*2.4. Whole-Genome Sequencing Data Analysis*

Genomic data analysis was executed with minor modifications as described previously by Dumetz et al. 2017 [43]. The paired-end (PE) raw reads obtained from the sequencer were checked for the quality of the reads using FastQCv0.11.8 (Babraham Institute, Cambridge, UK) and were further trimmed to improve the quality of the reads using the Trimmomatic tool v0.38 (Usadellab. org, RWTH Aachen University, Germany) [44]. The *L. donovani* strain LdBPK282A1 reference genome was indexed and high-quality pair-end reads were mapped using Burrows-Wheeler Aligner (BWA-MEM v0.7.5a algorithm, Broad Institute, Cambridge, MA, USA) [45]. The generated SAM file was converted into BAM format and duplicates were removed using Picard toolkit v1.119 (Broad Institute, Cambridge, MA, USA). Further, the BAM file was used for identifying SNPs and InDels using GATK Haplotype caller v3.4 (Broad Institute, Cambridge, MA, USA). Filtering of SNPs and InDels was performed using Bcf tools v0.1.18 (Sanger Institute, Cambridge shire, UK) from subdirectory of SAM tools (mapping quality cut off 25 and read depth of 15) and the variants were annotated with the SnpEff v4.3 tool (McGill University, Montreal, QC, Canada) [46,47].

Estimation of CNV along with chromosomal somy was done in accordance with the protocol designed by Downing et al. 2011 [29]. CNV estimation was done using CNVnator (https://github.com/ abyzovlab/CNVnator) [48], and for somy assessment, the median read depth of each chromosome (*d<sup>i</sup>* ) was computed first followed by median depth estimation of 36 complete chromosomes (*dm*). The somy state of an individual chromosome is determined as the ratio of (*d<sup>i</sup>* /*dm*) and the chromosome ploidy value is specified as 2 × (*d<sup>i</sup>* /*dm*), considered often for diploid species [38]. The full cell-normalized chromosome somy (S)-value: S < 1.5, 1.5 < S < 2.5, and 2.5 < S < 3.5, was assigned to monosomy, disomy, and trisomy, respectively [43].

#### *2.5. Functional Annotation and Classification of Unigenes*

To identify all the unigenes present in K133WT and K133AS-R, a homology search was performed against the NCBI non redundant (NR) protein database in accordance with BLASTx program (NCBI, Bethesda, MD, USA) using a cutoff E-value of 10−<sup>05</sup> and the maximal aligned results with the lowest E-value were chosen to annotate the unigenes [49,50]. The Gene Ontology (GO)-based annotation of the unigenes was carried out using Blast2GO version 3.0 (Biobam, Valencia, Spain) and Web Gene Ontology Annotation Plot (WEGO) was utilized to designate GO classification on the basis of the distribution of gene functions in different species [51–54]. The basis of the functional classification considered was biological processes, cellular components, and molecular functions.

#### *2.6. Total RNA Isolation from Parasites*

Early log-phase promastigotes (1 × 10<sup>8</sup> ) of both K133WT and K133AS-R were used to isolate total RNA using TRIzol reagent according to the manufacturer's instruction. Extracted RNA was cleaned up using a RNeasy Plus mini kit (Qiagen, Hilden, Germany). The absorbance of purified RNA was taken at 260 and 280 nm using a Nanodrop Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The quality and integrity of RNA were assessed on an RNA 6000 Nano Assay Chips on Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). RNA of good quality based on the 260/280 values (Nanodrop, Thermo Scientific, Waltham, MA, USA), rRNA 28S/18S ratios, and RNA integrity number (RIN) was used for further analysis [24].
