Introduction

Lignocellulose, as the most abundant biomass component on earth, composed of cellulose, hemicellulose and lignin, is hydrophobic and difficult to enzymatically degrade [1]. Bamboo is a low-nutrition food with about 70–80% of its composition being lignocellulosic components [2, 3]. In general, mammals lack enzymes that breakdown complex lignocellulosic carbohydrates [4, 5], and heavily rely on the complementary action of the lignocellulose-degrading enzyme repertoire of their microbial symbionts [6, 7]. Bamboo rats and giant pandas are two famous specificity bamboo-eating animals and the understanding of their adaptations and mechanisms for digesting lignocellulose-based diet evokes enormous interest.

The giant panda retains a typical carnivore digestive system and the adaptation of their gut microbiome to a specialized bamboo diet has been widely studied [2, 8,9,10,11]. Several studies have documented, through metagenomic analyses, the ability of the gut microbiome of giant pandas to digest cellulose and hemicellulose [12,13,14,15]. It has been reported that the abundance of microbial genes in the gut of giant pandas contributing to lignocellulose degradation is lower than that of ruminants and other herbivores [16, 17], and is affected by diet, habitat environment (captive or wild) and seasonal factors (e.g., climate and food availability) [3, 13, 18, 19]. Gut microbial composition differs between captive and wild giant panda populations, and the gut microbes of wild giant pandas seem to be more capable of lignocellulose digestion [13, 20, 21]. In contrast, only a limited number of studies have focused on the bamboo rats. Research has shown that an endoglucanase with cellulase activity was identified from a bacterium strain in the gut of the bamboo rat [22], and rich CAZy families have been detected in the cecum microbiome that might contribute to the degradation of bamboo fibers [23].

Few studies have systematically characterized the lignocellulose utilization from an integrated view of both the host and gut microbiome, with most mainly focusing on the microbiome originating from stool. Stool and gut samples present different microbial communities in both membership and diversity, and the proportions of aerobic genera are elevated in stool samples [24]. The fact that most microbes involved in lignocellulose degradation are anaerobic [25], generates the need for metagenomic analyses of gut instead of stool samples. Furthermore, a single metagenomic analysis does not completely reveal the gene products of the gut microbes and their impact on host health [26]. This necessitates multi-omics techniques that enable the efficient characterization of dynamic changes and functional activities of both the gut microbiome and the host [27, 28].

Bamboo rats serve as an excellent model to study the utilization of lignocellulose due to their specialized bamboo-eating dietary characteristics. Here, we explored the degradation of recalcitrant lignocellulosic materials by complementary and synergistic cooperation between the host and their gut microbiome. We performed a study, which involved stages from the unweaned period to the specific bamboo-eating period in bamboo rats, and integrated multi-omics measurements including microbiomes from several gut segments, serum metabolome, and the host transcriptome to identify lignocellulose digestion mechanisms. In addition, similarities and differences in the gut microbiome response to dietary switches in bamboo rats and giant pandas were compared.

Results

Comparative profiling of gut microbial patterns in bamboo rats and giant pandas

We grouped samples according to the dietary stage of the bamboo rats and giant pandas (Fig. 1a, b; Tables S1 and S2) and compared their gut microbial compositions and associated functional profiles. Clusters of gut microbiome were observed to be closer in the unweaned giant pandas and bamboo rats based on Bray–Curtis dissimilarities at the OTU level, whereas huge differences were evident after these animals started consuming bamboo-based diets (Fig. 1c). Linear discriminant analysis (LDA) effect size (LEfSe) was employed to assess the microbiome attributes that differed significantly by dietary status. Our analysis revealed that the family Bacteroidaceae and genus Bacteroides were enriched in the milk diet stage of both the bamboo rat and the giant panda (Kruskall–Wallis test, p < 0.05, LDA > 3.5) (Fig. S1a, b; Table S3). After switching to the lignocellulosic diet, microbial families including Muribaculaceae, Lachnospiraceae, Ruminococcaceae, and Desulfovibrionaceae were enriched in the gut of the bamboo rats, whereas the phylum Proteobacteria, and the families Streptococcaceae and Clostridiaceae were observed to be significantly elevated in the gut of the giant pandas (Fig. S1a, b; Table S3). Notably, the gut microbiomes of captive and wild pandas were quite different at the OTU, CAZyme, KO, and COG levels (PERMANOVA, p < 0.05) (Fig. 1c, d, Fig S1c, d), and the microbial compositions and CAZy family profiles of wild pandas appeared closer to those of the weaned bamboo rats (Fig. 1c, d).

Fig. 1: Comparison of the gut microbial patterns in bamboo rats and giant pandas.
figure 1

a Multi-omics study design for bamboo rats. The study yielded host and microbial data from the liver, blood, intestinal tissue, intestinal content, and stools of bamboo rats with different ages. Animal and sample counts are shown. b Metadata of samples from two representative bamboo-eating animals. Data for bamboo rats were generated in this study and data for giant panda were accessed from public databases (see “Materials and Methods”). Principal coordinates analysis (PCoA) plots based on Bray–Curtis dissimilarity at (c), OTU level and (d), CAZy family level. Triangles represent bamboo rats and circles represent giant pandas. Each color corresponds to a different dietary stage. Ellipses are at the 70% confidence level.

Full size image

Adaptation of the gut microbiome to dietary switch in bamboo rats

The within-sample diversities based at the OTU and gene levels were calculated for the duodenum, cecum, and colon of the bamboo rats. In general, OTU and gene diversities of the gut microbiome in the S2 stage were higher than in the S1 stage (Fig. 2a and Table S4). The OTU diversity in the hindgut was lower than that of the foregut, but gene diversity in the hindgut was higher than that of the foregut (Fig. 2a and Table S4). PCoA plot based on Bray–Curtis dissimilarities showed that their gut microbiomes separated into three distinct clusters (PERMANOVA, R2 = 0.06, p < 0.05) (Fig. 2b). Moreover, non-metric multidimensional scaling (NMDS) analysis of the genes (PERMANOVA, R2 = 0.08, p < 0.05), CAZymes (PERMANOVA, R2 = 0.08, p < 0.05), KOs (PERMANOVA, R2 = 0.12, p < 0.05), and COGs (PERMANOVA, R2 = 0.17, p < 0.05) also showed clearly different clusters between the S1 and S2 stage samples (Fig. S2). Our findings indicate that substantial changes in gut microbial composition and function occur after diet conversion.

Fig. 2: Differences in the gut microbiome of bamboo rats before and after their diet switch.
figure 2

a Boxplots show microbial α diversity (Shannon index) at the OTU and gene levels at different ages. Samples from different intestinal segments are marked in different colors. b PCoA plot based on Bray–Curtis dissimilarities at the OTU level. Ages are distinguished by color, and the samples from different intestine segments are distinguished by shape. Ellipses are at the 80% confidence level. c Co-occurrence networks based on the OTU level. Gut microbiome in the two dietary stages were analyzed separately. Each node represents one OTU, and each edge represents a strong (Spearman |r| > 0.6) and significant correlation (FDR p < 0.01) between the two nodes. The size of each node is proportional to the degree of the OTUs. Total sum scaling normalized topological metrics of the networks are shown as radar chart. d Bar chart shows the average relative abundance of the predominant microbes in the S1 and S2 stages. The bacterial order Bacteroidales has a shift between the S1 and S2 stages. e Boxplot shows the bacterial order Bacteroidales related to the dietary types detected by MaAsLin2 based on the OTU level. All the samples in the S1, S2 and adult stage were used in this analysis.

Full size image

Next, co-occurrence networks were established to estimate microbial co-existence in the different dietary stages. According to the topology of these networks, Desulfovibrionaceae-centered microbial communities in the S1 stage were noted to be replaced by a stronger inner-connected microbial community containing Lachnospiraceae, Ruminococcaceae, and Muribaculaceae in the S2 stage (Fig. 2c). In addition, the order Bacteroidales was found to be predominant in both the stages. Notably, Bacteroidaceae, which were dominant in the S1 stage (average relative abundance: 31.6%), were largely replaced by uncultured Muribaculaceae in the S2 stage (average relative abundance: 38.2%) (Fig. 2d). This bacterial switch in the order Bacteroidales was significantly affected by dietary types (MaAsLin2, FDR p < 0.05) (Fig. 2e).

Synergistic actions of multiple CAZy families are necessary for bacteria to perform lignocellulose degradation. Changes in the relative abundance of CAZy families that were significantly affected by multiple factors, including intestinal segment, age, and dietary type are displayed in Fig. 3a and Table S5 (MaAsLin2, FDR p < 0.05). Overall, the relative abundance of CAZymes in the hindgut (colon and cecum) was higher than in the foregut and produced two different patterns consistent with dietary stage (Fig. 3a). We then focused on the CAZymes which were up-regulated in the S2 stage. These enzymes mainly originate from Muribaculaceae, Lachnospiraceae, Clostridiaceae, and Bacteroidaceae, which indicated the rich carbohydrate enzymatic functions of these microbes. Among these microbes, Muribaculum gordoncarteris, Lachnospiraceae bacterium GAM79, Roseburia intestinalis, Duncaniella dubosii, and Flintibacter sp were five representative species with high abundance (Fig. 3b). According to the structure of lignocellulose, we identified CAZymes involved in all steps of lignocellulose degradation in the gut microbiome of the bamboo rats (Fig. 3c), and these enzymes were significantly elevated in the S2 stage (MaAsLin2, FDR p < 0.05) (Table S5). For example, some peroxidases and laccases, including AA4 and AA10, required for lignin degradation and release of cellulose and hemicellulose, were increased in the S2 stage (Fig. 3a and c). Some redox enzymes and carbohydrate esterases including CE2, CE3, CE4, CE6, and CE12, necessary for lignocellulose degradation were also enriched in the S2 stage (Fig. 3a and c). Moreover, the endo-glucanases (i.e., GH5, GH9, GH10, GH53, GH81, GH85), endo-hemicellulases (i.e., GH53, GH98), exo-hemicellulase (i.e., GH39) and β-glucosidases (i.e., GH128), which contribute to cellulose and hemicellulose degradation, were noted to be enriched in the S2 stage (Fig. 3a and c).

Fig. 3: CAZymes profile adapted to lignocellulose utilization.
figure 3

a Heatmap of the CAZymes profile. CAZymes associated with dietary type, age and intestinal segment detected by MaAsLin2 are shown. CAZy families with up-regulated abundance in the S2 stage are marked. Relative abundance of each CAZy family is colored according to the row z-score. b Sankey diagram depicting the distribution of the top 20 bacterial species associated with bamboo-based diet ranked by CAZymes abundance. CAZymes are colored according to CAZy classes and bacteria are colored according to microbial families. Height of the rectangles indicates the relative abundance of the CAZymes (left) and species (right). c Model for lignocellulose degradation in the gut of bamboo rat (modified from Bredon et al [6]). CAZy families enriched in the S2 stage with corresponding enzymatic functions are listed.

Full size image

In addition, we also found functional changes in vitamin and amino acid synthesis in the gut microbiome of bamboo rats that accompanied the food conversion. Biosynthesis of vitamin B family (VB1, VB2, B3, VB7, and B9) and amino acids (including arginine, BCAAs and AAAs) were enhanced, while retinol and vitamin B6 metabolism were decreased in the S2 stage (Wilcoxon rank sum test, p < 0.05) (Fig. S3).

Serum metabolite patterns at the different dietary stages

A total of 2530 metabolites (1431 and 1099 metabolites for positive and negative polarity mode, respectively) were detected as blood circulating metabolites from the 24 bamboo rats. Partial least-squares discriminant analysis (PLS-DA) was used to reduce the dimensionality of the high-dimensional metabolome data. The score plot showed a clear separation between each group (Fig. 4a). Accordingly, we compared the differences in metabolites between the groups in the two dietary stages (Fig. 4b). Overall, 91 shared differential metabolites (SDMs), including 52 for positive mode and 39 for negative mode, were identified between the two dietary stages (VIP > 1.0, FC > 2.0 or FC < 0.5, p value < 0.05) (Fig. 4b).

Fig. 4: Serum metabolome characteristics.
figure 4

a PLS-DA plots of metabolomic features at 4 ages. Positive mode and negative mode of the untargeted metabolome are displayed separately. b Venn plots of the differential metabolite numbers between the paired groups. Differential metabolites were calculated separately based on the positive and negative mode. Differential metabolites were defined if VIP > 1, FC > 2.0 or FC < 0.5 and p value < 0.05. c Correlation between shared differential metabolites and dietary stage (left), and metabolite attributes are shown (right). The attributes of the metabolites, including annotation information (HMDB, KEGG and LIPID MAPS database), metabolic information, and source information have been provided. Scatter plots representing the relationship between (d), bile acid derived, (e), glucose derived metabolites with animal ages.

Full size image

We then assessed the correlations between these SDMs and the dietary patterns to explore potential biomarkers at the different stages. In this study, the SDMs conformed to correlation analysis (Spearman r > 0.7, FDR p < 0.05), and ROC curve estimation (AUC > 0.9) were selected as the potential biomarkers for the specific dietary stages (Table S6). Metabolites that were potential drug components or hormones were excluded. Finally, 27 and 12 potential biomarkers were identified in the S1 and S2 stages (Fig. 4c). For example, 12(13)-DiHOME, the terminal product of linoleic acid metabolism, was observed to robustly correlate with bamboo-based diet (Spearman r = 0.87, FDR p < 0.0001). We classified these metabolites according to their database annotation status (KEGG, HMDB, and LIPID MAPS database) and the attributes of their metabolic derivatives. Specifically, some bile acid derived metabolites (glycocholic acid, 1β-hydroxycholic acid, sulfoglycolithocholic acid, and taurohyocholic acid), and amino acid derived metabolites [roseotoxin S, DL-α-aminocaprylic acid, saccharopine, α-glutamyl-4-hydroxyproline, and ε-(γ-glutamyl)-lysine)], were abundant in the S1 stage (Fig. 4c, d). On the other hand, some glucose derived molecules such as p-cresol glucuronide, 5′-O-β-D-glucosylpyridoxine, 17-β-estradiol-3-glucuronide, and estriol 16-glucuronide were significantly elevated in the S2 stage (VIP > 1.0, FC > 2.0 or FC < 0.5, p value < 0.05) (Fig. 4c and e).

Increased fatty acid metabolism during the lignocellulosic diet period

In order to further explore the mechanisms of host response to diet switch, we attained transcriptomes from 36 samples from 3 tissues (liver, duodenum and colon) across 4 ages (15d, 30d, 45d, and 80d). These tissues are known to be critical for metabolic performance and nutrition utilization. About 92.2% of transcripts represented the available vertebrata conserved orthologs, which indicates that the transcriptomes were well assembled and relatively complete (Fig. S4). Differentially expressed genes (DEGs) were identified in the liver, duodenum, and colon of the bamboo rats across the 4 ages. KEGG pathway and GO enrichment analyses were performed for the DEGs in each pair groups (Figs. 5a, S5). Concretely, pathways related to carbohydrate digestion and absorption (hydrolase activity acting on glycosyl bonds) were enriched in the duodenum in the S2 stage (FDR p < 0.05) (Figs. 5a, S5). Moreover, nutrients absorbed by the intestine at the two dietary stages were different. For example, pathways related to the digestion and absorption of proteins, minerals and fats were enriched in the S1 stage, while pathways related to metabolism and absorption of vitamins and carbohydrates were abundant in the S2 stage (FDR p < 0.05) (Fig. 5a).

Fig. 5: Host response to diet switch.
figure 5

a KEGG pathway enrichment in the duodenum, colon, and liver. Up-regulated pathways in the S1 stage are shown on the left, and up-regulated pathways in the S2 stage are shown on the right. b KEGG pathway enrichment of metabolites with significant and robust correlations in the bamboo-eating stage (Spearman coefficient r > 0.6, FDR p < 0.05). c Scatter plot representing the relationship between linoleic acid (up), arachidonic acid (down) with animal ages. The correlations between metabolites and bamboo-eating stage are shown as Spearman coefficients. d Schematic map of linoleic acid metabolic pathway in the liver. Up-regulated genes and metabolites in the S2 stage are marked in red.

Full size image

Diverse KEGG pathways such as linoleic acid metabolism, arachidonic acid metabolism, retinol metabolism, and gluconeogenesis were enriched in the S2 stage (FDR p < 0.05) (Fig. 5a). Linoleic acid metabolism was also observed to be the most enriched KEGG pathway, based on serum metabolites, which strongly correlates to the bamboo-based diet (Spearman r > 0.6) (Fig. 5b). Arachidonic acid, a product of linoleic acid metabolism, had an elevated level in the S2 stage (Fig. 5c). However, only linoleic acid correlated strongly to bamboo-based diet (Spearman r = 0.65) (Fig. 5c). In the liver, cytochrome P450 family members CYP2C, CYP2J, CYP2E1, and CYP3A4 showed significant upregulation in the S2 stage, which facilitate efficient conversion of linoleic acid to the end product 12,13-DiHOME (Fig. 5d). Our results revealed that the source and end products of linoleic acid metabolism were significantly elevated in the S2 stage, which might be associated with adaptation of the bamboo rats to lignocellulosic diet.

The associations of gut microbiome and serum metabolome

We next explored the associations of gut microbiome and serum metabolome in bamboo rats that might facilitate lignocellulose utilization. We assessed the effect of duodenum, cecum and colon microbiomes on the serum metabolome, using the Mantel test (Fig. 6a). Overall, the microbiomes of the hindgut (cecum and colon) were more tightly coupled with serum metabolome than the foregut (duodenum) (Fig. 6a). For example, the correlations between the duodenum CAZymes profile and serum metabolome were weak with an extremely low effect size (Mantel test, r = 0.05, p > 0.05; r = 0.04, p > 0.05 for positive and negative metabolome modes, respectively), while only the colon microbiome, both in taxonomy profile and function profile, was significantly associated with serum metabolites (Mantel test, p < 0.05) (Fig. 6a). This indicates that colon microbial metabolism and associated products contributed most to the serum metabolome.

Fig. 6: The relationship between gut microbiome and host serum metabolome.
figure 6

a Mantel test quantifying variance explained between gut microbiome and serum metabolome. Taxonomy and CAZymes profiles of the gut microbiome from different intestinal segments were compared with serum metabolome. Mantel statistic r is provided in the rectangular area in the plot. *p < 0.05; **p < 0.01; ***p < 0.001. b Network plot shows the correlation between microbes based on OTU level and circulating metabolites. Microbiome are represented by diamonds, and metabolites are represented by circles. Red line: similarity values ρ > 0.60; blue line: similarity values ρ < −0.60.

Full size image

To identify co-variation between the gut microbes and serum metabolites, we used an integrative approach to analyze the colon microbiome (based on the OTU level) and serum metabolome (Fig. 6b). Microbes closely related to host metabolites reflected two categories; (a) the Muribaculaceae family, which was predominant in the S2 stage, and (b) a group of co-occurring microbes (i.e., families Lachnospiraceae and Ruminococcaceae) with abundant metabolic functions (Fig. 6b). As expected, short-chain fatty acid (SCFA)-derived metabolites, including butyl acrylate and butanoate-derived molecules, displayed strong significant positive correlation with the second group of co-occurring microbes (ρ > 0.6). Biomarkers of bamboo-based diet (Fig. 4c), namely 5-amino-6-(D-ribitylamino) uracil, margaric acid, and p-cresol glucuronide, were positively correlated with Muribaculaceae (ρ > 0.6). Notably, some polyunsaturated fatty acids (linoleic acid and arachidonic acid) had weak correlations with the gut microbes (ρ < 0.6) indicating that they were mainly derived from food and lacked a microbial origin.

Discussion

Multi-omics covering stages from milk diet to bamboo diet revealed that the bamboo rat is well adapted to lignocellulose-based food. The cooperation between their gut microbiomes and their own metabolic systems is displayed in Fig. 7.

Fig. 7: Putative cooperate mechanisms of bamboo rat and gut microbiome in response to the low-nutrient bamboo diet.
figure 7

In the bamboo-eating stage, the dominant gut microbiota were switched from Bacteroidaceae to Muribaculaceae, which enriched diverse CAZymes associated with lignocellulose degradation and increased the production of glucose, amino acids, B vitamins (VB1, VB2, VB3, VB7 and VB9), and SCFAs. Meanwhile, circulating concentrations of some metabolites such as 12(13)-DiHOME, 2-Hydroxyhippuric acid and glucose derived metabolites were elevated, and associated metabolic function were enhanced in liver. The metabolism of fatty acids (i.e., linoleic acid and margaric acid from food) and volatile SCFAs (i.e., butyrate and acetate produced by the gut microbes) could be the main energy sources in bamboo rats.

Full size image

The high-fiber bamboo diet after weaning increased the diversity and richness of the gut microbiome in bamboo rats (Fig. 2a), which is consistent with similar observations in humans [29], mice [30], swine [31], and ruminants [32]. Moreover, the level of serum 2-Hydroxyhippuric acid, a predictor of gut microbiome α-diversity [33], was found to be increased in our study (Fig. 4c). This signifies the powerful ability of fiber in shaping the gut microbiome. We found that the phyla Bacteroidetes and Firmicutes were predominant in the gut of bamboo rats, which have been reported to substantially contribute to the taxonomic and metabolic variations in the gut microbiome of humans and ruminants [34, 35]. Notably, we observed a major switch in the order Bacteroidales from family Bacteroidaceae to family Muribaculaceae, after the consumption of bamboo-based food (Fig. 2d). Recent studies have shown that Muribaculaceae are functionally distinct from their neighboring families and are versatile degraders of complex carbohydrates [36, 37]. We detected that Muribaculum gordoncarteri, Duncaniella dubosii, Duncaniella sp. C9, and Sodaliphilus pleomorphus, representing the Muribaculaceae family, harbor diverse functional CAZymes involved in lignocellulose degradation (Fig. 3b). Although the ruminant gut microbiome can digest recalcitrant lignocellulose, the order Bacteroidales has not been reported to be enriched with members from the Muribaculaceae family [35], and seem to lack cellulosome to degrade the complex plant polysaccharides [38]. Our findings suggest that the mechanisms for lignocellulose degradation by gut microbiome in ruminants and non-ruminants are different.

Our data revealed that the colon microbiome in bamboo rats contributed the most to serum metabolome composition (Fig. 6a). Lignin breakdown is the first step of lignocellulose degradation by gut microbiome with concomitant release of cellulose and hemicellulose [39]. The enrichment of Enterocloster bolteae and other bacterial species, which harbor ligninolytic enzymes such as AA4 and AA10 in the bamboo-eating period mainly contributed to lignin breakdown in the bamboo rats (Fig. 3b, c). Multiple CAZy families with oxidative, hydrolytic and non-hydrolytic activities identified in Muribaculum gordoncarteri, Lachnospiraceae bacterium GAM79, Roseburia intestinalis, Blautia producta, Duncaniella dubosii, and Bacteroides thetaiotaomicron typically act on cellulose and hemicellulose metabolism in the gut of post-weaning bamboo rats (Fig. 3b, c). These findings indicate that the gut microbiome of bamboo rats is functionally adapted to the utilization of complex lignocelluloses. These bacterial members formed a strong inner-connected microbial community with richness in metabolic and fermentative functions (Fig. 2c). As a result, the levels of acetate- and butyrate-derived metabolites in the blood were significantly increased in the S2 stage (Figs. 4c and 6b). Meanwhile, metagenomic analysis revealed that the biosynthesis of vitamin B family (VB1, VB2, B3, VB7, and B9) and amino acids (BCAAs, AAAs and arginine) was significantly enhanced (Fig. S3). The level of 5-Amino-6-ribitylamino uracil, the intermediate in VB2 (riboflavin) biosynthesis, was correspondingly significantly elevated in the blood (Fig. 4c). We noticed that the gut microbes provide a variety of metabolites for the host, including volatile SCFAs, vitamins and essential amino acids during lignocellulose utilization (Fig. 7).

The relative abundance of retinol and VB6 metabolism in the gut microbiome was reduced (Fig. S3), while transcriptome showed that the retinol metabolism pathway was significantly enriched in the gut and liver of bamboo rats in the bamboo-eating period (Fig. 5a). Previous studies reported that VB6 biosynthesis by gut microbes appear to be negatively correlated with the relative abundance of butyrate producers such as Blautia and Roseburia [40, 41], which is consistent with our findings. In contrast, bamboo-based food is rich in β-carotene [42], which can be processed into retinoic acid by intestinal epithelial cells for further metabolism or storage [43]. Retinol, available from the conversion of carotene present in bamboos, might already meet the needs of the bamboo rat host, thus, limiting the requirement for retinol biosynthesis by gut microbiome.

The transcriptome and metabolome analyses revealed that the host linoleic acid metabolism was significantly enhanced in the S2 stage (Fig. 5). Although some bacterial genera like Propionibacterium, Lactobacillus, and Bifidobacterium have shown promising abilities to produce conjugated linoleic acids [44, 45], the level of linoleic acid in the serum of bamboo rats was found not to be closely related with gut microbiome, but rather to robustly correlate with specific bamboo dietary type (Spearman r = 0.65) (Fig. 5c). The plants of the Poacea family, including bamboos, possess linoleic acid or linoleic acid-derived compounds (i.e., linoleic-CoA, linoleic acid ethyl ester and linoleoyl chloride) and margaric acid [46, 47]. This could explain the elevated serum concentrations of linoleic acid and margaric acid in bamboo rats during the bamboo-eating period. Similar trends were observed in bamboo-eating giant pandas where the concentrations of metabolites related to the linoleic acid metabolism predominate in the serum [48]. The above indicates that the intake and utilization of fatty acids like linoleic acid and margaric acid might be one of the adaptive mechanisms of animals to respond to low-nutrient bamboo food (Fig. 7).

It is generally believed that increased intake of linoleic acid is associated with overweight and obesity in humans [49]. However, the body fat rate of bamboo rat with enhanced linoleic acid metabolism during the bamboo eating period is low [50]. The serum metabolite 12,13-diHOME, a terminal product of linoleic acid metabolism [51], was found to be significantly increased in concentration in the S2 stage and identified as an important biomarker of bamboo diet in our study (Fig. 4c). It has been reported that the level of 12,13-diHOME is negatively correlated with the body-mass index (BMI) and insulin sensitivity, which could stimulate fat consumption and reduce blood triglyceride levels [52]. We also noticed that Bacteroides thetaiotaomicron, previously reported to reduce fat accumulation by altering the levels of serum metabolites [53], was one of the main bacterial species involved in lignocellulose utilization in the bamboo rats (Fig. 3b). Moreover, in the bamboo-eating period, the circulating concentrations of glucose-derived metabolites were significantly elevated (Fig. 4e), which might be attributed to increased gluconeogenesis in the liver (Fig. 5a). Therefore, bamboo rats, which feed on high lignocellulose diet, might reduce lipid accumulation by increasing lipolysis and induce gluconeogenesis to maintain glucose levels in the blood. Our observations are consistent with similar findings in ruminant sheep [54]. In addition, volatile SCFAs (especially butyrate) produced by gut microbes, provide 60–70% of the energy needed by the colon, muscle and brain cells in the hosts [55, 56]. Overall, our data indicate that in the process of feeding on low-nutrition bamboo diet, fatty acids (i.e., linoleic acid and margaric acid) available from food and volatile SCFAs produced by the gut microbes (i.e., butyrate, acetate) could be the main energy sources in bamboo rats (Fig. 7).

We compared the gut microbiome adaptation and mechanisms of lignocellulose degradation in bamboo rats and giant pandas, as both are bamboo specialists. The phyla Bacteroides was enriched in the gut of both bamboo rats and giant pandas before weaning; however, aerobes and facultative anaerobes like Proteobacteria and Streptococcaceae were found to be elevated in the gut of giant pandas after weaning (Fig. S1), along with a decrease in the diversity or richness of gut microbiome [10, 13]. Our analysis and previous studies provide evidence that the gut microbes of giant pandas, especially wild pandas, have the ability to metabolize lignocellulose (Fig. 1c, d) [20, 21]. However, the digestion and utilization rate of bamboo fiber in the gut of giant pandas is low [57]. Research has shown that the source of energy for giant pandas may be protein, hemicellulose, or intracellular nutrients [2, 3, 58]; however, the role of fatty acids has been unexplored. In response to the bamboo diet, the genes of giant panda involved in the digestion and utilization of bamboo fatty acids experience adaptive convergence [2], and linoleic acid related metabolites are enriched in their serum [48]. Combined with the results achieved in bamboo rats, we infer that the pattern of using fatty acids (i.e., linoleic acid and margaric acid) from food could be another major source of energy for giant pandas.

In conclusion, we systematically explored the utilization of lignocellulose in bamboo rats for the first time, using a multi-omics approach with data from gut digesta, blood, and tissue samples throughout their different dietary periods. Our study reveals that the diversity of Bacteroidetes involved in the utilization of lignocellulose in bamboo rats is different from ruminants and giant pandas. We comprehensively characterized the bacterial species and functional genes associated with lignocellulose degradation in the intestine of bamboo rats, and highlighted the importance of fatty acid metabolism (i.e., linoleic acid from food and SCFAs from gut microbiome) on a low-nutrient bamboo diet. Further assays such as the isolation of functional bacteria and in vitro or in vivo tests are needed to draw causality relationships. Our study not only provides novel insight into the complex mechanisms for lignocellulose digestion in the non-ruminant bamboo rats, but also provides data to support the future bio-utilization of lignocellulose involving specific bacteria.

Materials and methods

Animals and specimen collection

All protocols and studies involving animals were conducted in accordance with the guiding principles of the Animal Care and Use Committee of South China Agriculture University (permit number 2021g007).

All bamboo rats were purchased from an artificial breeding farm in Guangxi province, China. The captive breeding of this animal was legal before the COVID-19 outbreak, while it is now banned. We defined the pre-weaning period as being the S1 stage (age < 40d), while the post-weaning period before adulthood was the S2 stage (Table S1). Sampling from bamboo rats was conducted on March 28th–29th, 2019. Tissue, blood and intestinal contents were sampled from 24 bamboo rats (6 animals each at the ages of 15d, 30d, 45d, and 80d). For adult bamboo rats (age > 1 year), we only collected fecal samples (Table S1).

Data from captive and wild giant pandas, which cover different dietary periods, including 230 16S rRNA gene sequencing datasets (accession no. PRJNA524253 and PRJNA358755) [10, 20] and 88 shotgun metagenomic sequencing datasets (accession no. PRJEB24780, PRJCA000366, and PRJNA356809) [3, 17, 20], were download from the Sequence Read Archive (SRA) database (Table S2).

16S rRNA gene data analysis

Barcode and primer sequences were removed from the raw reads using Cutadapt (v1.18) [59]. OTU clustering were performed by usearch (v11, -cluster_otus) [60]. Non-redundant reads were compared with the reference Gold database to remove chimera sequence (usearch, -uchime2_ref). The filter of sequences that were not assigned to bacteria, OTU estimation and taxonomy annotation were achieved by alignment against the SILVA database (v132) through the QIIME platform [61]. RDP classifier with 97% sequence identity was used for taxonomic classification.

In the joint analysis of bamboo rats and giant pandas, total OTU counts were rarefied to 5000 per sample due to the magnitude difference in raw counts among the samples. In other cases, the dataset were rarefied to same OTU counts according to sample with the lowest count tags. Shannon indices were estimated for within-sample microbial diversity analysis. Beta diversity was evaluated using principal coordinates analysis (PCoA) based on Bray–Curtis dissimilarity values. Multivariable association of microbial community features with metadata including dietary type, age and intestinal segment were assessed using MaAsLin2 [62]. LEfSe was used to identify and compare different bacterial communities between each group [63].

Co-occurrence network analysis was performed for the S1 and S2 stages, separately, with co-occurrence patterns between the microbes assessed by calculating the correlations based on the cumulative sum scaling normalized OTU tables. Spearman’s rank coefficients (r) between the OTUs were calculated in a pairwise fashion, and significant and robust correlations (FDR  p < 0.01, |r| ≥ 0.6) were used to construct the networks [64]. The R packages vegan, hmisc and igraph were used for network statistical analyses [65, 66]. Network visualization was performed by Gephi software with the layout algorithm of Fruchterman-Reingold [67].

De novo transcriptome analysis

Raw RNA-seq reads were trimmed to remove adapter sequences and excluded low quality reads with a score cutoff < 30, using fastp (v0.19.7) [68]. Sequences that mapped to microorganisms including bacteria, viruses, archaea, and fungi were filtered. Due to the absence of a high-quality genome for Rhizomys pruinosus, we generated a de novo assembly of all the high-quality reads from the 36 samples using Trinity (v2.3.2), with default parameter (except -min_kmer_cov = 2) [69]. Bowtie (v1.0.0) [70] was used to map the reads from each sample. The bam format files were delivered to Corset (v1.07) [71] to cluster and remove redundant transcripts. Transcripts with maximum expression TMM > 1 were retained. Open reading frames (ORFs) were predicted for all transcripts by Transdecoder (v5.5.0) [72]. Transcripts with poorly supported protein-coding evidence (ORF length < 50 amino acids) were discarded. Completeness of the assembled transcripts was assessed by BUSCO (v3.0.2) based on conserved euarchontoglires and vertebrate orthologs [73]. Transcriptome annotation was performed on the predicted protein sequences using Trinotate (v3.1.1) [74] with default parameters.

We used the R package DESeq2 (v1.30.1) [75] to obtain differentially expressed genes (DEGs) by comparing gene expression values between each group. Genes with a twofold (or greater) change and an adjusted p < 0.05 were considered as differentially expressed. KEGG and GO enrichment analyses of the DEGs were performed using an online resource (http://www.omicshare.com/) with default instructions.

Metagenomic assembly and gene profile construction

Clean reads were achieved by removing the adapter and low quality sequences (quality score < 30) from the raw reads using fastp (v0.19.7) [68]. We constructed a nucleotide dataset containing the genomes of giant panda (AilMel_1.0), hoary bamboo rat (RhiPru_v1_BIUU), rat (Rnor_6.0), mouse (GRCm38.p6), Cavia procellus (Cavpor3.0), maize (RefGen_v4), and our de novo assembled transcriptomes. The clean reads were filtered by mapping reads to this merged nucleotide dataset using BWA-MEM (v 0.7.17) [76].

Megahit (v1.0.3) was used to assemble the high-quality reads [77]. To improve the assembly result for less abundant species in the gut of bamboo rats, the final metagenomic assembly results were combined with the results of each sample assembly, and the results of the two sample groups (S1 and S2) assemblies. Contigs larger than 500 bp were predicted by Prodigal (v2.6.3) and the generated coding sequences with length < 102 bp were discarded. The initial non-redundant gene sets were clustered by CD-HIT-est (v4.7), with parameter “-c 0.95 -n 10 -G 0 -aS 0.9” [78]. The encoded protein sequences were aligned against the NCBI-NR database (November 2018) on the basis of DIAMOND (v0.8.28.90, diamond blastp -evalue 10 -max-target-seqs 250) by CARMA3 (carma -classify-blast -type p -database p) [79, 80]. The genes classified as “eukaryota” and “unknown” were excluded from the non-redundant gene set, and the final gut microbial gene catalog of bamboo rat included 2.76 M genes.

Functional annotation of metagenomic assembled genes

For functional annotation, the protein sequences were aligned against the KEGG (release 79) [81], eggNOG (v4.5) [82] and CAZy (available November 2018) [83] databases. KEGG and eggNOG annotations were aligned against the corresponding databases using DIAMOND (v0.8.28.90), by taking the best hit with the criteria of E value < 1e−5. Annotation of CAZymes families was aligned against the CAZy database using the hmmscan tool of HMMER (v3.2.1), by taking the best hit with the threshold of E value < 1e−18 and coverage >0.35 [84].

Species annotation of the genes coding for CAZymes

The species origin of the genes coding for CAZymes was identified according to the following pipeline [85]: similarity of the predicted assembled gene segments to known species was estimated by alignment against NCBI-NT (Jul 9, 2021), using BLASTn (v2.11.0, default parameters except that -evalue 1e-10 -outfmt 6 -word_size 16); and result hits >90% identity and >70% alignment coverage were used as the critical value for species assignment. For the genes with multiple best-hits, the species annotation with the highest frequency and the highest average similarity was defined as the annotation of the gene.

Quantification of genes and construction of species and function profiles

Sequenced-based abundance profiling was performed as previously described [86]. For the calculation of relative gene abundance, the high-quality reads from each sample were aligned against the gene catalog by BWA-MEM with the criteria of alignment length ≥50 bp and identity >95%. Copy number of each gene in one sample was defined as the quotient of the number of mapped reads and the length of the gene. The species and function profile were calculated by summing the relative abundance of genes belonging to each category per sample. Overall differences in the genes, CAZymes, KOs, and COGs were evaluated by NMDS or PCoA analysis based on Bray–Curtis dissimilarity values, using the phyloseq and vegan packages in R [66, 87].

Metabolomics data processing

Basic metabolomic data handing including statistical analysis, Pattern Hunter analysis, biomarker analysis, and enrichment analysis were performed with MetaboAnalyst (v5.0) [88]. To reduce baseline noise, the features were excluded if their relative standard deviation were higher than 25% in the QC samples or near-constant throughout all samples detected by interquartile range. Data normalization was achieved by log transformation and uv scaling (mean-centered and divided by the standard deviation of each variable). For screening differential metabolites, we considered the variable importance of projection (VIP) value of the first principal component of the partial least-squares discriminant analysis (PLS-DA), fold change (FC) and statistical significance of the variables calculated by t-test. Differential metabolites were defined if VIP > 1, FC > 2.0 or FC < 0.5, and p value < 0.05. Correlation analysis based on Spearman correlation coefficient (r) was performed against a given pattern consistent with the dietary types (S1 or S2 stage). The metabolites with significant and robust correlation (FDR p < 0.05, |r| ≥ 0.6) were used for KEGG enrichment analysis. Classical univariate ROC (receiver operating characteristic) curve analysis was performed to generate the ROC curve. The excellent predictive performance was evaluated by calculating the area under the ROC curve (AUC) > 0.9 [89].

Estimation of microbiome association with serum metabolome

Mantel tests were used to quantify the covariation between different measurement types [90]. Bray–Curtis dissimilarity matrices from the OTU and CAZymes profiles, and the Euclidean distance matrix from the serum metabolite profile were obtained based on the corresponding normalized matrix. The paired matrices were then compared using the mantel function implemented the R package vegan [66]. The effect size (mantel statistical r) and significance (p value) of the Mantel test were used to quantify the covariation between the microbial profiles of different intestinal segments (duodenum, cecum, colon) and the serum metabolome.

R package mixOmics [91] was used to integrate and explore the metabolome and microbial taxonomy profile. The different types of omics datasets were processed and normalized before integration. Individual and paired analysis with the sPLS-DA (sparse PLS-DA) model was executed to test whether a set of OTUs or metabolites can distinguish the groups. To circumvent spurious associations, the optimum number of components and variables was determined with 50 × 5-fold cross-validation. The similarity values between pair of variables were calculated by the network function. Variable pairs with a high similarity measure (|ρ| ≥ 0.6) are considered as relevant and used for built networks. Networks displaying correlations between the gut microbiome and serum metabolites were generated using Cytoscape (v3.7.2) [92].

Statistical analysis

Bacterial features showing differential abundance between the different dietary stages of bamboo rat and giant panda were identified using LEfSe threshold criteria of LDA score > 3.5 and p value < 0.05. The significance in the difference of community compositions between the groups was analyzed with the permutational multivariate analysis of variances (PERMANOVA) method based on 999 permutations. Wilcoxon rank sum tests were performed to assess the differences in the alpha diversity and functional profiles of the gut microbiomes. Pattern hunter analysis with a threshold criterion of Spearman r > 0.6 and FDR p < 0.05 was used for estimating the metabolites associated with dietary stages. Mantel test with 9999 random permutations were used to assess the contribution of the microbiome of different parts of the gut to the serum metabolome, two measurement types were considered correlated if the p value < 0.05, and the absolute value of mantel statistical r indicates the degree of covariation. The p values in the multiple comparisons were adjusted by the Benjamini–Hochberg false discovery rate (FDR) method.