Impact of demographics on human gut microbial diversity in a US Midwest population

Biobank/consent

Subjects were selected from the Mayo Clinic Biobank, an institutional resource composed of biological specimens accompanied by clinical data obtained from patient medical records and patient-provided risk factor data (Olson et al., 2013). Written informed consent was provided by all individuals enrolled. The study was approved by the Institutional Review Board of the Mayo Clinic (#13-003694).

Subject selection/questionnaire

As described below, representative fecal samples were chosen according to sex, age, race, body mass index (BMI), alcohol use, and tobacco use. A survey accompanied the invitation to participate in the study and inquired about bowel symptoms, such as diarrhea, constipation, blood in the stool, and nausea or vomiting; current medication use; current supplement use; antibiotic use in the past 2 weeks; current smoking status; current alcohol intake; current weight; and cancer history.

To define a sample truly representative of the Midwestern United States (Minnesota, Iowa, North Dakota, South Dakota, and Wisconsin), we estimated the distributions of selected behavioral risk factors among subjects included in the Center for Disease Control’s 2011 Behavioral Risk Factor Surveillance Survey (BRFSS; http://www.cdc.gov/brfss/annual_data/annual_2011.htm). Mayo Clinic Biobank participants were then selected based upon BRFSS distributions of age, sex, race, BMI, smoking status, and alcohol intake. Two rounds of subject selection were employed for this study. The first round took place in the summer of 2013. During this round, 380 subjects (ages 20–49 years) were selected. In the summer of 2014, a second group of 267 patients (ages 50–79 years) was selected. Stool samples were collected from 25% of the samples who consented in the first round and 31% of the subjects in the second round.

The final, combined sample was chosen on the basis of age (20–29, 30–39, 40–49, 50–59, 60–69, and 70–79 year groups), sex (10 males and 10 females in each age group), race (white or non-white), BMI (<30 kg/m² or ≥30 kg/m²), smoking status (current or former/never), and alcohol use (yes or no). The proportions of all inclusion criteria were chosen to reflect the underlying population in the Upper Midwestern states, as estimated by data from the BRFSS. Exclusion criteria, assessed based on survey information completed at the time of enrollment into the study, included use of medications, non-GI cancers, use of any antibiotic within the past 2 weeks, or any known GI disease or symptoms including, but not limited to, inflammatory bowel disease, irritable bowel syndrome, esophageal cancer, intestinal cancer, constipation, diarrhea, Clostridium difficile infection, and celiac disease/sprue.

Collection method

The collection process was reviewed in detail with study participants, and a stool collection kit was provided to each subject. Stool samples were collected by the subjects and returned to Mayo Clinic Rochester within 24 hours of passing. Specimens were stored at −80 °C until DNA extraction.

Sample preparation and sequencing

Fecal DNA was extracted using the PowerSoil kit (MoBio, Carlsbad, CA, USA) according to the manufacturer’s instructions. Genomic DNA was used as a template for the polymerase chain reaction (PCR), with 0.3 μM V3–V5 barcoded primers (Caporaso et al., 2012) targeting 357F and 926R of the bacterial 16S gene (5′AATGATACGGCGACCACCGAGATCTACACTATGGTAATTGTCCTACGGGAGGCAGCAG3′ and 5′CAAGCAGAAGACGGCATACGAGATNNNNNNNNNNNNAGTCAGTCAGCCCCGTCAATTCMTTTRAGT3′, respectively). PCR conditions were as follows: 95 °C/3 min; 35 cycles of 98 °C/30 s, 70 °C/15 s, and 72 °C/15 s; and finally 72 °C/5 min in a Bio-Rad T100 Thermocycler (Bio-Rad, Hercules, CA, USA) using Kapa Hotstart Hi-Fi DNA polymerase (Kapa Biosystems, Boston, MA, USA). PCR product sizes were verified using the Agilent TapeStation with reaction cleanup, and DNA was purified using an epMotion automated system (Eppendorf, Hauppauge, NY, USA) with the Agencourt AMPure PCR Purification System (Beckman Coulter, Brea, CA). Final quantitation was performed using a QuBit HS dsDNA kit and the QuBit 2.0 fluorimeter (Life Technologies, Carlsbad, CA, USA). Samples were pooled to equal concentration and sequenced on one lane of a MiSeq at the Mayo Genomics Facility using the MiSeq Reagent Kit v2 (2 × 250 reads, 500 cycles; Illumina Inc., San Diego, CA, USA). Pre-processed sequence files were then processed via the IM-TORNADO bioinformatics pipeline with the default parameter settings to form operational taxonomic units (OTUs) (Jeraldo et al., 2014). IM-TORNADO uses paired-end reads to form OTUs, and was shown to be more sensitive than methods using single-end reads based on synthetic mock community studies.

Statistical analyses

Statistical analysis proceeded in two steps. First, overall associations between demographic and health behavior–related variables and the microbiota were investigated. Next, specific associations at the level of taxa were investigated.

To perform overall association tests, we summarized microbiota data using both α-diversity and β-diversity. α-diversity reflects species richness and evenness within bacterial populations. Two α-diversity metrics, the observed OTU number and the Shannon index, were investigated. The observed OTU number reflects species richness, whereas the Shannon index measures both species richness and evenness. Before calculating the α-diversity metrics, we first replaced the observed singleton OTU count with a more robust estimate to reduce the influence of sequencing errors (Chiu & Chao, 2015). We then computed the α-diversity estimates at the minimum sample coverage among all samples (0.997) to standardize by sample completeness (“estimateD” in the R package iNEXT) (Chao et al., 2014). β-diversity reflects the shared diversity between bacterial populations in terms of ecological distance; different distance metrics provide distinctive views of community structure. Two β-diversity measures, unweighted and weighted UniFrac distances, were calculated using the OTU table and a phylogenetic tree (“GUniFrac” function in the R package GUniFrac) (Lozupone & Knight, 2005; Chen et al., 2012). The unweighted UniFrac reflects differences in community membership (i.e., the presence or absence of an OTU), whereas the weighted UniFrac captures this information and also differences in abundance. Rarefaction was performed on the OTU table before calculating UniFrac distances.

To test the overall association based on these diversity measures, we used a regression model: ${\displaystyle \text{Microbiota}\sim \text{Batch}+\text{X}+\text{X:Batch},}$ where microbiota is the outcome variable, summarized by the diversity measures described above, and X is the variable of interest. In the model, age and BMI were treated as continuous variables while the other variables as categorical. We included a batch variable to account for potential batch effects, since sequencing was performed in two batches, and a covariate-Batch (X:Batch) interaction term. The interaction term allowed different association strengths for the two batches. We assessed the overall effect (main + interaction effect), which is equivalent to testing the null hypothesis of no association in either batch. To assess the association with α-diversity measures, we performed regular linear regression analysis with the likelihood ratio test of regression coefficients, as the outcome was approximately normal. To assess the association with β-diversity measures, we used the PERMANOVA procedure (“adonis” function in the R package vegan), which is a multivariate analysis of variance based on distance matrices and permutation (McArdle & Anderson, 2001). Significance was assessed by 1,000 permutations. Ordination plots were generated using principal coordinate analysis (PCoA) on unweighted UniFrac-based distances as implemented in R (“cmdscale” function in the standard R package).

A distance-based coefficient of determination, R² (i.e., the percentage of overall microbiota variability explained by a variable), was used to quantify the effect size of the overall association (McArdle & Anderson, 2001). Distance-based R² is defined as ${\text{R}}^2=\frac{\text{tr}\left(\text{HGH}\right)}{\text{tr}\left(\text{G}\right)},$ where tr(.) is the trace of a matrix, H is the projection matrix into the column space spanned by the corresponding variable, and G is the Gower’s centered matrix, which is defined as $\text{G}=\left(\text{I}-\frac{{11}^{\text{T}}}{\text{n}}\right)\text{A}\left(\text{I}-\frac{{11}^{\text{T}}}{\text{n}}\right),$ where I is an identity matrix, 1 is a vector of 1s, and $\text{A}={({\text{a}}_{\text{ij}})}_{\text{nxn}}={(-{\text{d}}_{\text{ij}}^2/2)}_{\text{nxn}}$ is a matrix constructed using the pairwise distances d_ij. In the formula, the total variability of the microbiota was summarized using a distance metric. To quantify the overall effect size well, a good distance metric should capture the association signals. In this study, we used the unweighted UniFrac distance, as most of the associations investigated were significant only for this metric. To account for potential overestimation due to a small sample size, we calculated an adjusted R², which is defined as ${\text{R}}_{\text{adj}}^2=1-\frac{(1-{\text{R}}^2)(\text{n}-1)}{\text{n}-\text{p}-1},$ where n is the sample size and p is the degree of freedom of the covariate.

Next, we assessed associations between demographic and health behavior–related factors and the relative abundance of taxa. To do so, we square-root transformed relative abundance data and used a linear model adjusting for batches where appropriate. To address the non-normality of the taxa data, significance was assessed with 1,000 permutations; the F-statistic was the test statistic. False discovery rate (FDR) control was performed based on the Benjamini-Hochberg procedure to correct for multiple testing (“p.adjust” in R). To reduce the number of tests, we confined the analysis to taxa with prevalences greater than 10% and median nonzero proportions greater than 0.05%. An FDR-adjusted P-value (or Q-value) of less than 10% was considered significant. All statistical analyses were performed in R 3.0.2 (R Development Core Team, Vienna, Austria).

Gut microbiota profile of the Midwestern United States

As revealed by a comparison with BRFSS data, the distribution of important demographic variables within the MWRP cohort, which consists of 118 subjects drawn from the Mayo Clinic Biobank, generally reflects the demographic characteristics of the Midwestern United States, although alcohol users are slightly over-represented in subjects older than 50 years of age (Table 1).

Stool samples from this cohort were collected and deep sequenced. 16S rDNA–targeted sequencing yielded 110,997 reads/sample on average (range: 33,579–385,720). Clustering of these 16S sequence tags produced 1,745 non-singleton OTUs at a 97% similarity level. The distribution of OTU abundance and prevalence was typical, dominated by rare and low-abundance OTUs (Figs. 1A and 1B). An average of 402 OTUs was detected in each subject (range: 186–810). The median OTU abundance was 163 read counts per OTU (range: 2–1,456,439), and the median OTU prevalence was 11% (range: 0.8%–100%). Only 61 OTUs (3.5% of total OTUs) occurred in more than 95% of samples, indicating a relatively small core set of OTUs. The detected OTUs were classified into 11 phyla, 80 families, and 171 genera using the Ribosomal Database Project (RDP) classifier. At the phylum level, Firmicutes accounted for 48.3% of total reads and Bacteroides for 48.1% of total reads. At the family level, the most abundant families were Bacteroidaceae (33.76%), Ruminococcaceae (22.94%), and Lachnospiraceae (15.13%). At the genus level, the dominant genera were Bacteroides (33.76%), Ruminococcus (9.83%), and Faecalibacterium (9.13%). Consistent with previous studies, we observed large intersubject variability in taxa abundance (Figs. 1C and 1D) (Human Microbiome Project Consortium, 2012; Flores et al., 2014). Table 2 lists the “core” taxa (prevalence greater than 95%) identified in this study; we identified 4 core phyla, 18 core families, and 22 core genera. The overall pattern in this cohort was similar to that in the HMP, although specific values were very different (Huse et al., 2012; Human Microbiome Project Consortium, 2012). This was probably due to biological variability, including differences in genetics, demographics, and health behavior, as well as technical variability, including differences in sample collection, preparation, sequencing, and bioinformatics processing.

Demographic and health behavior–related variables shaping the gut microbiota

Because of the wide demographic and health behavior–related variation it captures, the MWRP data set provides a good opportunity to identify associations between these variables and the microbiota. We first performed overall association tests based on α- and β-diversity measures. α- and β-diversity measures provide a holistic view of the microbiota; distinctive yet interrelated, these measures focus on different aspects of microbiota structure. In this study, we chose the following α-diversity measures: the number of observed OTUs (after rarefaction) as a species richness measure, and the Shannon index as an overall diversity measure incorporating both species richness and abundance. In terms of β-diversity, which describes overall microbiota structure, we chose the phylogeny-based unweighted and weighted UniFrac distance metrics; the unweighted UniFrac focuses on community membership, whereas the weighted UniFrac reflects both membership and abundance (Chen et al., 2012).

The overall association results are summarized in Tables 3 and 4. Increased BMI was associated with decreased species richness (P = 0.017, Fig. 2A). BMI was also significantly associated with overall microbiota structure, as revealed by both unweighted and weighted UniFrac analysis (P = 0.006 and 0.030, Table 4). Race was associated with both species richness (P = 0.013, Fig. 2B) and overall microbiota structure (P = 0.003, unweighted UniFrac), with white subjects exhibiting greater species richness. Sex and alcohol use were both associated with overall structure (P = 0.049 and 0.041 respectively, unweighted UniFrac, Figs. 2C and 2D). In addition, alcohol use was associated with increased species richness (P = 0.040, Table 3). Some evidence of association of tobacco use with overall microbiota structure was observed (P = 0.091, unweighted UniFrac). Age, by contrast, was not significantly associated with measures of α- or β-diversity. Statistical significance was primarily observed for unweighted measures, such as species richness and the unweighted UniFrac distance, indicating that these demographic and health behavior–related factors chiefly affect community membership or the rare biosphere of the microbiota (Chen et al., 2012).

Next, we quantified the effect sizes of these variables with adjusted, distance-based R² values (see Methods). Using the unweighted UniFrac distance to summarize overall microbiota variability, the adjusted R² values were 0.95% for race, 0.88% for BMI, 0.53% for alcohol use, 0.46% for sex, 0.33% for tobacco use, and 0.0% for age (Fig. 3), indicating huge intersubject variability relative to the variability associated with these factors.

Microbial signature of demographic and health behavior–related factors

We next set out to identify the microbial signature of the demographic and behavioral factors that were most significantly associated with overall microbiota structure (BMI, sex, race, tobacco use, and alcohol use). We used permutation-based univariate association tests to identify associated bacterial taxa at the phylum, family, and genus level. FDR control was used to correct for multiple testing.

Consistent with overall association tests, a number of bacterial taxa were associated with BMI at an FDR of 10% (Fig. 4, Table 5). In particular, we observed an increase in Firmicutes and a decrease in Bacteroidetes in obese subjects, consistent with findings from some previous studies (Turnbaugh et al., 2009; Walters, Xu & Knight, 2014). In addition, we detected an increase in the genus Eubacterium, consistent with two previous studies that reported associations between Eubacterium rectale and obesity (Walters, Xu & Knight, 2014). Interestingly, we also found an increase in Fusobacteria, which have been shown to be associated with colon cancer (Rubinstein et al., 2014).

Race, sex, and alcohol use were also associated with changes in the abundance of specific taxa (Fig. 4, Table 5); here, we focus only on the most statistically significant associations (Q-value < 0.05). Relative to non-whites, whites had higher levels of Holdemania and lower levels of Megamonas (Table 5). Interestingly, although Holdemania is one of the 30 most abundant genera in the gut microbiome of individuals from high-income countries (Arumugam et al., 2011), it was absent in a recent study of children from Bangladesh. It is not clear why race might be associated with the abundance of this common but not universal bacterial genus. Megamonas species are associated with colon cancer (Weir et al., 2013), but they are also associated with normal glucose tolerance (vs type 2 diabetes) (Zhang et al., 2013). Again, the connection between Megamonas and race is unclear. Finally, current alcohol users had lower levels of the family Rikenellaceae, which is also depleted in end-stage liver disease (Bajaj et al., 2014), and higher levels of Christensenellaceae (Table 5), a bacterial family enriched in individuals with low BMI (Goodrich et al., 2014). As is true for race, the mechanisms linking sex and alcohol use with changes in these taxa are not yet clear. We were unable to identify specific taxa associated with tobacco use after FDR control.