Total microbial load of the duodenum across patients with GI symptoms is log-normally distributed

A digital PCR-based determination of total microbial load [9, 29] from 250 human duodenal aspirates revealed samples that spanned loads from our detection limit of ~ 5 × 10³ rRNA gene copies/mL up to nearly 10⁹ copies/mL. The overall distribution of total loads was log-normal with mean = 6.13 Log₁₀ copies/mL and standard deviation = 1.12 Log₁₀ copies/mL (Fig. ​(Fig.1A,1A, B). A quantile-quantile (QQ) plot was constructed to compare the sample distribution to a log-normal distribution (Fig. ​(Fig.1B).1B). Data from our samples aligning with the y = x line on a QQ-plot indicate a high similarity between the sample distribution and a theoretical log-normal distribution [32]. Neither age nor gender significantly correlated with total microbial load (Fig. S1). Total microbial load also did not correlate with patient reported intake of probiotics supplements or yogurts, smoking, or usage of proton pump inhibitors (Fig. S2, Table S2). Current antibiotic usage appeared to lower the average total microbial load, but antibiotic usage in the previous 6 months had no impact (Fig. S2, Table S2).

Open in a separate window

Fig. 1

Microbial load distribution across 250 human duodenal aspirate samples. A Histogram of the total microbial load in 250 duodenal aspirate samples overlaid with a kernel-density estimate. B Quantile-quantile plot comparing the sample distribution of the log₁₀-transformed total microbial load in duodenal aspirate samples to a normal distribution. C Kernel-density estimate plots showing the absolute abundance distribution for the taxa with greater than 50% prevalence in duodenal aspirates. Prevalence (defined as a taxon’s frequency of occurrence in our dataset) and number of samples with each genus are labeled next to the distribution. A legend indicates strict anaerobes (red line through O₂) and the location each genus is commonly found (saliva and/or stool) [30, 31]. Classification of taxa as common in stool or saliva was determined by prevalence of ≥ 50% (stool data are not included in this study) in the 16 participants for whom we had paired samples

Digital PCR anchored 16S rRNA gene amplicon sequencing [9] (hereafter quantitative sequencing) provided absolute taxon abundances in each sample and a statistical framework for differentiation between real and contaminant taxa (Methods). We first compared the culture counts from aerobic (MacConkey agar) and anaerobic (blood agar) plates to the total load of microbes expected to grow on these plates (Fig. S3). For aerobic plating, we observed a bimodal distribution of combined Escherichia-Shigella, Enterobacteriaceae, Enterococcus, and Aeromonas bacterial load from quantitative sequencing and culture and a high correlation between the two measurements (Spearman, 0.61, P < 0.001, N = 244). For anaerobic plating, we observed lower concordance (Spearman, 0.35, P < 0.001, N = 244) between quantitative sequencing and culture. This discrepancy could reflect the difficulty in culturing many intestinal microbes [33], especially anaerobes that are initially collected and processed in aerobic environments.

Next, we analyzed the log-transformed absolute-abundance distributions for the most prevalent genera in our dataset (Fig. ​(Fig.1C).1C). We define prevalence as a taxon’s frequency of occurrence in our dataset. Streptococcus was present in all 250 samples and followed an approximately log-normal distribution with a mean load that was half an order of magnitude below that of the mean total microbial load and an equal standard deviation. Other genera showed wide-ranging distributions that deviated from normality. For example, Porphyromonas appears bimodal with two local maxima whereas Haemophilus exhibits a long tail towards higher microbial loads. The 23 most prevalent genera in this study are also commonly found in the oral microbiota [30]. A subset of these genera (Streptococcus, Veillonella, Prevotella 7, Haemophilus) are also commonly found in stool samples, indicating possible survival of these genera throughout the entire GI tract [31]. The majority of prevalent genera are either strict or facultative anaerobes, indicating that parts of the duodenal environment are likely anoxic in this patient population.

Direct transmission of microbes from saliva to duodenum

To investigate whether many of the taxa found in the duodenum originated from the oral cavity we analyzed a subset of 21 patients for whom we had paired saliva and duodenum samples that were collected during the same hospital visit. Digital PCR revealed that the total microbial load in saliva was roughly 2.5 orders of magnitude higher than the total load in the duodenum (Kruskal-Wallis, P < 0.001).

Further, the range in saliva total loads was 3 orders of magnitude smaller than the range in total loads of the duodenum samples (Fig. ​(Fig.2A).2A). No significant correlation was observed between the total microbial loads in paired saliva and duodenum samples (Fig. ​(Fig.2B).2B). In this study, all samples were collected with a custom double-sheathed catheter via endoscope (see “Methods” section) that moves beyond the outer sheath before aspirating duodenal fluid. This custom catheter should limit oral microbiota contamination of the duodenum during the procedure. Additionally, the optimized sample-processing protocol (see “Methods” section) should eliminate extracellular DNA from swallowed dead bacteria.

Open in a separate window

Fig. 2

Relationship between saliva and duodenal aspirate microbiomes. A Total microbial load of 21 paired duodenal aspirate and saliva samples. B No significant correlation between the total microbial load of 21 paired duodenal aspirate and saliva samples. C Percentage of taxa in duodenal aspirate samples also present in paired (same patient) vs the average of all non-paired saliva samples (Kruskal-Wallis, P < 0.001). D Volcano plot showing the ratio of relative abundances of species in duodenum vs saliva samples. The red dashed line indicates a significance threshold at q = 0.1 (Kruskal-Wallis with Benjamini-Hochberg correction). Undefined Streptococcus sp. classified as S. pneumoniae with 80% confidence and one base pair mismatch to common Streptococcus taxon found in all samples

To evaluate the direct transmission of microbes from saliva to duodenum, we compared the shared taxa between paired (same patient) and randomly paired samples from the same dataset. On average, 89% (± 6% S.D.) of the taxa in the duodenum were also found in the paired saliva sample, whereas only 66% (± 9% S.D.) were found in the average of all non-paired comparisons (Fig. ​(Fig.2C,2C, Kruskal-Wallis, P < 0.001), suggesting direct transmission of oral taxa to the duodenum. We then looked for genera that were proportionally enriched in either saliva or duodenum samples. Campylobacter was present in 21/21 saliva samples but only 10/21 duodenum samples. The absence of Campylobacter in about half of the paired duodenum samples indicates the oral cavity may be the preferred niche of Campylobacter or that Campylobacter has a high sensitivity to the antibacterial properties of the stomach and small intestine [34] (Fig. ​(Fig.2D).2D). In contrast, an undefined species of Streptococcus was only found in duodenum samples (6/21) (Fig. ​(Fig.2D).2D). A breakdown of the difference between duodenal and saliva abundance of all taxa is provided in Table S3. These differences in the relative abundance of specific taxa of microbes between paired saliva and duodenum samples also provide evidence against oral contamination in the duodenal samples.

Taxa co-correlations reveal disruptor taxa

We assumed that the taxa with the highest absolute abundance would have the highest potential for impacting the host. Thus, we began by analyzing the relationships between the top 20 most abundant genera. A co-correlation heatmap of these taxa revealed several distinct motifs (Fig. ​(Fig.3A):3A): (1) taxa whose absolute loads had a high correlation with total load, (2) taxa whose absolute loads had a higher co-correlation with another taxon’s absolute load than with total microbial load, (3) taxa with a mutually exclusive relationship with almost all other abundant taxa. Examples of the first motif are in the first column/row of the co-correlation heatmap in Fig. ​Fig.3A.3A. Correlation with total load was often an indicator of a prevalent taxon because the variance in total microbial load was larger than the variance in relative abundance. When two taxa have a higher co-correlation with each other than with total load (motif 2), it potentially indicates these taxa share preferred environmental factors or directly cooperate. One group of these co-correlating taxa that included several Prevotella species and a species of Porphyromonas matches a known shared metabolic niche in the oral cavity [35, 36] (Table S4).

Open in a separate window

Fig. 3

Co-correlations reveal which taxa co-occur in high abundance and which can be considered disruptor taxa. A Co-correlation matrix of the top 20 most abundant genera and total microbial load. Only significant correlations (q < 0.1, Benjamini–Hochberg correction) are shown. Color of each marker is determined by the sign of the Spearman’s correlation coefficient and size of each marker is determined by the magnitude of the coefficient. Disruptor taxa labels are bolded. B Clustered co-correlation matrix of the top 16 genera ranked by the difference between their maximum abundance and mean abundance. Two common genera in the dataset are shown at the bottom for reference. The color of each square indicates the Spearman correlation coefficient from negative (blue) to positive (red). Disruptor taxa labels are bolded. Taxa with known relevance to human health are indicated. Enterobacteriaceae and Escherichia-Shigella are unique sequence variants from the Enterobacteriaceae family but only Escherichia-Shigella could be classified at the genus level. HAI=hospital acquired infection; IBS, irritable bowel syndrome; IBD, inflammatory bowel disease; HACEK, Haemophilus, Aggregatibacter, Cardiobacterium, Eikenella, Kingella

Several genera stood out as having no significant correlation with almost all other abundant taxa (motif 3): Enterobacteriaceae, Escherichia-Shigella, Clostridium sensu stricto 1, and Lactobacillus (Fig. ​(Fig.3A).3A). For clarification, throughout the manuscript our references to Enterobacteriaceae and Escherichia-Shigella refer to unique sequence variants from the Enterobacteriaceae family, but only Escherichia-Shigella could be classified at the genus level. Based on evidence from a previous study [18] using the REIMAGINE cohort that found Klebsiella in several samples, we decided to measure the abundance of Klebsiella via qPCR in all samples containing a high abundance (at least 10⁵ 16S rRNA gene copies/mL) of Enterobacteriaceae. We found that the majority (16/22) of the samples with a high abundance of Enterobacteriaceae contained Klebsiella (Fig. S4A). Furthermore, in the samples containing Klebsiella, there was a high correlation (Pearson, 0.88, P < 0.001) between Klebsiella load and Enterobacteriaceae load (Fig. S4B). These taxa appeared to disrupt the commonly observed microbial structure (i.e., the prevalent taxa that generally co-correlate with one another) of the duodenal microbiome. This pattern of mutual exclusivity can be represented algorithmically by sorting all taxa by the difference between their maximum abundance and their mean abundance. Practically, this means that these disruptors are relatively rare (i.e., present in a small fraction of samples), but when they are present they usually dominate, excluding other common taxa. A clustered heatmap of the top 16 taxa as ranked by the difference in their maximum and mean abundances reveals two taxonomic signatures (Fig. ​(Fig.3B).3B). The first signature in the top left of the heatmap contained the mutually-exclusive taxa from the co-correlation heatmap, along with Enterococcus, Romboutsia, Aeromonas, and Bacteroides. The second signature contained taxa that were generally found in lower abundance, many of which are from the HACEK (Haemophilus, Aggregatibacter, Cardiobacterium, Eikenella, Kingella) group of organisms associated with infective endocarditis [34]. However, the second group also clustered with more common taxa in this dataset, such as Prevotella and Fusobacterium. Thus, we initially labelled all eight of the taxa in the first taxonomic signature as “disruptors” (Fig. ​(Fig.3B,3B, bolded taxa) because their presence appeared to be mutually exclusive with many other common taxa.

Aerobic disruptor taxa displace strict anaerobes and decrease diversity

After performing the co-correlation analysis, we ran a principal component analysis (PCA) on the absolute taxon abundances to investigate the drivers of variance in the dataset (Fig. ​(Fig.4A).4A). Total loads spanned 5 orders of magnitude, accounting for most of the variance. Total load cleanly separated samples along the PC1 axis. The second most explanatory axis, PC2, strongly correlated with the Shannon diversity index of samples (Spearman, 0.74, P < 0.001, N = 250). Ranked feature loadings for PC2 (Fig. ​(Fig.4B)4B) indicated that many of the disruptor taxa (dark blue) are the main drivers of separation in the positive direction of PC2 whereas the five taxa driving most of the separation in the negative direction (light blue) of PC2 consisted of four strict anaerobes (Porphyromonas, Leptotrichia, Prevotella, Prevotella 7) and one obligate aerobe (Neisseria). It should be noted that many more taxa were strongly associated with the negative direction of PC2 than the positive direction. This separation matches well with the mutual exclusivity seen between the disruptor taxa and other organisms in the co-correlation analysis. The two disruptor taxa with the highest loads are aerobic pathogens from the Enterobacteriaceae family and the taxa most associated with the negative direction of PC2 were strict anaerobes, so we next took a closer look at the composition of strict vs facultative anaerobes in each sample. We found a nearly 1:1 correlation between the strict and facultative anaerobe loads across all samples (Fig. ​(Fig.4C).4C). Additionally, the fraction of strict anaerobes in a sample was strongly correlated (Pearson, 0.71, P < 0.001, N = 250) with Shannon diversity (Fig. ​(Fig.4D),4D), indicating that the disruptor taxa appear to be mutually exclusive with strict anaerobes and the “bloom” of absolute abundance of disruptors decreases Shannon diversity. Furthermore, in half of the samples containing the two most common disruptor taxa (Enterobacteriaceae and Escherichia-Shigella), the total microbial loads were greater than 10⁷ 16S rRNA gene copies/mL, indicating a clear enrichment of disruptor taxa in samples with higher than average total microbial loads (Fig. ​(Fig.4E).4E). This signature of higher than average total microbial loads and mutual exclusivity with other microbes has been observed in some pathogenic microbial species [37, 38].

Open in a separate window

Fig. 4

Strict anaerobes and disruptor taxa control diversity. A PCA plot of absolute microbial abundances at the genus level with the top two correlated metadata variables overlaid. B Feature loadings for principal component 2. Top five value-ranked genera in each direction (positive and negative) are highlighted and labeled. C Correlation between the strict anaerobic microbial load and facultative anaerobic microbial load. D Relationship between the percentage abundance of strict anaerobes and Shannon diversity index. E Empiric cumulative distribution function (ECDF) plot for Enterobacteriaceae (N = 33), Escherichia-Shigella (N = 24), Campylobacter (N = 59), Lactobacillus (N = 42), and the common taxa Prevotella (N = 104)

Human vs mouse small-intestine microbiome

Several findings from this study emphasize how different the small-intestine microbiome is between mice and humans. Our previous study revealed that the coprophagic nature of mice resulted in total microbial loads spanning approximately one order of magnitude from 5 × 10⁸–5 × 10⁹ 16S rRNA gene copies/mL [8] in the small intestine while our human cohort spanned 5 orders of magnitude with a median of 10⁶ copies/mL (Fig. S11). Additionally, neither the most common disruptor family, Enterobacteriaceae, nor any of the taxa with at least 50% prevalence in this study were commonly found in our previous study examining microbial loads in the mouse small intestine [8]. Instead, in that study we found that the mouse small intestine was dominated by Lactobacillus and, as a result of coprophagy, several stool microbes [8]. The total microbial load of stool is similar between mice and humans [39] and they both share several common taxonomic groups [40]. These differences should be considered when using mice to model human health or disease impacted by the small intestine.

Value of quantitative analysis

The nearly 5 orders of magnitude spread in total microbial loads in the duodenum of these patients revealed the value of utilizing an absolute abundance measurement technique when analyzing microbial communities. Analyzing absolute abundances of individual taxa also let us filter out likely contaminants using Poisson loading statistics, which is critical for samples with low microbial abundance, such as those sometimes found in the human small intestine [41, 42]. The range of total loads in saliva and in stool each appear to be smaller than in the duodenum, closer to two orders of magnitude, which likely relates to differences in residence times, nutrient availability, and host defenses at these two sites compared with the small intestine [39]. Another benefit of using absolute rather than relative abundance measurements is the improved accuracy of correlations between microbes and host phenotype. For example, the 10 patient samples with the greatest disruptor loads had the highest prevalence of severe GI symptoms, but these samples had relative abundances of disruptor taxa that ranged from 8 to 97%. This wide range of relative abundances made samples with high disruptor loads indistinguishable from samples with intermediate disruptor loads when analyzing relative abundances.

Microbial connection between oral cavity and small intestine

The majority (89%) of identified microbial taxa in the paired duodenum samples were also present in the paired saliva samples. Our data supports the hypothesis of oral-duodenal transmission of microbes but a larger paired study utilizing shotgun metagenomic sequencing techniques would provide stronger evidence for this claim. Survival of microbes after ingestion is likely dependent on many host factors, including stomach-acid levels, bile secretions, antimicrobial-peptide production, and GI motility. The bimodal taxon abundance distributions (Fig. ​(Fig.1C)1C) observed for some taxa, including Prevotella, may indicate two subsets of patients with distinct stomach and/or duodenal environments that allow for differential abundance of specific taxa. For example, Campylobacter concisus, one of the most common oral Campylobacter species, is known to be sensitive to both stomach and bile acids [34]. Therefore, one could hypothesize that if a patient had low levels of stomach or bile acids some C. concisus may survive ingestion. Low-acid conditions could also allow many other bacteria to survive transit to the duodenum, resulting in higher total microbial loads in the small intestine. We suspect we observed something similar in our samples; the Campylobacter genus was only found in samples with greater than average total microbial loads (Fig. ​(Fig.4E).4E). However, we did not observe a relationship between total microbial loads in the duodenum and the patients’ use of proton pump inhibitors (PPI), which are known to reduce acid production. PPI impact on survival of microbes between the oral cavity and duodenum may be dependent on how recently the PPI was taken, however this information was not collected from patients in the REIMAGINE study. A conclusive comparison of the relative importance of various factors affecting bacterial survival in the duodenum would require additional information on small-intestine secretions of bile acids and antimicrobial peptides in these patients.

Several common oral microbes have been implicated in GI diseases when present in stool [30, 43]. A high microbial load in the small intestine could increase the likelihood of these microbes surviving all the way down the GI tract. The shared taxa between the small intestine and oral microbiota in our paired saliva-duodenum samples provides evidence that blooms of opportunistic pathogens in the mouth could also lead to colonization in the SI [30]. In this study, only 1 of the 21 paired duodenum-saliva samples contained disruptor taxa in the duodenum, but these taxa were not present in the corresponding saliva sample. Several Enterobacteriaceae species have been identified in oral samples [44] but usually at a low frequency in healthy populations. Many Enterobacteriaceae species are introduced into the gut from contaminated food and water sources [45] which would likely result in only transient oral residence. However, persistent oral Enterobacteriaceae species have been linked to the use of dentures and the presence of periodontal disease [46]. All the taxa we classified as disruptors in this study are more frequently found in stool than in the small intestine or oral cavity [30, 31]. Further studies should be performed to determine the source of disruptor taxa in the upper GI tract.

A number of taxonomic groups we identified in the duodenum have members known to be opportunistic pathogens. Beyond disruptor taxa, several taxa from the HACEK group of organisms [47] associated with infective endocarditis were found in high abundance in the duodenum. The route that these and other opportunistic pathogens take to reach the blood stream is not clear but our data show that the HACEK organisms are not limited to the oral cavity. The same traits that allow them to colonize the mouth and heart (biofilm production [48], and general resistance to most host secretions) likely contribute to their ability to survive in the small intestine. Additionally, in mouse models, the transmission of opportunistic pathogens, like Klebsiella, from the oral cavity to the intestine has been shown to induce inflammation [30]. The oral cavity presents a potential reservoir for a wide range of opportunistic pathogens that have been linked to GI disorders.

Potential relationship between oxygen and disruptor taxa

Several colonic GI disorders are linked to increased oxygen levels in the lumen resulting from decreased epithelial integrity and inflammation [49]. However, the barrier properties of the small intestine, an absorptive organ, are different from those of the colon. To our knowledge, shifts in absolute abundance of microbes capable of aerobic respiration and anaerobes have not been quantitatively studied previously in the human small intestine. The highly correlated abundance of both strict and facultative anaerobes that we observed could be a function of the oxygen gradients in the gut from the epithelial surface to the center of the lumen [50]. In our study, when diversity collapsed and disruptor taxa bloomed, the microbial composition shifted away from strict anaerobes to taxa capable of aerobic respiration. One clear outlier was a Clostridium classified as C .perfringens, which is a strict anaerobe but was highly correlated with the Enterobacteriaceae genera classified as disruptors. Previous mutualistic relationships between aerobic and anaerobic species that could help facilitate colonization have been observed in other studies with Bacteroides fragilis and either Klebsiella pneumonia or Escherichia coli [51, 52]. We have previously hypothesized that the surprising coexistence of aerobe-anaerobe communities can occur in multi-stable systems, and that these communities can persist due to hysteresis [51]. Although multi-stability and hysteresis have not yet been documented in the gut microbiome, this phenomenon could explain the unexpected coexistence and persistence of aerobe-anaerobe communities in the small intestine.

Disruptor taxa predict SIBO classification and likelihood of GI symptoms

Clinically, SIBO is classified by culture of duodenal aspirates on aerobic MacConkey agar or measurement of exhaled hydrogen and methane after intake of a fermentable sugar solution [10, 53]. The main disruptor taxa (Enterobacteriaceae) grow well on MacConkey agar plates, which may explain the high correlation between SIBO classification and samples with disruptor taxa. It is commonly hypothesized that overgrowth of these taxa in the small intestine is responsible for the gas production detected during a breath test, and our study further supports this understanding because we found a correlation between bloating symptoms (attributable to gas production) and disruptor taxa. Future studies should determine whether individuals with and without high loads of disruptor taxa yield positive breath test results. Our findings support a strong relationship between overgrowth of specific disruptor taxa and GI symptoms in subjects with SIBO. High total microbial load alone in the small intestine was not associated with GI symptoms usually observed in subjects with SIBO and other GI conditions and diseases. Microbial culture is never perfect and will not capture all taxa associated with SIBO and GI conditions. However, our data suggest that SIBO diagnosis via microbial culture should focus on quantification of a specific group of disruptor taxa (Enterobacteriaceae) rather than total microbial load. Additionally, SIBO diagnosis via quantitative sequencing should focus on the absolute abundance of the seven disruptor taxa identified in this study.

Lactobacillus seemed to be an exception among the disruptor taxa in several ways. It commonly co-occurred with other disruptors; however, it was also present in many “normal” samples at low abundance. Additionally, when present at high total loads in the presence of other disruptor taxa, Lactobacillus load had a negative correlation with bloating score. However, Lactobacillus also dominated several samples that had no other disruptor taxa but had high symptom scores. It should also be noted that individuals taking probiotics (N = 49) did not have increased prevalence or abundance of Lactobacillus in the duodenum. Overall, finer taxonomic resolution may be required to decipher the role of different Lactobacillus species and strains. Their impact on human health is likely also dependent on the overall microbial community and host environment.

Although most patients in this study have various GI complications that could result in abdominal symptoms independent of a microbial component, patient samples with high loads of disruptor taxa had a substantially higher likelihood of having many severe GI symptoms. However, total microbial load alone did not associate with GI symptoms. Of the 13 cytokines and chemokines measured, only IL8 levels were significantly higher in the serum of patients with disruptor taxa, potentially indicating an associated local inflammatory process. Future studies that analyze biopsy transcriptomes would be needed to determine whether there is an associated host response, such as immune infiltration or epithelial stress responses in regions with disruptor taxa and/or high total microbial loads.

We initiated this study with four expectations, only one of which was supported by our data. Because mice are coprophagic and humans are not, we expected to see a dramatic difference between mouse and human small-intestine microbiomes. We indeed observed large quantitative and qualitative differences between the two. However, we were more surprised and educated by the three expectations that were shown to be incorrect. First, we expected microbial load in the human duodenum to have a bimodal distribution, with low microbial loads for non-SIBO patients and much higher load for SIBO patients, which our findings did not support (Fig. ​(Fig.5B).5B). Second, because stomach acid and bile acid secretions isolate the duodenum from the upper GI tract and because the unidirectional flow of digesta and the ileocecal valve isolate the small intestine from the colon, we expected to find a unique population of microbes in the duodenum. We were surprised by the extent to which the oral microbiota appeared to influence the small-intestine microbiota (Figs. ​(Figs.1C1C and ​and2).2). Third, we expected to see microbiomes dominated by taxa generally thought of as commensals like Lactobacillus and Bifidobacterium. We were surprised by the prevalence and abundance of taxa known to be human pathogens (Figs. ​(Figs.1C1C and ​and3B),3B), especially given that the small intestine is an immune-rich, absorptive organ with a loose mucus structure that likely permits substantial exposure to microbial cells and microbial-associated proinflammatory molecules.

Questionnaires

Prior to EGD, all subjects completed a study questionnaire documenting demographic information and family and medical history, including GI disease and bowel symptoms, medication use, use of alcohol and recreational drugs, travel history, and dietary habits and changes. Subjects also reported any known underlying conditions, such as GI diseases and disorders, neurologic disease, hematologic disease, autoimmune disease, kidney disease, heart disease, and cancer. All medical information provided by subjects was verified through audits of medical records. All data were de-identified prior to analysis.

Blood collection and analysis

After completing the study questionnaire, fasting blood samples were collected in BD Vacutainer SST tubes (Becton Dickson, Franklin Lakes, NJ, USA). Levels of circulating pro- and anti-inflammatory cytokines and chemokines were analyzed on a Luminex FlexMap 3D (Luminex Corp., Austin, TX, USA) using a bead-based multiplex panel that included: GM-CSF, IFNγ, IL10, IL12P70, IL13, IL1B, IL2, IL4, IL5, IL6, IL8, MCP1, and TNFα (EMD Millipore Corp., Billerica, MA, USA, cat. #HCYTOMAG-60K).

Saliva and small-intestine lumenal sample collection

Prior to EGD procedure, saliva was collected in a sterile 5 mL tube. During the EGD procedure, samples of duodenal lumenal fluid were procured using a custom-designed sterile aspiration double-lumen catheter (Hobbs Medical, Inc.) [28]. Duodenal aspirates (DA) were collected using a custom-designed sterile inner catheter which was pushed through a sterile bone wax cap only after the endoscopist entered the second portion of the duodenum, in order to reduce contamination from the mouth, esophagus, and stomach. After collection, samples were immediately placed on ice and transferred to the laboratory for further analysis.

Aspirate processing and microbial culture

Prior to microbial culture, an equal volume of sterile 6.5 mM dithiothreitol (DTT) prepared with RNase and DNase PCR-grade sterile water was added at a 1:1 ratio to each saliva and duodenal aspirate (~ 1 mL) and the samples were vortexed until fully liquified (~ 30 s) as described previously [28]. A 100-μl aliquot of each duodenal sample (DA + DTT) was then serially diluted with 900 μL sterile 1× PBS and plated on MacConkey agar (Becton Dickinson), and on blood agar (Becton Dickinson). Plates were incubated at 37 °C for 16–18 h under aerobic (MacConkey) or anaerobic (blood agar) conditions. Plates without bacterial growth after 18 h were re-incubated for an additional 18 h. Colony forming units (CFU) were then counted electronically using a Scan 500 (Interscience, Paris, France). Saliva + DTT and the remainder of each DA+DTT were centrifuged at maximum speed (> 13,000 RPM) for 5 min. The supernatant was removed, and 1 mL of sterile Allprotect reagent (Qiagen, Hilden, Germany) was added to the microbial pellet and then stored at − 80 °C.

DNA isolation

On the day of the DNA isolation, DA pellets were thawed on ice and processed as described previously [28]. Microbial DNA was isolated using the MagAttract PowerSoil DNA KF Kit (Qiagen) on a KingFisher Duo (Thermo Fisher Scientific, Waltham, MA, USA), and quantified using Qubit dsDNA high sensitivity and Qubit dsDNA BR Assay kits (Invitrogen by Thermo Fisher Scientific) on a Qubit 4 Fluorometer (Invitrogen, Carlsbad, CA, USA).

16S rRNA gene sequencing

Extracted DNA was amplified, barcoded, and sequenced as described previously [8, 9, 29]. Briefly, amplification of the variable 4 (V4) region of the 16S rRNA gene was performed in 20 μL duplicate reactions with: 8 μL of 2.5× 5Prime Hotstart Mastermix (VWR, Radnor, PA, USA), 1 μL of 20× Evagreen (VWR), 2 μL each of 5 μM forward and reverse primers (519F, barcoded 806R, IDT, CoralVille, IA, USA), 3.5 μL of water, and 3.5 μL of extracted DNA template. A CFX96 RT-PCR machine (Bio-Rad Laboratories, Hercules, CA, USA) was used to monitor amplification reactions and all samples were removed in late exponential phase (~ 10,000 FRU) to minimize chimera formation and non-specific amplification [9, 54, 55]. Amplification was performed under the following cycling conditions: 94 °C for 3 min, up to 50 cycles of 94 °C for 45 s, 54 °C for 60 s, and 72 °C for 90 s. Several samples were rerun after diluting the template as they showed non-exponential amplification in the undiluted sample, a sign of PCR inhibition. Amplified duplicates were pooled together and quantified with KAPA library quantification kit (Roche, Basel, Switzerland) and then all samples were pooled at equimolar concentrations with up to 96 samples per library. AMPureXP beads (Beckman Coulter, Brea, CA, USA) were used to clean up and concentrate libraries before final library quantification with a High Sensitivity D1000 Tapestation Chip (Agilent, Santa Clara, CA, USA). Illumina MiSeq sequencing was performed with a 2 × 300 bp reagent kit by Fulgent Genetics (Temple City, CA, USA).

Raw reads were demultiplexed by Fulgent Genetics. Demultiplexed forward and reverse reads were processed with QIIME 2 2020.2 [56]. Loading of sequence data was performed with the demux plugin followed by quality filtering and denoising with the dada2 plugin [57]. Dada2 trimming parameters were set to the base pair where the average quality score dropped below thirty. All samples were rarefied to the lowest read depth present in all samples (45,386 reads) to decrease biases from varying sequencing depth between samples [58]. The q2-feature-classifier [59] was then used to assign taxonomy to amplicon sequence variants (ASV) with the Silva [60] 132 99% OTUs references. Resulting read count tables were used for downstream analyses in IPython notebooks (see “Data availability” section).

Klebsiella-specific qPCR

Primers specific for the Klebsiella gltA gene [61] (F: 5′-CAGGCCGAATATGACGAATTC-3′, R: 5′-CGGGTGATCTGCTCATGAA-3′) were first informatically evaluated for coverage across Klebsiella pneumoniae, Klebsiella oxytoca, and Klebsiella aerogenes via Primer-BLAST [62]. This primer set was found to have a perfect match against strains from all three tested Klebsiella species. These primers were also evaluated in the lab for specificity against Escherichia coli. No amplification after 40 cycles was observed with a DNA equivalent of ~ 10⁶ E. coli cells from the Zymo microbial community DNA standard (Zymo Research, Irvine, CA, USA). Klebsiella qPCR was performed in 10 μL reactions with 5 μL of Ssofast Evagreen Supermix (Bio-Rad Laboraties), 0.5 μL of 10 μM gltA primers, and 3.5 μL of water. A CFX96 RT-PCR machine (Bio-Rad Laboratories) was used for amplification with the following cycling conditions: 95 °C for 3 min, 40 cycles of 95 °C for 15 s, 62 °C for 30 s, and 68 °C for 30 s. Estimated conversion of cycle threshold (Cq) to copies/μL was performed where a Cq of 22.4 equals 1000 copies/uL. Klebsiella load was then calculated by adjusting for dilutions and normalizing to the collected sample volume.

Absolute abundance

The total microbial load (bacteria and archaea) of each sample and the absolute abundance of each taxon in individual samples was determined as described previously [9, 29]. Briefly, the Bio-Rad QX200 droplet dPCR system (Bio-Rad Laboratories) was utilized to measure the 16S concentration in each sample with the following reaction components: 1X QX200 EvaGreen Supermix (Bio-Rad), 500 nM forward primer, and 500 nM reverse primer (519F, 806R) and thermocycling conditions: 95 °C for 5 min, 40 cycles of 95 °C for 30 s, 52 °C for 30 s, and 68 °C for 30 s, followed by a dye stabilization step of 4 °C for 5 min and 90 °C for 5 min. The final concentration of 16S rRNA gene copies in each sample was corrected for dilutions and normalized to the extracted sample volume.

For each sample, the input-volume-normalized total microbial load from dPCR was multiplied by each amplicon sequence variant’s (ASV) relative abundance to determine the absolute abundance of each ASV. No correlation between collected sample volume and measured bacterial load was observed. The average of all sample volumes for a specific sample type was used for a few samples (11 duodenum, 10 saliva) that were missing the starting volume information. A 95% confidence interval of input volumes for duodenum samples ranged from 0.18 to 1.93 mL indicating that the estimated input volume measurement would likely be up to 4× off in either direction while the total microbial load ranged 40,000X. Similarly, a 95% confidence interval of input volumes for saliva samples ranged from 0.36 to 1.28 mL indicating that the estimated input volume measurement would likely be up to 2× off in either direction while the total microbial load ranged 82X.

Poisson quality filtering

Two separate quality-filtering steps based on Poisson statistics were used to determine the statistical confidence in the measured values. First, a 95% confidence interval was calculated from the repeated measures of water blanks. Samples with a total microbial load below the upper bound of this confidence interval were removed from further analysis.

Second, the limit of detection (LOD) in terms of relative abundance was determined for each sample. Sequencing can be divided into two separate Poisson sampling steps. First, an aliquot of sample is taken from the extracted sample and input into the library amplification reaction. The LOD of the library amplification step was determined by multiplying the total microbial load from dPCR by the input volume into the library amplification reaction and then finding the relative abundance corresponding to an input of three copies. Poisson statistics tells us that the likelihood of sampling one or more copies with an average input of three copies is 95%. The second Poisson sampling step in sequencing arises from the number of reads generated from the amplified library. The accuracy of the second Poisson sampling step was previously shown [9] to follow a negative exponential curve, LOD=7.115read depth^−0.115, between the total read depth and relative abundance at which 95% confidence of detection is observed. The minimum of the two described LODs (first determined per sample by total load, and second by sequencing depth) was then determined for each sample. For each sample, the abundance of any ASV with a relative abundance below the LOD was set to zero. After filtering, data tables for each taxonomic level were generated.

Data transforms and dimensionality reduction

For PCA, all absolute taxon abundances were log-transformed. To handle zeros, a pseudo-count of 0.1 reads was added to all taxon relative abundances before multiplying by each sample’s total microbial load as determined by digital PCR. PCA was performed with the sklearn.decomposition.PCA function in Python. Ranked feature loadings for each taxon on a given principal component were determined by scaling the corresponding eigenvector by the maximum transformed value for that principal component axis.

We thank the Caltech Bioinformatics Resource Center for assistance with statistical analyses, Jenny Ji for related analyses and Natasha Shelby for contributions to writing and editing this manuscript. We acknowledge OpenMoji for use of the saliva and stool graphics in Fig. ​Fig.1.1. We thank Stacy Weitsman, Walter Morales and Maria Jesus Villanueva-Milan from MAST for assistance with sample processing and data curation from the REIMAGINE study. We also thank the Gastroenterology team at Cedars-Sinai Medical Center for assistance with patient recruitment and endoscopy procedures.