Population Structure and Evidence for Both Clonality and Recombination among Brazilian Strains of the Subgenus Leishmania (Viannia)

Background/Objectives: Parasites of the subgenus Leishmania (Viannia) cause varying clinical symptoms ranging from cutaneous leishmaniases (CL) with single or few lesions, disseminated CL (DL) with multiple lesions to disfiguring forms of mucocutaneous leishmaniasis (MCL). In this population genetics study, 37 strains of L. (V.) guyanensis, 63 of L. (V.) braziliensis, four of L. (V.) shawi, six of L. (V.) lainsoni, seven of L. (V.) naiffi, one each of L. (V.) utingensis and L. (V.) lindenbergi, and one L. (V.) lainsoni/L. naiffi hybrid from different endemic foci in Brazil were examined for variation at 15 hyper-variable microsatellite markers. Methodology/Principal findings: The multilocus microsatellite profiles obtained for the 120 strains were analysed using both model- and distance-based methods. Significant genetic diversity was observed for all L. (Viannia) strains studied. The two cluster analysis approaches identified two principal genetic groups or populations, one consisting of strains of L. (V.) guyanensis from the Amazon region and the other of strains of L. (V.) braziliensis isolated along the Atlantic coast of Brazil. A third group comprised a heterogeneous assembly of species, including other strains of L. braziliensis isolated from the north of Brazil, which were extremely polymorphic. The latter strains seemed to be more closely related to those of L. (V.) shawi, L. (V.) naiffi, and L. (V.) lainsoni, also isolated in northern Brazilian foci. The MLMT approach identified an epidemic clone consisting of 13 strains of L. braziliensis from Minas Gerais, but evidence for recombination was obtained for the populations of L. (V.) braziliensis from the Atlantic coast and for L. (V.) guyanensis. Conclusions/Significance: Different levels of recombination versus clonality seem to occur within the subgenus L. (Viannia). Though clearly departing from panmixia, sporadic, but long-term sustained recombination might explain the tremendous genetic diversity and limited population structure found for such L. (Viannia) strains.


Introduction
The species of the subgenus Leishmania (Viannia) Lainson and Shaw, 1987, are exclusively endemic in the New World (NW) and infections of humans with these protozoan parasites constitute a significant public health problem in at least 18 countries of Latin America [1]. Subgenus L. (Viannia) parasites are capable of causing a variety of clinical symptoms ranging from cutaneous leishmaniasis (CL) with single or few lesions that may heal spontaneously, disseminated CL (DL) with multiple lesions, to disfiguring forms of mucocutaneous leishmaniasis (MCL) that may occur concomitantly or after remission of CL [2]. The outcome of human infections by Leishmania parasites is thought to be influenced by the immune status of the host and virulence of the infecting parasite [3]. At present, multilocus enzyme electrophoresis (MLEE) is the reference technique for the identification of Leishmania and was employed in most of the classification schemes, although MLEE is likely to be partially superseded by multilocus sequence typing (MLST). The application of numerical taxonomy and cladistic techniques to electrophoretic data has resulted in the identification of two species complexes in the subgenus L.  [4]). This classification has been largely supported by a recent molecular study comparing hsp70 sequences of different Leishmania species [5].
In Brazil, CL is endemic in all federal states and an annual mean of 27,250 CL cases has been registered from 1990-2010 (http://portal.saude.gov.br/portal/saude/profissional/area.cfm?id _area = 1560). The disease is caused by six species of the subgenus L.  [7,8]. So far, DL has been also mainly associated with L. (V.) braziliensis, and only sporadic cases with L. (V.) guyanensis [1]. Severe anergic diffuse cutaneous leishmaniasis (DCL) may be a long term sequel in a minority of L. (L.) amazonensis infections [1].
Transmission of species of the subgenus L. (Viannia) involves different species of phlebotomine sand flies and a wide variety of wild and domestic animals have been implicated as reservoir hosts. Some species have a more restricted transmission cycle, whereas others are more complex with several different vectors and hosts in different ecological and geographical regions [1,3]. Sympatry of different subgenus L. (Viannia) species has been reported particularly in the Amazon region [9] where separate epidemiological patterns have been described, involving different sand fly species [10]. There is also increasing evidence that pathogenic Leishmania strains can be maintained in both sylvatic cycles, involving wild animals and sylvatic sand flies, and urban cycles involving domestic animals and peridomestic sand flies [11]. Subgenus L.
(Viannia) parasites are characterized by tremendous genetic diversity and the described species vary enormously in their degree of such diversity [12,13]. Intra-specific polymorphisms are, for example, very frequent in L. (V.) braziliensis and L. (V) naiffi, and it has been suggested that the genetic diversity of the parasites is most probably related to the sand fly vector(s) and/or animal reservoir(s) involved in the transmission cycles. On the other hand, L. (V.) guyanensis, particularly strains circulating in the Brazilian Amazon region, and L. (V.) shawi have been found to be rather homogenous by MLEE and ITS-RFLP typing [12,13].
Multilocus microsatellite typing (MLMT) is currently the method of choice for molecular epidemiological and population genetic studies of different species of Leishmania (reviewed in [14]). It combines the advantages of co-dominance and higher discriminatory power when compared to MLEE, RAPD and the PCR-RFLP approaches used in many studies. Different sets of microsatellite markers have been designed and successfully applied for discriminating strains of subgenus L. (Viannia) with special emphasis on L. (V) braziliensis and L. (V.) guyanensis [15][16][17][18][19]. In a preliminary study we have demonstrated that our microsatellite marker set is highly discriminatory at intra-species level. Moreover, this genotyping scheme allows the detection of mild genetic structures at different levels, and is thus, relevant for epidemiological and population genetic studies of strains within the subgenus L. (Viannia) [17].
In the present study we investigated microsatellite variation in strains of the subgenus L.

Ethics statement
Research in this study was subject to ethical review by the European Commission and approved as part of contract negotiation for Project LeishEpiNetSA (contract 01547): the work conformed to all relevant European regulations. The research was also reviewed and approved by the ethics committee of the London School of Hygiene and Tropical Medicine (approval 5092). The Leishmania strains isolated from human and animal hosts and analysed in this microsatellite analysis, were received from the ''Coleção de Leishmania 34), and from the collection hosted by the Universidade Federal de Minas Gerais in Belo Horizonte. Only previously gathered samples from animals have been used in this study. All human strains of Leishmania had been isolated from patients as part of normal diagnosis and treatment with no unnecessary invasive procedures and with written and/or verbal consent recorded at the time of clinical examination. Data on human isolates were coded and anonymised.

Parasite and DNA samples
Sources, designation, geographical origins, MLEE identification, if known, and clinical manifestation for the 120 subgenus L. (Viannia) strains from Brazil that were used in this study are listed in Table S1. These Most of the strains were isolated from human CL cases, three from DL cases, and three strains from patients suffering from MCL. The reference strains were cloned, but all other strains represented uncloned material. Seven strains were isolated from sand fly vectors, and 18 from different animals, such as opossums (4), rodents (3), armadillos (3), dogs (2), sloths (2), pacas (2), a capuchin monkey (1) and a

Author Summary
Cutaneous leishmaniasis (CL) constitutes a significant public health problem in all federal states of Brazil. Most cases are caused by parasites of the subgenus Leishmania (Viannia) which can cause a variety of clinical symptoms ranging from single or few lesions, disseminated CL with multiple lesions, to disfiguring forms of mucocutaneous leishmaniasis. This study has used a multilocus microsatellite typing approach for exploring the genetic diversity and population structure among 120 strains representing different subgenus L. (Viannia) species and different Brazilian CL foci. Genetic diversity within the subgenus was much higher than expected, especially within L. (V.) braziliensis, L. (V.) shawi, L. (V.) naiffi, and L. (V.) lainsoni which were all from the north of Brazil. These strains could not be assigned to well-defined populations, but presented a rather loosely associated group. Strains of L. (V.) braziliensis isolated along the Atlantic coast of Brazil and strains of L. (V.) guyanensis formed, however, two clearly separated populations exhibiting remarkable levels of sexual exchange. The latter finding is in contrast to previous studies suggesting clonal modes of propagation or inbreeding for natural populations of Leishmania parasites and might explain the genetic heterogeneity and limited population structure for Brazilian strains of subgenus L. (Viannia) observed in this study.
porcupine (1). Table 1 summarises the number of strains per species according to geographical origin, zymodeme and clinical picture.
Most strains were obtained from the FIOCRUZ Leishmania collection (Coleção de Leishmania do Instituto Oswaldo Cruz, CLIOC, WDCM731, http://clioc.fiocruz.br). Seventeen strains of L. (V.) braziliensis, 15 from Minas Gerais and two from Pará, were obtained from the collection of the Universidade Federal de Minas Gerais, Belo Horizonte. Sample preparation and MLEE typing, based on the electrophoretic mobility of 11 enzymes in agarose gel electrophoresis, were performed as previously described [12].
DNA was isolated using proteinase K-phenol/chloroform extraction [20] or the Wizard TM Genomic DNA Purification System (Promega, Madison, WI, USA) according to the manufacturer's protocol, suspended in TE-buffer or distilled water and stored at 4uC until use.

PCR amplification assays and electrophoretic analysis of the microsatellite markers
The standard set of 15 primer pairs (CSg46, CSg47, CSg48, CSg53, CSg55, CSg59, 6F, 7G, 10F, 11H, 11C, B3H, B6F, AC01R, and AC16R), specific for L. (Viannia), was used for amplification of microsatellite containing fragments, as previously described [17]. PCRs were performed with fluorescenceconjugated forward primers. Screening of length variations of the amplified markers was done by automated fragment analysis using the ABI PRISM GeneMapper (Applied Biosystems, Foster City, CA). After manual checking the microsatellite repeat numbers were calculated for all loci by comparing the sizes of the respective fragments to those of the strains MHOM/BR/00/LTB300 (L. braziliensis) and MHOM/SR/87/TRUUS1 (L. guyanensis), which were included as reference strains in every experiment, and for which the repeat numbers had been determined by sequencing. These repeat numbers were then multiplied by two, since we have used dinucleotide microsatellites throughout, after which the size of the flanking region was added, as determined by sequencing of the reference strains. This rigorous normalization process was applied to correct for small size differences that could occur due to the use of different sequencing machines and/or fluorescent dyes during the analyses. These normalized fragment sizes for all markers were assembled into a multilocus microsatellite profile for every strain under study. In 1.7% of all loci three or four peaks were observed of which only the two most prominent bands were included in the microsatellite profiles. The microsatellite profiles of all strains analysed in this study are given in Table S1.

Data analysis
Population structure was investigated using the STRUCTURE software [21], which applies a Bayesian model-based clustering approach. This algorithm identifies genetically distinct populations on the basis of allele frequencies. Genetic clusters are constructed from the genotypes identified, estimating for each strain the fraction of its genotype that belongs to each cluster. The following parameters were used: ''burn-in'' period of 20,000 iterations, probability estimates based on 200,000 Markov Chain Monte Carlo iterations. The most appropriate number of populations was determined by comparing log likelihoods for values of K between 1 and 10, with ten runs performed for each K, and by calculating DK, which is based on the rate of change in the log probability of data between successive K values [22]. Factorial correspondence analysis (FCA) implemented in GENETIX v 4.03 software [23] was performed, which places the individuals in a three-dimensional space according to the degree of their allelic state similarities.
Phylogenetic analysis was based on microsatellite genetic distances, calculated with the program POPULATIONS 1.2.28 (http:// bioinformatics.org/,tryphon/populations) for the numbers of repeats within each locus using the Chord-distance [24], which follows the infinite allele model (IAM). Neighbor-joining trees (NJ) were constructed with the POPULATIONS software and visualized with MEGA [25]. Additionally, phylogenetic networks were inferred from the distance matrix obtained from the microsatellite dataset by using the Neighbor-Net method in SplitsTree4 [26].
Microsatellite markers as well as populations were analysed with respect to diversity of alleles (A), expected (gene diversity) and observed heterozygosity (H e and H o , respectively) applying GDA (http://hydrodictyon. eeb.uconn.edu/people/plewis/software. php). F IS is a measure of heterozygosity that assesses the level of identity within individuals compared to that between individuals. It ranges between 21 and 1, where a negative value corresponds to an excess of heterozygotes, and a positive value to heterozygote deficiency. F IS = 0 indicates Hardy-Weinberg allele proportions. Mean F IS estimates over loci in each population were calculated with the software FSTAT (version 2.9.3.2) [27] using Weir and Cockerman's (1984) unbiased estimators [28]. Confidence intervals per locus were assessed by randomization and bootsraping procedures over loci and individuals, implemented in GENETIX [23] using 1,000 random permutations. We also analysed the data, by computing estimates and tests of significance for various population genetic parameters. Genetic differentiation and gene flow was assessed by F-statistics [28][29][30] with the corresponding P-values (confidence test) using the MSA software [31]. Linkage between all pairs of loci in populations 1 and 2 was tested using the software ARLEQUIN, version 3.5 [32] and FSTAT [27]. P-values for multiple tests were corrected using a sequential Bonferroni correction to minimize the likelihood of Type 1 errors [33]. Composite digenic disequilibrium values were estimated and their significance was tested using Chi-square statistics as described by Weir [34]. An exact test for association between alleles across loci based on permutation [35] was also employed.
To assess the level of multilocus linkage disequilibrium, the Index of Association (I A , multilocus) and the r d statistic were calculated in MULTILOCUS 1.3b [36,37]. P values were derived through comparison to a null distribution of 1,000 randomizations; median values were taken from 1,000 diploid resamplings of the multiallelic dataset.

Results
Genetic diversity of Brazilian strains of subgenus L.

(Viannia)
All 15 microsatellite markers were polymorphic in the 120 strains of the subgenus L. (Viannia) analysed here ( Table 2). Taking together all strains, the number of alleles varied between 7-29, with a mean value of 15, with markers CSg47 and CSg48 being the most variable. The overall observed heterozygosity per marker ranged from 0.117 to 0.706 and was lower than the expected (0.646-0.940) for all markers. Overall inbreeding coefficients varied between 0.245 and 0.855. The discrepancy between expected and observed heterozygosity and the high F IS values most probably reflects population substructuring (Wahlund effect) although the existence of a considerable amount of inbreeding cannot be ruled out. A total of 107 strains had unique MLMT profiles (Table S1). Five profiles, Lgua13, Lbra6, Lbra22, Lbra48 and Lsha3, were each shared by two strains. Nine of the 15 strains of L. braziliensis isolated between 1986 and 1992 from human CL cases in Minas Gerais presented indistinguishable MLMT profiles (Lbra26). Only two of the Minas Gerais strains were different from the predominating genotype ( Figure 1, Figure S1).
More than two peaks were found for 30 of the 1800 loci (1.7%) analysed in this study. Twenty-three of the 120 strains presented such possibly aneuploid loci. Strain L. braziliensis L 2516 had more than two peaks in four loci, and strains L. braziliensis L-0018, L. lainsoni L-2500 and L-2503, and L. naiffi L-991 in two loci each. One putative aneuploid locus was seen in five strains of L. guyanensis, eight of L. braziliensis, two of L. lainsoni, two of L. naiffi and in the L. naiffi/L. lainsoni hybrid L-2490.
Population structure of Brazilian strains of the subgenus L. (Viannia) The multilocus microsatellite profiles consisting of the repeat numbers for 15 markers were processed using both model-based and distance-based methods and the results of both analyses are compared in Figure 1 and Figure S1. STRUCTURE analysis assigned the 120 Brazilian strains of the subgenus L. (Viannia) to three main populations (Table S1, Figure  S2 (Table S1), and strain L-2446 from Pernambuco for all three populations indicating gene flow between the populations.
The existence of Populations 1 and 2 was supported by FCA. However, strains of Population 3 were shown to be highly heterogeneous compared to Populations 1 and 2 ( Figure 2). Two strains of Population 2, namely L-2516 and L-2446 both from Pernambuco, grouped within the cloud formed by Population 3. F-statistics revealed significant genetic differentiation between the three populations identified by STRUCTURE, especially between Populations 1 and 2 ( Table 3). The correlation between population assignment and the geographical origin of the strains is shown in Figure 3.
The   (Table S1). The distance analysis based on the MLMT profiles obtained for all individual strains and the inferred neighbor-joining tree ( Figure 1, Figure S1 (Table S1) and might represent hybrids or mixed infections as cloned isolates were not used.
The phylogenetic NeighborNet network (Figure 4) largely confirmed the results described above. It clearly showed the tremendous diversity of the strains assigned to Population 3 and its four sub-populations. However, conflicting splits represented by boxes can be seen between and within the three main populations. The same analysis was carried out for each of the populations separately. The obtained phylogenetic networks confirm most of the sub-structures found previously by Bayesian analysis ( Figures  S3, S4, S5).

Population genetics characterization of identified populations
The difference in the degree of microsatellite polymorphism between the three main populations was also reflected by mean number of alleles (MNA) which fluctuated from 4.1 to 4.9 in Populations 1 and 2, respectively, up to 14.2 in population 3 (Table 4). In Populations 1 and 2 the observed heterozygosities (Table 4) were close to those expected under Hardy-Weinberg equilibrium (HWE), whilst global F IS values were only moderately positive for Population 1 and Population 2 (Fig. S6). Significant (P,0.05) excess of heterozygosity indicated by negative F IS values was detected for only two loci, locus CSg47 in population 1 (20.551) and locus B3H in population 2. Another 8 loci presented significant deficit of heterozygosity in a least one of the two populations. Both global F IS values were significant (0.184 and 0.088 respectively; P = 0.0017 and P = 0.0167) and therefore panmixia must be rejected. Testing HWE within population 1 and 2 supported this scenario since both populations significantly departed from the null-hypothesis (P,0.002 and P,0.001, respectively). The Maynard Smith index of association, I A , [36,39], which assesses multilocus linkage disequilibrium, was calculated. Another estimator, r d , was implemented because the I A tends to increase with the number of loci, a trend corrected by this statistic. To summarise, both populations displayed positive I A and r d values (1.295 and 0.097 for population 1 and 1.586 and 0.123 for populations 2, respectively) departing significantly (P,0.001) from panmixia (I A and r d = 0). However, the F IS values gathered from populations 1 and 2 are consistently lower than those observed in previous reports [40,41], indicating that gene conversion or recombination may play a substantial role in Viannia species present in the Amazon Basin and along the Atlantic coast and to a lesser extent westward from the Andes. In contrast, in Population 3 the observed heterozygosity was much less than the expected resulting in a high F IS value, most probably due to population subdivision (Wahlund effect) although high rates of gene conversion or inbreeding cannot be excluded. This group was not tested for all population genetic parameters, since it represents a composite and artificial unit.
In order to test whether associations between the 15 microsatellite loci are in gametic equilibrium in populations 1 and 2, as expected for random-mating populations, we have applied both the composite disequilibrium test and the exact test to all 105 pairwise comparisons. For population 1, only one strain representing the microsatellite profile that was shared by two strains (Table S1) was included in the data matrix. Using Chi-square test and exact test, 24 and 37 significant associations were respectively detected. However, after Bonferroni correction, those numbers dropped to 6 and 10. In population 2 nine strains of L. braziliensis from Minas Gerais shared an identical genotype, Lbra26, and further four strains had highly related genotypes, Lbra 28, 29, 30 and 31. These strains most probably represent an epidemic clone and only one Lbra 26 strain was therefore included in the data matrix for linkage disequilibrium calculations. The Chi-square test and the exact test revealed 30 and 49 significant associations (28 and 45 after Bonferroni correction), respectively for population 2.
The three populations detected according to the Bayesian algorithm were also significantly separate from each other with highly significant F st values ranging from 0.249 to 0.521 (Table 3).

Discussion
Because of its high resolution potential, its reproducibility and the possibility of data storage, multilocus microsatellite typing  Figure S1) was calculated for the MLMT profiles of 120 strains of different species of subgenus L. (Viannia), based on 15 microsatellite markers and using the Chord distance measure. The assignment of these strains to three main populations by the Bayesian model-based clustering approach implemented in STRUCTURE is indicated by coloured circles: population 1 (green), population 2 (red) and population 3 (blue). Strains belonging to these populations are listed in Table S1 [40,41], and confirmed that the agent of VL in the NW is L. (L.) infantum, which has been recently imported multiple times from southwest Europe to the New World [42]. In the present study, a MLMT approach employing 15 microsatellite markers previously shown to be highly discriminatory for strains of the subgenus L. (Viannia) [17], was used to explore the genetic diversity of 120 strains from Brazil in order to unravel discrete populations. Amazonas. This last group seems, however to be rather artificial, since seven of the ten strains have hybrid genotypes sharing alleles that are specific for Population 2 and 3, respectively. Strains of L. (V.) braziliensis of population 3 are considerably different from each other and very distinct from those that were assigned to Population 2, except those having mixed memberships to Populations 2 and 3. Whether such strains represent outliers, as stated in a different study that used AFLP for typing strains of L. (V.) braziliensis and L. (V.) peruviana mainly from Peru and Bolivia [43], or different taxa requires further investigations including additional strains and using DNA sequence-based comparisons.  NeighborNet analyses provide a snapshot overview of the general structure in the data and are useful as a guide for further analysis [44]. The phylogenetic network obtained here for the full sample set is distinctly non-treelike and demonstrates marked ambiguity in the signal. Only Populations 1 and 2 formed distinct clusters in the network. The reticulate patterns seen in the network between (Figure 4), and within the three main populations ( Figures  S3, S4, S5) could result from hybridization, recombination events or gene conversion.  [12,13,38], strains of L. (V.) braziliensis from northeastern Brazil belonging to zymodeme Z75 were however, found to be related to L. (V.) shawi [46] which is in agreement with the results of our MLMT study. Assessment of the taxonomic status of L. shawi thus warrants further investigation with more extensive DNA sequence comparisons. All strains of population 3 and only those are found on long branches in the overall network. Because NeighborNet is prone to long-branch attraction, rapidly evolving lineages can be inferred as being closely related regardless of their evolutionary relationships. Whether these strains have a high mutation rate leading to numerous homoplasies or to convergence, which could be misinterpreted as having evolved once in a common ancestor, remains to be established.
Given its high genetic diversity, Population 3 could represent the ancestral lineage and might have given rise to two new populations through bottleneck events (Populations 1 and 2). The Amazon forest seems to be the central distribution area with secondary spreads to the northeast, east and south. This would be consistent with the hypothesis of an Amazon origin of CL in Brazil, with later spread to other regions, most probably through human migrations [6]. However, we cannot exclude that sampling biases are responsible for the weak resolution of the strains in Population 3. More extensive sampling in the north of Brazil, where most of these strains were isolated, is needed to address these questions.

Genetic diversity of Brazilian strains of the subgenus L. (Viannia)
Almost all strains investigated in this study presented unique microsatellite profiles, except 13 strains from Minas Gerais that had identical or highly similar microsatellite profiles and might have been isolated during an outbreak of CL in this area.
Previous studies using isoenzyme typing [12] or PCR-RFLP of the ribosomal internal transcribed spacer [13] have already demonstrated that L. (V.) braziliensis is much more polymorphic than L. (V.) guyanensis. Furthermore, strains of L. (V.) braziliensis from the Amazon region and from Pernambuco were shown to present the highest level of genetic diversity whereas those from Rio de Janeiro were more homogeneous [11,46]. In our study the MNA (Table 4)  The strains of L. (V.) guyanensis investigated herein were more diverse than expected considering that they all belonged to the same zymodeme Z23 which is in agreement the results of a recent MLST study [45]. Only two strains, MHOM/BR/1997/203P and MHOM/BR/1997/203G, isolated from the skin (P) and a lymph node (G) of the same patient, shared an identical MLMT profile ( Figure S1), all other strains presented different patterns of microsatellite variation. Previous studies using monoclonal antibodies had already pointed to the existence of two distinct sub-populations in the Brazilian Amazon region [47]. In our study, the existence of two L. (V.) guyanensis sub-populations in the Amazonas state (data not shown) was not supported by the genetic distance analyses. Our MLMT approach confirmed however, the existence of a new L. (V.) guyanensis genotype in Acre, where this species is not commonly found. This strain was very closely related to L. (V.) shawi as previously suggested [38]. In French Guyana, two distinct populations of L. (V.) guyanensis had been studied using a PCR-RFLP approach targeting ribosomal DNA sequences and found to have originated from two ecologically different regions and to differ in clinical manifestations of CL [48]. In a previous study comparing microsatellite variation in a limited number of strains of L. (V.) guyanensis from Brazil, Peru, Suriname and French Guyana, the strains grouped according to their geographical origin [17]. Future investigations should include strains sampled in different locations, since it has been speculated recently that strains from the eastern and southern Brazilian Amazon region might represent different genetic groups of L. (V.) guyanensis [8]. We found only weak correlations between the MLMT profiles and the results of previous isoenzyme typing (Table S1, Table S2). As already mentioned above, the strains of L. (V.) guyanensis, despite all being of zymodeme Z23, could be individualised by MLMT but were all grouped in Population 1, with the exception of the single strain from Acre. Strains of L. (V.) braziliensis with identical isoenzyme patterns were assigned to different populations or genetic groups by MLMT. The majority of the strains presenting the predominant zymodeme Z27 grouped in Population 2 but those with hybrid genotypes were found in Population 3B. This implies that zymodeme Z27 is paraphyletic and does not reflect the genetic diversity of the strains which was also shown by the MLST study [45].
We did not find any correlation between a particular MLMT profile and the clinical presentation of the disease. The two strains isolated from MCL patients in Bahia grouped together with strains from CL cases from the same area, although the long term outcome of those CL cases is not known. This is consistent with previous studies which suggested that the clinical outcome of the disease caused by L. (Viannia) parasites is also influenced by host genetic and/or immune factors [49][50][51] which could possibly be stimulated through pre-exposure to sand fly saliva [52]. In conclusion, the only correlation found for MLMT patterns of the L. (Viannia) strains studied herein was that to their geographical origin. Similar observations were published earlier for L. (V.) braziliensis [11] and L. (L.) infantum in Brazil [53] and might be associated with different transmission cycles with different sand fly vectors and/or animal reservoirs involved in those areas.

Reproductive strategies among Brazilian strains of the subgenus L. (Viannia)
Despite the fact that recombination in Leishmania has been proved to occur in vitro in the sand fly hosts [54,55] and the growing evidence of gene flow coming from different population genetic studies using MLST and MLMT approaches (reviewed in [14]), Leishmania species are still considered as predominantly clonal organisms [56]. Especially in the case of strains of the subgenus L. (Viannia) this hypothesis has been challenged by the frequent detection of hybrids involving different species of the L. (Viannia) subgenus indicating that recombination events are much more frequent in these parasites than previously thought [12,57] [62], and these hybrids had been mostly isolated from patients living in areas with sympatric circulation of both putative parental species. In our study, some strains had mixed membership of the different populations identified and were considered to be putative hybrids, although this will have to be confirmed by analysis of cloned parasites. A recently published MLSA study of subgenus L. (Viannia) strains has also provided evidence for recombination occurring in both L. braziliensis and L. guyanensis [45].
Strong linkage disequilibrium, or non-random association of genotypes at different loci, and a distinct phylogenetic signal are the criteria for the identification of clonality [56]. To start with the latter, strong phylogenetic signals are clearly absent in both the NJ distance tree (Figure 1 and Figure S1) and the NeighborNet network ( Figure 4). Calculation of linkage disequilibrium which does not depend on the ploidy status [56] revealed a higher number of significant associations between loci for both populations 1 and 2 than what would be expected due to chance for a random mating population. On the other hand, the majority of the pairwise comparisons were not significant and many more loci appear to be recombining than would be expected for a strictly clonal population. We can, however not exclude that population subdivision (Wahlund effect) accounts, at least partially, for the amount of disequilibrium found for populations 1 and 2. Taking in consideration the limited linkage disequilibrium, the absence of overrepresented genotypes and the weak phylogenetic signal observed in this study, we would conclude that recombination has an important impact in populations 1 and 2. The MLMT approach used here, however also identified an epidemic clone consisting of 13 strains of L. braziliensis isolated between 1986 and 1992 from human CL cases in Minas Gerais. Nine of these strains shared an identical microsatellite profile (Lgua26) and the profiles of the other four strains differed from profile Lgua26 in only one locus each.
Different levels of clonality versus recombination have been earlier suggested to occur within some bacterial and protozoan species due to variation in geographic sampling [36] and this might be also the case for species of the subgenus L. (Viannia). The substantial heterozygote deficiency and extreme inbreeding found by MLMT analysis of L. (V.) braziliensis from Bolivia and Peru was consistent with a predominantly endogamic mode of reproduction (mating with relatives) with occasional recombination events between individuals of different genotypes [40]. In contrast, significant homozygosity and only little linkage disequilibrium was observed for populations of L. (V.) guyanensis from French Guyana suggesting a high level of sexual recombination and substantial endogamy [41]. Ramirez et al. [63] stated that heterozygosity statistics at microsatellite loci has to be interpreted with caution in the context of parasite sexuality because strong linkage disequilibrium can be accompanied by negative and positive Fis values. In our study, mild linkage disequilibrium was observed together with relatively low Fis values, compared to those previously published [40,41]. Indeed the Fis values observed for the L. guyanensis strains from the Brazilian Amazon region (0.184) and for the L. braziliensis strains from Eastern Brazil (0.088) are below the values observed in former publications (0.278 and 0.307) and this is without partitioning steps to correct for a potential Wahlund effect, which might deflate this index. Consequently, the selfing rates for the considered populations are likely to be %0. 50. This implies that clonality, selfing and random union of gametes contribute to the shaping of Viannia's natural populations. However, one of the future challenges is to understand why the contribution of sex is more significant in the Amazon Basin.
The analyses used for calculations of heterozygosity are based on the assumption of diploidy. Recently, whole genome sequencing and FISH analyses have, however confirmed significant chromosomal copy number variations for different species of Leishmania [64][65][66][67][68]. In the only L. braziliensis strain, MHOM/BR/ 75//M2904, that has been fully sequenced so far, 30 of the 35 chromosomes were clearly triploid, three were tetrasomic and one hexasomic [66]. Whether other strains of L. braziliensis show similar or different ploidy patterns, as shown for all other Leishmania species examined so far [64,66,67], remains to be established. More than two alleles have been observed for only 1.7% of the microsatellite loci analysed in this study. This could possibly be due to aneuploidy [64,66], although other reasons such as mixed strains, duplication or stutter bands cannot be excluded. Sterkers at al. ( [64,65] were able to show that in L. major chromosomal content varies not only from strain to strain but also from cell to cell creating 'mosaic aneuploidy'. This leads to high karyotypic diversity and conserved intra-strain genetic heterogeneity combined with loss of heterozygosity per cell. The total number of alleles can, however be maintained in a strain. As a consequence, DNA-based typing methods, including the microsatellite typing approach used herein, cannot decide if a cell population (or strain) consists of heterozygous cells or of homozygous cells presenting different allelic and ploidy content [65].
In conclusion, this study showed the extensive microsatellite diversity present in the subgenus L. (Viannia) and indicated that L. lainsoni also isolated in northern Brazilian CL foci. All findings concerning the strains from northern Brazil may, however, be subject to bias due to an inadequate sampling strategy. More strains need to be sampled from this region in order to fine tune the population structure of these parasites and their mode of reproduction. Figure S1 Rectangular NJ tree showing the populations and subpopulations of 120 Brazilian L. (Viannia) strains. A midpoint rooted Neighbour-joining (NJ) tree (rectangular version) was calculated for the MLMT profiles of 120 strains of different species of the subgenus L. (Viannia), based on 15 microsatellite markers and using the Chord distance measure. The assignment of these strains to three main populations by the Bayesian modelbased clustering approach implemented in STRUCTURE is indicated by colored branches: population 1 (green), population 2 (red) and population 3 (blue). Strains belonging to these populations are listed in Table S1 Figure S3 NeighborNet network of Population 1 as inferred by STRUCTURE. Six strains were isolated from animal hosts, two from Choloepus didactylus and four from Didelphis marsupialis, and one from a sand fly, Lutzomyia anduzei, all other strains were isolated from human CL cases. (TIF) Figure S4 NeighborNet network of Population 2 as inferred by STRUCTURE. Four strains were isolated from animal hosts, two from dogs and one each from Nectomys sp. and Mesocricetus auratus, three from human MCL and three from human DL cases, all other strains were isolated from human CL cases. The assignment of the strains to the two sub-populations of Population 2, A and B, is indicated. Strains presenting mixed membership coefficients in two sub-populations are highlighted in grey. The two subpopulations are largely confirmed by the phylogenetic network albeit some strains occur at intermediate positions.

Supporting Information
(TIF) Figure S5 NeighborNet network of Population 3 as inferred by STRUCTURE. Eight strains were isolated from animal hosts, three from Dasypus sp., two from Cuniculus paca and one each from Coendou sp., Cebus apella and Rattus rattus, five from sand flies, three from Lutzomyia whitmani and one each from L. tuberculata, L.