Effects of Sample Handling Method on Individual Meta’omic Features.

To assess the contribution of individual features (e.g., specific microbial species or genes) to this strong global agreement, we performed two-way ANOVA tests on each metagenomic and metatranscriptomic feature, normalizing abundance data (Methods) to partition feature variance across the eight subjects and three sample handling methods (Fig. 2 D–F). Only features exceeding a minimum relative abundance of 10⁻⁴ (0.01%) in at least 3 of the 24 samples were considered. Relative to between-subject variation, no individual microbial species demonstrated statistically significant variation across sample collection methods after correction for multiple hypothesis testing (Benjamini–Hochberg α = 0.05) (18) (Fig. 2D). Similarly, sample handling method was not observed to have a statistically significant effect on the relative abundance of any individual genes (Fig. 2E). These findings are consistent with the strong within-subject, between-method agreements observed for DNA-level species and gene relative-abundance profiles in the correlation analyses described above (Fig. 2 A and B), and further suggest that—in addition to strong global agreement—individual metagenomic measurements are robust to subject-collected stool sample handling methods.

The effect of sample shipping on metatranscriptomics was comparably small, with only a very small minority (n = 84, <5% of total) of transcripts exhibiting statistically significant variation across sample handling methods (Fig. 2F and Dataset S1). The nature of these differentially abundant transcripts was consistent with a pattern of live cells responding to an altered environment via changes in gene regulation. For example, up-regulated genes in ethanol-fixed samples were largely involved in oxidative metabolic processes, a signal consistent with bacteria responding to a combination of oxygen exposure and a new carbon source. On the other hand, a subset of genes up-regulated in the RNAlater-fixed samples were involved in response to osmotic stress [e.g., the glycine betaine/proline transport system (19)], which is consistent with the high saline content of RNAlater solution. It is also possible that some transcripts experienced a sample handling method effect because of variation in RNA stability across the three storage conditions.

Comparison of Oral-Gut Microbial Ecology in the HPFS and Human Micobiome Project Cohorts.

After establishing the data quality of subject-collected samples, we next sought to meta-analyze the relationship between the oral and gut microbial communities in our own HPFS cohort and the larger healthy population of the Human Microbiome Project (HMP) (20) (Fig. 3A). In addition to characterizing potential microbial transit along the gastrointestinal tract, this contextualized our eight subjects within a broader population, as the HPFS participants were healthy but of both restricted geography (Boston metropolitan area) and age range (over 65). We identified a subset of 62 commonly occurring species in the eight frozen HPFS saliva and stool samples and compared their abundance profiles with 69 oral (tongue dorsum) and 81 stool samples from the HMP. As expected, differences in body-site specific ecology proved to be the largest effect in both cohorts, with HPFS and HMP stool samples forming a single, well-mixed cluster and HPFS saliva samples associated with but distinct from HMP oral metagenomes from the buccal mucosa and tongue dorsum (SI Appendix, Fig. S3). Whereas the compositions of the oral versus gut samples were largely distinct, we did observe a small number of species that occurred regularly at both body sites (Fig. 3A).

Open in a separate window

Fig. 3.

Oral-gut ecology in the HPFS and HMP cohorts. (A) We isolated species observed in the eight pairs of frozen stool and saliva samples from the HPFS cohort with relative abundance of at least 10⁻² (1%) in two HPFS samples. The taxonomic profiles of these species were compared with stool and tongue samples from the HMP cohort, with tongue representing the oral community. Samples were clustered by Bray–Curtis distance, and species were clustered by rank correlation. Note that samples cluster strongly by body site (oral vs. gut). Highly abundant oral species are more likely to be detected at low levels in the gut. Green numbers associate oral-gut co-occurring species with detailed abundance profiles in B and C. (B) Eight abundant oral species detected in the HPFS saliva samples were detectable at low abundance in the stool samples from the same individuals, but showed minimal transcriptional activity in the stool. Gray lines connect oral DNA (light blue), gut DNA (dark blue), and gut RNA (red) from the same individual. (C) D. invisus is an unusual example of a gut-dominant species that also occurs in the oral cavity.

Detection of Oral Bacterial Strains in the Gut Microbiota.

Although the oral and gut environments are anatomically linked, the degree of exchange between their resident microbiota is not completely understood (21). Bacterial species from the oral community are carried along with food into the stomach, but the degree to which they survive or remain biologically active in the lower gastrointestinal tract has not been systematically characterized, particularly whether oral microbes contribute to the stable commensal gut community as measured by the stool. We first examined our frozen saliva and stool samples for cases of bacterial species co-occurring in the oral and gut communities of each subject. We defined a species as co-occurring if, for at least two subjects, the species occurred with relative abundance greater than 10⁻² (1%) in a subject’s saliva sample and greater than 10⁻⁵ (0.001%) in the same subject’s frozen stool sample. Of 33 species meeting the first criterion (common oral species), eight met our criteria for detection in the stool: four members of the Streptococcus genus (Streptococcus salivarius, Streptococcus parasanguinis, Streptococcus australis, and Streptococcus sanguinis), along with Haemophilus parainfluenzae, Veillonella atypica, Veillonella parvula, and Actinomyces odontolyticus (Fig. 3B). For each of these species, the typical drop in relative abundance between the oral and gut communities was one-to-two orders-of-magnitude, with higher oral abundance generally corresponding to higher gut abundance. This finding suggests that, although DNA from these oral species does survive transit to the gut, it does not form a dominant component of that community.

If a small fraction of oral microbiome members do survive transit to the healthy gut, the gut metatranscriptome provides one means of assessing their active biological contribution to this community, albeit agnostic of their transcriptional activities at the oral site. Our metagenomic samples were assessed, among other methods, by mapping DNA reads to sets of species-specific marker genes using MetaPhlAn (16) (Methods). We applied the same procedure to RNA read data, in which case MetaPhlAn provides a profile of a subset of the metatranscriptome that can be unambiguously assigned to individual species. Based on this analysis, we found that when oral species were detectable in the gut at the DNA level, they rarely appeared to be transcriptionally active there. It was only in those cases where an oral species achieved its highest DNA abundance in the gut that its species-specific RNA was detected, and typically with relative abundance one-to-two orders-of-magnitude lower than the corresponding DNA abundance (some cases of RNA nondetection may correspond to exceptionally small, but still nonzero, RNA abundances that fell below our detection limit). As oral species’ abundances in the gut community were already multiple orders-of-magnitude lower than in the oral microbiome, this indicates that oral microbes that do survive transit to the gut are not stable, active contributors to its ecology.

Although this evidence is consistent with a pattern of abundant oral species passaging through the gut at low levels, co-occurrence could also be explained by separate strains of the same species adapted to the oral and gut environments. To evaluate whether oral and gut signatures of a species represent the same or different strains, we compared profiles of species-specific marker-gene presence and absence across subjects’ stool and saliva samples. These profiles capture gene gain and loss events in specific strains of a species and were used here as molecular “barcodes” for identifying or differentiating strains (Methods). Barcodes for common oral species often differed markedly across the eight saliva samples, indicating the presence of subject-specific strain variation (SI Appendix, Figs. S4–S11). However, in cases of within-subject oral-gut species co-occurrence, we rarely detected markers for a species in a subject’s stool sample that were not also seen in the subject’s saliva sample. This finding supports the hypothesis that oral-gut species co-occurrence is driven not by separate pools of niche-adapted strains of the same species, but instead by oral strains surviving transit to the gut in low quantities.

An even smaller number of abundant gut microbes occurred at appreciable levels in the oral community, the only significant example being Dialister invisus; this was true for samples from both the HMP and HPFS cohorts (Fig. 3C). There were three cases of D. invisus oral-gut co-occurrence among the eight HPFS subjects, and in each case the same strain was carried with high abundance in the stool and low abundance in the saliva (SI Appendix, Fig. S12). Curiously, despite its high DNA abundance in the gut, D. invisus made almost no contribution to the pool of species-specific transcripts, an indicator of reduced transcriptional activity (Fig. 3C and SI Appendix, Fig. S13). Of the remaining five subjects, four carried D. invisus exclusively in the oral community. This evidence resolves previously disparate ecologies for D. invisus, which was isolated from the human oral cavity (22) but also identifed as a marker for human stool (23), suggesting that it is atypically capable of persisting at high abundance in both the oral and gut communities and may freely transit between the two.

Relating Microbial Genes in the Gut at the DNA and RNA Levels.

We next investigated possible global models for metagenome vs. metatranscriptome regulation in the gut microbiota. Among the host-adapted microbes in this community, DNA and RNA abundances would be correlated if many genes were not differentially regulated and were transcribed at the same constant rate. This would be the case, for example, if typical gut microbe molecular activity was regulated by genome modifications over evolutionary time, as opposed to transcriptional regulation on a more rapid time scale. To test this theory, we quantified the relative abundances of genes and transcripts using HUMAnN (17) in the HPFS stool samples, with the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthogroups (KO) database as a reference set for gene families (Methods). Across all samples, a total of 3,292 KOs were observed with relative abundance of at least 10⁻⁴ (0.01%) at either the DNA or RNA level. Averaged over the eight frozen stool samples, we found that gene abundance and corresponding transcript abundance were well correlated (Spearman’s r = 0.76; two-tailed P 0.001) (Fig. 4).

Open in a separate window

Fig. 4.

Up- and down-regulated pathways and clades in the gut metatranscriptome. Gene and transcript relative abundances are generally well correlated (Spearman’s r = 0.76). (A–H) Each scatterplot illustrates the average gene (DNA) and transcript (RNA) relative abundance for 3,292 KOs from the eight frozen HPFS stool samples, highlighting a prominent over- or underexpressed functional module. Red circles correspond to KOs where RNA > DNA; blue circles correspond to KOs where DNA > RNA. Marks on the x or y axis margins represent KOs with zero measured abundance in one dataset but nonzero abundance in the other. The trends illustrated here were all of large effect (fold-change > 2) and statistically significant following FDR correction (Methods).

To identify differentially regulated transcripts, we computed the log RNA/DNA abundance ratio for each gene on a subject-by-subject basis and then tested whether the mean of the eight log ratios was significantly different from zero following false-discovery rate (FDR) correction (indicating a pattern of consistent over- or underexpression) (Dataset S2). Although a substantial fraction of KOs were not consistently differentially regulated in the gut (nonsignificant fold-change or fold-change < 2; n = 1,340, 41%), we also observed many transcripts whose relative abundances were an order-of-magnitude higher or lower than expected from the DNA abundance of their corresponding gene families (significant fold-change > 10; n = 724, 22%) (Fig. 4). For example, tetA, which encodes a transporter protein conferring tetracycline resistance, was on average 1,000-times more abundant at the RNA level than the DNA level, one of the strongest such effects (one-sample t test, two-tailed P < 10⁻⁵) (Fig. 4A). We then used these rankings of gene-level expression as input to a functional enrichment analysis, searching for KEGG pathways and modules whose member genes (KOs) were enriched for consistent over- or underexpression (Fig. 4 and Dataset S3). All of the over- and underexpression relationships discussed in the following sections were both of large effect (fold-change > 2) and statistically significant following FDR correction (Methods).

Microbial genes encoding ribosomal proteins were among the most strongly overexpressed (Fig. 4B). Note that these are, of course, distinct from the ribosomal rRNAs depleted from metatranscriptomic assays (Methods). Indeed, we observed distinct clusters of ribosomal protein-coding gene overexpression across three domains of life, with bacterial ribosomal genes having the highest overall abundance, followed by archaeal ribosomal genes, and finally eukaryotic ribosomal genes detectable at the low end of DNA relative abundance. Notably, these archaeal ribosomal genes occurred as part of a “burst” of other highly expressed archaea-associated functions, including methanogenesis (Fig. 4C) and the archaeal RNA polymerase. This signal can be explained predominantly by the presence of the archaeon Methanobrevibacter smithii in five of the eight HPFS subjects. In these five subjects, the relative abundance of M. smithii at the DNA level ranged from 0.005 to 0.053 (0.5–5.3%), whereas its relative contribution to the pool of species-specific transcripts ranged from 0.021 to 0.147 (2.1–14.7%) (SI Appendix, Fig. S13). This finding confirms that M. smithii tends to be highly transcriptionally active relative to other species in the gut (9), an activity associated with its unique functional role in maintaining gut energy balance (24) by removal of hydrogen and correspondingly increased rates of fermentation among its bacterial neighbors (25). The second carbon oxidation step of the citrate cycle was also found to be highly overexpressed (Fig. 4D). Given that the vast majority of gut-dwelling bacteria are anaerobic (26), the observed overexpression of the citrate cycle is best explained by its role in replenishing biochemical intermediates, rather than its role in aerobic metabolism. This finding is consistent with the observation that only genes involved in one part of the cycle (reactions yielding oxaloacetate) were significantly overexpressed and emphasizes the need to avoid the blind application of broad pathway names to specific biochemical functions observed in the microbiome.

Conversely, other metabolic routes were consistently underexpressed, including the metabolism of starch and sucrose (Fig. 4E). Moreover, many pathways involved in the biosynthesis of small metabolites tended to be expressed at relatively low levels (Fig. 4F), consistent with earlier findings (8). Because these compounds will have high bioavailability in the gut because of the host’s diet, it is logical that it would, in most cases, be more energetically favorable to import them rather than synthesizing them directly. Genes involved in sporulation were also strongly underexpressed (Fig. 4G). Sporulation is a stress-response mechanism used by Gram-positive bacteria to survive in nonideal conditions (e.g., nutrient depletion) and also plays a role in certain bacterial infections [e.g., Clostridium difficile (27)]. However, neither of these conditions is applicable to this cohort’s healthy human gut environments, with a presumably undisturbed baseline spore fraction, and hence it is not surprising that the sporulation pathway is underexpressed. Finally, we also observed peptidoglycan biosynthesis to be significantly underexpressed in the gut microbiome (Fig. 4H). Peptidoglycans are critical components of bacterial cells walls and have roles in the maintenance of cell shape and integrity (28). To perform those roles, peptidoglycan biosynthesis must be carefully regulated to coordinate with cellular growth (29). Consistent with this idea, we found that peptigoglycan biosynthesis transcription and ribosomal protein transcription were strongly correlated across species (Spearman’s r = 0.78; two-tailed P 0.001) (SI Appendix, SI Methods and Fig. S14). However, within species, ribosomal protein-coding genes tended to be overexpressed relative to their metagenomic abundance, whereas peptidoglycan biosynthesis tended to be underexpressed (consistent with the community-level trends seen in Fig. 4). This pattern may result from the fact that the end-products of ribosomal protein transcription are physical building blocks of the ribosome and need to be continuously synthesized, whereas the end-products of peptidoglycan biosynthesis transcription are enzymes that act on the building blocks of the cell wall (and may thus persist at smaller concentrations).

The Metatranscriptome Varies Significantly Among Subjects Despite Metagenomic Concordance.

Recent large-scale studies of the human microbiome have shown that, although taxonomic composition at a particular body site can differ dramatically from one individual to the next, functional composition tends to be much more conserved (20, 30). Indeed, we observed similar patterns in our metagenomic (DNA) data from the HPFS cohort. Although there was considerable variation in the 10 most abundant genera across the eight subjects (Fig. 5A), the 10 most abundant functional categories were very well conserved (Fig. 5B). However, when we examined the 10 most abundant functional categories as measured from metatranscriptomic (RNA) data, we found that they were noticeably more variable than their DNA-based counterparts, although still more conserved than genus-level taxonomic composition (Fig. 5C).

Open in a separate window

Fig. 5.

Functional diversity at the transcriptional level suggests a pattern of subject-specific metagenome regulation. (A) The 10 most abundant genera in each subject. (B) The 10 most abundant Enzyme Commission (EC) level-3 function categories in each individual, based on metagenomic measurements. (C) The 10 most abundant EC level-3 function categories in each individual, based on metatranscriptomic measurements. (D) Comparison of within-sample evenness for taxonomic and functional profiles across individuals using the Pielou metric. (E) Comparison of between-sample (β) diversity for taxonomic and functional profiles using the Bray–Curtis metric. Two-tailed P-values are based on the Wilcoxon signed-rank test. Although taxonomic profiles are highly variable among subjects, metagenomes are stable but differentially regulated at the transcriptional level.

To more rigorously quantify these observations, we computed within-subject evenness and between-subject dissimilarity for the taxonomic, DNA-level function, and RNA-level function composition for the eight frozen stool samples. Taxonomic composition was significantly less even than either DNA- or RNA-level functional composition, with a small number of (not always the same) genera dominating each individual’s gut microbiome (Wilcoxon signed-rank test, two-tailed P = 0.01) (Fig. 5D). Moreover, in line with our qualitative observation above, we found that between-subject variation in RNA-level functional composition was significantly stronger than between-subject variation in DNA-level functional composition (P < 10⁻⁵) (Fig. 5E). Thus, although subjects share a highly conserved gut metagenome derived from microbial contributors that are quite variable among subjects, its constituent functions are fine-tuned at the transcriptional level in a subject-specific manner.

HPFS Cohort and Subject Recruitment.

We recruited the eight individuals sampled in this study from the Health Professionals Follow-up Study, a prospective cohort of 51,529 United States male dentists (58%), veterinarians (20%), pharmacists (8%), optometrists (7%), osteopathic physicians (4%), and podiatrists (3%), aged 40–75 y at enrollment, who returned a mailed health questionnaire in 1986 (www.hsph.harvard.edu/hpfs). Every 2 y, participants have completed comprehensive questionnaires with updated information on their lifestyle and medical history, with a follow-up rate exceeding 90%. In the spring of 2010, 210 active HPFS participants with current addresses in the greater Boston area were mailed a short questionnaire to assess their interest in providing stool and saliva samples for a pilot investigation of diet and gastrointestinal microbiota. We received positive responses from 75% of participants. From among these participants, we selected 63 men who represented extremes of Western versus prudent dietary intake and sent them kits for self-collection of saliva and stool at home. Subject recruitment and protocols were approved by the Harvard School of Public Health Institutional Review Board.

Sample Collection.

Saliva samples were collected using the self-collection kit Oragene-DNA OG-500 (DNA Genotek) and accompanying instructions. For the stool collection, we took advantage of an established protocol well-received in the HMP cohort, providing us with starting sample material comparable to previous studies (14). Stool was self-collected by deposition in a commode specimen collection system (catalog 02–544-208, Fisher Scientific). Both the saliva and stool samples were placed by the subject in a styrofoam box with frozen icepacks for pick-up by a member of our study team and delivery to our laboratory within 24 h of sample collection. We independently evaluated the effects of this initial transport protocol on stool metagenomes and metatranscriptomes by sequencing DNA and RNA from aliquots of a fresh stool sample that had been (i) immediately flash-frozen versus (ii) incubated on ice for 24 h (SI Appendix, SI Methods and Fig. S1). Participants also completed a brief questionnaire to report their recent diet and other pertinent exposures that may affect specimen composition, including recent consumption of antibiotics or probiotics or the use of oral hygiene products. The sample collection period ran from February to July 2011, during which time we collected samples from 25 participants. From among these participants, we randomly selected complete samples from eight men for this analysis.

Laboratory Procedures.

Upon receipt at the laboratory, saliva tubes were frozen to −80 °C. The stool samples were homogenized using a spatula and separated into three aliquots (fresh-frozen control, mock-shipped RNAlater-fixed, and mock-shipped ethanol-fixed). The homogenization step was intended to reduce between-aliquot variability for the purposes of our methods evaluation; subjects using the self-collection and self-shipping protocol validated here would not be required to homogenize their stool. The first 100-mg subsample (control) was flash-frozen on dry ice and then stored at −80 °C until being thawed for extraction. The control aliquot represented the baseline composition of the sample upon arrival at the laboratory, following the initial HMP-based sample collection and transport protocol described above. To a second 100-mg subsample we added 700 μL of RNAlater (Ambion) and then allowed the mixture to fix at −4 °C overnight. A third 1-g subsample was collected using a feces tube with spatula (catalog 80.623.022, Sarstedt), to which we added 5 mL of 95% (wt/wt) ethanol and allowed the mixture to fix at room temperature overnight. The stool masses described above were acquired by taking a single scoop with a corresponding collection spatula, just as subjects would do when following the proposed self-collection protocol. Ethanol- and RNAlater-fixed subsamples were then stored at −80 °C until all samples were ready for the mock-shipping stage of the experiment. During mock-shipping, all ethanol- and RNA-later fixed samples were transferred to a cardboard box and subjected to natural environmental condition changes (at ambient temperatures) for 48 h to mimic conditions during carrier or mail transport. After this phase, mock-shipped samples were then stored at −80 °C until extraction; the frozen control and RNAlater-fixed subsamples entered the extraction pipeline directly, whereas the ethanol-fixed subsamples were first sliced to yield 100 mg of stool.

DNA and RNA Extraction and Sequencing.

For saliva samples, DNA was extracted using procedures established for the Human Microbiome Project (14). Stool samples in RNAlater and ethanol were centrifuged at maximum speed for 10 min and the supernatant removed. DNA and RNA were extracted and processed following the protocol established in ref. 35. Briefly, bacterial lysis buffer containing lysozyme and proteinase K was added to ~100 mg of stool sample, incubated for 10 min, followed by 3 min of bead beating to extract the nucleic acids. Qiagen AllPrep DNA spin columns were used for RNA and DNA isolation following the manufacturer’s protocol (Qiagen). Illumina DNA-seq libraries were generated as previously described (36) with the following modifications. For each sample, 100 ng of genomic DNA was sheared to 150–300 bp in size using a Covaris LE220 instrument (Covaris) with the following parameters: temperature: 7–9 °C; duty cycle: 20%; intensity: 5; cycles per burst: 200; time: 90 s; shearing tubes: Crimp-Cap microTUBES with AFA fibers (Covaris). DNA fragments were end-repaired, 3′-adenylated, ligated with indexed Illumina sequencing adapter, and PCR-enriched, as previously described (37). The resulting Illumina fragment sequencing libraries were size-selected to contain inserts of 180 bp ± 3% in length using a Pippen Prep system (Sage Science) following the manufacturer’s recommendations.

For RNA-seq libraries, 5 μg of initial RNA was depleted for ribosomal RNA using Ribo-Zero (Epicentre), subjected to DNase treatment to deplete remaining sample DNA, fragmented, and then used as a template for strand-specific cDNA synthesis by dUTP marking and degradation of second strand cDNA (35). This procedure has been known to introduce 1–2% Escherichia coli genomic DNA into the final cDNA library (a result of E. coli-derived DNA polymerase I and ligase being used in the cDNA generation steps). Including versus excluding E. coli sequences in downstream bioinformatic analyses did not affect the conclusions of this work. RIN scores for all metatranscriptomic samples are provided in SI Appendix, Table S1. The average RIN score over all samples was 6.9, whereas the averages for the control, EtOH-fixed, RNAlater-fixed batches were 5.7, 7.0, and 7.8, respectively.

All DNA-based samples were then sequenced by Illumina HiSeq, with each DNA sample run on a single sequencing lane, and each RNA (cDNA) sample run on a pair of sequencing lanes. Raw sequences were then filtered to remove human sequence contamination, quality trimmed, and depleted of duplicate reads also using HMP protocols (14). Sequences from paired RNA sample lanes were merged before downstream analyses. Fig. 1 reports average sequence statistics for all raw and postprocessed samples. Sequence data are available for download from the Sequence Read Archive with accession no. PRJNA188481.

Taxonomic and Functional Profiling.

We performed taxonomic profiling with MetaPhlAn (16). Briefly, MetaPhlAn operates by mapping raw sequence reads to a database of predefined clade-specific marker genes. Markers are those genes occurring in isolates from a particular clade but not outside of that clade. After mapping reads to clade-specific marker genes, the resulting raw counts are normalized for total marker-gene length and outliers, yielding profiles of (i) clade relative abundance, (ii) marker-gene presence/absence, and (iii) marker gene abundance (in RPKM units). Because of gene gain and loss events, an individual strain will not necessarily carry all of the markers associated with its corresponding species. A specific pattern of marker presence and absence can therefore be used as a molecular “barcode” to identify a strain across samples. Examples of MetaPhlAn-based strain barcodes are illustrated in SI Appendix, Figs. S4–S12. We defined “commonly occurring species” to be those detected with relative abundance of at least 10⁻² (1%) in at least two of our eight subjects.

When applying MetaPhlAn to RNA read data, the output is a profile of relative marker gene transcription across clades (focusing on species-level measurements in this work). A species’ value in this profile depends on two factors: (i) its overall genomic abundance in the sample, and (ii) the average transcriptional level of its unique marker genes. Notably, only the first criterion matters when quantifying species abundance at the DNA level, as each marker gene is expected to occur at equal frequency after normalization. If a species has higher RNA-based abundance than DNA-based abundance, this suggests that its marker genes are more actively transcribed than marker genes from other species, which we can interpret as an approximate measure of the species’ clade-specific transcriptional activity. Notably, this measure can only approximate a species’ overall transcriptional activity, as the species’ unique marker genes may be more or less actively transcribed than other pathways (which includes many nonunique genes; e.g., ribosomal RNAs).

We performed functional profiling of genes and transcripts using HUMAnN (17) based on the KEGG database of gene families and pathways (31). This procedure involved first mapping quality- and human-filtered DNA and RNA reads from a meta’omic sample to protein-coding sequences from the KEGG Orthology using USEARCH (38). We then used HUMAnN to collapse the output of this mapping step to produce a table of relative abundances for all KOs. When the input to USEARCH and HUMAnN are DNA read data, the resulting functional profile describes the relative abundance of gene families in the sample; when RNA read data are used as the input, the profile instead quantifies the corresponding transcripts of these gene families. If a gene family or its corresponding transcript is never observed with relative abundance above 10⁻⁴ (0.01%) across all samples, then it is merged into an “other” category and excluded from subsequent statistical analysis.

Samples were additionally profiled by mapping reads against a collection of reference genomes (SI Appendix, SI Methods). Raw mapping results for the eight control stool metagenomes and metatranscriptomes are provided as Dataset S6.

We thank the participants in the Health Professionals Follow-up Study who graciously participated in this research; Elizabeth Frost-Hawes, Jill Arnold, and Siobhan Saint-Surin at the Harvard School of Public Health, who assisted with recruitment for this study; and James Bocchichio, Doyle Ward, Kimberly Kanki, Scott Mahan, and Erin Euler at the Broad Institute for their assistance in project logistics, sample handling, and nucleic acid extractions. This work was supported by National Institutes of Health (NIH) National Institute of Allergy and Infectious Diseases Contract HHSN272200900018C (to the Broad Institute); Grant CA166150 (to J.I.); Grants NIH R01 CA137178, P30 DK043351, and P50 CA127003, and an American Gastroenterological Association Gut Microbiome Pilot Research award (to A.T.C.); Danone Research PLF-5972-GD (to W.S.G.); NIH Grant R01 HG005969, National Science Foundation Grant DBI-1053486, Army Research Office Grant W911NF-11-1-0473, and a Center for the Study of Inflammatory Bowel Disease Pilot Study grant (to C.H.); and in part by the generous support of the Eric Larson and Barbara Wu Trust. A.T.C. is a Damon Runyon Clinical Investigator.