Colonization and genetic diversification processes of Leishmania infantum in the Americas

Abstract

Leishmania infantum causes visceral leishmaniasis, a deadly vector-borne disease introduced to the Americas during the colonial era. This non-native trypanosomatid parasite has since established widespread transmission cycles using alternative vectors, and human infection has become a significant concern to public health, especially in Brazil. A multi-kilobase deletion was recently detected in Brazilian L. infantum genomes and is suggested to reduce susceptibility to the anti-leishmanial drug miltefosine. We show that deletion-carrying strains occur in at least 15 Brazilian states and describe diversity patterns suggesting that these derive from common ancestral mutants rather than from recurrent independent mutation events. We also show that the deleted locus and associated enzymatic activity is restored by hybridization with non-deletion type strains. Genetic exchange appears common in areas of secondary contact but also among closely related parasites. We examine demographic and ecological scenarios underlying this complex L. infantum population structure and discuss implications for disease control.

Fig. 1. Different read-depth profiles found in L. infantum isolates from Brazil.

Del isolates contain a >12 kb deletion between 1.122 Mb and 1.135 Mb on chr31 (e.g., Del_MT_3219 in the left graph). NonDel isolates do not contain the deletion, showing full read-depth at the locus (center graph). 8HTZ isolates are heterozygous for the deletion, with read-depth dropping to ca. 50% (right graph). Quantitative PCR confirmed heterozygosity at the deletion locus in monoclonal HTZ subcultures. MIX isolates appear to contain a mixture of NonDel and Del or HTZ profiles based on subclone PCR by Carnielli et al.. However, full read-depth is observed at the deletion locus in all MIX isolates, except in MIX_PI_05A and MIX_PI_08A (showing ca. 75% read-depth, see Supplementary Fig. 4). This suggests that NonDel cells are more abundant than Del and/or HTZ cells within MIX isolates. Circle radius indicates the number of isolates (each from a different canine or human host) representing the study site. Dotted circles represent study sites where multiple read-depth profiles occur (see table inset). Fill color indicates the majority read-depth profile at such study sites. The map was created in the open-source geographic information system Quantum GIS version 2.18.4 using Open Layers plugin access to Bing Aerial imagery. Microsoft product screen shot(s) reprinted with permission from Microsoft Corporation.

Fig. 2. Quantitative PCR confirms that intermediate read-depth profiles represent heterozygous deletions in L. infantum clones.

a HTZ_PI_2949 and HTZ_MT_3134 were selected as representatives of isolates for which read-depth drops to ca. 50% between 1.122 Mb and 1.135 Mb on chr31 (see Supplementary Fig. 4). DNA from monoclonal subcultures established from these two isolates was analyzed in qPCR targeting LinJ.31.2380 (within the chr31 deletion site) and LinJ.31.2330 (downstream of the chr31 deletion site). Differences in Ct values for LinJ.31.2330 between each HTZ sample and the NonDel reference (NonDel_MS_2666) were used to normalize a fold change estimate at LinJ.31.2380 based on the ∆∆Ct method by Livak and Schmittgen. Student’s t-test was applied to test whether fold change estimates obtained from n = 3 independent reactions differed significantly from the 1 : 1 ratio represented by the reference sample. Results were considered significant at *p < 0.05 and indicate that intermediate read-depth profiles represent abundant heterozygous deletions as opposed to mixtures of deletion-carrying and non-deletion-type cells within isolates. b Fold change was calculated the same way for monoclonal HTZ subcultures using the parental isolate as the reference. Results indicate that “unbalanced” heterozygotes also occur, e.g., subclone 2949 G1 appears to contain three chromosome copies with the chr31 deletion and one copy without.

Fig. 3. Homozygosity relative to Hardy–Weinberg expectations in New and Old World L. infantum isolates.

a The box plot shows median and interquartile ranges of genome-wide inbreeding coefficients (F_IS). Values are generally high for New World isolates. Values for HTZ isolates, however, all occur below the second quartile and strong excess heterozygosity is suggested in HTZ_MT_3134, HTZ_MT_3135, and HTZ_MT_3137. b Relatively low genome-wide F_IS in HTZ isolates is not driven by values from a subset of chromosomes. Values appear low throughout the genome. Circle fill color indicates New vs. Old World origin and read-depth profile on chr31.

Fig. 4. Phylogenetic relationships among New and Old World L. infantum isolates.

The maximum-likelihood tree was built using a general time-reversible substitution model with branch lengths corrected for ascertainment bias (i.e., the use of only nonvariant sites in sequence alignment). Pairwise genetic distances are haplotype-based, defined as the proportion of non-shared alleles across all SNP sites for which genotypes are called for all individuals (i.e., no missing data in alignment). Outlier isolates NonDel_MS_MAM, NonDel_FR_47, NonDel_PT_151, NonDel_PA_317, and NonDel_PA_85 are excluded. L. donovani strain MHOM/NP/03/BPK282/0 was temporarily included as an outgroup, to identify an L. infantum sample to subsequently root the tree. NonDel_ES_1345 became the outgroup. Circle fill color indicates New vs. Old World origin and read-depth profile on chr31. Font color specifies states sampled in Brazil. Isolates from other countries are labeled in black.

Fig. 5. Metric multidimensional scaling, simulated mating, and tree-to-graph conversion suggest admixture and hybridization between Del and NonDel L. infantum groups.

a Metric multidimensional scaling separates New and Old World (NW and OW) isolates on two axes of variation (goodness-of-fit = 0.40). NonDel isolates from Mato Grosso do Sul (MS) and Del isolates from Rio Grande do Norte (RN, see asterisk) and Mato Grosso (MT, see double-asterisk) position at opposite ends of axis 1, the primary axis of divergence within and between NW populations. HTZ isolates occur at intermediate positions (see pink circles) between these dissimilar groups. Other isolates with such intermediate positions are labeled and may also represent mating events between dissimilar groups. Gray, white, and cyan fill colors, respectively, indicate NonDel, Del, and MIX read-depth profiles found in the NW. Circles for OW (NonDel) isolates are green. Five outlier isolates are excluded as in Fig. 4. b Neighbor-joining positions of simulated hybrids (blue font, left tree) correspond to those of observed HTZ isolates (pink font, right tree) from MT. Hybrids were simulated in two steps. Random 50% haplotype contributions were first drawn from Del and NonDel isolates observed in MT and MS. The resultant offspring genotypes were then either let diversify through random mutation or subjected to a second round of Mendelian recombination as before. The same tree topology resulted in each of 100 simulation replicates. Trees are midpoint-rooted as opposed to outgroup-rooted as in Fig. 4. c Given that mating can create non-treelike divergence patterns within species, TreeMix was used to search iteratively for up to five migration edges that improve the fit of a maximum-likelihood tree built based on Gaussian approximation of genetic drift among isolates from MT, MS, RN, and OW groups. This input tree (black edges) suggests dichotomous differentiation into MT/RN and MS/OW clades and has a log-likelihood of 84.9206. Tree-to-graph conversion by addition of a migration edge from MS to MT increases log-likelihood to 84.9775. No other edges further increase the fit of the input tree. A four-population test also supports post-split admixture between MS and MT or RN, because differences in allele frequencies between MT and RN isolates correlate with those within the other population pair (F₄-statistic = 5 × 10⁻⁵, Z-score = 3.51).

Fig. 6. Ecto-3’-nucleotidase and ecto-ATPase activity correlates to read-depth profiles on chr31.

a Ecto-3′-nucleotidase activity was quantified by measuring the rate of inorganic phosphate (_Pi) release during adenosine 3’-AMP hydrolysis as described in Freitas-Mesquita et al.. Bar plots show mean and S.E. for at least three replicate assays (n = 3 independent experiments). Welch’s t-test was applied to test for statistical significance between pairs of samples at *p < 0.05. Results indicate a significant reduction of enzymatic activity in Del isolates relative to all NonDel and HTZ isolates. Activity appears significantly higher in NonDel_SC_3737 than in other NonDel and HTZ samples. Activity also appears significantly higher in NonDel_PI_2972 relative to HTZ_3134_B1 and HTZ_2949_B2 subcultures. b Ecto-ATPase activity was quantified with the same protocol except replacing 3’-AMP with equimolar ATP and Mg²⁺. T-tests between NonDel and Del isolates suggest higher ecto-ATPase activity in Del isolates than in all NonDel isolates, but larger samples sizes are required to substantiate the effect.