Searching for signatures of positive selection in cytochrome b gene associated with subterranean lifestyle in fast-evolving arvicolines (Arvicolinae, Cricetidae, Rodentia)
BMC Ecology and Evolution volume 21, Article number: 92 (2021)
Mitochondrial genes encode proteins involved in oxidative phosphorylation. Variations in lifestyle and ecological niche can be directly reflected in metabolic performance. Subterranean rodents represent a good model for testing hypotheses on adaptive evolution driven by important ecological shifts. Voles and lemmings of the subfamily Arvicolinae (Rodentia: Cricetidae) provide a good example for studies of adaptive radiation. This is the youngest group within the order Rodentia showing the fastest rates of diversification, including the transition to the subterranean lifestyle in several phylogenetically independent lineages.
We evaluated the signatures of selection in the mitochondrial cytochrome b (cytB) gene in 62 Arvicolinae species characterized by either subterranean or surface-dwelling lifestyle by assessing amino acid sequence variation, exploring the functional consequences of the observed variation in the tertiary protein structure, and estimating selection pressure. Our analysis revealed that: (1) three of the convergent amino acid substitutions were found among phylogenetically distant subterranean species and (2) these substitutions may have an influence on the protein complex structure, (3) cytB showed an increased ω and evidence of relaxed selection in subterranean lineages, relative to non-subterranean, and (4) eight protein domains possess increased nonsynonymous substitutions ratio in subterranean species.
Our study provides insights into the adaptive evolution of the cytochrome b gene in the Arvicolinae subfamily and its potential implications in the molecular mechanism of adaptation. We present a framework for future characterizations of the impact of specific mutations on the function, physiology, and interactions of the mtDNA-encoded proteins involved in oxidative phosphorylation.
Understanding mechanisms involved in the formation of adaptations at the molecular level is among the main challenges of evolutionary biology. The occurrence of similar adaptations in pairs of phylogenetically related or distant organisms provides a good opportunity to test hypotheses about the mechanisms of natural selection in evolution at the molecular level.
Subterranean rodents represent a good model for testing hypotheses on adaptive evolution driven by important ecological shifts. They live in a subterranean environment characterized by high levels of carbon dioxide, low levels of oxygen, and a relatively constant temperature and humidity . Being subjected to the drastic change in energy requirements  associated with colonization of the subterranean niche, in particular the transition from oxygen-rich to hypoxic atmospheres [3, 4], it was supposed that selective regimes of the proteins involved in respiration may have experienced positive directional selection in response to their entry into the subterranean habitat .
Voles and lemmings of the subfamily Arvicolinae are remarkable as it is the youngest group among rodents, is rapidly evolving, and is one of the most diverse groups that have colonized almost all landscapes and habitat types in the Northern Hemisphere. Arvicolinae display the fastest documented adaptive radiation among modern mammals. The earliest arvicolids are known from the Late Miocene (ex. gr. 7–8 Ma) both in Eurasia and North America [6, 7]. The subfamily consists of about 150 species grouped according to various opinions into 28–30 genera belonging to 8–10 tribes . The number of extant arvicoline rodents is eight times greater than in the sister taxon Cricetinae, the most recent common ancestors (MRCA) of both are known in the fossil record from the Late Miocene, ca 10 Ma . Arvicolinae emerged during the series of repeated events of rapid speciation with at least three “explosive” periods of rapid divergence during its evolutionary history .
Within this group, at least 5 phylogenetically distant lineages counting altogether nearly 10 species show independent transition to the subterranean lifestyle. The long-clawed mole vole, Prometheomys schaposchnikowi Satunin, 1901 represents the earliest evolutionary lineage among all recent arvicolids and is the only subterranean form among the so-called first radiation of voles . This early split between Prometheomys and all other voles according to molecular dating was estimated around 7 Ma [9, 10]. Other subterranean lineages appear much later and are related to the most speciose and recent radiation wave. Among these lineages are the highly specialized subterranean mole voles of the tribe Ellobiusini with the only genus Ellobius Fisher, 1814 that counts 5 species in two subgenera. Fossil remains of mole voles are known from the boundary of Pliocene–Pleistocene, approximately 2.5 Ma [11, 12], molecular dating estimates of the mole voles lineage split are approx. 4.5–4.8 Ma [9, 13]. Several species from different genera within the most diverse and species-rich tribe Arvicolini (encountering 60–62 species) also show the various extent of adaptation to subterranean lifestyle, particularly: Terricola subterraneus de Selys-Longshamps, 1836, Microtus pinetorum Le Conte 1830 and Lasiopodomys mandarinus Milne-Edwards, 1830. These species belong to different nodes within the tribe  and are not descendants of the MRCA. The sister taxa of each species are surface-dwelling which indicate the independent multiple transitions to subterranean lifestyle within this tribe. The basal radiation of all microtines (including and Lasiopodomys) crown lineages may be estimated as 2.1 Ma so far as according to known fossil record multiple rootless forms assigned to Microtus appear throughout the Northern Hemisphere at 1.9–2.1 Ma [15, 16]. Despite this ecological and phylogenetical diversity, molecular signatures of selection associated with mastering various ecological niches, including subterranean environment are poorly studied.
Mitochondrial genes have often been assumed to be under strong purifying selection because they encode proteins involved in oxidative phosphorylation (OXPHOS) that can directly influence metabolic performance. However, lifestyle shifts that imply an alteration in metabolism might be associated with changes in the selection pressure of those proteins that participate in the biochemical pathways of cellular respiration. Cytochrome b (cytB) is a key component of bc1, one of the protein complexes involved in oxidative phosphorylation in the mitochondrial membrane. It catalyzes the reversible electron transfer from ubiquinol to cytochrome c coupled to proton translocation (Q-cycle [17, 18]). Despite cytB having been extensively used for phylogenetic and phylogeographic studies as a neutral evolutionary marker, many studies provide support for instances of adaptive selection in mammalian mitochondrial protein-coding genes and cytB. After  several papers gave evidence in favor of signatures of positive selection in the evolution of this gene [19,20,21]. Subsequent studies have suggested that, despite strong functional constraints, mtDNA may be subjected to positive directional selection in cases, for example, of energy-demanding lifestyles and/or the limited availability of oxygen [22,23,24]. In turn, Nevo  correlated sequence variation of a portion of the cytB gene with ecological differences between chromosomal races of blind mole-rats Spalax ehrenbergi Nehring, 1897. To date, the influence of these variations on respiratory function remains unclear. Da Silva et al.  found a significantly higher estimated dN/dS ratio (ω)—the ratio of nonsynonymous to synonymous substitutions—in independent lineages of subterranean rodents concerning their non-subterranean counterparts, suggesting a link between the evolution of this gene and the colonization of a hypoxic environment. This observation was later confirmed  using a set of seven complete mtDNA genomes of octodontoids.
As cytB is an important component of the electron transport chain and thus respiratory processes of the cell, several attempts have also been implemented to find adaptive amino acid changes in partial and complete sequences of this gene. Some amino acid changes may improve aerobic capacity and adaptations to new thermal environments [22, 26,27,28,29]. For instance, mutations in mitochondrial genes have been implicated in exercise intolerance in humans . Da Silva et al.  detected many amino acid substitutions in the cytB sequence in the tuco-tucos, coruro, pocket gophers and African mole rats, in both unique and in the same positions for all species.
In this study, we analyzed the between and within species variation of the cytB gene among subterranean vole species, with special reference to signatures of positive selection on the background of a wider group of surface-dwelling rodents from the Arvicolinae subfamily.
Earlier, an acceleration of substitution rates in cytB in subterranean rodents was shown ; however, authors carried out their analysis in a phylogenetic framework at the family level and over a large evolutionary period. Here, we tested the hypothesis of whether the phenomenon  described is true at a significantly smaller evolutionary scale and within a lower taxonomic level. Comparing different genera within subfamily and different species within one genus that independently mastered the subterranean lifestyle, we aimed to determine whether the transition to the subterranean life is followed by positive selection in cytB and whether the substitutions occur at homologous or different sites.
In this study, we search the natural selection traces in subterranean species of the Arvicolinae subfamily. We analyzed sequences of protein-coding mitochondrial gene cytochrome b for 62 species shown on Fig. 1. Our study covers representatives of all main taxonomic groups and genera that allow us to study species with various times of divergence. We obtained a list of TreeSAAP significant substitutions (categories 6–8) for the studied species in a phylogenetic context since the sequences were analyzed with an account of the phylogenetic position of the species on the tree. From the list of significant substitutions, we selected those typical of at least 3 subterranean species. We found three substitutions that satisfied our criteria: Ser57Pro, Asp214Asn and Ile338Val (Fig. 1). The same substitution type Asp214Asn was also found in specialized subterranean rodents belonging to other families.
Comparison of amino acid frequencies per site revealed (amino acid patterns) more than 80 sites being significantly different between subsets of subterranean and non-subterranean species (Additional files 1, 2b). The substitutions 57 and 338 were among these.
The serine to proline substitution at residue 57 in subterranean rodents potentially removes a phosphorylation site. We used two different methods to predict the phosphorylation of this site. NetPhos 3.1 Server predicted phosphorylation with CDC2 kinase with score 0.518, GPS 5.0 predicted AGC, PKN, PKN1 kinases with score 65.363. The predictions of the type of kinase do not agree with each other, however, all predictions assign high probabilities to this site being phosphorylated. The same methods did not predict phosphorylation for Asp214, and to our knowledge, neither Ile nor Val can be phosphorylated.
A single nucleotide polymorphism at position 338 (ATT > GTT) corresponding to Ile338Val substitution was detected as likely pathogenic in the ClinVar NCBI database and is associated with cancer processes: www.ncbi.nlm.nih.gov/clinvar/variation/143898/.
We modeled the structure of cytB to examine the possible effect of three substitutions: Ser57Pro, Asp214Asn, and Ile338Val (Fig. 2a). Based on our model, Ser57 faces the intermembrane space of the mitochondria. It resides on an unstructured loop segment spanning residues 54–60. This loop contacts the same loop on the second cytochrome bc1 monomer in the complex (Fig. 2b). Unlike Ser57Pro, substitution Asp214Asn is located on a loop facing the matrix of the mitochondrion. It contacts the N-terminus of ubiquinol-cytochrome c reductase complex III subunit VII (UQCRQ) (Fig. 2c). Substitution Ile338Val is located on the interface between α-helices in the transmembrane region of the complex (Fig. 2d). The modeled structure shows that this substitution favors a different rotamer of Ile350, which neighbors residue 58 of UQCRQ.
We observe a tendency towards weakening purifying selection in subterranean rodents in branch analysis. All species showed significant differences compared with the neutral model (b_neut ) on branches by LRT (Table 1) that indicate the signature of relaxe selection in cytB gene. Significant differences between foreground branches (subterranean species) and background branches (surface-dwelling species) were indicated for Ellobius species, Lasiopodomys mandarinus and Terricola subterraneus according to LRT results comparing free-ratio (b_free) and one ratio (M0) models. The same result was obtained in the analysis total set of subterranean species.
Branch-site analysis implemented in the Datamonkey server (aBSREL) found no evidence of episodic diversifying selection in analyzed phylogeny. RELAX confirmed changes in natural selection level of subterranean rodents. So, K-values for three species (Ellobius sp., L. mandarinus and P. schaposchnikowi; Table 1) demonstrated values < 1, which could indicate relaxed selection. Two of these species: Ellobius sp. and L. mandarinus, showed increased ω-values in branch analysis. K for T. subterraneus demonstrated a value much greater than 1 that could be interpreted as “selection strength has been intensified and correlated with lower ω-value in comparison with surface-dwelling species.
Furthermore, the calculation of the distribution of nucleotide substitutions at each site of cytB separately and combined by domain coordinates manually showed the following: site analysis indicated four positions (Table 2) with significantly higher values of substitutions in cytB in subterranean species (sites 4, 236, 237, and 241). The search of sites under positive selection, applying sites and branch-sites models, with programs: PAML (site model M2), FEL and MEME, taking into account all species and several clades separately has shown both unique and similar sites in different species (Table 2). The detected positions were further checked for variation in amino acid substitution pattern at intraspecific level (Additional file 1).
Despite programs giving different results in many cases, there are several sites consistent in several analyses: site 241 was indicated both with FEL and manual search on full-species dataset and on separate clades with Ellobius sp. (FEL) and L. mandarinus (PAML). Position 329 was detected by FEL in L. mandarinus and P. schaposchnikowi clades. Site 238 indicated by PAML and MEME analysis in L. mandarinus clade. All these sites (238, 241, and 329) show a significant difference in amino acid patterns.
Comparison of substitutions in whole domains (Table 3, Fig. 4) uncovered significant differences in membrane domains 1, 2, 5, and 9 and transmembrane domains 5 and 7 for nonsynonymous substitutions (Fig. 4a); and membrane domain 6 and transmembrane domain 5 for synonymous substitutions (Fig. 4b). Its visualization (Additional file 2a) shows that the location does not correlate with special positions relative to other domains or complex components.
Purifying selection is a predominant force in the evolution of mtDNA. But it is possible that weak and/or episodic positive selection occurs simultaneously during the shift to a lifestyle with greater energy demands or reduced availability of oxygen [23, 24]. We examined the possibility of it in subterranean lineages bearing in mind that colonization of the subterranean niche results in exposure to environmental changes that may affect the function of mitochondrial genes. Among Arvicolinae, several independent, relatively recent, and non-simultaneous colonization of the subterranean niche can be examined for features consistent with convergent evolution under similar directional selection, using their non-subterranean relatives for comparison.
By implementing the standard tests we show that (1) several phylogenetically distant subterranean species show similar amino acid substitutions in cytochrome b, and these substitutions are plausibly important for the protein complex structure, (2) cytB showed an increased ratio between non-synonymous and synonymous substitutions (ω) in subterranean lineages compared to non-subterranean with evidence that selection strength has been relaxed, and (3) eight protein domains possess increased substitution ratio in subterranean species as well several nucleotide positions. Taken together, these results are consistent with the hypothesis that colonization of the subterranean niche promotes a new selective regime of positive, directional selection in protein-coding mtDNA genes. Below we discuss these findings to compare our results with previous findings on subterranean and surface-dwelling species from other taxa.
Amino acid substitutions in cytB and their impact on the protein structure
Subterranean and surface-dwelling species of Arvicolinae were examined for the presence of amino acid substitutions in the mitochondrial cytB that could separate both. When only one sequence per species was used in TreeSAAP analysis, three sites with similar amino acid substitutions in distantly related subterranean species compared to all of the analyzed non-subterranean species were found, particularly Ser57Pro, Asp214Asn, and Ile338Val (Fig. 1). Using the results reported in , we compared substitution patterns in subterranean species from different families. A substitution at site 57 was also detected in African mole rats (family Bathyergidae), and in 214 in African mole rats and tuco-tucos (genus Ctenomys). An analogous substitution at site 214 was found to be under positive selection in high-altitude subterranean zokors Eospalax fontanierii Milne-Edwards, 1867 . It is remarkable that among subterranean voles, substitution Asp214Asn was found in Prometheomys schaposchnikowi, Ellobius fuscocapillus and Ellobius lutescens and the same amino acid substitution was found in most of the specialized subterranean rodent families (Fig. 1). Among the Arvicolinae, the three mentioned species referred to the lineages with relatively earlier shifts and longer evolution periods compared to other subterranean forms in the subfamily. Prometheomys schaposchnikowi, the oldest lineage within the subfamily, represents the first wave of species radiation within Arvicolinae and is sister to all other recent lineages. The molecular time estimate for these split dates to approximate 7 Ma . The putative origin of Ellobiusini is estimated as ca. 4 Ma, and other subterranean arvicolids, T. subterraneus, M. pinetorum and L. mandarinus belonging to the crown lineages within the Arvicolini tribe could not be earlier than 2 Ma. Thus, the occurrence of the same substitution in highly specialized species that originated at 2–5 Ma intervals, may indirectly indicate the importance of this substitution for adaptation to subterranean life. Also, the substitution at site 57 from serine to proline removes the opportunity for phosphorylation. No substitution was detected at site 338 in Spalacidae or Bathyergidae families but could be associated with the cancer process according to the ClinVar database. The latter may confirm the importance of substitution, but its impact and adaptive significance remain unclear.
The wide distribution of cytB as a phylogenetic marker allowed us to test changes in amino acid usage across the full-length protein. We detected more than 80 positions with significant changes including part of TreeSAAP detected positions. Amino acids Pro and Val are not major in positions 57 and 338, respectively, but used more often. Our models of the structure of cytB suggest that the substitution Ser57Pro can alter the dimerization of the complex since it is located on a loop segment that resides on the interface of the two bc1 monomers (Fig. 2). The functional role of the substitutions Asp214Asn and Ile338Val is less clear. The loop spanning residues 54–60 is rich in residues that potentially can form hydrogen bonds: Ser57, Asp58, Thr59, Thr60, Thr61 in Lemmus sibiricus. The Ser57Pro reduces the number of available hydrogen bonds, potentially weakening the interaction between loops. Additionally, this substitution eliminates a putative phosphorylation site, which can modulate the inter-loop interaction.
Both residues 214 and 338 of cytochrome b are in the vicinity of UQCRQ. In our model of Lemmus sibiricus Kerr, 1792, Asp214 forms an ion pair with an Arg2 of UQCRQ. This UQCRQ residue is conserved across multiple rodents (NCBI ID: RLQ55034.1, AAH28519.1, NP_777230.1, NP_001020305.1), and experimentally obtained structures from different mammals also display this interaction [32, 34]. Curiously, the Asp214Asn substitution in E. lutescens no longer interacts with Arg2 of UQCRQ; however, a co-occurring substitution—Asn212Asp—forms the ion pair instead of Asn214. The Ile338Val substitution may alter the binding of UQCRQ indirectly, through Ile350. Greater insight into the function of these substitutions may come with the gene sequences of UQCRQ for surface-dwelling and subterranean rodents. Taken together, structural modeling suggests the altered dimerization of cytochrome bc1 and UQCRQ binding.
Signs of positive selection in the cytB gene
As OXPHOS is a conservative mechanism that is essential for energy metabolism, purifying selection dominates the evolution of mtDNA . MtDNA-encoded COX subunits and cytB are the most conserved genes in mitochondrial genomes . Due to the functional importance of mitochondrial genes, purifying selection is the dominant force in their evolution; however, weak and/or episodic positive selection may occur in the background of strong purifying selection if selective pressures shift, as might happen when oxygen availability decreases. Given the drastic change in energy requirements  and the transition from oxygen-rich to a hypoxic subterranean habitat [3,4,5], it is likely that genes involved in respiration experienced positive directional selection. Under this hypothesis, accelerated rates of functional (nonsynonymous) substitutions, relative to silent substitutions in these genes, are expected in subterranean organisms .
Signatures of positive selection were identified in a series of branch analyses implemented for a subset consisting of a subterranean species and a group of sisters surface-dwelling taxa and on the complete dataset. Both in comparing four Ellobius species against Arvicola, Eothenomys and Chionomys (Fig. 3) and L. mandarinus with other Lasiopodomys species and Neodon the ω-ratio was found to be significantly different between the foreground branch of subterranean species against the group of non-subterranean taxa. A significant difference can be observed by analyzing T. subterraneus against other Terricola and Microtus species, but the dN/dS ratio is conversely less for subterranean T. subterraneus compared with surface-dwelling. Surprisingly, comparison of P. schaposchnikowi against species of first Arvicolinae radiation (Ondatra, Dicrostonyx, Myopus, and Phenacomys) and M. pinetorum with other Microtus did not yield a significant difference between ω-values. When ω-values are smaller than one, the study of a single gene does not allow us to reject alternative, non-selection explanations for the rate variation, such as a relaxation of purifying selection, variations in metabolic rate, body mass, population size, and generation time between lineages [35,36,37]. However, Da Silva et al.  found a significantly higher ω in the cytB among phylogenetically distant subterranean family-level taxa of rodents (tuco-tucos, coruros, pocket gophers, and mole rats) compared to their above-ground relatives, suggesting a link between positive directional selection in this gene in subterranean lineages. Both the data of RELAX analysis consistent with codeml branch model results and show many sites under positive selection favor the latter suggestion.
Interestingly, the same adaptations in COX and cytB genes were linked to increased evolutionary rates in simian primates [38, 39] and in mammalian species adapted to unusual oxygen requirements [21, 23, 29]. This confirms the adaptive property of the cytB gene. The fact that we could not observe this difference for all species with subterranean lifestyle and less level of ω for T. subterraneus than for phylogenetically close non- subterranean species sites may be related to less pronounced fossorial in the latter species than in highly specialized subterranean forms such as Ellobius and L. mandarinus.
Substitution density across cytB
We examined the substitution frequencies by domains in subterranean and surface-dwelling rodents. Significant differences in frequencies were observed more often with nonsynonymous substitutions than with synonymous ones. This result is consistent with our data on dN/dS estimations (Fig. 4).
According to several studies [40,41,42], the three main structural domains of the cytB protein are characterized by pronounced differences in the levels of amino acid variation. Gering et al.  consider these to be an adaptation to high-altitude in Peromyscus maniculatus, as the matrix domain was the most variable and the intermembrane domain was the least variable. In the taxa that were studied and described here, we came across the opposite situation: from six domains with significant differences in nonsynonymous substitutions, four belong to the membrane (intermembrane) domains and only two to the transmembrane domains. The influence of nonsynonymous substitutions on the protein structure is unclear due to a substitution-compensating mechanism. However, the great difference in substitution ratio (especially nonsynonymous) provides evidence in favor of the relaxation of purifying selection revealed with ω estimations.
Our results indicate the signatures of positive selection in the evolution of mitochondrial DNA in Arvicolinae during colonization of subterranean environments. We observe relaxation of selection in cytB sequence using dN/dS calculation with branch, branch-site, and site models. Also, we detect significant differences in substitution distribution by domain structure and changes in amino acid usage for underground rodents. On top of this, we found similar amino acid substitution among phylogenetically distant subterranean lineages that could affect protein structure. Our data corroborate the recent findings that suggest that the evolution of mitochondrial protein genes, in particular cytB, could be associated with metabolic adaptations to environments with low oxygen availability.
In our analysis, we used 62 cytB sequences from Arvicolinae species that represent all major genera and tribes. Among these, several phylogenetically independent lineages including those adapted to existence in the subterranean environment. Here, we consider all species of the genus Ellobius, monotypic Prometheomys schaposchnikowi, and the species Lasiopodomys mandarinus, Terricola subterraneus and Microtus pinetorum as subterranean. We compared these taxa with non-subterranean representatives of 22 genera: Alexandromys, Alticola, Arvicola, Blanfordimys, Chionomys, Clethrionomys, Craseomys, Dicrostonyx, Eolagurus, Eothenomys, Lagurus, Lasiopodomys, Lemmus, Microtus, Myopus, Neodon, Neofiber, Ondatra, Phenacomys, Synaptomys, Terricola, Volemys and used Cricetulus, Mesocricetus, Peromyscus and Phodopus as an outgroup. The full list of species and GenBank accession numbers are given in Additional file 3.
DNA extraction, amplification and sequencing
Muscle and skin tissue samples were collected between years of expeditions and stored in 96% ethanol at − 20 degrees Celcius in a tissue and DNA collection of the Group of molecular systematics of mammals (Zoological Institute RAS). Genomic DNA was isolated from ethanol-preserved muscle tissues using a standard salt extraction protocol . For the better resolution of the Arvicolinae tree, alongside mitochondrial cytB, we used and seven nuclear genes: breast cancer 1 gene (BRCA1), exon 11; growth hormone receptor gene (GHR), exon 10; a fragment of the lecithin cholesterol acyltransferase gene (LCAT), exons 2–5 and introns 2–4; tumor suppressor protein gene (TP53), exons 5–7 and introns 5–6; interphotoreceptor retinoid-binding protein gene (IRBP); von Willebrand factor gene (vWF), exon 28; and acid phosphatase type V gene (Acp5), exons 2 and 3 and partial coding sequence, were amplified. All the primers used and references where the PCR conditions are given are listed in Table 4. All sequences obtained in this work are marked with bold in Additional file 3.
PCR cleanup was performed using the Omnix kit (Omnix, Russia). PCR products were sequenced in both directions using ABI BigDye version 3.1. Sequences were edited and aligned with Geneious R11 (https://www.geneious.com), we also checked that sequences obtained were coded correctly. Final alignments had the following lengths: cytB 1143 bp, BRCA1 1022 bp, GHR 869 bp, LCAT 607 bp, TP53 949 bp, IRBP 1267 bp, vWF 1251 bp, and Acp5 454 bp. The sequences obtained in the current study were deposited in GenBank (Additional file 3).
The full list of genes and GenBank accession numbers used are provided in Additional file 3. The best-fit of several substitution models for each gene (TVM:G:5 for Acp5, GTR:G:5 for BRCA1, HKY:G:5 for GHR, TN:G:5 for TP53, J3:G:5 for LCAT, J2:G:5 for IRBP and vWF) was assessed using Treefinder  under the corrected Akaike information criterion (AICc). Bayesian analysis based on the concatenated alignment of seven genes (partitioned by gene) was performed in MrBayes 3.2.6 . Each analysis started with random trees and two independent runs with 4 Markov chains Monte Carlo (MCMC) were performed for 5 million generations, with sampling every 1000th generation; the standard deviations of split frequencies were below 0.01, potential scale reduction factors were equal to 1.0. Stationarity and convergence of separate runs was examined in Tracer v1.6 . A consensus tree was constructed based on the trees sampled after the 25% burn-in.
Amino acid substitutions detection
Significant physicochemical amino acid changes between residues in cytB were identified using the modified MM01 model implemented in TreeSAAP v3.2 . Eight categories (1–8) were used to represent the magnitude of radical substitutions, of which categories 6–8 indicate the most radical substitutions, by setting a sliding window of 20 codons and analyzing the properties of 31 amino acids . Significant positive z-scores (categories 6–8, P < 0.001) were accepted as a sign of significant change in function.
The distribution of synonymous and nonsynonymous substitutions was calculated between subterranean and surface-dwelling species in each site separately and combined by domain coordinates. Domain coordinates used accordingly UniProt information of Mus musculus cytochrome b domains: https://www.uniprot.org/uniprot/P00158.
Amino acid patterns for each position were determined using all sequences for selected species to consider intraspecific variation. Altogether more than 6 200 sequences were analyzed: 131 for subterranean species and 6 059 for surface-dwelling. This dataset includes all available in Genbank cytB sequences in August 2020. Amino acid patterns were calculated using the script on Python 3. Statistical test was selected considering unequal sample sizes for all analyses.
The significance of the substitution frequency and amino acids patterns were estimated using Fisher’s exact test and Holm multiple adjustment. All computations were performed in R software v.3.4.4 .
Synonymous and nonsynonymous substitution estimations
Variation in the estimates of dN, dS, and ω (dN/dS ratio) was explored using the codeml approach, as implemented in ete-toolkit  For each subterranean species (or group of species for Ellobius) branch analysis was implemented with free-branch model (b_free, were ωfrg and ωbkg are free), neutral-branch model (b_neut, were ωfrg is fixed to one) and M0 model, where all branches evolve at the same rate. Subdivision into analyzed groups was carried out according to the principle of selection of phylogenetic nearest surface-dwelling taxa for ω comparison and more distant as outgroup. Were calculated likelihood-ratio tests to compare different models. The comparison between free-branch and M0 showed if foreground branches have an ω significantly different from the rest of the tree. And the comparison between free-branch and neutral-branch models detect if the value of ωfrg is significantly higher than 1.
Several programs from DataMonkey web-server (datamonkey.org) was used for search selection signatures: aBSREL (An adaptive branch-site REL test for episodic diversification ), FEL (Fixed Effects Likelihood ) and MEME (Mixed Effects Model of Evolution ). Also, we performed PAML  site analysis with M2 model for looking sites under positive selection .
We modeled homology-based structures of cytochrome b as part of the bc1 complex from L. sibiricus and E. lutescens. The overall architecture of the bc1 complex varies little in crystal structures from different organisms ranging from yeast to several mammalian species [60, 61], justifying homology modeling. We based the model on the crystal structure and protein sequence of cytochrome bc1 complex from Bos taurus with a resolution of 2.4 Å (Protein Data Bank ID: 1NTM ). First, the homodimeric structure was recovered by exploiting the symmetry of the crystallographic group. Next, we used modeller (release 9.22)  to create structures of the bc1 homodimer by using the sequence of the cytochrome b from the corresponding organism but keeping the other parts of the complex from Bos taurus. We adopted the auto model protocol with default settings. The models were superimposed to 1NMT and substitutions were analyzed visually in PyMOL v.2.0 (Schrödinger, LLC), which was also used to produce figures. Transmembrane regions of the complex were estimated using the OPM web server .
- cytB :
Cytochrome b gene
Million years ago
Most recent common ancestor
Amount of nonsynonymous substitutions
Amount of synonymous substitutions
The ratio of nonsynonymous to synonymous substitutions
Cell division control 2 kinase
Protein kinase N
Protein kinase N1
Corrected Akaike information criterion
Markov chains Monte Carlo
- BRCA1 :
Breast cancer 1 gene
- GHR :
Growth hormone receptor gene
- LCAT :
Lecithin cholesterol acyl transferase gene
- TP53 :
Tumor suppressor protein gene
- IRBP :
Interphotoreceptor retinoid-binding protein gene
- vWF :
Von Willebrand factor gene
- Acp5 :
Acid phosphatase type V gene
Ubiquinol-cytochrome c reductase complex III subunit VII
Buffenstein R. Ecophysiological responses of subterranean rodents to underground habitats. In: Lacey E, Patton J, Cameron G, editors. Life underground: the biology of subterranean rodents, vol. 62. Illinois: University of Chicago Press; 2000. p. 110.
Vleck D. The energy cost of burrowing by the pocket gopher Thomomys bottae. Physiol Zool. 1979;52(2):122–36.
Nevo E. Mosaic evolution of subterranean mammals: Tinkering, regression, progression, and global convergence. In: Subterranean Rodents: News from Underground. Berlin: Springer Berlin Heidelberg. 2007:375–88.
Luo Y, Gao W, Gao Y, Tang S, Huang Q, Tan X, et al. Mitochondrial genome analysis of Ochotona curzoniae and implication of cytochrome c oxidase in hypoxic adaptation. Mitochondrion. 2008;8(5–6):352–7.
Da Silva CC, Tomasco IH, Hoffmann FG, Lessa EP. Genes and ecology: accelerated rates of replacement substitutions in the cytochrome b gene of subterranean rodents. Open Evol J. 2009;3:17–30.
Martin RA. Biochronology of latest Miocene through Pleistocene Arvicolid rodents from the central great plains of North America. Coloquios Paleontol. 2003;1:373–83.
Fejfar O, Heinrich WD, Kordos L, Maul LC. Microtoid cricetids and the Early history of arvicolids (Mammalia, Rodentia). Palaeontol Electron. 2011;14(3):12.
Musser GG. Superfamily muroidea. In: Mammal species of the world. 2005: 894–1531.
Abramson NI, Lebedev VS, Tesakov AS, Bannikova AA. Supraspecies relationships in the subfamily Arvicolinae (Rodentia, Cricetidae): an unexpected result of nuclear gene analysis. Mol Biol. 2009;43(5):834–46.
Steppan SJ, Schenk JJ. Muroid rodent phylogenetics: 900-species tree reveals increasing diversification rates. Huchon D, editor. PLoS ONE. 2017;12(8):e0183070.
Lychev GF, Savinov PF. Late Pliocene lagomorphs and rodents of Kiikbay. Mater Hist Fauna Flora Kazakhstan, Almaty. 1974;6:39–57.
Tesakov AS. Biostratigraphy of Middle Pliocene–Eopleistocene of Eastern Europe (based on small mammals). 2004.
Lebedev V, Bogdanov A, Brandler O, Melnikova M, Enkhbat U, Tukhbatullin A, et al. Cryptic variation in mole voles Ellobius (Arvicolinae, Rodentia) of Mongolia. Zool Scr. 2020;49(5):535–48.
Martínková N, Moravec J. Multilocus phylogeny of arvicoline voles (Arvicolini, Rodentia) shows small tree terrace size. Folia Zool. 2012;61(3–4):254–67.
Tesakov AS, Vangengeim EA, Pevzner MA. Findings of the oldest rootless voles Allophaiomys and Prolagurus in Eastern Europe. Dokl earth Sci. 1999;366:452–3.
Martin RA, Peláez-Campomanes P, Honey JG, Fox DL, Zakrzewski RJ, Albright LB, et al. Rodent community change at the Pliocene-Pleistocene transition in southwestern Kansas and identification of the Microtus immigration event on the Central Great Plains. Palaeogeogr Palaeoclimatol Palaeoecol. 2008;267(3–4):196–207.
Mitchell P. Possible molecular mechanisms of the protonmotive function of cytochrome systems. J Theor Biol. 1976;62(2):327–67.
Trumpower B. The protonmotive Q cycle in mitochondria and bacteria. Crit Rev Biochem Mol Biol. 1994;29(3):165–97.
Tomasco IH, Lessa EP. Two mitochondrial genes under episodic positive selection in subterranean octodontoid rodents. Gene. 2014;534(2):371–8.
Shao Y, Li JX, Ge RL, Zhong L, Irwin DM, Murphy RW, et al. Genetic adaptations of the plateau zokor in high-elevation burrows. Sci Rep. 2015;5:1–11.
Di Rocco F, Parisi G, Zambelli A, Vida-Rioja L. Rapid evolution of cytochrome c oxidase subunit II in camelids (Tylopoda, Camelidae). J Bioenerg Biomembr. 2006;38(5–6):293–7.
Blier PU, Dufresne F, Burton RS. Natural selection and the evolution of mtDNA-encoded peptides: evidence for intergenomic co-adaptation. Trends Genet. 2001;17(7):400–6.
Shen YY, Liang L, Zhu ZH, Zhou WP, Irwin DM, Zhang YP. Adaptive evolution of energy metabolism genes and the origin of flight in bats. Proc Natl Acad Sci U S A. 2010;107(19):8666–71.
Tomasco IH, Lessa EP. The evolution of mitochondrial genomes in subterranean caviomorph rodents: adaptation against a background of purifying selection. Mol Phylogenet Evol. 2011;61(1):64–70.
Nevo E. Mosaic evolution of subterranean mammals: regression, progression, and global convergence. Oxford Univ Press; 1999.
Grossman LI, Wildman DE, Schmidt TR, Goodman M. Accelerated evolution of the electron transport chain in anthropoid primates. Trends Genet. 2004;20(11):578–85.
Ballard JWO, Whitlock MC. The incomplete natural history of mitochondria. Mol Ecol. 2004;13(4):729–44.
Mishmar D, Ruiz-Pesini E, Golik P, Macaulay V, Clark AG, Hosseini S, et al. Natural selection shaped regional mtDNA variation in humans. Proc Natl Acad Sci. 2003;100(1):171–6.
da Fonseca RR, Johnson WE, O’Brien SJ, Ramos MJ, Antunes A. The adaptive evolution of the mammalian mitochondrial genome. BMC Genomics. 2008;9:1–22.
Woolley S, Johnson J, Smith MJ, Crandall KA, McClellan DA. TreeSAAP: selection on amino acid properties using phylogenetic trees. Bioinformatics. 2003;19(5):671–2.
Yang Z, Nielsen R. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol Biol Evol. 2002;19(6):908–17.
Gao X, Wen X, Esser L, Quinn B, Yu L, Yu C-A, et al. Structural basis for the Quinone reduction in the bc 1 complex: a comparative analysis of crystal structures of mitochondrial cytochrome bc 1 with bound substrate and inhibitors at the Q site. Biochemistry. 2003;42(30):9067–80.
Lychev GF, Savinov PE. Late Pleistocene lagomorphs and rodents of Kiikbay. Mater po Istor Fauny i Flory Kazakhstana. 1974;6:39–56.
Wu M, Gu J, Guo R, Huang Y, Yang M. Structure of mammalian respiratory supercomplex I1III2IV1. Cell. 2016;167(6):1598-1609.e10.
Bromham L, Rambaut A, Harvey PH. Determinants of rate variation in mammalian DNA sequence evolution. J Mol Evol. 1996;43(6):610–21.
Martin AP, Palumbi SR. Body size, metabolic rate, generation time, and the molecular clock. Proc Natl Acad Sci U S A. 1993;90(9):4087–91.
Martin AP. Metabolic rate and directional nucleotide substitution in animal mitochondrial DNA. Mol Biol Evol. 1995;12(6):1124–31.
Andrews TD, Jermiin LS, Easteal S. Accelerated evolution of cytochrome b in simian primates: adaptive evolution in concert with other mitochondrial proteins? J Mol Evol. 1998;47(3):249–57.
Adkins RM, Honeycutt RL. Evolution of the primate cytochrome c oxidase subunit II gene. J Mol Evol. 1994;38(3):215–31.
Esposti MD, De Vries S, Crimi M, Ghelli A, Patarnello T, Meyer A. Mitochondrial cytochrome b: evolution and structure of the protein. Biochim Btophystca Acta. 1993;1143:243–71.
McClellan DA, Palfreyman EJ, Smith MJ, Moss JL, Christensen RG, Sailsbery JK. Physicochemical evolution and molecular adaptation of the cetacean and artiodactyl cytochrome b proteins. Mol Biol Evol. 2005;22(3):437–55.
Gering EJ, Opazo JC, Storz JF. Molecular evolution of cytochrome b in high- and low-altitude deer mice (genus Peromyscus). Heredity (Edinb). 2009;102(3):226–35.
Miller DN, Bryant JE, Madsen EL, Ghiorse WC. Evaluation and optimization of DNA extraction and purification procedures for soil and sediment samples. Appl Environ Microbiol. 1999;65(11):4715–24.
Lebedev VS, Bannikova AA, Tesakov AS, Abramson NI. Molecular phylogeny of the genus Alticola (Cricetidae, Rodentia) as inferred from the sequence of the cytochrome b gene. Zoolog Scr. 2007;36(6):547–63.
Ohdachi S, Dokuchaev NE, Hasegawa M, Masuda R. Intraspecific phylogeny and geographical variation of six species of northeastern Asiatic Sorex shrews based on the mitochondrial cytochrome b sequences. Mol Ecol. 2001;10(9):2199–213.
Bannikova AA, Sizhazheva AM, Malikova VG, Golenishchev FN, Dzuev RI. Genetic diversity of Chionomys genus (Mammalia, Arvicolinae) and comparative phylogeography of snow voles. Genetika. 2013;49(5):649–64.
Petrova TV, Tesakov AS, Kowalskaya YM, Abramson NI. Cryptic speciation in the narrow-headed vole, Lasiopodomys (Stenocranius) gregalis, (Rodentia: Cricetidae). Zoolog Scr. 2016;45:618–29.
Poux C, Chevret P, Huchon D, De Jong WW, Douzery EJ. Arrival and diversification of caviomorph rodents and platyrrhine primates in South America. Syst Biol. 2006;55(2):228–44.
Jobb G, Von Haeseler A, Strimmer K. TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics. BMC Evol Biol. 2004;4:1–9.
Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61(3):539–42.
Rambaut A, Drummond AJ, Suchard M. Tracer v1. 6. http://beast.bio.ed.ac.uk/Tracer. 2014.
McClellan DA, McCracken KG. Estimating the influence of selection on the variable amino acid sites of the cytochrome b protein functional domains. Mol Biol Evol. 2001;18(6):917–25.
R Core Team. R: A language and environment for statistical computing. 2017.
Huerta-Cepas J, Serra F, Bork P. ETE 3: reconstruction, analysis, and visualization of phylogenomic data. Mol Biol Evol. 2016;33(6):1635–8.
Smith MD, Wertheim JO, Weaver S, Murrell B, Scheffler K, Kosakovsky Pond SL. Less is more: an adaptive branch-site random effects model for efficient detection of episodic diversifying selection. Mol Biol Evol. 2015;32(5):1342–53.
Kosakovsky Pond SL, Frost SDW. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22(5):1208–22.
Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Kosakovsky Pond SL. Detecting individual sites subject to episodic diversifying selection. Malik HS, editor. PLoS Genet. 2012;8(7):e1002764.
Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.
Nielsen R, Yang Z. Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics. 1998;148(3):929–36.
Hunte C, Koepke J, Lange C, Roßmanith T, Michel H. Structure at 2.3 Å resolution of the cytochrome bc1 complex from the yeast Saccharomyces cerevisiae co-crystallized with an antibody Fv fragment. Structure. 2000;8(6):669–84.
Crowley PJ, Berry EA, Cromartie T, Daldal F, Godfrey CRA, Lee D-W, et al. The role of molecular modeling in the design of analogues of the fungicidal natural products crocacins A and D. Bioorg Med Chem. 2008;16(24):10345–55.
Webb B, Sali A. Comparative protein structure modeling using MODELLER. Curr Protoc Bioinforma. 2016;54(1).
Lomize MA, Pogozheva ID, Joo H, Mosberg HI, Lomize AL. OPM database and PPM web server: resources for positioning of proteins in membranes. Nucleic Acids Res. 2012;40(D1):D370–6.
Blom N, Sicheritz-Pontén T, Gupta R, Gammeltoft S, Brunak S. Prediction of post-translational glycosylation and phosphorylation of proteins from the amino acid sequence. Proteomics. 2004;4(6):1633–49.
Xue Y, Liu Z, Cao J, Ma Q, Gao X, Wang Q, et al. GPS 2.1: enhanced prediction of kinase-specific phosphorylation sites with an algorithm of motif length selection. Protein Eng Des Sel. 2011;24(3):255–60.
Authors are grateful to Evgeny Genelt-Yanovskiy for branch analysis demonstration, critical reading of manuscript final version and Stanislav Bondarev for help in amino acid pattern analysis. English language was edited by proofreading service: https://www.proof-reading-service.com.
This study was conducted in the Zoological Institute RAS, under research theme No. AAAA-A19-119020790106-0 (for OVB, TVP and NIA) with financial support from grant from the Russian Foundation for Basic Research (RFBR) N 18-04-00730 (for OVB, TVP and NIA), and the Program of Presidium RAS “Dynamics of Gene Pools in Natural populations.” (for OVB, TVP and NIA). The funding body did not exert influence on the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Amino acid sites with significant difference between subterranean and surface-dwelling rodents.
Analyzed position visualization. A. Domains with a significantly increased frequency of nonsynonymous substitutions. B Sites with significant changes in amino acid usage.
Species and genes used in the study. Sequences received in this work marked bold.
About this article
Cite this article
Bondareva, O.V., Potapova, N.A., Konovalov, K.A. et al. Searching for signatures of positive selection in cytochrome b gene associated with subterranean lifestyle in fast-evolving arvicolines (Arvicolinae, Cricetidae, Rodentia). BMC Ecol Evo 21, 92 (2021). https://doi.org/10.1186/s12862-021-01819-4
- Subterranean lifestyle
- Cytochrome b
- Natural selection