Associations between the human intestinal microbiota, Lactobacillus rhamnosus GG and serum lipids indicated by integrated analysis of high-throughput profiling data

Ethics statement

The trial and its protocol have been approved by the Ethics Committee of the Hospital District of Helsinki and Uusimaa (Ethical protocol no HUS 357/E0/05). The subjects provided written informed consent. The details of this trial have been published previously (Kekkonen et al., 2008a; Kekkonen et al., 2008b).

Subjects

The subjects were healthy Finnish adults (n = 25, 18 females, 7 males) with a mean age of 42 years (23–55) and a mean BMI of 24 kg/m² (18–30) from Helsinki region, representing a subset of a larger cohort (Kekkonen et al., 2008a; Kekkonen et al., 2008b). The CONSORT flowchart (Fig. 1) summarizes the subject enrollment, intervention allocation, follow-up, and data analysis during the study. The subject characteristics are provided in Supplemental Table S1.

Dietary intervention

The subjects constituted two treatment groups in a randomized, double blind intervention study to receive either L. rhamnosus GG (probiotic; n = 11) or placebo (n = 14). During the intervention the subjects consumed daily a 250 mL milk-based fruit drink containing either L. rhamnosus GG (ATCC 53103, 6.2 × 10⁷ cfu/mL) or a similar placebo drink without probiotic bacteria for three weeks. No other probiotic-containing products were allowed three weeks prior or during the intervention; a list of fermented foods and commercial probiotic-containing products was given to the subjects. The subjects were not separately questioned about their abstinence from probiotic products, but they filled a study diary recording the daily intake of the study product.

Fecal samples and microbiota analysis with the HITChip and quantitative PCR

Three fecal samples per individual were collected three weeks before, during and three weeks after the intervention, as previously described (Kekkonen et al., 2008b). The fecal DNA extraction with modified Promega method, and quantification of the probiotic counts with a strain-specific quantitative PCR (qPCR) assay have been described previously (Ahlroos & Tynkkynen, 2009). The phylogenetic analysis of the intestinal microbiota composition with the HITChip microarray was performed as previously described (Rajilić-Stojanović et al., 2009; Salonen et al., 2010). The phylogenetic HITChip microarray targets the V1 and V6 hypervariable regions of the 16S rRNA gene of over 1000 bacterial phylotypes that present the majority of the so far detected phylotypes of the human intestine. Phylogenetic organization of the microarray probes and data preprocessing have been explained in detail elsewhere (Rajilić-Stojanović et al., 2009; Salonen et al., 2010; Jalanka-Tuovinen et al., 2011). Quantification of the relative differences in the taxon abundance between samples was obtained by summarizing the probe signals to phylotype (species-like), genus and phylum levels. Genus-level taxa with ≥ 90% sequence similarity in the 16S rRNA gene are referred to as Species and relatives, the latter being shortened in the text as “et rel.” (Rajilić-Stojanović et al., 2009).

Previously, HITChip-derived microbial profiles have been shown to correlate well with those obtained with fluorescence in situ hybridization (FISH) (Rajilić-Stojanović et al., 2009) and pyrosequencing (Claesson et al., 2009). Samples were also used for absolute quantification of total bacteria, methanogenic Archaea, Lactobacillus group and Bifidobacterium spp. with previously described qPCR primers and reaction conditions (Salonen et al., 2010).

Blood samples and their biochemical analyses

Of the 25 study subjects used for the microbiota analysis, 22 (8 from probiotic group and 14 from placebo group) were available for the parallel analysis of serum lipid profiles before and after the intervention (44 samples in total). Venous blood samples from the antecubital vein were taken at baseline and after the three-week intervention. The blood samples were stored at −20 °C for further analyses. Total cholesterol, high-density lipoprotein (HDL) and low-density lipoprotein (LDL) cholesterol as well as triglyceride levels were enzymatically determined from the serum as previously described (Kekkonen et al., 2008a).

Global lipid profiling: sample preparation and analysis by UPLC-MS

Extraction of lipids from the serum samples, analysis of the lipid extracts with Waters Q-Tof Premier mass spectrometer combined with an Acquity Ultra Performance LC^TM (UPLC), data preprocessing and identification of the lipid molecular species is described in detail elsewhere (Kekkonen et al., 2008a). Lipids have been named according to Lipid Maps (http://www.lipidmaps.org) with the following abbreviations: Cer: ceramide; ChoE: cholesteryl ester; lysoPC: lysophosphatidylcholine; PA: phosphatidic acid; PG: phosphatidylglycerol; PC: phosphatidylcholine; PS: phosphatidylserine; SM: sphingomyelin; TG: triglyceride. Where the fatty acid composition could not be determined, the total number of carbons and double bonds is indicated. The first number indicates the amount of carbon atoms in the fatty acid molecule, followed by the number of double bonds.

Statistical analyses

The effects of probiotic intervention on the intestinal microbiota were investigated with a combination of linear models, sparse principal component analysis (PCA; Shen & Huang, 2008; Lê Cao, González & Dèjean, 2009), unsupervised hierarchical clustering, and significance testing. The log-transformed HITChip and lipid profiling values were approximately Gaussian distributed and fulfilled the general statistical assumptions underlying the selected computational approaches. Background variables, including age, body-mass index and gender were compared in the baseline samples to exclude potentially confounding effects associated with these variables. The Wilcoxon test was used with continuous variables and the Fisher exact test for categorical variables, followed by Benjamini–Hochberg p-value correction for multiple testing. No significant differences between the background variables were observed between the treatment groups (p > 0.05 in all comparisons). The two-group comparisons between the time points and between the treatment groups were quantified based on a linear model with group-wise fixed effects and sample-specific random effects, as implemented in the limma R package (Smyth, 2004), to identify bacterial taxa with significant changes induced by the probiotic intervention and to assess the magnitude and significance of the effects. The function lmFit was used to fit the linear model, followed by significance estimation by empirical Bayes as described in Smyth (2004). In addition, power calculation was carried out to assess statistical power of the current study. The original data was randomly permuted, one significant 2-fold alteration was inserted, and Gaussian noise was added using the same average standard deviation as in the original data. Empirical power calculation with 1000 random permutations showed that 2-fold and higher alterations that follow the noise levels in the original data were detected in > 99.8% of the cases with sample size of 8 or more, based on the same detection criteria than in the current study, confirming the high statistical power of the analysis at the present sample size. PCA is a linear dimension reduction method used to compress information in the high-dimensional phylogenetic (genus-level) and lipid profiles into few informative features that capture the main variation in the data and allow two-dimensional visualization of sample similarities. The temporal and inter-individual similarity of the microbiota and lipid profiles was assessed by average intra- and inter-individual Pearson correlations (r) between and within the time points, respectively. The differences in profile similarity between the probiotic and placebo groups were estimated with Wilcoxon test. The biweight midcorrelation, which is more robust to outliers than the standard Pearson correlation, was used to quantify correlation between the microbiota and lipid profiles across the individuals. High-throughput screening studies involve considerable multiple testing and the traditional multiple testing correction approaches are prohibitively conservative in this context due to their emphasis on estimating the probability of a single false positive finding. Hence, we have used q-values for multiple testing correction in the high-throughput screening tests that include parallel comparisons of large numbers of lipids and bacterial phylotypes.

Correlation analysis of the intestinal microbiota and serum lipids

The biweight midcorrelation measure was used to quantify associations between the microbiota and lipid profiling data sets, and between the microbiota and biochemically determined lipids. The correlations between genus-level bacterial groups and lipid species were calculated together with significance estimates that were corrected for multiple testing. The significantly correlated (q < 0.05) phylotypes and lipids were selected for further analysis. It is important to note, however, that while the present sample size of 22 samples across 2 time points (before and during the probiotic intervention) is sufficient for highlighting significant correlation patterns, the present experimental design cannot uncover causal relationships and the observed correlations may be partly associated with diet or other confounding variables that simultaneously affect both lipid and microbiota profiles. Two-way average linkage hierarchical clustering of the bacterial taxa and lipids was applied to highlight and visualize groups of phylotypes and lipids sharing similar correlation patterns. Analogous ‘second-order’ correlations have previously been applied for instance in the context of gene expression studies (Parmigiani et al., 2004; Lahti et al., 2013). A constant plaid model biclustering was used for systematic detection of significantly correlated bacterial and lipid groups on the correlation heatmap (Lazzeroni & Owen, 2002). To interpret the detected biclusters, statistical enrichment (over-representation) of the implicated bacteria was quantified. Moreover, enrichment analysis was carried out for lipids containing even or odd number of carbon atoms, as well as the enrichment of lipids with zero, one or more double bonds. The enrichment analyses were carried out with Fisher’s exact test (Rivals et al., 2007). All analyses were performed within the R statistical environment (R Development Core Team, 2010).

Impact of the L. rhamnosus GG intervention in the intestinal microbiota

This study characterized the impact of L. rhamnosus GG intervention on the stability and composition of the intestinal microbiota. The subjects in the probiotic group consumed daily approximately 10¹⁰ (10.2log₁₀) colony forming units (cfu) of L. rhamnosus GG. The compliance was verified with the quantification of L. rhamnosus GG in the feces with strain-specific qPCR (Kekkonen et al., 2008b). The average excretion of L. rhamnosus GG in the probiotic group was more than 1000-fold higher than that in the placebo group (8.52; sd 0.73log₁₀ versus 5.20; sd 1.09log₁₀ genome copies per gram of feces, respectively).

The stability of the microbiota during the trial was quantified by the similarity of the microbiota profiles between the three time points with Pearson correlation (Table 1). The average intra-individual correlations were high; 0.94–0.95 (sd 0.02–0.03). No significant difference in the temporal stability of the microbiota between the probiotic and placebo groups was observed, indicating that the probiotic intervention did not alter the overall microbial stability. Principal Component Analysis (PCA) visualization of the relationships between the intestinal microbiota profiles (Fig. 2A) and hierarchical clustering (data not shown) further supported the conclusion that there were no systematic differences in the microbiota between the intervention groups. The average inter-individual microbiota correlation was 0.76–0.78 with no significant differences between the probiotic and placebo groups (Table 1). The intra-individual microbiota correlations (r = 0.94–0.95) were notably higher than the inter-individual correlations (r = 0.76–0.78), stressing the subject-specificity of the microbiota.

Linear models were used to quantify the effects of the L. rhamnosus GG intervention on individual taxa using genus- and phylotype-level microarray data. A specific and transient increase of bacteria related to L. rhamnosus was detected in the probiotic group immediately after the intervention (Fig. 2). There were no intervention-related effects in the other lactobacilli, Bifidobacteria or any other taxa either in the probiotic or placebo group; however substantial individual variation was evident.

To complement the HITChip microarray analysis, we also determined the absolute counts of total bacteria, methanogenic Archaea, Lactobacillus group and Bifidobacterium spp. using real-time PCR. The ingestion of L. rhamnosus GG was reflected in the total lactobacilli that showed a highly significant increase in the probiotic group after the intervention, returning to baseline levels in the follow up (q < 0.05 in both comparisons). No other significant changes were observed in the amount of targeted microbes neither in the probiotic or placebo group. We detected substantial inter- and intra-individual variation in the methanogenic Archaea but that was independent of the time point or the treatment group (data not shown). These observations support the conclusion that the L. rhamnosus GG intervention did not change the overall microbiota composition.

High-throughput profiling of serum lipids

Global profiling of the serum lipids was performed using the UPLC-MS platform to assess the serum lipid profiles at the molecular species level (Kekkonen et al., 2008a). In total, 407 lipid species from 11 different classes were identified. The lipid profiles were dominated with triglycerides (TG; 37%), phosphatidylcholines (PC; 25%) and phosphatidyletanolamines (PE; 13%) while phosphatidic acids (PA), phosphatidylglycerols (PG), sphingomyelins (SM), cholesteryl esters (ChoE), lysophosphatidylcholines (lysoPC), phosphatidylserines (PS), ceramides (Cer) and lysophosphatidylethanolamines (lysoPE) each contributed with 7% to 0.3% to the total lipid pool.

The stability of the lipid profiles was quantified based on the same correlation analysis approach than in the HITChip microbiota profiling analysis. No significant differences in the lipid profile stability were observed between the treatment groups (average intra-individual r = 0.92 and r = 0.93 for the probiotic and placebo groups, respectively; Table 1). In contrast to microbiota profiles, the serum lipid profiles were remarkably similar not only within the subjects but also between the subjects (average inter-individual r = 0.89–0.91; Table 1). No significant differences were seen between the probiotic and placebo groups, or between the time points. The stability of the lipid and microbiota profiles showed a weak positive correlation (r = 0.27) across the subjects independent of a treatment group. The global lipid profiles did not separate according to the treatment group in PCA visualization (Fig. 2B), and no statistically significant lipid alterations were detected between the time points neither in the probiotic or placebo group (data not shown; Kekkonen et al., 2008b). The fold-changes for each lipid species between the two time points followed normal distribution in both intervention groups (data not shown), which further supports the conclusion that the intervention did not induce systematic alterations on the lipid profiles.

Associations between the intestinal microbiota and serum lipids

In addition to studying the effects of probiotic L. rhamnosus GG intervention in the intestinal microbiota, we investigated the co-variation of the microbiota and serum lipid profiles. As discussed above, our results (Table 1, Fig. 2) conclusively show that the probiotic intervention did not affect the lipid profiles or the microbiota beyond the L. rhamnosus that was ingested. Hence, the subjects from probiotic and placebo groups were pooled to increase the sample size and statistical power in exploring the associations between the intestinal microbiota and serum lipids.

Heatmap visualization of the microbiota–lipid correlations across the 22 subjects from 2 time points revealed intestinal bacteria and serum lipids with significantly correlated abundance patterns (Fig. 3). In total, 86 bacterial group-lipid pairs with notable correlations ( ± 0.5 or higher) were identified at q < 0.05 significance level (Fig. 3). Among the significant correlations, 23 of the 131 genus-level taxa detectable by the HITChip were represented (Supplemental Table S3). Most of these taxa belonged to Proteobacteria (6) and Firmicutes within Clostridium cluster XIVa (5) and Bacilli (3). From the lipid side, the vast majority of the significant correlations were attributable to the most dominant lipids, TGs (62%) and PCs (30%). When analyzed at the genus-level, the significant correlations were strongly dominated by the positive correlations between bacteria related to Ruminococcus gnavus and different TG species (31 of 40 significant positive correlations, average r = 0.61; Fig. 3A). Detailed analysis at the phylotype-level data indicated that within the R. gnavus group, two uncultured phylotypes, uncultured human gut bacterium JW1G3 and JW1H4a, dominated the correlations while the type species R. gnavus or related R. torques did not correlate significantly with TG or any other lipid species. Other representatives of the Clostridium cluster XIVa included bacteria related to Dorea formicigenerans that also correlated positively with four TG species (Fig. 3; average r = 0.58). The remaining TG correlations were negative and involved phylogenetically diverse bacteria (Supplemental Table S3). Also PCs correlated significantly with numerous taxa, most of which correlated negatively with ester- or ether linked PCs. Interestingly, a set of Proteobacteria (Campylobacter, Helicobacter and Moraxellaceae) with low signal intensities had a clear negative correlation with PC. Similarly, two Lactobacillus species (L. rhamnosus that was ingested during the trial and L. salivarius) together with another group of Bacillus, Enterococcus spp., also correlated negatively with PC while few Firmicutes exhibited a positive correlation with this lipid species. Among the less abundant lipids (<10% of the total lipid pool), a positive correlation was detected between cholesteryl ester and Collinsella (r = 0.59; Supplemental Table S3, Fig. 5A). On the phylotype level, the correlation was mainly attributable to an uncultured bacterium clone Eldhufec074 within genus Collinsella. Another group of Actinobacteria, namely Actinomycetaceae, showed a significant negative correlation to TG and PA.

Hierarchical clustering and heatmap visualization of the significantly correlated lipid–bacteria pairs revealed groups of lipids and bacteria sharing similar correlation patterns. A constant plaid model biclustering approach was applied to detect coherent groups of significantly correlating microbe–lipid pairs (see Methods). The analysis revealed three major clusters of significantly correlated lipid–microbe pairs (Supplemental Table S2). One of the observed biclusters highlights the above discussed positive association between bacteria related to R. gnavus and TGs. Another TG cluster, containing highly unsaturated long-chain fatty acids, consisted of negatively correlating bacteria related to Megamonas hypermegale (family Veillonellaceae, Clostridium cluster IX) or belonging to the Actinomycetaceae. The third cluster contained a set of Proteobacteria (Helicobacter, Oceanospirillum and Moraxellaceae spp.) and bacteria related to Eubacterium cylindroides within the family Erysipelotrichaceae, which all correlated negatively with ether lipids.

To explore the nature of the identified lipid–microbe pairs, we quantified the enrichment of specific, functionally relevant lipid categories in the biclusters. In particular, we investigated the enrichment of lipids with even and odd number of carbon atoms and the degree of saturation (0, 1, or more double bonds). The origin of lipid can be inferred from its chemical composition: Polyunsaturated fatty acids (PUFAs; more than one double bond) and lipids that have uneven number of carbon atoms originate from plants, bacteria or marine organism as humans cannot synthesize them. Significant enrichment of long acyl chain (p < 0.05; Fisher exact test) was detected among the TG lipids that significantly correlated with bacteria related to R. gnavus (Supplemental Fig. S1; Supplemental Table S2). Such a positive association suggests a role for these R. gnavus-related bacteria in the absorption of dietary lipids. In support of this interpretation, TGs containing fatty acids with odd number of carbons were also relatively common in this group (p = 0.1). In the other two biclusters, no significant enrichment of odd/even carbon count or saturation level was detected but within the bicluster 2 (Supplemental Table S2), Actinomycetaceae correlated exclusively with highly unsaturated PUFAs (Fig. 3).

Associations between the intestinal microbiota and biochemically determined serum lipids

The mean values (SD) for the major serum lipids were total cholesterol 5.10 (1.02), LDL cholesterol 3.00 (1.21), HDL cholesterol 1.50 (0.33) and TG 1.20 (0.71) mmol/L (Kekkonen et al., 2008a). No significant changes were detected in the lipids of the probiotic or the placebo groups during the intervention (Kekkonen et al., 2008b). Hence, we analyzed the potential associations of the enzymatically determined major blood lipids with the intestinal microbiota (Table 2). A positive correlation between bacteria related to R. gnavus and TG was observed (q < 0.01; r = 0.60; Fig. 4B), corroborating the association identified within the global lipid analysis. In line with the spectroscopic lipid analysis, representatives of Bacteroidetes and Uncultured Clostridiales correlated negatively and other implicated Firmicutes positively with enzymatically determined TG (Table 2). Representatives of Proteobacteria and Actinobacteria that correlated negatively with TG in the lipid profiling data did not correlate significantly with enzymatically determined TG. Such inconsistency may partly arise from methodological reasons as the enzymatic assay captures not only TG but also diacylglyceride, monoacylglyceride and free glycerol, while lipid profiling captures TGs at the molecular level.

Collinsella spp. and Eubacterium biforme et rel. showed statistically significant (q < 0.05) and positive correlations to enzymatically determined total and LDL cholesterol (Fig. 5B; Supplemental Table S2). No other significant correlations were identified for total or LDL cholesterol while HDL cholesterol correlated significantly with numerous taxa. Eight different Firmicutes including taxa related to Ruminococcus obeum and D. formicigenerans were found to correlate negatively with HDL, while Uncultured Clostridiales I was the only taxon showing positive correlation to HDL (Supplemental Table S2).