Comparative analysis of the subjects' intestinal microbiota over time

To obtain an overview of the similarity of the total GI microbiota between all the study subjects, we unsupervisedly grouped the preprocessed HITChip profiles (see Methods and Figure S1) using hierarchical clustering with correlation distance method (Figure 1). The largest differences in the microbiota profiles were the inter-individual differences (mean inter-individual Pearson correlation; r=0.78, standard deviation; ±0.04). All samples clustered in a subject-wise manner (Figure 1), and a high overall stability of the individuals' microbiota was observed (mean intra-individual r=0.96, ±0.02 over all time points).

Open in a separate window

Figure 1

Hierarchical clustering of the HITChip profiles of 15 subjects and their six timepoints.

The subject-wise clustering is highlighted with boxes and the temporal variation is displayed with the variable length of the branches. Red vertical line is drawn at Pearson correlation of 0.925, below which the intra-individual sample similarity was detected only in the five unstable subjects (5,6,8,9 and 11) discussed in the text.

Despite the overall stability we detected substantial variation in the temporal dynamics of the subjects' microbiota as visualized by the variable lengths of the branches in the clustering tree. Subject 9 was on an antibiotic course during the collection of the fourth sample, which separated clearly from the other samples (Figure 1). Also the preceding sample taken two weeks earlier branched deeply in the dendrogram implying a distinct microbiota profile. At that time the subject was not medicated but was ill and few days later diagnosed with streptococcal pharyngitis. Apart from the antibiotic treatment, we observed no correlations between the microbiota composition or stability and the recorded medication (subjects 1,3,7 and 11 received regularly thyroid or estrogen hormones, or hypertension treatment). A decreased temporal stability was also observed in subjects 5, 6, 8 and 11 (Figure 1). For all these subjects, with the exception of subject 5, the decreased stability could be explained in the light of the background information. Subjects 6, 8 and 11 had records of travelling during the first weeks of the trial, passing different time zones in a range of 1–7 hours. The similarity between the samples taken before and after the travelling was statistically significantly lower than between the samples taken two to four weeks after travelling (r=0.90, ±0.02 and 0.96, ±0.02, respectively; t-test p-value <0.05) (Figure S2). Accordingly, the travelers showed significantly lower intra-individual stability of the GI microbiota than the nine healthy subjects (r=0.93±0.01 and 0.96, ±0.02, respectively; t-test p-value <0.05), suggesting that the GI microbiota was affected by the travel. Therefore, although the travelers were healthy, we classified them in this context as compromised together with the three subjects affected by an antibiotic treatment or persistent intestinal symptoms (Table 1). The rationale for such classification was to focus the compositional microbiota analysis on those healthy subjects for which we did not observe any traceable destabilizers.

Compositional microbiota analysis

The GI microbiota in the nine healthy uncompromised subjects was analyzed in detail using HITChip signal intensities of species-like taxa that we hereafter term phylotypes (level 3 phylogeny in Methods and [13]). We detected on average 470 phylotypes per subject, representing 68.4% of the total number of the phylotypes that passed the signal treshold (Figure S1). The most abundant phylum in the study subjects consisted of Firmicutes with over 80% contribution to the total microbiota, followed by an average 10% of Bacteroidetes and 1.5% of Actinobacteria (Figure 2). Verrucomicrobia and Proteobacteria were detected with proportional share below 1% (data not shown). As expected, the most abundant Firmicutes were the members of Clostridium clusters XIVa and Clostridium cluster IV, comprising of an average proportion of 40% and 35% of the total microbiota, respectively.

Open in a separate window

Figure 2

Microbiota composition of the study subjects.

Relative abundance of the phylotypes that contributed over 0.5% to the total HITChip signal was summed up to phylum level, except Firmicutes, which were summed up to class or Clostridium cluster level (left y-axis). The qPCR-based quantification of the methanogens is shown on the right y-axis.

Based on Methanobrevibacter–specific quantitative PCR (qPCR), intestinal methanogens were present in 80% of the individuals (12/15). In one subject (subject 2; Figure 2) the methanogens were detected only in a single timepoint, while other 11 carriers had them in all the timepoints. The abundance of the methanogens varied over two orders of magnitude between the carriers, ranging from log10 counts of 7.1 to 10.2 g-1 feces. The relative proportion of the methanogenic archaea in the total GI microbiota was determined by the ratio of the 16S rRNA gene copy numbers of the methanogens and the total bacteria. The proportion of the methanogens was below 1% in most of the samples while in subject 14 the methanogens represented 1–5% of the total community across the timepoints.

Temporal variation of the healthy microbiota and the individual core

We analyzed the effect of time on the microbiota variation in the nine healthy subjects by determining the similarity of the HITChip profiles between the intra-individual samples. The similarities were determined by a correlation coefficient computed for all possible time intervals from one to seven weeks. The mean intra-individual Pearson correlations ranged from 0.96 (±0.02; one week interval) to 0.95, (±0.02; seven week interval), showing only a slight, non-significant, decline in the microbiota similarity over time. In order to address the individual core consisting of microbes that retain over time, we determined the taxon-specific temporal variability by calculating the Coefficient of Variation (CoV) within a subject for each phylotype. The CoVs ranged from 0.83% (uncultured Ruminococcus callidus-like bacterium) to 46.1% (phylotype within uncultured Clostridiales) with an average of 6.3% (±1.1%), reinforcing the high overall stability. Among the 50 most stable genus-level taxa per subject, the most common statistically enriched genus-level groups were Ruminococcus obeum and relatives (et rel.), affiliated to the recently renamed genus Blautia [33] with an average CoV of 3.7%, and Clostridium symbiosum et rel. with an average CoV of 4.4% (Table S1). No other stable genus-level taxon was shared by >50% of the healthy subjects, indicating subject-specificity of the stable species, i.e. of the individual core microbiota. Similarly, we determined the 50 statistically significant, temporally most unstable phylotypes per each healthy subject (Table S2). The most common unstable genus-like taxa were uncultured Clostridiales and Bacteroides vulgatus et rel. (CoV 9.9%), but even they were common to only few individuals. Not only the commonness but also the actual CoVs of the unstable groups had larger inter-individual variation than the stable ones, indicating that no genus-level group was generically unstable. Few genus-like groups, including F. prausnitzii et rel., were in the top list of both stable and unstable taxa (Tables S1 and S2), highlighting the large individual variation in the temporal behavior of these taxa. In contrast, the methanogenic archaea appeared temporally very stabile among the carriers, as the intra-individual CoVs calculated from the qPCR data were limited to a range of 1% to 6%.

Correlation between intestinal symptoms and the GI microbiota

Despite the good general health of the participants, only one out of the 15 subjects did not experience any intestinal complaints during the study period (subject 1; Table 1). However, occasional and mild intestinal problems dominated the answers as only two subjects had intestinal health score below 90 out of 100 (Table 1). To detect the bacteria associated to specific intestinal symptoms, we correlated the HRQoL–questionnaire scores and the microbiota using the HITChip data from all the 15 subjects. We could identify 170 phylotypes, which correlated statistically significantly (q<0.005) with 14 HRQoL-questions related to intestinal symptoms (Figure 3). The significant correlations were most frequently found for general intestinal discomfort, urgent need for defecation, postprandial fullness, and abdominal pain. About half (53%) of the 170 phylotypes correlated negatively with the symptoms, showing higher abundance in the absence of the symptoms. The negatively correlating genus-like taxa included Bifidobacteria, R. obeum et rel., Streptococcus bovis et rel., Clostridium difficile et rel. (mostly due to the nonpathogenic Clostridium irregularis) and Ruminococci. The other half showed a positive correlation with the observed symptoms so that the probe intensities were higher when the symptoms occurred. The positively correlating phylotypes represented mainly members of the Clostridium clusters IV and XIVa, among which uncultured members of the F. prausnitzii et rel. and C. symbiosum et rel. were the most predominant groups.

Open in a separate window

Figure 3

Correlation heatmap of the microbiota profiles and intestinal symptoms.

The heatmap visualizes all significant (q-value <0.005) Spearman correlations between the intestinal symptoms and the HITChip-measured abundance of bacterial phylotypes. Each row represents a single phylotype, whose genus-level assignments are listed below the rows, or in the middle when several rows refer to the same genus-level group. The negative correlations between symptoms are indicated in red and the positive correlations with yellow; non-significant correlations are shown in gray.

We next aimed to identify the bacteria related to the symptoms that were most frequently experienced in our survey, namely bloating (31 complaints from 8 subjects) and abdominal pain (28 complaints from 8 subjects). For bloating, the severity of the recorded complaints varied from strong discomfort (2 records) to mild symptoms (23 records). Several uncultured phylotypes from Clostridium clusters IV and XIVa had statistically significant positive correlation with bloating (Table S3). The abundance of phylotypes within Anaerotruncus colihominis et rel., Ruminococcus callidus and Lachnospira pectinoschiza et rel. were over ten-fold higher when bloating was recorded.

The subjects experienced abdominal pain in varying degree so that the severity of the recorded complaints ranged from very strong disturbance (1 record) to mild symptoms (20 records). The most striking finding related to the abdominal pain was its negative correlation with the abundance of Bifidobacteria (r=−0.45, ±0.03; Table S4 and Figure 3). The majority of the Bifidobacterium phylotypes that were above the detection limit (14/21) had statistically significant (q≤0.02) negative correlation with abdominal pain. The subjects who experienced pain had over five-fold (on average 5.4, ±1.3) less Bifidobacteria compared to the subjects without recorded pain (Figure 4A). Moreover, a genus-specific qPCR, covering most cultured members of this group, revealed a statistically significant difference and a fold change of two (Wilcoxon test p<0.05) in the Bifidobacterial counts according to the symptom status (Figure 4B).

Open in a separate window

Figure 4

Quantification of Bidifobacteria in relation to abdominal pain.

The samples were divided into two groups: with and without concurrent abdominal pain, and the abundance of Bifidobacteria in each group was determined using A. HITChip and B. qPCR. The difference between the groups was statistically significant (q≤0.02 in A; p<0.05 in B). The box extends from 25th percentile to 75th percentile, with a line at the median; the whiskers extent to the highest and lowest values.

In addition to Bifidobacteria, one phylotype within R. lactaris et rel., Clostridium cluster IV, was significantly decreased in pain-associated samples (Table S4). Among the positive correlations, uncultured, potentially pathogenic phylotypes within Uncultured Clostridiales II, Anaerotruncus colihominis et rel. and R. callidus, were increased over ten-fold when pain was recorded (Table S4). C. symbiosum and uncultured phylotypes affiliated to F. prausnitzii et rel. were more abundant in subjects experiencing postprandial fullness (Figure 3). The positive correlation between postprandial fullness and the abundance of F. prausnitzii was confirmed with qPCR (data not shown). The symptom scores of bloating and abdominal pain correlated strongly (ρ=0.72). Four genus-level Clostridial groups were significantly correlated with the both symptoms, albeit in a non-uniform direction (Tables S3 and S4).

The common core microbiota

In-depth phylogenetic microarray analysis provided the opportunity to explore the common core microbiota that is potentially shared between the healthy individuals. Most of the present efforts in determining and defining the microbial core have focused on the presence and absence of phylotypes in a defined proportion of the study subjects. Instead of using arbitrary criteria for the abundance and prevalence of the phylotypes, we addressed the common core microbiota using a novel approach that exploits adjustable values of these two parameters. Considering the relative abundance, we utilized the full dynamic range of the HITChip microarray and incorporated the phylotypes that at lowest comprised as little as 0.02% of the estimated total signal while the most dominant phylotypes had over 5% relative abundance. With regard to the prevalence, we accommodated all possibilities between 0–100% (n=0–9). The averaged HITChip data from the six timepoints per person were used. As a result, the size of the common core is not a single number but instead a continuum from zero to hundreds, depending on the selected abundance and prevalence values as visualized with a perspective plot (Figure 5A). At the prevalence of 50%, only 10 phylotypes remained in the common core when the minimum abundance was set as high as 0.5% of the estimated total signal. On the other hand, when all the phylotypes above the detection limit were accepted, changing the prevalence from 50% to 100% had much less effect in the number of the core phylotypes. In other words, the number of the shared phylotypes decreased dramatically when a phylotype was required to be in high abundance in all samples, emphasizing the contribution of the low-abundance bacteria in the core.

Open in a separate window

Figure 5

Definition of the common core microbiota.

A. Perspective plot is used to visualize how the number of phylotypes in the common core is a function of the selected abundance and prevalence. Light gray indicates the area where no phylotypes passed the given criteria. B. The common core microbiota in the healthy subjects. The y-axis represents the range of the phylotype count shared, thus forming the common core, in the nine healthy subjects. The line visualizes how the number of the phylotypes depends on the selected percentile abundance, and how it decreases drastically when the phylotypes with lowest abundance are excluded. List shows phylum-level summary of the phylotypes that contributed with over 0.03% relative abundance C. The common core phylotypes summarized to genus-level taxa. The boxplots visualize the intra-individual variation in the abundances of the core taxa. Bacteroides fragilis et rel. and Bifidobacteria (two leftmost boxes) are visualized for reference purposes. These taxa were not part of the common core due to their large intra- and inter-individual variation. For boxplot details, please see the legend for Figure 4.

By including all the phylotypes above the signal treshold in all the nine healthy subjects, we detected 288 phylotypes in the common core (Figure 5B and Table S5). The 288 phylotypes represented 41.9% of all the detected phylotypes, indicating a high degree of phylogenetic overlap between the individuals. The exact taxon numbers should be regarded with caution as they are subject to noise arising e.g. from unintended cross-hybridization with non-target sequences and differences in the probe binding affinities. To estimate such noise, and to exclude the possibility that related probe-level effects would have remarkably affected the common core composition and size, we used the Robust Probabilistic Averaging algorithm designed for probe performance analysis [34]. For the vast majority of the common core phylotypes (234; 81%) the correlation between the probabilistic and the standard method was >0.9, indicating a low proportion of potentially false positives among the detected core phylotypes. The preprocessing did not affect the phylogenetic composition of the core as the 22 phylotypes detected with most unreliable probes (r<0.8) were all assigned to the diverse and polyphylogenic Clostridium clusters XIVa and IV that vastly dominated the core with 145 and 111 phylotypes, respectively (Table S5). The majority of the core organisms were uncultured phylotypes, belonging to 43 different genus-like taxa from five phyla (Figures 5B and 5C). The most represented groups were R. obeum et rel. (37 phylotypes), F. prausntizii et rel. (25 phylotypes) and O. guillermondii et rel. (20 phylotypes). Bacteroidetes, Actinobacteria, Proteobacteria and Verrucomicrobia were also represented in the core (Figure 5 and Table S5), in line with their established inhabitance in the human gut. Most of the core phylotypes had a highly variable abundance among the healthy individuals (Figure 5C).

Fecal samples and DNA extraction

The fecal samples we collected at home, immediately frozen at −20°C, and transported within 12 hours to the study center where they were stored at −70°C until analyses. Community DNA was extracted from the fecal samples with a repeated bead beating method that has been recently validated [29]. The extracted DNA was quantified with a Nanodrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and analyzed by 0.8% agarose gel electrophoresis for integrity.

Microarray analysis

The HITChip microarray analyses were performed as previously described [13]. In brief, the HITChip consists of over 4 800 oligonucleotide probes targeting the V1 and V6 hypervariable regions of the 16S rRNA gene, spotted in duplo on the custom Agilent arrays (Agilent Technologies, Palo Alto, CA, USA). To prepare the fecal DNAs for the microarray analysis, full-length 16S rRNA gene was amplified using primers T7prom-Bact-27-for and Uni-1492-rev [13]. The PCR products were transcribed into RNA and labeled with Cy3 and Cy5 prior to fragmentation and hybridization on the array. The arrays were scanned with Agilent DNA Microarray Scanner and the data were extracted from the generated images using the Agilent Extraction Software 9.1.

Computational pre-processing of the microarrays

The array normalization was conducted as described earlier [13], [29] including spatial normalization with polynomial regression of the each scanner channel followed by outlier detection and quantile normalization of both of the Cy3 and Cy5 dyes for each sample. The duplicates that showed over 0.98 Pearson correlations between the dyes were considered for further analysis and the dyes were averaged, otherwise the hybridization was reproduced. The between-array normalization of the samples was conducted by using the assumption of normal-exponential (background-signal) distribution and quantile normalization [57]. The array probes were organized based on their 16S rRNA sequences on three levels of phylogeny as described before [13]. In brief, the HITChip signal intensity was analyzed using the following phylogenetic assignment levels: 1) the phylum-level, with the specification of Firmicutes down to Clostridium clusters, creating altogether 23 groups; 2) the genus-like level, including 131 groups of sequences with ≥90% sequence identity, and 3) the phylotype (species-like) level with 1033 distinct phylotypes with ≥98% sequence similarity to cultured species or clones corresponding to uncultured micro-organisms. “Genus-level taxa with ≥90% sequence identity distributed over multiple genera are termed “et rel. ”. To reduce the experimental and potential cross-hybridization noise, and to reliably assess whether certain phylotypes were present in a sample based on the hybridization intensity distributions of all the subjects, we estimated the signal intensity threshold from the data using the approach originally developed for gene expression microarrays [58]. The thus obtained signal intensity threshold was log10 intensity of 1.8 (Figure S1), corresponding to estimated 0.02% of the total signal as a detection limit for each phylotype. To increase the analysis reliability, we additionally determined the probe specificities for each phylotype by calculating the number of the targeted phylotypes per probe in order to specify taxa that are potentially subject to cross-hybridization (i.e. are detected with probes that have more than one target phylotype based on the hypervariable V1 and V6 sequences of the 16S rDNA). Multiple targeted phylotypes were essentially detected only within the same genus-level taxa. Additionally, the Robust Probabilistic Averaging algorithm was used to estimate the reliability of the individual HITChip probes in regard to noise attributable e.g. to unintended cross-hybridization with non-target sequences [34], potentially resulting in false positives in the common core analysis. Previously, a good correlation between the HITChip- and pyrosequencing-derived microbial profiles has been recorded [30], [59].

Quantitative PCR analysis

Triplicate quantitative PCR (qPCR) amplifications of the 16S rRNA gene were performed for all available samples using the MX3005P Real-Time PCR System (Stratagene, La Jolla, CA, USA) and CYBR green chemistry. Total bacteria were quantified by amplifying 0.5 ng of fecal DNA with universal primers [60]. The archaea were quantified using primers specific for the genus Methanobrevibacter [61] and 25 ng of the template DNA. The assays for Faecalibacterium prausnitzii and Bifidobacterium spp. were conducted as described previously [62]. The reaction conditions and construction of standard curves have been previously described [29], [62]. The data are expressed as 16S rRNA gene copies µg-1 template DNA.

We thank all the volunteers who participated in this study and Valio Ltd. for the collection of fecal samples and health questionnaires. We gratefully acknowledge the assistance of Dr. Mirjana Rajilić-Stojanović, Hans Heilig and the HITChip team (Wageningen University, NL) in the HITChip analysis, and MD Peter Sweins in the interpretation of the HRQoL data. The work was performed in the Centre of Excellence on Microbial Food Safety Research, Academy of Finland.