Dysbiosis in the Gut Microbiota in Patients with Inflammatory Bowel Disease during Remission.

Inflammatory bowel disease (IBD) is a chronic, relapsing, inflammatory disorder which comprises two main conditions: Crohn's disease (CD) and ulcerative colitis (UC). Although the etiology of IBD has not been fully elucidated, the gut microbiota is hypothesized to play a vital role in its development. The aim of this cross-sectional study was to characterize the fecal microbiota in CD or UC patients in a state of remission to reveal potential factors sustaining residual levels of inflammation and triggering disease relapses. Ninety-eight IBD patients in a state of clinical remission (66 UC, 32 CD) and 97 controls were recruited, and stool samples, as well as detailed patient data, were collected. After DNA extraction, the variable regions V1 and V2 of the 16S rRNA gene were amplified and sequenced. Patients with IBD had a decrease in alpha diversity compared to that of healthy controls, and the beta diversity indices showed dissimilarity between the cohorts. Healthy controls were associated with the beneficial organisms unclassified Akkermansia species (Akkermansia uncl.), Oscillibacter uncl., and Coprococcus uncl., while flavonoid-degrading bacteria were associated with IBD. Network analysis identified highly central and influential disease markers and a strongly correlated network module of Enterobacteriaceae which was associated with IBD and could act as drivers for residual inflammatory processes sustaining and triggering IBD, even in a state of low disease activity. The microbiota in IBD patients is significantly different from that of healthy controls, even in a state of remission, which implicates the microbiota as an important driver of chronicity in IBD. IMPORTANCE Dysbiosis in inflammatory bowel disease (IBD) has been implicated as a causal or contributory factor to the pathogenesis of the disease. This study, done on patients in remission while accounting for various confounding factors, shows significant community differences and altered community dynamics, even after acute inflammation has subsided. A cluster of Enterobacteriaceae was linked with Crohn's disease, suggesting that this cluster, which contains members known to disrupt colonization resistance and form biofilms, persists during quiescence and can lead to chronic inflammation. Flavonoid-degrading bacteria were also associated with IBD, raising the possibility that modification of dietary flavonoids might induce and maintain remission in IBD.

INTRODUCTION

The gut microbiota has been implicated to have a major role in health and disease. With an estimate of 4 × 1013 microbes (1) in the human intestinal tract comprising thousands of different species of bacteria, the abundant and diverse members of the gut microbiota collectively possess a large number of genes that enable them to have critical functions that contribute to the host’s health. These roles include, among others, fermentation of complex carbohydrates to produce short-chain fatty acids (SCFA) (2), synthesis of amino acids (3), and synthesis of certain vitamins (e.g., biotin and phylloquinone) (4). The microbiome is also involved in the development and function of both the local and systemic immune systems with the uncompromised gut microbiota protecting against pathogen infections via a phenomenon called “colonization resistance” (5–7).

The gut microbiome is continually exposed to external perturbations such as infections, antibiotics, and dietary changes. This ecosystem has the ability to respond to perturbations by reaching a limited number of “stable equilibrium states” that are beneficial to both the host and the microbiota (8). A few of these states, however, can disrupt the symbiosis between the host and the bacterial community. Such unfavorable states are termed “dysbiosis.” Dysbiosis can be defined as a “stable” and detrimental disturbance in the complex biological system of the microbiome. In the structure of the microbiome, this can have long-lasting compositional and functional implications for the microbiota and the host and can be an important factor in the pathogenesis or chronicity of certain diseases (8). One such disease attributed to the microbiota is inflammatory bowel disease (IBD). Ulcerative colitis (UC) and Crohn’s disease (CD), the two main types of IBD, are chronic conditions resulting in relapsing and remitting inflammation of the intestine. Patients with IBD can be either in the active phase of their disease or in remission. Remission is characterized by a decrease or absence of symptoms, together with lack of biochemical, endoscopic, or radiological evidence of inflammation. The recurrent and long-lasting nature of IBD has still not been completely elucidated, with studies postulating that environmental factors, genetics, immunological abnormalities, and particularly the intestinal microbiota are contributory to the onset and course of the disease.

Changes in complexity and composition of communities have been reported before in CD and UC (9–11), which can result in decreased abundances of anaerobic producers of SCFAs and increases in oxidative metabolism, both detrimental to the gut mucosa. In addition, amino acid biosynthesis and carbohydrate metabolism are decreased in favor of nutrient uptake by the microbiota (12). Studies have identified several species or clusters to be differentially abundant between controls and patients with IBD in remission (11, 13–15). Thus, this state of dysbiosis might persist during quiescence, despite the patient not having any clinically significant signs of inflammation, and may contribute to the relapsing and chronic nature of IBD.

We set out to compare the fecal bacterial microbiota in adult patients with a diagnosis of IBD (CD or UC) who are in the state of remission with a cohort of healthy controls. Analyzing patients in the quiescent stage of their disease allows us to assess whether significant and stable community differences and altered community dynamics persist, even after acute inflammation has subsided. This study, therefore, provides insight into potential microbial factors which can have a role in the maintenance of chronicity in IBD or which can act as possible triggers to new flare-ups.

RESULTS

Patients above the age of 18 years diagnosed with IBD according to the Copenhagen Diagnostic Criteria for either UC or CD who were in remission were recruited from Mater Dei Hospital, Malta. Ninety-eight (n = 98) patients were included in the study (Table 1, Fig. S1). The mean age was 46.1 years (standard deviation [SD] of 16.59) and 62.9% were male. In total, we analyzed 66 UC patients (67.3%) and 32 CD patients. With regard to smoking status, 64.3% of patients with IBD were nonsmokers (had never actively smoked cigarettes), 10.2% were smokers (still actively smoking), and 25.5% were ex-smokers (had stopped smoking for at least 6 months). Ninety-seven (n = 97) healthy controls were recruited, of which 51.5% were male. The mean age was 44.8 years (SD = 15.88). Gender distribution between the patient cohorts was not significant (P = 0.2622; Fisher’s exact test), while CD patients were slightly younger than UC patients (χ2 = 7.8428, degrees of freedom [DF] = 2, P = 0.01981; Kruskal-Wallis test) and body mass index (BMI) was not different between groups (χ2 = 5.7864, DF = 2, P > 0.050; Kruskal-Wallis test).

TABLE 1

Demographics of patients with UC and CD and of controls

Demographic	UC	CD	Control
Gender
Male	43 (65.2%)	18 (56.3%)	50 (51.5%)
Female	23 (34.8%)	14 (43.7%)	47 (48.5%)
Age
<30	11 (16.7%)	13 (40.6%)	30 (30.9%)
31–60	34 (51.5%)	14 (43.8%)	51 (52.6%)
>61	21 (31.8%)	5 (15.6%)	16 (16.5%)
Smoking status
Nonsmoker	45 (68.2%)	18 (56.3%)	81 (83.5%)
Smoker	4 (6%)	6 (18.8%)	6 (6.2%)
Ex-smoker	17 (25.8%)	8 (25%)	10 (10.3%)
BMI
18.5–24.9	21 (31.8%)	16 (50%)	41 (42.3%)
25–29.9	25 (37.9%)	11 (34.4%)	37 (38.1%)
>30	20 (30.3%)	5 (15.6%)	19 (19.6%)
UC classification
E1: ulcerative proctitis	20 (30.3%)
E2: left sided UC	24 (36.4%)
E3: pancolitis	22 (33.3%)
CD classification
Age at onset
A1: <16 yrs		6 (18.8%)
A2: 17–39 yrs		21 (65.6%)
A3: >40 yrs		5 (15.6%)
Disease behavior
B1: inflammatory		23 (71.9%)
B2/B3: stricturing and/or penetrating		9 (28.1%)
Disease location
L1: terminal ileum		11 (34.4%)
L2: colon		8 (25%)
L3: ileocolon		13 (40.6%)
Medication
5-ASA only	39 (59.1%)	6 (18.8%)
5-ASA + thiopurine	14 (21.2%)	6 (18.8%)
5-ASA + thiopurine + biological	8 (12.1%)	8 (25%)
5-ASA + biological	5 (7.6%)	4 (12.5%)
Thiopurine only	0	2 (6.3%)
Thiopurine + biological	0	6 (18.8%)

Disease-specific differences in community diversity and composition persist in the absence of active inflammation.

Alpha diversity describes the complexity of communities of a single environment by measuring the number of species (richness), their abundance patterns (evenness), and their phylogenetic complexity. The strongest factor influencing the complexity of the fecal microbial communities was the difference in health condition (healthy control, CD, UC). Patients with IBD displayed reduced average species richness and community evenness compared to those of healthy controls (Fig. 1; Table 2). Correlations of diversity and subject age were observed, which differed with respect to the underlying health condition. Simpson diversity appeared relatively stable across ages for healthy individuals and CD patients, while under UC, communities became more diverse with increasing age (Fig. 1B). This points toward a possible age-dependent decrease of colonization resistance or simply continuing succession of the bacterial communities (16, 17). Phylogenetic measures of alpha diversity include nearest taxon index (NTI) and net relatedness index (NRI), with positive values indicating phylogenetic clustering, values close to zero indicating neutral or random community assembly, and negative values indicating phylogenetic overdispersion, either over the whole phylogenetic tree (NRI) or across the closest related species/tips of the phylogenetic tree (NTI). This study showed an overall decrease in phylogenetic clustering and a trend toward random assembly to phylogenetically overdispersed communities with increasing age. However, on average, phylogenetic clustering could be observed, which was strongest in patients suffering from UC and/or CD compared to the healthy controls, hinting toward either loss of species or replacement by a phylogenetically restricted group of bacteria (Fig. 1C and D; Table 1). Thus, CD and UC patients show on average a less diverse and phylogenetically constrained microbial community compared to that of healthy controls.

FIG 1

(A) Alpha diversity via the Chao1 species richness in the different cohorts showing a lower alpha diversity in patients with IBD than in controls. (B) Correlation of Simpson diversity with age for the different cohorts. (C and D) Correlation of the phylogenetic measures net relatedness index (NRI) and nearest taxon index (NTI) with age for the different disease cohorts. The dashed line marks the cutoff between phylogenetic overdispersion and phylogenetic clustering (*, P ≤ 0.050; **, P ≤ 0.010; ***, P ≤ 0.001; #, P ≤ 0.1000). Correlation of community diversity and age within the respective patient groups has been assessed for Simpson diversity (R2Con. = 0.00110, R2UC = 0.09695, R2CD = 0.0266), NRI (R2Con. = 0.04094, R2UC = 0.02767, R2CD = 0.07245), and NTI (R2Con. = 0.03094, R2UC = 0.07355, R2CD = 0.01193).

TABLE 2

Linear model results of the analysis of the different alpha diversity metrics with respect to health condition and age

Index	Model factor	DF	F	P	Comparison	P pairwise
Chao	IBD	2,191	8.63883	0.00026	Contr.-CD	0.0041
	IBD				Contr.-UC	0.0016
	IBD				CD-UC	0.8882
Simpson	Age	1,188	2.24911	0.13537	Age:Contr.-CD	0.3515
	IBD	2,188	3.57301	0.02999	Age:Contr.-UC	0.0258
	Age:IBD	2,188	4.00957	0.01971	Age:CD-UC	0.0284
Net relatedness index	Age	1,190	6.19897	0.01364	Contr.-CD	0.9986
	IBD	2,190	2.90330	0.05728	Contr.-UC	0.0566
	IBD				CD-UC	0.2269
Nearest taxon index	Age	1,190	4.70107	0.03139	Contr.-CD	0.0499
	IBD	2,190	5.59379	0.00436	Contr.-UC	0.0094
	IBD				CD-UC	0.9995

Models were reduced by stepwise model selection, minimizing AIC. DF, degrees of freedom; F, F value (ANOVA model); P, P value. Age/BMI was employed as a covariate, and where appropriate, its interaction with IBD pathology was analyzed. “Age:IBD” refers to differences in slope of the age by diversity relationship between the different pathologies. To describe the patterns in more detail, the data were subsetted to describe all potential pairwise differences in slope or average values in the diversity by age relationships. The notation “Age:Contr.-CD” refers to potential differences in slope between control and CD patients, considering its relationship to alpha diversity.

Analyses were also done within the respective disease cohort to investigate the potential differences in diversity between IBD subtypes, following the Montreal classification for CD and UC (18) as well as the influence of medication (e.g., 5-ASA, thiopurine, anti-tumor necrosis factor [TNF]). Although BMI did not have significant effects on alpha diversity in a global analysis and was removed during model selection procedures, it showed minor effects in interaction with different disease subtypes (Table A1). Overall, subtypes had only minor effects on the diversity in CD and UC patients (Table A1), while medication (namely, 5-ASA) in CD patients increased evenness and decreased phylogenetic clustering with age (NRI and NTI). In UC patients (all of whom where on 5-ASA), evenness and richness of the communities increased with increasing age under thiopurine treatment and decreased without it (Table A2). This suggests that, with an advancing age, certain medications can have different effects on the microbiome depending on a person’s age and, thus, may differentially influence inflammatory processes.

When analyzing beta diversity, significant differences in community composition and structure on the taxonomic (Bray-Curtis dissimilarity and Jaccard distance) and phylogenetic levels (generalized UniFrac distance) were detected with respect to health condition (control, CD, UC; P = 0.00010, P = 0.00010, P = 0.00010), as well as to age (Fig. 2A; Table 3; Appendix Fig. A2A, A3A, A4). Even after correcting for differences in age, BMI, and gender among individuals, these effects were highly significant. In particular, differences between microbial communities of healthy individuals to either disease were comparably high, while the separation of UC communities was consistently highest in pairwise comparisons (see Table 3). Interestingly, both IBD pathologies differed significantly from each other, although less strongly than they differed from the healthy subjects. Gender significantly influenced beta diversity under Bray-Curtis dissimilarity (P = 0.04162) and Jaccard distance (P = 0.01050), while smoking status had no significant effect on overall community composition (Table 3).

FIG 2

(A) Principal coordinate analysis (PCoA) of Bray-Curtis dissimilarity with respect to IBD status, including the significant correlation of community distance with age (Table 2). The tips of the colored lines each characterize a sample of a specific IBD status (blue, orange, and red representing controls, CD, and UC, respectively). The point where the lines converge, shown as a dot, depicts the center of weight of the respective cohort (centroid). The distance between the centers of weight visualizes the average dissimilarity between the IBD pathologies. Arrows point in the direction of increasing values for that variable, so the lower the PCo1, the greater the age. The amount of variation captured by the respective dimension/axis is shown for the first two axes. (B) PCoA to visualize community differences among only CD patients and (C) UC patients (Appendix Table A3). Crosses indicate the centroid for the nominal significant parameters as detected via permutative ANOVA (PERMANOVA). See Fig. A4 for PCoA of controls.

TABLE 3

PERMANOVA results with respect to gender, smoking status, and IBD using different beta diversity metrics focusing on differential abundance of ASVs (Bray-Curtis dissimilarity), differential occurrence of ASVs (Jaccard distance, see Fig. S2), and phylogenetic differences between communities (generalized UniFrac distance, see Fig. S3)

Model or comparison	DF	F value	P value	adj. P value (FDR)	R 2	adj. R2
Bray-Curtis dissimilarity
IBD	2, 192	1.44428	0.00010	0.00054	0.00589	0.00458
Gender	1, 191	1.13670	0.01850	0.04162	0.01098	0.00071
Smoking status	2, 191	1.06048	0.06919	0.12509	0.01490	0.00063
Contr., CD	1, 126	1.31345	0.00010	0.00015	0.01032	0.00246
Contr., UC	1, 160	1.72652	0.00010	0.00015	0.01068	0.00449
CD, UC	1, 96	1.15183	0.01710	0.01710	0.01186	0.00156
Jaccard distance
IBD	2, 192	1.25200	0.00010	0.00045	0.00564	0.00260
Gender	1, 191	1.08846	0.00350	0.01050	0.01028	0.00046
Smoking status	2, 191	0.99202	0.66913	0.70993	0.01294	−0.00008
Contr., CD	1, 126	1.19526	0.00010	0.00015	0.00940	0.00154
Contr., UC	1, 160	1.42725	0.00010	0.00015	0.00884	0.00265
CD, UC	1, 96	1.03809	0.05249	0.05249	0.01070	0.00039
Generalized UniFrac distance
IBD	2, 192	2.07221	0.00010	0.00067	0.00629	0.01099
Gender	1, 191	1.21542	0.06159	0.10985	0.01099	0.00111
Smoking status	2, 191	1.06161	0.21498	0.27640	0.02124	0.00064
Contr., CD	1, 126	1.82581	0.00010	0.00015	0.01428	0.00646
Contr., UC	1, 160	2.77997	0.00010	0.00015	0.01708	0.01093
CD, UC	1, 96	1.22807	0.05329	0.05329	0.01263	0.00235

Contr, control.

Further analyses assessing the influence of various factors on beta diversity among patients with CD and UC showed that age also affected community composition in both CD (P = 0.00639, P = 0.00350) and UC (P = 0.00260) (Table A3). Among individuals suffering from CD, gender was also a factor influencing community composition (P = 0.00710). The location of inflammation (subtype L as per Montreal classification) left only weak traces in the taxonomic (Jaccard) and phylogenetic composition of the microbial communities, while age of disease onset (subtype A) significantly influenced the phylogenetic makeup of the communities in CD (Fig. 1B; Fig. A2B, A3B; Table A3). Among UC patients, gender, most medications, and disease subtype (E) did not significantly influence the microbiota composition, while age and smoking exerted some effect, although the associations were no longer significant after correction for multiple testing (Fig. 2C; Fig. A2C, A3C; Table A3).

Differential abundance and indicator taxa.

Differential abundance analysis was performed to detect differences in taxonomic groups between the different health conditions, correcting for potential influences of age, BMI, and gender. At the phylum level, abundances differed significantly with respect to health condition. Firmicutes and Actinobacteria were identified as being more abundant in UC patients (false-discovery rate P value [PFDR] of 0.00264 and 0.00005, respectively), while Fusobacteria was slightly increased among CD patients (PFDR = 0.05931). Only the phylum Verrucomicrobia was significantly more abundant in healthy controls (PFDR = 7.42962 × 10−11) (Fig. 3A; Table A4). Several of these patterns were also represented at the species/amplicon sequence variant (ASV) level. An ASV is a representative, theoretically true sequence derived from error-corrected next-generation sequencing (NGS) amplicon marker experiments, representing a group of bacteria carrying this sequence (19). Twenty-five ASVs that were differentially abundant between the different health conditions were identified (Fig. 3B; Table A5); 13 were overabundant in healthy individuals (e.g., ASV 883 Faecalibacterium uncl. and ASV 20 Akkermansia muciniphila), while 6 were the least abundant among the controls. Three ASVs showed higher abundances among CD patients (ASV 6 Escherichia/Shigella uncl., ASV 497 Dorea uncl., and ASV-709 Subdoligranulum uncl.), while 3 other ASVs had higher abundances in UC patients (ASV 156-Parabacteroides uncl., ASV 373-Faecalibacterium uncl., and ASV 836-Parasutturella uncl.) than in patients with other health conditions. Indicator species analysis, done as a complementary analysis to determine characteristic bacteria associated with the different health conditions, showed that healthy controls were characterized through significant associations with two butyrate producers: ASV 296 Oscillibacter uncl. (P = 0.00040, P < 0.1000) and the also differentially abundant ASV 420 Coprococcus uncl. (P = 0.00050, P < 0.1000; Fig. 3A; Table A6). Butyrate is vital in maintaining homeostasis within the intestinal epithelium, and a decrease in butyrate-producing bacteria has been linked with IBD (20). Another member of the Escherichia/Shigella cluster (ASV 24) was also detected to be significantly predictive of CD (P = 0.00090, P < 0.1000) (Table A6).

FIG 3

(A) Phylum abundances with respect to IBD status/health condition analyzed via DESeq2. Firmicutes (PFDR = 0.00264) and Actinobacteria (PFDR = 0.00005) are more abundant in UC patients, while Fusobacteria are increased in CD patients (PFDR = 0.05931) and Verrucomicrobia are more abundant in healthy controls (PFDR = 7.42962 × 10−11). (B) Heatmap showing significantly differentially abundant and indicative bacteria (Appendix Tables A5 and A6). *, P ≤ 0.050; **, P ≤ 0.010; ***, P ≤ 0.001; #, P ≤ 0.1000.

Differential abundance analyses with respect to disease subtypes as per Montreal classification (18) were also performed. Within the disease cohorts (CD and UC), a total of 181 ASVs were found to significantly differ in abundance depending on age of onset (A), disease behavior (B), or disease location (L, E) (Table A7; Fig. A5 to A8). Interestingly, several members of the Proteobacteria showed a higher abundance in patients with stricturing phenotype (B2 as per Montreal classification) and/or penetrating disease (B3) in CD Escherichia/Shigella uncl. (ASVs 13, 31, 282, and 422), Klebisella uncl. (ASVs 75 and 101) and Enterobacter uncl. (ASV 219; Table A7, Fig. A6). Previous studies have documented a relative increase in bacterial species belonging to the phylum Proteobacteria, such as Enterobacteriaceae, in patients with IBD (12, 21–23).

Not only did smoking affect beta diversity as described above but associations were also identified with different bacteria, like ASV 227 Anaerostipes uncl., which was more abundant in healthy subjects as well as in nonsmoking individuals, or ASV 373 Faecalibacterium uncl., which showed increased abundance in UC patients and people who actively smoked (Table A8). Stratified by health condition, the influence of three major treatment options, 5-ASA, thiopurine, and infliximab/adalimumab (anti-TNF), on the abundance and occurrence of bacteria was investigated (Table A9). Several bacteria appeared to be influenced in similar ways by the same treatment in both UC and CD, including ASV 1278 Veillonella parvula and ASV 14 Alistipes uncl., which were abundant under the influence of infliximab, or ASV 13-Escherichia/Shigella uncl. and ASV 8-Dialister uncl., which were more abundant under thiopurine treatment. Conversely, some bacteria were influenced in opposite directions by the same treatment, e.g., Collinsella aerofaciens (ASV 60). In CD patients, this ASV was associated with thiopurine and anti-TNF but not 5-ASA. In UC (where all patients were treated with 5-ASA), this bacterium was associated with individuals without thiopurine or anti-TNF (Table A9).

Correlation network analysis.

To investigate the bacterial community as a system of interacting agents whose interactions may change as a cause or consequence of IBD, an ASV coabundance network was constructed and analyzed (Fig. 4). Different bacterial network characteristics were measured to approximate the structural importance of the single bacterial taxa within the network, where the degree of a node refers to the number of connections (number of correlated bacteria) and betweenness measured its importance in connecting other nodes (being in the shortest path between nodes). Eigenvector centrality and PageRank index give importance to ASVs that relate with other important ASVs, with eigenvector centrality emphasizing highly connected groups (24–26). On average, network members which were indicators for IBD or differentially abundant among health conditions showed an increased centrality/importance within the correlation networks (Fig. 4B), suggesting that they play more important roles than other bacteria in the community.

FIG 4

(A) Species level ASV correlation network based on SparCC. (B) Comparison of different average node centralities (degree, betweenness, PageRank, and eigenvector centrality) between IBD indicators and nonindicators tested via permutative Wilcoxon test (PFDR ≤ 0.05). (C) Community modules as detected via the fast-greedy algorithm and (D) the association of these community modules to disease are demonstrated by showing the average completeness of the respective module in the respective cohort (*, P ≤ 0.050; **, P ≤ 0.010; ***, P ≤ 0.001; #, P ≤ 0.1000; Appendix Tables A10 and A11).

Directly testing the observed importance of single network members showed several indicators and differentially abundant species to be more important than expected by chance. These included ASV 87 Flavonifractor plautii (formerly known as Clostridium orbiscindens or Eubacterium plautii [27, 28]), a flavonoid-degrading bacteria which was associated with IBD patients more than with controls. Flavonifractor plautii showed a high network importance on several levels, with a very central position (betweenness) and many associations with other community members (node degree), which were themselves important mediators (PageRank). Intriguingly, this bacterium was identified before as being more abundant among CD patients than among controls (Table A10). Eggerthella lenta (ASV 648), another flavonoid-degrading bacterium, also appeared to be highly important in the network and associated with patients suffering from UC (Fig. A9; Table A10). Flavonoids, which are plant metabolites, are known to have multiple beneficial properties which include anti-inflammatory effects (29). In addition, several bacteria which occupied highly central positions in the community network were associated with controls, including Oscillibacter uncl. (ASV 296 and ASV 90), Alistipes shahii (ASV 16), and members of the order Clostridiales (ASV 321 and ASV 749).

Through the analysis of compositional differences of modules between health conditions, closely interacting bacterial groups, which structurally changed due to disease, were identified (Table A11). Four modules (including module 8) were influenced by IBD, which themselves also included several indicator taxa for the different disease conditions. A strongly correlated group of Proteobacteria ASVs was identified by having a very high eigenvector centrality and forming a clearly defined network module (module 8; fast-greedy algorithm [30], Table A11; Fig. 4C). This module was characterized by almost exclusive interactions among Enterobacteriaceae genera (Klebsiella, Escherichia/Shigella, Enterobacter, Citrobacter; 18/23) and some less dominant taxonomic groups. This particular network module further displayed a significantly low completeness in UC (Table A11; Fig. 4D), and its completeness was positively associated with SCCAI (simple clinical colitis activity index) (31) in UC patients and negatively associated with CDAI (Crohn’s disease activity index) (32) in CD patients (Table A11; Fig. A10), adding more evidence that this cluster of Enterobacteriaceae may play a role in the pathogenesis of IBD (33). Module completeness refers to how “complete” or “whole” a group of apparently associated bacteria is, as derived from the correlation network (i.e., module). The analyses are based on evaluating how complete the membership of these bacterial groups (ranging from 0, no group members present, to 1, all group members present) is within each individual subject. The resulting measure was then subjected to further analyses with regard to its association to patient characteristics (e.g., IBD condition, age). The Enterobacteriaceae have previously been linked to intestinal inflammation via production of lipopolysaccharide that can trigger TLR4 signaling (34). Thus, depending on the type of disease, this bacterial cluster may be able to influence disease severity. Other network modules showed patterns of reduced completeness in the cases of pathology and were also enriched in bacteria associated with changes in health condition (Table A12). This translated to specific subgroups within the microbiome which either drive changes in relation to inflammation or are particularly sensitive to a changing gut environment due to inflammation.

DISCUSSION

Alpha and beta diversity differ between patients with IBD in remission and healthy controls.

In this cohort of Maltese IBD patients in remission, significant differences were observed compared to healthy controls. A decrease in alpha diversity in IBD patients is in accordance with other studies (9, 35, 36). This reduction in diversity suggests that the microbiota has less functional capability and redundancy, making it more susceptible to detrimental effects of external perturbations (8). While alpha diversity, as measured by the Simpson index, was relatively stable with increasing age in the controls and CD patients, it increased with age in patients with UC. This points toward a possible age-dependent decrease of colonization resistance or simply continuing succession of the bacterial communities (16, 17). Significant differences were notable in the community composition (beta diversity) between the different cohorts, showing that the communities varied between disease conditions and among disease subtypes even in the absence of active inflammation. Age also influenced beta diversity, suggesting changes in composition with time due to natural progression of communities with host age as observed elsewhere (37). So far, most studies have been conducted with individuals in a state of active disease, which therefore does not reflect the condition of quiescence found between flare-ups as was demonstrated in this study.

Differential abundance analysis shows community differences in an inactive disease state.

Certain ASVs demonstrated higher abundances among IBD patients, like. Dorea uncl. (ASV 497) and Subdoligranulum uncl. (ASV-709), which were more abundant in CD patients. Studies have reported a decrease in these species in patients with IBD; however, either these patients were newly diagnosed and thus in a state of active disease or no information was available on the disease activity of these patients (12, 38). In UC, the Parabacteroides uncl. (ASV 156) was more abundant than in the other health conditions, as shown by Rajilić-Stojanović (11). Parasutturella uncl. (ASV 836), a member of the Proteobacteria, was also noted to be more abundant in UC, and this tallies with another study (39). Members of Verrucomicrobia, including Akkermansia muciniphila (ASV 20), were more abundant in controls. A. muciniphila degrades mucin, producing acetate and propionate (40). Levels of A. muciniphilia are decreased in obesity and metabolic disorders (41), while administration of this bacterium to mice resulted in partial protection against diet-induced obesity and was beneficial to the gut mucosal barrier (42). Akkermansia spp. have been documented to be decreased in IBD patients in remission (11) and also in active disease (43, 44). In mouse studies, A. muciniphila improved dextran sulfate sodium-induced UC (45), though another study showed an increase in the levels of A. muciniphila in mice with dextran sulfate sodium-induced colitis (46). The decreased abundance in our patients with IBD compared to that in our healthy controls supports the notion that A. muciniphilia or its components could be used as a therapeutic agent (41).

The observed differences in bacterial abundances show that although patients are in the quiescent phase of their disease (lack of symptoms, normal fecal calprotectin, normal endoscopy, normal imaging), a high degree of disease-related community differences remain. This pattern is well reflected in the significant differences in community composition between the different disease cohorts in this study. A study assessing whether the presence of certain microbes is associated with a greater risk of relapse when anti-TNF treatment is stopped found that lower levels of Bacteroides spp., Faecalibacterium prausnitzii, and Clostridium coccoides result in a higher risk of relapse (47). This suggests that further assessment of the microbiome during remission can influence treatment decisions.

Flavonoid-degrading bacteria have a central role in a correlation network.

Two flavonoid-degrading bacteria, Flavonifractor plautii (ASV 87) and Eggerthella lenta (ASV 648), showed a high and consistent network importance on several levels and are also strongly associated with CD and UC, respectively. Flavonoids are secondary plant metabolites with antioxidant properties, which have been shown to suppress inflammation in both in vivo and in vitro studies, decreasing the severity of various inflammatory diseases, including IBD (29, 48, 49). Furthermore, many flavonoids have been described to possess antibacterial, antifungal, and even antiviral properties (50). Flavonoids have been studied as a potential therapeutic option for IBD through various mechanisms, including preservation of the epithelial barrier, various immunomodulatory functions influencing the gut microbiota, and modulation of enterohormone secretion (51). Though not commonly recognized, Flavonifractor plautii was shown to be more abundant in controls than in IBD patients in remission (11) or cases of colorectal cancer (52). Thus, the role of flavonoid-degrading bacteria in IBD remains to be elucidated, but one can hypothesize that the increased breakdown or modification of dietary flavonoids might add to chronicity of IBD by reducing the antibacterial and anti-inflammatory effects of these dietary compounds (29, 50, 52).

Enterobacteriaceae cluster is linked to CD.

Differential abundance analyses revealed that patients with the stricturing and/or penetrating subtype of CD (higher severity) had a higher abundance of Escherichia/Shigella uncl. (ASVs 13, 31, 282, and 422), Klebsiella uncl. (ASVs 75 and 101), and Enterobacter uncl. (ASV 219) than patients with nonstricturing, nonpenetrating phenotype. In addition, ASV 6 Escherichia/Shigella was significantly reduced among healthy individuals. Network analyses revealed a strongly correlated group of Proteobacteria ASVs (Klebsiella, Escherichia/Shigella, Enterobacter, Citrobacter). Some of the ASVs in this group were strong indicators of CD and also showed repeatedly high network importance in other measures. Various studies have linked Enterobacteriaceae with inflammatory processes. Taxa in the Enterobacteriaceae have been shown to disrupt colonization resistance and form biofilms (7), which might be reflected in the tightly connected (strongly coabundant) module 8. Biofilm formation may be beneficial during inflammation and ecological disturbance and increases persistence in the gut environment (53), albeit at the expense of potential competitive abilities (54). Adherent-invasive Escherichia coli has been isolated from ileal mucosa of patients with CD and has been shown to disrupt the gut mucosal barrier by attachment to the epithelial cells or mucosal tissue invasion (33). Klebsiella spp. have been previously associated with IBD, especially in active CD (36), possibly by a TH1 cell induction capability (55) as evidenced by research on the microbiota, immunological investigations, and mouse studies (56). Klebsiella spp. have also been associated with primary sclerosing cholangitis (57), which in turn is associated with IBD (mainly UC). Citrobacter has been documented to trigger inflammation by induction of Th17 cells (6) and was also found to elicit colitis in mice (58). The persistence of the Enterobacteriaceae can therefore be part of the explanation for the chronicity of IBD, as the completeness of this particular bacterial group also correlates with measures of severity.

In conclusion, a high degree of dysbiosis remains in IBD patients in remission, further supporting the notion of long-lasting dysbiotic patterns supporting and maintaining inflammation and chronicity in IBD. The persistent decrease in alpha diversity and clear differences in community composition in patients with IBD may leave these individuals in a state more vulnerable to external and internal perturbations despite being in remission. In view of the association of Enterobacteriaceae and flavonoid-degrading bacteria to IBD, one could hypothesize that these bacteria, some of which are known pathobionts, can flourish during a perturbation, resulting in a relapse of active inflammation. Future studies might indicate whether abundance of certain pathobionts is associated with a shorter time to relapse and might influence decisions pertaining to management of IBD.

MATERIALS AND METHODS

Patient recruitment.

Ethical approval was obtained from the University of Malta Research Ethics Committee (Ref 32/2017) and from Mater Dei Hospital. Maltese patients above the age of 18 years diagnosed with IBD according to the Copenhagen Diagnostic Criteria for either UC (59) or CD (60) were recruited from the Gastroenterology Outpatients department at Mater Dei Hospital, Malta between September 2018 and September 2019. Patients had to be in the quiescent phase of the disease. Quiescence was defined as a fecal calprotectin (FCP) value below 250 μg/g (61) and/or normal recent (<6 months) endoscopic examination and/or normal recent (<6 months) small bowel imaging. All FCP tests were carried out at the Bioscientia Institute for Medical Diagnostics, while all endoscopies and small bowel imaging were carried out at Mater Dei Hospital.

Exclusion criteria were (i) usage of antibiotics, probiotics, bowel-cleansing medications, or systemic corticosteroids within the previous 3 months, (ii) acute illness, chronic autoimmune or inflammatory medical conditions, a history of malignancy, previous gastrointestinal surgery, or current other active gastrointestinal disease, and (iii) females who were pregnant or lactating.

Control samples were randomly collected from a healthy cohort. Controls had the same inclusion and exclusion criteria, with the exception that they were not known to suffer from IBD.

DNA extraction.

The microbiota analysis was carried out at the Institute of Clinical Molecular Biology in Kiel, Germany. DNA was extracted and purified using the QIAamp fast DNA stool minikit automated on the QIAcube (Qiagen, Hilden, Germany) as described by Trautmann et al. (62).

Bacterial 16S rRNA gene sequencing.

Variable regions V1 and V2 of the 16S rRNA gene were amplified using the primer pair 27F-338R in a dual-barcoding approach according to Caporaso et al. (63). DNA was diluted 1:10 prior to PCR, and PCR products were verified using the electrophoresis in agarose gel and normalized using the SequalPrep normalization plate kit (Thermo Fischer Scientific, Waltham, MA, USA). Subsequently, PCR products were pooled equimolarly and sequenced on the Illumina MiSeq v3 2 by 300 bp (Illumina Inc., San Diego, CA, USA). Demultiplexing after sequencing was based on 0 mismatches in the barcode sequences bcl2fastq.

Quality control, classification, and binning of sequences.

Data processing was performed using the DADA2 version 1.10 (64) workflow for big data sets resulting in abundance tables of amplicon sequence variants (ASVs). All sequencing runs were handled separately for error correction, read merging, and combined chimera detection. ASVs underwent taxonomic annotation using the naive Bayesian classifier implemented in DADA2 using the RDP database (v16) as a taxonomic reference and, if possible, classified sequences at the species level (65, 66). ASV sequences were aligned via NAST-alignment to the SILVA core database and filtered for informative sites (constant gaps, constant bases) in mothur (67–70). Phylogenetic tree construction on ASV alignment generated was carried out using FastTree 2.1 using the CAT substitution model with Γ-correction and improved accuracy, employing more minimum evolution rounds for initial tree search (-spr 4), more exhaustive tree search (-mlacc 2), and a slower initial tree search (-slownni) (71).

Statistical methods for microbiome analyses.

Microbiome data were rarefied to 11,800 reads/sample to ensure comparable and sufficient coverage across samples (average Good’s coverage of 99.87% ± 0.001 SD). Species diversity indices were calculated in vegan for R 3.5.1 (72), analyzed via multivariate linear models, and simplified via backwards selection. Phylogenetic measures of alpha diversity (NTI and NRI) were derived using the picante package, based on 999 permutations against a null model preserving relative species richness within the communities (73, 74). Multivariate analyses of alpha diversity measures were carried out with linear models using BMI and age as covariates in interaction with IBD and using model selection to derive the most optimal statistical model (minimizing Akaike information criterium [AIC] criterion, likelihood ratio tests) in R v3.5.3. For beta diversity analyses, taxon-based beta diversity metrics were calculated in the vegan package in R v3.5.3 (72, 75, 76). As measure of phylogenetic differentiation between bacterial communities, the generalized UniFrac measure (d0.5) was utilized (77, 78). Distances were ordinated via principal coordinate analysis (PCoA) avoiding negative eigenvalues by square root transformation. Testing for community separation between the different health conditions or other patient characteristics was performed by distance-based (conditional) redundancy analyses and permutative analysis of variance (PERMANOVA; 10,000 permutations) (79–81) including evaluation of fit via R2 and adjusted R2 (82). The fit of centroids and/or biplots of continuous variables was assessed via a goodness-of-fit test (10,000 permutations), based on the results of the PERMANOVAs.

Indicator species analysis was based on 10,000 permutations utilizing the generalized indicator value (IndVal.g) as implemented in the indicspecies package for R. Differential abundances were tested via general linear models (GLMs) using a negative binomial error structure as implemented in DESeq2 using bacteria present in at least 5% of samples (83), correcting for the covariates’ genders and ages in the analysis (model: bacteria ~ gender + age + IBD + error), and correcting for the covariates gender and age in the analysis (model: bacteria ~ gender + age + IBD). P values of differential abundance analyses and indicator analyses were adjusted via the Benjamini-Hochberg procedure (84).

To generate coabundance networks, the SparCC algorithm (85) as implemented in mothur was employed (50 samplings, 100 iterations, 10,000 permutations, P value adjustment via Benjamini-Hochberg procedure) based on ASV abundances (shared among 10% samples) at a P value cutoff of ≤0.05 and a correlation strength above the median of nonzero correlations. Network modules were detected by the fast-greedy algorithm (30), and node-based importance/centrality measures weighing different node characteristics were calculated as implemented in igraph 1.2.4.1 on the weighted networks (|R| values). (24–26). To assess whether bacteria are more important than expected by chance, bacterial vertices were permuted 10,000 times and their centrality within the random network was calculated. These individual statistics of random communities served as a baseline random distribution to compare the observed centralities against via a one-sided Z-test to determine significance. Module differentiation between health conditions was tested via redundancy analysis of Hellinger transformed community subsets based on module membership, comparable to so-called “eigengene analysis” (86). Binomial GLMs were employed to assess relative module completeness with respect to health condition and the indices CDAI (32) and SCCAI (31).

Data availability.

Raw sequence data and relevant metadata can be accessed online at the European Nucleotide Archive (https://www.ebi.ac.uk/ena/) under the accession number PRJEB44440 for CD and UC cases, as well as under PRJEB47162 for the healthy controls.