We downloaded the reference genome of Nannospalax galili (GCF_000622305.1) and pair-end (PE) reads from NCBI at https://www.ncbi.nlm.nih.gov/. PE reads were filtered by removing low-quality reads and adapters with fastp v0.20.1 (Chen S. et al., 2018). The clean reads were mapped to the reference genome with BWA-MEM v0.7.17-r118 15 (Li and Durbin, 2009), sorted and indexed with SAMtools v1.11 (Li H. et al., 2009). The duplicates were marked and variants were detected using GATK v4.0 16 (McKenna et al., 2010). SNPs were extracted and performed hard filtering by GATK with default parameters. Then SNPs were further filtered using VCFtools v0.1.16 (Danecek et al., 2011) with parameters: --maf 0.05. --minDP 5. --maxDP 50. --minGQ 20. --hwe 0.01 and --max-missing 0.9. A genetic distance matrix was generated by PLINK (Purcell et al., 2007). We next used FastME 2.0 (Lefort et al., 2015) to construct a phylogenetic tree with default parameters.

In total, 12 gut samples were collected for chalk (6 samples) and its abutting derivative basalt (6 samples) from “Evolution Plateau” (EP): in the eastern Upper Galilee, Israel, in January 2016 and were stored in liquid nitrogen immediately for DNA extraction. DNA was isolated using QIAmp DNA stool mini kit (Qiagen, United States). Qualified DNA was sonicated to about 350 bp fragments, after purification, end-repair, A-tailing, and adaptor ligation were carried out and followed by PCR amplification. After quality control, libraries of each sample were sequenced with 150 bp paired-end reads. We generated ~50 million 150 bp PE raw reads from the shotgun metagenomic sequencing for each sample. Paired-end raw sequencing reads were first trimmed by SOAPnuke (Chen Y. et al., 2018) to remove adapters and low-quality bases, and the host reads (S. galili) were removed by SOAPaligner (Li R. et al., 2009). We also collected 12 soil samples for shotgun metagenomic sequencing at the same location (Mukherjee et al., 2022), and the data of soil were processed as described below.

First, clean reads were uploaded into the MG-RAST (Meyer et al., 2019) server, and bacteria were found to be the most dominant kingdom (over 99.5%). Kraken2 v2.1.2 (Wood et al., 2019) was used for taxonomic profiling, and the relative abundance was estimated by Bracken v2.5 (Lu et al., 2017). Meanwhile, alpha and beta diversity were calculated using R packages “vegan”.¹ Species with significant abundance differences between basalt and chalk were identified using LEfSe (linear discriminant analysis effect size; Segata et al., 2011) on the Galaxy website (v1.0).² We generated a Bray-Curtis distance matrix based on the microbiome composition of the gut using R packages “vegan,” and the microbiota dendrogram was constructed based on the distance matrix by R packages “ape.” Mantel-test was performed to compare the relationship between the host genetic distance matrix and Bray-Curtis distance matrix based on the gut microbiome. Procrustes analysis and species accumulation curves were performed on the Tutools platform.³ Statistical analysis was conducted with R.

Contigs were assembled by MEGAHIT v1.1.3 (Li D. et al., 2015) with default parameters. After assembly, assembled contigs were binned, refined, and reassembled by metaWRAP v1.3.2 (Uritskiy et al., 2018). The completeness and contamination of each bin were evaluated by CheckM v1.0.12 (Parks et al., 2015). All bins were aggregated and then dereplicated using dRep v3.2.2 (Olm et al., 2017; parameters: -comp 50 -con 10 -pa 0.90 -sa 0.95 -nc 0.30 -cm larger). The taxonomy of bins was assigned using GTDB-tk v1.7.0 (Chaumeil et al., 2020) and these results were visualized in iTOL v6 (Letunic and Bork, 2007) as a phylogenetic tree. The relative abundance of bins was estimated by CoverM v0.6.1.⁴ Bray-Curtis distance matrix and PCoA were generated by R packages “vegan,” and the heatmap was plotted by https://www.bioinformatics.com.cn.

The coding sequences (CDS, > 100 bp) were predicted by Prodigal v2.6.3 (Hyatt et al., 2010), and we used CD-HIT v4.8.1 (Fu et al., 2012; parameters: -c 0.95 -d 0 -aL 0.9 -uL 0.05 -aS 0.9) to remove redundancy. The non-redundant gene catalog was searched against the Kofam database by KofamKOALA v1.3.0 (Aramaki et al., 2020). Because KEGG Orthology and KEGG pathway have a many-to-many relationship, we used ReporterScore (Bäckhed et al., 2015) and biostack-suits⁵ to analyze all KOs and used the overall trend to reflect the change in the pathway. Carbohydrate-active enzymes (CAZyomes) were identified by dbCAN2 (Zhang et al., 2018). Clusters of Orthologous Groups of proteins (COG) and Gene Orthology (GO) were determined by eggNOG-mapper v2.1.6 (Huerta-Cepas et al., 2017). Meanwhile, alpha diversity and beta diversity were also calculated by R packages “vegan.” Significantly different GO terms and CAZyomes between basalt and chalk were identified by STAMP v2.1.3 (Parks et al., 2014). Statistical analysis was conducted with R and significant outliers were removed. The relative gene abundance was calculated as follows (Qin et al., 2012; Liu et al., 2021).

Step 1: Clean reads were aligned to the non-redundant gene catalog using BWA (Vasimuddin et al., 2019), and genes with a number of mapped reads less than two were removed; the copy number of each gene was: b_i = x_i/L_i;

Step 2: Calculation of the relative abundance of gene i: a_i = b_i/∑b_i;

a_i: the relative abundance of gene i.

b_i: the copy number of gene i from sample N.

L_i: the length of gene i.

x_i: the number of mapped reads.

The calculations mentioned above were performed with custom Python scripts.

First, we used kraken2 v2.1.2 to identify the most abundant species and significantly different species were identified by LEfSe; second, we downloaded their genomes from NCBI as references, but the reference genome of Akkermansia muciniphila we used was generated in the present study. SNPs were detected the same as described above with the following parameters: “--maf 0.05. --hwe 0.01 and --max-missing 0.9.” Principal component analysis (PCA) was accomplished by PLINK v1.90b6.21 (Purcell et al., 2007). Neighbor-joining (NJ) tree was conducted by TreeBeST v1.9.2 (Vilella et al., 2009). Demographic history was analyzed by SMC++ (Terhorst et al., 2017). F_ST and nucleotide diversity (π) were calculated by VCFtools in 1 kb sliding windows.

Because the abundance of Flavonifractor plautii is also high in the host dwelling soil, so we compared both the gut and soil microbiome to investigate the differentiation of this species. Furthermore, admixture v1.3.0 (Alexander and Lange, 2011) was used to perform structural analysis with the number of clusters (K) ranging from 2 to 4, respectively. A genetic distance matrix was generated by PLINK, and visualization of the genetic network was finished by R packages (“netview,” “network,” “graph,” “sna,” “visNetwork,” “three,” and “networkD3”) with k = 7 (Steinig et al., 2016).

We collected eight liver samples of the host Spalax, including four individuals from basalt and four from chalk, for proteome comparison. Tissues were frozen in liquid nitrogen immediately and stored at −80°C. The samples were ground first, then add 0.4 ml protein cracking liquid, containing 100 mM Tris–HCl (pH 8.0), 10 mM DTT, 8 M urea, and 1X protease-inhibitor. The mixture was incubated on ice for 30 min and then centrifuged at 10,000 g and 4°C for 30 min. The supernatant was collected, and protein concentration was determined with a BCA protein assay kit. Protein quality was identified with SDS-PAGE.

Protein samples were diluted with NH₄HCO₃ (200 mM), and incubated with 10 mM DDT for 1 h at 56°C, after cooling, add 55 mM iodoacetamide (IAA) for 40 min in the dark to alkylate samples, then digested with 5 μg trypsin for 14-16 h at 37°C. Digested peptides were concentrated to about 1 mg/ml, then separated with chromatography using Eksigent 425 (AB SCIEX). Separated peptides were performed mass spectrometric analysis with Q-Exactive (Thermo Scientific).

Raw files were extracted by MaxQuant (Cox and Mann, 2008). MS data were searched against Uniprot- Heterocephalus glaber database with parameters: max missed cleavages were 2; Carbamidomethylation (C) were set as fixed modifications, and oxidization (M) were set as variable modifications; Peptide Mass Tolerance was ±15 ppm; fragment mass tolerance was 20 mmu; peptide length was >4. The cutoff of the global false discovery rate (FDR) for peptide and protein identification was set to 0.01. Then screening proteins based on the threshold with a minimum fold change of 1.5 (p < 0.05). Statistical analysis was conducted by Student t -test with Microsoft Excel. For proteins with missing values, data for that protein were kept if its lowest value in one group was higher than its highest value in another group. Finally, Gene ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genome (KEGG) pathway enrichment analysis were performed to characterize the properties of the proteins identified using STRING (Szklarczyk et al., 2015), these results were visualized by https://www.bioinformatics.com.cn.

A total of ~95 GB of raw reads (Supplementary Table 1) from the shotgun metagenomic sequencing was generated, and ~ 80 GB of clean reads (Supplementary Table 1) were retained for six basalt and six chalk gut samples after trimming the low-quality bases. We obtained an average of 570,905 contigs (>200 bp) for each metagenomic assembly for downstream analysis (Supplementary Table 2). In total, 339 bins (>50% completeness, <10% contamination) were obtained (Supplementary Figures 1, 2), with 139 high-quality bins as defined by Bowers et al. (2017; >90 complete, <5 contamination). All bins can be annotated at the family level (Supplementary Table 3), with the majority of bins (614 bins, 90.4%) also annotated at the genus level (Supplementary Table 3), indicating potential new species in the gut. One bin (Chalk3.37) can be annotated to the species level (Akkermansia muciniphia). PCoA based on the relative abundance of all bins showed a clear separation between the chalk and basalt populations (Figure 1E; ANNOVA, p = 0.001998). The heat map also showed similar results (Supplementary Figure 2), indicating that chalk and basalt differ at the strain level.

We checked the sample size saturation and found that as the sample number increased to nine the species accumulation curve (Supplementary Figure 3) began to plateau, indicating that our 12 samples were sufficient to cover most species. We characterized the overall microbiome composition variation by PCoA based on Bray-Curtis distance and found there was no significant difference between the chalk and basalt (Supplementary Figure 4; ANNOVA, p = 0.314). Additionally, chalk and basalt did not differ significantly in alpha diversity (Supplementary Figure 5). We found Akkermansia muciniphila, Clostridium innocuum, Stenotrophomonas sp. CW117, Enterobacter hormaechei, Clostridium autoethanogenum were overrepresented in chalk, and Christensenella minuta was enriched in basalt (Figure 1F; LDA score (log10) > 1.5).

The most dominant kingdom in both populations gut microbial communities were bacteria (over 99.5%), followed by archaea and eukaryota. For the bacterial community, 57 phyla, 99 classes, 214 orders, 469 families, 1721 genera, and 6,401 species were identified from 12 gut samples. We found 2 phyla, 8 classes, 16 orders, 24 families, 109 genera, and 313 species were significantly different between basalt and chalk (Welch’s t-test, p < 0.05) populations. In both populations, the most abundant phyla were Firmicutes that accounted for 59.6% of abundance, followed by Proteobacteria (17.4%), Bacteroidetes (10.7%), ActinobacteriaI (7.61%) and others (4.57%, Figure 1D). Paenibacillus (6.21%), Lachnoclostridium (5.83%), Flavonifractor (4.24%), Oscillibacter (3.99%) and Blautia (3.79%) were the most abundant genera. It is worth noting that Chalk-3 contained considerable Akkermansia muciniphila which belonged to the Verrucomicrobia phylum, Akkermansia genus (6.09%, Figure 1D). Moreover, we performed Welch’s t-test and found significant difference in relative phylum abundance of Proteobacteria (Figure 2D; p = 6.32e-3) and Candidatus Micrarchaeota (p = 0.013). At the genus level, Christensenella, Acidovorax, Bordetella, Bibersteinia and Faecalitalea were enriched in basalt; While in chalk, Providencia, Acidithiobacillus, Gimesia, Gemella and Tardiphage were enriched. The ratio of Firmicutes to Bacteroidetes (F/B ratio) was higher in the basalt compared to chalk (Supplementary Figure 6). We also compared the microbial diversity between the soil and gut and found significant differences in community composition (Figure 2A) and significantly higher Shannon diversity in the soil microbiome (Supplementary Figure 7). Meanwhile, Procruster analysis (Figure 2B; M² = 0.0761, p = 0.55) and mantel test (r = − 0.01518, p = 0.5317) showed there was no correlation between the gut and soil microbiomes.

We plotted the cumulative curves for KEGG, GO and CaZy and all curves began to plateau, indicating that our 12 samples were sufficient to reveal the function of both communities (Supplementary Figure 8). Principal coordinate analysis (PCoA) of the non-redundant gene catalog based on Bray-Curtis distance showed that the functional composition of basalt and chalk was different (Figure 3A; ANNOVA, p = 0.00995). The Shannon index of basalt was higher than that of chalk (Figure 3B) but was not significantly different. The richness index, however, was significantly higher in chalk (Figure 3C; p = 0.0039). For GO terms, the non-redundant genes of basalt were enriched in growth, oxidoreductase complex, transporter activity, cation transmembrane, etc. (Figure 3D), and enriched in mannan-binding, catalytic activity, and negative regulation of cardiac muscle cell differentiation in chalk. The significantly enriched KEGG pathways (Figure 3E) in basalt included: Oxidative phosphorylation, Purine metabolism, Glyoxylate and dicarboxylate metabolism, Methane metabolism, etc. CAZymes (carbohydrate-active enzymes) compositions were significantly different (ANNOVA, p = 0.014; Figure 2D) between basalt and chalk. The diversity (Supplementary Figure 9) and abundance (Figure 3F) of CAZymes were higher in chalk, which contains more (Figure 3G) Carbohydrate-binding module (CBM), carbohydrate esterase (CE), glycoside hydrolase (GH) and glycosyl transferase (GT). GH43, GH28, GH94, GH97 and PL12 were enriched in basalt, while GT4, GT10, GH92, PL9, GT26 and GT11 were enriched in chalk. For COG (Cluster of Orthologous Groups database) annotations (Supplementary Figure 10), the metabolism category was dominant in both chalk and basalt; cellular processes and signaling were the second most dominant. Furthermore, Replication, recombination and repair (COG L) was the most abundant COG type, followed by Carbohydrate transport and metabolism (COG G) and Amino acid transport and metabolism (COG E; Supplementary Figure 10). Recombination and repair (COG L) were enriched in basalt (Supplementary Figure 11C). Transcription (COG K), cell wall/membrane/envelope biogenesis (COG M) and energy production and conversion (COG C) were enriched in chalk (Supplementary Figures 11A,B,D).

The species identified by LEfSe were similar between the chalk and the basalt (Supplementary Figures 26–31) populations; however, many species showed separate genetic clusters between basalt and chalk measured by PCA (Figure 4A; Supplementary Figures 12A–31A), phylogenetic tree (Supplementary Figures 12A, 13B–31B) and genetic networking (Figure 4C) based on SNPs. Demographic analysis revealed that microbiome differed between chalk and basalt, but the differentiation was not high (Supplementary Figures 12B, 13C–31C). Interestingly, the differentiation time of Akkermansia muciniphila between the chalk and the basalt populations is estimated to be approximately 6,000 to 7,000 years. These changes have occurred more recent than the genetic differentiation of the two host species which is estimated to be approximately 228,000 years 33 by SMC++ analysis (Supplementary Figure 26C). Other microbial species also have a shorter differentiation time (Supplementary Figures 26–31). F_ST and nucleotide diversity (π) were also calculated for these species based on SNPs (Supplementary Table 4).

Furthermore, we explored the differentiation of the same bacterial species of Flavonifractor plautii between the chalk and basalt gut microbial composition and the two types of soils using PCA (Figures 4A,B) and genetic networking (Figure 4C). We revealed that it was separated into two clear-cut clusters in the soil but not in the gut (Figures 4C,D). F_ST was also calculated for both the gut and soil, and it was significantly higher in the soil than in the gut microbiome (Figure 4D). Compared to the genetic differentiation of the hosts between the chalk and basalt, the differentiation of the gut microbiome is later and smaller, indicating that the gut microbiome was not a major factor influencing host divergence. When the individuals were separated into two groups (K = 2), the basalt population is completely separated from the chalk population (Supplementary Figures 32, 33).

We performed a mantel test and found there was no correlation (r = 8.078e-3, p = 0.4528) between the host genetic matrix and the Bray-Curtis distance matrix based on the gut microbiome. The microbiome dendrogram did not reflect the host phylogeny (Figure 5A). PCoA based on protein abundance showed little difference in protein composition between the hosts from basalt and chalk (Figure 5B). However, we did identify 127 significant differential proteins (Figure 5C; p < 0.05 Student t-test) that were significantly enriched in metabolic pathways such as endocytosis, biosynthesis of amino acids, endocrine and other factor − regulated calcium reabsorption, arginine and proline metabolism, vasopressin−regulated water reabsorption, and pyruvate metabolism (Supplementary Figure 34). GO enrichment of Biological Process and Molecular Function also showed similar results (Supplementary Figure 34): differential proteins were mainly involved in the synthesis of proteins (e.g., Organonitrogen compound metabolic process, cellular metabolic process, translation, Peptide biosynthetic process), they may be secreting proteins (e.g., Regulation of multivesicular body size and Retrograde transport, endosome to the plasma membrane), and were followed by the synthesis and decomposition of various biological substances (e.g., Amide biosynthetic process, Negative regulation of cholesterol biosynthetic process and Fumarate metabolic process). Enrichment of GO Cellular Component (Supplementary Figure 35) also indicated that these proteins may be secretory proteins or have transport functions (e.g., Cytoplasmic vesicle, Late endosome, and Clathrin-coated vesicle membrane).

Although hundreds of samples are frequently present in many studies, we had only a small number from each population but were still able to demonstrate that richness was not increasing with sample size (Supplementary Figures 3, 8), suggesting the sample size was not a restriction for the analysis. Firmicutes microbes had the largest relative abundance in both chalk and basalt gut microbiome (Figure 1D), which supports evidence found in other studies (Qin et al., 2010; Xiao et al., 2015; Huang et al., 2018). Proteobacteria was significantly higher in the basalt digestion tract (Figure 2D). This particular phylum of bacteria is reported to be positively correlated with the fat intake of the host diet, and is significantly richer in populations with a high-fat diet than that of the malnutritional population (Méndez-Salazar et al., 2018). In this study, the basalt mole rat population was mainly feeding on geophytes (Jiao et al., 2021; Joyce et al., 2022), which have higher fat than that in Eryngium sp. roots from chalk; this is congruent with the functional enrichment pathway of fatty acid degradation and Glycerolipid metabolism (Figure 3E) in basalt population. Bacteroidota was higher in chalk than in basalt (Supplementary Figure 6A), which may be caused by the higher fat content in the basalt diet suppressing this phylum (Jeong et al., 2019). If the host is feeding on more proteins and fat, the ratio of Firmicutes/Bacteroidota would be higher (De Filippo et al., 2010). In the present study, the ratio was higher in basalt than in chalk (Supplementary Figure 6B), which may be because of the different host diets primarily geophytes on basalt and primarily roots in the chalk population.

Although the microbiome composition between the two populations is not clearly differentiated (Supplementary Figure 4), the functional compositions showed separated clusters (Figure 2A; Supplementary Figure 9), which may be due to short differentiation time (Supplementary Figure 26) and sharp divergent host diet (Nelson et al., 2013; Lövy et al., 2015; Zmora et al., 2019). Higher functional diversity and abundance of the basalt population (Figures 2B,C) reflect higher food diversity resources in basalt (Hadid et al., 2013). The more abundant (Figures 3F,G) and diverse (Figure 3G; Supplementary Figure 9) CaZymes in chalk populations may be caused by adaptation to poorer food root quality. Mannans are a major component of plant cell walls (La Rosa et al., 2019) and cannot be hydrolyzed by the host itself (Teng and Kim, 2018). We found genes enriched in Mannan binding (GO:2001065; Figure 3D) and catalytic activity (GO:0140097) in chalk, which may be due to higher fiber in the chalk root diet. In basalt populations, the enrichment of transporter activity (GO:0005215), cation transmembrane transport (GO:0022857) and glycogen binding (GO:2001069) reflected better food conditions (Figure 3D). Furthermore, growth (GO:0040007) and oxidoreductase complex (GO:1990204) were also enriched in basalt (Figure 3D), which implies gut microbes play an important role in supplying energy to the host (Tremaroli and Bäckhed, 2012; Zhang et al., 2016). These results correspond to higher metabolic rates and activity rates in basalt populations (Hadid et al., 2013; Zhao et al., 2016). KEGG pathway enrichment analysis also showed similar results (Figure 3E) in basalt, the enrichment of oxidative phosphorylation, photosynthesis and pyruvate metabolism (Figure 3E) revealed basalt populations required more energy and had a higher metabolic rate. The enrichment of fatty acid degradation and glycerolipid metabolism suggested the basalt population consumed a high-fat diet. The metabolic pathways of various other substances (Figure 3E e.g., purine metabolism, thiamine metabolism) also indicated the diversity of basalt food resources. These results suggest that the host diet is a main driver of the gut microbiome divergence.

The divergence of soil microbiome between the chalk and basalt is larger than gut microbiome. At the species level (Mukherjee et al., 2022), we found that the soil bacteria differentiation between basalt and chalk is larger than that of the same species in the gut (Figure 4D), which may be due to a combination of more contrasting edaphic stresses and longer divergence time of the soil microbiome. The divergence of the soil microbiome started when volcanic eruptions formed new habitats 1 million years ago, but the hosts of blind mole rats split much later only 0.228 Mya which hampered the gut microbial divergence (Supplementary Figure 4). As the bacteria from the same host individual are under the same intestine stresses, we can expect that the effects of host on the gut microbiome are similar. However, some species displayed significant divergence between basalt and chalk (Supplementary Figures 2, 12–33), which also echoed the results above (Figures 4A,C; Supplementary Figures 12–31).

Generally, the host blind mole rats showed separate clusters in protein between the basalt and chalk populations (Figure 5B). We identified 127 significantly differential proteins (Figure 5C). The enrichment results showed that these differential proteins were involved in the metabolism and synthesis of a wide range of species (Supplementary Figure 34) and were mainly secreted proteins (Supplementary Figure 35). These results correspond to the different diets of the two populations. One individual from the basalt population was clustered into the chalk group, possibly due to the proteome differentiating more slowly than the genome.

The community composition between the chalk and basalt gut microbiome was not significantly different (Supplementary Figure 4) when compared to its host phylogeny (Figure 1D), host diet (Figure 2) or the environmental soil metagenomes (Mukherjee et al., 2022), which may be due to the fact that divergence between the chalk and basalt microbiome occurred much later (Li K. et al., 2015). However, the main question still stands: which environmental factor is the main driver of the microbial divergence between the chalk and basalt, or have all of them shaped the gut microbiome together? For host phylogeny, we showed that the two mole rats populations diverged in the genome (Li K. et al., 2015), methylome (Li et al., 2020b), transcriptome, and genomic editing (Li et al., 2016) and even the proteins (Figure 5B,C; Supplementary Figures 34, 35). If the gut microbiome were shaped mainly by host phylogeny and changed synchronously, we would expect the host phylogeny to mirror that of the gut microbiome (Brooks et al., 2016; Perofsky et al., 2019). However, we found that the host phylogeny was not consistent with the gut microbiome (Figure 5A). Additionally, the mantel test showed there were no correlations (r = 8.078e-3, p = 0.4528) between the host genetic matrix and Bray-Curtis distance matrix based on the gut microbiome; this may be due to the short differentiation time of the gut microbiome (Supplementary Figures 12B, 13C–31C), or that the host phylogeny does not affect it. The significant differences in functional composition (Figure 3) illustrates that the gut microbiome had adapted to the host’s local ecology which is dramatically different between the calcareous chalk and siliceous basalt (Figures 3D–G; Supplementary Figures 9–11). Basically, the same species or similar species compositions in the gut microbiome of blind mole rats can perform different functions, which was consistent with our finding of differentiation of the probably same species between the two types of environments (Figure 4; Supplementary Figures 12–31). The microbe from the two contrasting soils was significantly different in both community composition and function (Mukherjee et al., 2022). The divergence of the gut microbiome was significantly smaller than that of the soil microbiome (Figure 4); this is likely because the soil difference is larger than that of the gut environment. Another possibility that should be explored in the future is that the bacteria in the microbiome underwent sympatric speciation following their hosts, This possibility should be studied in the future in representative dominant bacteria in the microbiome of both hosts.The chalk and basalt rocks and soils are abutting but contrasting with different chemicals (Grishkan et al., 2008, 2009). The community composition of the soil is significantly different from that of the gut (Figure 2A). Both the procruster analysis (Figure 2B) and mantel test (r = −0.01518, p = 0.5317) showed there was no significant correlation between the soil and gut microbiome, allowing us to reject the hypothesis that the local environment was the main factor shaping the gut microbiome.

The food store for S. galili chalk is from Eryngium sp., poaceae roots, and Ranunculus sp. leaves; while the S. galili basalt mainly feed on geophytes, including Hordeum bulbosum, Bellevalia sp., Iris histrio, and a very small amount of Eryngium sp. (Jiao et al., 2021). The geophytes from the basalt are rich in protein and fat, while it is rich in cellulose in bushlets roots in chalk (Lövy et al., 2015, 2017). The derivative S. galili basalt split from its ancestral S. galili chalk, and so does the gut microbiome, which may explain why there was no significant differences in species diversity and composition (Supplementary Figures 4, 5). Unless bacteria also underwent sympatric speciation like their hosts, a major issue that should be explored in the future. Nevertheless, in terms of species and function such as Proteobacteria (Figure 2D) and CAZymes composition, characteristics of the gut microbiome adapted to their host’s diet were also found. Differences in diet may also be involved in the divergence of blind mole rats (Martin, 2013; Nelson et al., 2013; Li K. et al., 2015; Ekimova et al., 2019). SS of the hosts, blind mole rat, and functional adaptation rather than the composition of the gut microbiome were expected to be synchronized. Additionally, here we demonstrated that function or species composition of the gut microbiome are evolving to adapt to their local ecology dramatically different between chalk and basalt ecologiest (Figures 1E, 2C,D, 3A, 4A). Although we can make a preliminary conclusion that host diet is the main driver of functional divergence of the gut microbiome, further experimental tests are required to preclude the roles of phylogeny and environment on the gut microbiome differentiation. Likewise as indicated above the likely possibility of sympatric speciation of the gut microbiome followed that of their hosts.