Degrees of convergent evolution in rodent adaptations to arid environments

A multipurpose data set to study convergent evolution

Our goal was to study convergence at two different phylogenetic scales, Rodents and Murinae, bringing together kidney expression and coding sequences data for the kidney. We combined publicly available data (cDNA from published genomes and kidney RNA-seq data) with newly generated RNA-seq data to obtain coding sequence and expression data in a variety of rodents, ensuring a reasonable overlap between the two types of data (12 families, 48 species out of which 33 species with expression data) (phylogeny built from our assemblies in Fig. 2; Supplemental Data S1). Most of our sampling effort was meant to increase the number of Murinae species taken from the wild, but we also included species that are more easily obtained from laboratory colonies or commercial strains (Supplemental Tables S1–S3; Supplemental Methods).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Phylogram showing the species used in the study, made using maximum likelihood on a concatenate of 4065 orthologous gene families. Expression (circles) and coding sequence assemblies (squares) with their source (colors) are indicated. Mus musculus strains (gray box, not included in the 33 species count), Murinae_seq data set (dashed lines box), and Total_seq data set (plain line) are delineated. Families and environments are color-coded, and precipitation is indicated. Transitions to arid environments are indicated by stars together with minimum transition ages (or maximum when not available) for data expression analysis.

We associated each species with a biological status corresponding to its arid (orange) or mesic (green) environment. We determined its geographical distribution area and extracted the corresponding bioclimatic variables (Methods) (Supplemental Table S4; Supplemental Methods). Because an annual average pluviometry can hide large differences between seasons, we decided to use the precipitation of the driest quarter of the year and an environment with <40 mm was defined as arid. We annotated 28 species living in arid environments and presumably adapted to arid environments (Fig. 2). For convenience, below we refer to them as “arid species,” whereas other species are named “mesic species.”

To assess the independence of the transitions to arid environments and date them, we used a semiautomated method, enlarging our annotation to 1898 species. Then we estimated ancestral states along a published rodent phylogenetic tree containing 2260 rodent species (Fabre et al. 2012; Supplemental Methods; Supplemental Data S2–S4; Supplemental Tables S1, S4; Supplemental Fig. S1), which was congruent with our phylogenomic tree presented in Figure 2. For species in our data set and their relatives, we further manually curated the transitions based on additional information from specific literature on the concerned species (i.e., more complete phylogenies by specialists of the considered group considering phylogeography) and bracketed them to distinguish recent from older transitions (minimum dates in Fig. 2, intervals and bibliography in Supplemental Fig. S3; Supplemental Table S5; Supplemental Methods).

This provided two data sets for expression analyses, so-called Total_expr and Murinae_expr, which encompassed 33 and 11 species, respectively (Figs. 2, 3A,E). Not all species for which RNA-seq data were available were selected for these data sets, as the selection needed to balance the many requirements for gene expression analysis (reproducible kidney dissection, RNA quality and batch effect, number of replicates; available expression data for an outgroup mesic species or a closely related species; minimizing overrepresentation of a genus; keeping a reasonable balance between mesic and arid species). In these final data sets, nine selected species had no replicates. They were nevertheless selected because (1) a closely related species from the same genus could serve as biological pseudoreplicates for the considered environmental transition, and/or (2) arid and corresponding mesic species were all dissected and sequenced by us at the same time, ensuring minimal confounding effects.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Detection of convergence in expression levels in rodents and in Murinae. (A–D) Detection in the Total_expr data set. (A) Phylogeny of species in the total expression data set. Species names are colored by environment (mesic: green; arid: orange). Rodent families are color-coded. Minimal dates of transitions to arid environments are indicated. (B) Whole transcriptome patterns extracted with a principal component analysis (PCA) made using batch-corrected expression values for 5443 genes and 80 individuals. PC1 is shown with PC5, because it is best correlated with the environment. Main patterns from a phylogenetic principal component analysis (Phylo-PCA); the components associated with most variance are represented. (C) Coexpression modules are represented by their eigen gene with the color showing expression levels in each species. Modules significantly down- or upregulated in arid species are depicted at the top of the panel, with grayscale colors indicating the strength of the correlation. Modules are named according to the sign and significance of this correlation and a number (tu: Total up; td: Total down; tn: Total not correlated). Barplots represent the numbers of genes in each module. (D). Venn diagram of differentially expressed genes found with DESeq2 and EVEmodel. The position of genes passing the EVEmodel FDR threshold is shown with a small bar. (E–H) Results for the Murinae expression data set. (E–G) Detection in the Murinae data set, for 7779 genes and 29 individuals. Colors as in A, and legend as in B–D. (mu) Murinae up, (md) Murinae down, (mn) Murinae not correlated. Code for species names in Figure 2.

In contrast, for coding sequence analysis, we maximized phylogenetic power by working on the largest possible data sets, named Total_seq and Murinae_seq, which encompassed respectively all the rodent species (48) and all the Murinae species (14).

Higher degree of gene expression convergence at gene, module, and whole-transcriptome level in the Murinae subdata set compared with the total data set

We investigated whether differences in the whole transcriptome that were recorded agnostically by principal component analysis (PCA) or by coexpression module analysis carried a signal associated with the mesic/arid environment. Next, we investigated the possibility of identifying genes that exhibit differential expression in relation to the environment.

First, we analyzed broad patterns of expression differences in a PCA. In the Total_expr data set, the difference between arid and mesic individuals only becomes apparent when the phylogeny is accounted for (Methods) (Fig. 3B; Supplemental Fig. S4). We found that the fifth principal component (PC), which accounted for only 4% of the variation, only partially separates arid from mesic species. We verified that this was not an artifact of technical batch effects (Supplemental Fig. S5; Supplemental Table S3). In the Murinae_expr data set, the first component carries most of the separation between arid and mesic individuals and accounts for a large proportion of the variation (25%) (Fig. 3F). However, on the first two PCs, three individuals of the arid species Mus macedonicus locate with mesic species of the genus Mus, probably because they are closely related and because M. macedonicus lives in the moderately arid mediterranean environment. To investigate the phylogenetic structures in these multivariate comparative data, we then used a phylogenetic principal component analysis (pPCA) (Jombart et al. 2010a,b; Supplemental Methods). We observed that in the Total_expr data set, the third component, which shows a phylogenetic signal, is also correlated with the environment (PC3) (Fig. 3B; Supplemental Fig. S6). In the Murinae_expr data set, the two components with the most variance are correlated with the environment, one with a phylogenetic signal (PC1) but the other in opposition with the phylogenetic signal (PC12) (Fig. 3F; Supplemental Fig. S7). This suggests convergence at the whole-transcriptome level.

By constructing dendrograms based on gene expression distances of DE genes (thus, without phylogenetic correction) we confirm that species adapted to arid environments tend to be closer, but the bootstrap supports are low. Hence, as in the PCA, there is a strong phylogenetic signal mixed with environment signal in expression divergence (Supplemental Figs. S4, S8).

Second, we hypothesized that convergence in gene expression could be detected in functionally related sets of genes that exhibit coordinated changes in expression throughout the phylogeny. We ran correlation network analyses to detect these modules of covarying genes (Horvath 2011). Seven modules in the Total_expr and 10 in the Murinae_expr data set showed a significant correlation with the arid habitat (Fig. 3C,G; Supplemental Table S11). Within each rodent family, the expression level in the Total_expr modules varied consistently with the habitat, confirming that a convergence signal is present alongside the phylogenetic signal (Supplemental Figs. S4, S9). The 10 Murinae_expr modules presented several overrepresented functional categories consistent with the literature mentioned in introduction, such as lipid metabolic processes, cholesterol metabolism, anion transport, and symporter activity (Supplemental Table S10; for top 25 terms, see Supplemental Fig. S10).

Third, we assessed differential expression by pairwise contrasts between arid and mesic conditions. Using the standard tool DESeq2 for differential expression analysis, we accounted for the family as a batch effect to identify genes that were significant at the 0.1-adjusted P-value threshold and exhibited a log-fold change greater than 0.4 (for lists and other thresholds, see Supplemental Tables S6, S7). We found 41 genes in the Total_expr data set and 744 genes in the Murinae data set. This is 30.0 and 31.7 times more genes than expected (P-values < 0.0001 estimated from permutations) (Supplemental Fig. S11). Among these genes, we found genes involved in kidney metabolism and kidney adaptation to aridity (e.g., Slc35b4 and Slc40a1 in Total_expr; Aqp2, Aqp7, Slc14a2, Slc8b1, and Slc27a2 in Murinae_expr) (Supplemental Material; Supplemental Tables S9, S10). In addition, we applied the EVE phylogenetic modeling method, an approach based on an Ornstein–Uhlenbeck model with two optima that allows statistical assessment of branch-specific shifts in expression while accounting for interspecies variation (Rohlfs and Nielsen 2015; Gillard et al. 2021) but works one gene at a time, contrary to DESeq2. After P-value adjustment for multiple testing, we found only two genes in the Total data set and four in the Murinae data set (0.1-adjusted P-value threshold). However, many genes found by DESeq2 are significant with EVE, suggesting that both methods are to some extent congruent (70.1% of DESeq2 significant genes have EVE P-value < 0.05 for Total_expr, 77.6% for Murinae_expr) (Fig. 3D,H; Supplemental Fig. S12A). The remaining genes may reflect different abilities of the two methods to detect phylogenetic structures of convergent expression signals (Fig. 1).

Because the various approaches are consistent with one another, after correcting for the phylogeny, we conclude that the Total_expr rodent data set exhibits a moderate signal for habitat-related convergent evolution of gene expression levels, whereas the Murinae_expr data set displays a stronger signal.

Transcriptome convergence is associated with convergent shifts in kidney cell composition

Tissue bulk RNA-seq data reflect variation in both expression per cell and cell-type composition (abundance). Here, different species may exhibit divergent kidney structures as part of their adaptation to arid environments. We investigated whether tissue composition distinguishes species living in mesic and arid environments by deconvoluting our bulk RNA-seq data with MuSiC (Wang et al. 2019). This method uses cell-type-specific gene expression from single-cell RNA-seq data in a reference data set to characterize cell-type compositions from target bulk RNA-seq data (Fig. 4A). We made several controls to choose the best reference data set (a rat data set from Balzer et al. 2023), and to evaluate the limits of the method (Supplemental Figs. S13–S17; Supplemental Methods). We estimate that the method can be used for moderately distantly related species, because markers of cell types are largely conserved between rats, mice, and hamsters.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Convergent tissue composition in species living in an arid environment. Cell proportions were estimated by deconvolution using rat scRNA-seq data. (A) Principle of the method. UMAP redrawn from selected samples from Balzer et al. (2023). All cells from the scRNA-seq data except immune cells were used. (B,D) PCA on cell-type proportions is presented for the Total_expr and Murinae_expr data sets. Centroids are represented for species living in a mesic (green) and arid (orange) environment. (C) Cell proportions for species of the Total_expr data set, segregated by habitat. P-values of linear models testing for an effect of habitat accounting for the rodent family are indicated. (E) Correlations of Murinae coexpression modules (shown in Fig. 3G) with cell proportions, colored according to key. Cell types are abbreviated as follows (see also schema of panel A): endothelial cells (Endo), stroma cells (Stroma), proximal tubule cells (Prox tub), loop of Henle (LOH), distal convoluted tubule (DCT), connecting tubule (CNT), collecting duct principal cells (PC), and collecting duct intercalated cells (IC), immune cells (Immune).

In both the Total_expr and the Murinae_expr data set, the first two components of a PCA calculated on these proportions separated species living in arid/mesic environments (Fig. 4B,D; Supplemental Table S12). This discrimination was significant (discriminant analysis [DA], P-values = 0.001 and 0.002 respectively, Monte-Carlo tests based on 1000 replicates). Visually the separation is noticeable but imperfect in the Total data set but is very clear in the Murinae data set, with the exception of the M. macedonicus samples (which resemble mesic species as seen above). Cell types that mostly contributed to this separation were connecting duct, distal convoluted tubule, stroma, intercalated cells, and endothelial cells on the arid side and proximal tubule and loop of Henle cells on the mesic side. Several cell-type-associated P-values were significant after controlling for the family effect (Fig. 4C, linear models). Notably, for certain cell types, we saw the same pattern across most families (e.g., proximal tubule cells decreased in proportion in species living in arid environments) but not for other cell types (e.g., loop of Henle cells tend to increase in proportion in Cricetidae and Heteromyidae but not Muridae) (Supplemental Fig. S18).

Two observations show that this convergent evolution in cell proportions globally impacts the transcriptome. First, variation in cell proportions correlates with coexpression modules: Modules upregulated in arid species correlate positively with cell populations whose proportions are enriched in these species and reversely for modules downregulated (Fig. 4E; Supplemental Fig. S19). Second, marker genes tend to be enriched in differentially expressed genes significantly up- or downregulated in species living in an arid environment, differently depending on the cell population they mark: For instance, intercalated cells marker genes contain upregulated and proximal tubule marker genes downregulated DE genes (Supplemental Fig. S15).

Taken together, these results imply that both Total and Murinae data sets exhibit a significant level of convergent evolution of cell proportions, again with a stronger signal in the Murinae data set. This convergence in cell proportions very likely shapes the global trend of transcriptome convergence pointed above with PCA and gene modules.

Convergent evolution in coding sequences is mosaic

Studies on convergent evolution in coding sequences have shown that convergence on the exact same amino acid occurs but is rare. Few convergent phenotypes, including echolocation or C4 carbon fixation in plants, exhibit convergent substitutions that can be deemed perfect, such that a transition to an identical amino acid is found on all branches where a transition to the new phenotype occurred (Besnard et al. 2009; Marcovitz et al. 2019). Moreover, the exact same change is not needed to accomplish a similar functional effect on a protein. Here, to characterize cases of convergent evolution in coding sequences, we looked for differences in amino acid usage at a given site between arid and mesic species in the Total_seq and Murinae_seq data sets. We used Pelican, a method that takes into account the phylogeny of the species and that proved to be the best of its kind in a recent benchmark (Duchemin et al. 2023). By simulating positive and negative control sites, we set up a false-positive threshold (0.001) to detect cases with at least two events of convergent evolution (Fig. 5A; Supplemental Fig. S20; Supplemental Methods). We also ensured that our detection is not severely affected by gene length (Supplemental Fig. S21).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Convergent evolution in coding sequences. (A) Distributions of log₁₀(P-value) of Pelican tests for real Murinae_seq data (black; all sites and excluding constant sites), simulated negative controls (blue; 1/1000 quantile is indicated), simulated positive controls (red; number of transitions is indicated; dashed lines show the minimum P-value). (B) Scatterplot displaying log₁₀(P-value) of Pelican tests for all sites common to Murinae_seq (x-axis) and Total_seq (y-axis). The density plot is in black. Blue and dashed red lines as in A for Murinae_seq and as in Supplemental Figure S20A for Total_seq. (C) Top 21 best sites found by Pelican with Murinae_seq data set. (D) Genes ordered according to Pelican's best site-wise P-value in the Murinae_seq data, with 1/1000 threshold in blue. DE genes indicated by triangles. Genes that belong to selected significant GO-term are indicated. Star indicates a significant enrichment (adjusted P-value < 0.1). (E) Same for Total_seq data.

In the Total-seq data set, we selected 2350 gene families with well-aligned single-copy orthologs found in at least 20 mesic and 25 arid species and used in the Total_exp data set to be able to cross the results. We used Pelican to score the genes based on their best site (with the lowest P-value) and retained 729 genes with a P-value below the threshold (Fig. 5, blue line). Six GO terms displayed a significant enrichment, but differentially expressed genes were not enriched, probably owing to their small number (adjusted P-value < 0.1) (Fig. 5D; Supplemental Fig. S20). Our simulated positive controls and the inspection of the best sites show that the lowest observed P-values correspond to about six independent events out of 22 transitions (Fig. 5B; Supplemental Figs. S20, S22). Repeated evolutions were not observed in all independent transitions but were limited to a subset of the families of rodents, and this subset was proper to each gene (but often correlated between sites of the same gene).

In the Murinae_seq data set, using the same method, we characterized 459 genes with convergent evolution in coding sequences out of 2308 gene families with more than seven arid and nine mesic species, corresponding to five transitions. Differentially expressed genes as well as five GO terms showed an enrichment (adjusted P-value < 0.1) (Fig. 5E). Here again, we did not observe sites with convergent evolution at all five transitions but several sites with convergent evolution at four out of five transitions (Fig. 5B,C).

Patterns of convergent evolution in coding sequences therefore appear similar in Murinae_seq and Total_seq data sets, with 20%–30% of genes carrying at least one site with two transitions, but no site undergoing convergent evolution at all transitions. Convergence seems stronger in Murinae as some sites converge at almost all transitions (four of five), which is far from being the case in the Total data set (about six of 22 transitions).

Altogether, the analysis of expression levels, tissue composition, and genomic convergence suggests a greater level of convergent evolution within a single rodent family, Murinae, than in the total data set. One possible reason is that Murinae have recently diverged and thus share a common genomic background. Another reason may be that their transitions all occurred in a more similar time frame, whereas in the total data set some species belong to lineages that have adapted from mesic ancestors tens of millions of years ago, but others adapted very recently.

Convergent evolution in expression levels is stronger for ancient transitions than for recent transitions

To explore the effect of the age of the transitions, we prepared two data sets of size similar to the Murinae, in terms of number of species and number of transitions to the arid habitat (four transitions), but these transitions are either relatively recent (recent_expr, within the same genus and less than 6 MY ago) (Fig. 6A), or more ancient (ancient_expr, at the base of a rodent family and/or more than than 6 MY ago) (Fig. 6E). In both cases, transitions happening in different rodent families were represented so that the “common background” effect was minimized (three different families for the older transitions, two families for the younger transitions; see later how we controlled for this effect).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Detection of convergence in data sets with ancient and recent transitions to arid environments. (A) Phylogenetic relationships of the subset with recent (within genera) transitions. The environments and families of the species are color-coded as in Figure 2, and the dates of transitions are indicated. (B) Phylo-PCA plot and PCA plot using family-corrected expression values. (C) Heatmap representing eigengenes per species and per coexpression module. (D) Differentially expressed genes called by DESeq2 and genes detected by the Eve model (genes passing the FDR threshold are only found in the nonoverlapping genes as shown with the small gray bar). (E–H) Same for the set with ancient (between genera) transitions.

These dates were assessed using time-calibrated trees and an extensive review of the literature (Fig. 2; Supplemental Methods; for chronogram, see Supplemental Fig. S2). It should be noted that some mesic species are present in both data sets, but this should not affect our conclusions, as it is the events happening in the arid branches that interest us. We then applied the same methods as in the previous sections to analyze convergent evolution in expression levels.

In the recent_expr data set, we found that the fourth axis of a PCA analysis of expression levels correlated best, albeit imperfectly, with arid/mesic status and explained 7% of the variation (Fig. 6B). In the ancient_expr data set, the second axis of the PCA clearly separated arid and mesic species and explained 16% of the variation (Fig. 6F; for main axes in the phylogenetic PCAs correlated with phylogeny, see Supplemental Figs. S23, S24). This concords with the number of genes with differential expression found by DESeq2, which was much lower in the recent_expr data set than in the ancient_expr data set (Fig. 6D,H). The number of genes with convergent evolution found by the EVE model was relatively larger in the recent_expr data set, but we show by simulations that this topology tends to favor detection with this method (Supplemental Methods; Supplemental Fig. S12).

Three gene expression modules were significantly correlated with the arid environment in the recent_expr data set, whereas seven were correlated in the ancient_expr data set (Fig. 6C,G). They were enriched for relevant functions in the ancient_expr but not in the recent_expr data set (Supplemental Fig. S25; Supplemental Table S10). This was also coherent with dendrogram analyses of DE genes (Supplemental Fig. S26).

Altogether, whole-transcriptome patterns, gene modules, and differential expressed genes detected by both methods show more convergent evolution in ancient transitions than in recent transitions.

Disentangling the contribution of environmental and temporal factors on expression level and cell composition evolution

In prior sections, we emphasized that convergent evolution in gene expression levels is more pronounced within a single family than throughout the entire data set and that it is also more pronounced for ancient transitions than for more recent transitions. This was noted for data sets that were equivalent in species count and transition count, albeit for a singular instance in each scenario. To generalize this discovery and separate these two effects, we constructed every possible species pair in the data set. Pairs in which both species inherited a change in habitat were eliminated. We quantified the correlation between transcriptomes and between cell compositions for each pair. Of note, we cannot run statistical tests on this analysis because the pairs are not phylogenetically independent.

We observed more variation between kidney transcriptomes of arid species compared with mesic species. This is explained if adaptation to aridity drives kidney diversification and/or drives faster genomic evolution, but mesic species tend to conserve more similar kidneys and/or have slower evolution rates Whatever its origin, this effect will limit the amount of detectable convergent evolution.

Kidney transcriptomes were more similar for pairs that belonged to the same family than for pairs that did not, as expected from phylogenetic relatedness (Fig. 7A). This effect is even greater for pairs of species independently adapted to the arid environment, with species from different families more dissimilar and species from the same family more similar on average than expected from mesic species. We made a similar observation when comparing cell proportions (Fig. 7B). This is explained if convergent evolution in expression levels and cell proportions is stronger when species evolve independently from ancestors that were more alike.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 7.

Similarity of expression levels and cell composition in pairs of species. Transcriptome similarities of expression levels (A) and cell-type composition (B) measured with Pearson correlation coefficients for pairs of species living in similar or different environments, belonging to the same or to different rodent families. (C) Transcriptome similarities of expression levels measured for pairs of species living in arid environments, belonging to the same or to different rodent families, after ancient or recent transitions (all independent pairs from Fig. 3A were taken). (D) Same but correlations measured on all differentially expressed genes detected in the study (N = 1325 genes).

In order to separate this effect from the age of transition, we split the pairs of species adapted to the arid environment depending on the age of the transition (Fig. 3A). As we saw earlier, transcriptome correlations were higher for pairs of species belonging to the same family but also for more ancient transitions (Fig. 7C). We also observed this when we calculated the correlations only on the union of all the differentially expressed genes found in the previous sections (1325 genes) (Fig. 7D). We observed high correlations for the pairs with ancient transitions, whether the species belong to the same rodent family or not, but observed lower correlations for the pairs with recent transitions. This is explained if convergent evolution in expression levels and cell proportions is stronger between ancient adaptations.