Characterization of the mitochondrial genomes of two toads, Anaxyrus americanus (Anura: Bufonidae) and Bufotes pewzowi (Anura: Bufonidae), with phylogenetic and selection pressure analyses

Sample Collection and DNA extraction

The specimens of A. americanus and Bufotes pewzowi were collected in Ottawa, Ontario, Canada (45°25.38′N, 75°43.11′W) and Aksu, Xinjiang, China (41°10.54′N, 80°16.81′E), respectively. Both were identified by J.Y. Zhang. Toe-clip samples of both species were stored at −80 °C in the Institute of Ecology, Zhejiang Normal University. Our experimental procedures complied with the current regulations on animal welfare and research in China and Canada. The Animal Research Ethics Committees of Zhejiang Normal University and Carleton University approved the experimental design (No.1082196).

PCR amplification, and sequencing

Total DNA was extracted from the toe-clips samples of A. americanus or Bufotes pewzowi using an Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech Company, Shanghai, China). The 11 universal primers for standard polymerase chain reaction (PCR) amplification of mitogenomes were slightly modified according to Yu et al. (2015) and Zhang et al. (2013) and 15 specific primers were designed based on the sequenced fragments from universal primers using Primer Premier 5.0 (PREMIER Biosoft International, CA, USA) (Tables S1 and S2). All PCR amplifications were carried out in a 50 µL reaction mixture and the procedures were performed using an Eppendorf Thermal Cycler (Mastercycle® nexus GSX1, Hamburg, Germany). We used both standard PCR and Long-PCR methods with TaKaRa Ex-Taq and TaKaRa LA-Taq Kits (Takara Biomedical, Dalian, China). These two methods were slightly modified from Yu et al. (2015) and Zhou et al. (2009). All PCR products were detected by electrophoresis on 1% agarose gels, and sequences were obtained in an automated DNA sequencer (ABI 3730) by Sangon Biotech Company (Shanghai, China).

Mitogenome annotation and sequence analyses

Sequences were checked and assembled using SeqMan (Lasergene version 5.0) (Burland, 2000). The 22 tRNAs were identified by their cloverleaf secondary structure using tRNAscan SE 1.21 (Lowe, 1997) (http://lowelab.ucsc.edu/tRNAscan-SE/) or determined by comparison with the available tRNA genes of closely related anurans downloaded from GenBank. Location of the 13 PCGs and 2 rRNA genes were determined by comparison with homologous sequences of mtDNA from other Bufonidae species using ClustalW (Thompson, Higgins & Gibson, 1994) and then PCGs were checked and translated to amino acids using the vertebrate mitogenome code by Mega 5.0 (Tamura et al., 2011). The mitogenome maps of A. americanus and Bufotes pewzowi were constructed using GenomeVx (http://wolfe.gen.tcd.ie/GenomeVx/) (Conant & Wolfe, 2008). The A+T and C+G content values, codon usage and relative synonymous codon usage (RSCU) of protein-coding genes were calculated using Mega 5.0 (Tamura et al., 2011). Nucleotide sequence skewness was calculated according to the following formulae: AT-skew = (A −T)/(A +T) and GC-skew = (G −C)/(G +C) (Perna & Kocher, 1995).

Phylogenetic analyses

To confirm the phylogenetic relationships among Bufonidae, 19 sequences of complete or partial mitochondrial genomes were used. The data set was inclusive of the ingroups of the 2 species from this study, 13 other species from Bufonidae (Cao et al., 2006; Dong & Yang, 2016; Igawa et al., 2008; Jacob Machado et al., 2018; Jiang et al., 2017b; Machado, Lyra & Grant, 2016; Wang et al., 2013; Yang et al., 2016; Zhang et al., 2005; Zhang et al., 2013) and the outgroups of 4 species from Mannophryne and Dendrobatidae (Lyra et al., 2017; Zhang et al., 2013). Accession numbers of all mitogenomes are summarized in Table 1. The amino acid and nucleotide sequences of the 13 protein-coding genes from all 19 species were employed to construct BI and ML phylogenetic trees according to the methods of Zhang et al. (2018) and Zhou et al. (2009). All of the 13 PCGs were aligned using Clustal W in Mega 5.0 (Tamura et al., 2011) and were analyzed with Gblocks 0.91b (Castresana, 2000) using default settings to select conserved regions. The best partition scheme and evolutionary model were determined with the PartitionFinder v.1.1.1 program (Lanfear et al., 2012), using the Bayesian Information Criterion (BIC) (Schwarz, 1978). For the ML analysis, we used the RAxML program (Stamatakis, 2014) under the GTRGAMMAI model with 1,000 bootstrap replications. For BI analysis, we used MrBayes 3.1.2 (Huelsenbeck & Ronquist, 2001) under the GTR+I+G model. Markov Chain Monte Carlo (MCMC) was run with four chains for 10 million generations, with sampling every 1,000 generations. The first 25% of generations were removed as burn-in, which was decided by checking convergences of -log likelihood (-lnL). After the average standard deviation of split frequencies in Bayesian was below 0.01, we judged that the Bayesian analysis had reached sufficient convergence.

Analysis of positive selection

The program EasyCodeML (Gao et al., 2019) was used to analyze the selective pressure on mitogenomes; this is an interactive visual tool for detecting selection in a molecular evolutionary analysis based on CodeML (Yang, 2007). The ω ratio is the rate of nonsynonymous (dN) versus synonymous (dS) substitution (dN/dS) and can indicate natural selection acting on the proteins. All of the concatenated 13 PCGs (Table S3) were used in the analysis and the values for ω ratio >1, =1 or <1 indicate positive selection, neutral evolution or negative selection, respectively (Yang, 2007). To investigate whether positive selection occurred on specific branches, branch models were run under the one-ratio model (M0) or the two-ratio model with A. americanus or Bufotes pewzowi as the foreground branch, respectively. M0 assumes that all branches have the same ω ratio values whereas the two-ratio model assumes one ω ratio value for the branches of interest and the other for the background branches. Because different topologies of trees will affect the results, both of the phylogenetic trees structured in BI and ML were used in the analyses. In addition, a likelihood ratio test (LRT) was performed to assess the significant difference between the results of the M0 and the two-ratio model (Yu et al., 2011). Variation happens mostly in several base pairings and may affect a few sites in some lineages (Yang, 2007; Yu et al., 2011). Consequently, the site model was applied to detect the potential selection among sites and allow for different ω ratios in different sites, codons or amino acids (Yang, 2007). Seven useful codon substitution models were taken into account in the calculations, including M0 (one ratio), M1a (Nearly Neutral), M2a (Positive Selection), M7 (β), M8 (β & ω) and M8a (β & ω = 1). We also used LRTs to assess these models and Bayes Empirical Bayes (BEB) to evaluate the posterior probability of positive selection sites. The three-dimensional (3D) structures of the amino acid positive selections in the ND6 protein were formed using SWISS-MODEL Workspace (Waterhouse et al., 2018).

Mitogenome organization and arrangement

The lengths of the complete A. americanus and Bufotes pewzowi mitogenomes are 17,328 base pairs (bp) and 17,551 bp, respectively (Table 2). Both mitogenomes are circular and contain the typical 37 genes (Tables 3 and 4). Most of the genes are coded on the H-strand, except for 8 tRNA genes and the ND6 gene on the L-strand. Gene structures are detailed in Figs. 1 and 2. In addition, the gene order and composition are identical with that of other mitogenomes of Bufonidae (Dong & Yang, 2016; Jiang et al., 2017b; Zhang et al., 2016). The different lengths of the mitogenomes are primarily caused by the different sizes of intergenic nucleotides (IGNs) (Tables 3 and 4), particularly the length of the CRs. The overall base composition, A+T and G+C content, as well as AT and GC skew of the A. americanus and Bufotes pewzowi genomes are listed in Table 2; these data show that A. americanus has the highest A+T content (62.4%) and a strong A+T bias. The H-strand of both mitogenomes showed a negative AT-skew and GC-skew, which is expected for most vertebrates (Fonseca, Froufe & Harris, 2006; Hao, Ping & Zhang, 2016; Zhang et al., 2019). Previous studies showed that the asymmetry of the nucleotides resulted primarily from mutations affecting the H-strand during its single-stranded state (Sahyoun et al., 2014).

Protein-coding genes and codon usages

All the typical 13 PCGs ranged from 165 bp (ATP8) to 1,789 bp (ND5) in both A. americanus and Bufotes pewzowi. In addition, the total size of the 13 PCGs in A. americanus and Bufotes pewzowi were identical (11,290 bp). Mitochondrial PCGs have no introns, but have several overlapping nucleotides with the adjacent gene. There are 5 and 6 reading frame overlaps in the mitogenome of A. americanus and Bufotes pewzowi, respectively (Tables 3 and 4). All PCGs begin with ATA or ATG as the start codons, except ND1 that begins with TTG. Although TTG is an uncommon start codon among PCGs, it is often found in amphibians (Wang, Cao & Li, 2017; Zhang et al., 2005). The stop codons for A. americanus and Bufotes pewzowi are mostly complete TAA, AGG, AGA, and TAG codons with some incomplete T codons. The functionality of these latter is believed to be restored by post-transcriptional polyadenylation (Liu, Wang & Su, 2005; Ojala, Montoya & Attardi, 1981). Furthermore, a comparative analysis indicated that the mitogenome of A. americanus had the highest A+T content (61.9%) (Dong & Yang, 2016; Jacob Machado et al., 2018). We also analyzed the relative synonymous codon usage (RSCU) of the A. americanus and Bufotes pewzowi mitogenomes, excluding stop codons (Table 5; Figs. 3 and 4). The results show that A and T at the third codon position are slightly overused when compared to the synonymous codons of C and G. Moreover, we found that the most frequently used amino acids in these two mitogenomes are Leu (CUN), Ile and Ala, with those encoding Cys and Ser (AGY) being rare.

Ribosomal and transfer RNA genes

We found that the length of rRNAs in Bufotes pewzowi (2,543 bp) is longer than other Bufonidae species (Cao et al., 2006; Igawa et al., 2008; Jiang et al., 2017b; Wang et al., 2013; Yang et al., 2016). The complete mitogenome of A. americanus or Bufotes pewzowi contains 22 typical tRNA genes, which ranged from 64 bp to 73 bp in length. In the H-strand, A. americanus showed a positive AT-skew and GC-skew of the tRNAs, and Bufotes pewzowi also showed a positive AT-skew and GC-skew (Table 2). Furthermore, all of the tRNA genes displayed the potential to fold into the typical cloverleaf secondary structure, excluding the trnC and trnS (AGY) (Figs. S1 and S2). In both species, the trnS (AGY) has an extra unpaired adenine (A) in the anticodon arm and the trnC has lost the DHU loop. This phenomenon is also found in other Bufonidae species (shown in Fig. S3) (Cao et al., 2006; Igawa et al., 2008; Zhang et al., 2005; Zhang et al., 2013). The trnS (AGY) in A. americanus had also lost the DHU arm, whereas in Bufotes pewzowi lost only the DHU loop. Hence, the unusual secondary structure of trnS (AGY) is in line with the molecular trend observed in metazoa (Wolstenholme, 1992). We also found 13 unmatched base pairs in A. americanus and 15 in Bufotes pewzowi. The putative origin of L-strand replications (O_L) are 29 bp and 28 bp long in A. americanus and Bufotes pewzowi, respectively. Both of these have the potential to fold into the characteristic stem and loop structure (Jiang et al., 2017a; Sahyoun et al., 2014) that has been demonstrated by Hixson & Brown (1986) to be involved in the transition from RNA to DNA synthesis.

Control region and intergenic regions

The control regions (CRs) in A. americanus and Bufotes pewzowi have lengths of 1,916 bp and 2,129 bp, respectively. The A+T content value (70.3%) of A. americanus is the highest known to date among Bufonidae species. In addition, A. americanus showed a negative AT-skew and GC-skew in its CR, whereas Bufotes pewzowi showed a positive AT-skew and negative GC-skew. Guanine (G) was the scarcest nucleotide at the third codon position of the H-strand due to a strong bias against guanine usage in the A. americanus mitogenome, which is common for mitogenome strands of vertebrates (Kan et al., 2010). The CRs can be divided into three regions that depend on the distribution of the variable nucleotide positions and differential frequencies of the nucleotides, and contain repeat regions at both 5′ and 3′-sides. We found 3 tandem repeats in A. americanus with consensus sizes of 162 bp, 21 bp and 2 bp. In Bufotes pewzowi there are only two tandem repeats whose consensus sizes were 104 bp and 91 bp. The intergenic spacers (IGNs) have a variable length within the mitogenomes of A. americanus (Table 3) and Bufotes pewzowi (Table 4). All lengths of IGNs in both species are smaller than 5 bp, except for one that was 37 bp in A. americanus and 39 bp in B. pewzowi, both located between the trnS (AGY) and ND5 genes. Comparative analysis showed that the longest IGN between A. americanus and Bufotes pewzowi had a high similarity (76.9%). In many vertebrates, there is usually a trnL (CUN) gene following the trnS (AGY) gene (Cheng et al., 2018; Lin et al., 2014; Ni et al., 2016; Ye et al., 2016; Yu et al., 2012; Yu, Zhang & Zheng, 2012). However, when the trnL (CUN) gene is translocated, an IGN is left in its original position (Zhou et al., 2009). Thus, Cao et al. (2006) speculated that this might reflect the evolution of the mitogenome arrangement in Anura.

Phylogenetic analyses

Phylogenetic relationships based on the nucleotide sequences of the 13 PCGs were obtained with BI and ML analyses (Fig. 5). Our BI analysis was relatively robust and provided resolution with high clade posterior probabilities, whereas the ML results showed some modest clade frequencies. Furthermore, the phylogenetic relationships deduced from BI and ML analyses showed somewhat different topologies. In the BI and ML analyses, the monophyly of Bufonidae was well recovered. In the BI analysis, Bufo japonicus was a sister group to the clade of (Bufo gargarizan DQ275350 + (Bufo tibetanus JX878885 + (Bufo gargarizan KM587710 + Bufo gargarizan KU321581))). However, in the ML analysis, Bufo japonicus was a sister group to Bufo stejnegeri and then the group of (Bufo japonicus + Bufo stejnegeri) was a sister group to the cluster of (Bufo gargarizan DQ275350 + (Bufo tibetanus JX878885 + (Bufo gargarizan KM587710 + Bufo gargarizan KU321581))). Coincidentally, Dong & Yang (2016) also reported that Bufo japonicus was a sister group to Bufo stejnegeri and the monophyly of family Bufonidae and genus Bufo were corroborated. Frost et al. (2006) and Van Bocxlaer et al. (2009) corroborated the monophyly of Bufonidae as well as the results of Portik & Papenfuss (2015). Our phylogeny within Bufonidae is generally similar to the results of Pramuk et al. (2008), Ron et al. (2015) and Dong & Yang (2016). The ML and BI analyses also confirmed the monophyly of Bufo, which has also been demonstrated by other researchers (Graybeal, 1997; Dong & Yang, 2016). By contrast, our results were different from Brandvain et al. (2014) and Jiang et al. (2018), who reported that Bufo was a paraphyletic group, a result that was also deduced by Pyron & Wiens (2011). However, we found that Dong & Yang (2016) confirmed the monophyly of Bufo using the complete mitogenomes, whereas Brandvain et al. (2014), Pyron & Wiens (2011) and Jiang et al. (2018) confirmed a paraphyletic relationship using partial mitochondrial genes (12S, 16S and CYTB) and several nuclear genes. Hence, the different results from these studies are likely caused by the different data sets and/or methods applied.

The family Bufonidae was divided into 3 major groups: (Leptophryne + (Bufotes + (Duttaphrynus + (Strauchbufo + (genus Bufo))))), (A. americanus + Rhinella acrolopha), and (Melanophryniscus simplex + M. moreirae). Among them, the clade of A. americanus was the sister group to R. acrolopha whereas Bufotes pewzowi was a sister group to the clade of (Duttaphrynus + (Strauchbufo + (genus Bufo))). Genus Anaxyrus was a sister group to genus Rhinella and then the clade of (A. americanus + R. acrolopha) was the sister group to (Leptophryne + (Bufotes + (Duttaphrynus + (Strauchbufo + (Bufo))))). Portik & Papenfuss (2015) also found that genus Anaxyrus was the sister group to genus Rhinella. Our phylogenetic tree showed that the clade of (M. simplex + M. moreirae) (BI: 1, ML: 100%) was recovered as the sister group to the other Bufonidae species. We also found a similar result in the research of Portik & Papenfuss (2015). The South American genera include Melanophryniscus, Nannophryne, as well as the group of Rhinella that is endemic to South America and is distributed from the southern tip of Texas to as far south as Argentina (Vallinoto et al., 2010). This phenomenon was consistent with the results from other researchers and implies that Bufonidae may have originated from South America. Coincidentally, nearly all recent phylogenetic analyses have asserted a South American origin for Bufonidae (Pérez-Ben, Gómez & Báez, 2019; Pramuk et al., 2008). In addition, discrepancies in the results between our work and prior studies were likely to be caused by the different composition of the data chosen (Gao et al., 2018). Therefore, more molecular data are needed to better comprehend the phylogenetic relationship of the family, especially from taxa that, in other studies, led to a conclusion that Bufo was paraphyletic.

Analysis of positive selection in 13 protein-coding genes

The results of positive selection are shown in Tables 6 and 7. In the analyses of the branch model, we reached a similar conclusion no matter which toad species was set as the foreground branch. The ω ratio value in the M0 model was 0.03375 when using the BI tree and 0.03373 using the ML tree, with ω ratio values all smaller than 1. This means that these two branches are under negative selection. When we set A. americanus as the foreground branch, the LRT of the comparison (M0 vs. the two-ratio model) was highly significant (p < 0.01), whereas when Bufotes pewzowi was set as the foreground branch, the LRT value was greater than 0.05. The site model was used to detect positive selection sites and we got similar conclusions regardless of which tree-building method was used. The LRT of M7-M8 comparison showed high significance (P < 0.01) and two amino acid positions were found have BEB values >0.90 (positions 3,590 and 3,594 in the concatenated amino acids of the 13 PCGs). Amino acid residues 3,590 and 3,594 correspond to amino acid positions 115 and 119 in the ND6 protein, respectively (Table S3). Residue glycine (G, 115) in the ND6 protein can be found in A. americanus, Bufo gargarizan, Bufo japonicus, Bufo stejnegeri, Bufo tibetanus and S. raddei, which are distributed across northern regions of the earth at similar latitudes (33°N to 44°N), and residues serine as well as glutamic acid can also be found in this position in other Bufonidae. Residue aspartic acid (D, 119) in the ND6 protein can only be found in A. americanus whereas most of the other Bufonidae species contain glycine at this site. Amino acid positions 115 and 119 (Fig. S4) are in a part of the ND6 protein that protrudes outside the mitochondrial inner membrane into the intermembrane space (as determined using the TMHMM Server v. 2.0; http://www.cbs.dtu.dk/services/TMHMM/) (Möller, Croning & Apweiler, 2001). In this study, A. americanus was obtained from the highest latitude and contained both of the two positive selection sites. At 115 and 119 positions of ND6 protein, aspartic acid has a negatively charged polar side chain, whereas glycine has an uncharged side chain with no polar group (Shiraishi & Kuwabara, 1970). The amino acid changes in the ND6 protein may have a role in modulating Complex I redox potential and ROS production (Zhuang & Cheng, 2010).

Consequently, the above evidence possibly implies that the mitogenomes of these northern species are under natural selection. The changes in the relative mass and electrical charges of amino acids that they encode may be related to adaptive evolution to low temperature environments in A. americanus. As the first enzyme complex of the respiratory chain, mitochondrial complex I consists of 45 subunits, its seven hydrophobic subunits being encoded by the mitogenome (ND1-6 and 4L) (Formosa et al., 2018; Zhou et al., 2014). Thus, the adaptive changes in several ND4 and ND6 gene sites may affect the proton-pumping process and metabolic performance (Yu et al., 2011). Sun et al. (2018a) found positive selection on ND4 from Tetranychus truncatus during adaptation to low temperature. Zhuang & Cheng (2010) proposed that the modification of ND6 gene probably improved complex I subunit interactions at low temperatures after analyzing positive selection results. Lamb et al. (2018) also found positive selection of the ND6 gene when tested for evolution caused by climate-linked selection. Furthermore, residue aspartic acid (D, 119) can also be found in the same position in the ND6 proteins of other anurans, such as in the family Ranidae (Rana dybowskii, R. cf. chensinensis and R. huanrenensis), Leptodactylidae (Leptodactylus melanonotus) and Microhylidae (Phrynomantis microps). We observe that the three frogs in the family Ranidae are mainly distributed in northern China (Li, Lei & Fu, 2014; Dong, Zhou & Yang, 2015), where temperatures are relatively low. L. melanonotus and P. microps are distributed near the equator, where temperatures are relatively high. More research is needed to reach definitive conclusions about low or high temperature positive selections in mitochondrial genes (e.g., Banguera-Hinestroza et al., 2018; Ben Slimen et al., 2018; Sun et al., 2018b; Zhou et al., 2014). According to this evidence, we can hypothesize that residue aspartic acid (D, 119) in the ND6 protein of A. americanus may be related to temperature adjustment and in different anurans may be under different climate-linked selection. However, we did not find evidence for positive selection in the mitogenome of Bufotes pewzowi and this suggests that gene adaptations that improve osmoregulation (if they occur) may be associated with nuclear genes. Hence, more information on the evolution of nuclear genomes in toads is needed to analyze the potential adaptation to salinity.