Fecal Microbial Composition of Ulcerative Colitis and Crohn’s Disease Patients in Remission and Subsequent Exacerbation

Background Limited studies have examined the intestinal microbiota composition in relation to changes in disease course of IBD over time. We aimed to study prospectively the fecal microbiota in IBD patients developing an exacerbation during follow-up. Design Fecal samples from 10 Crohn’s disease (CD) and 9 ulcerative colitis (UC) patients during remission and subsequent exacerbation were included. Active disease was determined by colonoscopy and/or fecal calprotectine levels. Exclusion criteria were pregnancy, antibiotic use, enema use and/or medication changes between consecutive samples. The microbial composition was assessed by 16S rDNA pyrosequencing. Results After quality control, 6,194–11,030 sequences per sample were available for analysis. Patient-specific shifts in bacterial composition and diversity were observed during exacerbation compared to remission, but overarching shifts within UC or CD were not observed. Changes in the bacterial community composition between remission and exacerbation as assessed by Bray-Curtis dissimilarity, were significantly larger in CD versus UC patients (0.59 vs. 0.42, respectively; p = 0.025). Thiopurine use was found to be a significant cause of clustering as shown by Principal Coordinate Analysis and was associated with decreases in bacterial richness (Choa1 501.2 vs. 847.6 in non-users; p<0.001) and diversity (Shannon index: 5.13 vs. 6.78, respectively; p<0.01). Conclusion Shifts in microbial composition in IBD patients with changing disease activity over time seem to be patient-specific, and are more pronounced in CD than in UC patients. Furthermore, thiopurine use was found to be associated with the microbial composition and diversity, and should be considered when studying the intestinal microbiota in relation to disease course.


Introduction
Crohn's disease (CD) and ulcerative colitis (UC) are chronic inflammatory diseases of the gastrointestinal tract, collectively referred to as inflammatory bowel disease (IBD). IBD is a heterogenous disease ith respect to disease location, disease course, occurrence of extra-intestinal manifestations and therapeutic response. The disease course is characterized by exacerbations and remissions. IBD is generally considered to arise from the interaction between host genetics, environmental factors, dysregulated immune responses and alterations in the intestinal microbiota composition [1]. IBD, especially active disease, is associated with a decreased quality of life and high health care costs, [2,3] especially due to the use of medication [2,4]. Treatment is mainly based on symptom reduction by nonspecific immune-modulating drugs and can be associated with serious side effects [5]. Further insight in causative factors associated with the development of exacerbations (i.e. active mucosal inflammation) may contribute to new specific treatment options for IBD.
The intestinal microbiota is considered to play a central role in the pathogenesis of IBD and numerous studies have corroborated evidence for intestinal dysbiosis in IBD patients compared to healthy controls [6][7][8][9]. The gut microbiota of healthy individuals is dominated by the bacterial phyla Firmicutes and Bacteroidetes, and to a lesser extent by Proteobacteria, Actinobacteria and Verrucomicrobia [10]. In IBD patients, members of the Firmicutes phylum appear to be reduced, [9,11,12] whereas members of Gammaproteobacteria seem to bloom [13][14][15][16]. In CD patients the Clostridia cluster IV group, in particular Faecalibacterium prausnitzii, has shown to be decreased [15,17]. Members of Clostridia group XIVa, belonging to the Roseburia genus, also seem to be decreased in all IBD patients [17][18][19]. Data on Bacteroidetes are more ambiguous; inconsistent findings have been reported for their presence in IBD compared to controls [20][21][22][23][24][25]. In addition to these differences in relative abundances of specific phylotypes, there appears to be a general decrease in biodiversity in IBD patients [25,26].
Altogether, these studies provide compelling evidence for changes in the gut microbial communities in IBD patients as compared to healthy controls. Only a few studies, however, investigated the intestinal microbiota in relation to disease activity. Differences in bacterial species or groups (e.g. F. prausnitzii, Clostridia, E. coli and Fusobacterium varium) were found, when comparing active with inactive UC and CD patients, [7,9,20,[27][28][29] but could not always be confirmed [17]. These studies are prone to confounding effects due to their cross-sectional designs and differences in medication use. Studies monitoring the microbiota composition within patients with changing disease activity over time and controlling for medication use are currently lacking. Therefore, the aim of the present study was to prospectively monitor the microbial composition in UC and CD patients during remission and subsequent exacerbation by means of 16S rDNA pyrosequencing.

Population and Design
The current study was conducted within the context of a prospective follow-up cohort of IBD outpatients, in which standardized demographic and clinical data, feces and blood samples were collected at study entry, at every subsequent outpatient visit and during an exacerbation [30].
The diagnosis of IBD was based on clinical and endoscopic or radiological findings. For case definition of CD, the Lennard-Jones criteria were applied [31]. UC was defined as continuous mucosal inflammation with granulomata, affecting the rectum and/or some or all of the colon in continuity with the rectum. Data of concomitant use of medication, body mass index (BMI), disease duration since diagnosis and disease phenotype (Montreal classification) were obtained using the computer-based medical registration databases. Furthermore at each visit, the occurrence of infections was checked by medical history and culture if indicated, and questions on overall changes in dietary habits were completed.
Fecal samples were collected at home on the evening before or the morning of each visit and stored at 4uC. Upon arrival at the hospital, the faecal samples were split directly. Part was sent to the laboratory of Clinical Chemistry for routine analysis of fecal calprotectine and the remaining part was frozen at 280uC within 24 hours after defecation for analyses of the microbiota.
Within the prospective cohort, patients with faecal samples at remission as well as after subsequent development of an exacerbation during a one-year follow-up period were eligible for the present study.
As clinical indices do not correlate very well with endoscopic scores [32,33] and as we aimed to assess the association of the fecal microbiota with the presence of mucosal inflammation, an exacerbation was defined by colonoscopy, a calprotectin value . 150 mg/g feces or .5-fold increase in calprotectin as compared to remission levels. For CD, only patients with colonic involvement were included as fecal samples are not likely to reflect microbial composition in the small intestine. For UC, eligible patients had pancolitis, left-sided colitis or proctosigmoiditis. Exclusion criteria were pregnancy, use of rectal enemas, use of antibiotics in the 3 months prior to inclusion or during follow-up, and/or a change in immunosuppressive medication between sampling at remission and subsequent development of disease exacerbation.

Ethics Statement
The study was approved by the Medical Ethics Committee of the Maastricht University Medical Center+ (NL31636.068.10), and written informed consent was obtained from all subjects.

16S rRNA V1-V3 Amplicon Library Preparation
Approximately 200 milligram feces was used for DNA isolation by the PSP SPIN Stool DNA plus kit (Stratec Molecular GmbH, Berlin, Germany) according to the manufacturer's instructions and eluted in a final volume of 200 mL.
PCR amplifications (in a volume of 50 mL) were performed using 16FastStart High Fidelity Reaction Buffer, 1.8 mM MgCl 2 , 1 mM dNTP solution, 5 U FastStart High Fidelity Blend Polymerase (from the High Fidelity PCR System (Roche, Indianapolis, USA)), 0.2 mM reverse primer, 0.2 mM of the barcoded forward primer (unique for each sample) and 1 mL of template DNA. PCR was performed using the following cycle conditions: an initial denaturation at 94uC for 3 min, followed by 25 cycles of denaturation at 94uC for 30 s, annealing at 51uC for 45 s and extension at 72uC for 5 min, and a final elongation step at 72uC for 10 min. Amplicons (20 mL) were purified using AMPure XP purification (Agencourt, Massachusetts, USA) according to the manufacturer's instructions and eluted in 25 ml 16 low TE (10 mM Tris-HCl, 0.1 mM EDTA, pH 8.0). Amplicon concentrations were determined by Quant-iT Pico-Green dsDNA reagent kit (Invitrogen, New York, USA) using a Victor3 Multilabel Counter (Perkin Elmer, Waltham, USA). Amplicons were mixed in equimolar concentrations to ensure equal representation of each sample. A one-region 454 sequencing run was performed on a GS FLX Titanium PicoTiterPlate with a GS FLX pyrosequencing system (Roche, Branford, USA).

Data Analysis
The V1-V3 16S rDNA bacterial sequences analyzed in this paper have been deposited in the MG-RAST database (project ID: 4728).
Raw pyrosequencing reads were initially passed through quality filters to reduce the overall error rate using Mothur version 1.23 [35]. Only those sequences with perfect proximal primer fidelity and a threshold quality score of $20, a read length between 200 and 540 nucleotides, a maximum of one ambiguous base call and a maximum homopolymer length of 9 were retained for further analysis.
Barcodes were used to identify sequences belonging to each patient sample. The UCLUST algorithm was used to cluster sequences into operational taxonomic units (OTUs) or phylotypes based on 97% similarity (species level) against the Greengenes reference set [37]. The following nondefault search parameters for UCLUST were applied: maxrejects = 100 and stepwords = 16. Creation of new clusters for sequences that did not cluster to reference sequences within the given similarity threshold was disabled, to further reduce the influence of pyrosequencing errors.
Species richness and diversity within communities (alphadiversity) was measured by means of the observed OTUs (observed richness), Chao1 index (estimated richness), and the Shannon diversity index. Beta-diversity or diversity shared across patient communities was determined by UniFrac distance and Bray-Curtis dissimilarity (BC). UniFrac distances are based on the fraction of branch length shared between two communities within a phylogenetic tree constructed from the 16S rRNA gene sequences from all communities being compared. A relatively small UniFrac distance implies that two communities are compositionally similar, harboring lineages sharing a common evolutionary history [38]. Both weighted and unweighted UniFrac were used, respectively taking into account the abundance of each bacterial species for creating branch length or not.

Statistical Analysis
BC and Unweighted UniFrac distances were calculated for each patient (remission versus active sample) and distances were subsequently compared between CD and UC patients using the Mann-Whitney U test. Distances between subsequent samples from the same patient (within subjects) and distances between remission samples of different patients (between subjects) were also compared using the Mann-Whitney U test. The above mentioned statistical analyses were conducted using SPSS version 20.
All subsequent statistical analyses were conducted within QIIME 1.5.1, unless stated otherwise.
Age was dichotomized in two groups with 52 years (median age) as cut-off. Differences in alpha-diversity with respect to IBD type, disease location, age, gender and medication use (thiopurines, TNF-a inhibitors, aminosalicylates, prednisone or methotrexate) were tested using a nonparametric two-sample t-tests (i.e. using Monte Carlo permutations to calculate the p-value) at a rarefaction intensity of .6,190 sequences per sample. To test for differences in alpha-diversity within patients with respect to disease activity, the Wilcoxon Signed-Rank test was applied (conducted within SPSS v. 20).
OTU coverage was estimated using the conditional uncovered probability method [39]. Analysis of similarities (ANOSIM) was used to test for differences in community structure (UniFrac) among the various groups (IBD type, disease location, age, gender and medication use) in remission state. ANOSIM is a permutation-based test of the null hypothesis that within-group distances are not significantly smaller than between-group distances [40].
To test for associations between OTUs (presence) and IBD type, disease location, age, gender and medication use, the G-test statistic was used. To test for differences in the relative abundance of OTUs, the paired-sample T-test was used with respect to disease activity (within group comparisons) and ANOVA was used to test for differences with respect to IBD type, disease location, age, gender and medication use between groups. The false discovery rate (FDR) method was used to correct for multiple testing at a cutoff of 0.25. The relatively high value was used so as not to miss possible associations.

Study Population
Out of the 323 patients from the prospective IBD cohort, with an average of three subsequent samples per subject, ten CD and nine UC patients that fulfilled our criteria could be selected. Patient characteristics are presented in Table 1. All patients, except one, had calprotectin levels above 150 mg/g feces at time of active disease. The single patient with a calprotectin level below 150 mg/g, showed a 9-fold increase (from 14 toward 126 mg/g faeces) during the flare as compared to remission.

Sequencing
In total, 230,026,772 bases and 551,607 sequences were recovered following pyrosequencing. After trimming, filtering and binning 324,085 sequences (mean length 6SD: 310672.7 bases), ranging from 6,194 to 11,030 sequences per sample, remained for downstream analysis. A total of 2,839 OTUs were found with an OTU coverage of 97.1% 61.1% (mean 6 s.d.).
Patient-specific shifts in bacterial composition were observed in all patients. In a subgroup of patients, major shifts were observed for specific bacteria. Two CD patients (labeled #10 and 11 in Figure 1) showed strong increases in the relative abundances of Bacteroides fragilis during active disease, i.e. from 0.1% to 10.1% and from 0.1% to 81.9% respectively, whereas another CD patient (#17) showed a large increase in relative abundance of Akkermansia muciniphila (1.6% to 62.9%) during activity. Similarly, one UC patient (#8) showed a strong increase in the relative abundance of Bacteroides dorei (2.6% to 18.3%) and another UC patient (#3) showed an increase (1.2% to 13.9%) of an undefined Lactobacillus species. Despite these shifts in the individual microbial profiles, no significant difference in presence or relative abundance of any particular species or group was detected based on disease activity on a group level neither in the CD nor in the UC patients. Furthermore, no sequences from pathogenic species such as Campylobacter spp., Helicobacter spp., Salmonella spp. or Mycobacterium spp. were detected.

Effects on Richness and Diversity
No significant differences in any alpha-diversity metrics were found comparing subgroups based on IBD type, or disease activity (Table 2), nor for disease location, age, gender or use of TNF-a inhibitor, aminosalicylate, prednisone or methotrexate (data not shown). Alpha-diversity seemed to be affected in only one CD patient (#11) during follow up. This patient showed a large relative increase of fecal B. fragilis abundance during relapse compared to remission as reported above.
Thiopurine was used by 7 patients both at remission and during exacerbation (median treatment duration 49.5 (2-182) months), and was found to be associated with a lower number of observed OTUs, Chao1 and Shannon index ( Table 2). The median daily dose expressed as 6-mercaptopurine (6-MP) equivalences [41] was 57.9 mg per day.

Effects on Community Membership and Composition
Unweighted UniFrac-based principal coordinate analysis (PCoA) did not show clustering according to disease activity, IBD type, TNF-a inhibitor, aminosalicylate, prednisone or methotrexate use (Figures S2, S3, S4). Age (R 2 = 0.16; p = 0.001), thiopurine use (R 2 = 0.41; p = 0.001), gender (R 2 = 0.35; p = 0.003), disease location of UC (R 2 = 0.34; p = 0.008) were related to clustering (Figure 3 and Figure S1), as tested by ANOSIM. For thiopurine use, Fusobacteria and Verrucomicrobia phyla appeared to explain clustering for users, whereas especially Lentisphaerae appeared to determine clustering for non-users ( Figure S5).   Unweighted UniFrac distances of patients were significantly closer (p,0.001) between their own remission and exacerbation samples than to samples of other patients (Figure 4). Only the remission sample of patient #11 showed more similarity to the other patients' remission samples than to his own exacerbation sample.
Interestingly, the mean distance in community membership between paired Crohn's disease samples (0.51; s.d. 0.12) was significantly larger as compared to UC samples (0.42; s.d. 0.04; p = 0.027). Using Bray-Curtis dissimilarity distances, the same trend was observed with CD at 0.59 (s.d. 0.18) and UC at 0.42 (s.d. 0.09) (figure 5). This difference between the means was also statistically significant as determined by Mann-Whitney U test with p = 0.027. Several associations were found between specific OTUs and patient characteristics as determined by G-test and ANOVA, but none of these withstood correction for multiple testing by means of FDR (Tables S2, S3, S4, S5).

Discussion
To our knowledge, this is the first prospective study using nextgeneration sequencing to examine the fecal microbiota composition in UC and CD during a quiescent disease phase and a subsequent exacerbation.
Similar to previous studies on the human fecal microbiota, the divisions of Firmicutes and Bacteroidetes predominated [42][43][44]. Although no overarching shifts in microbial communities could be  detected in the CD and UC populations, shifts in composition and decreased diversity were observed in individual patients. Moreover, it was demonstrated that the microbiota within CD patients was less stable compared to the microbiota from UC patients. Finally, thiopurines appear to affect the gut microbial composition and diversity.
In the present study, we were not able to identify overall patterns in microbial changes related to the presence of an exacerbation. An exacerbation was defined by endoscopic scores and/or fecal calprotectin levels as indicators for mucosal inflammation, as clinical scores are reported not to correlate very well with endoscopy findings [32,33]. Mucosal inflammation is increasingly recognized as outcome parameter for active disease and is important for optimization of therapeutic strategies and to prevent complications [45].
Previous studies did show microbial differences comparing active with inactive IBD patients. For example, a combination of leucocytes in the fecal-mucosal transition zone and the concentration of F. prausnitzii were found to be strong predictors of active CD and UC [7]. Diminished F. prausnitzii, as well as a lower Firmicutes/Bacteroidetes ratio was also associated with active disease in another study on 49 IBD patients [9]. Furthermore, a decrease in the Clostridium family was found in active UC and inactive and active CD compared to inactive UC and healthy controls [20]. A decrease of C. coccoides, C. leptum subgroup, Atophium cluster and B. ovatus was reported in feces of active compared to inactive UC, [8] while increased counts of bifidobacteria, E. coli and Clostridia were found in tissue specimens of active compared to inactive UC by others [27].
The above-mentioned studies had a cross-sectional design and are therefore more vulnerable to selection bias and confounding. None of these studies did match patients with quiescent disease to patients with active disease based on characteristics such as disease location, medication, diet and other possible confounders, nor did they adjust for these factors in their statistical analyses. In contrast, two more recent studies that used pyrosequencing did not show clustering based on disease status. A study by Willing et al. used a multivariate statistical approach (partial least-squares discriminant analysis) to distinguish disease phenotypes based on microbial composition in 40 twin pairs [17]. Six patients had active disease based on clinical scores, which could not be distinguished from 34 patients with inactive disease. In the second study, analyzing 231 biopsies and stool samples of IBD patients and healthy subjects, disease activity was not associated with specific shifts in the microbiome composition even after adjusting for other factors such as sample type (i.e. stool or biopsy), IBD type, age, smoking and medication [29].
A potential explanation for the lack of overarching shifts in microbial composition at time of exacerbation in our study may be due to the relatively small study population. Although we included only CD patients with colonic disease and excluded IBD patients with changes in medication over time, some heterogeneity in the study population cannot be excluded. IBD can be considered a set of diseases with overlapping phenotypes rather than a single disease [46,47]. This may contribute to the patient-specific microbial composition shifts, and lack thereof at group level. The fact that from a large cohort of 323 subjects, only 19 fulfilled our selection criteria, illustrates that it is very difficult to include a rather large homogenous subpopulation.
The main strength of the present study is the prospective design comparing the microbiota in association to disease activity within the same patients followed over time, particularly considering the interindividual differences in microbiota composition. Clear subject-specific changes have been observed. As modifications of medication use over time were not allowed, these changes were most likely associated with the changing disease activity. Major increases (up to 81.9%) in specific bacterial subspecies were observed in five patients, being Lactobacillus spp. (in one UC patient), B. dorei (in one UC patient), A. muciniphila (in one CD patient), and B. fragilis (in two CD patients (#10 and #11)). These two latter patients were both non-smoking females (age 35 and 50 years respectively) using a combination of thiopurines and TNFinhibitors continuously during the period of interest. Apart from developing an exacerbation, no major changes in clinical or environmental factors did occur during the follow-up period in any of these patients that could have contributed to the observed findings. Although larger patient numbers need to be studied before final conclusions can be drawn on the role of the individual commensal in the pathogenesis of IBD, the potential importance of B. fragilis is supported by findings from others. B. fragilis was found to be significantly more abundant in biopsies from CD patients compared to UC patients and healthy subjects [23]. In another study, the mucosal biofilm mass in IBD patients appeared to be dominated by B. fragilis [24]. Some strains of B. fragilis are enterotoxigenic by producing a zinc-dependent metalloprotease toxin, and appear to be more prevalent in active IBD patients as compared to inactive patients and control subjects [48]. Animal models have confirmed that enterotoxigenic B. fragilis alone is sufficient to induce colitis [49]. Furthermore, a potential role for the eukaryotic-like ubiquitin gene in B. fragilis in the pathophysiology of CD has recently been postulated [50].  In the present study, we found a different level of microbial variability contributing to exacerbations in either CD or UC. This may indicate that changes in gut bacterial composition play a more important role in CD than in UC exacerbation development. These findings are in line with findings by others. Willing et al. (2010) found different microbial profiles between CD and controls but not between UC and controls based on 454 sequencing [17]. In addition, Andoh et al. found inactive UC patients to cluster with healthy controls, while inactive and active CD clustered separately together with active UC [20].
Although the numbers were small, we found that thiopurine treatment was associated with the microbial composition and the diversity. Based upon the Principal Coordinate analysis, thiopurine use was found to be the most important factor responsible for clustering and was furthermore associated with a significant decrease in alpha-diversity. It has previously been found that the thiopurines 6-mercaptopurine (6-MP) and azathioprine inhibit growth of M. avium subspecies paratuberculosis in vitro [51,52]. We could not demonstrate an antimicrobial effect of 6-MP in vitro as tested on Bacillus subtilis by disc diffusion test (test range 1-64 mg, data not shown). This does however not exclude the possibility of an (indirect) effect of 6-MP or one of its metabolites on the intestinal microbiota in vivo.
As such, the present findings indicate that thiopurine use, especially changes in administration or dosage over time, should be considered as a potential confounding factor when studying the microbiota composition in IBD. Replication of our findings on the association between thiopurine use and shifts in microbiota composition as well as more mechanistic insight on the potential underlying mechanisms is warranted to prove a causal relationship.

Conclusion
Although we did observe patient-specific shifts in microbial composition, we could not demonstrate general changes in microbial composition or diversity in IBD patients at time of exacerbation as compared to a quiescent disease state. Patientspecific shifts in fecal microbial composition seemed to be associated with larger Bray-Curtis dissimilarity between remission and active disease for CD as compared to UC. Furthermore, thiopurine use was found to have significant impact on the microbial composition and diversity, and should be considered when studying the intestinal microbiota in relation to disease course.
Larger prospective studies that enable controlling for potential confounders such as medication, disease subtype and location, are needed to unravel potential general shifts in the microbiota related to exacerbations in IBD. Figure S1 Communities clustered using Principal Coordinates Analysis (PCoA) of the unweighted UniFrac distance matrix. (A) PC1 and PC2 are plotted on xand y-axes. Each point corresponds to a community colored according to age (red, ,52 years; blue $52 years). All samples are shown. The percentage of variation explained by the plotted principal coordinates is indicated on the axes. (B) PC1 and PC2 are plotted on x-and y-axis. Each point corresponds to a community colored according to gender (blue, male; red, female). (TIFF) Figure S2 Communities clustered using Principal Coordinates Analysis (PCoA) of the unweighted UniFrac distance matrix. (A) PC1 and PC2 are plotted on xand y-axes. Each point corresponds to a community colored according to aminosalicylate use (red, no aminosalicylate use, blue, aminosalicylate use). All samples are shown. The percentage of variation explained by the plotted principal coordinates is indicated on the axes. (B) PC1 and PC2 are plotted on x-and y-axis. Each point corresponds to a community colored according to IBD type (red, CD; blue UC). (TIFF) Figure S3 Communities clustered using Principal Coordinates Analysis (PCoA) of the unweighted UniFrac distance matrix. (A) PC1 and PC2 are plotted on xand y-axes. Each point corresponds to a community colored according to TNF-a inhibitor use (red, TNF-a inhibitor use; blue, no TNF-a inhibitor use). All samples are shown. The percentage of variation explained by the plotted principal coordinates is indicated on the axes. (B) PC1 and PC2 are plotted on x-and y-axis. Each point corresponds to a community colored according to disease activity (red, exacerbation; blue, remission). (TIFF) Figure S4 Communities clustered using Principal Coordinates Analysis (PCoA) of the unweighted UniFrac distance matrix. (A) PC1 and PC2 are plotted on xand y-axes. Each point corresponds to a community colored according to methotrexate use (red, methotrexate use; blue, no methotrexate use). All samples are shown. The percentage of variation explained by the plotted principal coordinates is indicated on the axes. (B) PC1 and PC2 are plotted on x-and y-axis. Each point corresponds to a community colored according to prednisone use (red, no prednisone use; blue, prednisone use). (TIFF) Figure S5 Communities plotted together with phyla using Principal Coordinates Analysis (PCoA) of the unweighted UniFrac distance matrix. PC1, PC2 and PC3 are plotted on x-, yand z-axes. Each point corresponds to a community colored according to thiopurine use (red, thiopurine use; blue, no thiopurine use). Note: Axes were deliberately skewed to most clearly depict relation of clustering to bacterial groups. (TIFF)