The Institutional Review Board at the University of Virginia approved this study (IRB-HSR# 21646) with waiver of consent, as samples were de-identified remnants. Diarrheal samples from CDI+ (n = 20) and CDI− (n = 20) patients collected by the University of Virginia Medical Center between September 2019 and August 2020. CDI was diagnosed based on the presence of a conserved sequence of the tcdB gene, as determined by Xpert® C. difficile assay (Cepheid, Sunnyvale, CA, United States). To be included in the study, patients had to be at least 18 years old. Samples were stored in a −80°C freezer after collection until further processing.

RNA was extracted from samples using the ZymoBIOMICS RNA Miniprep Kit (Zymo Research, Irvine, CA, United States) according to the manufacturer’s protocol with the following exceptions: 1 volume of lysis buffer was used, the optional DNase I treatment was completed, and extracts were eluted with 50 uL of DNase/RNase free water. After extraction, quantification was conducted using an Invitrogen Qubit 4 Fluorometer and Qubit RNA High Sensitivity Assay Kit (ThermoFisher Scientific, Waltham, MA, United States).

RNA was then processed to make MT libraries using the NEBNext Ultra II RNA Library Prep Kit for use with Illumina (New England Biolabs, Ipswich, MA, United States), which utilizes a random priming approach for reverse transcription. The protocol specifically for use with purified mRNA was followed. Depletion of rRNA was not performed.

Libraries were quality checked using an Agilent 2100 Bioanalyzer and DNA High Sensitivity kit. Results from these steps were used to pool samples in an equimolar ratio. The pool was then gel purified using a 2% agarose gel and the Qiagen QIAquick gel extraction kit (Qiagen, Germantown, MD, United States). Following purification, the pool was sequenced using an Illumina NextSeq 550 platform to produce 2 × 150 bp reads by Wright Labs LLC (Huntingdon, PA, United States).

In total, 40 samples were subjected to MT sequencing, yielding 759,555,378 raw sequences, and after quality filtration, 639,700,992 sequences remained. Of those samples, 39 (CDI + =20, CDI− = 19) had enough data for downstream analysis after quality filtering (range 942,194–82,265,454 sequences per sample, Figure 1).

Pielou’s evenness values within the expressed genes dataset tended to be higher in CDI− samples compared to the CDI+ cohort, though the difference was not significant (Table 1). Overall, alpha diversity was variable among samples but did not differ significantly between CDI+ and CDI− for the active microbial taxa or expressed genes datasets based on either the Pielou’s evenness or Observed Features metrics (Figure 2).

Significant clustering based on CDI status was observed for both active taxa and expressed gene datasets (PERMANOVA, p = 0.005 and p = 0.013, respectively; Figure 3). Subsequent LEfSe analysis identified the microbial taxa and KEGG Orthologs (KOs) driving differential clustering.

The commensal bacteria Clostridium butyricum, Lacticaseibacillus rhamnosus, and Roseburia intestinalis were the most differentially active microbial species in CDI− (Figure 4). However, variability in the expression of all three taxa within CDI− was relatively high (Supplementary Table S1). In contrast, though its corresponding LDA score was lower, Lactiplantibacillus plantarum was more consistently active in the CDI− samples (Supplementary Table S1).

Multiple butyrate producers were identified as being more active in CDI+, including Anaerostipes hadrus (Allen-Vercoe et al., 2012), Coproccocus comes (Louis and Flint, 2009), and Roseburia hominis (Louis and Flint, 2009), and all three were found to express genes within the Butanoate Metabolism pathway (ko00650) in this study. However, the taxa Clostridioides difficile, Enterocloster bolteae, and Ruminococcus gnavus were the three most differentially active species in CDI+ (Kruskal-Wallis, p ≤ 0.05, log(LDA) ≥ 3.0), with C. difficile being the most differential (Supplementary Figure S3). Clostridium scindens was also identified as being more active in CDI+.

When our microbial taxa table was subsetted to include only fungal taxa, the family Saccharomycetaceae was more active in CDI− samples (Figure 5). In contrast, 31 fungal taxa were more active in CDI+ (Supplementary Table S2). Several differential genes (Supplementary Table S3) and pathways (Supplementary Table S4) were also identified within the fungi functional gene datasets. However, we recovered relatively low numbers of fungal sequences (Supplementary Figure S4).

More differential genes were identified within the CDI+ cohort in the total expressed genes dataset (CDI+ genes = 99, CDI− genes = 31), including several that were especially of interest due to their involvement with biofilm formation or sporulation and their expression by C. difficile (Table 2). However, a greater number of pathways were expressed associated as differential in CDI− (Supplementary Table S5). Interestingly, 22 of the genes more expressed in CDI+ code for ribosomal subunits (Supplementary Table S6). The most significantly different contributors to those genes were taxa significantly more active in CDI+ (Supplementary Table S7).

Clostridioides difficile was one of the taxa found to be significantly more critical to the expression of differential rRNA genes in CDI+ (Supplementary Table S7). It was also found to be a significantly more important contributor to the expression of three differential KOs related to spore formation in CDI+ (Supplementary Table S8). Likewise, Enterocloster bolteae contributed to the differential expression of two of those KOs (Supplementary Table S8). Similarly, though it was not associated with differential expression of those three spore formation KOs, Ruminoccocus gnavus expressed genes relating to spore formation as well. Clostridioides difficile also constituted a greater proportion of the expression of genes relating to biofilm formation in CDI+, specifically K03073 (SecE), K03075 (SecG), and K03666 (hfq; Supplementary Tables S9, S10). Clostridioides difficile and R. gnavus were also significantly more important to the expression of genes within the differential Flagellar Assembly (ko02040) in CDI+ (Supplementary Table S11). Additionally, C. difficile was identified as expressing the tcdB gene in three CDI+ samples (13, 19, and 5), albeit very few sequences were associated with both (1, 2, and 5 sequences respectively).

Random forest modeling was used to help evaluate the consistency of features identified by LEfSe for differentiating samples based on CDI status, regardless of possible confounding factors. The performance of those models indicates that both microbial taxa and expressed genes identified by LEfSe differed consistently between CDI+ and CDI–. Both random forest models using differential features (genes or species) had average accuracies of 84.2% (Table 3). As expected, the differential species model excluding C. difficile performed worse, but it still had a 9% greater average accuracy than the random forest model that was generated using all species. Likewise, the genes model had a 7% increase in average accuracy between the dataset containing all genes and the dataset containing only genes identified as differentially expressed by LEfSe. Additionally, eight of the 10 best predictors for the full species dataset were also identified as being predictive of CDI status by LEfSe (Figure 6). The best predictors for the other models are shown in Supplementary Figures S5–S9.

Co-occurrence network analysis helped us model potential transkingdom interactions between bacteria and fungi with respect to CDI status. The networks revealed negative relationships between fungi and bacterial taxa in both CDI− and CDI+ networks. In the CDI− network, 18 fungi were present, with Nakaseomyces glabratus (formerly Candida glabrata) having negative correlations with 57 taxa, 55 of which were bacteria (Figure 7). In contrast, N. glabratus was also entirely absent in the CDI+ network (Figure 7). Furthermore, the CDI+ network only had two fungi, Saccharomyces cerevisiae and Debaryomyces hansenii.

Several predicted metabolites (n = 16) were found to be significantly more abundant (Kruskal-Wallis, p ≤ 0.05; Supplementary Table S12) in CDI+ patients, including caproic acid (Figure 8). In addition, a secondary bile acid, lithocholic acid, was also predicted to be more abundant in CDI+. Fewer predicted metabolites (n = 6) were associated with CDI–.

Although previous studies have noted significant differences in alpha diversity, those studies utilized 16S rRNA and ITS data to characterize microbial communities (Antharam et al., 2013; Lamendella et al., 2018; Zuo et al., 2018). To the best of our knowledge, this study is the first to examine possible alpha diversity differences based on CDI status using MT data. However, these results are aligned with another study that compared CDI− diarrheal communities to CDI+ communities that also did not observe significant differences between those groups (Antharam et al., 2013).

Still, similar to other CDI microbiome studies (Lamendella et al., 2018; Zuo et al., 2018; Stewart et al., 2019) overall community composition and gene expression differed significantly based on CDI status. The recapitulation of this finding is notable because these results indicate that both the active community members and the functions they express differ with respect to CDI status.

The species most strongly associated with CDI− samples are beneficial to human health. Specifically, both Clostridium butyricum and Roseburia intestinalis are butyrate producers (Hayashi et al., 2021; Nie et al., 2021). Similarly, the use of Lacticaseibacillus rhamnosus as a probiotic has been associated with increased butyrate production (Berni Canani et al., 2016; Lin et al., 2020; Carucci et al., 2022), and Clostridium and Roseburia genera are contained in a newly FDA-approved fecal microbiota capsule (SER-109) for use against recurrent CDI (Khanna et al., 2016). In the context of CDI, butyrate has been found to reduce inflammation and improve intestinal barrier function following C. difficile colonization (Fachi et al., 2019). Roseburia intestinalis was most abundant in the healthy cohort in a study examining the gut microbiome of children based on CDI and hospitalization (IC vs. non-IC) status (Mohandas et al., 2020). Additionally, multiple studies have found C. butryicum promotes resistance to C. difficile infection (Hagihara et al., 2021; Hayashi et al., 2021), and it has been investigated as a possible treatment for CDI (Oka et al., 2018; Lee et al., 2022). Likewise, studies have found evidence administering Lactiplantibacillus plantarum as a probiotic may help prevent CDI (Klarin et al., 2008; Kujawa-Szewieczek et al., 2015; Dudzicz et al., 2018). Therefore, the higher activity of C. butryicum and L. plantarum in CDI− samples could have helped protect those patients against acquiring C. difficile.

Multiple butyrate producers were also identified as more active in CDI+. Although previous research has shown that butyrate was significantly reduced in CDI cohorts compared to healthy controls (Pensinger et al., 2023), another study did not find a significant difference in butyrate producers (based on 16S rRNA analysis) in their diarrheal controls compared to CDI, while their healthy cohort had significantly more butyrate producers compared to their CDI and diarrheal control groups (Antharam et al., 2013). Therefore, since both CDI+ and CDI− groups in this study were comprised of diarrheal samples, differential gene expression attributed to butyrate producers in both is consistent with previous findings.

Interestingly, Clostridium scindens was also identified as being more active in CDI+. Clostridium scindens had previously been associated with resistance to CDI as a result of producing secondary bile acids (Buffie et al., 2015; Greathouse et al., 2015), but resistance may be contingent on the resulting concentration of deoxycholate (Dubois et al., 2019). Still, C. scindens has also been found to inhibit C. difficile growth independent of bile acid production (Aguirre et al., 2021). Instead, the authors of that study suggested resistance could be due to competition for nutrients with C. difficile (Aguirre et al., 2021). Therefore, though unexpected, if nutrient conditions were conducive to C. difficile growth in our CDI+ patients, increased activity of C. scindens could have also been favored.

Besides Clostridioides difficile, the other species most strongly associated with CDI+ are also known pathogens. Enterocloster bolteae has been previously noted to be an opportunistic pathogen and are also correlated with the expression of inflammatory genes (Pandit et al., 2021). Similarly, Ruminococcus gnavus has been found to produce a glucorhamnan that causes inflammation (Henke et al., 2019). Additionally, Ruminoccocus gnavus was previously associated with CDI in another study that found several operational taxonomic units (OTUs) belonging to this species were more abundant in their CDI IBD group compared to the non-CDI IBD group and the healthy subjects (Sokol et al., 2017). Notably, Ruminococcus gnavus is a mucin degrader, and a previous study found C. difficile grown in the presence of MUC2 with mucin-degrading taxa also present had increased expression of genes associated with flagella (Engevik et al., 2021).

The larger number of fungal taxa associated with CDI+ aligns with previous studies that found fungal taxa were more abundant in CDI patients (Sangster et al., 2016; Lamendella et al., 2018; Stewart et al., 2019). Therefore, this study provides further evidence that the mycobiome may have an important impact on the pathogenesis of CDI. For instance, Encephalitozoon cuniculi, a species of fungi enriched in CDI+ patients, is a known opportunistic pathogen (Moretto et al., 2015). Conversely, the use of Saccharomyces boulardii, which is a part of the family Saccharomycetaceae (more active in CDI− patients), as a probiotic seems to reduce CDI-associated diarrhea (Goldenberg et al., 2017).

The relatively large number of differential KEGG Orthologs associated with ribosomal subunits was likely due to the taxa expressing them being more transcriptionally active in CDI+ samples because of other reasons, rather than a driving reason for increased activity themselves. Higher expression of genes coding for ribosomes has previously been correlated with an increase in activity for the microbes expressing them, though not consistently (Blazewicz et al., 2013). Still, the relatively greater contribution of CDI+ differential taxa to these genes provides another line of evidence that they were more active in CDI+.

Clostridioides difficile is known to form spores during infection (Paredes-Sabja et al., 2014), which fits with the observation of it as a significantly more important contributor to the expression of three differential KOs related to spore formation in CDI+ in this study. Likewise, Enterocloster bolteae can form spores (Magdy Wasfy et al., 2023) and contributed to the differential expression of two of those KOs (Supplementary Table S8). Similarly, though Ruminoccocus gnavus had historically been considered to be non-spore forming, more recent work has demonstrated that it can form spores under certain conditions (Crost et al., 2023). Therefore, it is possible that members of all three taxa could form spores in response to antibiotics, which may help explain their presence in CDI+ samples, as antibiotics are often used to treat CDI.

Both SecE and SecG have previously been proposed to be involved with C. difficile biofilm (Poquet et al., 2018). SecE and SecG have also been associated with biofilm formation in Staphylococcus aureus (Resch et al., 2005) and Bifidobacterium longum (Zhang et al., 2023) respectively. Therefore, although we did not recapture our previous finding of differential biofilm pathways (Stewart et al., 2019), the observed overall increased expression of hfq, SecE and SecG and their expression by C. difficile in this study, as well as other genes relating to biofilm formation, provides further in vivo evidence that biofilm formation is involved in the pathogenesis of CDI. However, hfq, SecE, and SecG are also components of other pathways, and of the 21 KOs expressed by C. difficile that are associated with at least one biofilm formation pathway, only two are associated solely with biofilm formation (Supplementary Table S10). Therefore, future work is needed to confirm in vivo biofilm formation. Still, biofilm formation, if occurring, could help protect C. difficile from antibiotics and consequently, contribute to CDI’s high reoccurrence rate.

Two of the four pathways more expressed in CDI+ [Ethylbenzene Degradation (ko00642) and Flagellar Assembly (ko02040)] were also identified as being differentially expressed in our prior study (Stewart et al., 2019). Notably, the production of flagella has been associated with pathogen colonization (Rossez et al., 2015), and C. difficile, E. bolteae, and Parabacteroides distasonis were significantly more important contributors to genes within that pathway in CDI+ (Supplementary Table S11). Like E. bolteae, P. distasonis is also an opportunistic pathogen (Ezeji et al., 2021). Therefore, the enrichment of the flagellar assembly pathway here was driven by the increased activity of potentially pathogenic bacteria.

Correlations between bacteria and fungi within CDI+ samples have been observed in other studies (Lamendella et al., 2018; Stewart et al., 2019). However, Lamendella et al., 2018 found no correlations between fungi and bacteria in CDI− (Lamendella et al., 2018), and Stewart et al., 2019 had only one unidentified fungi with a significant correlation to an unidentified Bacteroides OTU (Stewart et al., 2019). Networks in both were based on ITS data, so the use of MT data here instead is a likely contributor to the differences in the prevalence of fungi in the CDI− network. Still, despite evidence of fungi’s importance to CDI pathogenesis, co-occurrence analysis itself did not identify a significant correlation between fungal species and C. difficile in this study or in either of those previous studies (Lamendella et al., 2018; Stewart et al., 2019).

Caproic acid was previously detected in 41% of CDI+ fecal samples according to a gas chromatography study used to predict C. difficile (Levett, 1984). In contrast, lithocholic acid, also associated with CDI+ in this study, has previously been found to inhibit the growth of C. difficile (Kang et al., 2019). However, it is possible that lithocholic acid could have been more abundant in our CDI+ samples but still below a concentration that would have been lethal to C. difficile, and if so, then it is also possible that lithocholic acid could have induced biofilm formation, as was the case with deoxycholate in a previous study (Dubois et al., 2019). This possibility is supported by our observance of C. difficile expressing multiple genes related to biofilm formation.

Still, it should be noted that tools like MelonnPan are predictive metabolic modeling tools, and it is unknown whether the lithocholic acid metabolite was truly more abundant in our CDI+ samples. Consequently, directly measuring metabolites of interest like lithocholic acid and caproic acid, in conjunction with microbiome analysis, would be valuable to further understanding how the abundance of metabolites previously associated with CDI relates to microbiome composition and function.

Unfortunately, no information about antibiotic usage or clinical history was reported for the individuals used in this study. Previous work has shown that antibiotic usage in conjunction with CDI has a large impact on the host’s overall microbial community (Lamendella et al., 2018). Specifically, in that study, the fidaxomicin and metronidazole groups had significant differences in alpha and beta diversity between CDI+ and CDI−, while the vancomycin group did not (Lamendella et al., 2018). Additionally, different taxa were detected as being significantly differential between CDI+/− depending on the antibiotic class (Lamendella et al., 2018). Still, another study did not find a significant difference in alpha diversity within their CDI+ and CDI− groups when split by antibiotics exposure (Antharam et al., 2013). However, it is possible that a differential usage of antibiotics could have impacted our results. Consequently, we limited much of our discussion to comparing the results obtained herein to previous CDI studies.

The failure to consistently capture sequences associated with both C. difficile and the gene used for diagnosis (tcdB) could also indicate our sequencing depth was not sufficient to fully capture all expressed genes in our samples. Additionally, due to the cost associated with MT, our sample size was smaller than some previous studies. Though, post hoc power analysis using SparseDossa2 showed we had reasonable (>80%) power to detect changes corresponding to at least 7.5 log2 fold change in the active species dataset and 4.5 log2 fold in the expressed genes dataset with LEfSe.

Despite antibiotic usage being a potentially large confounding variable, this work has contributed great value in characterizing active microbes and their associated functions in the context of CDI. To the best of our knowledge, MT analysis had only been applied to human samples associated with this disease state in a previous study by our group, in which it was used to identify differential genes and pathways (Stewart et al., 2019). Therefore, this study represents one of the first applications of MT to characterize differences in overall community diversity between CDI+ and CDI− cohorts. Although alpha diversity did not differ significantly, beta diversity for active species and expressed genes differentiated CDI+ and CDI− cohorts.

Subsequent differential feature analysis identified the specific genes and taxa driving those overall differences, revealing significant increases in the activity of several beneficial bacteria in CDI− and increased activity of taxa associated with inflammation in CDI+. Although the previous finding of differential biofilm formation pathways (Stewart et al., 2019) was not recaptured, several specific genes associated with biofilm formation were identified as differentially expressed in CDI+, and C. difficile was identified as a contributor to them, serving as additional evidence for the role of biofilm formation in CDI pathogenesis. Biofilm formation during infection would help explain CDI’s high recurrence rate, and if that is the case, treatments designed to target biofilms specifically could be effective for preventing reoccurrence.

Additionally, multiple fungal taxa were identified as being more transcriptionally active in CDI patients, which provides additional support for the potential importance of the mycobiome to CDI. Notably, E. cuniculi can cause inflammation, and increased inflammation due to it or other fungi could be one mechanism by which severe CDI could occur irrespective of bacterial toxin production. However, coverage for both E. cuniculi and fungi overall in this study tended to be low. Determining if specific fungal taxa are important for preventing or enabling CDI could facilitate the use of additional treatment methods. The low proportion of CDI+ samples that had tcdB identified within them suggests that greater sequencing coverage would also be valuable for characterizing the functional capabilities of these microbial communities. Therefore, future work should include deeper metatranscriptomics sequencing in order to increase coverage of the entire microbial community.

Increasing fungal sequences specifically, either through overall increased sequencing depth or targeting sequencing, would also aid better understanding if fungal activity contributes to the pathogenesis of CDI. Likewise, the development of methods to increase the proportion of fungal nucleic acids recovered would similarly help facilitate investigations of fungi’s role in CDI. Additionally, future iterations of the CDI status model could theoretically be trained to differentiate non-C. difficile diarrhea with colonization from clinically-relevant C. difficile infection, as a means to reduce unnecessary anti-C. difficile antibiotic use (Polage et al., 2015). Metabolomics data would also be invaluable for further exploring the role that previously identified metabolites of interest play in CDI.