Whole mitochondrial genome sequence and phylogenetic relationships of Williams’s jerboa (Scarturus williamsi) from Turkey

Sampling, gDNA extraction and long-range PCR amplification

A tissue sample (muscle) was obtained from one female individual of S. williamsi collected from Yeşilli Village, Sulakyurt, Kırıkkale (road-killed individual; collection number: 484). Genomic DNA (gDNA) isolation was conducted using the QIAGEN DNeasy^® Blood & Tissue Kit by following the manufacturer’s instructions. The quantity of gDNA was determined by using the Qubit^® 2.0 Fluorometer and Qubit^™ dsDNA BR Kit, and an aliquot of the extract was visualized on a 1% agarose gel. The mitogenome was amplified from gDNA by means of Long-Range PCR with NEB LongAmp^® Taq 2× Master Mix (M0287S, NEB) in two overlapping fragments (~11 kb and ~7.1 kb in length) using two specific primer pairs (Table S1). The Long-Range PCR mixture contained 1× NEB LongAmp^® Taq 2× Master Mix (12.5 µL), 0.6 μM of each primer (1.5 µL of each 10 μM primer) and ~30 ng gDNA (1 µL). The PCR conditions included a pre-denaturation at 94 °C for 1 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C (CrocAL1-2024L-CrocBH1-13002H)–55 °C (ScVu-11712L-LuLu-2503H) for 45 s, and extension at 65 °C for 11 min (CrocAL1-2024L-CrocBH1-13002H)–7.5 min (ScVu-11712L-LuLu-2503H). A final extension step was carried out at 65 °C for 10 min. A negative control (no gDNA) was included in the reaction for to detect possible contamination. A total of 5 µL of the PCR products were run on 1% agarose gel, and 4 µL of the PCR products was used to determine the DNA concentration with the Qubit^™ dsDNA BR Assay Kit according to the manufacturer’s protocol. The standardized PCR products were diluted with ddH₂O, and 1 ng of amplicon (5 µL) was taken to prepare the DNA sequencing library.

Preparation of libraries, sequencing and data analysis

The sequencing library was constructed by using the Nextera XT DNA Library Prep Kit (Illumina, San Diego, CA, USA) and Nextera XT DNA Library Preparation Index Kit v2 Set A (Cat. No: FC-131–2001, Illumina, San Diego, CA, USA) following the protocols of the manufacturer. The quantity of each library was normalized via bead-based normalization. The library loading volume and concentration were determined by using the Qubit^™ dsDNA HS Kit. Library sequencing was performed by using the MiSeq Reagent Kit v2 (500 cycles) (Illumina) and Illumina MiSeq sequencing platform (Genome and Stem Cell Centre, GENKOK, Erciyes University).

Analysis of the obtained 656,266 raw paired-end (2 × 250 bp) reads (average length, 215.7 bp) was carried out using Geneious Prime^® v2019.1.3 (Kearse et al., 2012). We used the BBDuk Trimming Tool in Geneious Prime to filter (remove) the following from the raw-read data: short reads (<50 bp), low-quality bases (Q-score <20), and adapters. The remaining 577,554 reads from S. williamsi were assembled to the O. sibirica mitogenome (National Center for Biotechnology Information (NCBI) accession number: KX058130) using the Geneious Mapper algorithm with the following parameters: Sensitivity; Highest sensitivity/Medium and Fine Tuning; iterate up to 25 times. Contig sequences were obtained, and gene annotations were performed based on the O. sibirica mitogenome. Moreover, gene borders were checked by using MITOS2 (Bernt et al., 2013) and manually curated. The Tandem Repeats Finder Web server (Benson, 1999) was used to determine sequence repeat motifs in the control region (D-loop). AT-skew and GC-skew analyses were calculated using the formulas (A − T)/(A + T) and (G − C)/(G + C), respectively. The graphical mitogenome map of Williams’s jerboa was drawn with the CGView Server (Grant & Stothard, 2008).

Phylogenetic analyses

The S. williamsi mitogenome obtained in this study and the available mitogenome data from taxa belonging to the suborder Myomorpha (n = 43) in the GenBank database (NCBI) were combined into a multiple sequence alignment using the MAFFT algorithm (Katoh et al., 2002). After the control region was removed from the alignments, a total of 16,169 bp remained for phylogenetic analyses. The most suitable substitution model was determined according to the AICc (corrected Akaike Information Criterion) and BIC (Bayesian Information Criterion) criteria by using jModeltest 2.1.10 (Darriba et al., 2012). The Maximum Likelihood (ML) and Bayesian Inference (BI) methods were used to reconstruct phylogenetic trees. The ML tree was generated with MEGA7 (Kumar, Stecher & Tamura, 2016) using the GTR+G+I model and 10,000 bootstrap pseudo-replicates, while the BI tree was constructed with MrBayes v3.2.6 (Ronquist et al., 2012) with four million generations of Markov Chain Monte Carlo iterations. Genetic distance values were calculated in MEGA7 using the Kimura 2-parameter (K2P) nucleotide substitution model.

Mitogenome characterization of Williams’s jerboa

The mitogenome assembly of Williams’s jerboa had a mean coverage of 5.800× (min: 1.780, max: 13.432) and a total length of 16,653 bp (NCBI accession number: MT079957). The mitogenome was similar in structure to those of other Dipodoidea mitogenomes (Yue et al., 2015; Ding et al., 2016; Luo & Liao, 2016; Luo, Ding & Liao, 2016; Reyes et al., 2016, https://www.ncbi.nlm.nih.gov/nuccore/AJ416890.1) and included 13 protein-coding (PCGs), 22 transfer RNA (tRNA), two ribosomal RNA (rRNA) genes, a light chain (L) replication origin (O_L), and a large non-coding control region (D-loop). The ND6 gene, O_L origin and eight tRNA genes (tRNA^Ala, tRNA^Asn, tRNA^Cys, tRNA^Glu, tRNA^Gln, tRNA^Pro, tRNA^Ser(UCN) and tRNA^Tyr) were located on the light chain (L chain), while 12 PCG, 14 tRNA and 2 rRNA genes were located on the heavy chain (H chain) (Fig. 1; Table 1). While the S. williamsi mitogenome (16,653 bp) was longer compared with the mitogenomes of S. concolor (16,492 bp), J. jaculus (16,546 bp) and E. setchuanus (16,573 bp), it was also shorter relative to the mitogenomes of D. sagitta (16,664 bp), O. sibirica (16,685 bp), S. telum (16,696 bp) and E. naso (16,705 bp). The nucleotide composition of the mitogenome was similar to that of typical mammalian mitochondrial genomes, guanine had the lowest frequency (A > T > C > G), with A+T content (59.72%) higher than G+C content (40.28%). The A+T ratio of the D-loop was relatively lower (57.45%) than those of the whole mitogenome (59.72%), PCGs (59.47%), rRNAs (60.40%) and tRNAs (62.88%) (Table S2). Similar to other vertebrate mitogenomes (Ding et al., 2019; Krause et al., 2008; Wada et al., 2010), annotation of the Williams’s jerboa mitogenome showed overlapping regions and intergenic intervals between protein-coding and tRNA genes. The S. williamsi mitogenome had nine overlapping regions with a total length of 66 bp (between 1 and 43 bp) and 13 intergenic intervals with a total length of 42 bp (between 1 and 7 bp) (Table 1).

Protein coding genes with a total length of 11,396 bp constituted 68.5% of the mitogenome and coded a total of 3,787 amino acids, except for the termination codons. While the open reading frame of 10 protein-coding genes started with the ATG codon (76.9%), the ND2 and ND3 genes started with an ATA codon and the ND5 gene started with the ATC codon. Eleven of the PCGs were terminated with full termination codons; TAA (ATP6, ATP8, COX1, COX2, ND3, ND4L and ND6), TAG for ND1, ND2 and ND5 genes and AGA for CYTB. Moreover, COX3 and ND4 were terminated with a T—incomplete termination codon (Table 1). In many metazoan taxa, mitochondrial protein-coding genes may contain an incomplete termination codon (Boore, 1999). This incomplete termination codon can form a complete termination codon through post-transcriptional polyadenylation (Ojala, Montoya & Attardi, 1981). When compared with the available Dipodoidea mitogenomes, the S. williamsi mitogenome showed was a deletion of a serine residue at the C-terminus of the ATP8 gene, similar to observations reported by Yue et al. (2015) for E. setchuanus and S. concolor. However, this deletion was not observed in the mitogenomes of D. sagitta, J. jaculus and E. naso. None of the remaining PCGs showed evidence of insertions or deletions. While 344 amino acids are encoded by ND2 in Muridae species and S. concolor, 345 amino acids are encoded by ND2 in Cricetidae species (Yue et al., 2015). By comparison, 347 amino acids are encoded by ND2 in S. williamsi (This Study) and in a great majority of the remaining rodents (Yue et al., 2015). Relative Synonymous Codon Usage (RSCU) and codon usage analyses of the mitogenome are presented in Fig. 2. Among the amino acids that coded on the mitogenome, those with hydrophobic characteristics were the most frequently encoded (n = 2,379, 62.6%) while those with acidic groups were the least frequently encoded (n = 165, 4.3%). The most frequently used codons were CTA (n = 284, 7.5%), ATA (n = 196, 5.2%) and ATC (n = 172, 4.5%); AGA (n = 1), TAG (n = 3), CGG (n = 3) and TGT (n = 4) were rarely used codons. AGG was not used at all.

The lengths of the tRNAs ranged from 58 bp to 75 bp (total length, 1,506 bp; Table 1). The secondary structures of all tRNAs indicated a typical clover shape; however tRNA^Ser(AGY) did not have the typical dihydropyridine (DHU) arm (Fig. S1). The missing clover structure in tRNA^Ser(AGY) could be activated through structural compensation mechanisms between the other arms (Steinberg & Cedergren, 1994). The two rRNA genes had lengths of 964 bp (12S rRNA) and 1,584 bp (16S rRNA), and were located between tRNA^Phe and tRNA^Leu(UUR); these rRNAs separated from each other by tRNA^Val. The putative origin of replication of the light strand (O_L), one of the non-encoding regions of the mitogenome, was within the WANCY tRNA region. The replication starting sequence at the base of the O_L stem (5′-GCCGG-3′), which is commonly found in Dipodoidea members, can show variations among species and genera (Ding et al., 2016) and was determined in the S. williamsi mitogenome. The control region (D-loop) of the mitogenome, the largest non-encoding region, was located between tRNA^Pro and tRNA^Phe and was 1,194 bp in length. The D-loop, which includes regions a participating in the replication and transcription of mitogenome, is characterized by separate and conserved sequence blocks and divided into three large groups as the extended termination-associated sequence (ETAS), the central conserved domain, and the conserved sequence block (CBS) (Saccone, Pesole & Sbisa, 1991). The S. williamsi D-loop included all three regions (ETAS, 449 bp; central conserved domain, 300 bp; and CSB, 445 bp). It also included twelve 22 bp-long (CATACACACGTACACGCATACG) and one 11 bp-long tandem repeat elements totaling 275 bp in length. A similar pattern, that is, 6–14 bp-long motifs repeated in various numbers was previously observed in the mitogenomes of other Dipodoidea members (Yue et al., 2015; Ding et al., 2016; Luo & Liao, 2016; Luo, Ding & Liao, 2016; Reyes et al., 2016, https://www.ncbi.nlm.nih.gov/nuccore/AJ416890.1). In this context, the longest tandem repeat motif was observed in the S. williamsi mitogenome. This finding supports the view that the tandem repeat sequences found in the control regions are common within Dipodinae species (Yue et al., 2015).

Genetic diversity and phylogenetic analysis

Molecular-based studies on S. williamsi samples are scarce and only partial mitochondrial sequences of this species are available in the NCBI database (Dianat et al., 2013; Kryštufek et al., 2013; Hamidi, Darvish & Matin, 2016; Olgun Karacan et al., 2019, https://www.ncbi.nlm.nih.gov/nuccore/MG255335.1). Kryštufek et al. (2013) established the mitochondrial CYTB (1,104 bp) and 16S rRNA (317 bp) sequences of three Williams’s jerboas from Turkey and found that their samples were grouped with samples of S. williamsi from north western Iran and samples of S. euphraticus from Lebanon. Kryštufek et al. (2013) suggested that the distributional area of S. williamsi extends to southern Lebanon, Turkey, and Iran; the authors also mentioned that S. euphraticus samples were different from those obtained from Syria and formed a different lineage. Kryštufek et al. (2013) emphasized the presence of a cryptic species within S. euphraticus. Olgun Karacan et al. (2019, https://www.ncbi.nlm.nih.gov/nuccore/MG255335.1) recorded the CYTB sequences (888 bp) of S. williamsi samples in the NCBI database. Besides these studies, two other molecular-based studies (Dianat et al., 2013; Hamidi, Darvish & Matin, 2016) used S. williamsi samples from Iran. In these studies, Dianat et al. (2013) used CYTB (894 bp) and COX1 (632 bp) sequences of nine samples from different regions of Iran and suggested that S. williamsi is a species complex depending on the genetic variability within the genus. Hamidi, Darvish & Matin (2016) reported the first distributional record of S. williamsi (S. cf. williamsi) from northeastern Iran, in the Kopet-Dag Mountains based on the CYTB (695 bp) sequence obtained from a single sample.

The S. williamsi mitogenome obtained in the present study and the mitogenomes of 43 Myomorpha species were combined into a multiple sequence alignment without the D-loop region and analyzed. The genetic distance between the mitogenomes of S. williamsi and O. sibirica was calculated as 0.176 (17.6%) according to the K2P nucleotide substitution model. The average genetic distance within Dipodinae is 0.178 (17.8%). Mean genetic distances (between-group mean distance) among subfamily Allactaginae and subfamilies Dipodinae, Euchoreutinae, Sicistinae, and Zapodinae were calculated as 0.232 (23.2%), 0.237 (23.7%), 0.304 (30.4%) and 0.268 (26.8%), respectively. BI (Fig. 3) and ML (Fig. S2) analyses produced phylogenetic trees with similar topologies that could be used to delineate subfamilies. All subfamilies within Muroidea formed a well-resolved monophyletic group with high nodal support values. Dipodoidea was located as a basal clade within Myomorpha. The superfamily Dipodoidea formed a sister group with Muroidea, and this relationship has been well demonstrated in previous studies (Holden & Musser, 2005; Steppan, Adkins & Anderson, 2004; Zhang et al., 2013; Michaux & Shenbrot, 2017; Steppan & Schenk, 2017). The results of the present study are in agreement with the results of previous studies. Subfamily Sicistinae, which includes S. concolor, formed the most basal branch of Dipodoidea, consistent with previous studies (Blanga-Kanfi et al., 2009; Churakov et al., 2010; Zhang et al., 2013; Yue et al., 2015). S. williamsi and O. sibirica were found within the subfamily Allactaginae (nodal support: BI, 1.0; ML, 100%), and the species were closely associated with subfamily Zapodinae (E. setchuanus) and subfamily Euchoreutinae (E. naso). D. sagitta, S. telum, and J. jaculus were clustered together (nodal support: BI, 1.0; ML, 99%) and formed subfamily Dipodinae. The species of subfamily Allactaginae were more basal than the species of Euchoreutinae and Dipodinae (Fig. 3). The results of the current study are largely similar to those of previous studies (Yue et al., 2015; Ding et al., 2016; Luo & Liao, 2016; Luo, Ding & Liao, 2016; Steppan & Schenk, 2017) and support the view that both species should be classified within different genera (Scarturus and Orientallactaga) (Michaux & Shenbrot, 2017).