Predictive regulatory and metabolic network models for systems analysis of Clostridioides difficile

Abstract

We present predictive models for comprehensive systems analysis of Clostridioides difficile, the etiology of pseudomembranous colitis. By leveraging 151 published transcriptomes, we generated an EGRIN model that organizes 90% of C. difficile genes into a transcriptional regulatory network of 297 co-regulated modules, implicating genes in sporulation, carbohydrate transport, and metabolism. By advancing a metabolic model through addition and curation of metabolic reactions including nutrient uptake, we discovered 14 amino acids, diverse carbohydrates, and 10 metabolic genes as essential for C. difficile growth in the intestinal environment. Finally, we developed a PRIME model to uncover how EGRIN-inferred combinatorial gene regulation by transcription factors, such as CcpA and CodY, modulates essential metabolic processes to enable C. difficile growth relative to commensal colonization. The C. difficile interactive web portal provides access to these model resources to support collaborative systems-level studies of context-specific virulence mechanisms in C. difficile.

Keywords: Biological networks; Clostridioides difficile; Commensals; EGRIN; Host-pathogen interactions; In vivo adaptive response; Integrated Regulatory and Metabolic Network Model; PRIME; Web Portal.

Figure 1.. Inference pipeline and general properties of the resulting Environment and Gene Regulatory Influence Network (EGRIN) model of C. difficile.

(A) Framework used to build the EGRIN model. (B) Distribution of residual values for the 406 detected co-regulated gene modules. 297 modules with residual ≤ 0.55 (shown in purple) were labelled as high quality. (C) Distribution of gene count for the high quality gene modules. (D) Coverage of all genes (4,018), the subset of metabolic genes (1,030) and TFs (309) by EGRIN modules for different residual thresholds. The red dashed line indicates the 0.55 residual cutoff.

Figure 2.. The EGRIN model of C. difficile recapitulates known biology of the pathogen.

(A) Co-regulated gene modules are enriched with functional terms derived from expert curated annotation of the C. difficile genome (Girinathan et al., 2021). The pie chart shows terms over-represented in three or more modules. Number of modules associated with each functional term is shown in parenthesis. (B) Enriched EGRIN modules among nine (out of 13) manually-defined and experimentally supported TF regulons (Table S2). (C) EGRIN identified the known DNA binding motif of CodY (Dineen et al., 2010). (D) EGRIN also identified the known DNA binding motif of SigL (Soutourina et al., 2020). (E) The EGRIN model recapitulated the previously reported influence of CodY on tcdA expression. The module #182 contains tcdA, it is enriched with the CodY regulon and contains a GRE (shown in panel C) similar to the experimentally determined CodY motif. (F) The EGRIN model captured the interaction between toxin expression and sporulation via module #397 that contains tcdB and is enriched with genes regulated by sporulation-related transcriptional regulators. (G) Expression profiles of tcdAB and tcdR (positive regulator of the PaLoc). Log2 ratios were computed with respect to dataset-specific references (Table S1), and grouped by condition blocks. Boxes cover the 25^th-75^th percentile range (median indicated by horizontal black line).

Figure 3.. The EGRIN model offers insights on potential functions of uncharacterized genes of C. difficile.

Hypotheses regarding the functions of 48 uncharacterized genes were generated based on their membership in high quality EGRIN modules significantly enriched with specific functional terms. (A) Barplot with the number of unknown genes associated with each functional term (from the C. difficile genome annotation in (Girinathan et al., 2021)). (B) The involvement of 10 uncharacterized genes in sporulation was supported by their significant downregulation (with respect to a wild-type control) in single deletion strains of sporulation regulators (Table S1). Boxes cover the 25^th-75^th percentile range (median is indicated by horizontal black line).

Figure 4.. The EGRIN model identifies TFs driving the in vivo response of C. difficile when interacting with gut commensals.

(A) Expression profile of module #48 across the transcriptional compendium used to build EGRIN. Each box represents log2 fold-change (with respect to a dataset-specific reference, Table S1) of members of the module in a single transcriptome. Transcriptomes are color coded according to their membership in the 11 condition blocks in the compendium and ranked based on their median log2 fold-change. The ‘Early infection’ condition block is statistically over-represented in the 20% of highest (indicated with the 80% dashed line) transcriptomes median log2 fold-change. (B) Expression profile of module #158. Only one condition block (‘In vivo vs. in vitro’) was found in the 20% of highest fold-change. (C) EGRIN modules enriched with genes differentially expressed (absolute log2 fold-change > 1 and adjusted p-value < 0.05) in C. difficile mono-colonized mice at 24-h vs 20-h of infection. X-axis shows module IDs. Modules were annotated according to their functional enrichment and overlap with manually curated TF regulons (Table S2). (D) Enriched EGRIN modules in C. sardiniense+C. difficile co-colonized mice vs C. difficile mono-colonized mice at 24-h of infection. Due to space constraint, only abbreviations of functional terms not shown in other panels are displayed. Modules enriched with Spo0A and sporulation sigma factors are indicated with ‘Spo0A/-’. (E) Enriched EGRIN modules in P. bifermentans+C. difficile co-colonized mice vs C. difficile mono-colonized mice at 24-h of infection. (F) Enriched EGRIN modules in P. bifermentans+C. difficile co-colonized mice vs C. sardiniense+C. difficile co-colonized mice at 24-h of infection. For all comparisons, only modules with absolute median fold-changes ≥ 0.5, and enriched with TF regulons or functional categories are displayed.

Figure 5.. Metabolic model predictions.

(A) Details of the in vitro metabolic models (icdf834 and icdf836) of C. difficile 630 (Kashaf et al., 2017). General properties of the icdf836 model (generated in this study) after adding the required in vivo exchanges, transports and reactions are shown. The ‘*’ indicates that there were two duplicated genes in the icdf834 model, reducing the total number of genes to 832. (B) ROC curves showing the accuracy of icdf834- and icdf836-predicted gene essentiality in nutrient rich medium evaluated against a Tn-seq functional screen (Dembek et al., 2015). (C) Details of the in vivo (mono-colonized) model derived using the GIMME algorithm (Becker and Palsson, 2008) on the icdf836 model where only the active reactions are included from in vivo transcriptome. (D) Venn diagram showing the number of model-predicted essential genes for growth of C. difficile 630 in vitro vs in vivo. Only genes predicted as essential for in vivo (mono-colonized) growth were considered in the analysis.

Figure 6.. PRIME-estimated metabolic fluxes and gene essentiality in mono- and P. bifermentans co-colonized conditions.

(A) BioTapestry visualization of in vivo gene regulatory network for C. difficile 630: 79 essential genes with Inferelator-predicted transcriptional regulators are shown. The genes and regulators shown as five digit numbers represent the nomenclature preceded by ‘CD630_’. All regulators with less than two essential targets were combined in the ‘Others’ meta-regulator. The ‘#’ symbol indicates the 10 genes predicted as essential in the mono-colonized condition but not essential in vitro. (B) PRIME-estimated metabolic fluxes for each reaction in the mono- and co-colonized conditions. Z-score transformation was independently applied to each reaction and condition. The ‘*’ symbol indicates that there is some contention about the existence of this pathway in C. difficile and it may need to be revised in future models.