Bothrops mitogenomes characterization and evolution
We present the first complete and structurally annotated mitochondrial genome of B. insularis reconstructed by a reference-guided assembly. B. insularis mitochondrial DNA molecule is similar to the mitogenome of B. jararaca (ALMEIDA et al., 2016), but differ in gene order and structure (presence or absence of components) when compared to B. diporus and B. pubescens (KLEIZ et al., unpublished). These divergences can be associated to duplication process of some regions: (i) a non-coding region (NC) followed by a duplicate tRNA-Phe (F*) in B. insularis and B. jararaca mitogenomes; and (ii) a duplicate tRNA-Pro (P*) in B. pubescens mitogenome.
Vertebrate mitogenomes tend to bear small non-coding regions and sometimes duplicates of tRNAs (BOORE, 1999; LANG et al., 1999; TANMAN, 1999), and their existence may be explained in part as a vestige of rearrangements that modified genes/tRNAs/control regions positions and orders (MACEY et al., 1997; PEREIRA, 2000; BERNT et al., 2013b). Duplications of tRNAs have been proposed for many taxonomic groups (see JÜHLING et al., 2012) and it is worth mentioning that these duplicated segments can become nonfunctional due to the accumulation of deleterious mutations in these molecules over time (MABUCHI et al., 2004).
Two plausible (and not mutually exclusive) hypotheses for the origin of duplicate regions in Bothrops mitogenomes can be suggested: (i) the occurrence of a duplication process involving the end of the control region I and the tRNA-Phe (F) with translocation of the duplicate regions prior the rRNA-12S subunit, which can explain the similarity of the non-coding region (NC) to the end of the control region I (< 90%) and the tRNA-Phe duplicate (F*) to the tRNA-Phe (F) (< 80%) on both B. insularis and B. jararaca molecules; and (ii) the occurrence of a duplication process involving the tRNA-Pro (P) with translocation of the duplicate region prior the CRII on B. pubescens molecule.
The methods implemented here for detection of gene-wise positive selection, or a shift in its strength (RELAX and aBSREL), did not show evidence for changes associated with either origin or diversification of the B. jararaca + B. insularis clade, neither for any of these individual taxa, compared to relatively more ancient branches in the phylogeny. This is in accordance with previous studies showing a lack of evidence for gene-wise positive selection in vertebrate mtDNA protein-coding regions (BALLARD & KREITMAN, 1995; BALLARD & RAND, 2005), albeit this molecule is pointed to evolve under purifying selection due to its role in the ATP synthesis and oxidative phosphorylation processes (SUN et al., 2018).
Snake mitogenome evolution perspectives and insights
The comparative analyses using 129 mitogenomes of snake species, distributed in 18 families, allowed the inclusion of 64 species and four families in the snake mitogenome evolution scenario. The improvement of new families was due the: (i) use of the taxonomic propositition suggested in the phylogenomic revision for the Colubroidea superfamily which were divided in 18 families (ZAHER et al., 2019); and (ii) the inclusion of a new specimen of the family Gerrophilidae.
The phylogenetic analyses partially recovered most evolutionary relationships found for Squamata in phylogenomic studies (for example, ZHENG & WEINS, 2016 – used 52 genes and 4162 species; and ZAHER et al., 2019–6 mitochondrial and 9 nuclear genes, and 1278 species), however some divergences should be highlighted. The first one is related to the Scolecophidia paraphyly that was recovered previously (MIRALLES et al., 2018; QIAN et al., 2018) but it disagrees with the monophyly found by Zheng and Weins (2019). Once again it is worth mentioning that the paraphyletic clades recovered herein (Clade 1: Gerrophilidae and Typhlopidae; Clade 2: Leptotyphlopidae) differ from the ones recovered previously (Clade 1: Leptotyphlopidae, Gerrophilidae, Xenotyphlopidae and Typhlopidae; Clade 2: Anomalepididae; MIRALLES et al., 2018). Moreover, we recovered Natricidae as the sister clade to Dipsadidae, whereas Natricidae was recovered as the sister clade to (Sybnophiidae + Colubridae) by Zaher and collaborators (2019). These divergences might be related to the mtDNA genes being representative of a single non-recombining haplotypic molecules, thus not being representative of the set of independent gene topologies present in studies embracing phylogenomic of different species (e.g., due to forces such as incomplete lineage sorting). Furthermore, such cases of inconsistent phylogenetic signal may also be due to a smaller set of genes compared to phylogenomic analyses, different parameters implemented in the analyses.
Aside from the phylogenetic analyses, comparative analyses focusing on gene order and structure allowed the recognition of 18 mitotypes within snakes. Although some of them had been previously described, others were revised after analyses involving unannotated regions (as described in the results section). The presence of unannotated regions might be related to automatic annotation of these molecules, which can lead to errors and bias on classification and identification of its components as observed in mammals by Prada and Boore (2019).
Rearrangements of gene order have been observed in tRNA clusters in the mitogenome structure and some authors pointed them out as rearrangement hotspots due to the higher probability of accumulating mutations (KUMAZAWA & NISHIDA, 1995; MACEY et al., 1997). Divergence involving mitogenome components found in this study are associated with rearrangements of non-coding regions or tRNAs, specifically in four clusters: (i) the WAN-Ol-CY; (ii) the CR (I or II) + adjacent tRNA; and (iii) the S2D.
Herein, we highlighted that these rearrangements may have led to: (i) the presence of a tRNA N duplicate in Hydrophis cyanocinctus; (ii) the inversion in the order of the OL region and the tRNA N in H. curtus, H. melanocephalus, H. platurus, and Emydocephalus ijimae; (iii) the presence of a tRNA-F duplicate (F*) and a non-coding region in B. insularis and B. jararaca (previously described for B. jararaca in ALMEIDA et al., 2016); (iv) the presence of a non-coding region in Daboia ruselli; (v) the presence of a tRNA-I duplicate vestige (I**) in Ophiophagus Hannah and the presence of a tRNA-I duplicate (I*) in Hypisiglena jani texana 2; (vi) the presence of a tRNA-Ser2 duplicate (S2*) in H. curtus and H. cyanocinctus; and (vii) the presence of a tRNA-Asp duplicate (D*) in Micrurus fulvius.
Even though the original paper that described the mitogenomes of O. hannah reported a presence of the duplicated tRNA I* sequence (CHEN & LAI; 2010), our analysis recovered an additional duplicated tRNA-I** vestige prior the duplicate tRNA I*. Furthermore, the description of the Hypsiglena jani texana 2 were not available when the last snake mitogenome revision were done (QIAN et al., 2018; MYERS & MULCAHY, 2020). Then, they are also included as reannotated or novelty in this revision.
Several groups of metazoans show variations in the WAN-Ol-CY cluster order (for example, see MUELLER & BOORE, 2005; MAURO et al., 2006; FUJITA et al., 2007; BERNT et al., 2013a) and in the CR (I or II) and adjacent tRNA cluster order (for example see BERNT et al., 2013a). In snakes, variations related to these clusters involved the lack of the OL, the translocation of the tRNA-Q, and the presence of duplicates of the tRNA-P (P*), tRNA-I (I*), and tRNA-F(F*) and their translocations (CHEN & ZHAO, 2009; ALMEIDA et al., 2016; QIAN et al., 2018; MYERS & MULCAHY, 2020). However, variations in the S2D cluster order had not been described in snake mitogenomes and are considered as one of this revision novelties.
Although the rearrangements highlighted in the last paragraphs occur in different regions of the mitogenomes, their origin might be explained by the occurrence of errors in the mtDNA replication process associated with the highest mutation rates in snake mitogenomes when compared to other animal taxa (KUMAZAWA & NISHIDA, 1995; MACEY et al., 1997; JIANG et al., 2007; CASTOE et al., 2009; BERNT et al., 2013a). These errors may lead to the occurrence of rearrangements, such as the tandem duplication and random loss (TDRL) process followed by the translocation of these duplicates or the original segment promoting the gene order variations, as observed.
Moreover, the mitogenome complete structure observed in some species (for example Hydrophis melanocephalus) has the presence of largest nonannotated segments in the WAN-OL-CY cluster that weren’t recognized as tRNAs duplicates due to the low similarity (< 50%). Therefore, they might be related to vestiges of TDRL processes involving the OL region.
Mitogenomes 1 and 2 diverge firstly and were associated with Scolecophidian species. As proposed by Yan and collaborators (2009), our results indicate that mitogenome 1 diverge firstly followed by the mitogenome 2, whereas differ from Qian and collaborators (2018), which recovered the mitogenome 2 diverging firstly followed by the mitogenome 1. Both mitogenomes found in Scolecophidia could have diverged independently from the same common ancestors that possessed a mitogenome whose composition is not clear. However, a similar characteristic of these mitogenomes is the lack of the OL region which might (or not) be associated with these snakes’ fossorial habits (MIRALLES et al., 2018).
Mitogenome 3 and its derivatives (3A-J) are distributed in Alethinophidia, with the mitotype 3 as the most frequent and observed in all families but Tropidophiidae, Viperidae and Homalopsidae. Due the dimness of the information about the common ancestor mitogenome structure, we proposed (as QIAN et al., 2018) that the mitotypes (1, 2 and 3) found for Scolecophidia and Aletinophia evolved independently from a common ancestor.
The dominance of the mitotype 3 in Alethinophidia might be associated with the diversification rates and the dispersal observed in snake evolution (KLEIN et al., 2021). Derivatives from the mitogenome 3 observed in Alethinophidia occur in specific species, except the mitogenome 3B, 3C, and 3D.
The mitogenome 3B and its derivatives are exclusive of the family Viperidae, these mitotype might be originated on the divergence episode of these family in the medium to late paleogene (KLEIN et al., 2021) and might be coevolved with specific characteristics of this family (such as its venomic traits).
The mitogenome 3C have multiple origins – firstly observed in Homalopsidae - and might be related to the divergence episode of this family in the late Paleogene (KLEIN et al., 2021). The other 3C mitotypes probably have emerged in episodes of the diversification of Colubridae and Dipsadidae in the Neogene (KLEIN et al., 2021).
At last, the mitogenome 3D and its derivatives are observed in exclusively marine snakes. Marine snake genomes showed chromosomal rearrangements, genomic and venomic adaptations related to the marine lifestyle adaptation (PENG et al., 2020; LI et al., 2021). Then, we could suggest that the mitogenome 3D recovered herein for marine species (E. ijimae and Hydophys sp.) could also be an adaptation to the marine lifestyle, even though it is not possible to demonstrate a direct association with these variables based on the data provided herein.