Characterizing the rhythmic oscillations of gut bacterial and fungal communities and their rhythmic interactions in male cynomolgus monkeys

In this work, eight male cynomolgus monkeys with ages ranging from 2 to 3 years were selected and used to collect feces every 6 h during a 24-h period (five time points) (Fig. 1a). The measurement of routine blood indexes indicated that all these monkeys were in good health (Fig. S1). Finally, 40 fecal samples were collected and used for 16S rRNA gene sequencing (Fig. 1a). To minimize the effects of sequencing depth on α- and β-diversity measure, the number of 16S rRNA amplicon sequencing from each sample was rarefied to 67,279, which yielded an average Good’s coverage above 99% (Fig. S2A). As for microbial richness and diversity, the Ace and Shannon index showed no obvious rhythmicity during a diurnal cycle (Fig. 1b). A Venn diagram showed that 734 amplicon sequence variants (ASVs) were shared by the 16S rRNA amplicon sequencing data of five time points, while 843, 829, 888, 817, and 843 ASVs exclusively belonged to ZT6, ZT12, ZT18, ZT0, and ZT6, respectively (Fig. 1c). Principal coordinate analysis (PCoA) of weighted UniFrac distances showed a significant difference between the bacterial microbes of ZT12 and the other time points, with the principal component 1 (PC1) axis explaining a substantial proportion of the variability (Fig. 1d).

Linear discriminant analysis (LDA) effect size (LEfSe) was performed to identify the distinguishing bacterial taxa between ZT6 (A), ZT12 (B), ZT18 (C), ZT0 (D), and ZT6 (E) (LDA score >2.0). As shown in Fig. S3, Alistipes and uncultured_bacterium_g_Solobacterium were enriched at ZT6 (A), Prevotella and Fournierella were abundant at ZT12 (B), Prevotellaceae_bacterium_DJF_LS10 and Enterococcus were enriched at ZT18 (C), Rhodocyclaceae was abundant at ZT0 (D), and Christensenellaceae and unclassified_g_Monoglobus were enriched at ZT6 (E).

To further explore the oscillation patterns of bacterial microbes, the relative abundance of microbiota at the phylum, genus, and species levels were analyzed. At the phylum level, Firmicutes and Bacteroidota showed obvious rhythmic oscillations (Fig. 1e). Specifically, the relative abundance of Firmicutes decreased at ZT12, while the relative abundance of Bacteroidota increased at the same time point. At the genus level, the genera Prevotella, norank_f_Eubacterium_coprostanoligenes_group, and Peptococcus showed obvious rhythmic oscillations (Fig. 1f). Correspondingly, the bacterial microbes unclassified_g_Prevotella, unclassified_g_norank_f_Eubacterium_coprostanoligenes_group, and uncultured_bacterium_g_Peptococcus exhibited similar rhythmic oscillations (Fig. 1g). By performing the bacteria–bacteria correlation analysis, we found that the average relative abundance of Firmicutes was negatively correlated with Bacteroidota (Fig. 1h). At the genus level, the average relative abundance of Prevotella was negatively correlated with the genera norank_f_Eubacterium_coprostanoligenes_group and Peptococcus (Fig. 1i).

Based on these results, we concluded that most of the bacterial microbes of male cynomolgus monkeys exhibited rhythmic oscillations and were closely correlated with each other during a diurnal cycle.

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) was used for the function analysis of bacterial communities at ZT6 (A), ZT12 (B), ZT18 (C), ZT0 (D), and ZT6 (E) (Fig. S4). By focusing on the bacterial functions that showed obvious rhythmicity, we found that the relative abundance of ABC transporters, two-component system, quorum sensing, propanoate metabolism, flagellar assembly, and valine, leucine, and isoleucine degradation were decreased at ZT12 and positively correlated with the relative abundance of norank_f_Eubacterium_coprostanoligenes_group during a diurnal cycle (Fig. 2a and b). Differently, the relative abundance of alanine, aspartate, and glutamate metabolism, fructose and mannose metabolism, folate biosynthesis, pentose and glucuronate interconversions, cyanoamino acid metabolism, and riboflavin metabolism were increased at ZT12 and positively correlated with the relative abundance of Prevotella during a diurnal cycle (Fig. 2c and d).

Based on the 40 fecal samples rhythmically collected from the eight male cynomolgus monkeys, here, we performed ITS amplicon sequencing to investigate the rhythmic oscillations of intestine fungal microbes (Fig. 3a). To minimize the effects of sequencing depth on α- and β-diversity measures, the number of ITS amplicon sequencing from each sample was rarefied to 26,513, which still yielded an average Good’s coverage above 99% (Fig. S2B). As to the microbial richness and diversity, the Ace and Shannon index showed no obvious rhythmicity (Fig. 3b). The Venn diagram showed that 138 ASVs were shared by the ITS amplicon sequencing data of five time points, while 277, 182, 178, 224, and 141 ASVs exclusively belonged to ZT6, ZT12, ZT18, ZT0, and ZT6 (Fig. 3c). PCoA of weighted UniFrac distances showed no obvious difference between the fungal microbes of five time points (Fig. 3d).

LEfSe was performed to identify the distinguishing fungal taxa between ZT6 (A), ZT12 (B), ZT18 (C), ZT0 (D), and ZT6 (E) (LDA score >1.0). As shown in Fig. S5, Aspergillus penicillioides, Oidiodendron cereale, Talaromyces islandicus, and Penicillium capsulatum were enriched at ZT6 (A); Piskurozymaceae, Solicoccozyma, and Phaeoacremonium were abundant at ZT12 (B); Herpotrichiellaceae was enriched at ZT18 (C); and Acremonium sordidulum and Exophiala were abundant at ZT6 (E). Of note, no fungal taxa abundant at ZT0 (D) was detected.

Next, the relative abundance of fungal microbes was analyzed to investigate their oscillation patterns. As shown in Fig. 3e through i, the relative abundance of Ascomycota, Aspergillaceae, Aspergillus, Talaromyces, Aspergillus penicillioides, and Talaromyces funiculosus was decreased at ZT18 and exhibited slight rhythmic oscillations.

FUNGuild was used for the function prediction of fungal communities at ZT6 (A), ZT12 (B), ZT18 (C), ZT0 (D), and ZT6 (E). By focusing on the top 9 fungal functions, we found that the relative abundance of undefined saprotroph decreased at ZT18 (C) (Fig. 4a and b) and showed a similar oscillation mode with the fungal genus Aspergillus (Fig. 3h). By performing the correlation analysis between the fungal genera and fungal function, we found that the relative abundance of Aspergillus was positively correlated with the relative abundance of undefined saprotroph (Fig. 4c). Apart from undefined saprotroph, we also found that the relative abundance of fungal parasite and plant pathogen exhibited different oscillation patterns and showed no obvious correlation with the fungal genera (Fig. 4d and e).

We performed Spearman correlation analysis to uncover the bacterial–fungal correlations during a diurnal cycle (Tables S1 to S5). As shown in Fig. S6, 11 fungal genera (unclassified_f_Bionectriaceae, Chaetomium, Bionectria, Botryotrichum, Kiflimonium, Acaulium, Chloridium, Clonostachys, Solicoccozyme, Cephalotrichum, and Pseudombrophila) formed a positive and high-magnitude correlation network at the A (ZT6) time point. Besides, among the five time points, the bacterial and fungal microbes formed the most complicated interactions at the B (ZT12) time point.

By calculating the number of total, negative, and positive interactions between gut microbes, we found that the most positive bacteria–bacteria interactions (205) were observed at the B (ZT12) time point, while the negative bacteria–bacteria interactions (100) were abundant at the A (ZT6) time point (Fig. 5a). For the fungi–fungi correlation analysis, the most positive (119) and negative (16) interactions were observed at the A (ZT6) time point (Fig. 5b). For the bacteria–fungi correlation analysis, most of the positive interactions (78) occurred at the B (ZT12) and C (ZT18) time points, while the negative interactions (114) were abundant at the B (ZT12) time point (Fig. 5c). Of note, the fungal or bacterial genera that positively or negatively correlated with the top 5 bacterial or fungal genera were analyzed at the five time points, respectively (Fig. 5d through h).

The limitation of our study is that we mainly focus on the male cynomolgus monkeys. Since the gut microbes underwent sex-dependent rhythmic alterations in rodents (13, 14), our results might not be applicable to the female cynomolgus monkeys, and further investigation is needed.

Our study uncovers the diurnal oscillation patterns of bacterial and fungal microbes in male cynomolgus monkeys and investigates their correlations during a diurnal cycle. Compared with the fungal microbes, the bacterial microbes exhibited obvious rhythmic oscillations, especially for the genera norank_f_Eubacterium_coprostanoligenes_group and Prevotella whose relative abundance varied significantly at ZT12 (19:00). Consistent with this result, the gut bacterial functions were also oscillated significantly and closely correlated with the abovementioned two bacterial genera. By characterizing the bacterial and fungal interactions during a 24-h diurnal cycle, we found that the bacterial and fungal microbes closely interacted with each other, and the most bacteria–fungi interactions occurred at ZT12 (19:00).

In this work, the eight captive male cynomolgus monkeys (Macaca fascicularis), with ages ranging from 2 to 3 years (the lifespan for cynomolgus monkeys was about 20–30 years), were housed at the Songjiang Non-human Primate Facility of Institute of Neuroscience. These monkeys were housed in an air-conditioned environment with controlled temperature (22°C ± 1°C), humidity (50% ± 5% RH), 12-h light/12-h dark cycle (lights-on time 07:00–19:00), and continuous access to municipal water. As for the diet, these monkeys were fed rhythmically with commercial monkey diet (Anmufei, Suzhou) twice daily (200 g per monkey at 8:00 and 15:00) and fruits and vegetables once daily (100 g per monkey at 10:00) to provide essential nutrition and vitamins.

Before the experiment, all the monkeys were housed singly for 2 months and had no antibiotic exposure. Then, the fecal samples of the single-caged monkeys were collected rhythmically at ZT6 (13:00), ZT12 (19:00), ZT18 (1:00), ZT0 (7:00), and ZT6 (13:00) during a 24-h diurnal cycle. The size of single cages used in this study was 0.7 m × 0.9 m × 1 m.

The fresh fecal samples of cynomolgus monkeys were collected by experienced working staffs at ZT6 (13:00), ZT12 (19:00), ZT18 (1:00), ZT0 (7:00), and ZT6 (13:00) during a 24-h diurnal cycle. The time used to collect the fecal samples at each of the time points was about 1 h. The collected samples were frozen immediately with liquid nitrogen in sterile tubes and stored at −80°C until use.

The total genomic DNA samples were extracted by using the Omega Soil DNA Kit (M5635-02) (Omega Bio-Tek, Norcross, GA, USA). The quality and concentration of DNA were determined by 1.0% agarose gel electrophoresis and a NanoDrop ND-2000 spectrophotometer (Thermo Scientific Inc., USA) and kept at −80°C prior to further use. For 16S rRNA amplicon sequencing, the hypervariable V3-V4 region of 16S rRNA was amplified using the primer pairs 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) by an ABI GeneAmp 9700 PCR thermocycler (ABI, California, USA). For ITS amplicon sequencing, the fungal ITS1-2 regions were amplified by PCR using the following primer pairs ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2R (5′-GCTGCGTTCTTCATCGATGC-3′). The PCR mixture included 4 µL 5 × Fast Pfu buffer, 2 µL 2.5 mM dNTPs, 0.8 µL each primer (5 µM), 0.4 µL Fast Pfu polymerase, 10 ng of template DNA, and ddH₂O to a final volume of 20 µL. The PCR amplification cycling conditions were as follows: initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 45 s, and single extension at 72°C for 10 min, and end at 4°C. Then, the PCR products were purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using Quantus Fluorometer (Promega, USA). Equal amounts of amplicons were used for paired-end sequencing using the Illumina MiSeq PE300 platform (Illumina, San Diego, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

Raw FASTQ files were de-multiplexed using an in-house Perl script, quality-filtered by fastp version 0.19.6 (48), and merged by FLASH version 1.2.11 (49) with the following criteria: (i) the 300-bp reads were truncated at any site receiving an average quality score of <20 over a 50-bp sliding window, the truncated reads shorter than 50 bp were discarded, and reads containing ambiguous characters were also discarded; (ii) only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence. The maximum mismatch ratio of the overlap region is 0.2. Reads that could not be assembled were discarded. (iii) Samples were distinguished according to the barcode and primers, and the sequence direction was adjusted, with exact barcode matching and two nucleotide mismatch in primer matching. The high-quality sequences were de-noised using DADA2 and viewed as ASVs. The taxonomic assignment of ASVs was performed using the Naive Bayes consensus taxonomy classifier implemented in Qiime2 and the SILVA 16S rRNA database (v138) and ITS database (UNITE 8.0). The bacterial and fungal function was predicted by PICRUSt2 (50) and FUNGuild based on ASV representative sequences.

Bioinformatic analysis of the gut microbiota was carried out using the Majorbio Cloud platform (https://cloud.majorbio.com↗). The rarefaction curves and α-diversity indices including ACE, Shannon index, and Good’s coverage were calculated with Mothur v1.30.2 (51). The similarity among the microbial communities in different samples was determined by PCoA based on weighted UniFrac distances using the Vegan v2.4.3 package. The PERMANOVA test was used to assess the percentage of variation explained by the treatment along with its statistical significance using the Vegan v2.4.3 package. A Venn diagram was generated using R package “VennDiagram.” The LEfSe (52) (http://huttenhower.sph.harvard.edu/LEfSe↗) was performed to identify the significantly abundant taxa (phylum to species) of bacteria and fungi among the different groups.

Gephi software (version 0.9.4) was used to perform the correlation network. Spearman correlation analysis was used to uncover the bacterial–fungal correlations at ZT6, ZT12, ZT18, ZT0, and ZT6. A connection is indicated for Spearman’s correlation with a coefficient >0.6 (positive correlation) or <−0.6 (negative correlation) and a significant (P < 0.05) correlation. In the correlation network, bacteria are labeled as a yellow node; fungi are labeled as a green node. The size of each node is proportional to the relative abundance. The red lines represent positive correlations between the nodes, and the green lines represent negative correlations, with the line width indicating the correlation magnitude.

The statistical analysis data were presented as the mean ± SEM. The statistical significance between two time points was analyzed using the Student t-test (unpaired, two-tailed). All statistical analyses were calculated by using GraphPad (version 7.0).