Chitwan district is located in the southwestern part of Bagmati Province in Nepal. It covers an area of 2,238.39 km² and had a human population of 719,859 in 2022. Nepal’s first national park, Chitwan National Park (CNP), is located in this district (Figure 1) and was established in 1973 (27). The park covers an area of 932 km² and is located in the subtropical Terai lowlands of south-central Nepal. It lies in a river valley basin or dun, along the floodplains of the Rapti, Reu and Narayani rivers. The Chitwan valley consists of tropical and subtropical forests. Sal forest covers 70% of the park. Sal leaves are used locally for plates in festivals and religious offerings, and are also known for their Ayurvedic use in the treatment of wounds, coughs and other ailments (28). Grasslands cover 20 per cent of the park. There are more than 50 different species of grasses, including elephant grass (Saccharum spp), known for its immense height. The climate is mainly dominated by the summer monsoon, and the valley experiences three distinct seasons each year: winter, summer and monsoon.

Chitwan National Park is regarded as one of the best national parks for wildlife viewing in Asia (29). It is home to 68 species of mammals, more than 576 species of birds, 49 species of reptiles and amphibians, 120 species of fish and several species of invertebrates that contribute significantly to ecosystem processes in the region. The reserve is home to many of Nepal’s charismatic wildlife species, including the one-horned rhinoceros (Rhinoceros unicornis), Bengal tiger (Panthera tigris), leopard (Panthera pardus) and Asian elephant (Elephas maximus), among others (29, 30). Of the various recreational activities available in the park, elephant safaris are the most popular, especially in and around the buffer zones. Elephants are also used by the government for transport and patrolling the park. Buffer zones in the CNP were established in 1996 with the aim of involving people in the management of park resources for the conservation of biodiversity and the livelihoods of buffer zone communities (31). The buffer zone covers an area of 750 km² and is home to approximately 45,516 houses and 2,60,352 people (32). Buffer zones are primarily designed to create human-wildlife coexistence by providing an ecological and socio-economic buffer for communities, and in these areas, where wildlife and humans share the same landscape in close proximity, their interaction is very obvious (33).

Samples were collected from two deceased adult (>30 years) elephants in Sauraha (Elephant PK and Elephant DG), Chitwan, owned by the government, and used for tourism, grass collection and patrolling. Sauraha is a village at the eastern entrance to Chitwan National Park, near the East Rapti River. These animals had died in captivity and at necropsy showed granulomatous lesions in the lungs compatible with tuberculosis. Tissue samples from lung biopsies (with lesions) were taken from both elephants. The human sputum sample (No. 52G111) collected in Chitwan as part of a previous study (20) was also genotyped, although it was not closely related to the diseased elephants. It was included to understand the evolutionary relationship with the M. tuberculosis strains isolated from elephants, as it was collected from the same region.

The elephant lung tissue samples (PK and DG) and human sputum samples (52G111) were sent to GENETUP (German Nepal Tuberculosis Project microbiology lab) for culture confirmation of M. tuberculosis as described in (34). The heat-killed cultures were then resuspended in 100 μL of PBS for DNA extraction. The DNA extraction was performed using a DNA-SorB kit according to the manufacturer’s instructions (Sacace Biotechnologies, Italy). Briefly, 300 μL of the lysis solution was mixed with 100 μL of the heat lysed M. tuberculosis PBS solution and incubated at 60°C for 5 min. The samples were centrifuged at 15000 g for 5 min to separate the cell debris and DNA in the supernatant. The supernatant was then mixed with 20 μL of the sorbent to bind the DNA to the sorbent. The sorbent was washed twice, dried at 65°C for 5 min and eluted with 25 μL of DNA eluent. Elution was further enhanced by incubating the tubes with DNA eluent at 65°C for 5 min and periodically vortexing. The extracted DNA was then quantified using Qubit (Thermo Fisher Scientific, USA) and M. tuberculosis infection was diagnosed using the specific real-time PCR (35).

Two complementary whole genome sequencing approaches were used in this study: Illumina short reads and MinION long reads. In the case of Illumina short reads, genomic DNA from the heat-killed cultures of M. tuberculosis was purified using AMPure beads (Beckman Coulter, United States) (1:1), quantified using Qubit (Thermo Fisher Scientific, United States) and normalised to 0.2 ng/μL. The normalised samples were tagged and indexed using the Nextera Indexing Kit (Illumina, United States). The indexed library was then purified using (0.8X) AMPure beads and quantified using HS kit (Thermo Fisher Scientific, United States). Each library was analysed using a bioanalyzer and normalised to a final concentration of 2 nM. The three normalised libraries were then pooled in equal volumes of 5 μL to give a final pooled library of 2 nM. The pooled library was spiked with 5% of 2 nM PhiX and denatured with NaOH. The denatured library was diluted to 20 pM, then 10 pM using hybridisation buffer, and the final 10 pM library was loaded onto the Illumina MiSeq platform using the MiSeq Reagent Kit V2 300-cycles (Illumina, United States), with an estimated coverage of approximately100× per sample.

For long-read sequencing in Nanopore MinION, normalised genomic DNA (0.2 ng/μL) from the heat-killed cultures of M. tuberculosis isolated from the two elephants was also used for long-read sequencing in Nanopore MinION, using the Nextera XT Library Preparation Kit (Illumina, USA). The 20 μL of tagged genomic DNA was end-repaired using Ultra II End Prep Buffer and Ultra II End Prep Enzyme (NEB, United Kingdom) according to the manufacturer’s instructions. The end-repaired DNA was then subjected to Nanopore DNA preparation using Ligation sequencing kit (Oxford Nanopore, United Kingdom) with PCR barcoding expansion kit (Oxford Nanopore, United Kingdom). The final library was purified using 0.4X AMPure beads and Long Fragment Buffer (LFB) to obtain the eluate in 15 μL elution buffer. The adapter-ligated library was quantified using the Qubit Fluorometer and approximately 15 fm of the total library was finally loaded onto the MinION for sequencing using 15 μL of sequencing buffer and 10 μL loading beads.

Tb-Profiler version 4.1.0, a stand-alone tool (36), was used for lineage identification and drug resistance assessment of M. tuberculosis isolates from whole genome sequences. The tool accepts raw data from both MinION Nanopore and Illumina MiSeq sequences as well as BAM and FASTA files. It is a command-line tool that can use whole-genome data to predict lineages and identify 21 types of drug resistance (based on small variants and large deletions associated with drug resistance).

Quality control and filtering of Illumina reads was performed using Fastp v0.20 (37), while NanoFilt v2.6.0 (38) was used to filter Nanopore reads. After quality filtering, the Illumina and Nanopore reads were mapped separately to the M. tuberculosis reference genome (GenBank accession NC_000962.3) using BWA v0.7.17 (39). SAMtools v1.15 (40) was used to convert the SAM files of the mapped reads from both Illumina and Nanopore to BAM files, and to sort and index the BAM file reads using the view, sort and index commands. SAMtools merge was then used to merge the index and sorted BAM files into a single BAM file containing both long and short reads for each sample (40).

To performe a comparative phylogenomic analysis of M. tuberculosis isolates, we selected a sequence dataset comprising 94 full-length M. tuberculosis whole genomes, including data from human isolates in Nepal (n = 8), and a genome of M. africanum to use as an outgroup. These genomes were from 33 different countries and the four different geographical lineages previously described from elephants worldwide (lineages 1 and 4) and from humans in Nepal (lineages 1, 2 and 3), and they were retrieved from the NCBI Sequence Read Archive (SRA) database (Supplementary Table S1). The prefetch command in the SRA toolkit v2.11.0 (41) was used to download all whole-genome sequence files of our dataset in SRA format, which were later converted to FASTQ format using the fastq-dump command. Fastp v0.20 was used for quality control and read filtering (37). Read mapping to the reference genome (GenBank accession: NC_000962.3) was performed using BWA v0.7.17 (39). SAMtools v1.15 was used to convert the SAM files of the mapped reads to BAM files, and then to sort and index the BAM file reads (40). The SAMtools mpileup command generated a text pileup output, summarizing the base calls of the aligned reads to the reference sequence. The consensus call was performed using the call command in BCFtools v1.10 (40), which produced a VCF file that was converted to a FASTQ file using the vcf2fq function in the vcfutils.pl script (42). Seqtk v1.3 (42) was used to convert the FASTQ files to FASTA files using the seq command with the masking of bases with quality less than 20.

Variant calling for all the downloaded M. tuberculosis sequences and three M. tuberculosis isolates from this study was performed using PhaME v1.0.4 (43), with the H37Rv genome (GenBank accession: NC_000962.3) as the reference. The PhaME analysis tool was used to remove repeats from each genome in the dataset, and genome alignment to the reference genome was performed using the nucmer2 tool from MUMmer v3.0, as described in the pipeline (43). The output genome alignments were filtered, and SNP tables were generated for each alignment using the delta-filter and show-snps utilities, respectively, both in MUMmer v3.0. The PhaME pipeline then collated the alignment files to produce a core genome alignment and output a multiple alignment FASTA file (~200 kbp) with the variant sites. This was further filtered with SNP-sites to remove non-informative synonymous sites, resulting in a final non-synonymous SNP alignment dataset of 2,596 base pairs.

Using the non-synonymous SNP alignment data, phylogenomic analysis was performed using maximum likelihood and Bayesian inference. For maximum likelihood, the best substitution models and likelihood trees were evaluated using IQTREE v1.6.11 (44) with 500 bootstrap replicates and default parameters. Similarly, for Bayesian inference, the best nucleotide substitution model, TVM, was selected based on the Bayesian Information Criterion score in jModeltest2 v2.1.8 (45, 46), followed by tree reconstruction in MrBayes v3.2.7 (47) with 2,000,000 iterations, discarding 25% as burn-in.

The two heat-killed culture samples from elephants yielded 1.07 ng/μL and 0.35 ng/μL DNA for DG and PK, respectively, allowing confirmation of M. tuberculosis by real-time PCR. The M. tuberculosis isolated from the human sample 52G11 was already a confirmed cultured M. tuberculosis with a DNA concentration of 2.34 ng/μL.

Lineage analysis using TB-Profiler showed that the three isolates belonged to three different M. tuberculosis lineages (Table 1). The elephant isolates DG and PK belong to lineage/sub-lineage 1.2.2.2 and 2.2, respectively, whereas the human isolate 52G111 belongs to lineage/sub-lineage 4.1.2.1. None of the three isolates had drug-resistant variants.