Integrating gut microbiome and host immune markers to understand the pathogenesis of Clostridioides difficile infection.

Clostridioides difficile (C. difficile) infection is the most common cause of healthcare-associated infection and an important cause of morbidity and mortality among hospitalized patients. A comprehensive understanding of C. difficile infection (CDI) pathogenesis is crucial for disease diagnosis, treatment, and prevention. Here, we characterized gut microbial compositions and a broad panel of innate and adaptive immunological markers in 243 well-characterized human subjects (including 187 subjects with both microbiota and immune marker data), who were divided into four phenotype groups: CDI, Asymptomatic Carriage, Non-CDI Diarrhea, and Control. We found that the interactions between gut microbiota and host immune markers are very sensitive to the status of C. difficile colonization and infection. We demonstrated that incorporating both gut microbiome and host immune marker data into classification models can better distinguish CDI from other groups than can either type of data alone. Our classification models display robust diagnostic performance to differentiate CDI from Asymptomatic carriage (AUC~0.916), Non-CDI Diarrhea (AUC~0.917), or Non-CDI that combines all other three groups (AUC~0.929). Finally, we performed symbolic classification using selected features to derive simple mathematic formulas that explicitly quantify the interactions between the gut microbiome and host immune markers. These findings support the potential roles of gut microbiota and host immune markers in the pathogenesis of CDI. Our study provides new insights for a microbiome-immune marker-derived signature to diagnose CDI and design therapeutic strategies for CDI.

Introduction

Clostridioides difficile infection (CDI) is the most common cause of healthcare–associated infection and an important cause of morbidity and mortality among hospitalized patients[1-3]. Exposure to toxinogenic C. difficile can lead to a range of clinical outcomes ranging from asymptomatic colonization to mild diarrhea and more severe disease syndromes such as pseudomembranous colitis, toxic megacolon, bowel perforation, sepsis, and death[4,5]. Asymptomatic C. difficile carriage is mainly characterized by C. difficile colonization in the absence of symptoms of infection. The diagnosis of CDI is based on clinical signs and symptoms in combination with laboratory testing, including enzyme immunoassays (EIA) for TcdA and TcdB, nucleic acid amplification tests (NAAT), selective toxinogenic culture, cell cytotoxicity neutralization assay, and glutamate dehydrogenase EIA[6-8]. However, currently available approaches do not accurately differentiate CDI from diarrhea with another cause in a patient colonized with toxinogenic C. difficile.

Current treatment strategies for CDI, including vancomycin and fidaxomicin, have inconsistent cure rates and treatment failure or CDI recurrence may occur in approximately one-third of cases[9,10]. Antibiotic exposure is considered the most important factor predisposing patients to CDI[11,12]. In fact, treatments with antibiotics have a tremendous impact on the composition and functionality of the gut microbiota, and accordingly are associated with reduced colonization resistance against pathogens such as C. difficile[13-15]. It has been reported that several gut commensal bacteria may contribute to the prevention of C. difficile colonization and infection[16,17]. Once colonized, C. difficile can produce toxins that mediate a robust inflammatory response. Toxin A (TcdA) and toxin B (TcdB) are the primary virulence factors of C. difficile[18] and act on intestinal epithelial cells first, inducing pro-inflammatory cytokines, loss of tight junctions, cell detachment and an impaired mucosal barrier[19-21] leading to further exposure of immune cells to toxins. The innate and adaptive immune responses to CDI play crucial roles in disease onset, expression, severity, progression, and overall prognosis[22,23]. The innate immune defense mechanisms against C. difficile and its toxins include the commensal intestinal organisms, mucosal barrier, intestinal epithelial cells, and mucosal immune system[24,25]. TcdA and TcdB have multiple effects on the innate immune system, including inducing expression of numerous pro-inflammatory mediators (e.g., cytokines, chemokines and neuroimmune peptides) and the recruitment and activation of a variety of innate immune cells[26,27]. Adaptive immunity is also sufficient to provide some protection from CDI, likely via antibody-mediated neutralization of TcdA and TcdB[28-31]. The role of the immune response combined with the knowledge that a balanced microbiota can prevent colonization and infection demonstrates the importance of combining both gut microbiota and host immune markers in understanding the pathogenesis of CDI.

Machine learning has a great impact in many areas of medical research, as it offers a principled approach for developing sophisticated, automatic, and objective algorithms for analysis of complex data. Indeed, previous studies indicate that supervised learning can be successfully employed for clinical disease assessment for diverse disorders[32-35]. In our previous work, we found that specific immune markers, particularly G-CSF, can be used to distinguish adults with CDI from other groups including asymptomatic carriers and NAAT-negative patients with and without diarrhea[36]. Here, we leverage machine learning tools to integrate the host immune marker data and newly obtained gut microbiome data from subjects of the same cohort to identify collections of bacteria and immune markers that can be associated with CDI. Our aim is to quantify the role of intricate interactions between gut microbiota and immune response in CDI pathogenesis, which can inform the design of microbiome-immune marker-based diagnostic test and therapeutic strategies.

Results

Baseline demographic and clinical characteristics of participants

Our clinical cohort consists of 243 well-characterized recruited participants, who were divided into four groups (see Methods)[36]: (1) Control (n = 47); (2) Non-CDI Diarrhea (n = 44); (3) Asymptomatic Carriage (n = 40); (4) CDI (n = 112). The first three groups can be combined as the Non-CDI group. The entire clinical cohort had a mean ± SD age of 63.66 ± 14.85 year and was 48.15% female. Demographic data of the cohort are summarized in Table 1. In total, 187 participants (76.95%) had both gut microbiome and immune marker data available (see Supplementary Table 1).

Table 1.

Demographic characteristics of the enrolled subjects

	NAAT negative		NAAT positive
Characteristics	Control (n = 47)	Non-CDI Diarrhea (n = 44)	Asymptomatic Carriage (n = 40)	CDI (n = 112)
Sex
Female	14 (29.79%)	22 (50.00%)	20 (50.00%)	61 (54.46%)
Male	33 (70.21%)	22 (50.00%)	20 (50.00%)	51 (45.54%)
Age, Avg ± SD	62.40 ± 12.33	63.07 ± 13.15	62.15 ± 17.25	64.99 ± 15.62
Ethnicity
Hispanic	1 (2.13%)	3 (6.82%)	1 (2.50%)	6 (5.36%)
Non-Hispanic	38 (80.85%)	37 (84.09%)	31 (77.50%)	96 (85.71%)
Unknown	8 (17.02%)	4 (9.09%)	8 (20.00%)	10 (8.93%)
Race
White	33 (70.21%)	28 (63.64%)	28 (70.00%)	89 (79.46%)
Other	4 (8.51%)	10 (22.73%)	3 (7.50%)	23 (20.54%)
Unknown	10 (21.28%)	6 (13.64%)	9 (22.50%)	0 (0.00%)

Microbial community structure

To compare the overall microbial community structure of the four groups, we first calculated the alpha diversity (i.e., the within-sample taxonomic diversity) of each sample at the genus level using four different measures: taxa richness (the observed number of different taxa present in the sample), Chao1 (abundance-based estimator of taxa richness), Evenness (the uniformity of the population size of each taxa present in the sample), and Shannon diversity index (estimator of taxa richness and evenness: more weight on richness). As shown in Figure 1(a-d), we found that taxa richness and Chao1 did not differ significantly among these groups. The gut microbiota of Non-CDI Diarrhea subjects showed lower evenness than that of the Control group. Shannon diversity was significantly lower in the Non-CDI Diarrhea and CDI groups than in the Control group.

Figure 1.

Comparing the diversity of the gut microbiota (and host immune markers) of subjects with different . (a) Taxa richness. (b) Chao1. (c) Evenness. (d) Shannon index. (e) Principal Coordinates Analysis (PCoA) plot based on Bray–Curtis dissimilarities of microbial compositions. (f) Boxplot of the gut microbiome Bray–Curtis dissimilarity between subjects within each group. (g) Principal component analysis (PCA) plot of host immune marker concentrations. (h) Boxplot of the Euclidean distance for the host immune markers of subjects within each group. Statistical significance was determined by Mann–Whitney test, *P < .05, *P < .05, **P < .01, ***P < .001

To determine whether the gut microbial compositions of participants are affected by C. difficile infection/colonization status, we performed Principal Coordinates Analysis (PCoA) at the genus level using Bray–Curtis dissimilarity (which is a beta diversity measure to quantify the between-sample compositional dissimilarity). We found no distinct clusters corresponding to the four different phenotype groups (Figure 1e). Interestingly, by directly comparing the beta diversity of each group, we did find that the CDI group displays higher beta diversity than other groups (Figure 1f), indicating that the microbial compositions of participants within the CDI group vary more prominently than other groups. PERMANOVA (permutational multivariate analysis of variance) showed that the overall bacterial composition differed significantly among different groups based on the CDI status (P < .001; Supplementary Table 2), whereas other host factors such as age, sex, race and ethnicity had no significant effect on the microbiome composition.

To identify microbiome markers (i.e., certain taxa with very high discriminatory ability) to differentiate those different phenotype groups, we performed differential abundance analysis. In particular, we used ANCOM[37] (analysis of composition of microbiomes) with a Benjamini-Hochberg correction, and adjusted for age and sex. We found that the abundances of 15 genera were significantly different between CDI and Asymptomatic Carriage groups (Figure 2a and Supplementary Table 3). Among the 15 genera, 4 of them (Veillonella, Enterobacter, Granulicatella, and Dialister) of these genera were enriched in the CDI group, while the other 11 genera (including Lactococcus, Dorea, Moryella, Stenotrophomonas, and Agathobacter) were enriched in the Asymptomatic Carriage group. We also found 16 differentially abundant genera between the Non-CDI Diarrhea group and the CDI group (Figure 2b and Supplementary Table 4). Of these, 10 genera (including Clostridioides, Enterobacter, Dialister, and Veillonella) were enriched in the CDI group, and the other 6 genera ([Eubacterium]_hallii_group, Collinsella, Agathobacter, Dorea, Stenotrophomonas, and Streptococcus) were enriched in the Non-CDI Diarrhea group. ANCOM analysis also enabled us to identify 40 genera (including Clostridioides and Veillonella) that have significant differential abundances between the CDI group and the whole Non-CDI group (Figure 2c and Supplementary Table 5). Note that a total of 6 differentially abundant genera were identified from all the three comparisons. Among them, Veillonella, Enterobacter and Dialister were enriched in the CDI group, while Dorea, Stenotrophomonas and Agathobacter were depleted in the CDI group.

Figure 2.

Relative abundances of differentially abundant genera identified by ANCOM in comparing different groups. (a) CDI vs. Asymptomatic Carriage. (b) CDI vs. Non-CDI Diarrhea. (c) CDI vs. Non-CDI. The top differentially abundant taxa were ranked based on their W statistics (a high “w score” generated by this test indicates the greater likelihood that the null hypothesis can be rejected, indicating the number of times a parameter is significantly different between groups) (from left to right). The relative abundance (%) are plotted on log10 scale. The notches in the boxplots show the 95% confidence interval around the median

Microbial correlation networks

As the gut microbiota is complex in both structure and function, this complexity can be well represented and modeled as networks. Network methods can be applied to microbiome studies to model the co-occurrence of microorganisms, find microbial relationships essential for community assembly or stability, and deduce the influence of various associations on the host health. To compare the microbial communities of the four groups at the network-level, we constructed the genus-level microbial correlation network for each group using SparCC[38] (sparse correlations for compositional data). We found that the microbial correlation network of the CDI group has quite different structure compared to other groups (Figure 3). More precisely, it has fewer nodes and edges, lower average degree, but higher modularity (Supplementary Table 6). These indicate that the overall microbial correlations in the CDI group are much weaker than those in other groups.

Figure 3.

Microbial correlation networks of different groups. (a) Control. (b) Non-CDI Diarrhea. (c) Asymptomatic Carriage. (d) CDI. Nodes represent genera and are colored based on their phylum. Edges represent microbial correlations: green/red means positive/negative correlations, respectively. Edge thickness indicates correlation strength, and only the high-confidence interactions (p-value < 0.05) with high absolute correlation coefficients (> 0.3) were presented. For each group, we further identified the top-three most connected genera/nodes. They are Ruminococcus_1, Roseburia and Lachnospiraceae_UCG_008 for the Control group, [Ruminococcus]_torques_group, [Eubacterium]_hallii_group and Blautia for the Non-CDI Diarrhea group, Ruminiclostridium_5, Enterococcus and Lachnospiraceae_UCG_008 for the Asymptomatic Carriage group, and Alistipes, Ruminiclostridium_5 and Lachnoclostridium for the CDI group

To analyze these patterns in more detail, we used NetShift[39] to identify potentially important “driver” taxa responsible for the change of microbial correlations. In the NetShift pipeline, the common taxa present in both ‘control’ and ‘case’ sample sets with a minimum abundance threshold are extracted to construct the microbial correlation networks. Then the ‘driver taxa’ are identified based on their Neighbor Shift (NESH) scores and the Betweenness Centrality (BC) measures. In a nutshell, the NESH score of a taxon/node is minimum when the associated partners of this node are the same in both ‘case’ and ‘control’ networks, intermediate when there is only a subset of the associated partners present in the ‘case’ network, and maximum when a completely new set of associated partners appear in the ‘case’ network. The Betweenness Centrality of a taxon/node quantifies its involvement in connecting other nodes in the network. A taxon with an altered set of edges (identified by a high NESH score), while still being increasingly important (i.e., with higher scaled BC in the ‘case’ network than in the ‘control’ network), necessarily holds a key significance in microbial interplay and is identified as a ‘driver’ taxon. This analysis revealed 24 potential driver taxa linked with the change of microbial correlations between CDI and Asymptomatic Carriage groups (Supplementary Figure 1). The top driver taxa were Alistipes, Clostridioides, Desulfovibrio, Eggerthella, Erysipelatoclostridium, Klebsiella, Odoribacter, Proteus, [Ruminococcus]_torques_group, Streptococcus, Vagococcus and Veillonella. We then identified 24 genera as potential driver taxa underlying the change of microbial correlations between CDI and Non-CDI Diarrhea groups (Supplementary figure 2). The top driver taxa were Alistipes, Buttiauxella, Citrobacter, Clostridium_sensu_stricto_13, Desulfovibrio, Klebsiella, Oscillibacter, Phascolarctobacterium, Streptococcus and Veillonella. Finally, Netshift analysis revealed 38 potential driver taxa underlying the change of microbial correlations between CDI and Non-CDI groups. The top driver taxa were Bifidobacterium, Clostridioides, Klebsiella, Oscillibacter, Streptococcus and Veillonella (Supplementary Figure 3). Together, these results suggested that certain bacterial taxa (e.g., Clostridioides, Klebsiella, Streptococcus and Veillonella) could play an important role in driving the changes of microbial correlations in subjects with different C. difficile infection/colonization status.

Host immune markers and CDI

To determine the systemic levels of proinflammatory cytokines in CDI, we measured the circulating levels of granulocyte-colony stimulating factor (G-CSF), interleukin-1β (IL-1β), IL-2, IL-4, IL-6, IL-8, IL-10, IL-13, IL-15, monocyte chemoattractant protein-1 (MCP-1), vascular endothelial growth factor-A (VEGF-A), tumor necrosis factor-alpha (TNF-α), and serum concentrations of immunoglobulin A (IgA), IgG, and IgM antibodies against C. difficile toxin A and toxin B as previously reported[36]. We previously demonstrated specific markers of these innate and adaptive immunity that can distinguish CDI from each of the other three groups[36]. In the current study, we are particularly interested in comparing the CDI group and the combined Non-CDI group. Based on the Mann-Whitney U test, we identified in total 11 immune markers that displayed significantly different concentrations in these two groups, including G-CSF, IL-4, IL-6, IL-8, IL-10, IL-15, TNF-α, MCP1, IgA anti-toxin A and B, and IgG anti-toxin A in blood (Supplementary Table 7). All of these immune markers had higher concentrations in the CDI group than in the Non-CDI group. Host immune marker variations between samples were evaluated using the Principal Component Analysis (PCA) (Figure 1g). PCA plot showed no clear clustering of those subjects based on immune marker concentrations. However, boxplot of Euclidean distance of immune marker profiles from CDI patients showed higher within-group variation than that in all the other three groups (Figure 1h). PERMANOVA analysis indicated that the immune homeostasis was significantly different among different groups based on the CDI status (P = .016; Supplementary Table 2).

Interactions between gut microbiome and host immune markers

To reveal the interactions between the gut microbiome and the host immune system, we calculated the correlations between microbial compositions and the circulating levels of host immune markers for each of the four groups. The results are shown in Figure 4 and Supplementary Figure 4. For the Control group, the most significant correlations were identified as Chiristensenellaceae R-7 group negatively correlated with TNFα, Bifidobacterium positively correlated with VEGFA and IL-13, Rothia positively correlated with IL-15, and Veillonella positively correlated with IL-4 (Figure 4a and Supplementary Figure 4). For the Non-CDI Diarrhea group, Ruminococcaceae UCG-011 was negatively correlated with IL-8 and IL-6, Defluviitaleaceae UCG-011 was positively correlated with IL-1β, and Blautia was negatively correlated with MCP1 levels (Figure 4b). For the Asymptomatic Carriage group, we found that Lactobacillus was negatively correlated with VEGFA, Akkermansia was positively correlated with IL-6, and Enterococcus was positively correlated with TNFα (Figure 4c). For the CDI group, negative correlations involved Akkermansia and IL-10, Lactococcus and G-CSF, while positive correlations involved Lactobacillus and IgG and IgA anti-toxin B (Figure 4d). Interestingly, none of these most significant correlations was universally present across different groups. This indicated that the interactions between gut microbiota and host immunological markers can be very sensitive to the status of C. difficile colonization and infection. Although the rudimentary correlation analysis cannot reveal any nonlinear interactions between gut microbiota and host immune markers, the result implies that the integration of gut microbiota and host immune markers might be quite useful for highly accurate classification of CDI.

Figure 4.

Correlations between gut microbial abundances and host immune markers in different groups, quantified by Spearman correlation with Benjamini-Hochberg correction. (a) Control. (b) Non-CDI Diarrhea. (c) Asymptomatic Carriage. (d) CDI. Rows represent genera; columns represent immune markers. The layout of the heatmap is followed the hierarchical clustering results of Control cohort (see Supplementary Figure 4). Red/blue represents positive/negative correlation, respectively. The intensity of the colors denotes the strength of the correlation. *α < 0.05, **α < 0.01, ***α < 0.001

Classification of CDI using host immune markers and gut microbiota

To determine whether host immune markers or gut microbiota could serve as biomarkers to classify subjects into different groups, we constructed a multi-class classifier based on random forests (RF). One of the most popular performance metrics of a classifier is the Area Under the receiver operating characteristic Curve (AUC). The performance of a multi-class classifier is measured by both micro-average and macro-average AUCs. We considered three different feature types: (1) host immune maker concentrations alone; (2) gut microbial compositions alone; and (3) the integration of (1) and (2) in our classification analysis. To eliminate confounding effects, we excluded the genus Clostridioides from our classification analysis. The immune marker-based classifier achieved macro-average AUC ~ 0.827 and micro-average AUC ~ 0.828 (Supplementary Figure 5a), which are quite comparable to the performance of microbiota-based classifier (Supplementary Figure 5b). Interestingly, integrating immune marker with gut microbiota showed much better classification performance (macro-average AUC ~ 0.926 and micro-average AUC ~ 0.869) (Supplementary Figure 5c).

Figure 5.

The performance of RF-based classification models based on various types of features in differentiating CDI from other groups. (a) CDI vs. Asymptomatic Carriage. (b) CDI vs. Non-CDI Diarrhea. (c) CDI vs. Non-CDI. For each classification task, we used different types of features: (1) the top-1 immune marker feature (based on mean decrease accuracy); (2) the top-1 genus feature; (3) all immune markers; (4) all genera; (5) integration of all immune markers and genera; (6) selected features from the set of all immune markers and genera. Error bars represent the standard errors of the means (SEM)

We further performed binary classifications to distinguish CDI subjects from Asymptomatic Carriage, Non-CDI Diarrhea, and Non-CDI subjects, using different feature types (Figure 5). The goal of this analysis was to assess whether any single taxon or immune marker could reliably differentiate CDI status. In the classification of CDI vs. Asymptomatic Carriage, we found that G-CSF and Moryella were the most important immune and microbial features based on mean decrease accuracy (MDA, the decrease in model accuracy from permuting the values in each feature), respectively (Supplementary Figure 6a-b). But the classification based on G-CSF (or Moryella) alone did not yield very high performance: mean AUC ~ 0.817 (or 0.701), respectively (Figure 5(a1, a2)). When we used all the immune markers (or all the genera) as features, we achieved mean AUC ~ 0.867 (or 0.805), respectively (Figure 5(a3, a4)). Interestingly, when we integrated all the host immune markers and gut microbial composition data together, we achieved a much higher performance with mean AUC ~ 0.900 (Figure 5(a5)). In order to select a subset of features that is as discriminatory as the whole set of features, we followed the “1-SE” rule (i.e., one chooses the model with fewest features such that its classification performance is less than one standard error away from that of the model with all the features), and selected the following 4 features: 2 bacterial genera (Moryella and Veillonella) and 2 immune markers (G-CSF and IL-6) in classifying CDI and Asymptomatic Carriage groups (Supplementary Figure 6:g-j). The RF classifier with those selected features displayed an outstanding classification performance, with mean AUC ~ 0.916 (Figure 5(a6)). Note that a significant negative correlation between Moryella and G-CSF was found in the Asymptomatic Carriage group (Figure 4c), which might contribute to the outstanding performance of the RF classifier with Moryella and G-CSF as selected features.

In the classification of CDI vs. Non-CDI Diarrhea groups, we found that G-CSF and [Eubacterium]_hallii_group are the top immune and microbial features, respectively (Supplementary Figure 6:c-d). But the classification based on G-CSF (or [Eubacterium]_hallii_group) alone did not perform very well: mean AUC ~ 0.747 (or ~ 0.630), respectively (Figure 5(b1, b2)). When we used all the immune marker (or all the microbial genera) as features, we achieved mean AUC ~ 0.851 (or ~ 0.884), respectively (Figure 5(b3, b4)). By integrating all features from both host immune marker and gut microbial genera, we further improved the classification performance to mean AUC ~ 0.918 (Figure 5(b5)). Following the “1-SE” rule, we selected the following 5 features: 3 genera: Enterococcus, Epulopiscium and [Eubacterium]_hallii_group; and 2 immune markers: G-CSF and IgA anti-toxin A (Supplementary Figure 6:h-k). The RF classifier with those selected features achieved mean AUC ~ 0.917 (Figure 5(b6)), which is quite comparable to that of using all the features. Note that Enterococcus was found to be significantly associated with G-CSF in the Non-CDI Diarrhea group (Figure 4b). This might partially explain the outstanding performance of the RF classifier with Enterococcus and G-CSF as selected features.

In the classification of CDI vs. Non-CDI groups, we found that G-CSF and Curvibacter are the top immune and microbial features, respectively (Supplementary Figure 6:e-f). Classification based on G-CSF (or Curvibacter) alone achieved mean AUC ~ 0.802 (or ~ 0.683), respectively (Figure 5(c1, c2)). When we used all the immune marker (or all the microbial genera) as features, we achieved mean AUC ~ 0.878 (or ~ 0.903), respectively (Figure 5(c3, c4)). Integrating all features from both host immune marker and gut microbial genera, we further improved the classification performance to mean AUC ~ 0.941 (Figure 5(c5)). Following the “1-SE” rule, we selected the following 10 features: 6 genera: Stenotrophomonas, Curvibacter, Enterobacter, Anaerobacillus, Fusobacterium and Veillonella; and 4 immune markers: G-CSF, IL-6, TNF-α and IgA anti-toxin B (Supplementary Figure 6:i-l). Classification with those well selected features achieved mean AUC ~ 0.929 (Figure 5(c6)).

Derive nonlinear interactions between gut microbiota and host immune markers using symbolic classification

As mentioned earlier, traditional correlation analysis cannot reveal any nonlinear interactions between gut microbiota and host immune markers. This fact and the outstanding classification results based on well-selected features prompt us to derive simple mathematical models to quantify the intricate interactions between gut microbiota and host immune markers. To achieve that, we leveraged symbolic classification (SC)[40,41], a genetic programming technique that automatically searches the space of mathematical expressions to find the model that best fits a given dataset. The fitness function in SC is a maximization function, and the number of generations is chosen based on the saturation of the fitness score (Supplementary SFigure 7). Using the same set of selected features and trained with the entire dataset, the SC model outperformed logistic regression (LR) in differentiating CDI (see Table 2).

Table 2.

Diagnostic scores derived from symbolic classification (SC) and logistic regression (LR). For each subject , we calculate his/her diagnostic score (or ) based on one of the following formulas derived from SC (or LR), respectively. For SC, the class of subject is CDI if ; or Asymptomatic Carriage (or Non-CDI Diarrhea, Non-CDI) if . For LR, the class of subject is CDI if ; or Asymptomatic Carriage (or Non-CDI Diarrhea, Non-CDI) if . Here, both and were learned from the entire dataset. Features used here include: GCSF; IgA_toxA; IgA_toxB; IL6; TNFα; Anaerobacillus; Curvibacter; Enterobacter; Enterococcus; Epulopiscium; [Eubacterium]_haillii_group; Fusobacterium; Moryella; Stenotrophomonas; Veillonella. In particular, for each classification task (regardless of using SC or LR), the following selected features were: (1) CDI vs. Asymptomatic Carriage: ,, and ; (2) CDI vs. Non-CDI Diarrhea: , , , , and ; (3) CDI vs. Non-CDI: , , , , , , , , and . Note that in the calculation of precision, recall and F1-score, we can treat either CDI (or Asymptomatic Carriage, Non-CDI Diarrhea, Non-CDI) as the true positive. Results shown in the parenthesis represent the latter case

Model	Diagnostic	Formula	Accuracy	Precision	Recall	F1-score
SC	CDI vs. Asymptomatic Carriage	fi=x_1x_15x13−0.2x_13+0.4+1.1x_1−0.1x_4−18.25	0.896	0.914 (0.840)	0.949 (0.75)	0.931(0.792)
	CDI vs. Non-CDI Diarrhea	fi=x_9x_20.5x_10−1+x_110.02x_11−x_1+x_21−10x_1−0.003x_9	0.900	0.946 (0.826)	0.897(0.905)	0.921(0.864)
	CDI vs. Non-CDI	fi=x_1x_3(0.2x_1x_5x_6x_14+0.04x_1x_7+0.3x_1x_15x84+x_1x_120.5x_7+x_1x_14+x_70.1x_4−x_6+x_14x_14−2	0.882	0.889 (0.878)	0.821 (0.927)	0.853 (0.902)
LR	CDI vs. Asymptomatic Carriage	logpi1−pi=0.66725−0.04442x_1+0.01022x_4+7.51484x_13−85.00213x_15	0.830	0.895 (0.667)	0.872 (0.714)	0.883 (0.690)
	CDI vs. Non-CDI Diarrhea	logpi1−pi=0.01974−0.002084x_1−0.02391x_2+1.895x_9−12740x_10+163.9x_11	0.800	0.814 (0.765)	0.897 (0.619)	0.854 (0.684)
	CDI vs. Non-CDI	logpi1−pi=2.122−0.01002x_1+0.01833x_3−0.006334x_4−0.009566x_5−4609x_6−8576x_7−40.75x_8−101.1x_12+32.84x_14−43.4x_15	0.813	0.841 (0.798)	0.679 (0.908)	0.752 (0.850)

Indeed, as shown in Table 2, we derived a simple SC model with selected features, reaching a very high accuracy (0.896) in distinguishing CDI subjects from Asymptomatic Carriage. Basically, for each subject , we calculate the diagnostic score that will be used for CDI diagnosis: the class of subject is CDI if ; Asymptomatic Carriage, if . Similarly, we derived a SC model with accuracy of 0.900 (or 0.882) in distinguishing CDI from Non-CDI Diarrhea (or Non-CDI) with the corresponding diagnostic score shown in Table 2. To ensure the SC models learned from the entire dataset are not overfitting, we performed cross-validation. With different training sets, SC will derive different mathematical formulas (i.e., diagnostic scores). However, those SC models learned from different training datasets demonstrated quite robust performance in terms of Accuracy (ratio of the total number of correct predictions and the total number of predictions), Precision (the number of positive class predictions that actually belong to the positive class), Recall (the number of positive class predictions made out of all positive examples in the dataset) and F1-score (a weighted average of the precision and recall) (see Supplementary Table 8). More importantly, even trained with less data, the SC models still outperformed LR models learned from the entire dataset.

As shown in Table 2, in the formulas of the diagnostic score, we colored the gut microbiota (or host immune marker) features in red (or blue), respectively. It is clearly seen that any potential interactions between gut microbiota and host immune markers are completely ignored in the formulas derived from LR. But for the formulas derived from SC, nonlinear interactions between gut microbiota and host immune markers can be clearly seen. We emphasize that those nonlinear interactions are not always pairwise. Those explicit interaction terms could inform further mechanistic studies to further reveal the role of intricate interactions between gut microbiota and host immune markers in CDI pathogenesis.

Discussion

Our studies suggest that the interactions between gut microbiota and host immune markers are very sensitive to the status of C. difficile colonization and infection. We demonstrated that incorporating both gut microbiome and host immune marker data into classification models can better distinguish CDI from other patient groups. Our classification models display robust diagnostic performance to differentiate CDI from Asymptomatic carriage. Using selected features to derive simple mathematic formulas we can explicitly quantify the interactions between gut microbiome and host immune markers, the two key components of CDI pathogenesis.

Consistent with previous studies[42-45], we found that the gut microbiota of CDI patients was characterized by lower Shannon diversity than that of the Control group. Interestingly, we observed an increased variation of both immune markers and gut microbial compositions in the CDI group with respect to other studied groups. This suggests that CDI is characterized by a significantly less stable microbiome and immune homeostasis. Our findings are in line with the Anna Karenina principle, which suggests that CDI linked changes in the microbiome and immune homeostasis are likely stochastic, leading to community instability[46-48].

We identified several candidate driver taxa (e.g., Desulfovibrio, Klebsiella, Streptococcus and Veillonella) that played a key role in driving the changes of microbial correlation networks between CDI and Asymptomatic Carriage (or Non-CDI Diarrhea, Non-CDI) groups. Among those driver taxa, Streptococcus has previously been shown to produce lactate thus impacting C. difficile TcdA expression to alleviate CDI[49]. Previous study indicated that Desulfovibrio has a pathogenic role in ulcerative colitis due to its ability to generate sulfides[50]. Klebsiella bacteria have been increasingly shown to develop antimicrobial resistance, most recently to the class of antibiotics known as carbapenems[51,52]. It is thus possible that the CDI pathogenesis is further enforced by the enrichment of antagonistic bacteria present in the gut microbiome of CDI subjects.

We developed classification models aimed at differentiating CDI status based on host immune markers and gut microbiome data. We were able to identify specific immune and microbial features that could accurately distinguish CDI subjects. In addition, most of the selected features identified by feature selection were also differentially abundant genera and differentially expressed immune markers. From the classification of CDI and Asymptomatic Carriage, we were able to select a few features with outstanding discriminability, including Veillonella and Moryella. Interestingly, a positive relationship between Veillonella and CDI has been identified in recent studies[53-56]. An important role for Veillonella in CDI is supported by the fact that Veillonella species were associated with low coprostanol levels that correlated strongly with CDI[53]. A similar negative relationship between Moryella species and CDI has previously been observed[57]. Enterococcus, a feature selected from the classification of CDI vs. Non-CDI Diarrhea, has been reported to be associated with CDI due to vancomycin resistance[58]. Consistent with the findings from previous reports[59,60], Epulopiscium was significantly enriched in the CDI group and played an important role in differentiating this comparison. Among those features selected from the classification of CDI and Non-CDI groups, Enterobacter and Fusobacterium have been considered as opportunistic pathogens involved in multiple diseases[61,62].

Machine learning method has the potential to identify biomarkers and aid in the diagnosis of many diseases. However, the learnt relationships between predictors and outcome are typically nontransparent, especially non-linear methods (i.e., decision tree learning)[63]. Classical logistic regression is one of the most common machine learning models in medicine. Yet, it fails to solve non-linear problems where there are multiple or non-linear decision boundaries[64]. Furthermore, the log odds scale in LR is hard to interpret[65]. Symbolic classification based on genetic programming is an automated technique to derive formulas from features[66]. Using the selected features from the random forests model, we demonstrated that the mathematical formulas automatically derived from symbolic classification have robust diagnostic accuracy to differentiate CDI patients from Asymptomatic Carriage (or Non-CDI Diarrhea, and Non-CDI groups). Specifically, symbolic classification provides explicit mathematic formulas as its output, which significantly improves the transparency of the learned relationship between predictors and outcomes. We previously demonstrated the potential clinical utility of a specific immunological biomarker (i.e., G-CSF) for diagnosis of CDI[36]. This study leverages the same unique and well-characterized study cohort, allowing us to study integrated host immune marker and microbial signatures associated with CDI. The fundamental differences between this study and the previous one are the clinical utilization of comprehensive immunological and microbial markers to explore the pathogenesis of CDI, and to generate clinical diagnostic models to detect CDI.

We acknowledge the following limitations of this study. First, the 16S rRNA gene sequencing may not have captured additional insights associated with CDI at the species or strain level. Second, observed associations do not prove causal relationship, and further studies are needed to validate the mechanism underlying the observed associations between these biomarkers and CDI. Our findings support the potential role of gut microbiome and immune markers in CDI and may serve as a starting point for future mechanistic studies. Third, diarrhea itself can affect the gut microbiome. When we compare the gut microbial compositions of groups with and without diarrhea, it is hard to disentangle the impact of diarrhea from that of disease on the gut microbiota, especially when the diarrhea is caused by the disease (e.g., in the CDI group). Completely resolving this limitation is out of the scope of the current study. Finally, further external validations of the classification models and derived formulas need to be performed on an additional cohort with same inclusion criteria as the current cohort.

In sum, utilizing this well-characterized cohort and leveraging machine learning tools, we proposed an effective computational framework to quantify the role of the intricate interactions between gut microbiota and host immune markers in CDI pathogenesis. We believe this framework offers the potential for microbiome-immune marker derived CDI diagnosis, as well as therapeutic strategies for its prevention and treatment.

Material and methods

Study cohort

The background and design of this cohort has been described in detail previously[67]. Basically, our clinical cohort consists of 243 well-characterized recruited participants, who were divided into four groups associated with different C. difficile infection/colonization statuses: (1) CDI (n = 112): Eligible patients were inpatients 18 years old with new-onset diarrhea, positive clinical stool NAAT result, and a decision to treat for CDI. The stool sample was captured as a discarded sample, and a discarded serum sample collected within 24 hours of that stool sample also captured. Patients were excluded if the diagnostic stool specimen was 72 hours old, if they had received CDI treatment for 24 hours prior to stool collection, or if they had a colostomy. Assessment for the presence of diarrhea included review of nursing input/output logs for number and consistency of stools, consultation with treating clinicians, and detailed chart review (requiring mention of “diarrhea”, “loose stools”, and/or increased frequency, in notes written by multiple providers). Patients for whom there was any doubt about the presence of diarrhea, or who had chronic diarrhea, were excluded. (2) Asymptomatic Carriage (n = 40): Eligible patients were inpatients 18 years old, admitted for at least 72 hours, who had received at least one dose of an antibiotics within the past 7 days, and did not have diarrhea in the 48 hours prior to stool specimen submission. Patients with 2 or more loose stools within 24 hours were excluded; patients who had 1 loose stool were included only if they had recently received a laxative. Patients were excluded if they had a colostomy; received oral or intravenous metronidazole, oral vancomycin, oral rifaximin, and/or oral fidaxomicin for 24 hours within the prior 7 days; had been diagnosed with CDI in the past 6 months; or had tested negative for C. difficile within the past 7 days. Stool specimen were collected prospectively under verbal informed consent. A discarded serum sample from within 24 hours of the stool specimen were also captured. NAAT (Xpert C. difficile/Epi) was performed on all samples, and positive samples retained as the Asymptomatic Carriage cohort. (3) Non-CDI Diarrhea (n = 44): patients with diarrhea (confirmed using the same definition used for the CDI cohort) but had NAAT-negative stool on clinical C. difficile testing; (4) Control (n = 47): patients without diarrhea who had screened as eligible for the Asymptomatic Carriage cohort but were NAAT-negative on research stool testing. In our previous study[36], the four groups were named as (1) CDI-NAAT = CDI; (2) Carrier-NAAT = Asymptomatic Carriage; (3) diarrhea NAAT-negative = Non-CDI Diarrhea and (4) no Diarrhea NAAT-Negative = Control. In this work, for simplicity we used the simpler and more clearly descriptive titles.

Serum immune marker measurement

The measurement of host serum cytokines concentrations of IL-2, IL-4, IL-6, IL-8, IL-10, IL-13, IL-15, IL-1β, G-CSF, MCP-1, VEGF-A, and TNF-α was performed using a Milliplex magnetic bead kit and Luminex analyzer (MAGPIX) (Millipore Sigma, Inc., Burlington, MA) as per the manufacturer’s instructions. Purified toxin A and B were separately prepared from C. difficile strain VPI 10463 (American Type Culture Collection 43255-FZ, Manassas, VA). Serum antibody (IgA, IgG, and IgM) levels against C. difficile toxins A and B were measured by semi-quantitative enzyme-linked immunosorbent assay (ELISA). All the experimental details have been reported previously[36,67].

Fecal DNA extraction and bacterial 16S rRNA sequencing data analysis

Stool DNA was extracted using the DNeasy PowerSoil Pro Kit (Qiagen, cat# 12888–100) in a QiaCube-automated DNA extraction system (Qiagen) according to instructions. Briefly, 250 mg stool was transferred into a PowerBead Pro Tube provided with the kit and 200 ug RNaseA and 800 μl of CD1 solution were added. Tubes were vortexed briefly, transferred into an adapter, and then vortexed at maximum speed for 10 min. Tubes were centrifuged at 15,000 xg for 1 min and about 500–600 μl supernatant was used for DNA extraction according to instructions. DNA were eluted in 70 μl elution solution C6 and stored at −80°C until use. 16S rRNA microbiome characterization was performed by sequencing the V4 region of the 16S rRNA gene using the Illumina MiSeq[68]. Each sample was amplified using a barcoded primer, which yielded a unique sequence identifier tagged onto each individual sample library. Illumina-based sequencing yielded greater than 15,000 reads per sample. CLC Genomics Workbench version 12 (Qiagen) was used for OTU clustering and generation of abundance tables. Analyses were performed using the tutorial “OTU Clustering Step by Step” updated September 2, 2019 and available at: https://resources.qiagenbioinformatics.com/tutorials/OTU_Clustering_Steps.pdf

Microbial diversity and differential abundance analysis

The diversity measures and permutational multivariate analysis of variance (PERMANOVA, compare groups of objects and test the null hypothesis that the centroids and dispersion of the groups as defined by measure space are equivalent for all groups) were calculated using the vegan package in R (see Supporting methods for details). For differential abundance analysis, we used ANCOM[37] (analysis of composition of microbiomes), with a Benjamini–Hochberg correction at 5% level of significance, and adjusted for age and sex. The Mann–Whitney U test was used to compare the difference of immune marker levels between different groups.

Microbial correlation network and microbiome-immune association analysis

The microbial correlation networks were constructed using SparCC[38] (sparse correlations for compositional data, https://github.com/luispedro/sparcc) (see Supporting methods for details). We also used the NetShift[39] (https://web.rniapps.net/netshift) to identify potential “driver” taxa that are responsible for the differences of microbial correlations between the CDI and Asymptomatic Carriage (or between Non-CDI Diarrhea and Non-CDI) networks. The driver taxa were identified based on their neighbor shift (NESH) scores, and betweenness centrality (BC) measures in the two networks. Associations between the gut microbiota and host immune markers were quantified by Spearman correlation coefficients in combination with Benjamini-Hochberg FDR correction to account for multiple hypothesis testing (significance threshold α ≤ 0.05). All included genera were required to be detected in ≥50% of all samples in each group.

Classification with Random Forests model

To build a classification model capable of testing the overall contribution of immunological or microbial data in distinguishing the CDI status, we developed a multi-class random forest (RF) classifier. The data is split into a training set and a test set, with 70% of the data forming the training data and the remaining 30% forming the test set. The performance of the multi-class model was measured by micro-average (aggregate the contributions of all classes to compute the average metric) and macro-average (compute the metric independently for each class and then take the average) AUC (the area under the receiver operating characteristic (ROC) curve, the measure of the ability of a classifier to distinguish between classes and is used as a summary of the ROC curve). To determine whether more specific host immune markers or gut microbial taxa could differentiate CDI subjects from Asymptomatic Carriage, Non-CDI Diarrhea and Non-CDI groups, we constructed the binary classifiers based on RF models with integrated immune markers and microbiome data (see Supporting methods for details).

Symbolic classification with genetic programming

We employed Karoo GP[69], a genetic programming application suite written in Python that support both symbolic regression (SR) and symbolic classification (SC) analysis, to derive simple formulas for CDI diagnosis. Due to the different training sets, SC will derive different formulas, but their classification performances are quite comparable (S8 Table). The formulas shown in Table 2 were derived based on the whole dataset (for details see supplementary methods). To demonstrate the advantage of SC, for each classification task (i.e., CDI vs. Asymptomatic Carriage, CDI vs. Non-CDI Diarrhea, and CDI vs. Non-CDI), we also performed logistic regression (LR) using the same set of selected features as used in SC (Table 2) (see Supporting methods for details).

Click here for additional data file.