Insights into Genetic Heterogeneity and Drug Resistance in Leishmania donovani of the Indian sub-continent from genomic data

doi:10.21203/rs.3.rs-7267613/v1

Download PDF

Research Article

Insights into Genetic Heterogeneity and Drug Resistance in Leishmania donovani of the Indian sub-continent from genomic data

https://doi.org/10.21203/rs.3.rs-7267613/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 19 Nov, 2025

Read the published version in Antonie van Leeuwenhoek →

You are reading this latest preprint version

Leishmaniasis is one of the neglected tropical diseases that is endemic to over 90 countries and it’s cases are being reported from non-endemic countries as well like Austria. The cases Visceral Leishmaniasis (VL) caused by Leishmania donovani are concentrated on the Indian Subcontinent. Several studies on the genetic heterogeneity, aneuploidy, and drug resistance emergence reported from Indian subcontinent as well as globally. However, no research has yet inspected the genomic data from the Indian subcontinent to depict a larger microscopic investigation at the genome level. We have considered whole genome sequence (WGS) data from the publicly available database i.e. ENA. Analysis has shown that there is tetraploidy in the chromosome 31, trisomy in chromosome 2 and 8 and trisomy was observed among chromosome 6 and 15 among some samples. This aneuploidy pattern is evolving over time observed in the present study. The pattern of aneuploidy variations of the Indian subcontinent sample differ from from other continents. Further, ATP-binding cassette (ABC) family, Amastin-like surface proteins, A2 genes, amino acid permeases, heat shock 70-related protein 1, mitochondrial precursor and sodium stibogluconate resistance protein were the proteins that exhibited maximum number of mutations among all the analyzed samples. The proteins showing highest number of mutations belong to membrane prtoeins that are involved in drug resistance mechanism. Most of these proteins are involved in the virulence and drug resistant mechanism. The present study provides the possible candidates which can be targeted to disarm the virulence of the protozoan and drug candidates for therapeutic interventions.

Whole Genome Sequences

Leishmania

Indian sub-continent

Drug Resistance

Single Nucleotide Polymorphism

Mosaic Aneuploidy

Leishmaniasis is caused by vector-borne, obligate intracellular, kinetoplastid parasites belonging to the genus Leishmania (Family Trypanosomatidae) (Torres-Guerrero, Quintanilla-Cedillo et al. 2017). Transmission is through the bite of an infected female phlebotomine sandfly(Schlein 1993). Leishmaniasis has three main types of clinical manifestations: Cutaneous Leishmaniasis, Mucosal Leishmaniasis, and Visceral Leishmaniasis. Disease manifestation and severity depends upon the host immune response, the causal pathogen and/or any comorbidities(Alvarenga, Escalda et al. 2010) (Awasthi, Mathur et al. 2004). The World Health Organization (WHO) has declared Leishmaniasis as a Neglected Tropical Disease (NTD) (Kumari, Lakhanpal et al. 2023). It is most prevalent in the subtropical and tropical regions of the globe; mainly Europe, Northern Africa, the Middle East, Asia and part of South America (Martínez, Guerrero et al. 2011). According to the CDC, current estimates of disease incidence for Cutaneous Leishmaniasis have ranged from 700,000 to 1.2 million cases per year and it's been estimated that the number of new cases per year for Visceral Leishmaniasis has decreased from 400,000 to less than 100,000 (Prevention, 2020).

Visceral Leishmaniasis (VL), caused by the Leishmania donovani complex which includes L. donovani and L. infantum, is the most lethal form of Leishmaniasis. Majority of disease burden (90%) has been reported in the following countries; India, Eiteria, Iraq, Somalia, Brazil, Ethiopia, Kenya, South Sudan, Nepal, China, and Sudan (Patil and Chatterjee 2023, Saha, Borah et al. 2023). India contributes 18% of the total reported VL cases. VL is endemic to 54 districts scattered in four states — Bihar (33 out of 38 districts), Jharkhand (4 out of 24 districts), Uttar Pradesh (6 out of 75 districts) and West Bengal (11 out of 23 districts). In recent time, the vector of leishmaniasis has been detected in the Western Himalayan region of India like Jammu and Kashmir, Gujarat, Assam, Himachal Pradesh, Haryana, Tamil Nadu, Puducherry, Madhya Pradesh, Kerala, Sikkim and Uttaranchal which raises a concern (Patil and Chatterjee 2023, Sofi, Godara et al. 2023). With the unconventional spread of Leishmaniasis to Himalaya, there is a need for the reassessment of Leishmaniasis endemicity in the Himalayan region of North India.

Pentamidine, Antimonial [Sodium stibogluconate (SSG) and meglumine antimoniate (MA)], Miltefosine, Amphotericin B and Paromomycin are the anti-leishmaniasis agents used in the first line and second line therapy (Bharadava et al., 2024; (Mann, Frasca et al. 2021)). Beauvericin (BEA) is a fungal secondary metabolite synthesized by Beauveria bassiana which is potent ABC transporter inhibitor. Its application along with anti-leishmanial drug may enhance therapeutic index of anti-leishmanial agents (Al Khoury et al., 2024). Emergence of drug-resistance strains has rendered traditionally employed treatment methods ineffective. The use of chemotherapy remains viable but it is expensive and has toxic side effects (Lira, Sundar et al. 1999). Hence there is an urgent need for novel drugs to combat the emerging peril of drug-resistant strains and effectively decrease the virulent effect of Leishmania (Bora and Jha 2020). Other mode of therapy includes thermotherapy, electric field stimulation and diverse topical treatment methods that are in clinical or preclinical trials, such as formulations containing paromomycin, imiquimod, miltefosine, morphine, and buparvaquone components (Pinheiro

et al., 2022). Currently new chemical classes which are in clinical phase 1 like anti-leishmanial classes nitroimidazoles and benzoxaboroles, with DNDI-0690 and DNDI-6148 as lead compounds, respectively, the pyrazolopyrimidine scaffold recognized as CRK-12 kinase inhibitors.

Advancements in high throughput sequencing technologies have greatly reduced the time and capital required to sequence whole genomes of organisms, making it possible to analyse complete genome sequences and investigate the biology in greater detail (Reuter, Spacek et al. 2015). WGS has enabled the genetic characterization of novel spp. of Leishmania ellisi from USA (Sapp, Low et al. 2024). The meta-analysis of transcriptomic data of Leishmania spp. has been carried out to understand the immune system’s response to its infection (Rezaei, Tahmasebi et al. 2024). Recent changes have been observed in the epidemiology of leishmaniasis with reports of it emerging in the non-endemic country such as Austria. The development highlights the impact of climate and habitat changes, urbanization, and globalization on the spread of the disease in regions previously unaffected (Bora and Nath Jha 2019, Riebenbauer, Czerny et al. 2024). There are reports of VL from certain villages of Nepal, where the patients did not have any travel history outside of their village. However, it is likely that traders and visitors might have contributed to the spread of the disease. The reported incidences of infection in new loci of Nepal can be a threat to the elimination effort of Leishmaniasis (Joshi, Banjara et al. 2024). L. donovani offspring obtained from a cross between two sand fly species has resulted in genetic diversity in HASPB antigen which has major impact on the therapeutic vaccine design (Kumari, Mamud et al. 2023, Sádlová, Yeo et al. 2024). The genomic studies of Leishmania spp. has explored the genetic plasticity and constantly ongoing rearrangement of the genome that leads to variation in copy number, gene clusters or the chromosome (Santi and Murta 2022). Several studies indicate that the clinical isolates of Leishmania spp. has shown increased resistance against pentavalent antimony (Santi and Murta 2022). In Indian samples, it has been observed that high levels of As(III) in the water has led to the selection of Sb-resistant parasites. Further, change in aneuploidy has been triggered by the drug resistance against antimony (Dumetz, Imamura et al. 2017). The changing epidemiology of leishmaniasis globally and in India is a matter of concern. To understand the genetic constitution of Leishmania spp., we have retrieved publicly available Leishmania donovani whole genome sequencing (WGS) data specifically isolated from the Indian subcontinent. We analyzed the genomes of the various strains and reported the variants which can be targeted for the drug development, aiming to reduce the severity of the disease. We tracked and reported the variants located within the functional domains of proteins. The study will help to explore the genome of Leishmania under different stress, environmental conditions and drug pressure.

Retrieval of WGS data of Leishmania donovani and alignment to reference genome

A total of 34 sample data of the Indian Subcontinent were downloaded through a combination of text-based and advanced search approach from the European Nucleotide Archive (https://www.ebi.ac.uk/ena/browser/) (Supplementary Table 1). The samples sequencing approach was whole genome sequencing and biological source of the samples was late passage, bone marrow and Atypical cutaneous Leishmania donovani. Fastq files were downloaded and quality check was performed using the FASTQC version: 0.11.9 (Andrews 2010). FASTQC generates a visual report of the sequencing files based on several analytical modules. These reports were used to select specific cut-offs for riding sub-quality sequence calls. Base calls with a Q score of less than 20 were filtered out using Trimmomatic Version: 0.39 (Bolger, Lohse et al. 2014). Leishmania donovani reference genome assembly ASM22713v2 (Downing, Imamura et al. 2011) was utilised and the GenBank (Clark, Karsch-Mizrachi et al. 2016) assembly GCA_000227135.2 was retrieved from NCBI. The reference genome was indexed using Burrows-Wheeler aligner (bwa) version: 0.7.17-r1188 (https://bio-bwa.sourceforge.net/bwa.shtml) ‘index’ tool while SAMtools Version: 1.10 (using htslib 1.10.2-3ubuntu0.1) (Danecek, Bonfield et al. 2021), PICARD version:2.27.5 (https://github.com/broadinstitute/picard) were used to create a dictionary of the reference sequence. The ‘faidx’ (Shirley, Ma et al. 2015) tool was used to create a bed file. Reads were aligned against the reference sequence using bwa mem that resulted in a sequence alignment map (SAM) files (Li, Handsaker et al. 2009). SAM files were sorted, duplicate marked, converted into binary alignment map (BAM) files, and read group information added using PICARD. SAMtools were used to get alignment, depth and coverage information, which were employed to discard low depth and low coverage samples from further analysis.

Variant discovery and annotation

Genome Analysis Toolkit (GATK v4.3.0.0) (Van der Auwera, Carneiro et al. 2013), a suite of tools was used to call variants following the best practices workflow for germline short variants – Single Nucleotide Polymorphisms (SNPs) and Insertion/Deletions (InDels) – discovery. Variant calling was performed using GATK HaplotypeCaller using default settings. The resultant variant call format (VCFs) were filtered for SNPs, removing Insertion/Deletions. The following thresholds were chosen to filter SNPs following the GATK hard filtering recommendations, QD < 20, QUAL < 30, SOR < 3, FS < 60, MQ < 4, MQRankSum < − 12.5 and ReadPosRankSum < − 8, to remove low quality or false positive variants which may contaminate downstream processing, annotation and analysis. Annotation was performed using ANNOVAR (Wang, Li et al. 2010) and snpEff (version:5.1d) (Cingolani, Platts et al. 2012) separately to obtain common annotated genes.

Variant Analysis

Multianno files obtained after annotating VCF files using ANNOVAR were imported into R and analyzed using the R package MAFtools from Bioconductor (Mayakonda, Lin et al. 2018). VCF files annotated with snpEff were converted to Mutation Annotation Format (MAF) using VCF2maf perl script (Mskcc/VCF2maf, n.d.) to extract variant information such as gene id, variant, variant location, variant type and effect, HGVS.p. SNPs and InDels were split into different maf files due to quality control and processing differences. Using the MAFtools R package, a summary of the MAF files was plotted and most mutated genes were discovered and reported using the plotmafsummary and oncoplot functions of the package. MAF files containing InDel information were analyzed using shell commands and various information associated with InDels were extracted. InDel affected genes common to all strains were investigated; top mutated genes with more than 15 InDel events were compared and projected with a heatmap using the pheatmap package.

Aneuploidy and copy number analysis

Aneuploidy and Somy profiles were calculated using a read depth based approach. Samtools was used to extract raw read depth information of the aligned bases from the BAM files. Normalised read depths were estimated by dividing raw depth by the mean raw read depth of specific chromosomes. The haploid depth was defined as the median of the normalised read depths and chromosomal somy was estimated by multiplying the haploid depth by two given the assumption that there exist two copies of chromosomes. Sample ploidy was defined as the median of the chromosomal somies calculated in the previous step. A heatmap was prepared to visualise the different chromosomal profiles using the pheatmap package in R. Copy Number Variations (CNV) were investigated with Control-FREEC (Boeva, Popova et al. 2012) using a sliding window approach, window size of 5000 and step size of 250. The following parameters were also included in the config file: breakPointThreshold = 0.8, minExpectedGC = 0.35, maxExpectedGC = 0.65, ploidy = 2, and inputFormat = BAM.

Chromosomal Somy and Copy Number Variant estimation

Chromosomal Somy was estimated using read depth statistics of the sequenced reads. Raw depth of a position (i _{Raw depth)} and mean read depth (Ch _{Mean depth)} of a chromosome were used to calculate the normalised read depth of a position (i_Norm) and the chromosomal Somy (S) was calculated using the following formula: 2 × Median( i _Norm = i _{Raw depth} / Ch _{Mean depth} ). The range of monosomy, disomy, trisomy, tetrasomy, and pentasomy was then used to define the full cell-normalised chromosome depth or S as S < 1.5, 1.5 ≤ S < 2.5, 2.5 ≤ S < 3.5, 3.5 ≤ S < 4.5, and 4.5 ≤ S < 5.5, respectively, as previously described. (Dumetz F, Imamura H, Sanders M, Seblova V, Myskova J, Pescher P, Vanaerschot M, Meehan CJ, Cuypers B, De Muylder G, Späth GF, Bussotti G, Vermeesch JR, Berriman M, Cotton JA, Volf P, Dujardin JC, Domagalska MA. Modulation of Aneuploidy in Leishmania donovani during Adaptation to Different In Vitro and In Vivo Environments and Its Impact on Gene Expression. mBio. 2017 May 23;8(3):e00599-17. doi: 10.1128/mBio.00599-17. PMID: 28536289; PMCID: PMC5442457.

)The results were projected onto a heatmap and dendrograms using the R package ‘pheatmap’. The isolates were clustered into groups based on the euclidean distances between their chromosomal profiles. Chromosomes were also clustered into groups based on the Euclidean distance of chromosomal copy numbers profile of a chromosome across isolates.

CONTROL-FREEC (Boeva, Popova et al. 2012) was employed to investigate the Copy Number Variation (CNV). Control-FreeC parameters were defined as follows:

Window size = 5,000, Step size = 250, breakPointThreshold = 0.8, minExpectedGC = 0.35, maxExpectedGC = 0.65, ploidy = 2, and inputFormat = BAM.

Functional domain analysis

For functional domain analysis, genes affected by missense mutations which were common to most samples were analyzed. Multiple missense variants were observed in most of the genes but variants which affected the majority of the samples were selected for downstream analysis. 70 genes were annotated using the Uniprot database (The UniProt Consortium, 2017); 19 genes were uncharacterized/ hypothetical and were removed from subsequent analysis. Protein sequence information for 51 remaining genes was retrieved from Uniprot, which included protein name, length, domain(s) and domain location(s). The domain annotations were provided by Interpro (Paysan-Lafosse, Blum et al. 2023). Only 30 protein sequences were found to have relevant Interpro domain annotations; variants found inside these domains were reported.

Sample details and quality check

The organism was isolated from human, biological replicate and laboratory-adapted strains. Leishmania sp. cultured under different conditions and drug resistance has shown different mutations and genetic makeup (Supplementary Table 1). The samples were quality checked, processed, and aligned to the reference genome. The Binary Alignment Map (BAM) files were used to assess the alignment and depth statistics of each sample (Fig. 1). 2 samples were discarded from further analysis because of insufficient alignment to the reference genome, due to the incompatibility of the read files with BWA alignment software. The 31 remaining samples were used for further analysis. The samples used in this study have been provided in Supplementary Table 2.

Aneuploidy change pattern in Leishmania donovani samples from Indian Sub-continent

All the analyzed samples have shown tetrasomic state as reported in the earlier studies. In addition to this, trisomy was observed in chromosomes 2 and 8 in BHU 1220, BHU 575, MHOM/IN/1983/AG83, MHOM/IN/2013, MHOM/IN/2014, K133 samples and tetrasomic states in BPK275/0/cl18 samples. Contradictorily, the Nepalese strain BPK275/0/cl18 did not show either trisomy or tetrasomy in chromosomes 2 and 8, and clustered away from the other Nepalese strains, with the K133 strain isolated from India. A few Indian strains along with one Nepalese strain, BPK275/0/cl18, K133PMM/R, K133AS/R, K133WT, MHOM/IN/2010/T3, MHOM/IN/IN/T8, MHOM/IN/2012/T9, MHOM/IN/2014/A2, MHOM/IN/2013/P2, MHOM/IN/2014/H1, BHU 575, showed the disomic state of chromosome 2 and clustered together separate from the other Indian cluster (Fig. 2). It should be noted that the strains are clustering into 4 major clusters on the similarity in their aneuploidy profile. These clusters encompass the two distinctly separate groups of Indian strains, a Nepalese cluster and a mixed group with Nepalese and Indian strains. Here, we observed that the Indian cluster of strains isolated from India BHU 1220, MHOM/IN/1983/AG83.i and ii, West Bengal MHOM/IN/2013/M1, Jharkhand MHOM/IN/2014/P6, and MHOM/IN/2014/T10, showed expected karyotypes as most chromosomes were disomic barring the trisomic chromosomes 6 and 15, and the tetrasomic chromosome 31. On the other hand, the second Indian cluster comprising strains isolated from Bihar MHOM/IN/2010/T3, MHOM/IN/2010/T7, MHOM/IN/2012/T9, MHOM/IN/2013/P2, West Bengal MHOM/IN/2012/T8, Jharkhand MHOM/IN/2014/H1, MHOM/IN/2014/A2 and BHU 575 strain showed a remarkable difference in aneuploidy profile; chromosome 31 and 8 showed an increase in copy number which reflected in their tetrasomic and trisomic states, chromosomes 2, 5, 9, 11, 13, 14, 16, 20, 22, 23, 26, 35 showed variable aneuploidy, ranging from disomy to trisomy, which was lacking in the previous Indian cluster. The Nepalese strains grouped together and surprisingly were closer to an Indian cluster due to similarity in aneuploidy profiles, except the chromosome 2 which distinctly showed tetrasomic condition in Nepalese strains BPK275/0/cl18. Finally, the remaining cluster, was a mixed group of one Nepalese strain and three Indian strains. This may have been due to the similarity in tetrasomic state of chromosome 5 along with trisomy observed in chromosomes 6, 9, 13, 15, 20, 23, 26, and 33. Chromosomes 2, 8 and 31 were also observed to be clustered together, possibly due to their peculiar nature of existing in unusually higher copy number states, mostly trisomic and/or tetrasomic with respect to the other chromosomes. The chromosomes 5, 9, 13, 14, 16, 20, 23, 26, 33, and 35 were grouped together, possibly because of their tendency to exhibit variable aneuploidies ranging from disomy to tetrasomy. Chromosomes 6 and 15 show occasional trisomy and hence grouped together. Rest of the chromosomes showed normal copy numbers with minimal deviation from disomy. The overall heatmap of the aneuploidy change in the samples from Indian subcontinent has shown trisomy or tetrasomy in the whole or partially in chromosome 2 and 31. This observation is supported by earlier study that has shown chromosome 31 is superanumerary and tetrasomic (Pilling, Reis-Cunha et al. 2023). The aneuploidy changes of Leishmania under the effect of paromomycin can be observed in the samples with BPK* id, most of the samples have shown trisomy in Chromosome 2. WGS samples of late passage has shown tetrasomy in chromosome 31 and trisomy in chromosome 5 and 6. The present analysis show a very high genetic modification in Leishmania donovani sequenced from patients of Indian sub-continent and drug resistance also alter the aneuploidy of Leishmania.

Distribution of variants among samples

The variant analysis has shown that the sample from the Himachal region of India has recorded the highest number of variants in the coding sequences, gene interval and chromosomal interval (Supplementary Figs. 1 and 2). Rest of the samples have shown variants in the same range i.e. 1000 to 10000 in number. We observed the genotype distribution of SNPs in different strains; Indian strains MHOM/IN/2010, MHOM/IN/2012, MHOM/IN/2013, and MHOM/IN/2014 showed an average SNP frequency of 1900. These strains were isolated from different regions of India mainly from Bihar, Delhi, Jharkhand, and West Bengal. Strains with unknown states as geographical locations were tagged as ‘Indian’. The strain MHOM/IN/2010/T7, isolated from Bihar, had the highest SNP occurrence of 1975 while Jharkhand isolate MHOM/IN/2014/H1 had the lowest SNP occurrence of 1868. Indian strains BHU1220, isolated from Bihar, showed 2089 SNPs, MHOM/IN/1983/AG showed a higher SNP frequency of 2707 and 2719, and BHU575/0/cl2 with 3800 SNPs was the highest. Indian strains K133WT, K133AS-R, and K133PM-R showed abundance with each having 3713, 3281 and 3488 SNPs respectively. We looked at the distribution of genotypes in these samples; 86–91% of SNPs had Heterozygous Reference Alternate (Het RA) genotype, 9–13% had Homozygous Alternate (Hom AA) genotype and a miniscule percentage, around 0–0.002% of SNPs represented the Heterozygous Alternate (Het AA) genotype, which represents two different alternate alleles in place of the reference (Fig. 3).

The Nepalese strains BPK275/0/cl18 showed a distinctly higher SNP abundance than Indian strains and an average SNP abundance of 3482 SNPs. BPK275/0/cl18.vii, isolated from Koshi province, showed the lowest SNP abundance of 2962 while BPK275/0/cl18 was the highest – 4003, among all strains from both India and Nepal. Most of the Nepalese strains were isolated from the Koshi province of Nepal, barring BPK275/0/cl18 of the study SRR15851294 which did not have an associated geographical location. The variant frequency distribution is similar among all the samples except the one obtained from the Himalayan region.

We also analyzed the Insertion/Deletions occurrence frequency in the samples and reported the different genotypes associated (Fig. 4). Quite interestingly, the Indian strains MHOM/IN/2010 T3, T7, MHOM/IN/2012/T8, T9, MHOM/IN/2013/M1, P2, MHOM/IN/2014/A2, H1, P6 and T10 had remarkably high insertion/deletion frequency compared to rest of the samples. The highest occurrence was in P2–29760, and the lowest – 18,153 – in T8. Among the strains, there existed unusually high counts of Het ALT/REF alleles nearing 20,000 to 30,000 in a few strains. The average insertion/deletions variants in these strains were 24,060. While the Het REF/ALT alleles were in abundance, the Hom ALT and Het ALT alleles showed sparse distribution. Allele

frequency ranged from strain MHOM/IN/2014/T10, from Jharkhand, having the most with 314 Hom AA and 16 Het AA alleles and strain MHOM/IN/2014/A2, from West Bengal, having the least with 165 Hom AA and 8 Het AA alleles. Other Indian strains such as the BHU1220, BHU575/0/cl2, MHOM/IN/1983/AG83 and K133 strains showed markedly higher frequency of Hom AA and Het AA while exhibiting negligible Het RA alleles relative to the strains mentioned previously.

The Nepalese strains BPK275/0/cl18 showed limited InDels relative to the Indian strains, while also exhibiting a higher frequency of Hom ALT and Het ALT alleles. The strain BPK275/0/cl18.vii showed lower amounts of insertion/deletions in comparison to other Nepalese strains. The sample of Himachal state of India has around 1 lakh missense and 40000 InDels.

To decipher the distribution of different variant types within the strains used in this study, we utilized snpEff for annotation. The method involves identifying and attaching the relevant sequence ontological terms to the variants, as well as predicting their effect types and prospective impacts. Most of the variants appeared to be impacting non-coding genomic regions which included the 5’Flank, 3’Flank, and IGR. We observed an abundance of Frame_Shift_Del events and Missense_Mutation which impacted the coding portion of the genome. Few variants which impacted the coding regions were associated with in-frame deletion or insertion, missense and silent mutations. A tiny proportion of the variants represented intron, nonsense, nonstop, splice region, splice site and targeted region classes (Fig. 5). We looked at the distribution of different variant types across strains. The majority of the variants were in 3’ Flank or 5’ Flank across the strains that constituted 75.1% of the total variants. The strains had an average of 7.43% IGR variants. A relatively small share of the variants affected the coding regions of the gene sequences, these were Frame_Shift_Del, Frame_Shift_Ins, and Missense Mutations, contributing to an average of 8.37%, 1.18% and 4.54% respectively (Fig. 6).

When we compared the average variant proportion between Indian and Nepalese strains we found there to be significant differences between Frame_Shift_Del variants, 3’Flank and IGR regions.

We observed that Indian strains MHOM/IN/2010/, MHOM/IN/2012/, MHOM/IN/2013/, and MHOM/IN/2014 exhibited substantial amounts of variants tagged as ‘Frame_Shift_Del’, except for the BHU 1220, BHU575/0/cl2, MHOM/IN/1983/AG83 and K133 strains. This might be due to the presence of a higher number of insertion/deletion events, potentially leading to an increased representation of these variants in the samples. This was supported by the higher frequency of these variants relative to other strains in the population. It ranged from 4364 in MHOM/IN/2013/M1 and 1497 in MHOM/IN/2014/H1. The proportion also remained high at 27 to 18% of the total variants in these samples. Frame_Shift_Del events in other Indian strains remained lower; only occurring as less than 1% of the total variants in the samples. The frequency relative to other Indian strains previously mentioned also remained significantly low i.e. 30, 41, 44, 53 respectively. This might be due to the lesser insertion/deletion events in these strains which could be contributing to remarkably higher frame shift events. Whereas the Nepalese strains exhibited less than 1% of these variants, except strain BPK/275/0/cl18.vii which had 5% of ‘Frame_Shift_Del’ variants, which could be possible due to relatively lower insertion/deletion events in these strains. Nepalese strains also showed a higher mean share of Missense_Mutations, around 5.59%, than the Indian strains which showed 3.78%. IGR also contributed to 9.74% in Nepalese strains while 5.76% in Indian strains. In the Himalayan region sample, downstream and upstream variants comprises of 36% and 35% of the total variants and 9.2% in the exon region.

Single nucleotide polymorphisms impact functional domains

We used the Mutation Annotation Format (MAF) files generated by the R package MAFtools to investigate the abundance of ‘Missense variants’ in different isolates. We observed a median of 310 missense variants in the isolates affecting an array of genes, mainly imparting function as proteins involved in biosynthesis pathways, transport, anchorage, and as receptor proteins involved in cellular uptake of essential amino acids, interaction with host immune system and drug resistance. Among the strains from the Indian subcontinent those from Nepal showed a relatively high frequency of missense variations. Although the Indian isolate K133WT exhibited the highest missense mutation frequency, also K133 and BHU575/0/cl2 strains had a higher frequency of mutations which could be due to the exposure from Paromomycin (PMM) drugs (Fig. 7). Fifty genes were found to be affected by missense variants in the entire population of

isolates (Fig. 8). These included genes of ATP-binding cassette (ABC) family, Amastin-like surface proteins, A2 genes, amino acid permeases, and others. ABC subfamily A, member 2 (ABCA2) and member 9 (ABCA9) genes showed the most mutations, 503 and 326 respectively. We observed an increase in copy number of ABC subfamily A, member 9 protein (LDBPK_270840) in all the isolates. ABCA9 protein has two ‘ABC transporter domains’ from amino acids (aa) 760–993 and 1547–1777; while none of the domains were affected by variants, several missense variants of ‘Moderate Impact’ were observed in the gene sequence. It has been reported that Leishmania modulates aneuploidy in response to drug pressure, as ABC family has been shown to be present in higher copy number. ABC family of proteins are membrane transporter proteins which have been observed to be involved in drug resistance against antileishmaniasis through efflux (Ouellette, Légaré et al. 1998). Three amastin genes that exhibit high mutation count; Amastin-like surface protein (LDBPK_341150), Amastin surface glycofamily protein (LDBPK_044340), and Amastin surface glycofamily protein (LDBPK_140500) were among that constituted high number of mutations. It has been proposed that Amastin-like surface glycoproteins play a crucial role in macrophage-parasite interaction and parasite virulence. As reported by in a study wherein they knocked down the δ-Amastin expression in L. braziliensis using RNA-interference and observed decreased intracellular parasites in amastigotes in-vitro and in-vivo. We observed 2 missense mutations and 3 silent mutations in LDBPK_341150, 2 missense and 1 silent mutation in LDBPK_044340 and 1 missense mutation in LDBPK_140500. The 3’a2rel (LDBPK_220560) and 5’a2rel (LDBPK_220660) related proteins of the A2-A2rel cluster showed an increased missense variant count, the A2rel genes have been predicted to be involved in the visceralization of the leishmaniasis infection, as reported in a study on BALB/c mice by infection of L. major Friedlin V9 strain which was transfected with A2 expressing plasmids .The 3’a2rel sequence showed 7 moderate impact missense variants, 2 silent mutations and the 5’a2rel sequence showed missense, a modifier 5’ Flank and silent mutations. Several surface protein encoding genes showed high frequency of missense mutations (Fig. 9) including the Phosphoglycan-ꞵ-1,3-Galactosyltransferase, partial protein (LDBPK_020150), Mannosyltransferase-like protein (LDBPK_311920) and Amastin-like surface protein (LDBPK_341150). Several studies have implicated the Lipophosphoglycan (LPG), which are glycolipid surface proteins, which play a key role in the propagation of the parasite, through the insect vector, by facilitating the interaction of promastigotes with midgut epithelial lining of sandfly. (Ilg T., Sacks D et al.) LPGs are involved in parasite defence against the host innate immune system by avoiding the formation of membrane attack complexes. These glycoconjugates have also been shown to have a role in promoting cellular uptake of Leishmania parasites into the host macrophages, through the complement system, where they take refuge and multiply inside the acidic phagolysosomes. (SM Puentes et al., PJ Green et al.) In response to the selective pressures by the host environment and immune system Leishmania parasites alter surface physiology by altering LPG side chains which contribute to enhanced survivability and pathogenicity. (Roditi I et al.) The geographical distribution of sandfly influences species-specific LPG alterations and contribute to the successful transmission of parasite from sandfly vector to vertebrate host. (Ilg T.) Differential host-species and host-strain interactions could be contributing to the observed mutational frequency in the Galactosyltransferase gene. Mannosyltransferases are a type of galactosyltransferase which act on mannose and produces Mannogen, a unique carbohydrate material soluble in the cytoplasm, a homopolymer of ꞵ-1,2-linked mannose, 4–40 residues long. (Ralton et al., Sernee et al.) Mannogen is constitutively produced by seven Mannosyltransferases/phosphorylases (MTPs). Infective stage amastigotes have accumulated Mannogen, as it has been reported that defective Mannogen biosynthesis can lead to loss of virulence and attenuated growth in macrophages. (Sernee, M. Fleur et al.) There exist 𝛼-mannosyltransferases that catalyse the synthesis of N-linked glycans and glycosylphosphatidylinositol (GPI) glycolipids which promotes anchoring of surface proteins, including LPGs, by GPI. (Hilley et al.) Among the multi-hit top mutated genes were Amino acid permease, putative (LDBPK_270530), Tryparedoxin peroxidase (LDBPK_151140), cAMP specific phosphodiesterases (LDBPK_151450). Arginine is an essential amino acid for Leishmania hence it needs to be sourced from exogenous sources. Arginine uptake and homeostasis is maintained by transmembrane proteins such as amino acid permeases which transport Arginine on depletion in the parasite. (Opperdoes FR et al., Darlyuk, Ilona et al.) Tryparedoxin peroxidase enzyme, which is upregulated in the presence of Reactive Oxygen Species (ROOs), helps mitigate cellular toxicity by eliminating toxic agents responsible for Ca2 + influx, also an increase in Tryparedoxin peroxidases leads to lesser responsiveness to anti-leishmanial Antimony. (Iyer JP et al.) Various membrane and surface proteins have been implicated in the pathogenicity of Leishmania which also showed high frequency of missense variations, such as the ATP-binding cassette protein subfamily A, member 2, putative (LDBPK_111210) and member 9, putative (LDBPK_111210), 3’a2rel-related protein (LDBPK_220560). The A2 gene family is essential for visceralization of the disease as it has been observed to be absent in Leishmania major which causes Mucocutaneous Leishmaniasis (MCL) but is present in L. donovani which results in the visceral manifestation of the disease. (Zhang et al.). Next, we investigated if any of the mutations affected the functional domains of protein products of the gene sequences. We found 20 genes which had a reported functional domain that happened to be impacted by the missense mutations. The findings are reported in the Table 1.

Indian and Nepalese L. donovani isolates have distinct mutational profiles

We checked for genes which were affected by indels and found 370 genes common to all isolates. Our analysis revealed differential mutation profiles for each gene in different isolates. A cluster analysis of isolates was performed based on Indel mutation profiles and the results are shown in the below. (Fig. 10)

We observed the clustering of Leishmania donovani strains into three main groups, one of the clusters showed exclusively Indian isolates, cluster A, while the other two showed an intermixing of both Indian and Nepalese isolates, cluster B and C.

Cluster A comprised the strains MHOM/IN/2010/T3, MHOM/IN/2010/T7, MHOM/IN/2012/T8, MHOM/IN/2012/T9, MHOM/IN/2013/P2, MHOM/IN/2013/P2/M1, MHOM/IN/2014/A2, MHOM/IN/2014/H1, MHOM/IN/2014/P6, and MHOM/IN/2014/T10. These strains showed increased mutations in Bridge-like lipid transfer protein family member 1 N-terminal domaining containing protein (LDBPK_160180), Mannosyltransferase like-protein (LDBPK_311920), Calpain-like cysteine peptidase, putative (LDBPK_270510, ATP-binding cassette protein subfamily A, member 2, putative (LDBPK_111210), and uncharacterized proteins (LDBPK_333030-LDBPK_333040). Strains M1, P2 and T10 isolated from West Bengal, Bihar and Jharkhand grouped together since they had frequency of variants in the above-mentioned genes.

Cluster B comprised Indian and Nepalese isolates; BPK275/0/cl18.i, BPK275/0/cl18.ii, BPK275/0/cl18.iii, BPK275/0/cl18.vi, BPK275/0/cl18.viii, BPK275/0/cl18.ix,BPK275/0/cl18.xi, MHOM/IN/1983/AG83.i, MHOM/IN/1983/AG83.ii, and K133AS-R, K133PMM-R and K133WT. This cluster exhibited higher mutations in Uncharacterized protein (LDBPK_310470), uncharacterized proteins (LDBPK_333030-LDBPK_333040), Calpain-like cysteine peptidase, putative (LDBPK_270500-LDBPK_270510), Phosphoglycan beta 1,3 galactosyltransferase (LDBPK_020150), Amastin-like surface protein, putative (LDBPK_341690), Amino acid permease, putative (LDBPK_272680), Amastin-like surface protein, putative (LDBPK_341720), Amastin surface glycofamily protein-CHR_END (LDBPK_044360-CHR_END), more than 30 variants were seen in these genes. A few more genes were observed that contained mutations, not as many as the above mentioned, these were; pseudo (LDBPK_310470), Uncharacterized protein (LDBPK_120661-LDBPK_120667), ATP-binding cassette protein subfamily A, member 2, putative (LDBPK_111210), Uncharacterized protein (LDBPK_093020), Tyrosine phosphatase family protein (LDBPK_321070), Amastin surface glycofamily protein (LDBPK_044340), and Uncharacterized protein (LDBPK_333030-LDBPK_333040).

Cluster C comprised Indian and Nepalese isolates; BHU 1220, BHU 575/0/cl2, BPK275/0/cl18.iv, BPK275/0/cl18.v, BPK275/0/cl18.vi, BPK275/0/cl18.x, BPK275/0/cl18.xii, and BPK275/0/cl18.xiii. Strains showed similar mutation profiles to the strains in cluster B but distinct to the cluster A strains. Isolates showed a peculiar abundance in mutation in the Amastin surface glycoprotein (LDBPK_044360-CHR_END) gene and the Amastin-like surface protein, putative (LDBPK_341720).

Indels aggregate in specific regions of the L. donovani genome

Here we provide a comprehensive report on the Indel variation in the population of Nepalese and Indian isolates (Table 2). 91 indels were seen to be occurring in all the samples, with different variant classifications, in noncoding and coding regions of the genome. Only 7.6% i.e. 7 indels affected the coding regions. These variants were classified as either ‘In_Frame_Ins’ or ‘In_Frame_Del’, meaning they existed inside the reading frame of the genes and were size multiples of three, thereby preserving the reading frame. Whereas 92%, i.e. 84 indels were in noncoding regions and mostly were size multiples of two. Indels in coding regions may have a severe impact on the structure and function of gene products and hence are negatively selected, which explains the lower frequency of reading frame disrupting indels. Non-coding indels may be better tolerated due to more lenient selection pressure. It has been proposed that non-coding variants may play a role in influencing the gene expression of upstream and downstream genes. 91 variants were spread over coding and non-coding regions of 54 genes. These variants were observed to be occurring in and around specific nucleotide positions of gene sequences. The ATP-binding cassette protein, family A, member 2 (ABCA2) showed 6 frameshift indels, 585_586insTA, 567_568insCTTGC, 557_558insA, 580_581delCA, 572delA, and 559_563delGGTGT, all of which push the reading frame of the coding sequence and are concentrated around nucleotide (nt) 550 to 590. Hypothetical predicted transmembrane protein (LDBPK_240440) had three in-frame indels and two frameshift indels concentrated at nt 1270 to 1290. IGR between two Amastin-like proteins, LDBPK_080760 and LDBPK_080780 had four indels, which could possibly be involved in influencing gene expressions. Indels were seen in the 5’ Flank region of the Histone H2 protein. 3',5'-cyclic nucleotide phosphodiesterase, putative indels were at 5’Flank upstream of the gene sequence. Sodium stibogluconate resistance protein, putative had indels which resulted in reading frame shift which were predicted to be of high impact. This trend continued in the rest of the genes, thereby indicating that indels tend to arise in specific regions of the gene sequences most prone to alterations.

Copy number alterations

We observed that all of the chromosomes were predicted to have local copy number alterations in them. Chromosomes 36 through 26, 24, 22, 21, 16, 12, 7 and 2 were observed to have developed variations in copy number in all 31 samples. While chromosome 1, 3 and 4 were predicted to have local copy number variations in 5, 14 and 17 samples respectively. We checked genes that were associated with loss in copy number and the frequency of samples they occurred in. We observed that most of the genes showing a loss in copy numbers, that occurred in ≥ 30 samples, existed on chromosomes 2, 3, 7, 8, 9, 16, 22, 26 through 36. We checked the same for genes showing a gain in copy number and observed that the same chromosomes, which showed increased loss in copy number alterations, also had increased gain in copy number events. We found 393 genes with a gain in copy numbers; the highest number of CNV alterations were predicted in the chromosome 8 and 31 which may be because of their increased chromosomal copy number.

In this study, WGS data of the Indian sub-continent has been analyzed to investigate the SNP, CNV and ploidy distribution in this region. There are few reports on Indian samples that deals with the genome-wide analysis of paromomycin (PMM) resistant L. donovani and genetic characterization of L. donovani variant from Himachal Pradesh (Ghosh, Kumar et al. 2022, Lypaczewski, Thakur et al. 2022). The aneuploidy analysis has shown that the samples BPK* that belong to Nepal has shown different aneuploidy patterns as compared to other samples. The possible reason is that L. donovani samples considered for analysis are PMM resistance that act as stimulus for aneuploidy modification. K133WT (Delhi, India) sample that has maximum number of variants has shown triploidy at chromosome 5. BHU.1220 sample collected during 1992 shows tetraploidy at the chromosome 31 only. Therefore, the analysis clearly indicates that Leishmania has evolved over the time i.e. from 1992 to 2019 and there has been change in the ploidy in response to drug resistance and other stress conditions. The samples from Bihar, Jharkhand and West Bengal have shown triploidy in chromosome 31, however two samples have shown triploidy. It is evident that L. donovani has different genetic constitution in Nepal area than Indian states (Fig. 10). Drug pressure in K133* sample and BPK* sample may lead to genetic drift in the parasite population that lead to change in ploidy. (Fig. 11)

L. donovani sample from Himachal has shown highest number of variants and this is attributed to its genetic hybridization between L. donovani parents from ‘Yeti’ clade reported from Nepal highlands. The SNPs that were present in either samples were removed and the novel SNPs present in this sample has been discussed in the following text. The details of the variants of the sample are given in Supplementary Table 3. It was observed that the following proteins were most mutated in the sample of Himachal namely, LdBPK_111210 encoding an ABC protein subfamilyA, LDBPK_060660 protein kinase and putative LdBPK_171090.1, LDBPK_281830, LDBPK_361500, LDBPK_091050, LDBPK_161130, LDBPK_262120, LDBPK_251150, LDBPK_191110 genes encode for hypothetical proteins (Supplementary Fig. 3). It is evident that the proteins which have shown more mutations are hypothetical proteins and need to be characterized. However, this sample follows the same pattern as other samples in case of higher mutation rate in transporter protein. Since this sample organism has been reported as an intraspecies hybrid two distinct clades of L. donovani from Ethiopia can be the reason that there are numerous novel mutations (Lypaczewski, Thakur et al. 2022). However in all samples of Indian sub-continent transporter protein, amastin like proteins, protein involved in drug resistance, heat shock proteins and hypothetical proteins have shown maximum variants. The details of the variants present in all samples are given in Supplementary Table 4 and the name of protein encoded by corresponding gene is given in Supplementary Table 5. There is genetic plasticity in Leishmania genome, CNV affects either some genomic region or whole chromosomes as observed for chromosome 8 and 31 that leads to aneuploidy. In recent studies, only chromosome 31 has shown high number of CNVs, however the present analysis shows high number of CNVs in chromosome 8 and 31 (Ghosh, Kumar et al. 2022, Lypaczewski, Thakur et al. 2022). This analysis of Indian sub-continent data is providing a wholesome image of CNV which is different from ealier obsevations. The proteins which have shown maximum missense mutations are transporter protein, proteins involved in carbohydrate metabolism, heat shock protein, amastin proteins, sodium stibogluconate resistant protein and hypothetical proteins. While in the case sample from Delhi, India, mutations were observed in transcription, translation (60S ribosomal protein L22, 60S ribosomal protein L40, RNA binding protein, eukaryotic translation initiation factor 4E-1 and 4E-2), virulence factor (protein-tyrosine phosphatase), and signaling molecule (protein kinase and protein-tyrosine phosphatase) (Ghosh, Kumar et al. 2022). 42 ABC genes namely ABCA to ABCH has been reported. In previous studies, membrane/transporter proteins like ABCG1, ABCG2, ABCG3, ABCA2, ABCA7, and ABCA8 has shown missense mutations (Ghosh, Verma et al. 2020). However, in the present study, ABC subfamily A, member 2 (ABCA2) and member 9 (ABCA9) has shown maximum missense mutations which are primarily involved in drug efflux and responsible for multi-drug resistance. Two transmembrane glycorotein i.e. Amsatin-like surface protein and amastin surface glycofamily protein has shown maximum mutations after transporter proteins. As reported in earlier studies, amastin-like protein is reported to be an important precursor for the phenotype resistance development against antimonial and paramomycin (Ghosh, Kumar et al. 2022). Mannosyltransferase protein has shown high number of missense mutations. The overexpression of Mannosyltransferase protein in promastigote and amastigote forms of L. braziliensis protects from trivalent antimony (Ribeiro, Rocha et al. 2019). Heat shock protein 70-related protein 1, putative, mitochondrial has high missense mutations and in the earlier studies, other types of HSP proteins has shown mutations or altered expression and are also linked to miltefosine resistance (Vacchina, Norris-Mullins et al. 2016, Alcolea, Alonso et al. 2023). Leishmaniasis treatment failure to sodium stibogluconate has been described in some studies which are related to expression of γ-Glumatylcysteine synthase, however there is no record of the missense mutations of sodium stibogluconate resistant protein putative. It was observed that sodium stiboglocunate resistant protein putative is among the top five mutated protein in the samples of Indian subcontinent. The mutated proteins observed from WGS data of L. donovani from ISC may act as initial step for further validation of these proteins in the virulence and understanding the drug resistance mechanism of these proteins.

Owing to the changing climate conditions, Indian sub-continent geographical area is more prone to these changes. The analysis of the clinical samples of L. donovani from Indian sub-continent has provided novel insights into the genetic makeup of the protozoan present here. The overall trend of the mutations in the protein of L. donovani from Indian sub-continent shows that there is increasing mutations in the proteins which are involved in the virulence and drug-resistant mechanism. The membrane proteins has shown most of the variants as they throw out the drugs from the pathogen by efflux mechanism. CNV analysis has shown different pattern of aneuploidy in the L. donovani from Indian sub-continent when compared to the earlier studies from other part of the world. Therefore, the proteins and variants obtained can be the founding steps to develop strategic drugs for the treatment of Leishmaniasis.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Author Contribution

AV has curated the data, carried out the WGS analysis, and written the first draft of the paper. IK: Streamlined the analysis workflow, written the results and discussion. ANJ: has checked, edited, and given suggestions in the manuscript.

Acknowledgement

ANJ gratefully acknowledges the ICMR project for providing computing facilities.

Alcolea, P. J., A. Alonso, F. García-Tabares, J. Larraga, L. T. Martins, F. J. Loayza, S. Ruiz-García and V. Larraga (2023). "An insight into differential protein abundance throughout Leishmania donovani promastigote growth and differentiation." International Microbiology26(1): 25-42.
Alvarenga, D. G. d., P. M. F. Escalda, A. S. V. d. Costa and M. T. F. D. Monreal (2010). "Leishmaniose visceral: estudo retrospectivo de fatores associados à letalidade." Revista da Sociedade Brasileira de Medicina Tropical43: 194-197.
Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data, Babraham Bioinformatics, Babraham Institute, Cambridge, United Kingdom.
Awasthi, A., R. K. Mathur and B. Saha (2004). "Immune response to Leishmania infection." Indian Journal of Medical Research119: 238-258.
Boeva, V., T. Popova, K. Bleakley, P. Chiche, J. Cappo, G. Schleiermacher, I. Janoueix-Lerosey, O. Delattre and E. Barillot (2012). "Control-FREEC: a tool for assessing copy number and allelic content using next-generation sequencing data." Bioinformatics28(3): 423-425.
Bolger, A. M., M. Lohse and B. Usadel (2014). "Trimmomatic: a flexible trimmer for Illumina sequence data." Bioinformatics30(15): 2114-2120.
Bora, N. and A. N. Jha (2020). "In silico metabolic pathway analysis identifying target against leishmaniasis–a kinetic modeling approach." Frontiers in Genetics11: 179.
Bora, N. and A. Nath Jha (2019). "An integrative approach using systems biology, mutational analysis with molecular dynamics simulation to challenge the functionality of a target protein." Chemical biology & drug design93(6): 1050-1060.
Cingolani, P., A. Platts, L. L. Wang, M. Coon, T. Nguyen, L. Wang, S. J. Land, X. Lu and D. M. Ruden (2012). "A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3." fly6(2): 80-92.
Clark, K., I. Karsch-Mizrachi, D. J. Lipman, J. Ostell and E. W. Sayers (2016). "GenBank." Nucleic acids research 44(D1): D67-D72.
Danecek, P., J. K. Bonfield, J. Liddle, J. Marshall, V. Ohan, M. O. Pollard, A. Whitwham, T. Keane, S. A. McCarthy and R. M. Davies (2021). "Twelve years of SAMtools and BCFtools." Gigascience 10(2): giab008.
de Paiva, R. M. C., V. Grazielle-Silva, M. S. Cardoso, B. N. Nakagaki, R. P. Mendonca-Neto, A. M. C. Canavaci, N. Souza Melo, P. M. Martinelli, A. P. Fernandes and W. D. daRocha (2015). "Amastin knockdown in Leishmania braziliensis affects parasite-macrophage interaction and results in impaired viability of intracellular amastigotes." PLoS pathogens 11(12): e1005296.
Downing, T., H. Imamura, S. Decuypere, T. G. Clark, G. H. Coombs, J. A. Cotton, J. D. Hilley, S. de Doncker, I. Maes and J. C. Mottram (2011). "Whole genome sequencing of multiple Leishmania donovani clinical isolates provides insights into population structure and mechanisms of drug resistance." Genome research 21(12): 2143-2156.
Dumetz, F., H. Imamura, M. Sanders, V. Seblova, J. Myskova, P. Pescher, M. Vanaerschot, C. J. Meehan, B. Cuypers and G. De Muylder (2017). "Modulation of aneuploidy in Leishmania donovani during adaptation to different in vitro and in vivo environments and its impact on gene expression." MBio8(3): 10.1128/mbio. 00599-00517.
Ghosh, S., V. Kumar, A. Verma, T. Sharma, D. Pradhan, A. Selvapandiyan, P. Salotra and R. Singh (2022). "Genome-wide analysis reveals allelic variation and chromosome copy number variation in paromomycin-resistant Leishmania donovani." Parasitology Research 121(11): 3121-3132.
Ghosh, S., A. Verma, V. Kumar, D. Pradhan, A. Selvapandiyan, P. Salotra and R. Singh (2020). "Genomic and transcriptomic analysis for identification of genes and interlinked pathways mediating artemisinin resistance in Leishmania donovani." Genes11(11): 1362.
Joshi, A. B., M. R. Banjara, M. L. Das, P. Ghale, K. R. Pant, U. R. Pyakurel, G. Dahal, K. P. Paudel, C. L. Das, A. J. T. A. J. o. T. M. Kroeger and Hygiene (2024). "Epidemiological, Serological, and Entomological Investigation of New Visceral Leishmaniasis Foci in Nepal." The American Journal of Tropical Medicine and Hygiene110(1): 44.
Kumari, I., D. Lakhanpal, S. Swargam and A. Nath Jha (2023). "Leishmaniasis: Omics approaches to understand its biology from molecule to cell level." Current Protein and Peptide Science24(3): 229-239.
Kumari, P., A. Mamud and A. N. Jha (2023). "Review on the Drug Intolerance and Vaccine Development for the Leishmaniasis." Current Drug Targets24(13): 1023-1031.
Li, H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin and G. P. D. P. Subgroup (2009). "The sequence alignment/map format and SAMtools." bioinformatics 25(16): 2078-2079.
Lira, R., S. Sundar, A. Makharia, R. Kenney, A. Gam, E. Saraiva and D. Sacks (1999). "Evidence that the high incidence of treatment failures in Indian kala-azar is due to the emergence of antimony-resistant strains of Leishmania donovani." The Journal of infectious diseases 180(2): 564-567.
Lypaczewski, P., L. Thakur, A. Jain, S. Kumari, K. Paulini, G. Matlashewski and M. Jain (2022). "An intraspecies Leishmania donovani hybrid from the Indian subcontinent is associated with an atypical phenotype of cutaneous disease." Iscience 25(2).
Mann, S., K. Frasca, S. Scherrer, A. F. Henao-Martínez, S. Newman, P. Ramanan and J. A. Suarez (2021). "A review of leishmaniasis: current knowledge and future directions." Current tropical medicine reports 8: 121-132.
Martínez, F. V., E. T. Guerrero, R. Arenas and M. Q. Cedillo (2011). "Leishmaniasis en México." Medicina Cutánea Ibero-Latino-Americana 39(4): 163-183.
Mayakonda, A., D.-C. Lin, Y. Assenov, C. Plass and H. P. Koeffler (2018). "Maftools: efficient and comprehensive analysis of somatic variants in cancer." Genome research 28(11): 1747-1756.
Ouellette, M., D. Légaré, A. Haimeur, K. Grondin, G. Roy, C. Brochu and B. Papadopoulou (1998). "ABC transporters in Leishmania and their role in drug resistance." Drug Resistance Updates 1(1): 43-48.
Patil, R. R. and P. K. Chatterjee (2023). "Epidemiology of Visceral Leishmaniasis in India."
Paysan-Lafosse, T., M. Blum, S. Chuguransky, T. Grego, B. L. Pinto, G. A. Salazar, M. L. Bileschi, P. Bork, A. Bridge and L. Colwell (2023). "InterPro in 2022." Nucleic acids research 51(D1): D418-D427.
Pilling, O. A., J. L. Reis-Cunha, C. A. Grace, A. S. Berry, M. W. Mitchell, J. A. Yu, C. R. Malekshahi, E. Krespan, C. K. Go and C. Lombana (2023). "Selective whole-genome amplification reveals population genetics of Leishmania braziliensis directly from patient skin biopsies." PLoS Pathogens 19(3): e1011230.
Reuter, J. A., D. V. Spacek and M. P. Snyder (2015). "High-throughput sequencing technologies." Molecular cell 58(4): 586-597.
Rezaei, Z., A. Tahmasebi and B. Pourabbas (2024). "Using meta-analysis and machine learning to investigate the transcriptional response of immune cells to Leishmania infection." PLOS Neglected Tropical Diseases 18(1): e0011892.
Ribeiro, C. V., B. F. B. Rocha, D. d. S. Moreira, V. Peruhype-Magalhães and S. M. F. Murta (2019). "Mannosyltransferase (GPI-14) overexpression protects promastigote and amastigote forms of Leishmania braziliensis against trivalent antimony." Parasites & vectors12: 1-7.
Riebenbauer, K., S. Czerny, M. Egg, N. Urban, T. Kinaciyan, A. Hampel, L. Fidelsberger, F. Karlhofer, S. Porkert and J. Walochnik (2024). "The changing epidemiology of human leishmaniasis in the non-endemic country of Austria between 2000 to 2021, including a congenital case." PLOS Neglected Tropical Diseases 18(1): e0011875.
Sádlová, J., M. Yeo, D. S. Mateus, J. Phelan, L. A. Hai, T. Bhattacharyya, S. Kurtev, O. Sebesta, J. Myskova and V. Seblova (2024). "Comparative genomics of Leishmania donovani progeny from genetic crosses in two sand fly species and impact on the diversity of diagnostic and vaccine candidates." PLOS Neglected Tropical Diseases 18(1): e0011920.
Saha, D., N. J. Borah and A. N. Jha (2023). "Molecular scaffold recognition of drug molecules against essential genes of Leishmania donovani using biocomputing approach." South African Journal of Botany 162: 52-63.
Santi, A. M. M. and S. M. F. Murta (2022). "Impact of genetic diversity and genome plasticity of Leishmania spp. in treatment and the search for novel chemotherapeutic targets." Frontiers in Cellular and Infection Microbiology12: 9.
Sapp, S. G., R. Low, G. Nine, F. S. Nascimento, Y. Qvarnstrom and J. L. Barratt (2024). "Genetic characterization and description of Leishmania (Leishmania) ellisi sp. nov.: a new human-infecting species from the USA." Parasitology Research 123(1): 52.
Schlein, Y. (1993). "Leishmania and sandflies: interactions in the life cycle and transmission." Parasitology Today 9(7): 255-258.
Shirley, M. D., Z. Ma, B. S. Pedersen and S. J. Wheelan (2015). Efficient" pythonic" access to FASTA files using pyfaidx, PeerJ PrePrints.
Sofi, O. M. U. D., R. Godara, R. Katoch and A. Yadav (2023). "Morphological and molecular evidence of sandfly (Diptera: Psychodidae: Phlebotomine) and its relevance to recent cases of leishmaniasis from Jammu region of North India." Journal of Vector Borne Diseases 60(3): 328-332.
Torres-Guerrero, E., M. R. Quintanilla-Cedillo, J. Ruiz-Esmenjaud and R. Arenas (2017). "Leishmaniasis: a review." FResearch 6.
Vacchina, P., B. Norris-Mullins, E. Carlson and M. Morales (2016). "A mitochondrial HSP70 (HSPA9B) is linked to miltefosine resistance and stress response in Leishmania donovani." Parasites & vectors 9: 1-15.
Van der Auwera, G. A., M. O. Carneiro, C. Hartl, R. Poplin, G. Del Angel, A. Levy‐Moonshine, T. Jordan, K. Shakir, D. Roazen and J. Thibault (2013). "From FastQ data to high‐confidence variant calls: the genome analysis toolkit best practices pipeline." Current protocols in bioinformatics 43(1): 11.10. 11-11.10. 33.
Wang, K., M. Li and H. Hakonarson (2010). "ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data." Nucleic acids research 38(16): e164-e164.
Zhang, W. W. and G. Matlashewski (2001). "Characterization of the A2–A2rel gene cluster in Leishmania donovani: involvement of A2 in visceralization during infection." Molecular microbiology 39(4): 935-948.
Bharadava, K., Upadhyay, T. K., Kaushal, R. S., Ahmad, I., Alraey, Y., Siddiqui, S., & Saeed, M. (2024). Genomic Insight of Leishmania Parasite: In-Depth Review of Drug Resistance Mechanisms and Genetic Mutations. ACS omega, 9(11), 12500-12514.
Al Khoury, C., Thoumi, S., Tokajian, S., Sinno, A., Nemer, G., El Beyrouthy, M., & Rahy, K. (2024). ABC transporter inhibition by beauvericin partially overcomes drug resistance in Leishmania tropica. Antimicrobial Agents and Chemotherapy, 68(5), e01368-23.

Tables 1 and 2 are available in the Supplementary Files section.

No competing interests reported.

Download PDF

Journal Publication

published 19 Nov, 2025

Read the published version in Antonie van Leeuwenhoek →

Editorial decision: Revision requested
22 Aug, 2025
Reviews received at journal
22 Aug, 2025
Reviews received at journal
07 Aug, 2025
Reviewers agreed at journal
07 Aug, 2025
Reviewers agreed at journal
06 Aug, 2025
Reviewers invited by journal
06 Aug, 2025
Editor assigned by journal
01 Aug, 2025
Submission checks completed at journal
01 Aug, 2025
First submitted to journal
01 Aug, 2025

You are reading this latest preprint version

Insights into Genetic Heterogeneity and Drug Resistance in Leishmania donovani of the Indian sub-continent from genomic data

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methodology

Results

Discussion

CONCLUSION

Declarations

Declaration of competing interest

Author Contribution

Acknowledgement

References

Tables

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1