*3.10. Correlation of Whole-Genome Sequencing Analysis with Transcriptomic Data*

Advancement in genomics and transcriptomics technologies have conferred considerable enhancement of our knowledge related to the set of changes that occur within the parasite at the molecular level resulting in the evolution of drug resistance. Comparative analysis of the genome and transcriptome data of the two distinct strains of *L. donovani* (K133WT and K133AS-R) provided a correlation of the mechanism behind the development of artemisinin resistance. Eight upregulated and 10 downregulated genes in the transcriptome data matched with the NGS data. Out of eight upregulated genes, two amastin-like surface protein genes (LinJ08\_V3.0700/LdBPK\_080710) and (LinJ34\_V3.0700/LdBPK\_43111505) showed one and five mutations respectively; the remaining six genes had upstream gene variants. The types of mutations observed were insertion, deletion, and transition. Genes that correlated in both included the amino acid transporter ATP11 and cathepsin—L like cysteine protease. Out of 10 downregulated genes, ABCA2 (LinJ11\_V3.1230/LdBPK\_111210) and receptor-type adenylatecyclase b (fragment) (LinJ17\_V3.0140/LdBPK\_170120) displayed a frameshift mutation at two sites. The genes that were observed to be downregulated included phosphoacetylglucosamine mutase-like protein (LinJ07\_V30930/LdBPK\_070930) and pteridine transporter (LinJ06\_V3.1320/LdBPK\_061320). Cathepsin L-like protease, a type of lysosomal endopeptidases, is present in both the promastigote and amastimogote stage of *Leishmania* species and involved in crucial biological process of parasites, such as evasion of the host immune system [63–65].

#### **4. Discussion**

Sesquiterpene, artemisinin, a secondary metabolite extracted from *Artemisia annua*, is an important antimalarial drug that has shown antimicrobial and antiviral activities [66,67]. Several in vitro and in vivo studies suggested potential antileishmanial activity of this drug [8–10,68]. However, the possibility of the emergence of resistance following the use of artemisinin as antileishmanial treatment cannot be denied. In our previous study, we reported that in vitro-selected artesunate-resistant *Leishmania* parasites were more virulent, successfully modulating the host cell defense mechanism, and exhibited altered expression of genes involved in the unfolded protein response, as compared to sensitive parasites [41]. The present study aimed to explore the genome and transcriptome of artesunate-resistant *Leishmania* parasites in order to understand the mechanism of resistance and to safeguard this drug for future use. Next-generation sequencing (NGS) platforms have advanced to provide a precise and comprehensive means for the detection of molecular mutations. Genomic and transcriptomic analyses would help in the advancement of our understanding of the biology of *Leishmania*. This comparative analysis of whole-genome sequences attempted to explicate genetic factors responsible for drug resistance in *L. donovani*. Here, we demonstrated that the in vitro-selected artesunate-resistant (K133AS-R) parasite was quite distinct from the sensitive wild-type (K133WT) at the genome and transcriptome level.

Major findings of the study are summarized in three sections. Firstly, from the genomic landscape, we found a high number of SNPs and InDel, many of them having a pronounced influence (stop codon gained/lost and frame-shifts) on essential biological functions. Briefly, in K133AS-R, upstream gene variants were higher followed by intergenic region gene variants. Out of the total number of gene mutations, SNPs were high as compared to InDel. The highest number of SNPs was observed on chromosome number 12, 31, 34, and 35. Among InDel mutations, insertions were greater than deletions. Non-coding mutations, such as upstream gene variants and downstream gene variants, affect regulatory elements and lead a to loss of function that results in reduced gene expression, or a gain of function resulting in differential gene expression [69]. In this study, we analyzed that selective forces are majorly acting on non-coding regions of the genome. Secondly, the major changes observed were concerned with local copy number variations (CNVs). In K133AS-R, higher deletion occurred as compared to duplication and CNV lengths in the range of 1–5 kbp were either deleted or duplicated, depicting that changes occurred at small sequences rather than larger sequences. The highest number of CNVs were observed in K133AS-R on chromosome no 31, 29, 20, and 18. In the absence of regulation of gene expression at the initiation site, duplication/deletion of specific genes in a genomic sequence modulates the transcript level and its products [69,70]. Complex chromosomal copy number variation is often observed in *Leishmania* parasites due to their asexual mode of replication [26]. Thirdly, second-generation sequencing data obtained with the Illumina analyzer sets out remarkable read depth coverage throughout the chromosomes of *Leishmania*, and both the K133AS-R and K133 WT exhibited a uniform read depth in all chromosomes that is disomic except chromosome number 12 in

K133 AS-R. A read depth greater than two-fold was observed in case of chromosome 12, suggesting that the chromosome is present in the trisomy condition. Chromosomy variation in *Leishmania* is a well-known adaptive strategy in response to experimental drug resistance selection [71]. Aneuploidy is mostly influenced by the environmental condition and is more prevalent in promastigotes under in vitro conditions than in amastigotes present inside the vertebrate host. It arises through unlicensed replication due to a lack of proper cell cycle regulation and/or mitotic non-disjunction [43]. Further, GO terms based functional annotation of genes lead to classification into different categories, including metabolic, cellular processes, and biological regulation, which include the response to stimulus, cell signaling, and growth. WGS data analysis gives ample information regarding genetic variation compared to other sequencing approaches that include SNPs, InDel, as well as structural variants *viz* CNVs, inversion, translocation, and ploidy variation in chromosomes [72].

Analysis of transcriptome data by microarray and further experimental validation of differentially expressed proteins resulted in several important findings. To maintain cellular homeostasis, the eukaryotic cells have developed specialized mechanisms, such as lysis of intracellular proteins and organelles, which regulate cellular functions like enzymatic activity, removal of toxic or misfolded proteins, and the production of free amino acids to ensure cell survival under stressful conditions. The eukaryotic cells are known to perform these functions by the process of autophagy, which is believed to have originated at a later point during evolution [73]. Artesunate causes high levels of ROS generation within the cell. Further, it has been reported in cancer cells that autophagy plays a cytoprotective role within cells by inhibiting ROS. In K133AS-R, downregulated expression of Atg8 suggested reduced inhibition of ROS. In addition, deceased expression of HSP70, both at transcript and protein levels, suggested an accumulation of a higher number of misfolded proteins, resulting in higher ER stress and finally higher ROS production in K133AS-R parasites. On the other hand, upregulated expression of lipoate protein ligase leads to upregulated expression of GSH and thus plays an important role in glutathione biosynthesis and response to oxidative stress. This may be a compensatory approach of K133AS-R parasites to survive under oxidative stress.

Downregulated expression of gene phosphoacetylglucosamine mutase, which is eventually involved in the formation of fructose-6-phosphate, an important component of glycolysis/ gluconeogenesis, suggested downregulation of these pathways in K133AS-R parasites. In addition, the downregulated expression of various glucose transporters suggested that K133AS-R parasites may not depend on carbohydrate metabolism for energy requirements. Hence, artesunate-resistant parasites may depend on amino acids and lipids for energy generation as inferred by the upregulated expression of methylmalonyl CoA mutase (involved in isoleucine, valine, and leucine metabolism) and glutamine aminotransferases (involved in glutamine metabolism) while myo-inositol-1-phosphate synthase and a hypothetical protein are involved in lipid metabolism.

Leucin-rich AMA1 protein secreted by many *Leishmania* species, including *L. donovani*, helps them to interact with cholesterol present in the host cell membrane and thereby assist the internalization of parasites [74–76]. Amastin, a transmembraneglycoprotein, encoded by a large gene family initially reported in the amastigote stage of trypanosomes and later observed as a surface protein expressed in *Leishmania* species (encoded by six copies of genes) plays an important role in visceralization [77]. The importance of amastin in the pathogenesis of *Leishmania* species is well documented in a previous study [55]. The data analysis showed that the parasites may undergo genomic alterations to express certain genes differentially to adapt to the drug-induced selection pressure.

K133AS-R parasites showed downregulated expression of several genes involved in DNA synthesis and translation machinery. Reduced DNA/protein synthesis leads to an arrest of parasites in a quiescent state, which may be responsible for drug resistance as reported in case of artemisinin resistance in malaria [78]. Further, there was a downregulated expression of metallo- and carboxy-peptidase involved in protein degradation, which may be an adaptive approach of K133AS-R parasites to overcome reduced protein synthesis.

In *Leishmania*, AQP1 plays an important role in providing nutrients from the host organism, mainly glucose, amino acids, and fatty acids. These may also be responsible for discarding waste and metabolic end-products, such as lactate, from the parasite's cytosol [79]. In the presence of AQP1 inhibitor, drug-resistant mutants showed a significant increase in susceptibility towards artesunate at the promastigote stage; however, no significant alteration in drug susceptibility was observed in drug-sensitive parasites, indicating an important role of AQP1 in the selection of artesunate resistance in *Leishmania*. AQP1 has been reported to be involved in the uptake of antimonial drugs and its downregulated expression has been found to be associated with drug resistance [80,81]. Interestingly, in artesunate resistance, higher expression of AQP1 both at the mRNA and protein levels was observed to be associated with drug resistance. Another interesting observation was the decrease in the susceptibility of drug-sensitive parasites towards artesunate at the intracellular amastigote stage, whereas no significant alteration in drug susceptibility was observed with drug-resistant parasites.

Higher expression of ABC transporters has been widely reported in drug resistance in *Leishmania* [24,81–83]. Upregulated expression of ABCG1 (ABCG subfamily) was observed in artesunate resistance. Further, the use of the ABC transporter verapamil resulted in a significant increase in the susceptibility of K133AS-R parasites towards artesunate both at the promastigote and amastigote stage, suggesting an important role of ABCG1 in the selection of drug resistance. However, functional characterization of ABCG1 needs to be carried out in order to establish its role in artesunate resistance. LPG5B (UDP-galactose transporter) plays diverse roles in parasite survival, like the control of parasite binding to the sand fly midgut wall, resistance to lysis by complement, protection from oxidative damage, and delayed fusion of phagolysosomes. Upregulated expression of LPG5B in K133AS-R may be helpful to parasites for their survival under drug pressure.

Comparative genome as well as transcriptome data analysis resulted in several major findings, such as upregulation of cathepsin L-like protease, amastin-like surface protein, and amino acid transporter at both the genome as well as the RNA level. Downregulated genes that were observed to be in sync with NGS data were ABCG2, Pteridine receptor, receptor-type adenylatecyclase, phosphoaceylglucosamine mutase-like protein, and certain hypothetical proteins.

Our data explicate a better insight in genomic and transcriptomics alteration that occurs during artemisinin stress under in vitro conditions and would act as a baseline for further studies involving the applicability of genomic changes encountered in the study of the clinically resistant and sensitive *L. donovani* isolated from patients of leishmaniasis. Overall, this study highlights genes and interlinked pathways contributing to artemisinin resistance using *Leishmania* as a model and highlights putative mechanisms that have applicability not only to malaria but also other diseases against which the drug is found to be effective.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/11/11/1362/s1, Figure S1: Comparative transcriptional responses following ART adaptation in *L. donovani*, Figure S2: Cytotoxicity of (A) AQP1 inhibitor and (B) Verapamil to host macrophages (mice PECs). Percentage cell viability ± SD with the increasing drug concentration has been plotted here. The data was obtained from three independent experiments, Table S1: List of genes validated for their modulated expression by Quantitative real time-PCR, Table S2: Pattern of up-regulated and down-regulated gene expression in K133 AS-R parasite, Table S3: Genes differentially modulated with their functional categories in K133 AS-R *L. donovani* parasites.

**Author Contributions:** Conceptualization, R.S., P.S., A.V. and S.G.; Methodology, S.G. and A.V.; Software, R.S. and D.P.; Validation, S.G., A.V., V.K. and D.P.; Formal Analysis, A.V., S.G. and V.K.; Investigation, S.G. and A.V.; Resources, R.S., P.S., A.S. and D.P.; Writing—Original Draft Preparation, A.V., S.G., V.K. and R.S.; Writing—Review & Editing, R.S., P.S. and A.S.; Supervision, R.S., P.S. and A.S.; Project Administration, R.S.; Funding Acquisition, R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Indian Council of Medical Research (ICMR, India) grant no. 53/16/2012-CMB/BMS and 6/9-7(188)2018-ECDII to RS. SG and AV are grateful to ICMR for financial support.

**Conflicts of Interest:** The authors declare no conflict of interest.
