Reprograming of gut microbiome energy metabolism by the FUT2 Crohn's disease risk polymorphism.

Fucosyltransferase 2 (FUT2) is an enzyme that is responsible for the synthesis of the H antigen in body fluids and on the intestinal mucosa. The H antigen is an oligosaccharide moiety that acts as both an attachment site and carbon source for intestinal bacteria. Non-secretors, who are homozygous for the loss-of-function alleles of FUT2 gene (sese), have increased susceptibility to Crohn's disease (CD). To characterize the effect of FUT2 polymorphism on the mucosal ecosystem, we profiled the microbiome, meta-proteome and meta-metabolome of 75 endoscopic lavage samples from the cecum and sigmoid of 39 healthy subjects (12 SeSe, 18 Sese and 9 sese). Imputed metagenomic analysis revealed perturbations of energy metabolism in the microbiome of non-secretor and heterozygote individuals, notably the enrichment of carbohydrate and lipid metabolism, cofactor and vitamin metabolism and glycan biosynthesis and metabolism-related pathways, and the depletion of amino-acid biosynthesis and metabolism. Similar changes were observed in mice bearing the FUT2(-/-) genotype. Metabolomic analysis of human specimens revealed concordant as well as novel changes in the levels of several metabolites. Human metaproteomic analysis indicated that these functional changes were accompanied by sub-clinical levels of inflammation in the local intestinal mucosa. Therefore, the colonic microbiota of non-secretors is altered at both the compositional and functional levels, affecting the host mucosal state and potentially explaining the association of FUT2 genotype and CD susceptibility.

Introduction

The human intestinal microbiome contributes vital biological functions to healthy hosts, including maintenance of immune homeostasis, modulation of intestinal development and enhanced metabolic capabilities (Li ; Turnbaugh ; Lee and Mazmanian, 2010; Qin ; Human Microbiome Project C, 2012). Dysbiosis, which refers to perturbations of the normally stable intestinal microbiota, has been associated with the development and progression of many conditions, including inflammatory bowel diseases (IBDs) (Frank ; Willing ; Lepage ; Morgan ), type 2 diabetes (Qin ) and obesity (Turnbaugh ). The reasons for such associations are not yet clear and may reflect either causal or secondary processes due to the impact on microbial composition and function of inter-individual variability and the contributions of environment and host genetics (Spor ). The contributions of such factors on human microbial composition are beginning to emerge through dietary and environmental studies (Wu ; Claesson ; Morgan ; Smith ), as well as twin studies (Zoetendal ; Stewart ; Turnbaugh ; Smith ). However, much work is still necessary to fully understand the extent of host genetic influence on the composition and function of the gut microbiome and the mechanisms linking these genetic traits with microbial function and disease biology.

A recent genome-wide association study published by our group identified Fucosyltransferase 2 (FUT2) gene as a Crohn's disease (CD) risk locus (McGovern ), a finding that has been validated in a meta-analysis of CD and ulcerative colitis genome-wide association scans (Jostins ). However, the molecular mechanism of the association between non-secretor status and CD remains unknown. Mucin 2 (muc2), the predominantly secreted mucin in the colon, has an important barrier role in intercepting and excluding bacteria from the mucosal cell surface, thereby reducing host susceptibility to colitis (Van der Sluis ; Johansson ; McGuckin ). Research from our group has shown that aberrant glycosylation of muc2 core proteins causes spontaneous colitis in mice (Fu ). Both the core 1- and core 3-derived O-glycans of mucin core proteins are terminally fucosylated, which serve as interceptive binding structures for bacteria (Linden ). Moreover, a subset of human intestinal microbiota produce glycosidases capable of hydrolyzing α-1,2-fucosyl linkages present in various mucin-type glycoproteins, as well as mucus glycan structures that are not capped by fucose (Katayama ). A mass spectrometry-based analysis of insoluble colonic mucin of both Fut2-null and wild-type mice (Hurd ) identified 17 different oligosaccharides with up to eight sugar residues of which 11 were neutral, five sulfated and one sialylated. The most abundant structures were composed of core 2 (Galβ1-3(GlcNAc β1-6)GalNAc-) glycan sequence with some based on core 1 (Galβ1-3GalNAc-) glycan structures. The primary difference in oligosaccharides was the presence of terminal fucose residues forming the blood group H-type epitope in most of the oligosaccharides in wild-type mice. In contrast, all the peaks for oligosaccharides carrying blood group H-type epitopes were absent in the Fut2-null mice. Therefore, FUT2 deficiency may alter the composition of intestinal microbiota by affecting either microbial adhesion and/or utilization of host-derived glycans, potentially leading to dysbiosis.

The phylogenetic composition in non-secretor individuals have been characterized in two studies recently (Rausch ; Wacklin ), showing that the FUT2 genotype was associated both with deviations of overall community ecology and with altered abundances of specific microbes. However, these descriptions did not address the degree to which these alterations were functional nor their potential mechanisms of action in IBD risk. Both questions are of particular interest, because microbial composition can exhibit large inter-individual variations compared with function-based analyses even in healthy individuals (Qin ). This may also be one of the reasons for existing between-study discrepancies (Rausch ; Wacklin ), other than the difference between measurements of the luminal/fecal microbiota and those at the mucosal surface (Eckburg ). As bacterial colonization largely occurs in the outer mucous layers (Johansson ) where the residual glycans that fuel bacterial growth are degraded, lavage sampling of the mucosal surface compartment coupled with functional and metabolic profiling is arguably more biologically relevant to host–microbial glycan metabolism.

We present here a comprehensive description of the mucosal luminal interface of healthy individuals distinguished by secretor status, capturing multiple aspects of the microbial ecosystem including microbiome composition, imputed function, metabolome and proteome. 16S rRNA gene sequencing can be used to characterize the composition and diversity of the microbiota, and with recent advances, it allows us to impute functionality of the microbiome. Deeper insight into microbial functionality can be provided by combining 16S rRNA gene sequencing with proteomic and metabolomic data (Erickson ). We detailed the phylogenetic and functional profiles of the mucosal microbiome associated with FUT2 polymorphism, indicating a strong effect of host genetics on the re-programing of energy metabolisms into dysbiotic setting. The combination of multi'omic analysis also provided us with unprecedented understanding of the dynamics of host–microbial interaction.

Materials and methods

Subject cohort and lavage sample collection

A subject cohort of 39 individuals (Supplementary Table S1) was recruited from patients presenting for screening colonoscopy at Cedars-Sinai Medical Center. Following institutional review board approval, subjects were consented and then included in the study if colonoscopy did not reveal any mucosal abnormalities. Enrolled subjects were prepared for colonoscopy by taking Golytely the day before the procedure. The mucosal lavage samples representing the mucosal luminal interface were collected from cecum and sigmoid colon as described previously (Li ).

Animals

All animal protocols were in accordance with Administrative Panel on Laboratory Animal Care, the Stanford Institutional Animal Care and Use Committee. Conventionally housed Fut2 mice (B6.129 × 1-Fut2tm1Sdo/J; backcrossed 12 generations with C57BL/6J) were re-derived as germ free and maintained in gnotobiotic isolators. Eight-week-old non-littermate germ-free Fut2-deficient Fut2 (n=10), wild-type Fut2 (n=10; C57BL/6J) and heterozygous Fut2 (n=8) mice were colonized with feces obtained from a healthy human donor (secretor) by oral gavage of 200 μl of human fecal sample (male, age 38 years, American diet). The sample was prepared by mixing stored frozen human fecal sample with filter-sterilized pre-reduced phosphate-buffered saline. Mice were singly housed and maintained in gnotobiotic isolators on a strict 12-h light cycle for the experiment. Mice were fed a standard autoclaved mouse diet (Purina LabDiet 5K67, LabDiet, St Louis, MO, USA). Fecal samples were collected 4 weeks after humanization for 16S rRNA gene sequencing using the 454 titanium platform.

Genotyping

The single-nucleotide polymorphism rs601338 (G>A) defines secretor status in Europeans and Africans (Ferrer-Admetlla ). We used rs516246, which is in strong linkage disequilibrium with rs601338, to infer secretor status. Estimate of linkage between rs516246 and rs601338 is 100%. These two single-nucleotide polymorphisms tag each other perfectly, as they are in perfect LD with one another, with R2 of 1.0 and D′ of 1.0, according to Hapmap3 release 2, CEU population. The individuals with the homozygote A/A genotype are defined as non-secretors (Supplementary Materials and Methods). In this cohort, 97% (38/39) of the subjects are Caucasian. One subject is African American, who is heterozygous G/A genotype for rs516246 and therefore categorized as Sese. Mouse genomic DNA was prepared from ear tissue obtained by ear punch. PCR amplification using three primers (F: 5′-CCTGCCATGCTTTCTTTCCTG-3′, R: 5′-ATTCCTTCTCTGACAGGGTTTGG-3′ (WT), 5′-TGGGTAACGCCAGGGTTTTC-3′ (KO)) yielded either a 191-bp band (Fut2) or 154-bp band (Fut2) or both (Fut2).

16S rRNA gene sequencing and microbial composition analysis

Genomic DNA was extracted as previously described (McHardy ). The V4 region of 16S ribosomal RNA genes were amplified and sequenced on an Illumina HiSeq 2000 (Illumina, Inc., San Diego, CA, USA) as previously described (McHardy ). HiSeq reads were processed using QIIME v1.5.0 (Caporaso ) with parameters of: minimum Q-score considered high quality: 20, maximum number of consecutive low-quality base calls allowed before truncating: 3, and maximum number of N characters allowed: 0. All filtered reads had a length of 101 bp. The number of reads per sample ranged from 326 481 to 1 021 473, with a mean of 646 140 and totaling 48 460 491. Sequence sub-sampling was performed for each sample at the depth of 300 000 reads per sample. This normalized data set was used for all the following analysis, including alpha-diversity analysis, beta-diversity analysis and imputed metagenomic analysis. For mouse fecal pellets, after DNA isolation (MoBio fecal DNA kit, MO BIO Laboratories, Inc., Carlsbad, CA, USA), 626 bp amplicons spanning 16S variable regions 3′–5′ (V3–V5) were generated using barcoded forward primer (338F, 906R) (Kashyap ). Samples were sent for pyrosequencing to Duke ISGP (Durham, NC, USA) using the Roche 454 titanium platform (454 Life Sciences, Branford, CT, USA). Operational taxonomic units (OTUs) were picked against the 4 February 2011 version of the Greengenes database (http://greengenes.lbl.gov/cgi-bin/nph-index.cgi) (DeSantis ), pre-filtered at 97% identity. For quality control, all the singletons were removed. After reference-based OTU picking, 97.5% of the total reads were successfully mapped to the reference Greengenes database. These steps were performed using QIIME v1.5.0 (Caporaso ). Alpha rarefaction was performed using the number of observed species, Chao1 and phylogenetic diversity. The comparison of alpha diversity between the two groups at certain sampling depths was performed using a two-sided Student's t-test. Beta diversity of 16S rRNA gene and imputed metagenomic data sets were estimated by computing unweighted UniFrac and Bray–Curtis distances between samples, respectively, using QIIME. Ordination of the resulting distance matrix was performed using principal coordinate analysis (PCoA). Pair-wise comparisons between SeSe, Sese and sese individuals were conducted using the Kruskal–Wallis test to identify differentially abundant bacterial phylotypes at phylum and 97% OTU levels. Multiple hypothesis tests were adjusted to produce a final Benjamini and Hochberg false discovery rate (FDR; Benjamini and Hochberg, 1995), and significant association was considered below a FDR q-value threshold of 0.25.

Identification of functional microbial communities (FMCs)

We used 97% OTUs with a minimum of 30 reads in the data set that were present in>50% of samples. We first defined a co-occurrence similarity measure, which was used to define the network. Assume that the vector x specifies the abundance of the i-th OTUs across the samples, the pair-wise Sparse Correlations for Compositional data (SparCC) ρ was inferred from the abundance profile of each OTU x and x as the measurement of co-occurrence relationship (Friedman and Alm, 2012). A signed weighted adjacency matrix (network) was defined by raising ρ to a power a=(0.5+0.5ρ), with β=4 (Zhang and Horvath, 2005). Modules were defined as branches of a hierarchical clustering tree based on the topological overlap measure. The modules were detected after applying the dynamic tree cut method (Langfelder ). These network modules (clusters) were interpreted as functional microbial communities (FMCs). The classical multidimensional scaling of the topological overlap distance matrix was plotted using the ‘cmdscale' function in R (v.2.15.1). To summarize the profiles of co-occurrence modules, we calculated the eigenOTU, which is defined as the first right-singular vector of the standardized module abundance matrix. The computation is implemented in the R function module Eigengenes (Langfelder and Horvath, 2008). As a weighted average abundance profile, it provides a mathematically optimal way of summarizing the co-occurrence patterns of all OTUs belonging to each module (Tong ). To identify modules (FMCs) that were correlated with clinical traits, we used correlation tests to relate each eigenOTU to the clinical traits. These steps were performed using the WGCNA package (version 1.13) in R (version 2.13.1) (Langfelder and Horvath, 2008). R tutorials explaining the analysis steps can be found on the webpage: http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/ WGCNA/Tutorials/.

Imputation of microbial gene content and metagenomes

This study takes advantage of PICRUSt, a program that infers the metagenome of a sample from its phylogenetic composition and was recently validated against conventional deep-sequencing metagenomics (Langille ). The OTU table was used as the input file for metagenome imputation of individual human and mouse samples. For the metagenomic profiling of FMCs, the OTU table of the six FMCs was generated with one count for each 97% OTU in a given FMC. The gene content of 2590 KEGG (Kyoto Encyclopedia of Genes and Genomes) reference genomes was used to infer the approximate gene content of the detected phylotypes using PICRUSt (v0.1) (http://picrust.sourceforge.net/) (Langille ). The program output the inferred metagenome represented by KEGG Orthology (KO) for each FMC. Taking the PICRUSt KO gene abundance inferences as inputs, the metabolic pathways were re-constructed using HUMAnN (v0.98) (Abubucker ). We restricted our analysis to the KEGG pathways that were present in at least 90% of the samples. Pair-wise Kruskal–Wallis tests between SeSe, Sese and sese individuals were performed to identify imputed KEGG pathways with differential relative abundance. Multiple hypothesis tests were adjusted to produce a final Benjamini and Hochberg FDR (Benjamini and Hochberg, 1995), and significant association was considered below a FDR q-value threshold of 0.05.

Mass spectrometry analysis

For metabolomic analysis, each human lavage samples was subjected to solid-phase extraction to eliminate a polymeric contaminant believed to originate from the lubricant used during colonoscopy preparation. The eluate was dried and reconstituted in 2% acetonitrile in water before mass spectrometry analysis. A 5-μl aliquot of extracted metabolites from each sample was injected onto a reverse-phase 50 × 2.1 mm ACQUITY 1.7-μm C18 column (Waters Corp., Milford, MA, USA) using an ACQUITY UPLC system (Waters Corp.). A Waters Q-TOF Premier was operated in negative-ion (ESI−) or positive-ion (ESI+) electrospray ionization mode with a capillary voltage of 3200 V and a sampling cone voltage of 20 V in negative mode and 35 V in positive mode. Data were acquired in centroid mode with a mass window of 50 to 850 m/z and processed using MassLynx software (Waters Corp.) (Supplementary Materials and Methods). To profile the meta-proteome of lavage samples, matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS) was performed using the soluble fraction of samples as previously described in Li ). The abundances of metabolomic and proteomic peaks were compared using analysis of variance (ANOVA) to identify features associated with FUT2 genotype. Multiple hypothesis tests were adjusted to produce a final Benjamini and Hochberg FDR (Benjamini and Hochberg, 1995), and significant association was considered below a FDR q-value threshold of 0.25. The relatively high FDR was used to avoid excessively strict filtering of metabolomic and proteomic features.

Results

Whole-community microbiome ecology differs according to secretor status

To study the host–microbial interaction at the mucosal luminal interface, 75 lavage samples were collected from the cecum and sigmoid colons of 39 healthy subjects (Supplementary Table S1). We assessed differences in overall microbial ecology between secretors (both homozygous SeSe and heterozygous Sese for the functional allele) and non-secretors (sese).

We first examined the microbial composition in these samples, to affirm that the present cohort matched previously reported differences between secretors and non-secretors in microbial composition (Rausch ; Wacklin ). The microbiota from these samples was profiled by multiplex sequencing, and a total of 47 171 132 reads (628 948±130 744 s.d. reads per sample) were generated after quality control. A total of 4074 OTUs were then identified by grouping reads at a 97% sequence similarity threshold. Compared with SeSe individuals, both Sese and sese individuals exhibited lower alpha-diversity based on the number of observed species (t test, P=0.012 and 0.085, respectively), although the difference between SeSe and sese individuals was not statistically significant (Supplementary Figure S1a). We also measured other diversity indexes, including Chao1 and phylogenetic diversity. Compared with SeSe individuals, Sese individuals exhibited significantly lower alpha-diversity as indicated by Chao1 and phylogenetic diversity indexes at the depth of 300 000 reads per sample (t test, P=0.019 and 0.02, respectively). The same trend was observed in sese individuals, although not statistically significant (t test, P=0.10 for Chao1 and 0.18 for phylogenetic diversity) (Supplementary Figures S1b and c).

The beta-diversity measured by unweighted UniFrac distance matrix was calculated for each sample to evaluate the similarity between microbial communities. PCoA demonstrated that the phylogenetic compositions of SeSe microbiomes were significantly different from those of Sese (Adonis test, P=0.016) but not of sese individuals (Adonis test, P=0.092) (Supplementary Figure S2a). The significant difference in phylotype abundances reported previously, namely, the increase in Bacteroidetes among non-secretors, was confirmed at the phylum level (Rausch ) (Supplementary Figure S2b).

To analyze at lower taxonomic levels, we filtered out low-abundant 97% OTUs based on the criteria of (1) minimum total observation count of 30 across all samples and (2) being observed in at least 60% of the samples, reducing the number of OTUs from 4074 to 419. Among these OTUs, 19 (4.5%) of them were depleted in Sese and sese compared with SeSe individuals (Kruskal–Wallis, FDR q<0.25) (Supplementary Data set S1). In summary, the FUT2 polymorphism was significantly associated with selected phylotypes of colonic microbiota in Sese and sese individuals, and the alterations in Sese individuals resulted in a significant shift of microbial composition compared with SeSe. These data revealed the gardening effect of FUT2 polymorphism on phylogenetic composition of the colonic microbiota.

Non-secretor-associated functional changes revealed by imputed metagenomes

We hypothesized that these compositional changes result in selectively augmented or deficient functional capabilities that may be relevant to CD susceptibility. To test this idea, we inferred the metabolic capacities of mucosal microbiota associated with secretor status, using a recently developed bioinformatic pipeline centering on the PICRUSt (Langille ) and HUMAnN tools (Abubucker ). In the ‘gene content inference' step, the gene contents and 16S rRNA gene copy number of the detected phylotype were predicted based on its evolutionary similarity with the 1119 KEGG reference genomes. In the subsequent ‘metagenome inference' step, the resulting gene content predictions for all microbial taxa with the relative abundance of 16S rRNA genes in each samples are combined and corrected for expected 16S rRNA gene copy number, to generate the expected abundances of gene families in the entire community represented by KOs. The prediction accuracy of PICRUSt has been validated using human and mammalian gut microbiome with paired 16S rRNA gene and metagenome sequencing data (Langille ). The relative abundances of KEGG pathways in each sample were then reconstructed by mapping KOs to these pathways using HUMAnN. At our routine sampling depth, both Sese and sese individuals harbored 15% fewer microbial genes on average than SeSe individuals (Figure 1a), which is consistent with the significant lower compositional diversity observed in Sese individuals compared with SeSe individuals.

Figure 1

Imputed metagenomes reveal the significant enrichment of KEGG pathways in secretors and non-secretor individuals. (a) Distribution of bacterial genes in SeSe, Sese and sese individuals. The proportions of individuals having a given number of genes are shown. (b) Communities clustered using PCoA of the Bray–Curtis distance matrix. Each colored point corresponds to a sample. The clustering of SeSe was significant compared with both Sese and sese individuals (Adonis test, P=0.004 and 0.004, respectively). (c) Relative abundance of KEGG metabolic pathways in microbiome samples was colored by secretor status. Only the 13 pathways showing concordant alterations in both human and murine data sets were plotted.

The similarity of functional states of the microbiomes from secretors and non-secretors was evaluated by the composition of imputed metagenomes. PCoA using Bray–Curtis distance demonstrated separation of the samples from SeSe, Sese and sese individuals along PC1. The clustering of SeSe was significant compared with both Sese and sese individuals (Adonis test, P=0.004 and 0.004, respectively), suggesting that variations at the metagenomic level were more profound than those at the compositional level and that FUT2 exhibited haploinsufficiency in programing the metagenomic functions (Figure 1b). In support of the distinction between the phylogenetic and imputed metagenomic data sets, Procrustes analysis of the Bray–Curtis PCoA plots derived from 16S rRNA gene and imputed metagenome data sets showed that the clustering of samples across data sets was not significant (P=0.510; Supplementary Figure S3). Among the 154 imputed metabolic pathways, 23 (14.9%) were differentially abundant between SeSe and Sese individuals. The alterations in sese individuals were greater, as shown by the changes of abundances of 43 (27.9%) of the imputed pathways compared with SeSe individuals (Kruskal–Wallis, FDR q<0.05) (Supplementary Data set S2, Supplementary Figures S4 and S5).

The FUT2-associated changes were more robust at the imputed metagenomic versus the phylotypic level. As compared with SeSe individuals, a diverse consortium of metabolic functions represented by 27 KEGG pathways were depleted in sese individuals, including amino-acid metabolism pathways, cofactors and vitamins metabolism pathways and genetic information processing pathways (Supplementary Figure S6). A broad-based decrease in amino-acid biosynthesis was observed, including lower abundances of lysine (KO00300), valine, leucine and isoleucine (KO00290), and phenylalanine, tyrosine and tryptophan (KO00400) biosynthesis pathways. Accompanying the depletion, 16 microbial pathways were enriched in these individuals, highlighted by aspects of energy metabolism, including carbohydrate and lipid metabolism, cofactors and vitamins metabolism and glycan biosynthesis and metabolism (Supplementary Figure S7). These data suggest that FUT2 gene polymorphism acted in a haploinsufficient manner to perturb metabolic pathways such as amino-acid biosynthesis encoded by the gut microbiome at the mucosal interface.

To validate these findings, we performed the same analysis on a 16S rRNA gene data set of the fecal samples collected from humanized FUT2−/− mice (germ-free FUT2−/− mice colonized with human feces from a healthy secretor) (Kashyap ). Among the 146 metabolic pathways reconstructed, 47 (32.2%) of them were differentially abundant between the FUT2+/+ and FUT2−/− mice (Supplementary Data set S3), comparable to our findings with the human samples. After cross comparing the two data sets, we identified 13 pathways that were consistently enriched or depleted with FUT2 haploinsufficiency (Kruskal–Wallis, FDR q<0.05) (Figure 1c). Specifically, the carbohydrate and lipid metabolism and glycan biosynthesis-related pathways were over-represented in Sese/sese individuals and FUT2+//FUT2−/− mice, whereas the relative abundances of five amino acid and vitamin metabolism-related pathways were enriched in SeSe individuals and FUT2+/+ mice (Figure 1c). These findings suggest that FUT2 genotype had a similar impact on imputed microbial metagenomic functions in humans and mice, notably in reduced amino-acid synthesis capabilities.

FMCs associated with non-secretor status

It is well understood that even healthy individuals differ remarkably in their fecal and mucosal surface gut microbial composition, especially at the genus and species level. Mucosal microbiota can be further assessed for functional relatedness based on their co-occurrence patterns (Faust ; Koren ). To determine whether such ecological structures can be observed in this data set, we developed a methodology to infer microbial co-occurrence networks (see Materials and methods). Nodes (OTUs) of these networks were grouped based on their topological overlaps using hierarchical clustering. Using this approach, six modules, ranging from 39 to 97 OTUs, were identified (Supplementary Figure S8). These modules of OTUs represent putatively key ecological features, which we term as functional microbial communities (FMCs). Multi-dimensional scaling was used to depict module structure and network connections (Figure 2a). Phylogenetically related OTUs were clustered into the same FMC, presumably because preferences for ecological niche are more likely to be shared between more closely related microbes. However, each FMC also included phylogenetically distinct OTUs from different phyla, suggesting that in addition to phylogenetic relatedness the formation of FMC depended upon additional ecological affinities, which could range from syntrophic dependences to convergent functionality between distinct phylotypes (Supplementary Data set S4).

Figure 2

FMCs associated with non-secretor status. (a) Classical multi-dimensional scaling plot in which OTUs in each FMC represented by colored dots tend to form distinct clusters. (b) FMC-trait correlations and P values. Each cell reports the Pearson correlation coefficient (and P-value) derived from correlating FMC eigenvectors (rows) to traits (columns). For the association with non-secretor status, the SeSe and Sese individuals were grouped together as secretor. For the association with FUT2 genotype (rs516246), additive genetic model was used. The table was color-coded by correlation according to the color legend.

The ‘abundance' of each FMC can be quantified by defining the OTU abundance profiles of each FMC (its module eigenvector; defined in Material and methods). One can thus correlate the abundance of the FMCs with the metadata, including host genotype, disease phenotype, age and so on. When using an additive genetic model, the abundances of turquoise and blue FMCs significantly associated with the copy number of the FUT2 loss-of-function allele reciprocally (P=0.04 and 0.05, respectively, with rs516246 by Pearson correlation) (Figure 2b). This observation was concordant with results at the individual OTU level: when examining the membership of the FMCs, we found that 12 of the 19 OTUs enriched in SeSe individuals were assigned to the blue FMC. Gender was another subject phenotype that significantly associated with microbial composition: the turquoise (P=0.004, Pearson correlation) FMC was enriched in females, whereas the blue and red FMCs (P=0.04 and 0.006, respectively, Pearson correlation) were more abundant in males. The gender effect on intestinal microbiome has also been reported previously in humans and murine model (Mueller ; Markle ). Because of the difference of male/female ratio in SeSe, Sese and sese individuals (Supplementary Table S1), the gender effect could potentially contribute to the association of turquoise and blue FMCs with FUT2 genotype.

To determine whether these co-occurring microbial communities represented distinct functional units at the mucosal surface, we profiled the metabolic capabilities of FMCs using the approximate gene contents imputed previously. After aggregating the individual inferred genomes according to module membership, the relative abundances of metabolic pathways in each FMC were re-constructed. The functional profiles of FMCs were highly variable (Figure 3). The pathways associated with FUT2 clustered into two groups that were over-represented in the FUT2 loss-of-function allele-associated turquoise FMC and functional allele-associated blue FMC, respectively. Moreover, the pathways enriched in SeSe individuals tended to cluster into amino-acid metabolism class and cofactors and vitamins metabolism class, whereas the pathways enriched in Sese and sese individuals highlighted amino-acid metabolism, lipid metabolism and biosynthesis of secondary metabolites (Figure 3). These metabolically specialized microbial communities were therefore responsible for the imputed metagenomic alterations associated with FUT2 polymorphism.

Figure 3

Variations of KEGG metabolic pathways in the functional microbial communities. The heatmap shows the functional profiles of FMCs (columns) based on the relative abundance of FUT2-associated metabolic pathways (rows) after z-score transformation. The color bar on top shows module membership. The dendrograms show the hierarchical clustering of columns and rows, respectively, using Euclidean distance. The two pie-charts show the number of pathways in each functional class for the cluster associated with turquoise and blue FMC, respectively.

Non-secretor-associated metagenomic changes reflected by metabolomic and proteomic profiling

To determine whether the imputed metagenomic alterations associated with FUT2 polymorphism correlate with changes in metabolic activities of mucosal microbiota, we profiled the soluble metabolites of the same lavage samples using Q-TOF MS. The analysis generated a rich metabolomic data set consisting of 649 and 576 spectral features in the cecum and sigmoid regions, respectively. Putative IDs were assigned to 372 ions by comparing their m/z values to those available in online databanks using a predefined mass error window of 20 p.p.m. The putative IDs were then used to map out the ions to various metabolomic pathways in the KEGG data set. In accord with our previous study,∼50% of all metabolites were located at the terminal end of metabolic pathways, suggesting enrichment for end-products (McHardy ). In the cecum, 48 metabolites were mapped to the 13 KEGG pathways associated with non-secretor individuals in the human and murine data sets and were present in>90% of the samples. Of these, 13 (27.1%) were differentially abundant among secretors (SeSe and Sese) and non-secretors (sese) (ANOVA, FDR q<0.25). In the sigmoid, 30 metabolites were mapped to the non-secretor-associated pathways; 4 (13.3%) were differentially abundant (ANOVA, FDR q<0.25; Figure 4a and Supplementary Data set S5). A less stringent q-value, up to 0.25, was used to avoid missing significant associations, as shown in recent comparable study designs (Wu ; Morgan ). When using more stringent q-value threshold (FDR q-value<0.1), these metabolites did not show significant association, which is comparable with the results from study with similar design (Wu ). Thus, differences in imputed microbial metagenome content corresponded to abundances of metabolic end-products directly detected in the same samples.

Figure 4

Meta-metabolomic and meta-proteomic features that differentiate secretors and non-secretors. Relative abundance of meta-metabolomic (a) and meta-proteomic (b) features in lavage samples is colored by secretor status.

To determine how the host reacted to the changes of the functional state of the microbiota, we also profiled the proteomic features of the same samples using MALDI-TOF MS. We focused our analysis on peaks of human origin. Of the 453 peaks included, the abundances of 16 (3.5%) were significantly different among secretors and non-secretors (ANOVA, FDR q<0.25). Although only 17 (3.8%) molecular features could be identified, we found that the expression levels of human neutrophil peptides 1 and 2 (HNP-1 and -2) were significantly higher in non-secretors, consistent with a sub-clinical inflammatory state at the mucosal surface (Figure 4b). This difference could be driven by one particular sese sample that had high expression of HNP-1 and -2 relative to the other samples (Figure 4b). To exclude such possibility, we repeated the analysis excluding this sample, and found that the levels of HNP1 and HNP2 were still significantly higher in non-secretors based on nominal P-values (ANOVA, nominal P=0.007 and 0.047, respectively), although FDR q-values were>0.25 (FDR q=0.27 and 0.49, respectively). Our clinical records did not indicate that this individual had gastrointestinal symptoms or chronic inflammation in the intestine at the time of sampling. These data suggest that the non-secretor state of the mucosa, via alteration of the mucosal surface ecosystem, changes the inflammatory state of the human intestinal mucosa.

Discussion

We combined 16S rRNA gene sequencing, metagenome imputation, meta-metabolomic and meta-proteomic profiling to delineate the integrated landscape of the mucosal surface ecosystem. The phylogenetic diversity and composition of intestinal mucosal microbiota in non-secretor individuals were significantly different from that of secretors. Compared with SeSe, the metabolic functions encoded and expressed by the gut microbiome in non-secretors were enriched for carbohydrate and lipid metabolism, cofactors and vitamins metabolism and glycan biosynthesis and metabolism and depleted for seven pathways related to amino-acid metabolism. These alterations in humans were highly consistent with analogous changes in the murine genetic counterpart. Changes in the imputed metagenomes were reflected by concordant metabolite pools as determined by meta-metabolomic assays, providing validation for certain imputed metagenomic changes as functionally consequential. Moreover, these microbial functional changes were accompanied by sub-clinical intestinal inflammation. FUT2 therefore appears to have a role in shaping the functional state of the mucosal surface by affecting not only microbial composition but also the resulting functional state of the microbiota at the human intestinal mucosal surface.

It was surprising that difference with the SeSe group were in several cases (for example, alpha-diversity) greater in the Sese rather than in the sese group. This raises two issues. One is the extent of glycan difference produced by haploinsufficiency. As there was a strong haploinsufficiency phenotype in all the facets of this study, we surmise that a substantial glycan change is produced by haploinsufficiency. However, to our knowledge, the glycan profile in heterozygous individuals has not been well described in the literature. The other issue is why the Sese group had a larger and more significant difference than the sese group. One possible explanation is that the inter-individual variation of gut microbial phylogenetic composition is inherently large (Human Microbiome Project C, 2012; Yatsunenko ). In this context, it is notable that the sizes of the test groups were modest (9 for sese, 18 for Sese). It is possible that the particularly small size of the sese group made it prone to outliers and potentially underpowered for establishing mean phenotypes and robust statistical comparisons (for example, the 95% confidence interval for alpha-diversity was larger for the sese (4.4) than for the Sese (1.8) group).

The imputed metagenome represents an accurate but approximate inference of the reference microbial genomes currently available. Potential bias may result from unmappable 16S rRNA gene reads and lack of sufficient reference genomes. After the reference-based OTU picking, 97.5% of the total reads were successfully mapped to the reference Greengenes database. The 2.5% unmappable 16S rRNA gene reads might cause the loss of metagenomic content that can be captured by shot-gun sequencing. Also, only the 2590 microbial genomes that had identifiers in the Greengenes reference tree were used as the reference to predict unknown genomes. Despite these bias and limitations, PICRUSt predictions usually reach high agreement with metagenomically measured gene content (Spearman r=0.8–0.9) (Langille ). In this study, the findings of imputed metagenomic analysis in human subjects were validated by metabolomic and proteomic data as well as by comparison to an independent 16S rRNA gene data set from humanized FUT2 mice. In the future, adding meta-transcriptomic data will further enrich our understanding of microbial functional capability at any given time.

Inconsistencies of FUT2-associated imputed metagenomic changes were also observed between human and murine data sets. In the humanized mouse gut microbiota, there are changes of microbial compositions as compared with the donor (Turnbaugh ). The fecal sample used for humanization was from only one healthy donor, which might cause inherent bias due to the limited sample size. Also, taxonomic biases between the two different 16S rRNA gene data sets may exist due to different PCR primer sequences, amplicon lengths and sequencing technologies (Claesson ).

One of the key drivers of gut microbiota composition and function is the type and quantity of complex carbohydrates, which are typically derived from either diet or host mucus (Wu ; Koropatkin ). These polysaccharides serve as a primary metabolic input for the abundant carbohydrate-fermenting bacteria within the microbiota (Martens ). However, different microbes within the gut are differentially endowed with abilities to use specific types of glycans (Sonnenburg ), and so differences in carbohydrate availability, such as the presence or absence of fucose in mucosal glycans, translate into selective and regulatory events that result in discrete alterations in the microbiota's functional properties (Kashyap ). Individual members of the microbiota can alter gene expression to accommodate the absence of fucose, while other members may be lost, and others recruited (Kashyap ). Additionally a change in carbohydrate utilization by relatively few members can cascade into ecosystem-wide alterations, given the interconnectedness of metabolic functions within the microbiome. Fucose processing by a gut-resident symbiont allows expansion of pathogenic species, such as Salmonella typhimurium. Similarly, increased sialic acid release by certain bacterial species (Bacteroides thetaiotaomicron) allows expansion of Clostridium difficile (Ng ). Differences in host glycan fucosylation result in distinct microbial ecosystems at the mucosal interface, the compositions and metabolic activities of which are likely precursors and predictors of ensuing disease phenotypes. This concept has been validated by the upstream role of microbial fucose processing for the expansion of the enteric pathogen S. typhimurium (Ng ).

Among the phylotypes that are depleted in Sese and sese compared with SeSe individuals (Supplementary Data set S1), Roseburia and Faecalibacterium are both short-chain fatty acid-producing bacteria (Scheppach, 1994; Wong ) and reported to be anti-inflammatory (Sokol ; Atarashi ). Moreover, depletion of Firmicutes and expansion of Proteobacteria members are also characteristic of changes associated with the IBD microbiome (Frank ; Sokol and Seksik, 2010). Both Sese and sese individuals harbored 15% fewer microbial genes on average than SeSe individuals (Figure 1a). The lower metagenome diversity is an unfavorable feature for the host, which has also been observed in IBD (Qin ) and obesity (Le Chatelier ). Among the pathways that were consistently depleted or enriched in Sese and sese individuals, the increase in glutathione metabolism and decrease in amino-acid biosynthesis (particularly lysine) pathways have been reported as a feature of the metagenome in IBD patients (Morgan ). These data indicated that the non-secretor-associated changes of microbial composition and imputed metagenomic functions are also characteristics of other chronic inflammatory conditions and therefore are unfavorable for the host.

FUT2−/− mice have a marked alteration in gastric mucosa glycosylation, characterized by diminished expression of alpha(1,2)fucosylated structures (Magalhaes ). As gut microbes have developed the ability to degrade host-derived glycans (Katayama ; Sonnenburg ), the deprivation of terminal fucosylation may affect the metabolic activity of the gut microbiota and thus its fermentation products potentially available to the host. Recent work reported that enterohemorrhagic Escherichia coli encodes a two-component system, termed FusKR, which responds to fucose and controls metabolic gene expression (Pacheco ). The imputed metagenomic changes in non-secretors highlighted the depletion of indispensable amino-acid biosynthesis. This group of metagenomic functions complements that encoded by the host genome (Metges, 2000; Gill ; Qin ). In the IBD metagenome, amino-acid biosynthesis and carbohydrate metabolism are reduced in favor of nutrient uptake (Morgan ). Such changes might reflect compensation by the microbiota for the lower availability of carbon sources. Amino-acid starvation can lead to host stress response and the induction of autophagy of intestinal epithelial cells (Tattoli ), which may increase the risk of IBD. Although the current imputed metagenomic analysis is limited to the KEGG pathway level, further insights could be gained by extending the analysis to individual KEGG module or enzyme. It would be helpful to further identify the individual genes or reactions that are differentially abundant in these pathways, which would serve as potential candidates for therapeutic manipulation.

Low richness of gut microbiota is a well-known feature of patients with IBD (Manichanh ; Lepage ), and other chronic conditions such as obesity (Turnbaugh ), and of elderly patients with inflammation (Claesson ). A recent study defined two groups of individuals who differed by the number of gut microbial gene as low gene count (LGC) and high gene count (HGC) (Le Chatelier ). LGC individuals exhibited an imbalance of pro- and anti-inflammatory bacterial species and evidence of low-grade inflammation. We have shown that both Sese and sese individuals harbored 15% fewer microbial genes on average than SeSe individuals and therefore exhibited low metagenomic richness. Similarly, in Sese and sese individuals, genera associated with LGC individuals, including Faecalibacterium and Coprococcus, were also more dominant. Moreover, the sub-clinical inflammatory state at the mucosal surface was reflected by the higher expression levels of HNP-1 and -2 in non-secretors.

The data presented here supports the hypothesis that the FUT2 loss-of-function allele increased the risk of CD by shaping the functional states of mucosal microbiota. Meta-analysis of genome-wide association studies has increased the number of confirmed IBD (both CD and ulcerative colitis) susceptibility loci to 167 (Jostins ), indicating that IBD is biologically heterogeneous. The analysis presented in this study focused on individuals without clinical symptoms. It will be important to extend the same analysis to patients with CD to determine to what extent the changes we present are recapitulated in the disease setting. In addition to FUT2, other risk genes have also been shown to affect the gut microbial composition, such as NOD2 (Petnicki-Ocwieja ; Frank ) and several defensin genes (Ivanov ; Salzman ). It is currently unclear whether the microbiota associated with genes of similar functions has the same compositional and functional signatures. The stratification of gut microbiome by host genetics is a crucial step for elucidating the pathogenic mechanism of IBD as well as the design of personalized therapeutic interventions.

To achieve unprecedented understanding of the ecological structures and biomolecular activities of the gut microbiome, it is necessary to extend the analysis to multiple levels of biological organization—genome content, gene expression, protein expression and metabolism (Raes and Bork, 2008; Nicholson ). In this study, we used multiple 'omic approaches to disentangle the complex host–microbial metabolic interplay. Meta-proteomic analysis in this case focused on the host side, but metabolomics could be extended to incorporate richer microbial analysis in the future. Although only a limited number of integrative ‘omics profiles of the gut microbiota currently exist (Erickson ; Perez-Cobas ), multi'omic studies have shown great potential in providing a holistic picture of the metabolic status of the gut microbiota and the host response to functional changes.