Snake mitochondrial genomes: phylogenetic relationships and implications of extended taxon sampling for interpretations of mitogenomic evolution

Mt genomes of Anilius scytale, Rhinophis philippinus and Charina trivirgata

The entire mt genome of C. trivirgata, except for the repeat regions within the control regions, was sequenced. All mt genes in A. scytale and R. philippinus were sequenced except for the tRNAs flanking the duplicate control region situated within the IQM tRNA cluster (Fig. 1). As previously described for other alethinophidian genomes [1–4, 9], the three new mt genomes possess two control regions. This was confirmed in all three species despite incomplete sequencing because it was possible to sequence both the 5' and 3' ends of both control regions in C. trivirgata and the 5' ends of both control regions for A. scytale and R. philippinus. In all three genomes prominent C-rich regions - characteristic of the control region 5' end in the majority of vertebrate mt genomes - were found. We failed to sequence the 3' ends and tRNAs adjacent to CRII for both A. scytale and R. philippinus, possibly because of extensive repeat regions typically found at the 3' end of the control region. While this paper was in preparation Castoe et al. [10] sequenced the mt genomes of A. scytale and Tropidophis haetianus. The mt genome of T. haetianus also has a duplicate control region and a translocated tRNA^LEU gene.

As in other alethinophidian mt genomes, the tRNA^LEU (UUR) gene in C. trivirgata was found translocated from its typical vertebrate position between the genes 16S and ND1 to downstream of CRII. The same can be assumed for A. scytale and R. philippinus because this gene was not found between 16S and ND1. As in other alethinophidian snakes, the origin of light strand replication is present in all three newly sequenced mt genomes, with the stem being 12 bp long.

Phylogenetic analyses

The total number of nucleotides in the alignment was 12316 (8874 excluding 3rd codon positions), with 3442 amino acid sites. The maximum likelihood (ML) and Bayesian phylogenetic analyses performed in this study yielded trees with unequivocal support for major snake taxa regardless of the method or model used. As we explain below, many of the differences in analytical results were relatively minor and did not conflict strongly in that the contentious nodes were weakly supported. The only strongly conflicting results are shown in Fig. 2 and 3, where boids are sister to either Python + Xenopeltis + Cylindrophis + Rhinophis (Fig. 2) or to all alethinophidians except Anilius and Tropidophis (Fig. 3).

The maximum likelihood (ML) GTR+I+Γ tree based on the 4-partition nucleotide model is shown in Fig. 2. This is identical to a tree recovered by running the 2nd codon positions under the TN+I+Γ model, and the Bayesian analysis in which nucleotide data was split into 26 partitions (see Methods). The Bayes factor for the 26- partition model (H₁) against the 4-partition model (H₀) was B₁₀ = 0.989 (2lnB₁₀ = -0.02), indicating that there was no significant difference between the two models.

Only node C (Fig. 2) is strongly supported or conflicted in different analyses. MrBayes analyses and amino acid ML analyses recovered node C with very high or maximal support (Fig. 2, Table 3), whereas this node received <70% bootstrap/LR-ELW support in ML nucleotide analyses. CAT-GTR and CAT-Poisson analyses recovered the alternative topology shown in Fig. 3 with ≥ 95 BPP. Interestingly, the mtREV analysis run with PhyloBayes also recovered the boid relationship shown in Fig. 3, albeit with only 0.65 BPP (see Additional file 2). Statistical tree comparisons with nucleotide data could not reject the tree shown in Fig. 3 as significantly suboptimal compared to that in Fig. 2 (Table 2), the difference in likelihood between the two trees only negligible (ΔlnL = 4.8). With amino acid data however, the CAT topology (ΔlnL = 23.47; Fig. 3) was rejected by all statistical tests at the 0.05 significance level (Table 2).

Relative rates of molecular evolution of mt genes

Relative nucleotide branch lengths of rRNA genes and protein-coding genes can be visualized in Figs. 4, 5 and 6. As expected, rRNA branches (Fig. 4) are substantially shorter overall than protein-coding branches (Fig. 5). CAT analyses also produced the same trend (data not shown). Branches leading to Serpentes, Alethinophidia, Typhlopidae and L. humilis show especially accelerated amino acid change relative to rRNA nucleotide branch lengths (Fig. 6). Our phylogenetic analyses produced two alternative resolutions of basal nodes among core henophidians (Figs. 2 and 3) and so the rates analyses were carried out twice. The different trees barely affected relative evolutionary rates, and the results shown here are based on the phylogeny in Fig. 2.

Patterns of evolution for different genes along different branches found in this study are broadly concordant with those presented by Jiang et al. [3]. Branch lengths of rRNA and protein-coding genes are generally positively correlated, positive deviances from this correlation being most prominent for CO genes, ATP genes, CYTB and ND6 (Fig. 7), with only slight deviances in other genes/clusters. These same trends were evident under both of the nucleotide partitioning schemes used (see Methods section 2.6), and in the plot of rRNA vs. amino acid branch lengths (Fig. 8), the latter suggests that markedly raised protein-coding gene evolutionary rates are not simply due to large numbers of synonymous changes. The most dramatic acceleration in relative rate of evolution occurred along the branches leading to Serpentes and to Alethinophidia, for which there was an acceleration in most protein-coding genes. Relative rate accelerations were found along the branch leading to A. scytale + T.haetianus in COI and ND6 (Fig. 7). There were no dramatic relative rate shifts in overall protein-coding gene evolution along the branch leading to core alethinophidians, although there was a slight relative acceleration in CYTB evolution here (Fig. 8). The other branch that showed the most dramatic elevation in protein-coding gene relative rates was that leading to Typhlopidae, which had elevated relative rates in COI, CO2, ATPs, ND2, ND4, ND5 and ND6. However, relative rates of evolution in most protein-coding genes along terminal typhlopid branches and the internal branch leading to Ramphotyphlops are moderately decelerated, also reflected in the amino acid plot (Figs. 7 and 8). In the rRNA tree (Fig. 4), typhlopids have the longest branches of all snake taxa except A. granulatus and L. humilis, the latter having elevated rates in all genes. The branch leading to Scolecophidia has relatively accelerated rates for COI, CO3, CYTB, ND3, ND4L and ND6, but dramatically decelerated for ND2 and ND4. The branch leading to Colubroidea occurs below the scatter distribution in the graphs for most genes, ND6 being a notable exception (Fig. 7). This deceleration is also seen in the amino acid plot (Fig. 8). The branch leading to Elapidae shows slight accelerations in rate in genes ND1, ND4 and ND6, which is also reflected in the amino acid plot (Fig. 8). Rates of evolution in ND1 show indications of relative acceleration in several terminal branches within Henophidia.

Relationships among the major lineages of snakes

The phylogenetic relationships among non-caenophidian snakes has fundamental importance for the interpretation of early snake evolution [7, 12, 13], but despite increased geno- and phenotypic character sampling over recent years, explicit quantitative phylogenetics has yet to deliver an entirely robust tree. Our analyses confirmed strong support for many nodes that have been robustly resolved in other studies, for example Colubroidea, Caenophidia, Colubridae, and nodes that were supported by other molecular studies but found to be equivocal or generally not recovered by morphological studies (Python+Xenopeltis, core alethinophidians, Cylindrophiidae + Uropeltidae). However, as with previous molecular studies, most deeper henophidian nodes could not be resolved robustly in this study. Low node support is often associated with short internal branches that can result from incongruence within data and/or too few character changes (e.g. [14, 15]). As is evident in Figs. 2 and 3, there are many short internal branches among henophidians, but that this is not a feature of mitochondrial data alone is indicated by relatively short branches and problematic nodes in analyses using nuclear data also [15].

Prior to the first broad-scale molecular snake phylogenies, morphological phylogenies had reached a consensus on several relationships among major snake lineages [7]. Morphologists were not entirely surprised by some of the new molecular findings because they mostly affected the less robustly inferred regions of morphological phylogenies, including the possible close relationship between pythons, Xenopeltis, and Loxocemus, the non-monophyly of macrostomatan alethinophidians and the non-monophyly of dwarf boas [16]. However, the sister-group relationship between Anilius and tropidophiids proposed by several molecular studies [7, 17–20] remains a notable sticking point between molecular and morphological studies [7], especially considering that recovery of this clade in larger nuclear gene studies [15] negates dismissal of this hypothesis as an aberration of mt data alone. It might be noted that this relationship was not very strongly supported in our analyses.

Mitogenomic evolution in snakes

Previous studies (e.g. refs [3, 4]) have mapped the duplication of the control region and translocation of tRNA^LEU to the internal branch leading to extant alethinophidians, but incomplete taxon sampling (especially the absence of Anilius) had left this proposition incompletely tested. Our identification of these same derived mitogenomic features in A. scytale provides empirical support for this. The same features are present in all alethinophidian genomes sequenced to date.

Our analyses of relative rate variation in mt gene evolution across the snake tree produced several results in concordance with the study of Jiang et al. [3], including overall dramatically accelerated (relative to rRNAs) protein-coding gene evolutionary rates along the internal branches leading to Serpentes and to Alethinophidia though with notable gene- and branch-specific variation (e.g., dramatic relative acceleration in evolution of COI, COII, ATP8, ATP6, ND6 and CYTB along branch leading to Alethinophidia). Notable differences between our results and those of Jiang et al. [3] are: no snake branches disproportionately longer for rRNA (Fig. 4) than combined protein-coding (Fig. 5) genes that would suggest dramatic relatively accelerated rRNA evolution; no notable acceleration in relative rate of CYTB evolution (and less dramatic increase in ND2) along branch leading to Serpentes; internal branch leading to Colubroidea with relatively decelerated COI and CYTB evolution (Fig. 7). The latter is potentially dependent on the position of Achalinus meiguensis, which was not resolved robustly here, but analyses with this taxon excluded produced the same patterns (data not shown).

This study has generated new findings, notably: an overall dramatically accelerated protein-coding gene evolution on the branch leading to Typhlopidae; terminal branches within Typhlopidae and the internal branch leading to Ramphotyphlops generally falling below the main distribution in gene-gene plots (Fig. 7); relative acceleration in COI and ND6 along the branch leading to A. scytale + T. haetianus; the terminal branch leading to A. scytale slightly accelerated for CYTB, ND4, ND5 and ND6; further relative acceleration in CYTB and ND6 evolution along branch leading to core alethinophidians; relative acceleration in ND1, ND4 and ND6 along the branch leading to Elapidae. Although we did not carry out the in-depth analytical tests of positive selection or protein structural modelling performed by Castoe et al. [6], the rRNA vs. protein-coding amino acid branch length plot (Fig. 8) demonstrates that some of the acceleration in molecular evolution along the branches leading to Serpentes, Alethinophidia and Typhlopidae is likely to be adaptive in that it leads to changes in protein sequences, and is not simply an increase in 'silent' substitutions.

Our greater taxon sampling enabled us to identify additional trends in mt gene evolution for some snake clades. Within Scolecophidia, we found dramatic relative acceleration in protein-coding gene evolution along the branch leading to Typhlopidae but subsequent relative deceleration (relative acceleration in rRNA evolution) along all terminal typhlopid branches. The terminal branch joining Leptotyphlops humilis to the rest of the tree appears above the main distribution in most gene-gene plots (for nucleotide and amino acid data), but comparison with other lineages is complicated because only one leptotyphlopid mt genome has been sequenced. Additional sampling of leptotyphlopids is required to clarify mitogenomic evolution along this branch. The inclusion of L. humilis is not responsible for the rapid protein evolution identified along the branch leading to Typhlopidae because this branch is similarly long in amino acid trees that exclude L. humilis (data not shown). The generally long branches within Scolecophidia and seemingly strongly fluctuating rate dynamics suggest that mitogenomic evolution within this group is worthy of greater attention. Although scolecophidians represent one half of the basal divergence among living snakes they are often overlooked, perhaps because of their superficial morphological and ecological homogeneity, which is, however, deceptive [21].

Castoe et al. [6] observed that dramatic adaptive shifts in the evolution of metabolic proteins occurred along the branch leading to Alethinophidia and so were potentially temporally associated with the origin of a duplicated control region, and they suggested that these shifts along this internal branch were causally associated with major anatomical, ecological, and behavioural changes such as an ecological niche shift to a non-burrowing lifestyle, increased body size, increased skull kinesis and gape and prey size, development of specialized venom proteins, and the ability of individuals to dramatically remodel their organs and physiology. The scenario of a single major switch from a burrowing to non-burrowing (and macrostomatan) condition early in alethinophidian history has been eroded by some of the more recent phylogenetic results for henophidian snakes, so that identifying broad suites of phenotypic adaptive change along the branch leading to Alethinophidia is not as trivial as previously assumed. However, the rapid remodelling of metabolic proteins before the burrowing Anilius diverged from other alethinophidians and, importantly, evidence of elevated relative rates of metabolic protein evolution in the branch leading to the exclusively burrowing Scolecophidia (which lack a duplicated control region) rejects the exact correlations that underpin some of the precise aspects of the hypotheses previously proposed to explain the observed mitogenomic evolutionary patterns.

Taxon sampling, DNA extraction, PCR and sequencing

Data assembly and alignment

Sequences were edited and assembled using EditView 1.0.1 (Perkin-Elmer). Gene organization within the mtDNA sequences was inferred by comparisons with other snake genomes using the program SeAl (v2.0a11). The twelve H-strand encoded mt gene sequences from all taxa were aligned using MAFFT [23] and concatenated. 12S and 16S ribosomal RNA genes were aligned using the program T-COFFEE [24] incorporating information from secondary structure [25]. All mt gene alignments were subsequently inspected manually. Alignment-ambiguous characters (where the program was equivocal in the placing of certain gaps) were excluded from the dataset.

Data partitioning and modelling

The data were analysed separately as nucleotide and amino acid sequences. The nucleotide data were subdivided into three partitions (rRNAs; protein-coding 1st codon positions; 2nd codon positions). The phylogenetic signal in 3rd codon positions was investigated by running preliminary parsimony and neighbour joining analyses using PAUP* [26], and found to contain substantially more signal when analysed under RY coding than as raw nucleotide data. In RY coding, pyrimidine bases C and T are analysed as one character, as are purine bases A and G. RY-coded third positions were thus included as a fourth partition in ML and MB analyses. Third codon positions were not included in PhyloBayes analyses because RY coding is not possible on this platform (N. Lartillot, pers. comm.).

Whether further partitioning of the data would yield an improvement over the 4-partition dataset was investigated by calculating the Bayes Factor (B₁₀), which is robust against Type I error ([27], J. Brown, pers. comm.). A dataset with 26 partitions (12S rRNA, 16S rRNA, ND1, ND2, CO1, CO2+ATP8+ATP6, CO3+ND3+ND4L, ND4, ND5, CYTB, with all protein-coding genes being split into three partitions comprising 1st, 2nd and 3rd codon positions) represented H₁ and the 4-partition dataset represented H₀. Third codon positions were coded as RY. Both datasets were analysed using the program MrBayes (see below). The harmonic means were calculated using the program Tracer (v 1.4) and incorporated into the formula: B₁₀ = harmonic mean L₁/harmonic mean L₀. The Bayes Factor was interpreted as in Kass & Raftery [28].

The model of choice for each partition was determined using the program Modeltest [29]. Both the Likelihood Ratio Test (LRT) and Akaike Information Criterion (AIC) found the GTR+I+Γ model [30, 31] as best fitting for all partitions except 2nd codon positions, for which LRT found TN+I+Γ [32] as best fitting.

Amino acid sequences were analyzed using the mtREV+I+Γ model [33]. The categories model with GTR (CAT-GTR) and CAT-Poisson models [34, 35] were also used to analyse nucleotide and amino acid data, respectively.

Phylogenetic analyses

Maximum likelihood (ML) analyses were performed in Treefinder [36]. For nucleotide sequences each analysis was started from ten different initial trees to increase the amount of tree space explored. Amino acid analyses were initiated from five different start trees. Support for nodes was determined using expected likelihood weights of local rearrangements (LR-ELW) with 1000 replicates for each node. Both nucleotide and amino acid data were also bootstrapped with 200 and 100 replicates, respectively. Bayesian analyses run under GTR+I+Γ and mtREV+I+Γ models were performed using MrBayes [37, 38] freely available through Bioportal cluster [39]. Nucleotide and amino acid data were run for 10,000,000 and 5,000,000 generations, respectively. Convergence was checked using Tracer. Analyses using CAT models were run with the PhyloBayes 2.3 program package [34] also available on Bioportal. Here amino acid data were analyzed using both CAT-Poisson and mtREV+I+Γ models. Convergence was checked with the bpcomp program, whereby convergence was reached if the maxdiff value = <0.1.

Because our primary concern was the relationships of snake lineages, the main phylogenetic analyses were not conducted with (distant) non-snake outgroups. Instead, we rooted Scolecophidia with Alethinophidia and vice versa, based on this being a well-supported relationship in previous morphological and molecular estimates of snake phylogeny [3, 7, 13, 15, 40].

Tree comparisons

Alternative topologies were compared statistically using the program Treefinder. Both nucleotide and amino acid data were analyzed. Six statistical tests were carried out simultaneously: Bootstrap Probability [41], Expected-Likelihood Weights (ELW) [42], Kishino-Hasegawa (KH) [43], Shimodaira-Hasegawa (SH), Weighted SH (WSH) [44] and Approximately Unbiased (AU) [45], all implemented in Treefinder [36].

Relative rates of evolution of mt genes

In accordance with Jiang et al. [3], we analysed temporal patterns of molecular evolution, measured by nucleotide branch length, among the different branches of the snake tree for individual protein-coding genes relative to branch lengths for 12S + 16S rRNAs. We also carried out similar analyses based on amino acid sequences in order to gain some insight as to whether notably rapid nucleotide change in protein-coding genes was associated with (potentially adaptive) change at the protein level. Because we wanted to include estimates of the length of the branch connecting snakes to other squamates, we included the following selection of lizard outgroups (Genbank accessions in parentheses): Gekko gecko (NC 007627), Coleonyx variegatus (NC 008774), Cordylus warreni (NC 005962), Eumeces egregius (NC 000888), Iguana iguana (NC 002793), Sceloporus occidentalis (NC 005960), Heloderma suspectum (NC 008776), Shinisaurus crocodilurus (NC 005959). The phylogeny used to calculate branch lengths was fixed, with the relationships among snakes based on our phylogenetic results, and those among lizards following refs [46] and [47].

The branch lengths of rRNA and protein-coding genes were compared by firstly estimating branch lengths separately for rRNAs and all protein-coding genes using both nucleotide and amino acid data. Because of the erosion of signal at 3rd codon positions (see section 2.3) 1st, 2nd and 3rd codon positions were weighted 2, 1 and 5, respectively. ND6 was analysed as a separate partition because of its aberrant base composition. The amino acid analysis was run without partitioning. Secondly, rRNA branch lengths were compared with those for individual protein-coding genes in bivariate plots. Following Jiang et al. [3] the short genes CO2, ATP8 and ATP6 were concatenated into one partition to reduce stochastic error, as were genes CO3, ND3 and ND4L. Two different partitioning schemes were used for each protein-coding gene/cluster: a) 1-partition with 1st, 2nd and 3rd codon positions assigned different weights as described above and b) equal-weighted 3-partitions based on codon position, with RY coding of 3rd codon positions. This second partitioning scheme was used as an alternative to amino acid data that, because of the shorter gene lengths, would increase the risk of stochastic error. The use of RY coding reduces variation at 3rd codon positions considerably, which means that the branch lengths would more accurately reflect evolution at 1st and 2nd codon positions. Nucleotide and amino acid branch lengths were estimated in Treefinder. Branch lengths for each protein-coding gene/cluster were plotted against the respective branch lengths for rRNAs in Microsoft Excel.