Bidirectional regulation of the brain-gut axis in Macaca mulatta : implications for wildlife conservation and experimentation

To establish an animal model of psychological stress, we selected wild-born, unhabituated rhesus monkeys and provided them with a captive breeding environment in collaboration with a wildlife conservation organization. Following captivity, these monkeys exhibited significant depression or emotional instability (Fig. 1A). We assembled the quality-controlled data separately for each sample (Table S1; Fig. S1A). There were no significant differences between the wild and captive groups in terms of Chao1, Shannon, and Simpson indices, indicating that both groups had high gene richness, evenness, and gene diversity (Fig. S1B). Additionally, there were no significant differences between the samples of the wild and captive groups, and the dispersion of the samples was within an acceptable range (Fig. S1C). Further analysis revealed that the species accumulation curves of both the wild and captive groups showed an upward trend at the end and tended to flatten, indicating that the sampling volume was sufficient to reveal the species diversity of the community (Fig. S1D). The Chao1 index of the wild group was significantly higher than that of the captive group (P = 0.041), indicating that wild macaques had more microbial species in their feces compared to captive macaques (Fig. S2A). However, there were no significant differences across groups in the Shannon and Simpson indices, indicating that the differences in species evenness between wild and captive macaques were not significant (Fig. S2B). Moreover, it was found that there were significant differences in the microbial community composition between the wild and captive groups (Fig. S2C and D), with extremely significant differences between the two groups (R² = 0.1752, P = 0.008, Fig. S2E). Taken together, captive-induced psychological stress led to behavioral disturbances and significant gut microbiota alterations, including reduced microbial richness and distinct community composition.

To clarify stress-induced gut microbiota changes, we characterized intergroup microbial composition differences. After species annotation, we obtained 8,702 species annotations, including 7,653 bacterial species, 409 archaeal species, 142 fungal species, and 377 viral species. We analyzed the community structure of the fecal microbiota of macaques in the wild and captive groups at the phylum, genus, and species levels. At the phylum level, the fecal microbiota of macaques in both the wild and captive groups was mainly dominated by Firmicutes, Bacteroidetes, and Proteobacteria (Fig. 1B). At the genus level, the fecal microbiota of macaques in both groups was primarily composed of Faecalibacterium, Prevotella, and Treponema (Fig. 1C). At the species level, the fecal microbiota of macaques in both groups was mainly dominated by Faecalibacterium_prausnitzii, Treponema_succinifaciens, and Dysosmobacter_welbionis (Fig. 1D). Analysis of the differences in fecal microbial species of the groups revealed that there were 4,293 common microbial species, with 491 unique species in the wild group and 250 unique species in the captive group (Fig. 1E). Using the linear discriminant analysis (LDA) effect size (LEfSe) test (P < 0.01, LDA > 2), we analyzed the differences in abundance between the wild and captive groups at the genus and species levels of the gut microbiota. The results showed that there were 100 genera and 160 microbial species with significant differences between the two groups. The dominant genera in the captive group were Odoribacter, Ornithobacterium, and Enterobacter, while the dominant genera in the wild group were Pasteurella and Succinivibrio (Fig. 1F). At the species level, the wild group exhibited significantly higher abundance of bacteria involved in host dietary digestion and energy metabolism, particularly cellulose-degrading and short-chain fatty acid-producing taxa (e.g., Faecalibacterium prausnitzii and Succinivibrionaceae Aeromonadales). In the captive group, the bacteria with significantly higher abundance, which were related to the host's dietary digestion and energy metabolism, capable of fermenting various organic substances and producing short-chain fatty acids, included Lactobacillus animalis and Intestinimonas_butyriciproducens (Fig. 1G; Fig. S2F). Collectively, our findings indicate distinct gut microbiota profiles between wild and captive macaques, with captive specimens showing reduced microbial diversity, while wild populations maintained higher abundances of fiber-degrading species linked to plant-rich diets.

To investigate the impact of different psychological stress levels on the functional distribution of gut microbiota in macaques, we performed functional annotation on the non-redundant genes predicted from the fecal microbiota of captive and wild macaques. Functional enrichment analysis revealed that the metabolic pathways in both groups were primarily enriched in global and overview maps, carbohydrate metabolism, amino acid metabolism, membrane transport, and metabolism of cofactors and vitamins (Fig. 2A). Based on Bray-Curtis distance, the genes in the captive and wild groups formed distinct clusters (Adonis test, R² = 0.08, P = 0.04), indicating significant differences between the two groups (Fig. 2B). The abundance of metabolic pathways based on functional annotation changed, with particularly noticeable differences in the pathways with the highest relative abundance. High-abundance pathways such as K07133 and K21572 exhibited higher abundance in the captive group (Fig. 2C). To further explore the specific functional differences between the two groups, we conducted enrichment analysis using the module database for the captive and wild groups. The results showed that the gut microbiota of the wild group had higher gene enrichment in 11 pathways, including sphingolipid biosynthesis; and globin and globin-like proteins; methane metabolism; O-antigen nucleotide sugar biosynthesis; alanine, aspartate, and glutamate metabolism; cysteine and methionine metabolism; and amino sugar and nucleotide sugar metabolism, compared to the captive group. In contrast, the gut microbiota of the captive group had higher gene enrichment in five pathways, including photosynthesis, phosphate and phosphonate metabolism, butanoate metabolism, atrazine degradation, and secondary bile acid biosynthesis, compared to the wild group (Fig. 2D). Therefore, captivity significantly affects the metabolic functions of the gut microbial community in macaques, with pronounced differences in abundance observed in several key metabolic pathways. Together, captivity induces significant metabolic reprogramming in macaque gut microbiota, with wild populations showing enrichment in sphingolipid and amino acid metabolism pathways, while captive specimens exhibit upregulation of butanoate metabolism and secondary bile acid biosynthesis.

To investigate how psychological stress influences host metabolism and biological functions through gut microbiota, we conducted further research. Carbohydrate metabolism, a vital pathway for host energy metabolism, exhibits high gene enrichment in related metabolic pathways. Principal coordinate analysis (PCoA) based on Bray-Curtis distance revealed that the composition of carbohydrate-active enzymes (CAZymes) in captive and wild macaques formed distinct clusters, indicating differences across groups (Fig. 3A). Using the Carbohydrate-Active Enzymes (CAZy) database, we predicted the composition of CAZymes in the gut microbiota of macaques. After classification and annotation, the major CAZymes were identified as glycoside hydrolases (GHs), glycosyltransferases (GTs), carbohydrate-binding modules (CBMs), carbohydrate esterases (CEs), polysaccharide lyases, and auxiliary activities (Fig. 4B). Statistical Analysis of Metagenomic Profiles (STAMP) analysis was employed to screen for the top 30 most abundant CAZymes with significant differences between the captive and wild groups. In the captive group, the abundance of genes encoding 21 glycoside hydrolases (e.g., GH13, GH3, GH97), 2 carbohydrate-binding modules (CBM20, CBM32), 2 carbohydrate esterases (CE0, CE11), and 3 glycosyltransferases (GT4, GT30, GT3) was significantly higher than in the wild group (P < 0.05; Fig. 3B). Conversely, the abundance of genes encoding two glycoside hydrolases (GH1, GH43_24) and one carbohydrate-binding module (CBM13) was significantly higher in the wild group than in the captive group (P < 0.05, Fig. 3C). This suggests that captive macaques possess a richer gene repertoire for carbohydrate-active enzymes compared to wild macaques.

PCoA based on Bray-Curtis distance also showed that the composition of antibiotic resistance genes in captive and wild macaques formed distinct clusters, indicating differences between the two groups (Fig. 3D). Comparison with the Comprehensive Antibiotic Resistance Database (CARD) revealed that the most abundant functional genes annotated were related to tetracycline, macrolide, aminoglycoside, cephalosporin, and fluoroquinolone antibiotics (Fig. 3E). Statistical analysis of the relative abundance (transcripts per million [TPM] values) of these genes showed that the abundance of tetracycline and macrolide resistance genes was significantly higher in the captive group than in the wild group (P < 0.05). In contrast, the abundance of aminoglycoside, penicillin, and fluoroquinolone resistance genes was significantly higher in the wild groups than in the captive groups (P < 0.05, Fig. 3F). In summary, captive macaques exhibit distinct gut microbial carbohydrate metabolism profiles with significantly enriched glycoside hydrolase genes and elevated tetracycline/macrolide resistance genes compared to wild counterparts, reflecting diet-driven microbial adaptation.

We further conducted untargeted metabolomics analysis on fecal samples from captive and wild macaques to investigate the differences in metabolite composition under different stress environments. Qualified samples were subjected to metabolite profiling (Fig. 4A). All detected metabolites were enumerated(Fig. 4B) and categorized by their primary metabolic pathways (superpathway). Lipid metabolites were further classified by sub-pathway (Fig. 4C). The results showed that among all fecal metabolites in macaques, the highest number of metabolites belonged to amino acid metabolism, with a total of 121 species, accounting for 22.92% of the total metabolites. Biosynthesis of other secondary metabolites had 98 species (18.56%), and lipid metabolism had 74 species (14.02%). Statistical analysis of the relative abundance of metabolites in the captive and wild groups revealed that the most abundant metabolites were amino acids, peptides and their analogs, and fatty acyls (Fig. 4D). The findings indicate that amino acid metabolism and biosynthesis of other secondary metabolites are the most dominant metabolic pathway categories in macaque fecal metabolites. This likely reflects the characteristics of gut microbial metabolism and host-microbiota interactions, providing an important basis for further research into their physiological functions and metabolic features.

To further clarify whether psychological stress alters gut metabolites, we analyzed the fecal metabolic composition using metabolomics. Firstly, we established principal component analysis (PCA) models for the wild and captive groups and validated them using permutation tests, finding that the samples in both groups clustered tightly and had high correlation, indicating minimal systematic error and good reproducibility of the experiments (Fig. 5A and B). Subsequently, we conducted univariate analysis to compare differential metabolites in the feces of the wild and captive groups. The results showed that, compared with the wild group, the captive group had 999 metabolites with significantly increased levels and 959 metabolites with significantly decreased levels (Fig. 5C and D). To further analyze the expression patterns of these differential metabolites, we performed clustering analysis on their expression levels. The results revealed that a total of 11,958 differential metabolites in the fecal samples of the two groups exhibited significant intergroup differences in expression levels, and these differential metabolites showed consistent expression patterns within the samples of the wild and captive groups, respectively (Fig. 5E). To more intuitively reveal the synergistic or antagonistic relationships between different types of differential metabolites, we constructed a chord diagram of differential metabolite correlations based on the Spearman correlation coefficients between metabolites (|r| > 0.8 and P < 0.05). The results showed that the expression levels of coumarins and derivatives and 1,4-dioxanes were negatively correlated, while the remaining groups generally exhibited positive correlations (Fig. 5F). To further explore the differences in high-abundance metabolites across groups, we analyzed the 20 metabolites with the largest differences between the groups. The results showed that the metabolites with high expression in the wild group included the human metabolome database (HMDB) 0029218 (urolithin C), 9.062_417.29982, and 9.028_315.19611 (15-deoxy-12,14-prostaglandin A2), while the metabolites with higher expression levels in the captive group included 3.512_206.11780 (2-piperidinobenzoic acid), PB00581 (angoletin), and HMDB0002085 (syringic acid) (Fig. 5G and H). These results demonstrate profound gut metabolite alterations, including stress-associated compounds like syringic acid (captive-enriched) and anti-inflammatory urolithin C (wild-enriched).

We performed metabolic pathway enrichment analysis on the differential metabolites, and the results showed that many differential metabolites were enriched in pathways such as dopaminergic synapse, retrograde endocannabinoid signaling, bile secretion, and steroid hormone biosynthesis (P < 0.01, Fig. 6A). To further analyze the overall changes of these differential metabolites in the pathways, we analyzed the top 10 significantly enriched metabolic pathways in the two groups. The results showed that all metabolites in the alanine, aspartate, and glutamate metabolism, alcoholism, amphetamine addiction, and cocaine addiction pathways were significantly upregulated in the captive group. In contrast, all metabolites in the linoleic acid metabolism and retrograde endocannabinoid signaling pathways were significantly downregulated in the captive group. In addition, the overall expression levels of metabolites in the dopaminergic synapse pathway were upregulated, while those in the steroid hormone biosynthesis, tyrosine metabolism, and bile secretion pathways were downregulated (Fig. 6B; Table S2). Finally, we conducted receiver operating characteristic (ROC) analysis on the differential metabolites, and the results showed that the area under the curve (AUC) value of metabolite C10757 was 0.999, indicating that this metabolite has high discriminability and is suitable as a biomarker to distinguish between wild and captive macaques (Fig. 6C). To reveal the sources of differential metabolites and the metabolic pathways through which microbes produce the relevant metabolites, we used the MetOrigin platform to trace the sources of differential metabolites. The results showed that 126 metabolites were annotated to the KEGG and HMDB databases: four were from the macaque host, 25 were derived from gut microbes, 33 were co-metabolites of microbes and host, and the remaining 207 metabolites were from other sources (Fig. 6D). Based on metabolic source pathway enrichment, there were 5, 24, and 50 enriched pathways for host, microbial, and microbial-host co-metabolism, respectively. Among them, 15 pathways were significantly enriched (P < 0.05). The microbial-specific metabolic pathway was polycyclic aromatic hydrocarbon degradation, while the host was only significantly enriched in the steroid hormone biosynthesis pathway. In the host-microbe co-metabolism pathways, tyrosine metabolism, alanine, aspartate, and glutamate metabolism and others were significantly enriched (Fig. 6E). In summary, metabolic pathway analysis revealed captivity-induced perturbations in neuroactive and endocannabinoid pathways, with microbial-host co-metabolism driving 50 enriched pathways, including tyrosine metabolism, while the high-accuracy biomarker C10757 effectively discriminates between the wild and captive states.

Wild macaques were studied at Nanwan Monkey Island in the Lingshui Macaque Nature Reserve, Hainan Province, where 100 fresh fecal samples Ling Shui (LS) were collected from 30 individuals. In consideration of wildlife welfare and conservation ethics, captive rhesus macaques were selected from a wildlife rescue center, matching the wild macaques in age, weight, health status, and diet (plants, fruits, etc.), collecting 90 fresh fecal samples BaoTing (BT) from 30 individuals. Each sample was collected within 2 min of defecation, with the uncontaminated internal portion placed into sterile cryovials, then preserved in liquid nitrogen and transported to the laboratory for storage at −80°C. Throughout the sampling process, no methods potentially harmful to the macaques were employed, and the study complied with relevant laws and regulations, including the Wildlife Protection Law and the Nature Reserve Management Regulations.

Data filtering was performed using SOAPnuke (v1.5.0), data alignment with Bowtie2 (2.2.5), and data processing with Samtools (1.2). Metagenomic assembly was conducted using MEGAHIT, functional annotation with DIAMOND, and species annotation with Kraken2. MetOrigin platform (http://metorigin.met-bioinformatics.cn↗) was utilized for metabolite tracing analysis. Data analysis and visualization were carried out using R (4.1.2), with reference databases including KEGG, HMDB, CARD, CAZy, Compound Discoverer v.3.3, mzCloud, ChemSpider, and BioSecure Multimodal Database (BMDB).

The central portion of fecal samples, uncontaminated by the external environment, was extracted and stored at −80°C. Genomic DNA was extracted following the MagPure Stool DNA KF Kit B (MFDP-0101) protocol, and qualified samples were used for library construction. The DNBSEQ platform was employed for metagenomic analysis of macaque fecal microbiota.

Library construction workflow: (i) DNA concentration and integrity were assessed, with qualified samples proceeding with library preparation. (ii) A Covaris sonicator was used to fragment a specified amount of metagenomic DNA. (iii) To focus sample bands around 200 bp–400 bp, magnetic bead-based fragment selection was performed on the fragmented samples. (iv) A reaction system was configured to repair double-stranded cDNA ends and append an A base to the 3´ end, followed by adapter ligation to the DNA. (v) A PCR reaction system was set up and programmed to amplify the ligation products, with magnetic beads used for product purification and recovery. (vi) PCR products were denatured into single strands, and a circularization reaction system was configured and reacted at an appropriate temperature to obtain single-stranded circular products. Linear DNA molecules were not circularized, yielding the final library. (vii) Pre-sequencing concentration detection was performed on the circularized products. (viii) DNBSEQ0T7, PE150 sequencing was conducted. DNA extraction and sequencing were completed on the BGI Genomics proprietary sequencing platform.

Raw sequencing data underwent filtering to obtain clean data. MEGAHIT software was used to assemble the metagenome from the quality-controlled clean data. Metagenomic genes were predicted de novo using Meta Gene Mark, and redundant gene prediction results from each sample were processed using CD-HIT software. Gene expression levels were quantified using Salmon software, with TPM values representing normalized gene abundance.

Non-redundant genes were functionally annotated using DIAMOND (BLASTP), with species profiling performed via Kraken2/Bracken. Functional analyses employed KEGG (metabolic pathways), CAZy (carbohydrate-active enzymes), and CARD (antibiotic resistance genes). Alpha (Chao1/Shannon/Simpson) and beta (Bray-Curtis/Jensen-Shannon) diversity metrics were computed in R, while differential abundance was assessed using DESeq2 and LEfSe permutational multivariate analysis of variance (PERMANOVA for PCoA validation). Pathway enrichment (reporter scores) and CAZy/antibiotic resistance genes (ARG) comparisons (STAMP; Wilcoxon/Kruskal-Wallis, P < 0.05) were visualized via extended bar charts and boxplots.

Corresponding samples for metabolomics were thawed at 4°C, then processed through grinding, ultrasonication in a 4°C water bath, static placement at −20°C, centrifugation, and resuspension. An appropriate amount of supernatant was taken as a quality control sample to assess the repeatability and stability of liquid chromatograph mass spectrometer (LC-MS) analysis.

Metabolites were separated and detected using a Waters UPLC I-Class Plus (Waters, Milford, MA) coupled with a Q Exactive high-resolution mass spectrometer (Thermo Fisher Scientific, Bremen, Germany). Chromatographic conditions were as follows: for positive ion mode, mobile phase A (0.1% formic acid in water) and mobile phase B (0.1% formic acid in methanol); for negative ion mode, mobile phase A (10 mM ammonium formate in water) and mobile phase B (10 mM ammonium formate in 95% methanol). Solvents were eluted in a gradient, with a flow rate of 0.35 mL/min, column temperature at 45°C, and injection volume of 5 µL. Mass spectrometry conditions included a scanning mass-to-charge ratio range of 70–1,050, with the top 3 precursor ions selected for fragmentation at energies of 20 eV, 40 eV, and 60 eV. Ion source parameters were as follows: sheath gas flow rate of 40, auxiliary gas flow rate of 10, spray voltage of 3.80 in positive ion mode and 3.20 in negative ion mode, ion transfer tube temperature of 320°C, and auxiliary gas heater temperature of 350°C.

Raw mass spectrometry data were imported into Compound Discoverer 3.3 (Thermo Fisher Scientific, USA) software. The exported results were further processed in metaX for data preprocessing and subsequent analysis. Preprocessing included probabilistic quotient normalization to obtain relative peak areas and batch effect correction using locally weighted polynomial regression fitting QC-based Robust LOESS Signal Correction (QC-RLSC). Compounds with relative peak area coefficients of variation greater than 30% in QC samples were removed.

Principal component analysis was performed in R to reduce data dimensions and ensure data quality. Metabolites were classified and annotated using the HMDB and KEGG databases. Variable importance in projection (VIP) values were calculated to assist in the selection of metabolic biomarkers. Data were log2-transformed to establish partial least squares discriminant analysis (PLS-DA) models for intergroup comparison and differential metabolite screening. Differential metabolites were identified using VIP values >1 in combination with univariate statistical analysis (P < 0.05). Results were visualized using R. Differential metabolites were traced using the MetOrigin platform. Spearman and Pearson correlation analyses were conducted in R to examine the relationships between differential microbes and metabolites in Hainan macaque feces.

Differences were considered statistically significant at P < 0.05 and highly significant at P < 0.01. Spearman and Pearson correlation analyses were performed in R to analyze the relationships between differential microbes and metabolites in Hainan macaque feces, with results visualized in plots.