Experimental Selection of Paromomycin Resistance in Leishmania donovani Amastigotes Induces Variable Genomic Polymorphisms

The relatively high post-treatment relapse rates of paromomycin (PMM) in visceral leishmaniasis treatment and the swift emergence of experimental drug resistance challenge its broad application and urge for rational use and monitoring of resistance. However, no causal molecular mechanisms to Leishmania PMM resistance have been identified so far. To gain insights into potential resistance mechanisms, twelve experimentally selected Leishmania donovani clonal lines and the non-cloned preselection population, with variable degrees of PMM resistance, were subjected to whole genome sequencing. To identify genomic variations potentially associated with resistance, SNPs, Indels, chromosomal somy and gene copy number variations were compared between the different parasite lines. A total of 11 short nucleotide variations and the copy number alterations in 39 genes were correlated to PMM resistance. Some of the identified genes are involved in transcription, translation and protein turn-over (transcription elongation factor-like protein, RNA-binding protein, ribosomal protein L1a, 60S ribosomal protein L6, eukaryotic translation initiation factor 4E-1, proteasome regulatory non-ATP-ase subunit 3), virulence (major surface protease gp63, protein-tyrosine phosphatase 1-like protein), mitochondrial function (ADP/ATP mitochondrial carrier-like protein), signaling (phosphatidylinositol 3-related kinase, protein kinase putative and protein-tyrosine phosphatase 1-like protein) and vesicular trafficking (ras-related protein RAB1). These results indicate that, in Leishmania, the aminoglycoside PMM affects protein translational processes and underlines the complex and probably multifactorial origin of resistance.


Introduction
Leishmania parasites belong to the family of Trypanosomatidae and are unicellular flagellates that are transmitted between different vertebrates via the bites of infected sand flies. Depending on the species, Leishmania can cause either cutaneous (CL), mucocutaneous (MCL) or visceral leishmaniasis (VL). These parasites are known for their high plasticity and rapid adaptation to different environmental conditions, such as drug exposure [1].
The efficacy of the aminoglycoside antibiotic paromomycin (PMM) against both CL and VL has been demonstrated in monotherapy and in combination with other antileishmanials, with cure rates up to 80% [2][3][4][5][6]. The relatively high post-treatment relapse rates challenge the broad application of PMM in VL treatment, and urge for rational use and monitoring of the emergence of resistance [7], particularly because several in vitro and in vivo laboratory studies demonstrated fairly easy selection of resistance both on promastigotes and intracellular amastigotes [8][9][10][11]. A stepwise increase of PMM pressure  50 values of the pre-selection population and the 12 different post-selection clones against paromomycin (PMM) and amphotericin B (AmB). Results are expressed as the average IC 50 (in µM) ± the standard error of the mean (SEM). The clones that were considered resistant (resistance index (RI) > 5) are indicated in bold. Those considered susceptible (RI < 2) are indicated in italics. All reported IC 50 values are the result of at least two independent experiments run in duplicate.
Promastigotes from twelve clones and the non-cloned population in exponential growth phase were used for DNA extraction by the QIAamp DNA Mini Kit according to the manufacturer's instructions (Qiagen, Antwerp, Belgium). Genomic DNA was used to generate paired-end libraries using the NEBNext Ultra II FS DNA library prep kit for Illumina sequencing (New England Biolabs, Ipswich, MA, USA), according to the manufacturer's protocol, and amplified using NEBNext multiplex oligos for Illumina (dual index primers, New England Biolabs). Pooled and barcoded libraries were sequenced to 2 × 150 bp on an Illumina Hiseq 3000 platform. Reads were mapped to the L. donovani BPK282A1 genome version 39 available from TriTrypDB using BWA mem [13,14]. The genome coverage was estimated based on the mean read depth of genes with non-outlier  [15] and GFF coordinates. Detailed descriptions of each read library are summarized in Supplementary Table S1.
Short-nucleotide variant (SNV) distribution and allele frequency variation among L. donovani isolates were used to estimate populational diversity and potential association of mutations with PMM resistance. SNVs were identified with a focus on singlenucleotide polymorphism (SNPs) and short insertion-deletion mutations (INDELS). Calls for each sample were generated with FreeBayes version 1.2.0 (FreeBayes, version 1.2.0; Haplotype-based variant detection from short-read sequencing; Eric Garrison and Gabor Marth; arXiv:1207.3907; 2012) [16], calling only on positions where the depth was at least 5 reads. SNVs were filtered by: quality > 2000, minimum global allele frequency (MAF) > 0.05, number of different alleles = 1, mapping quality and base quality of the reference and alternate allele > 40, and proportion of observed reference and alternate alleles which were supported by properly paired read fragments > 0.9, using VCFtools version 0.  [18]. To identify SNVs that occurred during selection for PMM resistance, from the total of 7448 SNVs, we excluded the 5400 variants that presented MAF < 0.05 or low-quality calls, representing unreliable calls or variants that were present in all isolates, as these are simply differences between the isolates and the reference assembly. The sharing of SNV positions among the 13 isolates was estimated using BCFtools isec, and plotted in R using ggplot2. The populational genetic diversity (π) was estimated with VCFtools window-pi, in 50 kb windows. The images of the diversity values of these windows along the 36 L. donovani chromosomes were generated in R, using ggplot2 (https://cran.r-project.org/web/packages/ggplot2/index.html; accessed on 19 July 2021).
To estimate fluctuation in the 'within-clone allele frequency' for each SNV, the frequency of the alternate allele was estimated by dividing the read depth of the alternate allele by the total number of reads that mapped in the position. The heatmap representing the fluctuations in alternate allele frequency along the chromosomes was generated in R, with the heatmap2 function, from the gplots library (https://cran.r-project.org/web/packages/ gplots/index.html; accessed on 19 July 2021). L. donovani isolates were clustered based on the Manhattan distance of the allele frequency of all SNVs, by UPGMA (unweighted pair group method with arithmetic mean). SNV positions were drawn according to their chromosomal positions, or reordered based on clustering by UPGMA of their Manhattan distance. The genome-wide distribution density of the alternate allele frequencies of each L. donovani sequenced isolate was estimated in R, with the density function, and used to infer the parasite ploidy.
The potential impact of fluctuation in allele frequency of each SNV with PMM resistance was estimated by the Pearson correlation of the allele frequency and the isolate PMM IC 50 , in R, and the images were generated with ggplot2.
To evaluate the impact of gene and chromosomal copy alterations on PMM resistance, the copies of all L. donovani 8135 genes and 36 chromosomes were estimated for the 12 sequenced clonal lines and the non-cloned population, and compared with their PMM resistance profiles. To estimate gene copy numbers, the median read depth of each gene was calculated with BEDtools genomecov version 2.27.1 (BEDTools, version 2.27.1; a suite of utilities for comparing genomic features; Aaron Quinlan and Ira Hall; Bioinformatics: University of Oxford, Oxford, England 2010), using gene coordinates from L. donovani GFF version 39 from TriTrypDB. Read depths were transformed to relative copy number values (copies per haploid genome copies) by dividing the gene coverages by the genome coverage, estimated as the median coverage of non-outlier genes (Grubb's tests, with p < 0.05). The somy of each chromosome from each L. donovani clone was estimated based in the median coverage of non-outlier genes, as described by Reis-Cunha et al. [19]. Briefly, mapped reads were filtered by mapping quality 30 using SAMtools version 1.10 (SAMtools, version 1.10; Suite of programs for interacting with high-throughput sequencing data; Heng Li; Bioinformatics: University of Oxford, Oxford, England 2009) [20] and genes with outlier coverages were excluded, based on iterative Grubb's tests, with p < 0.05. The median read depth of all non-outlier genes in each chromosome was normalized by the genome coverage and considered as the chromosomal somy. The potential impact of each L. donovani gene copy number and chromosome somies on PMM resistance was estimated by Pearson correlation of the gene copy number or somy against clone PMM IC 50 values. All statistical analyses were carried out using R v3.6.2 (R, version 3.6.2; The R project for Statistical Computing; R core team; Vienna, Austria, 2019) [21], and graphical representations were generated with ggplot2 [22].

Amastigote Susceptibilty to Paromomycin and Amphotericin B
As described elsewhere [8], all in vitro selected clones showed lower susceptibility to PMM than the parent strain, in several cases achieving 5-10-fold IC 50 values (Table 1). To evaluate potential cross-resistance to other drugs, the intracellular drug susceptibility to amphotericin B was determined ( Table 1). The in vitro susceptibility of this naturally antimony-resistant strain against pentavalent and trivalent antimonials remained unaltered, as did the susceptibility against MIL [9]. PMM and amphotericin B susceptibility were not significantly correlated, which reinforces that the developed PMM resistance is specific to this drug, and not a result of generic drug resistance mechanisms.

SNPs and Indels
We detected 7448 single-nucleotide variants (SNV), where 2048 (1222 single-nucleotide polymorphisms (SNP) and 826 indel positions) vary among the evaluated Leishmania isolates (Supplementary Table S1). Overall, the majority of these positions were present in all post-selection clones (Figure 1a,b), which is expected given their similar origin. When looking at the genetic diversity (π) between the different clones and the pre-selection sample (Supplementary Figure S1), chromosomes 6, 12 and 34 (mostly euploid) and chromosome 31 (aneuploid) had relatively high π spikes, which were generally a result of co-localization of SNPs and indels ( Figure 1c). Even though the SNV positions were generally shared amongst most post-selection clones, variations in within-clone allele frequency patterns were observed ( Figure 2). These inconsistent counts of reads from the two alleles may be due to gene duplications/losses. Therefore, SNVs in these clones were further compared at the allele frequency level in the downstream analyses. It is noteworthy that the genome-wide within-clone allele frequency patterns suggest that the progenitor isolate was triploid, as are all post-selection clones. We expect that the variation in withinclone allele frequencies is due to loss of heterozygosity due to fluctuating ploidy [23,24] and/or the difficulty in estimating true allele counts in a triploid sample.
To identify SNVs that associate with PMM resistance, we have compared withinclone allele frequencies to PMM IC 50 values. A total of 11 SNV positions were correlated with PMM resistance, from which 9 were located in intergenic regions and 2 resulted in synonymous mutations ( Figure 3, Supplementary Figures S2 and S3). These SNVs were not clustered in close chromosomal regions and were dispersed over different regions and chromosomes. The intergenic SNPs on chromosomes 18 and 35 showed the strongest correlations with resistance and are located upstream of the chaperone DnaJ protein (LdBPK_181410.1) and aspartate aminotransferase genes (LdBPK_350840.1). The 2 SNPs resulting in synonymous mutations were found in genes LdBPK_020680.1 (ATPdependent Clp protease subunit, heat-shock protein 78 (HSP78)) and LdBPK_362510.1 (sterol 24-c-methyltransferase) (Supplementary Table S2). For both genes, a negative correlation was found between the alternate allele frequency and PMM resistance, suggesting that loss Microorganisms 2021, 9, 1546 5 of 14 of the alternate alleles towards a reference-homozygous genotype correlates with PMM resistance. The synonymous SNP in LdBPK_020680.1 (UCG → UCC) results in a lower relative synonymous codon usage of the alternate codon (RSCU decrease from 1.439 to 1.115). For the sterol 24-c-methyltransferase, the synonymous mutation results in a higher frequency codon (UUU → UUC) with a substantial RSCU increase from 0.72 to 1.28 [25]. not clustered in close chromosomal regions and were dispersed over different regions and chromosomes. The intergenic SNPs on chromosomes 18 and 35 showed the strongest correlations with resistance and are located upstream of the chaperone DnaJ protein (LdBPK_181410.1) and aspartate aminotransferase genes (LdBPK_350840.1). The 2 SNPs resulting in synonymous mutations were found in genes LdBPK_020680.1 (ATP-dependent Clp protease subunit, heat-shock protein 78 (HSP78)) and LdBPK_362510.1 (sterol 24-c-methyltransferase) (Supplementary Table S2). For both genes, a negative correlation was found between the alternate allele frequency and PMM resistance, suggesting that loss of the alternate alleles towards a reference-homozygous genotype correlates with PMM resistance. The synonymous SNP in LdBPK_020680.1 (UCG → UCC) results in a lower relative synonymous codon usage of the alternate codon (RSCU decrease from 1.439 to 1.115). For the sterol 24-c-methyltransferase, the synonymous mutation results in a higher frequency codon (UUU → UUC) with a substantial RSCU increase from 0.72 to 1.28 [25].

Gene Copy Number Variations (CNV)
A significant correlation could be noted between PMM resistance and an increased or decreased gene copy number of 33 and 6 genes, respectively ( Figure 4, Supplementary  Figures S3 and S6). The genes with the highest correlation to PMM resistance are summarized in Table 2. The gene with the largest difference between the maximum and minimum gene copy values is a small hypothetical protein (LdBPK_020370.1), with no identifiable functional motifs. Although 14 genes are hypothetical genes, the other 25 are putatively enrolled in several functions, such as: transcription, translation and protein turn-over (transcription elongation factor-like protein, RNA-binding protein, ribosomal protein L1a, 60S ribosomal protein L6, eukaryotic translation initiation factor 4E-1, proteasome regulatory non-ATP-ase subunit 3 and RING-H2 zinc finger, with possible involvement in protein ubiquitination), virulence (major surface protease gp63, protein-tyrosine phosphatase 1-like protein), mitochondrial function (ADP/ATP mitochondrial carrier-like protein), signaling (phosphatidylinositol 3-related kinase, protein kinase putative and protein-tyrosine phosphatase 1-like protein) and vesicular trafficking (ras-related protein RAB1). Pearson correlation plots of the 11 significant SNVs with strong correlation to PMM resistance. The box title contains the SNV chromosome, position and correlation, r. The color of the dots reflects the level of resistance of the individual clones to PMM, where yellow, grey and purple correspond respectively to RI < 2, 2 < RI < 5 and RI > 5.

Chromosomal Copy Number Variation (CCNV)
Of the 36 L. donovani chromosomes, two had significant correlations between chromosome expansion and PMM resistance, Chr33 (r = 0.64, p-value = 0.02) and Chr34 (r = 0.56, p-value = 0.04) (Supplementary Figure S4). However, the correlation, r, was lower than 0.7 and the amplitude of variation in these chromosome copies was low (Chr33, minimum = 1.29, maximum = 1.39; Chr34, minimum = 0.89, maximum = 0.94). Within their complex pattern of CCNV, the "chromosomal copy number drift" that could be observed for all of the clones (Supplementary Figure S5) might be potentially explained by adaptation to in vitro culture [26,27].

Discussion
Although PMM resistance mechanisms are well-established for bacteria [28,29], several hypotheses have been postulated for Leishmania, such as interference with ribosomal protein synthesis and inhibition of respiration [11,[30][31][32][33][34][35]. However, no genetic resistance markers have been identified as yet. To try to identify these markers, the present study performed whole genome sequencing and comparison of 12 amastigote-selected clones and the pre-selection population, that had a large range of PMM resistances (Table 1). Consistent with previous analyses, we were not able to identify a specific genetic change that unequivocally explains PMM resistance. While the 12 post-selection clones contained numerous new SNV mutations that were fixed in all post-selection clones, these common mutations cannot explain the large variation in PMM resistance within these clones. As clone IC 50 values range from 57 µM (only marginally greater than the pre-selection strain at 45 µM) to 417 µM, we hypothesized that multiple other non-fixed genetic differences contribute to PMM resistance. To determine which of these genetic differences may be associated with resistance, we compared within-clone allele frequencies, gene CNV and chromosome CNV with PMM resistance, using a similar approach to the genome-wide association study (GWAS) approach.
The majority of the SNV positions whose within-clone allele frequency correlated with PMM resistance were observed in intergenic regions, including notably two upstream of the chaperone DnaJ protein and aspartate aminotransferase genes. Aspartate aminotransferase has been described to be significantly elevated in patients receiving PMM, indicating hepatotoxicity or another drug effect [36]. Further, only two synonymous SNPs were identified within genes, in HSP78 and sterol 24-c-methyltransferase. Although synonymous mutations are often not phenotypically relevant, they may affect translation by differences in codon usage. Both affected genes are of potential relevance to PMM resistance, as the HSP78 gene in L. donovani is an important virulence factor, needed to suppress immune activation and escape from NO-mediated toxicity in macrophages [37], whereas sterol 24-cmethyltransferase in L. major has been linked to mitochondrial function and virulence [38]. For both genes, a negative correlation was found between the alternate allele frequency and PMM resistance. The lower relative synonymous codon usage of the alternate codon in the HSP78 gene (RSCU decrease from 1.439 to 1.115 for UCG → UCC) would indicate a higher expression to be linked to PMM resistance. A negative correlation was also found for sterol 24-c-methyltransferase alternate allele frequency and the PMM IC 50 values. The reduced frequency of the alternate codon (codon UUC → UUU, RSCU decrease from 1.28 to 0.72) is anticipated to decrease the expression levels in favor of increased resistance. Genomic instability of sterol 24-c-methyltransferase has already been associated with amphotericin B resistance in L. mexicana, L. major, L. donovani and L. infantum [39,40]. In this study, however, no correlation could be demonstrated between the PMM resistance profiles of the different clones and their susceptibility towards amphotericin B. In a study with promastigoteselected PMM-resistant lines, no differential expression of sterol 24-c-methyltransferase was noted [40], but it should be mentioned that in these strains, PMM susceptibility of the amastigote stage remained unaltered [12,41], whereas selection on the amastigote level triggers a completely different outcome [9].
Despite stability of the PMM resistance upon parasite passage in mice or sand flies [9,42], suggesting a stable genomic modification, this comparison mostly identified gene CNVs correlating with resistance. Chromosomal somy changes could not be unequivocally correlated to the variable PMM susceptibility of the different clones. Due to the extremely high genome plasticity of Leishmania, somy changes should be cautiously interpreted, as they are readily affected by culture conditions, prolonged cultivation and the associated virulence changes [43]. Nonetheless, the CNV of 39 genes was found to be statistically associated with PMM resistance, where 33 were a result of gene expansions and 6 of gene losses. While 14 of these code for hypothetical proteins, the other 25 are known to be involved in transcription, translation and protein turn-over, virulence, mitochondrial function, signaling and vesicular trafficking. Although the exact mechanisms of action of PMM in Leishmania are yet to be elucidated, past research on experimentally selected resistant L. donovani already suggested that PMM not only impairs protein synthesis and modifies membrane fluidity and permeability, but also causes respiratory dysfunction with associated alterations of the mitochondrial membrane potential [11,32,35,40,[44][45][46].
When only considering those genes with a copy number difference amplitude of at least 0.3 haploid genome copies, 11 genes could be identified, of which 6 have been annotated (i.e., ADP/ATP mitochondrial carrier-like protein, ras-related protein RAB1, 60S ribosomal protein L6, proteasome regulatory non-ATP-ase subunit 3, an RNA-binding protein and ribosomal protein L1a). ADP/ATP mitochondrial carrier-like proteins are known to transfer molecules across mitochondrial membranes [47] and could therefore potentially play a role in drug transport into a main target organelle. The Ras-related protein RAB1 is involved in vesicular trafficking and is a known regulator of secretory pathways of several important Leishmania virulence factors, including gp63 and acid phosphatase [48]. Proteasome regulatory non-ATP-ase subunit 3 is important in the regulation of proteosomal functions, of which inhibition hampers specific stages of morphological differentiation in various protozoan parasites, including Leishmania [49]. In L. infantum, proteasomal functions were shown to be essential for replication and intracellular survival of amastigotes in the host cell [50]. Its upregulation upon PMM selection thus supports the increased intracellular parasite fitness observed by our group and others [12,43,51]. The increased virulence may also relate to the detected CNVs for the virulence genes major surface protease gp63 and protein-tyrosine phosphatase 1-like protein [52]. The identification of the ribosomal protein L6 gene is also noteworthy as it has been shown to be involved in bacterial resistance against aminoglycosides in the past [53]. Moreover, its close homologue, the 60s ribosomal protein L23, is known to be overexpressed in antimony-resistant parasites [54]. Additionally, RNA-binding protein overexpression has already been described in antimony-resistant strains [55], and could be linked to an increased parasite replication efficacy [56]. A negative correlation was found for the ribosomal protein L1a and quiescin sulfhydryl oxidase. Although not much information can be found in the literature on the role of both gene products in antileishmanial drug resistance, quiescin sulfhydryl oxidase has been associated with the induction of reduced cell proliferation and quiescence, for example in fibroblasts [57,58]. Its decreased expression may increase cell proliferation, which in theory fits the increased parasite fitness observed for PMM-resistant parasites [12]. Finally, the presence of several hypothetical proteins in which gene copy number alteration was potentially associated with PMM resistance reinforces the importance of characterizing unknown protein functions in Leishmania.
Other research groups have evaluated PMM resistance markers using genomic, metabolomic or transcriptomic comparisons [43,44,59]. Metabolic and lipidomic changes were shown to be strain-specific [43], and up until now, no suitable PMM resistance biomarkers have been identified. None of the post-selection variants that we found to be associated with PMM resistance were previously identified by earlier research as associated with Leishmania's PMM resistance. However, this is not completely unexpected as the present study was the only one to have selected PMM resistance directly on intracellular amastigotes. Our results indicate the likely multifactorial nature of PMM resistance, which complicates the identification of a single genetic marker for application in clinical settings. In-depth research and validation of the identified genetic variations in a range of clinical isolates will be necessary to understand their biological importance and use as potential resistance markers.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.3 390/microorganisms9081546/s1, Figure S1: Genetic diversity (π) in the 36 L donovani chromosomes, Figure S2: Correlation between SNVs and PMM resistance in all L donovani chromosomes, Figure S3: QQ plot of the genome variations association with phenotypes, Figure S4: Correlation between CCNVs from Chr33 and Chr34 and PMM resistance, Figure S5: CCNV among the L donovani isolates, Figure S6: Correlation between Gene CNVs and PMM resistance in all L donovani chromosomes, Table S1: WGS read libraries and SNV counts in the 12 clones and the polyclonal original population (Pre S57), Table S2: SNV correlated to PMM resistance.