Host circadian behaviors exert only weak selective pressure on the gut microbiome under stable conditions but are critical for recovery from antibiotic treatment

To test our hypothesis that the gut microbiome is under active selection by the host’s daily feeding patterns, we developed a novel model that combines a clock-knockout strain with constant dim red light (RR) for disrupted rhythmicity of feeding under conditions in which fecal samples could be easily collected without disturbing the light/dark conditions of the experiment. Our protocol for fecal collection comprised LD entrainment for 5 weeks, then RR for 26 weeks, and concluded with a return to LD for an additional 8 weeks (Fig 1A). We used mice in which 2 central circadian clock genes Period1 and Period2 were both knocked out (Per1/2-dko [21]). In the standard circadian mouse C57 background, these Per1/2-dko mice express robustly rhythmic nocturnal activity patterns in LD and arhythmicity and/or highly disrupted activity patterns in constant darkness (DD; [22]), whereas wild-type (WT) mice maintain robust nocturnal activity rhythms in LD and free-running rhythmicity in constant DD conditions with a period approximately 24 h. Per1/2-dko mice are fertile and healthy, unlike another commonly used clock-disrupted mouse strain, the Bmal1-ko mouse, that suffers from a variety of health problems (short lifespan, arthropathy, infertility, etc. [23–26]), and has been used in some other microbiome studies [13]. Moreover, the microbiome responses in Per1/2-dko mice to time-restricted feeding and jet-lag paradigms have already been characterized [14,27].

This study is a longitudinal analysis where each mouse genotype serves as its own control on the basis of comparisons pretreatment versus posttreatment. The treatments include both alterations in the light/dark conditions and/or antibiotics included in the drinking water. We tested the temporal patterns of locomotor activity and feeding of both strains in both LD cycles and uninterrupted constant conditions. For constant conditions, we used dim red light (RR) conditions rather than constant darkness (DD) to visually facilitate fecal sample collection (constant white light is known to suppress rhythmicity, even in WT). As shown in Fig 1B, the wheel-running locomotor activity of WT mice is strongly rhythmic with a 24-h period in LD 12:12, and sustains approximately 24 h but slightly longer period free-running rhythms in RR (period approximately 24.5 h, S1 Table). On the other hand, while the locomotor activity of Per1/2-dko mice remains highly rhythmic in LD (period = 24 h), it is dramatically disrupted in RR, where a variety of fragmented patterns are observed (often in the same mouse) from arhythmic to ultradian to very short period (ostensibly averaging 13.3 h, but note the dramatic loss of power that indicates arhythmicity and/or highly disrupted patterns, S1 Table).

In parallel experiments, we used infrared video recordings to quantify temporal feeding patterns of the WT and Per1/2-dko mice and found that feeding behavior under the LD and RR regimens corresponded closely to the locomotor activity patterns. Under a regimen of LD for 6 weeks, dim RR for 23 weeks and then back to LD for the remaining 10 weeks, we confirmed strong daily rhythms of both locomotor activity and feeding in LD for both strains (Fig 1C); in RR, WT mice retain robust daily rhythms of both kinds of activity, while Per1/2-dko mice exhibit highly disrupted patterns of both activity and feeding behavior in RR that are ultradian or arhythmic (Fig 1C and further analyzed in S2 Table; note the pronounced loss of power in the Per1/2-dko mice in RR). There were no obvious differences in the locomotor activity or feeding behavior between the antibiotic-treated and the untreated groups of mice.

Based upon the temporal patterns of feeding and locomotion in the 2 mouse strains, we tested whether temporally disrupted versus rhythmic behaviors (feeding, etc.) select for differing microbiotal composition and communities using the LD/RR/LD protocol shown in Fig 1A. Prior to treatments, there were no significant differences in the microbiome diversity, neither between the WT versus Per1/2-dko genotypes nor among the pretreatment time points (Table 1; based on the Shannon diversity index that assesses richness and evenness of species [28]). Thus, the pretreatment samples within each genotype serve as controls for the treatment- and time-dependent changes in the microbiomes in this longitudinal design. Quantification of the microbiota using 16S sequence analyses indicated that the relative abundances at the phylum level were not obviously changed by the LD to RR to LD transitions alone (the antibiotic-untreated or “UT” samples, Fig 2A and 2D). The total bacterial load in the UT stool samples did not change significantly among these LD/RR/LD transitions, although there was a downward trend in bacterial load over approximately 8 month timecourse (S1 Fig). To refine our assessments of microbiotal diversity, evenness, and richness of operational taxonomic units (OTUs; from phylum to species), we applied the commonly used Shannon and Chao1 alpha-diversity indices [28–30]. The Shannon diversity index is a measure of biodiversity that deals with richness (the breadth of differences among OTUs) and evenness (abundance of specific OTUs), while the Chao1 index focuses upon richness [28–30]. The alpha-diversity Shannon and Chao1 indices suggest a slight destabilization of microbiome diversity at the beginning of the RR treatment that stabilizes later in RR and continues in the final LD segment (“UT” samples, Fig 2B and 2E for Shannon, Fig 2C and 2F for Chao1). Assessments of diversity with the Chao1 metric (Fig 2C and 2F, plotted on a linear time scale) largely agrees with the Shannon diversity metric (Fig 2B and 2E, plotted on an equidistance timescale). Importantly, while there are some minor differences in the relative abundances at the phylum level between the WT(UT) and Per(UT) samples (Fig 2A versus 2D), there is no obvious difference in the long-term response of the microbiome to rhythmic {WT(UT)} versus arhythmic {Per(UT)} behaviors in RR as indicated by the alpha-diversity indices (compare Fig 2B and 2C with 2E and 2F).

We then proceeded to a complementary assessment, namely beta-diversity, which quantifies the compositional dissimilarity between treatments/environments that in our case are the temporal changes in the gut environment. Differences in microbiotal composition can be visualized by the Bray–Curtis beta-diversity index [29,31] and graphed by principal coordinates analysis (PCoA) dimensionality. Fig 2G and 2H summarize the PCoA analysis for 5 key time points; all 16 time points are shown in S1 Video and S2 Fig (S1 Video is organized as a video animation). By this complementary beta-diversity assay, the long-term microbiotal compositions are not changing significantly in the WT(UT) samples as the mice transition from LD to RR and back to LD (rhythmic behavior throughout; Figs 1, 2G, and S2). Surprisingly, neither is the microbiotal composition changing in the Per(UT) samples that are undergoing a transition from rhythmic behavior (first LD) to temporally disrupted behavior (RR) and back to rhythmic behavior (second LD) (Figs 1, 2H, and S2). Based on the results of prior short-term studies with other arhythmic mouse models, “jet-lag” simulations, and other tests of circadian disruption [13–17], this is an unexpected result, the implications of which will be addressed in the Discussion. Note that these alpha- and beta-diversity indices do not necessarily indicate that microbiotal composition was completely unchanged by the LD/RR/LD transitions. This issue will be addressed and confirmed by our shotgun sequencing analyses below.

In addition, we performed a perturbation analysis by testing whether rhythmic versus disrupted feeding alters the recovery of the microbiome from a brief antibiotic treatment. This perturbation analysis was accomplished by treating approximately half of the mice in each group (n = 11 for WT, n = 13 for Per1/2-dko) for the first 5 days of RR with an antibiotic cocktail that dramatically knocks down both gram+ and gram- bacteria (vancomycin, streptomycin, and metronidazole [32,33]), hereinafter referred to as the WT(T) and Per(T) groups. This knockdown is clearly indicated by the decrease of total bacterial load in response to the antibiotic treatment (S1E and S1F Fig). The 16S sequence analyses confirmed that the antibiotic treatment dramatically changed the relative proportions among the bacterial phyla (Fig 3A and 3D). For example, in both WT(T) and Per(T) samples, Bacteroidetes bacteria were strikingly diminished by antibiotics and thereafter gradually recovered over the next weeks in RR. Interestingly, Verrucomicrobia bacteria became a significantly larger proportion of the microbiome in the WT(T) samples after the antibiotic knockdown, but this effect was much less prominent in the Per(T) samples. Tables that detail the specific genera and species that respond to the antibiotic treatment are included in the Supporting information (S3 Table for genera and S4 Table for species). As in the case for the UT samples (Fig 2), we then applied Shannon and Chao1 alpha-diversity indices [28–30] to assess microbiotal diversity, evenness, and richness of OTUs. These indices highlighted the dramatic perturbation by the antibiotics in the treated samples (Fig 3B, 3C, 3E and 3F). By these alpha-diversity analyses, both the WT(T) and the Per(T) samples appeared to recover over 6 months from the antibiotic knockdown during the RR treatment. Changes in the taxonomic composition will be addressed by our shotgun sequencing analyses below.

The complementary beta-diversity index [29,31] was enlightening towards visualizing the differences between the WT(T) versus Per(T) samples as they recovered from the brief perturbation. Fig 3G and 3H shows 5 key time points in which arrows track the recovery trajectories (and all 16 time points are animated in S1 Video and graphed in S3 Fig). Just before the antibiotic/RR treatment (e.g., Day 0), the beta-diversity of both WT and Per1/2-dko mice is essentially the same (the WT(T) samples are somewhat more distributed in the PCoA1 dimension than the Per(T) samples), and their response to the antibiotic knockdown is comparable (Day 11, which is 6 days after the antibiotic treatment ended). Remarkably, the microbiome diversity of WT mice that feed with approximately 24-h rhythms in RR (Fig 1) largely recovers from the antibiotic perturbation by Day 182 in RR, and finally recovers to its original PCoA values in the final LD {WT(T), Day 238 in Figs 3G and S3}. In contrast, the beta-diversity of the microbiome from the temporally disrupted Per1/2-dko environment never entirely recovers to its initial starting diversity; even 8 further weeks in LD (microbiomes in now rhythmically behaving hosts) allowed minimal further progress towards recovery {Per(T), Day 238 in Figs 3H and S3}. Therefore, even though the alpha-diversity metrics suggested that both WT(T) and Per(T) microbiomes recover over 6 months from the antibiotic knockdown during the RR treatment, the beta-diversity metrics reveal that this “recovery” is associated with major compositional changes in the communities cultivated in the Per(T) environment.

It is as though the disrupted feeding patterns lock the microbiotal composition into a new steady state that remains largely intransient. There appear to be no major differences between the 2 mouse strains in terms of how their microbiotal diversity responds to LD/RR transitions (Fig 2G and 2H) or initially to antibiotic perturbations (Fig 3B, 3C and 3E–3H, Day 11; also see S3 and S4 Tables). However, in the mice treated with antimicrobials, the longitudinal recovery patterns of the microbiota differ significantly between mice that consistently maintain rhythmicity {WT(T)} versus temporally disrupted behavior {Per(T)}. Here, the abundance and presence/absence of microbial species in Per(T) mice failed to recover to its initial state following the return to LD (PERMDISP: Bray–Curtis: P = 2.38 × 10⁻¹⁰; Jaccard: P = 3.08 × 10⁻⁹), whereas the microbiota of WT(T) mice largely recovers in this timeframe (PERMDISP: Bray–Curtis: P = 0.46701; Jaccard P = 0.49837) [29,34–36]. Note that antibiotic treatment not only affects microbial abundance and diversity, but also has short-term effects on the gut epithelium, including the characteristics of circadian rhythmicity in intestinal epithelial cells [19]. However, as the perturbation inflicted in this study is a brief 5-day treatment with antibiotics, and the mucosal epithelial layer is reported to recover from antibiotic treatment within 3 weeks [37], any effects on the epithelium should not significantly affect the microbiome on the long-term time scale that we are analyzing. An unavoidable consequence of our long-term protocol is that the mice are aging for approximately 9 months during the procedure and this might have influenced the recovery from antibiotics, but this factor is also true for the untreated groups that exhibit minimal changes over time (Fig 2) and therefore aging appears to not be a major factor by itself.

In the WT(T) samples, we observed an interesting bifurcation in the kinetics of the return in RR to the initial PCoA coordinates. As seen in S4A Fig, at Day 42 the 4 mice in one of the group-housed cages (Cage#4) have dissociated from the main group and are practically back to the initial conditions, whereas it takes until Day 182 for the singly housed mice and the other group-housed cage to eventually return to the PCoA region of the initial LD samples (note that the microbiome of the singly housed mouse in Cage#1 was unexpectedly not influenced by the antibiotics in the drinking water). This bimodality of return is consistent with the conclusion that the microbiomes of animals that are housed together tend to be similar [38,39], and that the microbiomes of 1 or more mice in Cage#4 reverted to initial LD conditions more quickly and influenced the more rapid restoration of the microbiome diversity of the other mice in that cage. Our experimental design of including both types of housing therefore confirms that our perturbation analysis yields the same final results independently of the housing conditions. Interestingly, no such bifurcation was observed in the Per(T) samples—both singly and group-housed mice recovered partially (but only partially) towards the initial conditions in coherent distributions (S4B Fig).

To further refine our analyses to changes in the gut microbiome at the species/genera levels, we undertook shotgun sequencing of samples at 4 key time points that we identified as pivotal nodes in the beta-diversity data, namely Days −14, 11, 154, and 238. These time points correspond to initial conditions (Day −14), 11 days after transfer to RR (and antibiotics for treated mice, Day 11), 154 days after transfer to RR (Day 154), and 55 days after transfer back to LD (Day 238, see amber marks in Fig 1A). As expected, there was a large variety of differential responses among the identified species/genera to the changing environmental conditions (S5 Table). The data were analyzed at both genus and species levels by 4 different methods that have been applied to microbiome data and the results were consistent among all tests: (i) non-parametric Mann–Whitney U test; (ii) linear regression; (iii) linear mixed-effects model; and (iv) feature volatility [40–47]. We predicted that bacteria that adapt optimally to a temporally cyclic environment—perhaps due to resonance with an internal timekeeper [5]—might be differentially affected by arhythmic conditions as compared with rhythmic conditions. In the untreated (UT) mice, such a hypothetical bacterium would be predicted to be unaffected by the LD/RR/LD transitions in WT(UT) mice (consistent rhythmic behaviors), but show a disruption by RR in Per(UT) mice (arhythmic in RR) that recovers in LD (Fig 4A). Interestingly, that pattern was observed in the species Corynebacterium casei and Corynebacterium stationis (both species increase in RR in Per(UT) mice; S5 Table). A variant on these predictions was observed in 6 different species of Lactobacillus (L. acidophilus, L. amylovorus, L. acetotolerans, L. amylolyticus, L. helveticus, and L. crispatus), namely that the transfer to RR significantly knocked down the abundance of these Lactobacillus species (Day 11) and they thereafter recovered partially by Day 154 in RR, and recovered more completely in LD on Day 238 (S5 Table).

For the antibiotic-treated samples, Fig 4B depicts a “volcano plot” of the abundance and significance changes at the species level of the WT(T) and Per(T) microbiomes analyzed by Mann–Whitney U test followed by a p-value adjustment for false discovery rate [40,45,46]. These plots show changes in abundance where each datum is for a different species as identified by the shotgun sequencing. For the WT(T) samples, the abundance of many species is knocked down by the antibiotic treatment by Day 11 (755 species are underrepresented, blue points in Fig 4B), but there is considerable recovery in RR (Day 154) and essentially complete recovery by Day 238, as indicated by the absence of significant differences in the upper right panel of Fig 4B, which depicts the species representations on Day 238 (at which time only 27 species are underrepresented). On the other hand, the pattern for the bacteria residing in Per1/2-dko mice is quite different. Soon after antibiotic treatment, there are many bacterial species that are knocked down by the antibiotics and transfer to RR (363 species become underrepresented) as well as a few that become overrepresented (54 species, red data points in Fig 4B), as shown in the lower left panel (Day 11) of Fig 4B. At Day 154 in RR, some species remain repressed (blue data points), but there are species that were previously (Day 11) not knocked down by the antibiotics that now have become overrepresented (black data points in the upper right quadrant of the middle bottom panel of Fig 4B). This effect appears to be even more obvious on Day 238 (lower right panel of Fig 4B, where 116 are still underrepresented and 199 are overrepresented). Some of these species may be “opportunists” that have taken advantage of the antibiotic knockdown and the RR conditions to enhance their abundance relative to the total composition. Another way of depicting this volcano plot is to color-code on the basis of the final distribution (S5 Fig); this alternative depiction makes even more clear that most of the microbiota that are overrepresented on Day 238 in the Per(T) samples were not significantly knocked down by the antibiotic treatment and grew to an overrepresented abundance progressively during the recovery period. These volcano plots for the Per(T) samples imply that a new steady state microbial composition has been achieved in LD after the RR conditions that differs from the initial LD conditions, which is another way of visualizing the effects shown in Fig 3 for the Per(T) samples.

An alternative statistical test for longitudinal timecourse samples of microbiomes uses a linear mixed-effects model [42–44], for which the results of 258 most abundant species in the WT(T) and Per(T) samples are visualized in the volcano plot of Fig 4C. The same linear mixed-effects model was applied to genera, and the results of that analysis identified 182 most abundant genera that are depicted in S6 Fig. The overall trends between the WT(T) and Per(T) samples are similar among Figs 4B, 4C, and S6, but in addition, several genera of particular interest are highlighted in S6 Fig. As explained above, we predicted that bacteria that adapt optimally to a temporally cyclic environment—perhaps due to resonance with an internal timekeeper [5]—might recover more slowly in arhythmic conditions after perturbation than in rhythmic conditions (Fig 4A). Such a hypothetical bacterium would show a pattern after perturbation by antibiotics of (i) recovery in RR in WT(T) but more slowly in Per(T) samples (rhythmic environment in 1 situation but not the other), and (ii) recovery in the second LD in both WT(T) and Per(T) samples (rhythmic environments in both situations). Interestingly, the species Bacteroides_sp._PHL_2737, Neobacillus mesonae, Mucilaginibacter rubeus, Porphyromonas gingivalis, and Labilibaculum antarcticum conform to that predicted pattern (S5 Table), but Klebsiella aerogenes that has been suggested to harbor a circadian system [9,10] does not conform to that pattern. Most species of Clostridium seem to be relatively unaffected by all the conditions. Prevotella species are knocked down by antibiotics and generally recover in RR more rapidly in WT(T) than in Per(T), while species of Alistipes are knocked down by antibiotics, but mostly recover in RR in both WT(T) and Per(T).

Akkermansia muciniphila abundance is enhanced after the antibiotic treatment but eventually recovers to initial conditions in the second LD for WT(T) {but does not recover in Per(T)}, suggesting that it develops an “opportunistic” reaction to the knocked-down genera and grows in relative abundance (S6 Fig and S5 Table). Another particularly robust example of an opportunistic pattern that stood out from the other genera after the “volatility test” [41] was the genus Turicibacter. This genus was knocked down by the antibiotics, but recovered to much higher relative abundance in RR for both WT(T) and Per(T) (off the scale for the Per(T)/RR plot in S6 Fig), and was finally suppressed to initial conditions by the second LD in WT(T) but remained dramatically elevated in Per(T) (S5 Table). The mixed model results for species show this same large effect for species of Turicibacter, now plotted on an expanded scale in Fig 4C. One interpretation of these data is that Turicibacter is not strongly suppressed by the antibiotic treatment but exploits the altered gut environment to expand in an opportunistic pattern. This expansion is ultimately contained in the case of the rhythmic WT(T) mice, but is not brought back to preconditions in the temporally disrupted Per(T) mice; the volatility test indicated that Turicibacter is indeed an important contributor to the lack of full recovery of the microbiome of feeding-disrupted Per1/2-dko mice. Interestingly, Turicibacter shows the same effects in the untreated samples in RR, suggesting that Turicibacter is particularly adept in exploiting temporally disrupted gut environments, even in the absence of an antibiotic perturbation (S5 Table). Finally, while Figs 4 and S6 show the results of our Mann–Whitney U and linear mixed-effects model tests, we have also analyzed the data by a linear regression test [47], which yielded equivalent conclusions.

As a complementary approach, we analyzed whether these rhythmic versus disrupted daily patterns of environment have consequences for the functional capacities of the intestinal microbial community based on our shotgun metagenomic sequencing reads. We first filtered the shotgun sequencing reads with Bowtie2 to remove the mouse DNA and then we used HUMAnN3 to analyze the microbial metabolic pathway contributions. The data from Days 11, 154, and 238 are compared relative to the baseline (Day −14) data in S7 Fig, which depicts 24 pathways based on strength of signal and that include a variety of metabolic pathways, including carbohydrate, fatty acid, amino acid biosynthesis, and nucleotide metabolism. These 24 pathways were chosen for inclusion in S7 Fig on the basis of greatest significance and for illustrating (i) changes that occurred between T-samples and (ii) UT samples, and (ii) differences between WT(T) and Per(T) samples. As was the case for microbial abundance and diversity, the microbiota of untreated mice (either WT(UT) or Per(UT) mice) had stable pathway abundances that were not significantly affected by transitions between rhythmic and temporally disrupted conditions. This lack of response in metabolic pathways given the stability of richness and diversity of composition is not trivial, as other studies have found changes in metabolic pathways even when OTU diversity was not significantly affected—for example, in the case of disrupted and phase-shifted sleep/wake cycles in humans [48]. For the microbiomes from the treated mice, many pathways were dramatically down- or up-regulated by the brief 5-day antibiotic perturbation as expected (e.g., data for Day 11) [49], and under continuous rhythmic conditions {WT(T)}, almost all of the pathways had recovered by Day 238. However, for the pathways recovering in the temporally disrupted conditions {Per(T)}, many of these pathways had not returned to initial conditions, even by Day 238 (S7 Fig). Not surprisingly, the antibiotics knocked down bacterial species so that fundamental metabolic pathways such as aerobic respiration and gluconeogenesis were reduced, but these recovered completely only under rhythmic conditions but not disrupted conditions. Amino acid and nucleotide biosynthetic pathways exhibited blunted recovery from antibiotics under disrupted conditions. Sulfate reduction and seleno-amino acid biosynthesis were up-regulated by the brief antibiotic treatment (S7 Fig); increases in sulfur metabolism pathways are associated with antibiotic treatments [48,50]. The microbiotal recovery from this up-regulation is similarly protracted in the temporally disrupted environment (S7 Fig). Together, this analysis of functional pathways reflects the common themes observed with the richness and diversity analyses; the temporally disrupted environments (untreated groups) have limited consequences, whereas after perturbation the stably rhythmic environment significantly facilitates recovery.

All animal experiments were approved by the Vanderbilt University Institutional Animal Care and Use Committee and were conducted according to that committee’s guidelines under Animal Welfare Assurance Number A3227-01 (IACUC protocol # M1600221) that adhere to those published in “The Guide for the Care and Use of Laboratory Animals (8th Edition).”

Adult C57/BL6J background WT or Per1/2 double knockout (“Per1/2-dko” or “Per”) mice were bred to create age-matched cohorts of 3 to 4 weeks of age at the beginning of the experiment. Our Per1/2-dko mice were generated as described by Bae and colleagues [21] and had been introgressed into the C57 background (>8 generations of backcrossing). This mouse strain entrains to LD 12:12 cycles with a consolidated nocturnal activity pattern and exhibits arhythmicity and/or highly disrupted activity patterns in DD [22]. Prior to the study, the parental mice were housed in the same facility and bedding was periodically switched (once every week) between WT and Per1/2-dko mice cages to promote a similar microbiome composition. Bedding exchanges between and among cages of the 2 genotypes were also performed on the experimental mice up until the start of the experiment to homogenize the microbiomes prior to transfer to RR. Mice were provided with a standard lab diet (LabDiet Rodent 5001) ad libitum throughout the study. Because there are sex differences in the gut microbiome [13], we focused upon a single sex to reduce experimental variability, and males were chosen because their circadian wheel-running behavior tends to be more precise because it is not perturbed by the estrous cycle.

A total of 22 WT mice and 28 Per1/2-dko mice were used for the experimental protocol depicted in Fig 1A under both singly and group housed conditions so that their responses could be systematically compared (see S6 Table for housing conditions organized by cage number). Briefly, mice were under an LD 12:12 (12 h light/12 h dark, lights-on at 6:00 AM, cool-white light intensity approximately 400 lux) cycle for 5 weeks (35 days). Mice were then placed under constant dim red light (“RR” approximately 30 lux, see S8 Fig for spectrum of the light source) for 26 weeks (182 days) and then reintroduced to the original LD cycle for the remaining 8 weeks (56 days) of the experiment. Regular chow was provided ad libitum throughout the entire experiment. Fresh fecal samples were collected from all mice weekly during the first LD cycle, twice a week during the first 2 weeks of RR, once every other week for the remainder of the RR treatment, and once a week for the final LD treatment. All fecal collections occurred only between 2:00 PM to 4:00 PM clock time (LDT 8 to 10 for mice under LD conditions). All fecal samples were immediately put on ice after collection and placed in a −80°C freezer for long-term storage. Samples for quality control included: (i) empty tubes that were exposed to the air in the animal care facility during stool collection but which did not receive stool samples; (ii) samples of the chow; and (iii) samples of the water provided to the mice. These quality-control samples were also extracted for DNA and processed along with the stool samples to confirm that bacteria were not introduced to the analysis by artifactual exposure during sample collection.

During the first 5 days of constant dim red light vancomycin (0.5 g/L), streptomycin (2 g/L), and metronidazole (0.75 g/L) were added to the drinking water for the “Treated” mice in both the WT and Per1/2-dko groups (n = 11 for WT, n = 13 for Per1/2-dko) [32,33]. The antibiotic treatment did not cause diarrhea or any other detectable ill effects, but significantly decreased total bacterial load (S1 Fig) and affected many identifiable genera and species (S3 and S4 Tables). The remaining “Untreated” mice had no antibiotic treatment and experienced no significant effect of the LD/RR/LD transitions on total bacterial load (S1 Fig).

Locomotor activity and feeding activity were measured in mice in singly housed cages that underwent the same LD/RR/LD regimen as the mice from whom feces were collected (Fig 1 and S1 and S2 Tables). Both locomotor and feeding activity were monitored simultaneously for 1-week intervals during the second LD 12:12 (LD2) and the RR conditions. These 1-week intervals were dispersed throughout the light condition. Activity of WT and PER1/2 double knockout mice with or without antibiotic treatment was recorded (n = 4 for each mouse strain). Locomotor activity was measured as wheel-running activity and analyzed with Clocklab software by Actimetrics. Feeding activity was recorded using an infrared motion detecting camera (PYLE PLCM22IR Flush Mount Rear View Camera with 0.5 lux night vision, Pyle Audio, Brooklyn, New York, United States of America) directed at the chow feeder and recorded on a DVR (AMCREST 4 Channel Video Security Recorder, AMDV10814-H5) under infrared illumination. The video recordings were assessed for feeding activity by an independent observer who scored the behavior in accordance with the criteria outlined by Pendergast and co-researchers [62]. Specifically, feeding behavior was defined as a mouse reaching into the chow feeder and visibly eating a pellet with its mouth or holding a chow pellet in its paws and chewing for at least 3 s. Feeding behavior was logged as feeding events in 1-min intervals (multiple feeding events within 1 min were still only considered as 1 feeding event), and the data were then collected in 1-h bins. Therefore, the maximum number of feeding events that could be logged in a 1-h bin was 60. To better compare rhythmic activity between treatment groups, feeding and activity data of individual mouse data were then normalized by dividing by the hour of maximum recorded activity followed by smoothing with a 4-point (= 4 h) moving average and plotted as “counts” in Fig 1C. The feeding activity and locomotor activity were aligned by time using the timestamp on the video recording for comparison. There were no obvious differences in the locomotor activity or feeding behavior between the antibiotic-treated and the untreated groups of mice. Moreover, there were no obvious differences in the weight gain of the mice among any of the groups although this was not quantified.

DNA from fecal samples was extracted using a DNeasy PowerSoil Kit from Qiagen (Hilden, Germany). A total of 800 DNA samples (stool samples from the same 50 mice collected at 16 key time points) were processed for 16S rDNA sequencing by the Vanderbilt Vantage Genomic Resource Core. The 16S Amplicon sequencing was performed with Illumina MiSeq (PE 250) on V4 regions [63]. Approximately 200 DNA samples (from the 800 total) at 4 nodal time points (Days −14, 11, 154, and 238) were sent for shotgun sequencing by the Vanderbilt Vantage Genomic Resource Core with the Illumina NovaSeq 6000 (PE150) up to 10 M-reads coverage for each sample. Except as indicated otherwise, for data analyses, we used 5,000 reads per taxa as the cutoff for both 16S amplicon sequencing data (5,000 total reads/sample) and the shotgun sequencing data (5,000 total reads/species or genus) after summing up the total reads per taxa for all 200 samples.

QIIME2 was used for analysis of 16s rRNA sequencing [64] (QIIME2 was run online from https://cgc.sbgenomics.com/↗ with the embedded GreenGenes 16S database). Kraken2 was used for Metagenomic shotgun sequencing [65]. The shotgun sequencing reads were first screened against the mouse reference genome (GRCm38↗) with Bowtie2 to remove the host DNA [66]. Then, HUMAnN3 was applied to analyze the microbial metabolic pathway contributions [67]. Kraken2, Bowtie2, and HUMAnN3 analyses were all performed with the use of the Advance Computing Center for Research and Education (ACCRE) at Vanderbilt University.

Period analyses were performed with the Lomb–Scargle periodogram in Clocklab (Actimetrics, Illinois, USA) and R. A Wilcoxon rank sum test was performed on the power values of WT and Per1/2-dko mice to determine significant changes in power of the period for each condition. The OTU tables obtained from QIIME2 were uploaded to MicrobiomeAnalyst (https://www.microbiomeanalyst.ca/MicrobiomeAnalyst/home.xhtml↗) for statistics [29]. The raw sequencing reads for the 800 samples span from 950 to 47,573 reads. Samples with sequencing depths below 5,000 total reads were excluded from data analyses. For any feature retained, at least 20% of the values contained a minimal of 4 counts. Data were rarefied and scaled through total sum scaling (TSS). Shannon and Chao1 indices were applied to calculate alpha-diversities followed by Kruskal–Wallis tests for significance [28–30]. In order to identify whether the early stages of recovery after antibiotic treatment were different between WT and Per1/2-dko mice, we selected the alpha index data in which we observe the fastest rate of antibiotic recovery (Days 7 to 42) and fit quadratic regressions. This was done using the “stats” R package [68]. After this, we compared the coefficients obtained for the WT and Per1/2-dko regressions [69] in order to obtain the z-score between the 2 means. The p-values for the WT and Per1/2-dko regressions were then calculated through this z-score, and we found that there was no significant difference between the 2 strains, both for the Shannon and the Chao1 alpha indexes (p = 0.22 and p = 0.59, respectively). Bray–Curtis and Jaccard indices were used for beta-diversity calculations [29,31,36] followed by PERMANOVA and PERMDISP tests for significance [34,35]. PCoA was performed on all the samples together, but plotted separately for samples at 5 representative time points in Figs 2 and 3, and then for all samples as a gif video file in S1 Video and as graphs in S2 and S3 Figs. The “feature volatility” analysis with the shotgun data was applied with the q2-longitudinal plugin of QIIME2. The machine learning regressor (random forest by default) was used to learn the data structures and then to identify the important features [41]. For shotgun sequencing data, volcano plots at both genera and species levels were generated by 3 different statistical methods that are commonly used to analyze microbiome data [45] and the results were consistent among all 3 tests: (i) Mann–Whitney U test followed by a p-value adjustment for false discovery rate ([40,45,68], used by ref. [46]); (ii) a linear regression model ([70] used by ref. [47]); and (iii) a linear mixed-effects model (used by refs. [42–44]). In addition to the volcano plots, we used the data analyzed by the Mann–Whitney U test to identify species that conform to the predicted patterns of abundance changes described in Fig 4A. In order to avoid undefined values when calculating log2-fold changes, any zeroes that would end up in the numerator or denominator were changed to be effectively zero (specifically, to 1e-10, a value one order of magnitude smaller than the smallest value we obtained for our relative abundances). The data used for this analysis can be found in S5 Table, where we highlight log2-fold changes that are significant, so that patterns can be more readily identified. For the linear mixed-effects modeling of functional pathways, the samples (mice) were considered as the random effects, whereas the time points and cages were considered to be the fixed effects. The pathway relative abundances of mice on Days 11, 154, and 238 were compared with their values on Day −14 (reference time point) after a natural logarithm transformation and plotted as a heatmap in S7 Fig.

No human experiments, tissues, or data were included in this study.

All animal experiments were approved by the Vanderbilt University Institutional Animal Care and Use Committee and were conducted according to that committee’s guidelines under Animal Welfare Assurance Number A3227-01 (IACUC protocol # M1600221) that adhere to those published in “The Guide for the Care and Use of Laboratory Animals (8th Edition).”

In addition to the data included within the paper and itsfiles, the data generated during and/or analyzed are included inFile. The raw sequencing data are available from the Sequence Read Archive (SRA) at: PRJNA761116. Supporting Information S1 Data