Mitochondrial genomes assembled from non-invasive eDNA metagenomic scat samples in the endangered Amur tiger Panthera tigris altaica

Panthera tigris altaica scat samples

The raw sequence data used to assemble the mitochondrial genome of P. tigris altaica from eDNA was generated by He et al. (2018) and employed to describe the effect of two different diets on the gut microbiome of five different captive cubs. In brief, fresh fecal samples from five tiger cubs were collected at a single time point from the Heilongjiang Siberian Tiger Park, Heilongjiang Province, China, Next, genomic DNA (gDNA) was extracted from the fecal samples using the E.Z.N.A.® Stool DNA Kit (D4015–02; Omega, Inc., Norwalk, CT, USA) according to the manufacturer’s instructions. Then, paired-end (PE) shotgun libraries were constructed using a TruSeq Nano DNA LT Library Preparation Kit (FC-121-4001) following the manufacturer’s instructions. Finally, libraries were sequenced in an Illumina sequencer (Illumina, San Diego, CA, USA) using a 2 × 150 cycle. For further details, see He et al. (2018). The five different sets of sequences were retrieved from GenBank (SRA Accession numbers: SRR7429862–SRR7429866) (Fig. 1). The amount of data available in FASTQ format varied between 42,819,376 (SRR7429862) and 78,803,225 (SRR7429864) PE reads per sample (Table 1). All of the available reads generated for each faecal sample was used for mitochondrial genome assembly.

Mitochondrial genome assembly of Panthera tigris altaica from metagenomic faecal samples

Assembling of the mitochondrial genome from the different scat samples was attempted using the target-restricted-assembly pipeline GetOrganelle v1.2.3 (Jin et al., 2020). GetOrganelle uses a seed-and-extend algorithm that assemble organelles, including mitochondrial genomes, from whole genome sequencing (WGS), including metagenomic datasets, starting from a related or distant single ‘seed’ sequence (Jin et al., 2020) (Fig. 1). For each metagenomics sample, GetOrganelle was run twice. In the first run, the complete mitochondrial genome of P. tigris altaica (MH893763, from China) retrieved from GenBank was used as the ‘seed’ while in the second run, only a single gene fragment (cox1) from the same mitochondrial genome of P. tigris altaica retrieved from GenBank was used as a seed. These two runs per sample permitted us to explore the effect of seed length on the quality (i.e., completeness) and coverage of the retrieved mitochondrial genome assemblies. All the different runs used k-mer sizes of 21, 55, 85, and 115. Reads were not quality trimmed prior to assembly using the GetOrganelle pipeline following the developer’s guidelines (Jin et al., 2020). The software Bandage (Wick et al., 2015) was used to visualize the assembly graph generated by the pipeline above and any assembled contigs were compared to the nucleotide non-redundant database in NCBI’s GenBank in order to determine if these contigs belonged to the mitochondrial genome of the target species. I predicted that a circularized (or not) sequence ~17 kpb in length would be observed among the contigs if the pipeline above successfully assembled the mitochondrial genome of P. tigris altaica.

Annotation and analysis of mitochondrial genome or contigs

Complete or partially assembled mitochondrial genomes were first annotated in-silico using the web servers MITOS (http://mitos.bioinf.uni-leipzig.de) and MITOS2 (http://mitos2.bioinf.uni-leipzig.de) (Bernt et al., 2013) with the vertebrate genetic code (code 2). Manual curation of the in-silico annotations, including start + stop codons corrections, were performed using the web server ExPASy (https://web.expasy.org/) and the software MEGA7 (Kumar, Stecher & Tamura, 2016). Genome visualization was performed in OGDRAW–Draw Organelle Genome Maps (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) (Lohse, Drechsel & Bock, 2007) (Fig. 1).

Phylogenetic position of mitochondrial genomes retrieved from feces

The phylogenetic position of the mitochondrial genomes assembled from scat eDNA was examined among other representatives belonging to the genus Panthera. The five newly assembled and annotated mitochondrial genomes of P. tigris altaica plus those (complete or nearly complete) of other 67 specimens of Panthera, all of them available in the GenBank database (consulted 11 07 2021), were used for the phylogenetic analysis conducted with the program MitoPhAST (Tan et al., 2015). One mitochondrial genome belonging to the cloud leopard Neofelis nebulosi was used as an outgroup. In MitoPhAST, all 13 protein-coding gene nucleotide sequences from the species used in this analysis were retrieved from GenBank files, translated to amino acids, and aligned using Clustal Omega (Sievers et al., 2011). Next, poorly aligned regions with large numbers of indels were removed with trimAl (Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009). Partitions and best fitting models of sequence evolution for each partition were selected with ProtTest (Abascal, Zardoya & Posada, 2005). The concatenated and partitioned protein-coding gene amino acid alignments were used to perform a maximum likelihood phylogenetic analysis in the software IQ-TREE (Nguyen et al., 2015). The robustness of the ML tree topology was assessed by bootstrapping the observed data 1,000 times. In this study, we used amino acids instead of nucleotide characters for phylogenetic inference considering that the former markers have a greater ratio of phylogenetic information to noise, compared to nucleotides, for resolving basal and derived nodes in a phylogenetic tree (Page & Holmes, 2009).

Assembly of mitochondrial genomes from Panthera trigis altaica field scats

The program GetOrganelle assembled and circularized the mitochondrial genome of P. tigris altaica in four of the five tested scat eDNA metagenomics samples. Furthermore, the two strategies (full mitochondrial genome vs only cox1 seed) employed to assemble the mitochondrial genome of P. tigris altaica resulted in identical sequences in each scat eDNA sample. The length and coverage of these circularized mitochondrial genomes varied between 16,862 bp (coverage = 322.7x per nucleotide) and 17,095 bp (17.6x) in samples SRR7429865 and SRR7429863, respectively. In one sample (SRR7429866), Getorganelle assembled but did not circularize the mitochondrial genome of P. tigris altaica. The length and coverage of this assembled nearly complete mitochondrial genome was 16,964 bp and 101x, respectively (Table 1). Importantly, this assembled mitochondrial genome was only 16 bp shorter compared to the reference genome used for its assembly (MH893763).

Annotation of mitochondrial genomes assembled from scats of Panthera trigis altaica

All mitochondrial genome assemblies, either complete or partial, of Panthera tigris altaica retrieved from scat eDNA samples were composed of 13 protein-coding genes, two ribosomal RNA genes (rrnS (12S ribosomal RNA) and rrnL (16S ribosomal RNA)), 22 transfer RNA (tRNA) genes, and a non-coding region of variable length (Fig. 2, Table S1). Most of the protein-coding and tRNA genes were encoded on the L-strand, while only a single protein-coding gene (nad6) and eight tRNA genes (5′ to 3′ order: trnQ, trnA, trnN, trnC, trnY, trnS2, trnE, trnP) were encoded in the H-strand (Fig. 2, Table S1). The gene order observed in these complete or nearly complete mitochondrial genomes was identical to that reported before in Panthera tigris altaica as well as the whole genus Panthera (Wei, Wu & Jiang, 2009; Kitpipit & Linacre, 2012; Dou et al., 2016). A detailed characterization of the features present in all complete mitochondrial genomes is provided below in order to explore the accuracy of these assemblies retrieved from scat eDNA samples.

Nucleotide usage (overall % base composition) of the mitochondrial genomes assembled from the different scat samples (SRR7429862–SRR7429866) is presented in Table 2. All assembled mitochondrial genomes were compact with only a few intergenic spaces and overlaps among gene junctions (Table S1). In all circularized mitochondrial genomes, a single relatively long intergenic space varying in length between 1,413 bp in sample SRR7429865 to 1,531 bp in sample SRR7429863 was assumed to be the Control Region/D-loop that is involved in replication of this genome (Bernt et al., 2013).

In all the complete mitochondrial genome assemblies retrieved from scat eDNA samples, all of the 13 PCGs exhibited conventional vertebrate mitochondrial start codons; ATA, ATC, and ATG (Table 1). Invariably, eleven PCGs terminated with a complete and conventional stop codon (TAA or TAG) while the genes nad4 and cox3 terminated with the incomplete stop codons T and TA, respectively. Truncated stop codons are thought to be completed via post-transcriptional poly-adenylation (Ojala, Montoya & Attardi, 1981; Rorbach & Minczuk, 2012 and references therein).

In the retrieved mitochondrial genome assemblies of P. tigris altaica, the most frequently used codons found in the PCGs were AT-rich, and included CTA (Leu, n = 241–242 in 3 samples, 270 in one sample), ATA (Met, n = 176 in all samples), and ATC (Ile, n = 170–171 in all samples). In turn, discounting stop codons, least frequently used codons were GC-rich, and included CGG (Arg, n = 3 in all samples), CCG (Gln, n = 7–8 in all samples), GCG (n = 7 in all samples), and CGT (n = 7 in all samples). This aforementioned codon usage pattern is akin to that reported before for protein-coding genes in the mitochondrial genome of P. tigris altaica and other congeneric species whose codon usage pattern has been examined (Wei, Wu & Jiang, 2009; Kitpipit & Linacre, 2012; Dou et al., 2016).

In each one of the mitochondrial genome assemblies, tRNA genes ranged in length from 59 (trnS1) to 75 bp (trnL2) (Table S1). All but one tRNA gene (trnS1) exhibited a standard ‘cloverleaf’ secondary structure. In the trnS1 gene of all assemblies, the stem and loop of the pseudouridine arm (T-arm) was missing. This truncated trnS1 represents a conserved trait among vertebrates and eumetazoans (Bernt et al., 2013) and has also been found in all other representatives of the family Felidae in which the secondary structure of tRNA genes have been examined (Wei, Wu & Jiang, 2009; Kitpipit & Linacre, 2012; Dou et al., 2016). Whether or not truncated tRNAs interact with other molecular factors while decoding mRNA into protein remains to be investigated in species belonging to the family Felidae as well as in many other vertebrates and invertebrates (Watanabe, Suematsu & Ohtsuki, 2014).

In all assemblies, the length of the rrnS (12S) and rrnL (16S) genes was 962 and 1,574 bp, respectively. These two ribosomal genes were located close to each other; the rrnS gene was located between the trnF and trnV genes while the rrnL gene was located between the trnV and trnL2 genes. The two genes exhibited an AT-skew. Nucleotide usage of the rrnS and rrnL genes are presented in Table 2.

Overall, the annotation and detailed characterization above indicates that whole or nearly complete mitochondrial genome assemblies retrieved from eDNA are accurate.

Phylogenetic position of mitochondrial genomes retrieved from feces

The ML phylogenetic tree (terminals: 73, amino acid characters: 3,792, informative sites: 286) indicated monophyly of the genus Panthera with high support values (bv = 100) (Fig. 3). Within the genus Panthera, all representatives of P. onca (n = 3 terminals) clustered into a single well supported clade (bv = 97) that was sister to all other specimens used in the analysis. In the tree, the five mitochondrial genomes assembled from P. tigris altaica scat eDNA together with other 34 complete or nearly complete mitogenomes of the same species P. trigris retrieved from Genbank clustered into a very well supported clade (bootstrap value [bc] = 98). Within this monophyletic P. trigris clade, the different specimens did not segregate according to subspecies (Fig. 3). This is not necessarily unexpected considering that ancient hybridization events among species and historic admixture and gene flow between isolated populations are apparently widespread in felids, including within the genus Panthera (de Manuel et al., 2020). The other well-supported clade (bv = 99) was composed of all representatives belonging to P. leo and the extinct cave lion P. spelea. Within this P. leo + P. spelea clade, the different specimens did not segregate according to species, in contrast to a previous study that relied on a smaller number of mitochondrial genomes to reveal the phylogenetic position of the cave lion (Barnett et al., 2016). Two moderately-supported clades included those formed by all specimens belonging to P. uncia (bv = 59) and P. pardus (bv = 69). The latter clade was sister to the well-supported clade composed of P. leo + P. spelea (Fig. 2). Overall, this analysis demonstrates that the mitochondrial genomes assembled from eDNA scat samples can and do reliably identify specimens from which scat samples were obtained as belonging to P. tigris. Furthermore, our analysis shows that mitochondrial genomes assembled from scat samples can be distinguished from closely related congeners.

Our mitogenomic Panthera phylogeny differs from the most recent Panthera phylogeny estimated from genome-wide data (Figueiró et al., 2017). In contrast to that observed in our analysis, Figueiró et al. (2017) placed P. tigris as sister to P. uncia, as well as P. leo and the jaguar P. onca sister to P. pardus. The contrast in topologies based on mitochondrial genomes alone (this study) vs. genome-wide markers might be due to ancient hybridization events and historic admixture (de Manuel et al., 2020).