Phylogenetic relationships of closely-related phlebotomine sand flies (Diptera: Psychodidae) of Nyssomyia genus and Lutzomyia subgenus

Phlebotomine sand ?ies (Diptera: Psychodidae) are involved in the transmission of arboviruses, bacteria and protozoans. Among these pathogens, the trypanosomatids of the genus Leishmania spp. Ross are the etiological agents of the leishmaniases, a zoonotic disease considered as neglected with special relevance in public health.⁽¹⁾ Over 500 phlebotomine sand flies species distributed among 23 genera are known to occur in the Americas.⁽²⁾ Inside Nyssomyia genus, Nyssomyia intermedia (Lutz & Neiva), Nyssomyia whitmani (Antunes & Coutinho) and Nyssomyia neivai (Pinto) are Leishmania braziliensis Vianna vectors.^{(1, 3)} Nyssomyia whitmani is also a Leishmania shawi Lainson, de Souza, Póvoa, Ishikawa & Silveira vector and is suspected to transmit Leishmania guyanensis Floch.^{(3, 4)} Nyssomyia intermedia and Ny. whitmani are widely distributed in Brazil, including the five geographic regions of the country, but the latter is also registered in Suriname, French Guiana, Perú, Bolivia, Paraguay and northeast of Argentina. Whereas Ny. neivai is restricted to colder and drier areas at Bolivia, Paraguay, North, West-Central, Southeast and South regions of Brazil, and north of Argentina.^{(2, 3)} Another medical significant species is Lutzomyia longipalpis s.l. (Lutz & Neiva), the most important neotropical vector of Leishmania infantum Nicolle, causative agent of visceral leishmaniasis.⁽⁵⁾ Its geographical distribution comprises a wide but discontinued range from central Mexico to Argentina and Uruguay, with a recently southern expansion and urban environments colonisation.^{(5, 6)} A variety of evidence including chemical and behavioral suggests that it is a complex of at least seven species in Brazil, while molecular analyses identified eight haplogroups throughout its distribution.^{(5, 7, 8)} Regarding Lutzomyia cruzi (Mangabeira), its sibling, is also incriminated in Brazil as a vector of Le. infantum and suspected for transmitting Leishmania amazonensis Lainson & Shaw.⁽⁹⁾ While Lutzomyia alencari Martins, Souza & Falcão, has not been incriminated as a vector of Leishmania species.⁽³⁾ The geographic distribution of Lu. cruzi includes Bolivia and West-Central region of Brazil, while Lu. alencari is restricted to southeast of Brazil.^{(2, 3)}

Closely-related species but reproductively isolated that appear morphologically identical are often referred to as cryptic or sibling species, and a group of cryptic species is referred to as a species complex.⁽¹⁰⁾ The existence of sibling or closely-related species can hamper the process of species determination. Inside Nyssomyia genus, both males and females of Ny. intermedia, Ny. whitmani and Ny. neivai are closely-related species that can be determinate by subtle but sufficient morphological differences.⁽¹¹⁾ On the other hand, inside Lutzomyia subgenus, the males of Lu. longipalpis complex, Lu. alencari, Lu. cruzi, Lutzomyia pseudolongipalpis Arrivillaga & Feliciangeli and Lutzomyia matiasi Le Pont & Mollinedo presents very subtle details for their differentiation, while the female of Lu. longipalpis complex, Lu. cruzi and Lu. matiasi are morphologically indistinguishable together with Lutzomyia gaminarai (Cordero, Vogelsang & Cossio).⁽¹¹⁾ Therefore, females of these species have been identified based on their site of collection and through their association with males collected from the same location.⁽¹¹⁾ Molecular and biochemical findings suggest that reproductive barriers do not prevent gene flow among them, and evidence of introgression has been reported between Ny. intermedia and Ny. whitmani, and between Ny. neivai and Ny. intermedia.^{(12, 13)} Similarly, Lu. longipalpis complex and Lu. cruzi are not fully reproductively isolated with introgression occurring between them; the same is hypothesised for Lu. longipalpis complex and Lu. alencari in sympatric areas.^{(14, 15, 16)}

Mitochondrial DNA (mtDNA) markers, as COI gene, have been used for phylogenetic relationships inference, phylogeography, population genetics and molecular taxonomic identification (DNA barcoding) in phlebotomine sand flies from the new and old world.^{(15, 16, 17, 18, 19)} Given the medical importance of some phlebotomine sand ?ies within Nyssomyia genus and Lutzomyia subgenus, the studies of phylogenetic relationships will increase our understanding of the evolutionary relationships within these closely-related species, which to date has not yet been explored. Gene flow between these phlebotomine sand ?ies may influence vector competence and capacity and could have important epidemiological consequences such as insecticide resistance or susceptibility to Leishmania infection.⁽¹³⁾ Therefore, the aim in the present study was to infer the phylogenetic relationships among closely-related phlebotomine sand flies within Nyssomyia genus (Ny. intermedia, Ny. whitmani and Ny. neivai) and within Lutzomyia subgenus (Lu. longipalpis complex, Lu. cruzi and Lu. alencari) using a COI fragment.

MATERIALS AND METHODS

Study area, collection and taxonomic identification - The Leishmaniasis Research Network in Argentina (REDILA) captured phlebotomine sand flies as part of different epidemiological projects in 11 localities from six Argentinean provinces: Corrientes, Entre Ríos, Formosa, Misiones, Salta and Tucumán during 2015-2018 [Fig. 1, Supplementary data (Table I)]. Phlebotomine sand flies captured using REDILA-BL ultraviolet light traps, were separated from other insects and were preserved in 90% ethanol until being processed. For identification purposes, the head and the last three segments were dissected individually, clarified and mounted on slides while the rest of the body was stored at -20ºC until DNA extraction. Species determination was based on morphological characters according to Galati´s taxonomic key, and genera abbreviations follow the nomenclature proposed by Marcondes.^{(11, 20)} Reference specimens have been deposited at the Instituto Nacional de Medicina Tropical, Puerto Iguazú, Argentina.

DNA extraction, ampli?cation and sequencing - Genomic DNA was extracted individually from the rest of the body of the phlebotomine sand flies identified as Lu. longipalpis complex, Ny. neivai, Ny. whitmani and one Migonemyia migonei (França). We used DNA PuriPrep-S kit HighWay (INBIO, Argentina), following manufacturer’s instruction, with the addition of an incubation period at 56ºC that consisted in 4 h with vortex and spin every 20 min. The COI fragment amplification was done by a polymerase chain reaction (PCR) using the universal primers.⁽²¹⁾ The amplification conditions were an initial denaturalisation step of 94ºC for 4 min, followed by 35 cycles of 94ºC for 30 s, 56ºC for 30 s, 72ºC for 1 min, and a final extension of 72ºC for 10 min. PCR products were separated in 1.5% agarose gel electrophoresis, stained with 0.5 µL/mL SYBR Safe DNA (Invitrogen, USA) and isualized under blue light. The PCR products were purified with DNA PuriPrep-GP HW kit (INBIO, Argentina) following manufacturer’s instruction. The sequencing was carried out using both forward and reverse primers on an ABI 3730XLs (Macrogen, Korea) or was sent to the Canadian Centre for DNA Barcoding.

Genetic analysis - Cytochrome c oxidase subunit 1 (COI) sequences generated in this study of Ny. whitmani, Ny. neivai and Lu. longipalpis complex, (accession numbers in GenBank: MN857507-MN857540) in addition to sequences available from GenBank: Ny. intermedia (KP112715-KP112761), Ny. whitmani (KP112762-KP112770; MG234392-MG234412), Lu. longipalpis complex (JN845539- JN845542; KP112581-KP112595; KT806402-KT806405, KT806411, KT806413-KT806415, KT806417, KT806420, KT806426, KT806427, KT806430, KT806431, KT806437, KT806440, KT806443, KT806447, KT806451, KT806460, KT806468; GU909502-GU909505; HNLUZ004-HNLUZ005; MG234349-MG234370. BOLD Systems: MEXSM002-12.COI-5P, MEXSM003-12.COI-5P), Lu. cruzi (KP112572-KP112577) and Lu. alencari (KP112569-KP112571) were included in the analysis [Supplementary data (Table I)]. All sequences were manually aligned and edited using MEGA v.10 software. The multiple alignment of all sequences, was carried out in MEGA v.10 software. The number of polymorphic sites by species, the number of haplotypes and nucleotide divergence (Da) with the Jukes and Cantor (JC) correction among closely-related species were estimated using DnaSP v.6 software. Similarly, Kimura 2-parameter (K2P) genetic distances were estimated with 1000 replicates using MEGA v.10 software. The values of Da and K2P between Lu. longipalpis complex clades (see results) were also estimated.

Haplotype network of closely-related species was constructed using Network v.4.6.1.3, with the Median Joining algorithm to evaluate the relationships among haplotypes shared and the frequency distribution among species. The best-fitting model of evolution was estimated using the Bayesian Information Criterion (BIC) as implemented in JModeltest v.0.1.1 software. Bayesian inference (BI) was analysed with MrBayes v.3.2 software, four Markov Chain Monte Carlo (MCMC) were run for 10 million generations (sampled every 1000 generations) to allow adequate time for convergence (all ? 0.004369). The ?rst 25% of sampled trees were considered as burned in. The tree was visualised using FigTree v.1.4 (http://tree.bio.ed.ac.uk/software/?gtree/). Sequences of Mg. migonei (access number: MT430890) and Phlebotomus papatasi Scopoli (access number: KY848828) were used as outgroups.

RESULTS

Fifty-three phlebotomine sand flies were submitted to DNA extraction and COI amplification. Forty-eigth resulted in enough amplified product for sequencing (25 Nyssomyia genus and 23 Lutzomyia subgenus). Unfortunately, some of the sequences downloaded from GenBank and BOLD System were shorten than ours (658 bp). Then, during the multiple alignment all sequences were cut to the same length (543 bp for Nyssomyia genus and 483 bp for Lutzomyia subgenus).

Nyssomyia genus: Ny. intermedia, Ny. whitmani and Ny. neivai - A total of 103 COI sequences of 543 bp were analysed: 47 sequences of Ny. intermedia, 45 of Ny. whitmani and 11 of Ny. neivai. The species with more polymorphic sites was Ny. whitmani (44 sites, 8.1%). The values of K2P genetic distance ranged from 3.1 to 3.5%, while, the nucleotide divergence (Da) ranged from 0.017 to 0.025. The maximum and minimun K2P and Da values were between Ny. intermedia vs Ny. neivai and Ny. whitmani vs Ny. neivai, respectively (Table I).

Seventy-one haplotypes were identified from 103 Nyssomyia sequences, none haplotype was shared among them. Twenty-two haplotypes were identified in Ny. intermedia from localities of Brazil, 38 in Ny. whitmani from localities of Brazil and Argentina and 11 in Ny. neivai from Argentinean localities. The median-joining network identified at least 11 mutational steps between Ny. intermedia vs Ny. neivai and Ny. whitmani vs Ny. neivai,while between Ny intermedia vs Ny. whitmani there were at least 22 mutational steps. Argentinean Ny. whitmani (H61) from the northeast (Gobernador Lanusse and Piñalito Sur, Misiones provinces) and Pampeana (Concordia, Entre Ríos province) regions were identified by at least three mutational steps from Brazilian specimens (H58, Maranhão, Mato Grosso and Espírito Santo states) (Fig. 2). Bayesian inference (BI) was constructed using the HKY+I+G model as the most appropriate for the data (-lnL = 2152.4292; BIC = 5243.1277; Delta BIC = 0) with gamma of 0.1080 and p-inv of 0.6680. The analysis provided support with 0.77 posterior probability (PP) for the clade of Ny. intermedia and Ny. whitmani. Nyssomyia intermedia specimens grouped monophyletically with 0.96 PP, while some Ny. whitmani specimens also grouped monophyletically with 0.92 PP and another group presented polytomies (a comb-shaped phylogeny). Likewise, the basal clade with Ny. neivai specimens presented polytomies (Fig. 3).

Lutzomyia subgenus: Lu. longipalpis complex, Lu. cruzi and Lu. alencari - A total of 103 sequences of a 483 bp COI fragment were analysed: 94 sequences of Lu. longipalpis complex, six of Lu. cruzi and three Lu. alencari. Unfortunately, there were no sequences available in GenBank of Lu. pseudolongipalpis, Lu. matiasi and Lu. gaminarai, which are also closely-related species within the Lutzomyia subgenus.

Lutzomyia longipalpis complex showed the highest number of polymorphic sites (107 sites, 22%). The K2P genetic distance (ranged from 2.5 to 3.4%) and the nucleotide divergence (Da) (from 0.005 to 0.016) between closely-related species were low. The maximum and minimun K2P values were between Lu. longipalpis complex vs Lu. cruzi and Lu. longipalpis complex vs Lu. alencari, respectively. In contrast, the maximum and minimun Da values were between Lu. cruzi vs Lu. alencari and Lu. longipalpis complex vs Lu. cruzi and Lu. alencari (Table II). The intraspecific K2P distance and Da between clades from Lu. longipalpis complex (1 vs 2) (K2P = 9.8%, Da = 0.08) were high. Likewise, the values between clade 2 (samples from Colosó, Sucre) vs Lu. cruzi (K2P = 9.8%, Da = 0.081) and vs Lu. alencari (K2P = 9.6%, Da = 0.089) also were high [Supplementary data (Table II)].

Seventy-eight haplotypes were identified from 103 sequences analysed of the three closely-related species of Lutzomyia subgenus: 72 haplotypes of Lu. longipalpis complex, six of Lu. cruzi and three of Lu. alencari. The minimum spanning haplotype network showed three haplotypes shared among species: H5 was shared between Lu. longipalpis complex(Mato Grosso, Brazil and Tartagal, Argentina) and Lu. cruzi (Mato Grosso, Brazil); H12 was shared between Lu. longipalpis complexand Lu. cruzi (both from Mato Grosso, Brazil); and H59 was shared between Lu. longipalpis complex(Maranhão, Brazil) and Lu. alencari (Espírito Santo, Brazil). The specimens of Lu. longipalpis complex from Salta H6-H7, Tartagal H69 (Argentina) and Cáceres H11 (Brazil) are at least at one mutational step of Lu. cruzi H73 and H75 from Cáceres, Mato Grosso (Brazil). On the other hand, Lu. longipalpis complex specimens from Chiapas H71-H72 (Mexico) and Valle H26-H27 (Honduras) are separated by at least 24 and 26 mutational steps from Cundinamarca H28-H31 (Colombia, Andean region) and São Domingos do Maranhão H56 (Brazil) specimens, respectively (Fig. 4). Bayesian inference was constructed using the HKY+G model as the most appropriate for the data (-lnL = 2271.2317; BIC = 5852.6268; Delta BIC = 0) with gamma of 0.1040. The tree showed similar results to the network, its topology did not present structuration by species either by geography. Two monophyletic clades were identified: the first grouped specimens of Lu. longipalpis complex, Lu. cruzi and Lu. alencari (1.0 PP), while the second clade grouped only Lu. longipalpis complex specimens from Colosó, Colombia with 1.0 PP (Caribbean region) (Fig. 5). Indeed, Lu. longipalpis complex specimens from Colombia (Ricaurte, Cundinamarca and Colosó, Sucre) grouped within two monophyletic separated clades, the haplotype network identified at least 68 mutational steps between them (Figs 4-5).

DISCUSSION

The findings of the present study evidence low values of both K2P genetic distance and nucleotide divergence (Da) between closely-related species within Nyssomyia genus and within Lutzomyia subgenus. The genetically closest species within Nyssomya genus were Ny. whitmani and Ny. neivai; while within Lutzomyia subgenus were Lu. longipalpis complex with Lu. alencari andwith Lu. cruzi. Inside Nyssomyia genus, Mazzoni et al.^{(12, 22)} reported lower divergence values than us between Ny. intermedia and Ny. whitmani from southeast Brazil using ten nuclear loci. This can be due to: (a) Ny. neivai specimens were not analysed, (b) differences in nucleotide substitution rate between nuclear and mitochondrial genes, (c) the number of samples analysed, (d) the smaller geographical area they assessed or (e) a combination among them. Currently, de Melo et al.⁽²³⁾ reported even lower interspecific values inside Bichromomyia flaviscutellata complex between Bichromomyia olmeca olmeca (Vargas and Díaz-Nájera) and B. olmeca bicolor (Fairchild and Theodor). These authors concluded that these closely-related species represent the same species, supporting the morphological data. On the other hand, Scarpassa and Alencar,⁽¹⁸⁾ reported higher interspecific values (4.1-4.6%) between Nyssomyia anduzei (Rozeboom) and Nyssomyia umbratilis (Ward & Fraiha). Both are closely-related species with extensive overlapping areas and high morphological similarity. Arrivillaga et al.⁽¹⁷⁾ reported higher values of K2P distance ranging 1.3-17.7% within Lu. longipalpis complex while between Lu. longipalpis haplogroups (Hg) were 10.1-13%. In addition, K2P distance between Colombians specimens was comparable with the distance found between Lu. longipalpis complex clades (1 vs 2) in the present study. This could be explained by the geographical barriers in Colombia, such as the Andean Cordillera.^{(8, 17)}

The haplotype network and BI of Nyssomyia genus showed structuring by species although presented some polytomies, probably due to lack of samples. None haplotype was shared among species, contrary to mitochondrial introgression reported by Marcondes et al.⁽¹³⁾ between Ny. intermedia and Ny. whitmani as well as between Ny. intermedia and Ny. neivai using cyt b gene. Moreover, nuclear introgression using period gene was reported between Ny. intermedia and Ny. whitmani from southeast Brazil.^{(12, 22)} In the present study, haplotype network separated these three closely-related species by various mutational steps. Indeed, Ny. intermedia and Ny. whitmani have 22 mutational steps between them. All Ny. intermedia specimens from Espírito Santo (Brazil) presented few mutational steps among haplotypes, may be due to the restricted geographical range covered. Similarly, Meneses et al.⁽²⁴⁾ reported in Ny. intermedia population from Rio Janeiro (Brazil) few intra-population level polymorphism, using COI, 12S and 16S. In the BI, Ny. whitmani presented one monophyletic group (0.92 PP) and the rest of the specimens did not group (presenting polytomies). These not-grouped specimens from different localities of Brazil and Argentina are at least at seven mutational steps from the monophyletic group; suggesting that may form another monophyletic clade although more samples are needed to confirm this hypothesis. Indeed, haplotype network showed missing intermediate haplotypes between the two Ny. whitmani groups. Interestingly, Argentinian Ny. whitmani haplotypes from Concordia (its southernmost distribution until now) H61 and H62/H63 are separated by at least 17 mutational steps.⁽⁶⁾ Similarly, Ny. whitmani haplotypes from Puerto Iguazú (Northeast region, Argentina) have at least 11 mutational steps between them. These several mutational steps between the Ny. whitmani haplotypes from Argentina suggest that it could be a complex species occurring in sympatry. Ready et al.⁽⁴⁾ using cyt b gene suggested the presence of at least three phylogenetic lineages of Ny. whitmani in Brazil. Ishikawa et al.⁽²⁵⁾ adding samples of Ny. whitmani from the State of Rondônia (Brazil), identified two clades, one from eastern Amazônia and one from the Atlantic forest zone of the northeast, although with low support value (< 50.0% bootstrap). Unfortunately, in the present study, Ny. whitmani specimens did not cover the same locations previously studied by Ready et al.⁽⁴⁾ and Ishikama et al.⁽²⁵⁾ to compare the lineages. Regarding Ny. neivai, there are no studies available for comparisons, being these COI sequences the first available in GenBank.

Lutzomyia subgenus haplotype network and BI showed topologies without structuration by species neither by geography. The genetically closest species were Lu. longipalpis complex with Lu. alencari and with Lu. cruzi. Lutzomyia longipalpis complex shared haplotypes with both species or was found at least at one mutational step of them, compared with Lu. cruzi vs Lu. alencari that are separated by at least 16 mutational steps. Three haplotypes were shared among these species: (1) H5, Lu. longipalpis complex specimens from Argentina (Ar-Bra Hg) with Lu. longipalpis complex and Lu. cruzi from central-western Brazil, (2) H12, Lu. longipalpis complexwith Lu. cruzi (both from central-western Brazil) and (3) H59, Lu. longipalpis complex(northeast Brazil) with Lu. alencari (southeast Brazil). In this sense, the shared haplotypes detected both in symphatry and allophatrycally may be explained by introgression among these species. The introgressive hybridisation seems to be a recurrent phenomenon that shaped the genome of Lu. longipalpis complex.⁽²⁶⁾ This was suggested for Lu. longipalpis complex and Lu. cruzi based on molecular analyses using microsatellites markers⁽¹⁴⁾ and morphometrics analyses in sympatric areas.⁽²⁷⁾.

ACKNOWLEDGEMENTS

To Elizabet Vilacoba and Dario Lijtmaer from the Museo Argentino de Ciencias Naturales Bernardino Rivadavia. We are also grateful to Javier Liotta, Magali Giuliani for their assistance in the laboratory; Daniela Lamattina for providing samples (Instituto Nacional de Medicina Tropical - ANLIS Dr Carlos G Malbrán), and the workgroup of Concordia municipality.

AUTHORS’ CONTRIBUTION

Conceptualisation - AP-M and SLM; formal analysis - AP-M; funding acquisition - ODS, SLM, MM-A and MGQ; methodology - AP-M, SLM, MGQ and MM-A; writing/original draft - AP-M and SLM; writing/review and editing - AP-M, SLM, MGQ, MM-A and ODS.