In the current study we obtained new material for eight Arvicolinae species (Alticola lemminus Miller, 1898, Chionomys nivalis Martins, 1842, Terricola subterraneus de Selys-Longshamps, 1836, Dicrostonyx torquatus Pallas, 1778, Lemmus sibiricus Kerr, 1792, Myopus schisticolor Lilljeborg, 1844, Ondatra zibethicus Linnaeus, 1766, and Prometheomys schaposchnikowi Satunin, 1901), see details for sampling localities in Biosamples description. Tissue mix (muscles, liver, heart, lungs, and testes for males) was fixed in an intactRNA buffer (Evrogen, Moscow, Russia) to avoid RNA degradation directly at the moment of catching the animal in the field. The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of Zoological Institute Russian Academy of Sciences (permission #2-17/14-07-2022).

For nine species, SRA data were downloaded from the NCBI SRA database. Additionally, two already assembled transcriptomes of Mynomes ochrogaster Wagner, 1842 and Cricetulus griseus Milne-Edwards, 1867 were included in the analysis. Thus, in total, we analyzed transcriptomes of 19 Cricetidae species (see Supplementary Table 1 for details).

We categorized species as either subterranean or above ground based on their behavior and morphological traits by classical description of subterranean rodents (Lacey, 2000; Lessa et al., 2008) as previously (Bondareva et al., 2021b). A species was considered subterranean if it is known to spend the most time underground, perform regular digging activities, and come to the surface only incidentally and within the vicinity of burrow openings (Lessa et al., 2008; Smorkatcheva and Lukhtanov, 2014). Two species fit these definitions in our dataset: Ellobius lutescens Thomas, 1897 and Prometheomys schaposchnikowi. Following the data available in literature, Terrícola subterraneus and Lasiopodomys mandarinus Milne-Edwards, 1871 also could be assigned to subterranean or semi-subterranean species (Ognev, 1950; Kryštufek and Shenbrot, 2022). Taking this into account we categorized all four species as subterranean. They are marked by color on Figure 1. We were limited only by these species being unable to consider other subterranean Arvicolinae forms (ex gr. Microtus pinetorum Le Conte, 1830) in the analyses due to the lack of data in the SRA database at the period of analysis and the poorly described ecology habits for several potential subterranean species.

Total RNA was isolated with an RNeasy mini kit (Qiagen) under animal cells/spin protocol. RNA quality was quantified using a Bioanalyzer 2100 (Agilent Genomics) with a minimum of 7.0 for the RNA Integrity Number score. The combined protocol of the NEBNext Poly(A) mRNA Magnetic Isolation Module and the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina¹ was used to isolate polyA-RNA from the total RNA fraction and prepare DNA libraries. Preparation was carried out according to the protocol with the following modifications: (1) sample fragmentation was conducted at 94°C for 10 min; (2) the synthesis of the first chain was hold using the program with the following parameters: 25°C - 10 min, 42°C - 30 min, 70°C - 15 min; (3) DNA purification and concentration was carried out according to the protocol using Ampure XP magnetic particles. Elution was transferred with bidistilled water; (4) the number of cycles in PCR consisted of 10 or 15, depending on the initial RNA concentration.

Sample concentrations were measured on a Qubit fluorometer. The quality of the obtained libraries was checked using the Bioanalyzer 2100 Agilent using the DNA High Sensitivity kit. RNA extraction and library preparation were performed using resources of the Skoltech Genomics Core Facility².

Sequencing was performed on the Illumina HiSeq4000 with read length 75 bp and paired-end tracks. Demultiplexing and conversion of data to the fastq format was performed using the bcl2fastq2 Illumina software. Sequencing was performed using resources of the Skoltech Genomics Core Facility (see text footnote 2).

The quality of the obtained paired readings was assessed by the FastQC (Andrews, 2010). Trimmomatic 0.39 (Bolger et al., 2014) was used to remove low quality reads and Illumina adapters. We took only reads with a quality higher than 28–29.

Transcriptomes were assembled with a standard Trinity package (Haas et al., 2013). Base transcriptome statistics were also estimated with Trinity Stats. The ORFs in assembled transcriptomes were predicted by Transdecoder³. The resulting sequences were purified from chimeric ones using the DIAMOND program (Buchfink et al., 2015) with the reference from the NCBI protein base (nr). Transcriptome completeness was evaluated by BUSCO (Manni et al., 2021). Orthologous genes were determined using Proteinortho (Lechner et al., 2011). We searched for single-copy orthologs that were present in all analyzed transcriptomes for future analyzes.

Separate orthological genes were aligned using prank (Löytynoja, 2014) with the correction for the protein-coding genes tripletity. Aligned orthologs were concatenated into one meta-alignment with Python-based script.

Phylogenetic reconstruction for 19 species was based on found orthological genes and carried out in the IQ-tree⁴. Conserved blocks from multiple alignments were previously removed by Gblocks (Castresana, 2000).

Phylogenetic trees used for tracing signatures of the selection pressure on each subterranean species separately (Figure 2) were obtained with MrBayes 3.2.2 (Ronquist et al., 2012) using seven nuclear genes (6,421 bp), see methods in Bondareva et al. (2021a). The following analysis parameters were: nst = mixed and the distribution of the substitution rates between sites, the dataset was divided into partitions by genes. Each analysis was started with a random tree and had two replicates with four Markov chains (MCMC) and 1 million generations, with the results recorded every hundredth generation. Stationarity and convergence of separate runs were assessed using ESS statistics in Tracer v1.6 (Rambaut et al., 2014). The trees were visualized using the FigTree v1.6 program (accessed on 26 November 2021)⁵.

The evaluation of both convergent and parallel amino acid substitutions was carried out by the ProtParCon program (Yuan et al., 2019) with a further search for amino acids characteristic only for subterranean rodents: we looked for replacements occurring at least in two subterranean taxa using a script on Python 3. The significance of the substitution frequency was estimated with the Evolver program included in the PAML4 (Yang, 2007) using Monte Carlo simulation algorithm to obtain event probability of the substitution in each position. Simulated sequences evolving according to the following parameters: an evolutionary model, the branching pattern and branch lengths of the tree, amino acid frequencies and sequence length. Additionally Holm multiple adjustment used to obtain p-value in R software v.3.4.4 (R Core Team, 2020). Detected significant amino acid substitutions were subsequently checked for influencing protein structure by the Sift program⁶ and modeled in SWISS-MODEL (Waterhouse et al., 2018). Visualization and comparison of modeled structures was carried out in PyMol 2.2.0 (Schrödinger and DeLano, 2020).

To test for possible convergent occurrence of the detected substitutions in other subterranean rodents, we used protein sequences of homologous genes from NCBI Protein database (Heterocephalus glaber—XP_004848222.1, XP_021097910.1; Nannospalax galili—XP_008820604.1, XP_008852875.1; Fukomys damarensis—XP_010611512.1, XP_010637758.1).

Selective pressure was estimated as variation in the levels of non-synonymous (dN) and synonymous (dS) substitutions, as well as ω (dN/dS ratio) using several approaches. We applied the classic codeml program (Yang et al., 2000; Yang, 2007) with the comparisons between free-branch and M0, and free-branch and neutral-branch models with the ete-toolkit interface (Huerta-Cepas et al., 2016).

A branch models assumes significant differences between ω values for marked branches (foreground or ωfrg) and the rest branches of the tree (background branches or ωbkg). For branch analysis the total set of species was subdivided into several subsets so that each subset contained a subterranean species and its phylogenetically closest or sister aboveground species (Figure 2). For each subterranean species analysis was implemented with free-branch model (b_free, where ωfrg, and ωbkg are free), neutral-branch model (b_neut, where ωfrg is fixed to one) and M0 model, where all branches evolve at the same rate. The comparison between free-branch and M0 shows whether foreground branches have ω significantly different from the rest of the tree. The comparison between free-branch and neutral-branch models detects if the ωfrg value is significantly higher than 1. The values 999 and 0.001 were regarded as errors. The calculated p-values of all likelihood-ratio tests (LRTs) were corrected using Holm multiple adjustment in R software v.3.4.4.

In addition, we calculate a selection level for each species by free-ratio method. Free-ratio method does not imply the selection of branches and calculates across the whole three giving a unique value for each species. The Wilcoxon test was used to compare the differences between the mean values of dN/dS ratio for subterranean and terrestrial species groups. Test was conducted using R v.3.4.4.

Newly obtained raw sequence RNA-data for nine Arvicolinae species listed above (see “section 2 Materials and Methods,” Supplementary Table 1) are available at the NCBI Sequence Read Archive: SRR21996854—SRR21996864, SRR21996867—SRR21996869 (BioProject No. PRJNA590630).

After the de novo assembly, from 34,750 to 371,688 “Trinity genes” per species were generated (Supplementary Table 2). The N50 length was 1,644–3,571 bp, whereas the mean contig length was about 525 bp. Complete BUSCO values (C in Supplementary Table 2) varied from 48 to 80.1% except already assembled Mynomes ochrogaster and Cricetulus griseus with mean values 98 and 93%, respectively. Since we did not make a functional comparison of transcriptomes, even such a low percentage (48%) of complete BUSCO for several species was sufficient for further analysis.

The search for orthologs by Proteinortho resulted in more than 1,000 orthologous genes for the studied species. But only 112 single-copy orthologs were present in all analyzed transcriptomes. They were used in further analyzes. Such a small number of intersections is explained by our assembly’s incompleteness and limited amount of initial raw data. Since the already assembled and published transcriptomes of Mynomes and Cricetulus were also included in the analysis, we were confident that detected universal single-copy genes were not chimeric or erroneously generated by Trinity annotation.

The total length of the concatenated single-orthologs alignment was 214,696 nucleotides. The 98,595 nucleotides remained after cleaning with the Gblocks (45% of the original alignment). The topology of trees for nuclear (Figure 1) and previously published mitochondrial gene trees (Abramson et al., 2021) practically does not contradict each other. The only difference concerns the position of the taxa in the first basal radiation. According to the data on nuclear genes the earliest split is represented by the Prometheomys schaposchnikowi, whereas the tribe Lemmini is the earliest derivative according to mitochondrial genes (Abramson et al., 2021).

The calculation of free dN/dS ratio model values across the whole tree for each taxa was carried out for each of 112 single-copy orthologs independently. Free dN/dS ratio model values across all genes showed that median dN/dS meanings turned out to be significantly (p-value = 0.002) higher for subterranean rodents as compared to above ground forms: 0.078 and 0.037, respectively.

To test for branch-specific positive selection in each of the four subterranean rodents, we ran a branch model for each of these species independently in comparison only with phylogenetically close terrestrial taxa (Figure 2). For each of 112 single-copy orthologs we failed to find any genes with significant changes in selection level. Despite the observed trend of higher dN/dS values for subterranean species and possible selection relaxation for these genes (Supplementary Table 3), all LRT p-values were insignificant after Holm correction. This tendency is especially pronounced for Ellobius lutescens, where for several genes p-values were significant before Holm correction, in contrast to other subterranean Arvicolinae analysis.

We compared respective sets of substitutions generated by ProtParCon on 112 single-copy orthologous protein sequences in search of parallel molecular selective signatures among the four subterranean species. Parallel amino acid substitutions shared in at least two subterranean species were found in several genes: ERP29, RAD23B, HIKESHI, ZADH2, MRPS14, PYCR2, CCDC86, GTPBP2, SNAPC2, and TTLL12 (Table 1). We additionally checked manually that there was no variability in these positions for above ground rodents sequences. Subsequent verification by the Evolver simulation process showed the significance of non-random occurrence of substitutions in the RAD23 and PYCR2 genes. These genes are associated with DNA repair processes (Pohjoismäki et al., 2012) and response to oxidative stress (Kuo et al., 2016), respectively.

Substitution Ala314Thr in PYCR2 gene was predicted as affecting the protein function with a significant score of 0.00 by Sift program, while Thr143Ala substitution in RAD23 gene was shown as tolerated with a score of 0.81. The RAD23 tertiary structure modeling showed that the replacement site in the protein is located at a non-structural fragment and therefore does not have visible effect on the protein (Figure 3). It was not possible performing such modeling for the PYCR2 gene, because the length of the template turned out to be less than the length of our original sequence. By coincidence, the site of interest was located in the lacking part of the template. Other templates offered by the SWISS-MODEL program were much worse and could not be used.

Search of convergent protein substitution across ancient and highly specialized subterranean species revealed substitution in RAD23 gene of Heterocephalus glaber and Fukomys damarensis. For the same species, we found a substitution in the position 314 of the PYCR2 gene, but with an amino acid different from those found in subterranean Arvicolinae: Val instead of Thr, respectively. We were unable to find similar substitutions in another subterranean rodent subfamily Spalacidae—Nannospalax galii. Amino acids in the studied positions coincided for this species with those in terrestrial Arvicolinae rodents.

Within our transcriptome-based data set we uncover several genes with parallel amino acid substitutions. Verification showed the reliable substitutions in two genes: the RAD23B (Nucleotide Excision Repair Protein) and PYCR2 (pyrroline-5-carboxylate reductase 2). These substitutions were found in both genes in the genome of Ellobius lutescens, in PYCR2 in Prometheomys schaposchnikowi and in the RAD23B in Terricola subterraneus respectively (Table 1).

The found genes were not previously reported in studies of molecular adaptations of subterranean rodents. At the present stage it is impossible directly to link the functional impact of these substitutions to adaptive processes. However, the biochemical pathways and processes in which these genes are involved may play a certain role in the adaptation to subterranean lifestyle. The active work of the mitochondria potentially causes the formation of reactive oxygen forms that represent a damaging factor for the cell (Turrens, 2003). The RAD23B and PYCR2 genes are associated with DNA repair processes (Pohjoismäki et al., 2012) and response to oxidative stress (Kuo et al., 2016), respectively. The RAD23 homolog gene was found in the study of drought adaptation in plants, since clear signs of positive selection were observed (Zhang et al., 2013). The change in expression level of this gene was also revealed in the analysis of cold tolerance of Thinopyrum intermedium (Jaikumar et al., 2020). Moreover, the presence of the same amino acid substitutions in highly specialized subterranean rodent species from the family Bathyergidae may also indicate the importance of these substitutions in the adaptive processes.

The small number of genes with parallel substitutions typical for subterranean forms is not surprising. Even the wide scale study of Kalina Davies revealed (Davies et al., 2018) only 35 genes with parallel substitutions in all analyzed species which was less than 1% of the total number of analyzed genes. Another wide-scale study performed on subterranean rodents (Fang, 2015) discovered parallel substitutions for Heterocephalus glaber and Fukomys damarensis only in two genes: ARG1 and Na(V)1.7 across the whole-genome comparison. However, interesting to note that in the analysis of CYTB gene sequences across more than 60 Arvicolinae species we identify only three parallel substitutions occurring in subterranean forms (Bondareva et al., 2021a).

The previously observed tendency of the selection relaxation in mitochondrial genomes was not fully confirmed for nuclear data. The median dN/dS value is significantly higher across all analyzed genes for subterranean voles compared with above ground using free ratio model. However, branch-model analysis failed to detect genes with significant changes in selection level for individual subterranean lineage compared to its sister one at the tree. Opposite to this in the analogous branch-model analysis of mitochondrial genomes we observed dramatic differences in a number of genes with clear selection signatures in subterranean lineages compared to sister above ground ones (Bondareva et al., 2021b).

This difference between mitochondrial and nuclear genes may be explained both by the slower rate of nuclear genes evolution (Lin and Danforth, 2004) compared to mitochondrial ones and the insufficient number of analyzed genes. It is worth mentioning that Arvicolinae is the youngest group among order Rodentia and even its oldest subterranean lineage, Prometheomys schaposchnikowi, is not older than 7 Ma (Abramson et al., 2021). However even in the very old and highly specialized subterranean groups from different mammal superorders signs of positive selection were revealed only in 10% out of all protein-coding nuclear genes studied (over 8,000) (Davies et al., 2018).

Our screening of nuclear genes showed despite the small number of them in the analysis we nevertheless may see not strong but clear adaptive signs. These signs can be seen in the increase of dN/dS ratio in the free ratio analysis and in the presence of parallel substitutions in subterranean forms. These may from the first look indicate adaptive signals to subterranean life, but we cannot exclude any other explanation for these findings. Analogous results detected in other subterranean rodents were linked to a variation in metabolic rate, body mass, population size, and generation time among lineages (Martin and Palumbi, 1993; Martin, 1995; Bromham et al., 1996). However, we assume that in case of decrease of effective population size, that is characteristic for most subterranean forms, we would observe, in the first turn, the decrease of overall genetic diversity but not the increase of non-synonymous substitutions as we see in our case. Similarly, the other factors listed above may explain the difference in overall genetic diversity between subterranean forms and above ground ones, but could not be responsible for the significant increase in non-synonymous substitutions.

We aimed first to test the hypothesis on convergent molecular signatures of selection in several phylogenetically distant lineages of subterranean Arvicolinae and second, compare if found these signs with published data on other subterranean mammals. According to published data, we can identify two main trends: selection relaxation and occurrence of parallel amino acid substitutions. Combined results of analysis of mitochondrial genomes and transcriptomes showed both these tendencies for subterranean Arvicolinae (Figure 4).

Despite the limited number of nuclear genes in the analysis we have found parallel substitutions for subterranean forms in RAD23B and PYCR2 genes in the current work. Studies performed on representatives of the families Bathyergidae and Spalacidae revealed parallel substitutions in completely different genes, e.g., Na(V)1.7 and ARG1 (Fang et al., 2014b), EPAS1 and AJUBA (Shao et al., 2015). Thirty-five loci with parallel substitutions were common across subterranean pairwise comparisons (Golden moles, African mole-rats, blind mole rats, Star-nosed mole) by Davies et al. (2018). They include FREM1, LAMA3, MFI2, and PIEZO2 (Safran et al., 2010). Additionally, at least three of the loci (SLC26A8, TDRD6, and TEX15) have roles in sperm development, and at least eight (e.g., C3, CLCA1, and TLR4) are involved with an immune response (Safran et al., 2010). As mentioned above, RAD23B and PYCR2 genes were not listed in previous studies, but their function could be associated with adaptive processes. Our own analysis of the protein sequences of these genes in highly specialized Heterocephalus glaber and Fukomys damarensis from Bathyergidae family shows convergent substitution in the same position with Arvicolinae subterranean rodents.

Several previous studies have investigated convergent molecular evolution based on dN/dS ratio estimation in subterranean mammals. Da Silva et al. (2009) found a significant difference for its values in CYTB sequences of subterranean rodents from families Ctenomyidae, Geomyidae and Bathyergidae compared to their aboveground closely related species. Later, signatures of selection were found in several mitochondrial genes of subterranean rodents (Tomasco and Lessa, 2011, 2014; Tavares and Seuánez, 2018). These have revealed convergent evolution in genes relating to hypoxia between plateau zokors (Spalacidae) and naked mole-rats (Bathyergidae), and comparisons of subterranean and closely related terrestrial rodents found significant differences in substitution rates of protein-coding genes between the two groups (Shao et al., 2015). Additionally, two genome-wide screens for changes in evolutionary rates identified enrichment in vision- and skin-related loci (Prudent et al., 2016; Partha et al., 2017). Our previous studies on the Arvicolinae subfamily also showed the relaxation of selection strength in CYTB and mitochondrial protein-coding genes reflecting an increased ω ratio in subterranean lineages compared to aboveground ones (Bondareva et al., 2021a,b). In this work we observed the trend to selection relaxation by free-ratio analysis, but didn’t find individual genes with significant changes in selection level. We can associate it with the fact that Arvicolinae is the youngest group among order Rodentia and even its oldest subterranean lineage, Prometheomys schaposchnikowi, is not older than 7 Ma (Abramson et al., 2021). It is quite likely that this time span is insufficient for the pronounced selection change in slowly evolving nuclear genes. In addition, our sample dataset was limited, which could also affect the result of sensitive branch-based analysis.

The greatest number of selection signs within Arvicolinae subterranean species was found in the genus Ellobius both in mitochondrial genes (Bondareva et al., 2021a,b) and in the nuclear protein-coding genes. This genus is younger than the oldest representative of the subfamily, subterranean P. schaposchnikowi, by at least ca. 2 Ma. The oldest fossil remains of the true Ellobius are known from the Late Pliocene (Lytchev and Savinov, 1974; Zazhigin, 1988; Tesakov, 2004) but the morphological features in all species of the genus are very similar with those observed in the “model” subterranean rodents of the families Bathyergidae and Spalacidae: protruding incisors, very small eyes, and isolation of the mouth region by lips (Gromov and Polyakov, 1977).

So, on RNA-seq data we additionally confirmed our previous conclusion that the number of selection signatures is independent of the evolutionary age of the lineage but fits the degree of specialization to the subterranean niche and adaptive mechanisms to the subterranean lifestyle are convergent subterranean rodents from different families.