Acute Sleep-Wake Cycle Shift Results in Community Alteration of Human Gut Microbiome

Fecal samples were collected by the 22 study subjects 1 day before the experiment, to determine the baseline microbiome (baseline, T0), on the morning after the sleep-wake cycle shift (shift, T1), and 2 days after the subjects returned to their regular routine (recovery, T2) (Fig. 1A). The samples were sequenced by 16S rRNA gene amplicon sequencing, and 1,671 microbiota operational taxonomic units (OTUs) were identified by using Qiime2 (16). By comparing the OTUs that were detected repeatedly (detected in more than 2 samples) in the three time points, we observed only a marginally significant difference in the numbers of OTUs between baseline and shift (Mann-Whitney U test, P = 0.057) and between shift and recovery (Mann-Whitney U test, P = 0.065) (Fig. 1B). The OTUs were then collapsed to the genus level, and a significantly higher number of detectable genera (with at least 10 reads in no fewer than 3 samples at each time point) was observed at shift compared with the numbers at both baseline (Mann-Whitney U test, P = 0.006) and recovery (Mann-Whitney U test, P = 0.005). The number of genera at the recovery stage had decreased and was comparable to the baseline level (Fig. 1C). This result indicated a tendency for the number of genera to increase after the sleep time shift; however, this increase was quickly reversed after participants returned to their normal routines (recovery stage) (Fig. 1C). The alpha diversity was not significantly different among the three groups (Fig. 1D). Principal coordinate analysis based on the weighted UniFrac index exhibited no significant difference in bacterial composition among groups (Fig. 1E). Additionally, we calculated the Bray-Curtis similarity values between T0 and T1, T0 and T2, and T1 and T2 in each individual, and no significant differences were observed (Fig. 1F). The phyla Bacteroidetes, Firmicutes, and Proteobacteria dominated across the three time points (Fig. 1G). Given that changes in the ratios of the phyla Firmicutes and Bacteroidetes (F/B) have been associated with obesity (17) and type 2 diabetes (11), we compared the F/B ratios among any two time points. An increase in the F/B ratio at T1 (after sleep-wake time shift) compared with that at T0 (baseline) was observed (paired, one-sided Mann-Whitney U test, P = 0.049) (Fig. 1H).

There were no significant changes in taxa observed by statistical test among the baseline, shift, and recovery samples in our data set, which is in accordance with the results from a previous study of sleep restriction in rats and humans (2). Thus, we focused on the microbes whose average relative abundance changed more than 2-fold between any two of the time points to characterize the tendency of alteration, instead of identifying statistically differential abundance. Distinct patterns in changes of abundance were observed (Fig. 2). The phyla Fusobacteria and Tenericutes and classes Fusobacteriia and Mollicutes were increased after sleep shift and decreased to be comparable to the baseline level at recovery (Fig. 2A and B), representing the taxa that are resilient in response to acute sleep-wake cycle shift. The order Pasteurellales and family Clostridiales_Peptostreptococcacea gradually decreased from baseline to recovery, indicating a long-term impact of acute sleep-wake shift on them (Fig. 2C and D).

Next, to explore the alteration in function of the fecal microbiota at different time points, PICRUSt2 (phylogenetic investigation of communities by reconstruction of unobserved states) was used to infer the abundances of the MetaCyc pathways from the 16S rRNA gene amplicon sequencing. There was no alteration in pathways between time points T0 and T2. However, in the comparison between T0 and T1, five pathways were significantly enriched at T1, including ADP-l-glycero-β-d-manno-heptose biosynthesis, allantoin degradation IV (anaerobic), acetyl-coenzyme A (CoA) fermentation to butanoate II, acetylene degradation, and the superpathway for purine deoxyribonucleoside degradation (Fig. 3A). Four pathways were decreased at T2 compared with their abundances at T1, i.e., purine nucleotide degradation II (aerobic), the superpathway for l-methionine biosynthesis (by sulfhydrylation), acetyl-CoA fermentation to butanoate II, and the superpathway for N-acetylneuraminate degradation (Fig. 3B). Collectively, purine metabolism was increased during acute circadian rhythm shifts. Deregulated purine metabolism has been reported in multiple diseases, including gout, colitis, autism, and Alzheimer disease (AD) (18–21). The pathway for acetyl-CoA fermentation to butanoate II was significantly increased at T1 and significantly decreased at the recovery stage. Butanoate metabolism is related to short-chain fatty acid (SCFA) metabolism, which has been reported to be involved in energy metabolism (22), inflammation (23), and psychological functioning (24) of the host. Additionally, we also observed that most of the pathways that were increased at T1 had a tendency to be decreased at T2, though not statistically significantly (Fig. 3C).

Inspired by the idea that a slight change in abundance can introduce microbiome structural alteration (25), analysis of the cooccurrence network of genera at the three time points was conducted. The correlations with an r of ≥0.6 are shown in Fig. 4. The cooccurrences among Alistipes, Barnesiella, and Odoribacter were preserved at baseline (T0), shift (T1), and recovery (T2), and the interaction between Bacteroides and Parabacteriodes remained intact during the acute circadian rhythm shifts. Notably, the cooccurrences between Butyricimonas and other genera were mainly observed at T0 and T2 but not at T1. Members of Butyricimonas are butyric acid-producing bacteria. Butyric acid is an important energy source for intestinal epithelial cells and plays a role in the inflammatory process and maintenance of colonic homeostasis (26). A previous study in mice has implicated circadian rhythm disturbance as affecting the intestinal microbiota, which may have implications for inflammatory diseases (27). However, the taxa that cooccurred with Butyricimonas were different at T0 and T2, suggesting the existence of alterations in microbiome organization even after the return to the regular routine.

Next, the NetShift method (28) was used to analyze the driver genera that hold key significance in microbial interactions. The genus Odoribacter was predicted as a driver of the community change from T0 to T1 (Fig. 5A and D). Odoribacter is a producer of short-chain fatty acids like acetate, butyrate, and propionate (29), and studies have reported that increased Odoribacter may partly contribute to decreased inflammation (30) and may be associated with colon cancer (31). In addition, Odoribacter was also predicted to drive the community change from T1 to T2 (Fig. 5B and D), suggesting that the acute circadian rhythm shifts may introduce transient alteration in the gut microbiome interaction network through Odoribacter. Bacteroides is predicted as a driver to promote community change between T0 and T2 (Fig. 5C and D). Bacteroides is a well-known genus that is implicated in metabolic diseases, including obesity and diabetes (32, 33). For example, the depletion of species from the Bacteroides genus in obese individuals was proposed to be related to the higher concentrations of aromatic amino acids and branched-chain amino acids in circulation, which are known risk factors for type 2 diabetes and cardiovascular disease (34). This result suggested that the acute sleep-wake cycle shift may exert an effect on the interactions between gut microbes even after the return to the normal routine.

This study was approved by the Ethics Committee of Yixing People’s Hospital of Jiangsu University (ethics no. 2016005). All participants were provided a written informed consent upon enrollment. Once the consent forms were signed, we screened them for eligibility criteria and sent questionnaires and sample collection kits to participants.

Twenty-two healthy adults aged from 22 to 35 years were recruited from Nanjing Medical University and Yixing People’s Hospital (). All subjects had received no antibiotics, probiotics, prebiotics, or antifungal medications for 3 months prior to fecal sample collection, and none of them were on anti-inflammatory or antioxidant drugs. The participants were asked to maintain their regular lifestyle throughout the study, including the timing and type of meals, except for the sleeping time as we required. Table S1

Fecal samples of each subject were obtained at home as directed in the instructions provided and frozen immediately until delivered to the laboratory. Fecal samples were delivered together with dry ice and were stored at –80°C at the laboratory until extraction.

About 100 mg of sample was weighed in the laboratory from each fecal sample and used to extract total-genome DNA, following the protocol of the DNA extraction kit (product number DP328; Tiangen Company, Beijing, China). The concentration and purity of the extracted bacterial DNA were detected using a Qubit 2.0 fluorometer (Thermo Scientific, USA). DNA quality and quantity were determined by agarose gel electrophoresis.

The genomic DNA in fecal samples was extracted using the DNA extraction kit (product number DP328; Tiangen Biotech Co., Ltd., Beijing, China). The DNA concentration was measured using a Qubit 2.0 fluorometer (Thermo Fisher Scientific, USA). The 16S rRNA gene amplification procedure was divided into two PCR steps. In the first PCR, the hypervariable V4 region of the 16S rRNA gene was amplified using the conserved primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Then, amplification was performed in 96-well microtiter plates. Reactions were run with the following program: 2 min of denaturation at 94°C, followed by 30 cycles of 20 s at 94°C (denaturing), 30 s at 56°C (annealing), and 40 s at 68°C (elongation), with a final extension at 68°C for 5 min. The PCR products were quantified using the Quant-iT PicoGreen quantification system (Life Technologies), and samples with a concentration above 6 ng/μl were diluted to ∼3 to 6 ng/μl prior to further analysis. In the second PCR step, sequencing primers and adaptors were added to the amplicon products. The PCR was run according to the cycling program described above except with a reduced cycle number of 15. The amplification products were purified with Agencourt AMPure XP beads (Beckman Coulter Genomics, MA, USA) according to the manufacturer’s instructions. Equimolar amounts of the amplification products were pooled in a single tube. The pooled DNA samples were concentrated using the DNA Clean & Concentrator-5 kit (ZymoResearch, Irvine, CA, USA) according to the manufacturer’s instructions. The concentration of the pooled libraries was determined using the Quant-iT high-sensitivity DNA assay kit (Life Technologies) following the specifications of the manufacturer. Amplicon sequencing was performed on the Illumina HiSeq system (BGI, China). Automated cluster generation and 2 × 250-bp paired-end sequencing with dual-index reads were performed. The sequencing output was generated as demultiplexed fastq-files for downstream analysis.

Bioinformatics analysis of 16S rRNA gene amplicons was performed by using Qiime2 (version 2018.6.0) (16). Briefly, fastq reads were processed by using the dada2 program, and the dada2 denoise-paired command was used to delete the low-quality ones. Dada2 generates unique features that could be compared between different studies. The taxonomy of these features was assigned via the Greengenes reference database (version 13-2) classifier with 99% similarity. Determinations of alpha and beta diversities were also conducted in Qiime2.

The functional capacity of the gut microbial community was predicted using PICRUSt2. The OTU table generated from Qiime2 was rarefied according to the lowest sample sequence (10,214) and then supplied to PICRUSt2. Predicted functional genes were categorized into MetaCyc pathways (45) and compared across each of the two time points using LEfse. Statistical differences in MetaCyc frequencies were determined with the default parameters (linear discriminant analysis [LDA] score of >2).

The SparCC algorithm (46) was used to estimate the correlations from the compositional network. The network was demonstrated by cytoscape with edges connecting nodes (bacterial taxa) with a correlation coefficient of over 0.6 or less than −0.6. NetShift (28) analyses were conducted in the web server using the networks generated by SparCC.

All raw 16S rRNA gene sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) under accession number SRP227414↗.