Comparative mitochondrial genomics of snakes: extraordinary substitution rate dynamics and functionality of the duplicate control region
- Zhi J Jiang†1, 2,
- Todd A Castoe†3, 4,
- Christopher C Austin5,
- Frank T Burbrink1, 6,
- Matthew D Herron3, 7,
- Jimmy A McGuire5, 8,
- Christopher L Parkinson3 and
- David D Pollock9Email author
© Jiang et al; licensee BioMed Central Ltd. 2007
Received: 06 September 2006
Accepted: 26 July 2007
Published: 26 July 2007
The mitochondrial genomes of snakes are characterized by an overall evolutionary rate that appears to be one of the most accelerated among vertebrates. They also possess other unusual features, including short tRNAs and other genes, and a duplicated control region that has been stably maintained since it originated more than 70 million years ago. Here, we provide a detailed analysis of evolutionary dynamics in snake mitochondrial genomes to better understand the basis of these extreme characteristics, and to explore the relationship between mitochondrial genome molecular evolution, genome architecture, and molecular function. We sequenced complete mitochondrial genomes from Slowinski's corn snake (Pantherophis slowinskii) and two cottonmouths (Agkistrodon piscivorus) to complement previously existing mitochondrial genomes, and to provide an improved comparative view of how genome architecture affects molecular evolution at contrasting levels of divergence.
We present a Bayesian genetic approach that suggests that the duplicated control region can function as an additional origin of heavy strand replication. The two control regions also appear to have different intra-specific versus inter-specific evolutionary dynamics that may be associated with complex modes of concerted evolution. We find that different genomic regions have experienced substantial accelerated evolution along early branches in snakes, with different genes having experienced dramatic accelerations along specific branches. Some of these accelerations appear to coincide with, or subsequent to, the shortening of various mitochondrial genes and the duplication of the control region and flanking tRNAs.
Fluctuations in the strength and pattern of selection during snake evolution have had widely varying gene-specific effects on substitution rates, and these rate accelerations may have been functionally related to unusual changes in genomic architecture. The among-lineage and among-gene variation in rate dynamics observed in snakes is the most extreme thus far observed in animal genomes, and provides an important study system for further evaluating the biochemical and physiological basis of evolutionary pressures in vertebrate mitochondria.
The vertebrate mitochondrial (mt) genome has been an important model system for studying molecular evolution, organismal phylogeny, and genome structure. Despite extensive molecular studies, little is known regarding the ways in which genome architecture might affect the various aspects of genome function and evolution (including replication, transcription, and RNA/protein function, as well as rates and patterns of nucleotide evolution). Nevertheless, patterns linking mt genome structure, function, and nucleotide evolution have begun to emerge [1–3].
Among the most direct demonstrated links among genome architecture, function and nucleotide evolution is that relating the asymmetrical genome replication process with gradients of transition substitutions in vertebrate mitochondrial genomes [1–3]. Gradients of transition mutations, arising from deamination mutations, are observed due to the differential time regions of the mt genome spend in an asymmetric mutagenic state during genome replication (T AMS alternatively referred to as the time spent in a single-stranded state, D SSH , [4–6], but there is some controversy about this: see Additional file 1). Thus, gradients of transition biases are dependent upon the relative position of the functional origins of heavy and light strand replication. In vertebrate mt genomes, the origin of heavy strand replication (OH) is thought to be within the control region (CR), and the origin of light strand replication (OL) in the tRNA cluster referred to as the WANCY region (named for the five amino acids coded for by these five tRNAs). Among transition classes in vertebrate mt genomes, T→C light strand substitutions at degenerate 3rd codon positions increase linearly with increasing T AMS and C/T nucleotide frequencies at degenerate 3rd positions are good predictors of T AMS .
The mt genomes of snakes contain a number of characteristics that are unusual among vertebrates, and represent an ideal model for exploring potential links among genome structure, function, and evolution. Snake mitochondrial genomes appear to have the highest evolutionary rates among vertebrates and contain truncated tRNAs and other shortened genes [7, 8]. All snake species sampled to date, except the scolecophidian snakes Leptotyphlops dulcis, Ramphotyphlops australis, and Typhlops murius, have a duplicated control region (CR2) between NADH dehydrogenase subunit 1 (ND1) and subunit 2 (ND2), in addition to a control region (CR1) adjacent to the 5'-end of the 12s rRNA as it is in other vertebrates [7–11]. These two control regions appear to undergo concerted evolution that acts to homogenize the nucleotide sequence of each duplicate copy within a given genome [7–9]. The functionality of these two control regions in transcription and initiation of heavy strand replication is not clear, but given that the nucleotide sequence of each is nearly identical, any functional features that are not dependent on surrounding sequences should be similar. In contrast, recent evidence suggests that initiation of heavy strand replication may be distributed across a broad zone, including cytochrome b (CytB) and NADH dehydrogenase subunit 6 (ND6) , indicating that CR2 may not function as effectively in this role.
A number of interesting questions arise that might be addressed through comparative analysis, including: (1) does one or the other, or do both control regions function as origins of heavy strand DNA synthesis? (2) does the altered genome structure affect patterns of snake mt genome molecular evolution? (3) when during snake evolution did various features arise, and were any changes synchronous? (4) do patterns of mt molecular evolution vary at different depths of phylogeny? and (5) is there any evidence or plausible rationale for selection as a causative agent in generating differences in genomic structure and molecular evolutionary patterns?
To investigate outstanding questions regarding snake mitochondrial genome evolution, structure, and function, we analyzed a dataset consisting of three new complete snake mitochondrial genomes together with all eight previously published snake mitochondrial genomes that were available at the time of this study, and 42 other vertebrate mitochondrial genomes for comparative purposes. The new snake genomes were obtained from one Pantherophis slowinskii (Colubroidea: Colubridae; a corn snake from Louisiana; previously Elaphe guttata), and from two Agkistrodon piscivorus (Colubroidea: Viperidae; the cottonmouth or water moccasin; one specimen from Florida and the other from Louisiana).
Brief summary of the new complete snake mitochondrial genomes
Within a mt genome, the two copies of the CR in each newly sequenced species are nearly identical (e.g., Api1 CR1 and CR2), as is typical for alethinophidian snakes [7, 8]. In P. slowinskii there is a single point mutation and four extra nucleotides at one end of CR1, in Api1 there is one indel plus 14 extra nucleotides on one end of CR1, and in Api2 there are seven indels and two base changes between the two control regions. Between Api1 and Api2, CR1 differs by five indels and 19 point mutations, whereas CR2 differs by three indels (two at the 5' end) and 18 point mutations.
Comparison of Agkistrodon piscivorus genomes
Within A. piscivorus, the control regions (e.g., CR1 in Api1 vs. CR1 in Api2) are as similar to each other as are the rRNA genes, and more similar than the protein coding genes (Figure 2A). This is in strong contrast to the normal pattern of divergence between vertebrate species, for which control region similarity is far less than that of protein-coding or rRNA genes, e.g., [13, 14]. Between A. piscivorus and the other viperid, O. okinavensis, the control regions have 30% more differences (with indels included) than the rRNAs, and are on par with divergence in the protein-coding genes (Figure 2B). If indels are included, the control regions of these two species are nearly as different as the average 3rd codon position (Figure 2B). The high degree of similarity (low divergence) observed between the CRs of the two A. piscivorus individuals is surprising, and contrasts sharply with the high relative divergence of CRs between O. okinavensis and A. piscivorus (Figure 2).
The snake and overall amniote phylogeny are strongly supported by our analysis of this dataset, and we henceforth treat this phylogeny as accurate. We wish to emphasize, however, that the consistency of the phylogenetic results do not guarantee that they are, in fact, accurate. Although for simplicity we present a single nucleotide substitution model for the entire dataset we have analyzed an expanded version of this dataset (with additional unpublished snake and lizard mt genomes; with and without inclusion of the rRNA genes) using complex partitioned models for each gene and codon position. The results of this expanded (highly partitioned-model) phylogeny estimate (not shown) were essentially identical to those presented here in terms of the placement of snakes within squamates, and relationships among squamates. We provide evidence below for extremely complex non-stationary patterns of nucleotide substitution across branches and mt genome regions, and have previously identified asymmetric substitution gradients in mt genomes  that may vary among species (e.g., primates ). These latter patterns cannot be modeled using available phylogenetic programs (e.g., MrBayes ). We expect our phylogenetic estimates here to represent a good estimate of the relationships among mt genomes sampled, and if minor inaccuracies in the topology have occurred in our estimates, these changes should not substantially impact the qualitative conclusions of further analyses (e.g., sliding window analysis, SWA) because a majority of these later estimates are averaged over many branches of the tree, and the dynamics we concentrate on are quite dramatic and are likely to be obvious and qualitatively similar even with slight changes in the topology estimate.
Nucleotide frequencies and control region functionality
In A. piscivorus and P. slowinskii mt genomes, as in other vertebrates , nucleotides A and C are favored on the light strand, particularly at 3rd codon positions. This bias is probably related to elevated rates of deamination mutations on the heavy strand incurred during replication (see Background), and there is considerable variation in nucleotide content among individual mt genomes (see Additional file 2). Variation in snakes, even at 3rd codon positions, is not exceptional compared to other groups, and there is no clear snake-specific nucleotide bias evident (see Additional file 2) or strong branch-specific, or gene-specific nucleotide bias shifts across squamate mt genomes that would explain our findings of dramatic branch-specific and gene-specific rate dynamics.
Due to the simple linear relationship in most vertebrate mt genomes between C/T ratios and T AMS predicted based on the location of the (functional) control region, it is of interest to determine whether there has been any clear genetic effect of the duplicated control region in alethinophidians. Exclusive use of one control region or the other would be most strongly observable in ND1 (the only protein-coding gene located between the two control regions in alethinophidian snake mt genomes) because it is the only protein-coding gene that would spend a substantially different amount of time in the asymmetric mutagenic state (T AMS see Additional file 2) depending on which control region is functional. Since the nucleotide sequence of duplicate control regions is nearly identical within each genome, however, it is also reasonable to consider the possibility that both control regions are functional.
Results of mitochondrial genome replication model analyses
Gene length and stability of truncated tRNAs in snakes
In snakes, all mt protein-coding genes (except COX1), ribosomal RNAs, tRNAs, and individual CRs are shorter than their counterparts in most lizards and most other vertebrates (see Additional file 3). An exception to this is Sphenodon punctatus, for which the control region, ATP8 (ATP synthase subunit 8) and the 12s rRNA are all shorter than in snakes. With the increased sampling in this study, it appears that while the tRNAs and proteins became shorter prior to the divergence of all snakes, the tRNAs became shorter still within the Colubroidea (Figures 4 and Additional file 3). Additionally, the rRNAs did not become shorter in L. dulcis or the Acro-Heno clade, but are dramatically shorter in the Colubroidea (Figures 4 and Additional file 3).
The shorter length of tRNAs in snakes results mainly from a truncated T-arm in the secondary structure (see also [8, 9]). In some tRNAs, the D-arm is also shorter, but to a lesser extent than the T-arms. Although short tRNAs are typically less stable than long ones, there is only a minor effect of sequence length on secondary structure stability (ΔG) in snake tRNAs. The cloverleaf structures of most snake tRNAs are slightly less stable than their lizard counterparts (see Additional file 2), but two tRNAs (tRNAIle, tRNAMet) are actually more structurally stable in snakes than in other squamates with longer tRNAs.
Spatio-temporal substitution rate dynamics across mt genes and genomic regions
Cross-referencing results from Figures 5, 6, 7, we can summarize the apparent nucleotide evolutionary rate dynamics in snake mt genomes as follows (see also Figure 4). The branch leading to all extant snakes appears to have experienced accelerated evolution in the region starting near the end of COX1 through COX2, ATP8, and somewhat into ATP6, and also in the region including the end of ND5, ND6, and CytB (and a rise in ND1). The COX1, COX2, ATP8, and ND6 accelerations increased and were stronger in the branch leading to Alethinophidia, while the ND5 acceleration decreased, and a notable acceleration of CytB also occurred. In the branch leading to the Colubroidea, only the ND6 acceleration continued, but new rate peaks arose in ND5, 12s rRNA, and the first part of the 16s rRNA, followed by a strong dropoff in all gene-specific acceleration in terminal colubroid lineages, except in the end of 16s rRNA in O. okinavensis. In the branch leading to the Acro-Heno clade, the accelerated rates of evolution (in COX1, COX2, ATP8, and ND5 genes) observed along the branch leading to the alethinophidians diminished (except for ND6 as in the Colubroidea), but new rate peaks arose in ATP6, COX3, ND3, ND4L, and the latter half of the 16s rRNA. These punctuated gene-specific accelerations were followed by the complete elimination of all gene-specific signals of atypical relative rate in terminal lineages within this Acro-Heno clade. We find no evidence for a constant accelerated rate of snake mt genome evolution. Instead, our analyses of rates and patterns of substitution underscore both the spatial (gene-specific) and temporal (branch-specific) nature of molecular evolutionary relative rate dynamics in snake mt genomes.
The three new complete snake mt genomes presented here, together with previously existing vertebrate mt genomes, provide a preliminary perspective on a complex history of potentially adaptive mt genomic change in snakes. Unusual changes in gene size and nucleotide substitution rates are associated with changes in mt genomic architecture (Figure 4). Nevertheless, the changes in substitution dynamics cannot be directly explained by the changes in mt genome architecture. Snake mt genome evolution is most consistent with some type of broad selective pressure on the efficiency and function of oxidative metabolism in snakes early in their evolutionary history.
In mt genomes (particularly in vertebrates), the processes of replication and transcription are not entirely functionally independent, and genome structural organization plays a prominent role in both processes. The CR acts as the origin of heavy strand replication, in addition to its role as the promoter for both heavy and light strand transcription . Genome replication also depends on the processing of light strand transcripts to produce short primers required for heavy strand initiation of genome replication (originating from the CR ). The regular distribution of the tRNA genes throughout the mt genome is functionally significant, and these play an important role in RNA processing of polycistrons to yield mature RNAs, transcription initiation and termination, as well as initiation of light strand replication . Collectively, many functional ramifications are linked tightly to mt genome architecture in vertebrates.
Mitochondrial genome size reduction due to gene shortening in alethinophidians is more than offset by the retention of their duplicate control regions. If size reduction is caused by selective pressure, the long term retention of dual CRs suggests that having both copies provide some selective advantage. Although the duplicate control region appears to function in heavy strand replication in at least some snakes, there is considerable variation in CR usage across snake lineages (Table 1). Thus, if the duplication has been maintained by selection, control of replication may not be the singular or primary selective driving force.
The possession of two functional control regions in most snake mt genomes might be advantageous by increasing the rate at which genome replication proceeds, and/or increasing the overall number of genome copies per mitochondrion. Since the dual CRs essentially flank the rRNA genes, they (along with adjacent tRNAs) could also plausibly function to independently control rates of protein-coding and rRNA gene transcription. Across snake species, variation in the tRNAs flanking the CRs includes the translocation of tRNALeu (3' of CR2) and the duplication/translocation/truncation of tRNAPro. In vertebrates, tRNALeu has been shown to decouple rates of rRNA and mRNA transcription by acting as a terminator of ~95% of heavy strand transcripts (leading to ~20-fold higher rRNA vs. mRNA levels; ). Considering the ectothermy of snakes, transcriptional decoupling via independent control regions could provide a more direct means of countering thermodynamic depression of enzymatic rates at low temperatures.
Independent CR duplications have also been identified in eels , frogs , birds [22, 23], and lizards [24, 25]. Our results (and additional unpublished data) suggest that the dramatic shifts in rates and patterns of molecular evolution in snakes represent a unique phenomenon that we do not expect to be necessarily associated with CR duplication, but rather more likely associated with selection for mitochondrial function. Nevertheless, these independent duplications may be useful to test the consequences of duplication on mutational processes.
Concerted evolution in and around the duplicate control regions
The two control regions clearly undergo concerted evolution to maintain reciprocal homogeneity between control regions within a genome [7–9], presumably through gene conversion. Interestingly, an apparently nonfunctional partial (or pseudo) proline tRNA (Ψ-tRNAPro) in colubrid mt genomes also appears to be maintained by concerted evolution (Figure 1). The gene conversion process that homogenizes the control region may also occasionally pick up extra DNA, making tRNAPro, or part of it, prone to duplication at this location. The existence of a duplicate tRNAPhe between CR2 and tRNALeu in the viperid O. okinavensis  suggests that frequent gene duplication adjacent to the CRs may occur (these two tRNAPhe differ by only 3 of 64 bp; implying either concerted evolution or recent duplication). The concerted evolution of these tRNAs could be explained by a tendency for gene conversion events involving the duplicate control regions to extend into the homologous flanking tRNA regions.
Another point of interest concerning gene conversion that arises from this study is a preliminary indication of differential evolutionary processes operating on the CRs within versus between species. Vertebrate mitochondrial control regions typically evolve very rapidly, and this is the case in a comparison of the two viperid species (O. okinavensis and A. piscivorus) in which CRs from these species are (on average) approximately as divergent as the fastest evolving positions within the mt genome, third codon positions (Figure 2B). In contrast, the two A. pisvicorus genomes, Api1 and Api2, have surprisingly similar CRs between individuals (Figure 2A; Additional file 2), comparable to the similarity between rRNA genes, among the slowest evolving regions in the mt genome. A previous study on viperid snakes also showed slow within-species CR evolutionary rates , and other studies have demonstrated particularly slow intra-species rates and differential rates of CR evolution operating within versus between species in birds  and fish .
In this study we have found a great deal of rate heterogeneity among genes, so it is certainly possible that the normally unconserved control regions have become suddenly critical and conserved in A. piscivorus. Alternatively, it is plausible that the complex (and poorly understood) process of gene conversion of CRs within a genome may also alter rates of CR evolution within species through a yet unknown process of gene conversion that may involve intragenomic (or even intergenomic) recombination.
Comparative rates of molecular evolution
Previous studies have suggested that snake mt genomes have an accelerated rate of evolution [7, 8]. Our results suggest this general conclusion is an oversimplification of a much more complex scenario, and that rates of snake mt genome evolution incorporate broad temporal (branch-specific) and spatial (gene and gene region-specific) dynamics. Branches early in snake evolution appear to be associated with dramatically elevated evolutionary rates and extreme relative rate dynamics across the mt genome (Figure 4). In contrast, terminal branches appear to have patterns of mt genome evolution that are strikingly similar to other (non-snake) vertebrates.
In support of a hypothesis involving selection for overall oxidative metabolic function, the accelerated rates of molecular evolution in snakes appear to depend greatly on gene function, with most ND subunits accelerating only slightly and occasionally, while COX, ATP, CytB, and rRNA evolutionary accelerations are dramatic and punctuated. The roles of these proteins (and the mitochondria in general) in energetics via oxidative phosphorylation are well known, and it may be that a single causative agent accompanying the diversification of snakes that dramatically altered metabolic demand, or led to a fluctuation in metabolic demand, was responsible for large-scale changes in selective pressure on these proteins.
Snake mitochondrial genomes present a rare opportunity to investigate the evolutionary interactions and ramifications that link genome architecture, molecular evolution, and multi-level molecular function. Available evidence points to selective pressures acting at many hierarchical levels within snake mt genomes, and at different times during snake evolution, leading to diverse, dramatic, and broad-scale changes in the genome. Interestingly, some consequences of this adaptive shift appear to have diminished over time (e.g., accelerated evolutionary rates of COX and other genes), whereas others appear to continue in extant snakes (i.e., the effects of control region duplication on mutation gradients, replication, and potentially transcription, and remnant functional consequences of short and highly substituted genes). Although the precise causes are unknown, this outstanding example of an apparent punctuated adaptive shift involving multiple aspects of genome architecture evolution provides an important comparative tool for the study of vertebrate mt genome evolution.
Sampling, sequencing and annotation
DNA was extracted from vouchered specimens available at the Louisiana State University Museum of Natural Science (LSUMZ) and the University of Central Florida (CLP). The Agkistrodon piscivorus (cottonmouth or water moccasin; Viperidae) specimens were from Louisiana, USA (LSUMZ-17943) and from Florida, USA (CLP-73). We refer to these as Api1 (Louisiana specimen) and Api2 (Florida specimen). The Pantherophis slowinskii (corn snake; Colubridae) specimen was from Louisiana, USA (LSUMZ- H-2036). The genus Pantherophis  was recently erected to contain a clade of species formerly allocated to Elaphe. The species P. slowinskii was formerly considered Pantherophis (Elaphe)guttatus, and was recently recognized as a distinct species . Details of molecular laboratory methods (e.g., PCR, cloning, sequencing), genome annotation , and accession numbers are provided in Additional files (see Additional files 2 and 4).
Phylogenetic and sliding-window analyses
In addition to the three new snake mt genome sequences, the sequence dataset used included all eight snake mt genomes available at the time of the study, and 42 additional taxa for comparative purposes, including heavy sampling of birds, mammals (mostly primates), and lizards (species scientific names and access numbers are given in Additional file 2). Sequences of protein-coding and rRNA genes were aligned using ClustalX , followed by manual adjustment. Protein-coding genes were first aligned at the amino acid level, and then the nucleotide sequences were aligned according to the corresponding amino acid alignment. The alignment of rRNAs contained a small number of sites (corresponding to the loop-forming structures of the rRNAs) with somewhat ambiguous alignments only among major tetrapod lineages. Since we wanted to compare estimates of mitochondrial gene evolutionary rates and patterns, we chose not to exclude any sites of the alignment. This was also justified by preliminary phylogenetic estimates that suggested the incorporation of these few potentially ambiguous sites did not affect phylogenetic results. The main phylogeny presented here was inferred using the concatenated nucleotide sequence of all 13 protein-coding and two rRNA genes by maximum-likelihood (ML) analysis in PAUP 4.0 beta10 . This analysis used the GTR + Γ + I model of evolution, the best-fit model under all criteria in ModelTest .
Support for this topology was evaluated in two ways: (1) based on 1000 NJ bootstraps (in PAUP) with ML distances calculated under the same model as above, but with down-weighted synonymous sites to avoid saturation problems (rRNAs relative weight = 5 and 1st, 2nd, and 3rd codon positions relative weights = 4, 5, and 1) and (2) based on Bayesian posterior probability support estimated by conducting two simultaneous independent MCMC runs conducted for 106 generations (with the first 400,000 generations of each run discarded as burn-in) using a GTR + Γ + I model of evolution (in MrBayes 3.1 ). The burnin period was determined by visual assessment of stationarity and convergence of likelihood values between the chains. To analyze nucleotide substitution rate variation in different lineages and different genes, branch length estimates were separately calculated under the GTR + Γ + I model for different genes (COX1, ND1, ND2, ND4, ND5, CytB) and gene clusters (COX2 + ATP8 + ATP6, and COX3 + ND3 + ND4L; each comprising groups of individually short genes adjacent along the mt genome) using the ML topology and PAML . To further analyze fluctuations in nucleotide substitution rates, we conducted sliding window analyses (SWA) on the phylogenetic dataset. The program Hyphy  was used to estimate branch lengths (estimated numbers of substitutions) for 1000 bp windows. SWA was conducted using the GTR model with global parameter estimation and topological relationships specified based on the ML tree estimate, with a window slide of 200 bp. Based on preliminary trials, the size of the window and slide length were chosen to minimize noise observed with shorter windows, but to allow differentiation of patterns in different regions. To compare patterns of substitution across the mitochondrial genome for select branches or groups of branches, we first divided substitution estimates for each window by the median substitution rate across all windows. Since branch lengths are estimates of δ b t b (the branch-specific substitution rate times divergence time) this procedure estimates a ratio of substitution rates, , where is the branch- and window-specific substitution rate, and is the branch-specific substitution rate in the median window. To evaluate whether the windows had relative rates that were slower or faster than expected, we took the substitution rate ratio from the set of all branches in the non-snakes (NS) as a standard. This was then subtracted from the branch-specific ratio to obtain a "standardized substitution rate",. When relative rates of substitution are distributed similarly across the mt genome, in comparison with NS, this standardized rate comparison approaches zero.
The secondary structures of squamate tRNAs were determined under the guidance of the mammalian tRNA cloverleaf structures  and the tRNAscan program , and then used to modify tRNA alignments by hand (tRNASer [AGY] was not included in these analyses because it does not form a cloverleaf structure). To determine the relative stabilities of the tRNA secondary structures, we calculated the energy (ΔG) of the cloverleaf structure using the Vienna Package version 1.4 .
Analysis of control region functionality
The calculation of T AMS differs depending on whether CR1 or CR2 is functional, but only for the genes that are positioned between the two control regions, the two rRNAs and ND1 (see Additional file 2). Based on previous work, the light strand C/T ratio at synonymous two-fold and fourfold redundant 3rd codon positions is expected to increase linearly with T AMS , so we used this prediction to determine whether there was any evidence for activity of CR1 or CR2 in initiating heavy strand replication. We implemented a slightly modified version of the MCMC approach in  to estimate the most likely slope and intercept of the C/T ratio gradient depending on the calculated T AMS at every site. We applied these calculations using T AMS from CR1 and CR2, and also separately calculated the slope and intercept for the most likely weighted average T AMS for the two control regions. Other than the addition of the weighting parameter, all details of the Markov chain were as in . Relative support for alternative hypotheses was determined using Akaike Information Criterion (AIC) and Akaike weights [40, 41].
tRNA : ribosomal RNA, transfer RNA
origin of heavy strand replication
CR1, CR2: control region, control region 1, control region 2
origin of light strand replication
NADH dehydrogenase subunit #
Cytochrome C oxidase subunit #
Duration of time spent single-stranded by the heavy strand during replication
Time spent in an asymmetric mutagenic state during replication
T, A, G: cytosine, thymine, adenine, guanine
ATP synthase subunit #
Met, Pro, Thr, Leu, Phe, Ser: isoleucine, methionine, proline, threonine, leucine, phenylalanine, serine
sliding window analysis
million years ago
Louisiana State University Museum of Natural Science specimen tag
University of Central Florida specimen tag
Api2: Agkistrodon piscivorus specimen #
We thank Sameer Raina for modifying his program on Bayesian analysis of mitochondrial genome gradients to apply to this project, Jeremiah Faith for help with the preliminary annotations, Wanjun Gu for running some gradient analyses for us, and Judith Beekman for critical comments on the manuscript. This work was primarily supported by grants to D.D.P. from the State of Louisiana Board of Regents (Research Competitiveness Subprogram LEQSF (2001-04)-RD-A-08) and to C.L.P from a UCF startup package and a National Science Foundation grant (DEB-0416000). It was also in part supported by grants to D.D.P. from the National Institutes of Health (GM065612-01 and GM065580-01), the National Science Foundation through Louisiana EPSCOR and the Center for Biomodular Multi-scale Systems, and the State of Louisiana Board of Regents Millennium Research Program's Biological Computation and Visualization Center and Governor's Biotechnology Initiative.
- Krishnan NM, Raina SZ, Pollock DD: Analysis of among-site variation in substitution patterns. Biol Proced Online. 2004, 6: 180-188. 10.1251/bpo88.PubMed CentralView ArticlePubMedGoogle Scholar
- Krishnan NM, Seligmann H, Raina SZ, Pollock DD: Detecting gradients of asymmetry in site-specific substitutions in mitochondrial genomes. DNA Cell Biol. 2004, 23 (10): 707-714. 10.1089/dna.2004.23.707.PubMed CentralView ArticlePubMedGoogle Scholar
- Raina SZ, Faith JJ, Disotell TR, Seligmann H, Stewart CB, Pollock DD: Evolution of base-substitution gradients in primate mitochondrial genomes. Genome Res. 2005, 15 (5): 665-673. 10.1101/gr.3128605.PubMed CentralView ArticlePubMedGoogle Scholar
- Faith JJ, Pollock DD: Likelihood analysis of asymmetrical mutation bias gradients in vertebrate mitochondrial genomes. Genetics. 2003, 165 (2): 735-745.PubMed CentralPubMedGoogle Scholar
- Reyes A, Gissi C, Pesole G, Saccone C: Asymmetrical directional mutation pressure in the mitochondrial genome of mammals. Mol Biol Evol. 1998, 15 (8): 957-966.View ArticlePubMedGoogle Scholar
- Tanaka M, Ozawa T: Strand asymmetry in human mitochondrial-DNA mutations. Genomics. 1994, 22 (2): 327-335. 10.1006/geno.1994.1391.View ArticlePubMedGoogle Scholar
- Dong S, Kumazawa Y: Complete mitochondrial DNA sequences of six snakes: Phylogenetic relationships and molecular evolution of genomic features. J Mol Evol. 2005, 61 (1): 12-22. 10.1007/s00239-004-0190-9.View ArticlePubMedGoogle Scholar
- Kumazawa Y, Ota H, Nishida M, Ozawa T: The complete nucleotide sequence of a snake (Dinodon semicarinatus) mitochondrial genome with two identical control regions. Genetics. 1998, 150 (1): 313-329.PubMed CentralPubMedGoogle Scholar
- Kumazawa Y, Ota H, Nishida M, Ozawa T: Gene rearrangements in snake mitochondrial genomes: Highly concerted evolution of control-region-like sequences duplicated and inserted into a tRNA gene cluster. Mol Biol Evol. 1996, 13 (9): 1242-1254.View ArticlePubMedGoogle Scholar
- Kumazawa Y: Mitochondrial DNA sequences of five squamates: phylogenetic affiliation of snakes. DNA Res. 2004, 11 (2): 137-144. 10.1093/dnares/11.2.137.View ArticlePubMedGoogle Scholar
- Douglas D, Janke A, Arnason U: A mitogenomic study on the phylogenetic position of snakes. Zool Scr. 2006, 35 (6): 545-558. 10.1111/j.1463-6409.2006.00257.x.View ArticleGoogle Scholar
- Reyes A, Yang MY, Bowmaker M, Holt IJ: Bidirectional replication initiates at sites throughout the mitochondrial genome of birds. J Biol Chem. 2005, 280 (5): 3242-3250. 10.1074/jbc.M411916200.View ArticlePubMedGoogle Scholar
- Pesole G, Gissi C, De Chirico A, Saccone C: Nucleotide substitution rate of mammalian mitochondrial genomes. J Mol Evol. 1999, 48 (4): 427-434. 10.1007/PL00006487.View ArticlePubMedGoogle Scholar
- Ingman M, Kaessmann H, Paabo S, Gyllensten U: Mitochondrial genome variation and the origin of modern humans. Nature. 2000, 408 (6813): 708-713. 10.1038/35047064.View ArticlePubMedGoogle Scholar
- Gower DJ, Vidal N, Spinks JN, McCarthy CJ: The phylogenetic position of Anomochilidae (Reptilia: Serpentes): first evidence from DNA sequences. Journal of Zoological Systematics and Evolutionary Research. 2005, 43 (4): 315-320. 10.1111/j.1439-0469.2005.00315.x.View ArticleGoogle Scholar
- Lawson R, Slowinski JB, Crother BI, Burbrink FT: Phylogeny of the Colubroidea (Serpentes): New evidence from mitochondrial and nuclear genes. Mol Phylogenet Evol. 2005, 37 (2): 581-601. 10.1016/j.ympev.2005.07.016.View ArticlePubMedGoogle Scholar
- Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19 (12): 1572-1574. 10.1093/bioinformatics/btg180.View ArticlePubMedGoogle Scholar
- Fernandez-Silva P, Enriquez JA, Montoya J: Replication and transcription of mammalian mitochondrial DNA. Exp Physiol. 2003, 88 (1): 41-56. 10.1113/eph8802514.View ArticlePubMedGoogle Scholar
- Clayton DA: Replication of animal mitochondrial DNA. Cell. 1982, 28 (4): 693-705. 10.1016/0092-8674(82)90049-6.View ArticlePubMedGoogle Scholar
- Inoue JG, Miya M, Tsukamoto K, Nishida M: Evolution of the deep-sea gulper eel mitochondrial genomes: Large-scale gene rearrangements originated within the eels. Mol Biol Evol. 2003, 20 (11): 1917-1924. 10.1093/molbev/msg206.View ArticlePubMedGoogle Scholar
- Sano N, Kurabayashi A, Fujii T, Yonekawa H, Sumida M: Complete nucleotide sequence of the mitochondrial genome of Schlegel's tree frog Rhacophorus schlegelii (family Rhacophoridae): duplicated control regions and gene rearrangements. Genes Genet Syst. 2005, 80 (3): 213-224. 10.1266/ggs.80.213.View ArticlePubMedGoogle Scholar
- Abbott CL, Double MC, Trueman JWH, Robinson A, Cockburn A: An unusual source of apparent mitochondrial heteroplasmy: duplicate mitochondrial control regions in Thalassarche albatrosses. Mol Ecol. 2005, 14 (11): 3605-3613. 10.1111/j.1365-294X.2005.02672.x.View ArticlePubMedGoogle Scholar
- Eberhard JR, Wright TF, Bermingham E: Duplication and concerted evolution of the mitochondrial control region in the parrot genus Amazona. Mol Biol Evol. 2001, 18 (7): 1330-1342.View ArticlePubMedGoogle Scholar
- Amer SAM, Kumazawa Y: Mitochondrial genome of Pogona vitticepes (Reptilia; Agamidae): control region duplication and the origin of Australasian agamids. Gene. 2005, 346: 249-256. 10.1016/j.gene.2004.11.014.View ArticlePubMedGoogle Scholar
- Kumazawa Y, Endo H: Mitochondrial genome of the Komodo dragon: Efficient sequencing method with reptile-oriented primers and novel gene rearrangements. DNA Res. 2004, 11 (2): 115-125. 10.1093/dnares/11.2.115.View ArticlePubMedGoogle Scholar
- Ashton KG, de Queiroz A: Molecular systematics of the western rattlesnake, Crotalus viridis (Viperidae), with comments on the utility of the D-loop in phylogenetic studies of snakes. Mol Phylogenet Evol. 2001, 21 (2): 176-189. 10.1006/mpev.2001.1013.View ArticlePubMedGoogle Scholar
- Crochet PA, Desmarais E: Slow rate of evolution in the mitochondrial control region of gulls (Aves: Laridae). Mol Biol Evol. 2000, 17 (12): 1797-1806.View ArticlePubMedGoogle Scholar
- Tang Q, Liu H, Mayden R, Xiong B: Comparison of evolutionary rates in the mitochondrial DNA cytochrome b gene and control region and their implications for phylogeny of the Cobitoidea (Teleostei: Cypriniformes). Mol Phylogenet Evol. 2006, 39 (2): 347-57. Epub 2005 Oct 4.. 10.1016/j.ympev.2005.08.007.View ArticlePubMedGoogle Scholar
- Utiger U, Helfenberger N, Schätti B, Schmidt C, Ruf M, Ziswiler V: Molecular systematics and phylogeny of Old World and New World ratsnakes, Elaphe Auct., and related genera (Reptilia, Squamata, Colubridae). Russian Journal of Herpetology. 2002, 9 (2): 105-124.Google Scholar
- Burbrink FT: Phylogeographic analysis of the cornsnake (Elaphe guttata) complex as inferred from maximum likelihood and Bayesian analyses. Mol Phylogenet Evol. 2002, 25 (3): 465-476. 10.1016/S1055-7903(02)00306-8.View ArticlePubMedGoogle Scholar
- Slack KE, Janke A, Penny D, Arnason U: Two new avian mitochondrial genomes (penguin and goose) and a summary of bird and reptile mitogenomic features. Gene. 2003, 302 (1-2): 43-52. 10.1016/S0378111902010533.View ArticlePubMedGoogle Scholar
- Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG: The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 1997, 25 (24): 4876-4882. 10.1093/nar/25.24.4876.PubMed CentralView ArticlePubMedGoogle Scholar
- Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony (* and Other Methods). 1997, Sunderland, Massachusetts , Sinauer AssociateGoogle Scholar
- Posada D, Crandall KA: ModelTest: testing the model of DNA substitution. Bioinformatics. 1998, 14 (9): 817-818. 10.1093/bioinformatics/14.9.817.View ArticlePubMedGoogle Scholar
- Yang ZH: PAML: a program package for phylogenetic analysis by maximum likelihood. Comput Appl Biosci. 1997, 13 (5): 555-556.PubMedGoogle Scholar
- Pond SLK, Frost SDW, Muse SV: HyPhy: hypothesis testing using phylogenies. Bioinformatics. 2005, 21 (5): 676-679. 10.1093/bioinformatics/bti079.View ArticlePubMedGoogle Scholar
- Helm M, Brule H, Friede D, Giege R, Putz D, Florentz C: Search for characteristic structural features of mammalian mitochondrial tRNAs. RNA. 2000, 6 (10): 1356-1379. 10.1017/S1355838200001047.PubMed CentralView ArticlePubMedGoogle Scholar
- Lowe TM, Eddy SR: tRNAscan-SE: A program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997, 25 (5): 955-964. 10.1093/nar/25.5.955.PubMed CentralView ArticlePubMedGoogle Scholar
- Hofacker IL, Fontana W, Stadler PF, Bonhoeffer LS, Tacker M, Schuster P: Fast folding and comparison of RNA secondary structures. Monatshefte Fur Chemie. 1994, 125 (2): 167-188. 10.1007/BF00818163.View ArticleGoogle Scholar
- Akaike H: Information theory and an extension of the maximum likelihood principle. Second International Symposium on Information Theory. Edited by: Petrov BN, Csake F. 1973, Budapest , Akademia Kiado, 673-681.Google Scholar
- Akaike H: Information measures and model selection. Int Stat Inst. 1983, 22: 277-291.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.