Comparable adaptative mechanism however divergent demographic historical past of 4 sympatric desert rodents in Eurasian inland

Species distribution modeling, spatial climate segregation and niche width

To explore the selective regimes of the four species on external environmental factors, we first constructed species distribution modeling (SDM). We obtained a dataset including 22 environmental factors represented by climate, relief, and vegetation variables from 620 localities for DS, 1028 localities for OS, 581 localities for MM and 332 localities for PR, covering most of the species’ distribution ranges (Supplementary Fig. 1). The distribution areas of the four species overlapped widely. The contributions of environmental factors to SDMs showed similarities among the four species. The summer NDVI made important contributions for DS (41.0), OS (44.8), MM (32.5) and PR (8.1), and sand cover contributed significantly to PR (72.7) and DS (16.0) (Fig. 1c). Then, we assessed which set of environmental variables was most closely associated with species distribution via principal component analysis. The bioclimatic space occupied by the four species revealed a large overlap (Fig. 1d), which was consistent with SDM (Supplementary Fig. 1). The distribution of OS was more closely associated with higher-precipitation areas, whereas MM seemed to prefer areas with higher temperatures. Finally, we evaluated the macrohabitat niche breadth of each species. The breadths of environmental space occupation were similar for DS (0.527), MM (0.571), and PR (0.548) and slightly higher for OS (0.622), which suggests that niche selection among the four species is partially overlapping. In total, the four species are mostly similar in the selection of external environmental factors.

High-quality genomic landscapes of the four desert rodents

To investigate the genetic mechanism for desert adaptation of the four sympatric desert rodents, we generated four high-quality de novo genomes (Supplementary Fig. 2). The DS was sequenced using a combined strategy and generated 377.67 Gb of data from Illumina reads, 261.01 Gb from PacBio long reads, 299.51 Gb from 10X Genomics reads, and 389.13 Gb from Hi-C reads (Supplementary Table 1). The final genome size was 2.81 Gb with contig N50 of 31.41 Mb and ~472X mean coverage (Table 1, Supplementary Fig. 3, and Supplementary Tables 2, 3). The contigs for DS were further assembled into pseudochromosomes with lengths on the order of full chromosomes and a scaffold N50 size of 147.24 Mb (Fig. 2a, b, Table 1, and Supplementary Fig. 4). The OS, MM and PR were sequenced using the same hybrid strategy and generated 162.58 Gb, 172.22 Gb, and 214.34 Gb Illumina reads and 183.09 Gb, 161.34 Gb, and 186.45 Gb Oxford Nanopore Technologies long reads, respectively (Supplementary Table 1). The final assembly of OS, MM and PR was 2.83 Gb, 2.43 Gb, and 2.16 Gb with contig N50 of 25.87 Mb, 24.08 Mb, and 42.68 Mb, respectively (Table 1, Supplementary Fig. 4, and Supplementary Tables 2, 3).

Table 1 Genome assembly statistics of the four desert rodents.Fig. 2: High-quality assembly of Dipus sagitta genome and genomic elements of the four sequenced desert rodents.

a Hi-C heat map of Dipus sagitta genome assembly. b CIRCOS plot showing the distribution of GC content, transposable elements (TE), and coding sequences (CDS) in the D. sagitta genome. c Orthologous coding sequences composition inferred for thirteen rodents’ genomes. Mcar Mus caroli, Mmus Mus musculus, Mpah Mus Pahari, Mmer Meriones meridianus, Mung Meriones unguiculatus, Cgri Cricetulus griseus, Prob Phodopus roborovskii, Sgal Spalax galili, Osib Orientallactaga sibirica, Dsag Dipus sagitta, Jjac Jaculus jaculus, Hgla Heterocephalus glaber, Cpor Cavia porcellus. d Proportion of transposable elements (TEs). The barplots show the proportions of different types of TEs in corresponding species on the phylogenetic tree.

Analyses of the four draft genomes showed that 92.9–95.9% of mammalian BUSCOs were complete, and the GC content was 41.38–42.16% (Table 1 and Supplementary Table 3). Whole-genome annotation was performed via three complementary methods: ab initio prediction, homology-based prediction and RNA-seq based prediction. A total of 23,482, 22,859, 22,533, and 22,314 protein-coding genes were annotated for DS, OS, MM, and PR, respectively (Fig. 2c, Supplementary Fig. 5, Supplementary Table 4). Approximately 98.8–99.1% of genes were functionally annotated for the four species (Supplementary Table 4). Transposable elements (TEs) accounted for 31.38–53.02% of genome assemblies, which predominantly consisted of long-terminal repeats (LTRs), long interspersed nuclear elements (LINEs) and other unknown TEs (Fig. 2d). DS and OS displayed significant LTR expansion of 47.39% and 50.88% in four sequenced genomes, while MM showed an unexpectedly high LINE expansion of 28.99% and sharp LTR contraction to 9.38% (Supplementary Table 5).

Phylogenetic relationship and evolutionary history

Using 5,102 single-copy orthologous groups, we constructed a high-confidence phylogenetic tree using the maximum-likelihood algorithm, including time calibrations based on fossil records and previous studies (Figs. 1b, 2c)22. The phylogenetic tree strongly supported nodes uniting the subfamilies Murinae and Gerbillinae, which together represented the family Muridae (Supplementary Fig. 6). This group was sister to a clade containing cricetids. Spalacidae was recovered as the earliest divergent lineage from Muridae and Cricetidae in the superfamily Muroidea. The split of the most recent common ancestor of Dipodoidea and Muroidea dated to ~56.5 Mya (Fig. 1b, Supplementary Fig. 7). In the Miocene epoch (23 Mya–5.3 Mya), accelerated global geotectonic movement aggravated global climate drying and cooling23. Geological disruptions that modified landscapes and offered new habitats favored the early adaptive radiation of extant desert rodents. The ancestors of four sequenced species emerged separately during this period (Supplementary Note 1). Our phylogenetic tree is consistent with previous evolutionary research on rodents22 and supports the independent evolution of desert adaptations in Jerboas, Gerbils and Hamsters.

Expanded and contracted gene families

Comparative genomic analysis revealed 23/32, 4/22, 39/73, and 22/83 gene families exhibiting significant expansion/contraction in the genomes of DS, OS, MM, and PR, respectively (Fig. 1b and Supplementary Fig. 8). Genes belonging to the expanded/contracted families were functionally enriched (Fisher Exact < 0.05) in 127/44, 12/27, 48/81, 55/35 Gene Ontology (GO) terms and 7/15, 7/17, 3/18, and 9/5 KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (p < 0.05) in DS, OS, MM, and PR, respectively (Supplementary Data 1–2). Expanded genes were enriched in pathways including vasopressin-regulated water reabsorption (mmu04962, DS/PR), biosynthesis of unsaturated fatty acids (mmu01040, DS), microtubule-based movement (GO:0007018, DS/PR), kidney development (GO:0001822, DS) and so on, while contracted genes involved in pathways such as drug metabolism-cytochrome P450 (mmu00982, DS), microtubule-based process (GO:0007017, OS/PR), lipid catabolic process (GO:0016042, MM/PR), potassium ion transmembrane transport (GO:0071805, DS) and so on.

Analysis of selection and convergence

We analyzed the positively selected genes (PSGs) and rapidly evolving genes (REGs) in four desert rodents. The phylogenetic tree used for selection analysis has the same topology as the ML tree. When performing branch and branch-site model tests on one desert species, all other species without three other sequenced species were treated as background. Our analyses yielded 1972/1993/1547/1383 PSGs and 694/751/664/583 REGs for DS/OS/MM/PR (Supplementary Data 3–10). Among these genes, 55 PSGs and five REGs were shared by all four desert rodents (Supplementary Fig. 9, Supplementary Table 6).

The convergent/parallel evolution can be evidenced by identifying the same or similar amino acid substitutions of protein-coding sequences, transcriptional regulation of gene expression and similarly enriched metabolic pathways10,24,25. To test whether the three sympatric representative desert rodent taxa have convergent genetic changes to cope with similar challenges, we used the JTT-Fgene model26 and the PCOC method27 to identify genes that are under convergent evolution. PR represented the hamster group, the ancestral clades of MM and Mung represented the gerbil group, and the ancestral clades of DS, OS and Jjac represented the jerboa group (Fig. 1b). We first compared the observed number of convergent amino acid substitutions examined with neutral expectations between the three representative clades in the JTT-Fgene model. A total of 575 genes were identified under the JTT-Fgene model (p-value < 0.05, FDR < 0.05, Passion test, Supplementary Table 7). Three genes (ABTB1, SUSD4, and SYNE3) were identified to be under convergent evolution among the three taxa. The PCOC method considered shifts in amino acid preference instead of convergent substitutions28. We used the PCOC method to verify the convergent gene set obtained by JTT-Fgene model and detected 442 genes undergoing convergent evolution by both methods (posterior probability > 0.95, Supplementary Table 7).

Adaptations in the Arachidonic acid metabolism pathway

The AA metabolism pathway (AAMP) has been shown to be the primary adaptation pathway in almost all desert species12,13,14,15. We first investigated the selective signatures in AAMP across the four species and whether the evolved genes were in line with previous studies in some large desert-adapted animals. We examined 86 genes that belong to AAMP (KEGG pathway: mmu00590) to seek genetic signatures of adaptive evolution. Although the KEGG functional enrichment of AAMP could be observed for PSGs, REGs (except in DS), and convergent evolution genes, they were not significant (p > 0.01) (Supplementary Data 11). We screened the PSGs and REGs and found that the adaptive evolution of this pathway mainly focused on the conversion of lecithin to AA and the production of AA’s cyclo-oxygenase and lipoxygenase products (Fig. 3). DS showed positive selection in seven loci (PLA2G2F, PLA2G4C, PLA2G6, GPX4, ALOX5, PTGES2, HPGDS), OS showed positive selection in five loci (PLA2G1B, PLA2G3, GGT5, ALOX5, PTGES) and rapid evolution in two (GGT5, ALOX5), MM showed positive selection in four loci (PLA2G1B, TBXAS1, HPGDS, PTGIS) and rapid evolution in two (CYP2E1, PTGDS), and PR showed positive selection in two loci (PTGES, PTGIS) and rapid evolution in two (PLA2G4F, CYP2E1) (Fig. 3b).

Fig. 3: The enriched KEGG pathways and diagram of arachidonic acid metabolism pathway evolution in the four desert rodents, Dipus sagitta (DS), Orientallactaga sibirica (OS), Meriones meridianus (MM), and Phodopus roborovskii (PR).

a The enriched KEGG pathways of combined positively selected genes (PSGs) and rapidly evolving genes (REGs) in the four desert rodents. The x axis shows the -Log (P-value) of PSGs and REGs in each pathway. The size of the dot indicates the number of genes. The depth of dot color indicates the enrich ratio. b Alterations in the arachidonic acid metabolism pathway of the four desert rodents. Genes in yellow background are PSGs, genes in violet background are REGs, and red genes are under convergent evolution. The symbol next to the gene indicates the species in which the gene has evolved.

The phospholipase A2 (PLA2Gs) is involved in the hydrolyzation of cell membrane phospholipids, thereby releasing the polyunsaturated fatty acid (AA), which is one of the important sources of AA in organisms29. PLA2Gs showed positively and rapidly involved genetic signatures in all four species. In particular, PLA2G2F has also undergone convergent evolution between hamsters and jerboas (Fig. 3b; Supplementary Table 7). Similar selection signatures of PLA2Gs were also identified in desert adaptation animals such as dromedary (PLA2G2F)2, Tarim red deer (PA2GF and PA2GC)14, and Saudi Arabian indigenous chicken (PLA2G12B)30. After the formation of AA, it could be further metabolized into different eicosanoids, such as hydroxyeicosatetraenoic acid, thromboxanes and prostaglandins, and leukotrienes31. Genes of the cytochrome P450 (CYP) family can help to transform AA into 19(S)-HETE, which is a potent vasodilator of the renal preglomerular vessels that stimulates water reabsorption32. Although 19(S)-HETE is thought to play an important role in the survival of desert animals such as camels12, sheep15, and Tarim red deer14, only CYP2E1 in the cytochrome P450 pathway showed rapidly evolving selection in MM and PR, and the cytochrome P450 gene family showed contraction in all four desert rodents.

Thromboxanes and prostaglandins can influence endothelial and vascular smooth muscle cell function and then indirectly regulate blood pressure, water retention and reabsorption, and immediate hypersensitivity reactions33. Therefore, it is thought that selective signals related to vascular smooth muscle regulation or bronchoregulation may be associated with desert adaptation, such as asthma caused by airborne dust2,14. Cyclooxygenase metabolizes AA to prostaglandin H2 (PGH2), which is further converted into biologically active prostaglandins33. We detected selection signatures for some genes (TBXAS1, PTGES, PTGES2, HPGDS and PTGDS) involved in the bioactivation of PGH2. For example, PTGES/PTGES2 convert PGH2 to prostaglandin E2 (PGE2), which plays a key role in thermoregulation, inflammation response, and bronchodilatation34,35, as well as increases renal blood flow and provokes dieresis, natriuresis, and kaliuresis36. Although prostaglandins have a similar role in vascular smooth muscle tone regulation as angiotensin and 19(S)-HETE, adaptative signatures of genes in the cyclooxygenase pathway has rarely been mentioned in previous studies of desert animals12,13,15. In addition to biologically active prostaglandins, leukotrienes, which are lipid mediators of inflammation and chemotaxis, are also potent bronchoconstrictors that play an important role in immediate hypersensitivity reactions37. We detected selection signatures in lipoxygenase pathways (ALOX5, GGT5, GPX4) that produce leukotrienes, which are different from selective signatures such as TRAF2 in Tarim red deer14 or functional enrichment such as lung development (GO:0030324) in camels2.

The entire AAMP undergoing convergent evolution in desert-adapted mammals and birds has become a consensus3,12,13,15. However, in the selection of downstream AAMP, the specific genes involved differed among species and in the mechanisms by which they putatively affected adaptive phenotypes. For example, desert rodents seem to show more selective pressure on the prostaglandins and leukotrienes pathway (Fig. 3b) rather than selecting on cytochrome P450 in desert-adapted large herbivores or chickens12,14,15,30. Therefore, we hypothesized that regulation of biologically active prostaglandins might be another adaptative way to desert environments in desert organisms. There are also differences in selection preferences among the four sequenced rodents, such as the lipoxygenase downstream pathways showing selection signatures and functional enrichment “inflammatory response” (GO:0006954) only in jerboas (Supplementary Fig. 10, Supplementary Data 12). This may suggest that similar adaptive phenotypes in different desert organisms may be achieved through a combination of convergent and idiosyncratic adaptations.

Adaptations in thermostasis and energy homeostasis

Most desert rodents are nocturnally activity and hide in burrows during the day to avoid overheating. Although the daytime temperatures were extremely high in the four desert rodent habitats, the Eurasian inland where they were distributed was a cold desert with extreme diurnal ambient temperature variation. Adaptive thermogenesis is essential for them to ensure normal cellular and physiological function under conditions of environmental stresses38. Among the four species, we found similarly significant enrichment of functional pathways (p < 0.05) associated with thermic and energy processes, including 65 PSGs/REGs in the thermogenesis pathway (mmu04714), 34 PSGs/REGs in the oxidative phosphorylation pathway (mmu00190), 12 PSGs/REGs in the citrate cycle (TCA cycle) pathway (mmu00020) (Figs. 3a, 4a, Supplementary Data 11–13). In the presence of changes in ambient temperature, alterations in food intake, and/or sympathetic nervous stimulation, brown adipose tissue (BAT) activity increases resulting in heat generation (Fig. 4a). Upon activation by long-chain fatty acids, UCP1 transports protons through the inner mitochondrial membrane and uncouples mitochondrial substrate oxidation from ATP production thereby releasing the energy of fatty acid oxidation as heat39. Therefore, the oxidative phosphorylation pathway produces ATP to maintain physiological metabolism while also producing protons to maintain the proton gradient on both sides of the mitochondrial inner membrane (Fig. 4b). It is noteworthy that DS and OS showed more selective genes than the other two desert rodents (Supplementary Data 11), which may be because jerboas undergo hibernation during winter. As hibernating animals require intense thermogenesis to arouse from the deep torpor when they are hypothermic (as low as 4–6 °C), both shivering and nonshivering mechanisms need to be invoked rapidly (within a few hours) to raise the body temperature during arousal40. In addition to the four species, sweep signatures associated with cold-induced thermogenesis was also detected in other desert rodents10. In contrast, none of the small desert rodents10,13 showed enrichment on melanogenesis pathway compared with desert-adapted herbivores or chickens14,30,41.

Fig. 4: Schematic mechanisms of thermostasis and energy homeostasis pathways in the four desert rodents.

a Alterations in the thermogenesis pathway, insulin related pathway, and non-alcoholic fatty liver disease pathway of the four desert rodents. Genes in yellow background are PSGs, genes in violet background are REGs, and red genes are under convergent evolution. The symbol next to the gene indicates the species in which the gene has evolved. DS, Dipus sagitta, OS, Orientallactaga sibirica, MM, Meriones meridianus, PR, Phodopus roborovskii. b Evolved genes in the oxidative phosphorylation of the four desert rodents. Genes in red are PSGs, genes in green are REGs, and violet genes are under both positive selection and rapidly evolving. Different colored bars represent different species. The five sections of the same colored bar correspond to the five complexes. The underlying graph of the mitochondrial electron transport chain is from the KEGG website (https://www.genome.jp/pathway/mmu00190).

Moreover, the uptake and homeostasis of energy is important for desert animals living in food-deprived deserts. The scarcity of water and resources drives metabolism to rely on endogenous nutrients (carbohydrates, lipids, and proteins). The genome-wide features of the four species indicate that both “proteins-for-water”42 and “lipids-for-torpor”43 metabolic regulation contribute to desert adaptation. The transition to protein catabolism during early resource restriction may be an important water source for desert species10. Lipid acquisition and storage are important to extreme environmental-induced torpor because they enable desert species to maintain a low metabolic state for a long time10,43. For example, ACOT1, which catalyzes the hydrolysis of long-chain fatty acyl-CoAs into free fatty acids and coenzyme A44, was identified by KEGG pathways in fatty acid elongation (mmu00062, DS) and biosynthesis of unsaturated fatty acids (mmu01040, DS) (Supplementary Data 1), while functional enrichments of contracted gene families were identified in lipid-related GO categories (lipid catabolic process, GO:0016042, MM/PR; medium-chain fatty acid metabolic process, GO:0051791, MM) and protein-related pathways (cysteine and methionine metabolism, mmu00270, DS/MM; tryptophan metabolism, mmu00380, DS/MM; proteolysis, GO:0006508, PR) (Supplementary Data 2). The positively selected and rapidly evolving functional enrichments associated with endogenous nutrient metabolism of the four desert rodents included insulin signaling pathway (INSP) (mmu04910, DS/OS), non-alcoholic fatty liver disease (NAFLD) (mmu04932, all four species), lipid metabolism (mmu00561, DS/OS/MM; mmu00564, DS/OS/MM; mmu04979, OS/MM/PR; mmu00600, OS/PR; GO:0006629, all four species), lipid and cholesterol homeostasis (GO:0042632, OS/MM/PR; GO:0055088, OS), amino sugar and nucleotide sugar metabolism (mmu00520, OS/PR), proteolysis (GO:0006508, all four species), and so on. The convergent feature genes also demonstrated an enrichment of KEGG pathways and GO terms associated with metabolism or homeostasis of proteins or amino acids (mmu00330, mmu00270, GO:0006547, GO:0006548), lipids (mmu00600, GO:0090336), and carbohydrates (mmu00520, mmu00051, GO:0042593, GO:0045721) (Supplementary Data 13).

Many genes with selection signals or convergent features were functionally enriched in the INSP and NAFLD pathways in the four desert rodents, although the evolved genes in each species were not fully identical (Figs. 3a, 4a). The insulin-related pathway has been considered genetic bases for previously identified adaptative phenotypes, including energy storage, low energy expenditure, adaptive tolerance to starvation and dehydration2,3,10,45. Homeostasis of carbohydrates and lipids via INSP is achieved by regulating downstream pathways such as glucose uptake, glycolysis, glycogenesis, antilipolysis, lipogenesis and fatty acid synthesis. For example, G6PC (OS), hydrolyzing glucose-6-phosphate to glucose in the endoplasmic reticulum (ER), is the key enzyme in the homeostatic regulation of blood glucose levels46; LIPE (OS), with broad substrate specificity, catalyzes the hydrolysis of several lipids47; FASN (PR) catalyzes the de novo biosynthesis of long-chain saturated fatty acids in the presence of NADPH48. NAFLD, for which selection signals have been observed only in desert rodents, also plays a role in lipid homeostasis (Supplementary Data 11)10. The induction of insulin resistance activates key enzymes of lipogenesis and increases the synthesis of free fatty acids in liver leading to excess lipid accumulation49. For instance, NRI3H (OS) induces phospholipid remodeling in hepatocytes50; MLXIP and MLXIPL (DS and OS) play roles in the transcriptional regulation of glycolytic target and glucose-responsive genes (Supplementary Data 3–10)51.

Adaptations in stress response

Although desert rodents are mostly nocturnal, they are still not immune to the stress of heat and aridity from the regional environment. Thermal stress suppresses the transcription and translation machinery, increases DNA breaks and protein oxidation, alters cell morphology and cell adhesion, causes cell cycle arrest, and eventually leads to apoptosis and cell death52. Hyperosmotic stress caused by water shortage and water loss indirectly accelerates the rapid accumulation of cellular damage and exacerbates heat killing53. Repair of the genetic expression process has always been an important mechanism for animals to adapt to extreme environments2,11,13,54. Our results indicated that many genes with the strongest selection signatures were enriched in the GO terms and KEGG pathways related to various signaling pathways and different stages of genetic expression processes (Supplementary Data 11–13).

In the DNA replication and repair stage, many PSGs/REGs were enriched in multiple GO categories related to different types of DNA repair, as well as in KEGG pathways such as base excision repair (mmu03410) and homologous recombination (mmu03440) (Fig. 3a, Supplementary Fig. 10, Supplementary Data 14). Among these PSGs/REGs, four genes (LIG3, MUTYH, POLL, POLM) were also under convergent evolution and separately involved in base excision repair (LIG3, MUTYH), DNA double-strand break repair (POLM) or both (POLL) (Supplementary Data 14)55,56. Although DNA repair-related pathways are enriched in all four species and partly overlap among different species, genes under selection pressure in these pathways are not the same in the four species. The number of PSGs/REGs was significantly higher in jerboas, which also shared more PSGs/REGs.

In the transcription stage, a large number of pathways related to transcription and splicing were enriched in gene family evolution, selective evolution and convergent evolution, such as regulation of transcription, DNA-templated (GO:0006355 in four species), regulation of transcription by RNA polymerase II (GO:0006357 in DS, OS and MM), spliceosome (mmu03040 in four species), mRNA processing (GO:0006397 in DS and OS), RNA splicing (GO:0008380 in DS and OS), mRNA transport (GO:0051028 in OS and PR) (Supplementary Data 11–13). Similarly, all analyses indicated that traits associated with the synthesis and degradation of proteins have evolved under the influence of natural selection. The extended or contracted gene families were enriched in various GO categories and KEGG pathways, such as proteasome (mmu03050 in PR), phagosome (mmu04145 in DS, OS, and PR), apoptosis (mmu04210 in OS), and protein processing in endoplasmic reticulum (mmu04141 in MM) (Supplementary Data 1–2). The PSGs/REGs and convergent genes also demonstrated an enrichment of terms and pathways associated with this process, such as protein processing in endoplasmic reticulum (mmu04141 in DS, OS and MM), lysosome (mmu04142 in OS, MM and PR), ubiquitin mediated proteolysis (mmu04120 in DS and OS), ribosome (mmu03010 in OS and MM), proteolysis (GO:0006508 in four species), protein ubiquitination (GO:0016567 in DS, OS and MM) (Supplementary Data 11–13).

In the regulation of genetic expression processes, various signaling pathways are commonly linked by some of the same node genes and cascade reactions, which jointly regulate individuals’ responses to extreme environmental stress. For example, CREB is a phosphorylation-dependent transcription factor that stimulates transcription upon binding to the DNA cAMP response element (CRE)57. It is enriched in at least 20 pathways in each species, most of which are signaling pathways such as the insulin signaling pathway, PI3K-Akt signaling pathway, and AMPK signaling pathway (Supplementary Data 11–13). These signaling pathways did not directly affect some adaptive phenotypes, but indirectly affected several key traits of survival in hot arid environments by regulating downstream genes10,14,15,41, including energy metabolism, vascular smooth muscle tone regulation, repair of genetic expression processes, thermotolerance, or thermogenesis.

Similar results but divergent adaptative routes

Although all kinds of desert animals are adapted to desert environments, it seems that desert rodents have chosen different evolutionary routes in phenotypes and genotypes2,13,14,15. Taking heat regulation as an example, some large desert animals have adaptive evolution in storing heat without increasing body temperature3, the melanogenesis pathway14,28,41, photoreception and visual protection2,14 and other characteristics, while desert rodents have adaptative selection in thermogenesis (Supplementary Data 11–13). Regarding the water retention and utilization, previous studies and our results have found that more traits associated with the synthesis and degradation of proteins in desert rodents have evolved under the influence of natural selection10,13, while some large desert animals have relatively more evolved genes and functionally enriched pathways related to water balance and reabsorption, as well as osmoregulation2,14,15,41.

The three desert rodent taxa from the Eurasian inland show a potential framework of convergent or parallel evolution as well as similar selection on external environmental factors. However, our results suggest that the same adaptative molecular pathways could be achieved through several genomic routes instead of identical replacements of single amino acids within the encoded product of a protein-coding gene58. In macroecology, species with overlapping niches may come into competition with each other when they exploit the same limiting resources59, while niche segregation facilitates coexistence60. Here, we propose a hypothesis: if metabolic pathways are treated as kinds of resources (genomic niches), whether ecologically similar sympatric species are subjected to different evolutionary selection in the same metabolic pathway can be interpreted as genomic niche separation promoting their coexistence.

Demographic history

Quaternary climatic cycles have a profound influence on the formation and current distribution of species61. During the Middle Pleistocene (787–130 ka), the climate of the Eurasian inland gradually tended to be dry-cool, even though there were climate cycles. In the late Pleistocene (130–10 ka), the climate reached maximum dry cooling at 21 ka and then began to change rapidly to warm and humid characteristics (Fig. 5). To determine the relationships between the population demographic size and environmental dynamics over time, we estimated the effective population size (Ne) and climate temporal shifts for the four species over approximately the last one million years using the pairwise sequential Markovian coalescent (PSMC) method62. Our PSMC results show that the DS’s expansion began at ~1.0 Mya and presents an overall declining trend (Fig. 5). The OS’s effective population size showed drastic fluctuation with two expansions and two declines. The population expansion of OS significantly coincided with the expanded dry-cooling climate (Fig. 5). Moreover, during the period from 21 Ka to 12 Ka, there were significant decreases in cold temperatures (<−5 °C) and arid regions (< 400 mm) and sharp population reductions. This revealed that OS is more sensitive to warmer temperatures and wetter environments than the other species. MM almost reached the first expansion peak at ~1.1 Mya, followed by another two expansions. The PR’s expansion apparently began at ~500 Ka, followed by another two expansions. Although the macrohabitat choice and genomic features of the four species showed similar evolved selection in desert adaptation, they exhibited different dynamic patterns and relationships with temporal climate changes. We infer that the diversified genomic adaptative features lead to the different demographic histories. In total, all four species experienced repeated fluctuations in effective population size during the Quaternary climatic cycles, which may be related to the historically repeated expansions and contractions of desert–grassland–forest systems in the Eurasian inland. DS, OS, and PR appear to have declined in population size since the Last Glacial Maximum. In contrast, MM had a much larger overall effective population size than other species, which may be due to a broader gene pool of variation for selection10. This could be beneficial for MM to adapt to more diverse habitat types in the desert environment, which is consistent with our observations in the field. The population sizes of the two jerboas appeared to be more sensitive and vulnerable to fluctuations in environmental factors than MM.

Fig. 5: Demographic history of the four desert rodents during historical climate changes.

a Effective population sizes of the four desert rodents reconstructed using the PSMC model (g, generation time, years; μ, neutral mutation rate per generation). The period of the Penultimate Glaciation (PG, 330–130 ka), and the Last Glaciation (LG, 75–10 ka) are shaded in light blue, and the period of the Last Interglaciation (LIG, 130–800 ka) is shaded in light orange. b The relative temporal changes in annual precipitation and annual mean temperature from 787 Ka to present within a defined region (Lat: 33 N, 55 N; Long: 70E, 130E) including the distributions of the four desert species. The annual mean temperature and annual precipitation layers during the Quaternary were downloaded from PaleoClim database.