Patients with confirmed diagnoses of IBD who met all the inclusion criteria and none of the exclusion criteria (see Supplementary Information) were recruited from gastroenterology departments, IBD clinics and endoscopy units based at two tertiary referral hospitals in Sydney, Australia, between Oct 2015 and August 2017. Study participants were in complete clinical remission based on their disease activity indices; partial Mayo score for UC <2 (22), Crohn's disease activity index (CDAI) for CD <150 (23) and/or Harvey Bradshaw Index for CD <5 (23–25) confirmed by their gastroenterologists, and supported by endoscopic and histological results, if available. Baseline data were collected, and longitudinal data accumulated monthly. Data comprised scores related to symptoms of psychological state including perceived stress (PSQ) (26, 27), depression- anxiety and stress (DASS) (28, 29), depression in medically ill (DMI) (30, 31), personality characteristics, i.e., negative affectivity (NA) and social inhibition (SI) traits (32), wellbeing scores (33–35) and sleep quality (PSQI) (36) with clinical course and disease activity as a measure of outcome. To build the outcome, all variables including severe disease symptom/s and flare events were considered through formalized follow up assessments by study investigators. We studied the longitudinal dynamics of multiple contributing factors to disease activity from a cohort of 50 IBD patients (CD, n = 26; UC, n = 24; Tables 1, 2). IBD participants were grouped in disease types (CD or UC) and subdivided based on use or not of biological treatments. Monthly blood and stool samples were collected for assessment of serum C-reactive protein (CRP), fecal calprotectin (FC) levels and microbiome dynamics. The microbiome composition in each sample was determined by sequencing the V4 region of the 16S rRNA gene and a total of 3.3 million contigs were retained after quality control and subsampling. To determine the links between the gut microbiome and clinical components, we collected clinical data including inflammatory biomarkers, clinical indices for disease activity, and self-administered validated questionnaires to quantify psychological state on a monthly basis for a period of 12 months while participants still met the inclusion criteria for the study. Data related to routine medications, use of antibiotics and probiotics and disease activity scores were registered at monthly follow up visits. The end point was at 12 months follow-up or at confirmation of the onset of relapse; in the latter case the last assessment reflected the state at the time of disease relapse. The study was approved by the Human Research Ethics Committee of South Eastern Sydney Local Health District (Ref: 15/094 HREC/15/POWH 245−20 Aug).

Fecal samples were self-collected by participants and were aliquoted and stored at −80°C together with original collection pots for DNA extraction and quantifying of FC. The concentration of FC was assessed by commercially available enzyme-linked immunosorbent assay (ELISA; Calprotectin Elisa Buhlmann Laboratories, S100A8 and S100A9), according to the manufacturer's protocol. Twenty ml peripheral blood samples were drawn and centrifuged (at 2,000 g for 15 min). Serum samples were aliquoted and stored at −80°C to be used for CRP quantitation by high sensitivity ELISA (37). For microbiome assessment, stool samples were taken from collections at months 1, 4, 8, and 12; where there were any missing samples, the sample from the preceding available month was used.

Genomic DNA was extracted from 0.11 to 0.12 g of fecal material of each sample using the Allprep Power Viral DNA/RNA Kit (Qiagen) DNA extraction protocol. Briefly, the concentration of the extracted DNAs was quantified by Nanodrop mass spectrophotometry (38) and Qubit 2.0 fluorometer (Invitrogen) according to the manufacturer's instructions before dilution to 20 ng/μl. For DNA concentrations of <5 ng/μl (39 extractions), SYBR-Green-based qPCR assay was performed to quantify the absolute amount of a target sequence, to compare relative amounts of a target sequence between samples and to analyse whether they were amplifiable (Supplementary Figure 1; Supplementary Table 1 qPCR plot and table). DNA was submitted for sequencing at the Ramaciotti Center for Genomics (Australia).

Barcoding PCR for bacterial and archaeal 16S rRNA genes was carried out using a mix of 10 μL of HotMasterMix (5 PRIME), 0.2 μM of each primer and 1 μL of DNA template. Barcoded PCR primers based on 515F, and 806R (39). Reactions were kept at 94°C for 3 min for denaturing, followed by 35 cycles of denaturation at 94°C for 45 s, annealing at 50°C for 1 min and elongation at 72°C for 1 min 30s, ending with a final elongation at 72°C for 10 min and final hold at 4°C. All PCRs were carried out in 25 μL volumes. PCR concentrations were normalized and pooled using SequalPrep^™ Normalization Plate Kit (ThermoFisher) according to the manufacturer's instructions. The library was purified using Axygen^® AxyPrep^™ Mag PCR Clean-Up Kit (Fisher Biotec) as per the manufacturer's instructions. Concentration and quality of the pooled library were checked with Qubit^® and the library size on an Agilent 2,200 TapeStation instrument.

Raw sequencing data were processed with the OTUreporter pipeline (40), based on mothur v1.39.5 (41) and according to the MiSeq SOP. Samples with a length between 228 and 278 bp were retained and those with homopolymers longer than 8 bp were removed. Sequences were grouped into OTUs based on 97% similarity using the OptiClust algorithm. From each patient, quarterly microbiome samples with matched FC concentrations were sub-selected for use in downstream analysis.

Following the power analysis estimation (80%) to detect significant (two sided p-values ≤ 0.05) correlations (42) cluster analysis was used to examine whether categories of respondents (IBD patients) share common characteristics within clusters at the baseline and maintain the same properties over time (43). Utilizing SPSS 25 (SPSS Inc), general linear models, non-parametric methods, time series analysis and mixed model analysis were used to examine longitudinal data related to psychological and biological measures, their linear/quadratic trends and differences between two groups of treatment options over time (44). Furthermore, a regression analysis using panel data was conducted applying Stata (Stata V16; Stata Statistical Software: Release 16. College Station, TX: StataCorp LLC. 2019) to apply corrections for multiple analysis. Microbial profile analysis was carried out using the phyloseq (45) and microbiome (https://github.com/microbiome/microbiome) R packages to import and graph data, and vegan was used to perform differential abundance testing (46–48). Wilcoxon signed-rank test, and Permutational multivariate analysis of variance (PERMANOVA) test, were applied to examine statistical significance of differences in some bacteria abundances at family level and to evaluate compositional differences between the CD and UC groups, and within treatment modalities (49, 50). Bray Curtis dissimilarity matrices were used as input for PERMANOVA to evaluate compositional differences between the CD and UC and treatment modalities. Analysis of Variance (ANOVA) was used onto beta dispersion test output to validate significant result obtained by PERMANOVA. Community and species diversity were estimated using the Shannon diversity index (51) while species richness estimates were generated using Chao1 (52). Pielou's evenness index was used to examine species evenness (53). Categorical variables were used in subgroup analysis (alpha diversity, abundance testing with LEfSe) (54) (detailed information on statistical methods can be found in Supplementary Table 2).

The entry disease remission data related to bio-psychological state of all IBD patients were employed as a benchmark for assessment of their longitudinal dynamics. Sample analysis included 50 stool and 50 blood samples from CD (n = 26) and UC (n = 24) individuals who were in clinical remission. As expected, CD and UC participants had similar baseline distributions of inflammatory biomarkers (Figure 1) which were strongly and positively associated with one another, but none were significantly related to any of the psychological measures and sleep quality at study entry. There was an analogous baseline distribution for most psychological factors, while there were significant differences in baseline scores related to depression (DASS, Dep, p = 0.026) and negative affectivity (DS-NA, p = 0.001) between CD and UC patients at reference point with UC patients showing higher baseline depression and negative affectivity scores compared to CD patients (Supplementary Table 2). Baseline anxiety scores were significantly associated with sleep disturbances (p < 0.001, r² = 0.52) and baseline stress scores were strongly related to sleep duration (p < 0.001, r² = 0.46). Both outcomes suggest strong relationships between anxiety and stress, and sleep — Supplementary Table 2.

Both biologic and non-biologic treatment groups were similar in most baseline psychological scores with the exception of higher baseline stress scores (p = 0.004) and negative affectivity scores (p = 0.004) in non-biological group. At baseline, UC patients had higher mean psychological scores compared to CD patients in both treatment groups (except for DS-SI score which was higher in CD group who received biological treatment). These findings overall were not surprising given the mandated clinical remission at baseline—Supplementary Table 2.

Data showed no significant baseline differences in the Shannon index (p = 0.67), Pielou's evenness (p = 0.80) and Chao1 (p = 0.429) between CD and UC patients. Comparison of the baseline phyla abundance is shown in Figure 2. Results showed a strong negative baseline relationship between sleep latency and Shannon index (p = 0.001, r² = −0.425), Pielou's evenness (p = 0.002, r² = −0.401) and Chao1 index (p < 0.001, r² = −0.455), indicating that lower intestinal microbial diversity, richness, and uneven microbial composition was associated with longer time to fall asleep.

Detailed comparison of the baseline phyla abundance showed that Verrucomicrobia was present in 23/25 CD participants (4.0% ± 8.2), but only in 17/24 UC participants (2.3% ± 5.8). CD participants had more baseline abundance of some disease associated bacteria when compared to UC including Bacteroidetes (CD = 8.1% ± 9.8, UC = 6.2% ± 7.9) and Proteobacteria (CD = 4.2% ± 7.4, UC = 2.6% ± 6.2), but UC participants had more Firmicutes compared to CD at baseline assessment (CD = 75.7% ± 13.5, UC = 77.4% ± 11.8). Notwithstanding these data, none of the baseline differences were statistically significant between the CD and UC (Wilcoxon test, Figure 2; Table 3).

We used cluster analysis to study the bio-psychological behavior of IBD participants over time. Here clusters were identified based on systematic relationships found in psycho-biological and microbiome dynamics over time and across all study participants. Three stable clusters were identified which persisted in their baseline group categories during follow-up. One potential explanation is most participants remained in clinical remission during the 12 months follow up period (Table 4A). The first cluster (No 1) was the youngest IBD cluster which earned the highest anxiety, depression, and stress scores over time and revealed worst sleep quality sores, higher measures related to depression in medically ill state during follow up (Supplementary Tables 3, 4). This cluster also had lowest number of IBD patients on biological therapy. Second cluster (No. 2) is the oldest cluster with higher number of participants in this group (both male and female and mostly UC, Table 4B—detailed information on cluster analysis is included in Supplementary Tables 20, 21).

Longitudinal bio-psychological data of all remission IBD patients (CD, n = 25; UC, n = 17) were analyzed (460 blood and 460 stool samples) using mixed models. Results did not show significant shift in measures of inflammatory biomarkers within IBD cohort which was expected due to persisting remission. Results identified significant longitudinal coefficient of change in psychological scores (linear, quadratic or both; Supplementary Table 2) including negative affectivity showing greater magnitude of change in the remission UC group and sleep quality showing larger magnitude of change in the remission CD group. Microbial diversity and richness displayed larger magnitude of linear coefficient of change in the remission UC group. Further examination of complete study cohort including independent variables (psychobiological factors) with outcome variables (wellbeing scores and inflammatory biomarkers) was applied. In CD patients, results suggested a statistically significant negative relationship between wellbeing scores with depressive scores (p = 0.018) and positive relation with sleep quality (p = 0.027) over time. Outcomes did not suggest any significant interdependence between longitudinal psychological scores, sleep, mode of treatment and microbial indices in CD patients. In UC patients, longitudinal wellbeing scores retained a significant and positive relationship with sleep quality (p = 0.040), and a significant and unexpectedly negative relationship with stress (DASS stress, p = 0.009). FC had a positive and significant association with sleep quality (PSQI, p = 0.006) and a strongly negative relationship with stress (DASS stress, p = 0.023) (Supplementary Table 17). Assessment of longitudinal psycho-microbial dynamics in UC cohort did not suggest any significant interdependence between longitudinal psychological scores, sleep, and mode of treatment with microbial dynamics in UC patients except for a strong negative association between microbial diversity and depression (DASS depression, p = 0.011), and a negative association between microbial evenness and depression (p < 0.001, Figure 3) but not with microbial richness, Supplementary Table 19.

Analysis demonstrated a significant interplay between wellbeing scores with Shannon index and with Pielou's evenness but not with Chao1 (Supplementary Table 7). To examine microbial composition shift over time, microbial dynamics were explored in remission samples of IBD subtypes in addition to changes in FC concentration longitudinally. CD cohort in clinical remission showed greater microbiome fluctuations mainly by trading off different microbial families as well as in relative abundance of existing microbial profile (Figures 4, 5).

Longitudinal trends of psycho-biological factors and their interactions with modes of treatment (biological vs. no biological) were examined using mixed models of analysis. In CD and UC cohorts such factors were similarly distributed in both treatment groups over time (Supplementary Table 2). At the baseline there were similar ecological indices between the two treatment options in both disease classes which was also suggested by previous study (55) (Shannon, p = 0.380; Pielou's, p = 0.246; Chao1, p = 0.934), although the age effect was significant between the two groups for all three ecological indices (56) (information related to sample demographics for microbial analysis and output of PERMANOVA for both biologics and bio-flare variables in CD and UC cohorts are in Supplementary Table 9).

At the family level (Figure 6A), the CD group showed significant differences (Wilcoxon, p ≤ 0.05) between the treatment received (biologics vs. non-biologic) in the abundance of Barnesiellaceae, Bifidobacteriaceae, unclassified Clostridia, and Clostridiaceae (Figure 7). Results showed similar ecological index values in samples from the two treatment options. Linear discriminant analysis (LDA) with LEfSe was performed to characterize the differences between two treatment options in CD group (both in remission and during relapse). A total of 12 microbial biomarkers (OTUs) characteristic of CD under biological treatment, and 22 in CD with non-biological treatments, were identified (Figure 8).

Samples in UC groups on biological and non-biologic treatment modes (Figure 6B) were found to have significantly different microbial communities based on Bray-Curtis dissimilarity (PERMANOVA, p = 0.0031). Microbial family composition in UC group and on two treatment modalities showed a total of 12 families that were significantly different (Wilcoxon, p ≤ 0.05): Acidaminococcaceae, unclassified Bacteroidales, Barnesiellaceae, Christensenellaceae, Clostridiales vadinBB60 group, Coriobacteriaceae, Defluvitaleaceae, Eggerthelaceae, Fusobacteriaceae, Muribaculaceae, Prevotellaceae, and Streptococcaceae (Figure 9). Results identified significant differences in microbial diversity (Shannon index, p = 0.041), and evenness (Pielou's evenness, p = 0.045) between the two treatment modes (Figure 10) (see Supplementary Table 22 for microbial dynamics across the three clusters). No differences were identified in richness (Chao1) in UC participants (Figure 10) (57). LDA of UC group between two treatment options identified 11 microbial biomarkers (OTUs) which were more abundant in UC group with no-biological treatments and 35 microbial biomarkers were identified in UC group under biological treatments (Figure 11).

This study was the first prospective longitudinal study designed to evaluate bio-psychological interdependence, their complex orchestrated interplay, and their influence on the clinical course of IBDs including the treatment modalities. Time series analysis was used to examine bidirectional and longitudinal interplay between multiple disease contributing factors to measure the rebound effect. A limitation of this study was lack of healthy controls. Diet is a very significant factor affecting gut health and microbial profile (95, 96), but dietary intake was not assessed in this study. None of the participants were exmained for potential pre-existing psychological conditions. Clearly, replication of the current study is needed to further test the conceptual framework with a larger sample size, a better sample representation and by including a large cohort of healthy controls.