Leishmania infantum Genetic Diversity and Lutzomyia longipalpis Mitochondrial Haplotypes in Brazil

Leishmania infantum is the etiological agent of visceral leishmaniasis (VL) in the Americas with domestic dogs being its major reservoir hosts. The main VL vector is the sandfly Lutzomyia longipalpis, while other Lutzomyia species may play a role in disease transmission. Although the genetic structure of L. infantum populations has been widely evaluated, only a few studies have addressed this subject coupled to the genetic structure of the respective sandfly vectors. In this study, we analyzed the population structure of L. infantum in three major VL endemic areas in Brazil and associated it with Lutzomyia longipalpis geographic structure.


Introduction
Leishmaniases are parasitic diseases caused by protozoans from the genus Leishmania, which are transmitted by the bite of female sandflies from the family Psychodidae. The clinical manifestations of leishmaniases are particularly diverse and present different characteristics: visceral leishmaniasis (VL), the most severe one; mucocutaneous leishmaniasis, characterized as a mutilating disease; diffuse cutaneous leishmaniasis, caused by a deficient cellular immune response; and cutaneous leishmaniasis, which causes single or multiple lesions on the skin. The epidemiology of leishmaniasis is highly complex: there are 20 known species of Leishmania pathogenic to humans and at least 30 species of sandflies vectors. Furthermore, this disease can be designated as a zoonosis, which involves animals as the reservoir hosts or as an anthroponosis, when humans are the only source of parasites for sandflies. Leishmaniasis is widely spread in 98 countries and 3 territories, from which more than 70% are developing countries and 13 are among the least developed ones [1].
Visceral leishmaniasis can be either an anthroponosis (e.g., in the Indian subcontinent) or a zoonosis (e.g., in the Mediterranean or in the Americas), and it is characterized by chronic evolution and systemic involvement, which if untreated may result in death. In the Americas, Leishmania infantum is the etiological agent of the disease and Brazil accounts for over 90% of the cases in the continent [1,2]. Domestic dogs are the proven reservoir hosts in rural and urban areas, while the role of naturally infected wild mammals (canids and marsupials) as L. infantum reservoir hosts is still controversial [3]. The main sandfly vector is Lutzomyia longipalpis, but other Lutzomyia species might play a role in disease transmission; for example, in Corumbá, Mato Grosso do Sul, naturally Leishmania-infected Lu. cruzi have been discovered and because there is still no evidence of Lu. longipalpis in this region, that sandfly is considered the main vector [4,5].
In Brazil, VL typically occurred in rural settings, but since 1980 its incidence has been changing due to widespread urban outbreaks. The first major VL urban epidemic in the country happened in Teresina, Piauí State. Since then, epidemics occurred in Natal (Rio Grande do Norte) and São Luís (Maranhão), and the disease subsequently spread to other regions of the country. Autochthonous cases were recently described for the first time in the southernmost State of Rio Grande do Sul. The current epidemiological scenario of VL leaves no doubt regarding the severity of the situation and the unchecked geographic spread of the disease. In the 1990s, only 10% of the cases occurred outside the Northeast Region, but in 2007 the proportion reached 50% of cases. From 2006 to 2008, autochthonous transmission of VL was reported in more than 1,200 municipalities in 21 states [6].
The broad spectrum of leishmaniasis-associated symptoms, coupled with the wide diversity of vertebrate and invertebrate host species, suggests that both parasites' and hosts' genetic backgrounds determine the patterns of the disease [7]. On the other hand, clonal diversity and genetic heterogeneity, which can cause variability in parasite virulence, are quite common in Leishmania [8].
Several studies showed that genetic variability of L. infantum in Brazil is low, with restricted diversity and limited population clustering. In a recent work assessing parasite populations distributed over 18 states, three major clustered populations could be inferred using microsatellite typing. When the analysis is performed in parasites from closely related geographic regions, the overall diversity is even lower [9][10][11].
When we look at sandfly genetic variability, there is compelling evidence that the Lutzomyia population structure in Brazil is complex, with different genotypes identified depending on the geographic region assessed and also the species involved in parasite transmission [12][13][14].
Based on these studies, it is logical to hypothesize that the interactions of L. infantum genotypic variants with different hosts and vector populations may ultimately influence the transmission dynamics and severity of eventual outbreaks. Hence, assessing the genetic structure of both vector populations and parasites may help us to understand the dynamics of vector-parasite interactions and the epidemiological aspects of American visceral leishmaniasis. Here, we used PCR-RFLP of kinetoplast minicircle DNA (kDNA) to identify L. infantum genotypic variants from three VL endemic areas in Brazil: Teresina in Piauí State, Campo Grande in Mato Grosso do Sul State, and Bauru in São Paulo State. kDNA-RFLP analysis when compared to microsatellite genotyping has proven to be more sensitive to examine genetic data of closely related sympatric L. infantum strains [15]. In addition, in order to identify different haplotypes of Lu. longipalpis and Lu. cruzi sandflies from those three VL endemic areas, we used mitochondrial 12S rDNA sequencing. As a maternal inheritance, rapidly evolving, nonrecombinant and haploid molecular marker, 12S rDNA is suitable to trace genealogy and evolutionary history. To our knowledge, this is the first study that seeks to compare genetic variability of Leishmania infantum parasites to the genetic structure of its vectors in Brazil.

Ethics Statement. For insect collections in Mato Grosso
do Sul State, we obtained a permanent license for collecting and transporting zoological material N ∘ 25592-1 on behalf of Dr. Alessandra Gutierrez de Oliveira, issued by the System of Authorization and Information on Biodiversity of the Brazilian Institute of Environment and Renewable Natural Resources (Sisbio/IBAMA). For insect collections in São Paulo State and Piauí State, no specific permissions were required since the specimens were kindly provided by the Center for the Control of Endemic Diseases (SUCEN) and Federal Piauí State University, respectively. The collections were performed at private residences, whose owners personally granted permission to enter their backyards to capture the sandflies. All of these residences were located in urban areas and no endangered or protected species were collected in this study.   All reactions were carried out in a thermal cycler, with 35 cycles of 95 ∘ C for 20 s, 50 ∘ C for 15 s, and 60 ∘ C for 2 min. The amplified DNA was precipitated with 80 L of 65% isopropanol, washed with 200 L of 70% ethanol, and airdried for 5 min. Before injection, samples were resuspended in 10 L of HI-DI formamide (Life Technologies) and heated at 95 ∘ C for 3 min for DNA denaturation and immediately cooled on ice. Sample processing occurred in an ABI377 automatic sequencer.

Sequencing
Analysis. The forward and reverse 12S rDNA sequences were manually checked for quality and the polymorphisms confirmed and then matched using the online EMBOSS GUI tool package (http://imed.med.ucm.es/ cgi-bin/emboss.pl? action=input& app=merger). The obtained consensus sequences were aligned using Clustal X2 software. Polymorphisms in each sequence were identified and a haplotypic diversity test (Table 1) was performed with the DnaSP 5.10 software. Haplotype diagram generation was performed by TCS: phylogenetic network using statistical estimation parsimony software.  Table 2). The DNA from the promastigotes (from all cultured samples used and for the two parasite samples obtained from sandflies) was isolated with Chelex (Bio-Rad). Briefly, 1 mL aliquots of the cultures were transferred to 1.5 mL centrifuge tubes and spun down for 1 min at 10,000 g. The supernatant was discarded and the pellet resuspended in 300 L of 10% Chelex (w/v). Following, the samples were incubated for 15 min at 95 ∘ C and then centrifuged again for 1 min at 10,000 g. The supernatant containing the DNA was then carefully recovered and stored in a new tube at −20 ∘ C. For the two parasite samples isolated from sandflies, the whole insect was grinded in 300 L of 10% Chelex with the help of a motorized tissue grinder, following the same steps above. We had an average of 200 ng of DNA per culture sampled and 20 ng per sample for the two sandfly-derived parasites measured with NanoDrop 1000 (Thermo Scientific).
The DNA of L. infantum amastigotes was extracted following two different approaches. For dog bone marrow aspirates we used the Illustra Blood GenomicPrep Mini Spin kit (GE Healthcare) according to the manufacturer's recommendations. For slide-fixed human bone marrow aspirates we used the same protocol after scraping the contents of each slide into a 1.5 mL tube, as previously described [17]. We had an average of 100 ng per dog bone marrow sample and 25 ng of DNA per slide measured with NanoDrop 1000 (Thermo Scientific).     [10]. However, only one marker (Li45-24) was polymorphic and, due to its low variability, only two alleles could be identified. For this reason, we decided to perform only PCR-RFLP of kinetoplast DNA (kDNA) and RFLP analysis. For the analysis of the kinetoplast minicircle DNA, 157 L. infantum isolates were used ( Table 2): 98 cultured samples initially isolated from human patients by sternal bone marrow aspiration (44 from Teresina and 54 from Campo Grande), 42 samples from dog bone marrow aspirates from Teresina, 2 samples from sandflies blood-fed on L. infantuminfected dogs from this same study in Teresina, and 15 slidederived samples originated from bone marrow aspirates of human patients in Bauru, São Paulo State. PCR reactions were performed with primers LINR4 and LIN19 [18] and generated a 720 bp amplicon, which covers almost the entire minicircle. The 50 L reactions contained 1 mM MgCl2, 10 mM Tris-HCl (pH 8.3), 0.3 pmol of each oligonucleotide, 0.1 mM dNTPs, 1 unit of Taq polymerase (GE Healthcare), and 5 L of sample DNA. The amplification conditions were as follows: 3 min at 94 ∘ C, 33 cycles of 30 s at 95 ∘ C, 30 s at 58 ∘ C, and 1 min at 72 ∘ C, followed by a final extension step of 10 minutes at 72 ∘ C. The PCR products were then precipitated with ethanol, resuspended in water and digested with the restriction enzymes RsaI and HpaII (Promega) as previously described [15]. Approximately 1 g of each PCR product was used per digestion in order to ensure that all reactions had the same initial amount of DNA. Since the products smaller than 100 bp can be confused with primer dimers and the ones larger than 700 bp can be misidentified as undigested products, only the fragments within this range were used in our RFLP analysis.
Data analysis was performed using R software environment. A binary matrix was constructed based on the profile of fragments generated by each digestion, where 1 represents the presence of a fragment and 0 represents its absence. This matrix was converted into a similarity matrix using the package "proxy" and used for cluster analysis. After, -means partitioning method was used to infer the number of clusters using the package " -means" and Agglomerative Hierarchical Clustering dendrogram was built using the binary distance method and ward cluster method with the package "hclust".   (Figure 1). PCR reactions generated a mitochondrial 12S ribosomal DNA fragment of approximately 360 bp, as previously described [19], which was then partially sequenced (263 bp). Sequences were screened for significant polymorphisms, and 10 variable sites were found (Table 3). When polymorphisms were assessed with DnaSP 5.10 program, 13 haplotypes were generated: six haplotypes (H8, H9, H10, H11, H12, and H13) containing only individuals from Teresina (PI); five haplotypes (H3, H4, H5, H6, and H7) containing only Araçatuba individuals (SP); one haplotype (H1) containing one individual from Corumbá and one individual from Campo Grande (MS); and one haplotype Haplotypes SNPs H1

Parasites RFLP Analysis.
The kDNA fragments of interest were successfully amplified from the LinR4 and Lin19 oligos used in this study. RFLP analysis of kinetoplast minicircles DNA was also efficient in detecting restriction patterns between different samples. From the 157 tested samples, we could observe 55 unique genotypes in the cluster analysis dendrogram illustrated in Figure 3. -means partitioning identified 6 major clusters; there was a clear distinction between samples from Teresina, which grouped in two almost exclusive clusters, and all other samples; an exclusive Bauru cluster was also found. Two clusters presented with Teresina and Campo Grande samples, and one cluster presented with Bauru and Campo Grande samples. It is noteworthy that Campo Grande is distributed over 3 major clusters, one that groups together with one Teresina major cluster and the other two that are closer to Bauru clusters in a separate branch of the dendrogram. There was no clustering differentiation related to the years of collection.

Discussion
During the past 20 years, the epidemiology of VL has been constantly changing due to a continuous urbanization process, an increasing incidence of HIV/Leishmania coinfections, and syringe sharing by intravenous drug users [20] and the identification of novel L. infantum mammalian hosts/reservoirs [21]. This highlights the necessity of molecularly tracking the geographic distribution of different parasite and vector populations in order to enhance the knowledge on basic epidemiological aspects of the disease, such as its natural history and transmission. Several molecular approaches have been used in the characterization of genetic variants in the genus Leishmania: amplified polymorphic DNA (RAPD) markers [22], analysis by size polymorphism of restriction fragments (RFLP) of the ITS regions ribosomal DNA [23], and kinetoplast DNA [24]; analysis confirmed sequence amplified regions [25]; and analysis of regions of DNA with microsatellite markers [19,[26][27][28][29]. We then decided to proceed with PCR-RFLP analysis of minicircle kDNA because it has a higher resolving power when applied to population genetics studies involving either genetically or geographically closely related strains [24,30,31]. Our data revealed a clear distinction between samples from Teresina, which grouped in two almost exclusive clusters, and all other samples; an exclusive Bauru cluster was also found. Two clusters presented with Teresina and Campo Grande samples, and one cluster presented with Bauru and Campo Grande samples. These results allowed us to draw a relationship between genetic distance and geographic origin. Interestingly, geographic origin related to diverse genetic background was also found for L. infantum parasites in Brazil in the study performed by Segatto et al. [10].
Our data is partially in accordance with a previous microsatellite based genotyping study performed with parasite populations from all 5 Brazilian regions. In the study, three well-defined populations could be identified; one that  TER83  TER80  TER33  TER25  TER39  TER84  CGR91  TER07  TER35  TER22  TER10  TER01  TER12  TER04  TER25  TER88  TER06  TER09  TER85  TER86  TER87  TER76  TER77  TER81  TER73  TER74  TER75  TER70  TER71  TER72  TER67  TER68  TER69  TER64  TER65  TER66  TER61  TER62  TER63  TER59  TER58  TER60  TER55  TER56  TER57  TER52  TER53  TER54  TER49  TER50  TER51  TER46  TER47  TER48  TER43  TER44  TER45  TER14  TER15  TER24  TER02  TER11  TER16  TER21  TER30  TER03  TER36  TER08  CGR89  TER26  TER32  TER31  TER34  TER27  TER37  TER41  TER38  TER05  TER29  TER42  TER17  TER13  TER19  TER40  TER28  TER20  CGR109  CGR138  CGR90  TER82  TER79  CGR120  CGR125  CGR100  CGR97  CGR106  CGR123  CGR135  CGR137  CGR136  CGR117  BAU144  BAU145  BAU146  CGR93  CGR92  BAU143  CGR96  CGR95  CGR94  CGR107  TER16  TER78  CGR101  CGR102  CGR114  BAU155  BAU156  CGR111  BAU149  BAU151  BAU152  CGR131  BAU147  BAU148  CGR121  CGR113  CGR122  CGR118  CGR130  CGR104  CGR105  CGR103  CGR116  CGR124  CGR127  CGR115  CGR141  CGR140  CGR126  CGR99  CGR96  CGR110  CGR134  CGR142  CGR119  CGR112  CGR133  CGR132  BAU154  BAU157  CGR108  CGR139  BAU150  BAU153 CGR128 CGR129 was present mostly in Northeast region, (including Piauí State that was sampled in our study) and the other two present in Midwest region (including Mato Grosso and Mato Grosso do Sul States that were sampled in our study). On the other hand, parasites typed in Southeast region (including São Paulo State that was sampled in our study) are closely related to northeastern parasites while in our study they are closely related to Midwestern parasites [9]. Our findings corroborate the use of this technique in Leishmania genotyping studies and reinforce the idea that in some cases, especially when analyzing strains of very close geographical origin, it is the only molecular marker capable of producing detectable patterns of polymorphism [24,32]. All these genotyping studies on L. infantum suggest that the nuclear genomic variability of this species is likely to be low. Our hypothesis is that the kinetoplast genome can serve as a source of genetic variability for these parasites. The kDNA minicircles are essential for the function of the trypanosomatid's mitochondrial genes, as minicircles code for guide RNAs, which play an essential role in editing messenger RNA (mRNA) from the maxicircles that contain genes for essential mitochondrial proteins [33]. Therefore, this DNA is more prone to a rapid response to diverse ambient conditions and stress situations, and parasite fitness conferring different selective advantages might depend on which minicircle classes prevail in different Leishmania strains.
A similar phenomenon, known as transkinetoplastidy, has been described in Leishmania and is responsible for changes in minicircles classes when the parasites are challenged with increasing concentrations of drugs that are normally lethal. This will in turn cause a dramatic change in the abundance of certain minicircles classes, which during transkinetoplastidy will be increased or reduced and replaced by a previously less frequent class [34].
When we look at sandfly genetic analysis we can clearly observe a main haplotype (H2) comprising all individuals from Andradina and Birigui, 13 out of 14 individuals from Campo Grande, 6 out of 7 individuals from Corumbá, 20 out of 30 individuals from Araçatuba, and 11 out 30 individuals from Teresina. There is also a major haplotype (H8) comprising only individuals from Teresina (13 out of 29) and minor haplotypes from Araçatuba. From the 12S rDNA sequencing data, it was not possible to differentiate Lu. longipalpis from Lu. cruzi (Corumbá) since there was no haplotype clustering among Corumbá sandflies. This may suggest that the process of speciation is recent or still occurring. A microsatellite based study assessing the genetic variability of Lu. longipalpis and Lu. cruzi populations in Mato Grosso do Sul State showed evidence of introgression and limited gene flow between the two species, corroborating our findings [12].
In general, we can summarize the data obtained from haplotyping as follows: a major haplotype composed of 111 individuals (comprising 89% of SP, 90% of MS, and 38% of PI individuals); a main haplotype composed of 13 individuals exclusively from Teresina and giving rise to other 4 Teresina exclusive haplotypes (62% of individuals from Teresina with exclusive haplotypes); minor haplotypes comprising only individuals from SP (11% total) and from the same locality (Araçatuba).
When we compare data from parasite genotyping with sandfly 12S rDNA sequencing, the correlation of the two datasets is remarkable. Both show most samples from PI clearly separated from the MS and SP ones which are in turn much more related to each other when compared to PI that presented the highest haplotypic diversity ( Table 1). The exception comes from the minor vector haplotypes only found in Araçatuba samples. Araçatuba represents an important landmark in the natural history of VL in SP given the fact that the first VL outbreak registered in the state occurred in this location [35,36]. This could be a possible explanation to its greater number of unique haplotypes as one can assume that coevolution between parasites and vectors happens for a longer time in this area; this is supported by the high haplotypic diversity found for this population (Table 1). Taken together, these data corroborate that the sandfly vector probably plays an important role in shaping the genetic structure of L. infantum in Brazil as described by Ferreira et al. [9].
This work presents new insights towards the understanding of the population structure of L. infantum and Lu. longipalpis from VL endemic areas in Brazil. Further analyses will be needed to elucidate how different vector populations shape the genetic variability of L. infantum.

Conclusions
Taken together, our data indicate that the sandfly vector might play a role in selecting specific parasite strains at a regional level and therefore contributing to the genetic structure of L. infantum in Brazil. Assessing the genetic structure of both vector and parasite populations may help us to understand the evolution process surrounding vectorparasite interactions and shed light on a fundamental aspect of the ecoepidemiology of American visceral leishmaniasis.