key: cord-0077409-hfhybrir
authors: Sudhakar, Padhmanand; Andrighetti, Tahila; Verstockt, Sare; Caenepeel, Clara; Ferrante, Marc; Sabino, João; Verstockt, Bram; Vermeire, Severine
title: Integrated analysis of microbe-host interactions in Crohn’s disease reveals potential mechanisms of microbial proteins on host gene expression
date: 2022-02-22
journal: iScience
DOI: 10.1016/j.isci.2022.103963
sha: 47996231f0129ab41ecb4ec471051c874ad577df
doc_id: 77409
cord_uid: hfhybrir

Inflammatory responses of the intestinal epithelial barrier in patients with Crohn’s disease (CD), a chronic inflammatory bowel disease (IBD), are associated with gut microbial alterations. At a community level, there is scarce mechanistic evidence on the effects of gut microbial alterations on host mucosal barrier responses. We used a computational microbe-host interaction prediction framework based on network diffusion and systems biology to integrate publicly available paired gut microbial and intestinal gene expression datasets. The ileal signaling network potentially modulated by the microbiota was enriched with immune-related pathways such as those associated with IL-4, IL-2, IL-13, NFkB, and toll-like receptors. We identified bacterial proteins eliciting post-translational modifications on host receptors, resulting in the de-repression of pro-inflammatory cytokines via critical hub proteins such as NFkB. The signaling networks were over-represented with CD associated genes and CD drug targets. Using datasets generated from our validation cohorts, we confirmed some of the results.

Crohn's disease (CD) is a sub-type of inflammatory bowel disease (IBD), characterized by intestinal microbial dysbiosis (DeCruz et al., 2015; Lloyd-Price et al., 2019; Pascal et al., 2017; Schaubeck et al., 2016) . Although a complex set of factors including host genetics, lifestyle, diet, and exposure to antibiotics contribute to the development of CD (Ananthakrishnan et al., 2018) , the microbiome-host axis is considered as a fulcrum influencing the manifestation of CD symptoms (Ananthakrishnan et al., 2017; Erickson et al., 2012; Gevers et al., 2014; Halfvarson et al., 2017; Jostins et al., 2012; Pé rez-Brocal et al., 2015) . This is reflected by alterations in the gut microbial compositions of IBD patients at various taxonomic levels, including enhanced abundances of particular species and epithelial barrier responses (Erickson et al., 2012; Lloyd-Price et al., 2019) . Some of these responses result in hyperinflammation (Craven et al., 2012; Shaw et al., 2016) and degradation of the mucus layer which makes the epithelial barrier more vulnerable (Johansson et al., 2014; Roberton and Corfield, 1999; Wu et al., 2011 Wu et al., , 2014 Yu et al., 2014) to microbial components. These observations are also supported by changes in the expression of genes and/or proteins involved in these processes in CD patients exhibiting microbial dysbiosis .

Although several lines of evidence point to the importance of microbiome-host interactions in CD, the underlying molecular mechanisms that mediate the cross-talk between the microbiome and the host are largely unknown. This could be due to (1) the lack of comprehensive datasets with coherent/parallel sampling and high resolution profiling of the gut microbiome and intestinal gene expression from the same samples/subjects; (2) practical difficulties associated with the retrieval of enough biomaterial from luminal samples/biopsies to be able to measure multiple read-outs; (3) the inability of in vitro systems such as organoid cultures, cell-lines, and microfluidic systems to truly reproduce in-vivo conditions and the ensuing microbiome-host interactions; and (4) high costs and practical limitations associated with high-throughput technologies to perform live profiling of microbiome-host interactions.

Computational approaches provide an alternative option to investigate and explore the role of microbiome-host interactions (Guven-Maiorov et al., 2017) . Various methods have been applied in the last decade to infer and predict the interactions between microbial and host proteins (Cui et al., 2016;  Motivated by the availability of synchronous high-resolution datasets profiling the gut microbial alterations and gene expression measurements in CD patients via the IBDMDB/HMP2 project (Lloyd-Price et al., 2019) , we applied a modified workflow (Figure 1 ) of the MicrobioLink pipeline to the IBDMDB/HMP2 dataset to infer the cumulative effects of an altered CD gut microbiome on intestinal gene expression in CD patients. Topological and functional analysis of the resulting networks aided in the identification of the prominent hubs and signaling pathways mediating the effect of the altered microbiome on differentially expressed genes in the ileum and rectum of CD patients. Using an independent cohort composed of CD patients, we validated some of our findings. We also investigated the disease relevance and druggability of the hubs as well as the gene expression signatures modulated by the altered CD microbiome.

By combining information about enhanced abundances and transcriptional activity of bacterial species in CD as reported in the HMP2 IBDMDB dataset with other independent studies reporting the same trend, we identified eight bacterial species (Figure 2A , Table S1 ) relevant to CD. This approach based on multiple lines of evidence (enhanced abundance, transcriptional activity from one same study and additional confirmation of enhanced abundance from another study) imparts greater reliability in making inferences from iScience Article noisy microbiome datasets ( Figure 2A ). The eight identified species include Adhesive invasive Escherichia coli (AIEC), Bacteroides fragilis, Clostridium clostridioforme, Clostridium hathewayi, Clostridium symbiosum, Klebsiella pneumoniae, Ruminococcus gnavus, and Veillonella parvula ( Figure 2A , Table S1 ). By restricting the microbe-host protein-protein interaction prediction (see STAR Methods section for more details) to the homologs (in the eight CD relevant species) of enhanced microbial genes and transcripts (Table S2 ) identified in CD patients in the HMP2 IBDMDB cohort, we inferred 739 orthologous groups with the potential ability to bind to human proteins ( Figure 2B , Table S3 ). Limiting the microbe-host protein-protein interaction prediction to microbial genes and transcripts significantly enhanced in CD proteins further ensures the specificity of the microbial proteins with respect to their role in CD pathogenesis. The 739 orthologous groups thus identified make up about 3.7% of the total number of orthologous groups identified by clustering all the proteins from the eight bacterial species ( Figure 2B ). About 73% (539) of the 739 orthologous groups had protein members from either one or two species only which could indicate that a large proportion of protein-clusters are specific to particular species. In contrast, only ten orthologous groups (1.3%) contained proteins from all the eight species, suggesting that generic effects at the level of microbe-host PPIs are minimal compared to species specific effects ( Figure 2B ). Functional analysis (hypergeometric test; FDR % 0.05) of the 739 orthologous groups with potential binding capabilities identified several metabolic functions in addition to processes involved in post-translational modifications. (E and F) The segregation of the ileal signaling network into sequential layers based on their position in the signaling chains mediating the effect of the bacterial proteins on ileal gene expression. The size of the nodes in E denotes the clustering co-efficient, whereas in F it stands for the betweenness centrality. The top three ranked nodes according to the clustering co-efficient are indicated in E while the top three transcription factors (TFs) ranked by the betweenness centrality are shown in (F).

iScience 25, 103963, May 20, 2022 3 iScience Article Network diffusion highlights the potential modulation of host gene expression in CD by bacterial proteins

Even though we identified bacterial proteins capable of binding to human proteins, it does not necessarily imply that these interactions have an influence on disease development or progression. To fill this gap, we harnessed the power of network diffusion to investigate if and how the inferred bacterial-host PPIs can have downstream effects as measured by gene expression from the host. Using this approach, our analysis revealed some of the potential mechanisms and key regulators that could mediate the effect on host gene expression of the altered bacterial proteins in CD. The signaling network ( Figure 2C ), resembling a scalefree network of 901 nodes and 1600 interactions ( Figure 2D ), captures the potential effect of the altered bacterial proteins on ileal gene expression in CD patients. The prominent components as inferred by network topological metrics such as betweenness centrality (the extent to which a node is a bridge between two other nodes) and clustering co-efficient (the extent to which the neighbors of a node are connected to each other) ( Figure 2F , Table S4 ) point to proteins known to play important functional roles in intestinal inflammation, immune interactions/signaling and often with expression/activity modulations in CD patients (Tables 1 and S4 (Tables 1 and 4 ).

The human signaling network (i.e., excluding the bacterial proteins) capturing the effect of altered bacterial proteins on ileal gene expression in CD patients was over-represented with Reactome pathways such as signaling mediated by interleukins (IL-4 (UniProt: P24394), IL-2 (UniProt: P60568), and IL-13 (UniProt: P35225)), receptor tyrosine-kinases, multiple toll-like receptors including TLR4 (UniProt: O00206), MyD88 (UniProt: Q99836), NFkB (UniProt: P19838) (TRAF6 mediated induction of NFkB upon TLR7/8 or 9 activation), in addition to FCERI mediated MAPK activation (Figures 3A, 3B and Table S5 ). Meanwhile, functional over-representation analysis confined to genes differentially expressed (all upregulated - Figure 3C ) in the ileum concurrently point to interleukin (IL-4,10,13) signaling and immune system functions among other generic functions such as gene expression and transcriptional control ( Figures 3C and S1 ). Using an independent cohort of CD patients, we also validated the microbiome-induced gene expression in the ileum of CD patients in our test cohort. We observed that 78% (92/117) of the DEGs in our test cohort were also differentially expressed in the validation cohort (R2 = 0.907) ( Figure 3D , Table S6 ).

The ileal signaling network with both bacterial and host proteins could be broken down into three distinct modules based on the overall connectivity with other proteins in the network ( Figure 3E ) with common and specific functions ( Figure 3F , Table S7 ).While modules 1 and 2 shared several over-represented pathways such as TRIF-mediated TLR4 signaling, IL-4/IL-13 signaling among many others, several module specific pathways could also be identified. For example, module 1 was specific to the enrichment of interferon and PI3K/AKT signaling while FCERI (UniProt: P12319) and TGF-beta receptor signaling were confined to module 2. Module 3 was, however, peculiarly distinct in terms of its composition by harboring genes corresponding to exclusively enriched pathways such as response to infection by parasites, FCGR-3A mediated phagocytosis, cell-cell communication, BRAF/RAF1 UniProt: P15056) signaling, ERBB2 (UniProt: P04626) signaling among others ( Figure 3F , Table S7 ).

As potential drivers of the ileal gene expression changes described above, several bacterial proteins (which could be clustered into groups based on sequence homology) could be identified to have an influence as shown in Figure 3G . These include proteins with acetyltransferase, chaperone, peptidase, peptidyl-prolyl cis-trans isomerase and proteolytic functions among others. For example, the orthologous group OMA01692, representative of a peptidase S8 domain containing group of proteins present in 4 (C. clostridioforme, C. hathewayi, C. symbiosum, R. gnavus) of the 8 CD relevant bacterial species, potentially influences the expression of 43% of the differentially expressed genes in the ileum of CD patients.

Regulatory control which could be exerted by bacterial influences can potentially induce sitespecific inflammatory responses at the level of gene expression in CD

We also observed that the bacterial species identified as being relevant in CD tend to modulate ileal gene expression through prominent TFs such as SPI1, STAT1 and NFKB1 (UniProt: P19838) as identified by the network topology analysis ( Figure 4A , STAT3 Activation confined to actively inflamed colons in CD patients (Musso et al., 2005) , increased total STAT3 in CD (Mudter et al., 2005) , increased pSTAT3 correlated with Figure 4D ). When all the ileal DE-Gs under the control of the top TFs among the top 20 betweenness centrality ranked nodes were considered together, we observed a consistent trend between the test and validation cohorts (R2 = 0.908).

We also compared the regulatory footprints across the ileal and rectal signaling networks (Figures S4A and S4B) and found several commonalities and differences ( Figure 4E , Table S8 ) at the level of regulation of DE-Gs. Of the 87 DE-Gs regulated by a TF from among the set of top 20 nodes (based on betweenness centrality) of the ileal and rectal signaling networks, 54% (47/87) were common to both the ileum and rectum. In other words, ileum and rectum not only harbor common gene signatures in CD but could also be controlled by a set of shared regulators that could be driving the synchronous gene expression. With an analysis of regulators with at least 10 target DE-Gs in either ileum or rectum, we observed different levels of specificity (ranging from 35 to 61%) not only quantitatively but also qualitatively. For example, the DEG targets (COL1A1, MMP13, MMP2, PLAUR (UniProt: Q03405), CCL2, SOD2 (UniProt: P04179), TFPI2 (UniProt: P48307), TIMP1 (UniProt: P01033)) of NFKB1 exclusive to ileum are associated with extracellular matrix remodeling and organization. Conversely, the targets of NFKB1 common to both ileum and rectum are predominantly associated with immune and inflammatory functions. Furthermore, we also identified pathways and biological processes among the genes in the potential microbiota-mediated signaling networks in a site-specific manner ( Figure 4F , Table S5 ).While processes associated with anti-bacterial defense, cellular fate and cellular adhesion among others were over-represented in the ileal signaling network, the rectal network was characterized by the enrichment of metabolic processes. At the level of Reactome signaling pathways, ileum-specific ones include those related to collagen formation, assembly and degradation as well as platelet adhesion to exposed collagen structures.

Disease relevance and druggability of the potential microbiota-modulated signaling networks

To reveal the relevance of the inferred signaling networks modulated by the bacterial proteins, we checked if the proteins in the networks are associated with CD as reported in the OpenTargets database (Table S9 ) based on at least two sources of evidence (see Table S12 for the full filtered list of CD associated OpenTargets proteins). Among the OpenTargets annotated CD proteins in the ileal and rectal networks, 46/125 (37%) and 40/115 (35%) respectively are derived from evidence sources with a genetic basis although information on host genotypes was not integrated in the current study. In addition, by using drug target information from OpenTarget, some of the proteins in the ileal and rectal signaling networks are known targets of therapeutic molecules used to treat CD (Table 2) . Cumulatively, 24 proteins from the ileal and rectal signaling networks are existing targets of CD drugs. These include prominent cytokines such as IL23B (specific to the rectal network), IL6 and CXCL10 and chemokine receptors (CSF2RA (UniProt: P15509), CSF3R (UniProt: Q99062)) (common to both networks). The therapeutically targeted proteins also cover other functional categories such as those associated with JAK-STAT signaling (JAK1 (UniProt: P23458), JAK2 (UniProt: Q60674), JAK3 (UniProt: P52333), TYK2 (UniProt: P29597)), matrix organization and remodelling (MMP1 (UniProt: P03956), MMP7 (UniProt: P09237), MMP13, ITGA4 (UniProt: P13612)), T-cell proliferation and activation (CD80 (UniProt: P33681)) and enzymes involved in the production of inflammatory precursors (PTGS2).

In addition to the above demonstrated druggability of the potential microbiota modulated signaling network components (irrespective of status of differential gene expression), we also investigated if the Table S10 ) of CD patients potentially by the altered microbiota. Drugs which could induce a reversal of the expression patterns of the DE-Gs in the ileum or rectum include among others gefitinib, canertinib, dasatinib, trimpiramine, enalapril, and wortmannin. As an example, wortmannin influences the expression of genes involved in various phenotypes associated with inflammatory response, susceptibility to bacterial infection, wound healing and macrophage physiology/chemotaxis ( Figure 5C ).

Multiple studies have demonstrated that dysbiosis and alterations of the gut microbiome are associated with phenotypic responses including inflammation of the gut mucosal barrier in CD patients. However, at a meta-level, there is a gap in terms of the molecular mechanisms that mediate the effect of the altered gut microbiome on host responses. In this study, we used a computational pipeline (harnessing the power of microbe-host protein-protein interaction predictions, biological networks and network diffusion concepts) to integrate heterogeneous-omic datasets including whole genome shotgun sequencing of the fecal microbiome and transcriptomic read-outs from ileal and rectal biopsies of CD patients. We identified site-specific signaling networks for ileum and rectum, which may explain how the altered CD gut microbiome could modulate host gene expression. To start with, we identified eight bacterial species (AIEC, B. fragilis, C. clostridioforme, C. hathewayi, C. symbiosum, K. pneumoniae, R. gnavus and V. parvula) associated with the gut microbiome of CD patients compared to non-IBD controls. Although strain level resolution was not possible, species level resolution is expected to provide interesting insights, particularly given that (1) the inferences are based on three different evidence types and (2) only a subset of relevant genes/proteins with increased activities in CD patients are considered for the microbe-host interaction analysis. All of the eight species are either known pathogens or support the growth and survival of pathogens. They are mainly confined to the gut although non-gut niches have also been documented for some of them (Elsayed and Zhang, 2004; Knapp et al., 2017) . Conversely, pathogens such as K. pneumoniaewhich are predominantly found in non-gut niches like the respiratory tract have also been identified as being relevant to CD. This is understandable since K. pneumoniae is known to be a pathogen which has mastered the art of evading host immune defense mechanisms in different immune cells as well as a wide range of organ systems (Bengoechea and Sa Pessoa, 2019). However, despite the fact that depletion of commensals could have an impact on intestinal homeostasis and IBD (Koboziev et al., 2014) , since the goal of the study was to infer the functional effects of the microbiome on the host, we considered, for our global analysis, only microbial species and genes/transcripts with enhanced abundances and activities respectively.

Next, we identified bacterial protein groups involved in modulating the activity of human proteins in CD. The small fraction (3.7%) of such protein groups as a proportion of all possible protein clusters could possibly be attributed to various underlying reasons such as (1) the prevalence of mechanisms along other-omic layers such as metabolism in mediating microbe-host interactions (Metwaly et al., 2020) (Franzosa et al., 2019) (2) the ability of singular bacterial proteins or a small number of them to render noticeable phenotypic changes on the host (Carroll and Maharshak, 2013; Cohavy et al., 2000; Dalwadi et al., 2001) and (3) the lack of structural annotation (especially the lack of domain annotations) given that the interaction prediction between the microbial and host proteins is heavily reliant on such annotations. Nevertheless, based on their functional annotations, we identified several relevant host-interacting bacterial proteins capable of inducing post-translational effects which play a role in the pathogenesis of many chronic diseases including IBD (Caruso et al., Gordon et al., 2020; Schweppe et al., 2015; Zhang et al., 2018 Zhang et al., , 2020a .

The analytical pipeline uncovered the corresponding site-specific (ileum/rectum) signaling networks that captured the effects of the CD relevant bacterial proteins on differentially expressed host genes. Using an independent validation cohort of CD patients carrying the pathogenic B2 enterotype, we could verify our findings as exemplified by overall consensus (78%, R2 = 0.908) between the ileal DE-Gs in the test and validation cohorts. Ranking the proteins using topologically relevant parameters such as betweenness centrality (indicative of the number of shortest paths passing through the protein) (Kincaid and Phillips, 2011) resulted in the identification of central TFs and post-translational signaling proteins. The prevalence of TFs among the BCranked nodes is understandable given that upstream signals modulate TF activity which thus subsequently impacts the expression of their immediate targets downstream (Karin and Smeal, 1992) . Thus, TFs act as signal amplifiers by modulating gene expression and are involved in complex regulatory circuits which are required for quick and adaptive responses in biological systems (Gao et al., 2018) . Some of the top ranked TFs such as STAT3, STAT1, TP53, ESR1, STAT1, SPI1, SMAD3 and MYC are involved in modulating intestinal inflammation and/or known to be playing a role in CD/IBD. STAT3 and STAT1 for example are of interest due to the multiple mechanisms by which STATs impact IBD pathogenesis. STAT1 is differentially expressed in ileal biopsies of CD patients and its deficiency enhances predisposition to gut inflammation (Leon-Cabrera et al., 2018) . Adding to the more complex role of STAT3 in CD is its transient/dynamic involvement independent of genetic aberrations (Ferguson et al., 2010; Willson et al., 2012) .STAT3 is involved not only in mediating inflammatory responses to infections by pathogens (Backert et al., 2014) and generation of inflammatory Th17 cells (Durant et al., 2010; Liu et al., 2008) but also in modulating the return of the epithelium to steady-states following inflammatory insults (Goldsmith et al., 2011) . Synonymous with our observation of STAT3 as a major regulatory hub is the identification as a signaling hub of SOCS3 -a partner in crime with STAT3 in mediating intestinal inflammation. SMAD3, another hub, is involved in the TGF-beta driven inflammatory responses in CD (Monteleone et al., 2001) , mucosal immunity and T-cell activation (Yang et al., 1999) . In addition to the hubs, many of over-represented signaling pathways in the microbiome-modulated ileal and rectal signaling networks were also identified in independent experimental studies. For instance, TLR3 (UniProt: O15455)/TLR4 associated signaling pathways which were over-represented in the ileal network, have been demonstrated to be active in CD patients (Cario and Podolsky, 2000) . TLR4 is also over-expressed in ileal IECs of CD patients (Silva et al., 2008) and upregulated during intestinal inflammation (Hausmann et al., 2002) . This is intriguing since iScience Article TLR pathways, especially those associated with TLR4 are generally activated by bacterial lipopolysaccharides (LPS) in disease states and are also involved in mediating host response to microbial effects (van Bergenhenegouwen et al., 2016; Cerovic et al., 2009; Kløve et al., 2020) . However, our results suggest novel modes of TLR activation mediated by bacterial proteins in the context of CD.

Concomitant with observations from previous studies about the role of NFKB1 in mediating intestinal inflammatory responses in healthy states (Ramakrishnan et al., 2019) and CD (Han et al., 2017; Schreiber et al., 1998; Segain et al., 2000; Sorrentino, 2017) , our results indicated that (1) NFKB1 was among the central regulators in the ileal signaling network and (2) the expression of many inflammatory cytokines was upregulated via NFKB1 by bacterial proteins including those with proteolytic functions as demonstrated by the S8 domain containing peptidase example. Although many stimuli such as butyrate are known to inhibit inflammatory responses via NFKB1 (Segain et al., 2000) , the exact mechanisms have remained elusive, particularly in whether and how bacterial proteins can influence inflammatory responses via NFKB1. Under homeostatic conditions, NFKB1 under the influence of active upstream repressors such as MAP3K7 tends to have reduced activity and thereby results in the suppression of pro-inflammatory molecules. However, bacterial proteins such as the S8 domain containing peptidase present in four of the eight CD relevant species (C. clostridioforme, C. hathewayi, C. symbiosum, R. gnavus) potentially inactivates MAP3K7 by binding to the PCSK cleavage motif CLV_PCSK_PC1ET2_1 and thereby inactivating MAP3K7. This results in downstream effects culminating in the de-repression of TFs including NFKB1which is then able to activate the expression of pro-inflammatory cytokines such as IL6, IFNG, CCL2 among others. Interestingly, the S8 domain containing peptidase in absent in all but one proteome (Faecalibacterium sp. CAG:1138; protein uniprot id -R5FTW4) within Faecalibacterium sp and absent in all of the proteomes in Akkarmensia sp. Thus, bacterial proteins with enhanced abundances and activity in CD and mostly absent in commensals such as Faecalibacterium sp. and Akkarmensia sp which are depleted in CD, can ''switch-off'' mechanisms that under homeostatic conditions would have silenced inflammatory cytokines. That being said, the above described mechanism of action could only be one among many other microbiome-mediated mechanisms modulated by other bacterial proteins and several other players such as small RNAs and metabolites (Valter et al., 2020; Wlodarska et al., 2015) .

By comparing the ileal and rectal signaling networks mediating the possible microbial influences on ileal and rectal gene expression, we identified site specific commonalities and differences in the gene expression as well as the underlying regulatory programs. At the level of genes potentially modulated by the microbiota, 49 and 31% of the rectal and ileal DE-Gs were unique to the corresponding sites, possibly indicating specific effects of the CD altered microbiota on ileum and rectum. This is partially in line with observations from previous studies which report global differences in the transcriptional activity of genes in IBD patients between the two sites (Bogaert et al., 2010; Zhang et al., 2020b) . However, these studies did not consider the plausible microbial influences and hence could be capturing a range of unknown variables on ileal and rectal gene expression. Central regulators such as SPI1, NFKB1 and STAT1 displayed different levels of specificity thus pointing to site-specific regulatory programs mediating the possible effects of the microbiota. Although this is a novel finding with therapeutic implications, it needs further independent validation.

By integrating different levels of information such as collection of genes associated with CD pathogenesis, targets of known CD drugs as well as gene expression signatures in response to drugs, we interpreted the disease relevance and druggability of the ileal and rectal networks. Enrichment of our networks with CD associated genes and drug targets not only assign biological significance to the inferred networks but also impart translational value. Eventhough we did not integrate mutation information from CD patients, only a third of the OpenTargets annotated CD proteins in the ileal and rectal networks are known risk loci for CD. Furthermore, expanding this to the entirety of the two signaling networks, only 6.5% on average are known IBD susceptibility loci, a subset of which include prominent components of the inflam- (Table S14 ). This suggests that extrinsic factors such as potential microbial influences could affect the transient activity of the proteins in a mutation independent manner. Although further studies are required to ascertain this observation, it is concordant with the conclusions of many genetic studies that there is missing heritability (Chen et al., 2014; Jostins et al., 2012) . From a therapeutic viewpoint, functionally and topologically relevant nodes in our networks are known targets of approved IBD drugs such as ll OPEN ACCESS tofacitinib (JAK1, JAK2, JAK3, TYK2) and vedolizumab (ITGA4).ITGA4 interacts with vascular cell adhesion protein 1, a key protein involved in cellular adhesion and communication, and plays a major role in inflammation by recruiting memory and effector T cells to inflamed intestinal tissues (Elices et al., 1990; Reinhardt et al., 1997; Rü egg et al., 1992) . In total, 24 proteins from the ileal and rectal networks could be identified cumulatively as drug targets of 27 different CD therapeutics as inferred from completed or ongoing clinical trials (Tables 2 and S11) . From a translational perspective, though, topological relevant proteins in the signaling networks which are not known CD drug targets could be very handy in designing novel therapeutics since they are important mediators of plausible microbial influences on CD pathogenesis.

Despite the above-mentioned findings, our study is limited by some of the assumptions and existing bottlenecks in terms of data availability as listed below: (1) The lack of cell-type specific read-outs such as gene expression. This is relevant due to the recent studies reporting differences in cell-type compositions and activities between IBD patients and non-IBD controls (Martin et al., 2019; Smillie et al., 2019) (2) non-availability of strain-level resolution of the microbiome in both the test and validation cohorts, especially given that differences in pathogenicity are attributed to strain identities within the same species (Fang et al., 2018; Kalan et al., 2019; Zhang and Zhao, 2016) (3) lack of microbial read-outs from the intestinal mucosal samples (4) interaction prediction between microbial and host proteins are limited to predictions due to lack of large-scale bacterial-host PPIs (5) the effects of other regulatory layers such as miRNAs (Su et al., 2010) or long non-coding RNAs (Chen et al., 2013; Zacharopoulou et al., 2017) have not been considered during the network diffusion process (6) microbial components have been limited to proteins and excluded other molecules such as small RNAs, and metabolites that have been known to modulate host responses (Ahmadi Badi et al., 2020; Huang et al., 2019; Valter et al., 2020; Wlodarska et al., 2015) (7) studies have shown vast differences in composition and abundance between the rectal/colonic and ileal microbiomes (Zoetendal et al., 2008) . Hence, the lack of site-specific read-outs from the ileal and rectal compartments is another drawback of the study.

To conclude, using a combinatorial computational workflow enabling the integration of high-resolution microbial datasets representing the altered gut microbiome with the gene expression from ileal and rectal biopsies of CD patients, we inferred novel mechanisms which could mediate the plausible microbial effects on host responses in CD. By performing customized analyses on the inferred signaling networks representing the signal transduction between the CD altered microbial proteins and gene expression, we identified topologically relevant proteins with functional significance as exemplified by their involvement in CD pathogenesis.

Using specific examples, we pointed out potential mechanisms of action by which microbial proteins with enriched abundances or enhanced transcriptional activity can activate inflammatory responses in CD. Using gene expression data from a paired independent validation cohort of CD patients, we could validate the microbiome-modulated gene expression signatures in the ileum of CD patients. By comparing the inferred ileal and rectal signaling networks, we identified site specific differences in terms of individual genes as well as the regulatory rewiring occurring between the two sites. We also highlight functionally relevant proteins that are already known drug targets or those which could be targeted de novodue to their topological and functional relevance. Taken together along with the validation using an independent cohort of CD patients, our results provide deeper insights into site-specific mechanisms used by the microbiome in CD patients and provide opportunities for development of more tailored therapies.

Detailed methods are provided in the online version of this paper and include the following: 

PS and TA performed the analysis. PS wrote the manuscript. SaV and CC provided data for the validation cohort. SaV performed the differential expression analysis for the validation cohort. MF, JS and BV provided concrete feedbacks and discussions for the manuscript and edited the manuscript. SV supervised the study and edited the manuscript. All authors read and approved the final version of the manuscript.

BV reports financial support for research from Pfizer; lecture fees from Abbvie, Ferring, Takeda Pharmaceuticals, Janssen, and R Biopharm; consultancy fees from Janssen and Sandoz. JS reports lecture fees from Abbvie, Takeda Zhang, X., Ning, Z., Mayne, J., Yang, Y., Deeke, S.A., Walker, K., Farnsworth, C.L., Stokes, M.P., Couture, J.-F., Mack, D., et al. (2020a) . Widespread protein lysine acetylation in gut microbiome and its alterations in patients with Crohn's disease. Nat. Commun. 11, 4120. Zhang, Y., Shen, B., Zhuge, L., and Xie, Y. (2020b) . Identification of differentially expressed genes between the colon and ileum of patients with inflammatory bowel disease by gene co-expression analysis. J. Int. Med. Res. 48, 300060519887268. Zoetendal, E.G., Rajilic-Stojanovic, M., and de Vos, W.M. (2008) . High-throughput diversity and functionality analysis of the gastrointestinal tract microbiota. Gut 57, 1605-1615. iScience Article

The annotated proteomes corresponding to the identified bacterial species were downloaded from Uniprot (The UniProt Consortium, 2018). Where available, the reference proteomes were retrieved. The stand-alone version of the ortholog group mapping software OMA (Altenhoff et al., 2019) was used to identify orthologous groups. Default settings were used to carry out the ortholog analysis. Orthology tests against multiple members of Faecalibacterium sp. and Akkarmensia sp were performed using the NCBI blastp with default parameter settings (E-value % 0.05).

The microbe-host interaction prediction based on the inferred interactions and their downstream effects was carried out along the lines of the general framework laid out in the manuscript by Andrighetti et al. (2020) . To maintain the context specificity (i.e. prevalence in CD), interactions corresponding to bacterial genes/proteins with enhanced abundances in CD (as against the non-IBD controls) (Table S2) were considered for further analysis. Interaction prediction between microbial and host proteins were performed by using the domain-domain and domain-motif interaction prediction principle applied in the manuscript by Andrighetti et al. (2020) and used elsewhere in previous studies (Sudhakar et al., 2019) . Briefly, the interaction prediction assumes that protein-pairs with known interacting domain-motif or domain-domain pairs also interact. With the domain-motif method, in order to capture the desired direction of signal transduction (i.e., the modulation of host proteins by the bacterial proteins), we imposed the condition that the interacting domains need to be harbored by the microbial proteins and the motifs on the host proteins. Pfam domain annotation for the microbial proteins was retrieved from Uniprot (The UniProt Consortium, 2018) whereas information about the occurrence of motifs and the gold standard list of domain-motif interactions were obtained from the Eukaryotic Linear Motif database (Dinkel et al., 2014; Gouw et al., 2018) . The resulting set of predicted interactions between microbial and host proteins were then subjected to a disordered region based structural filtering procedure as described by Andrighetti et al. (2020) . Disordered region predictions were performed by using IUpred2A (Mé szá ros et al., 2018) . In short, only protein-protein interactions corresponding to motifs that lie within disordered regions of the host proteins were considered.

Where available, the reference proteomes corresponding to the bacterial species were downloaded from Uniprot (The UniProt Consortium, 2018).Since we did not find evidence for the majority of inferred bacterial species to be intracellular (with the exception of Adherent-Invasive E. coli(AIEC) to host cells, we considered them to be extracellular. Host proteins were also selected based on their localization information. iScience Article For compiling the extracellular list, selection was done by retaining only proteins located in ''plasma membrane'' and ''extracellular region'' according to the databases LocDB (Rastogi and Rost, 2011) and Human Protein Atlas (HPA) (Thul and Lindskog, 2018) . Proteins corresponding to all the remaining localization terms were considered as intracellular.

The set of differentially expressed genes (DE-Gs) between CD and non-IBD subjects were retrieved from Lloyd-Price et al. (2019) . Since biopsies were taken from two different sites (rectum and ileum), we obtained independent lists of differentially expressed genes for each of the two sites.

To infer the most probable signaling chains which capture the functional effects of the bacterial proteins on host processes, we used the network diffusion algorithm TieDie (Paull et al., 2013) . TieDie not only discounts for topological biases in networks such as hubs but also mines out logically consistent paths between a given set of start and stop nodes. Host proteins which are predicted to be bound by the bacterial proteins were used as the start nodes. As for the stop nodes, the transcription factors (TFs) of the differentially expressed genes in CD patients were used. TFs of the DE-Gs were retrieved from the set of directed and signed transcriptional regulatory interactions from DoRothEA (Garcia-Alonso et al., 2019) . For each of the TFs in the set of stop nodes, weights were calculated using the following formula (Equation 1). The equation considers both the fold change of the DEG as well as the sign (i.e stimulatory or inhibitory) of the relationship between the TF and the DEG. The set of directed and signed protein-protein interactions from OmniPath(Tü rei et al., 2016) were used as a source of interactions for inferring the signaling chains between the start and the stop nodes.

Given that there exist upper limits for the diameter of real-world biological networks, we set out path length limitations in order to capture the most relevant signaling chains inferring from the network diffusion methodology described above. In the case of signaling chains mediated by the set of extracellular bacterial proteins, we set a path length of four (i.e., the bacterial protein followed by the host receptor, an intermediary host protein, the TF and the DEG). As for the set of intracellular bacterial proteins (all from AIEC), we set a path length of three (i.e., the bacterial protein followed by the host receptor, the TF and the DEG). Paths of the desired lengths (as described above) were extracted using the command line version of PathLinker (Ritz et al., 2016) .

Network diffusion and chain extractions were performed separately for every bacterial species and the two sets of DE-Gs corresponding to both the ileal and rectal biopsies. Only chains not corresponding to interactions between bacterial and host proteins mediated by structural complex formation events were prioritized in order to retain the chains specific to post-translational modifications elicited by the bacterial proteins on the host proteins. Furthermore, chains corresponding to post-translational modifications elicited by bacterial histidine kinases (HKs) were discarded since the phosphorylation domains of HKs are predominantly embedded in the cytoplasmic-membrane boundary of bacterial cells (thus rendering the HKs incapable of inducing phosphorylation events on host proteins).

The R packages ReactomePA (Yu and He, 2016) and clusterProfiler (Yu et al., 2012) were used to perform the enrichment analysis. An adjusted P-value cut-off of % 0.1 was used to infer the statistically significant enrichment terms. For the enrichment analysis of the signaling networks, all the nodes in the directed PPI network were set as the background. For the corresponding analysis of the DE-Gs in the ileal signaling network, the entire set of DE-Gs were used as the background. For the identification of site-specific enrichment terms, the enrichment was performed with different gene sets separated by site and gene categories. All human genes (receptors, intermediary proteins, transcription factors (TFs) and DE-Gs) were considered for the site-specific enrichment analysis.

iScience 25, 103963, May 20, 2022 23 iScience Article Validation cohort A cohort of paired 16S read-outs from stool samples (n = 20) and gene expression data (n = 20) from ileal biopsies in CD patients, cross-sectionally collected at the University Hospitals Leuven was available for independent validation. The ethics committee of the University Hospitals Leuven approved the study (IRB approvals, B322201213950/S53684 and B322201627472/S57662). All individuals gave written informed consent.

Cell counting for all samples, was performed in duplicate. Briefly, 0.2 g frozen (À80 C) aliquots were dissolved in physiological solution to a total volume of 100 mL (8.5 g/L NaCl; VWR International, Germany). Subsequently, the slurry was diluted 1,000 times. Samples were filtered using a sterile syringe filter (pore size of 5 mmm; Sartorius Stedim Biotech GmbH, Germany). Next, 1 mL of the microbial cell suspension obtained was stained with 1 mL SYBR Green I (1:100 dilution in DMSO; shaded 15 min incubation at 37 C; 10,000 concentrate, Thermo Fisher Scientific, Massachusetts, USA). The flow cytometry analysis was performed using a C6 Accuri flow cytometer (BD Biosciences, New Jersey, USA) based on Prest et al. (Doherty et al., 2018) Fluorescence events were monitored using the FL1 533/30 nm and FL3 >670 nm optical detectors. In addition, also forward and sideward-scattered light was collected. The BD Accuri CFlow software was used to gate and separate the microbial fluorescence events on the FL1/FL3 density plot from the faecal sample background. A threshold value of 2000 was applied on the FL1 channel. The gated fluorescence events were evaluated on the forward/sideward density plot, as to exclude remaining background events. Instrument and gating settings were kept identical for all samples (fixed staining/gating strategy (Doherty et al., 2018) ). Based on the exact weight of the aliquots analysed, cell counts were converted to microbial loads per gram of faecal material. Faecal Moisture content was determined as the percentage of mass loss after lyophilization from 0.2 g frozen aliquots of non-homogenized faecal material (À80 C).

Faecal DNA extraction and microbiota profiling was performed as described previously (Gevers et al., 2014) . Briefly, DNA was extracted from faecal material using the MoBio PowerMicrobiome RNA isolation kit with the addition of 10 min incubation at 90 C after the initial vortex step. The V4 region of the 16S rRNA gene was amplified with primer pair 515F/806R (Caporaso et al., 2011) . Sequencing was performed on the Illumina MiSeq platform (San Diego, California, USA), to generate paired-end reads of 250 bases in length in each direction. After de-multiplexing, fastq sequences were merged using FLASH (Mago c and Salzberg, 2011) software with default parameters. Successfully combined reads were filtered based on quality (>90% of nucleotides with quality score of 30 or higher for every read) using Fastx tool kit (http:// hannonlab.cshl.edu/fastx_toolkit/). Chimeras were removed with UCHIME (Edgar et al., 2011) . Faecal samples were processed altering the protocol above to dual-index barcoding. Pre-processing was performed using the DADA2 (Callahan et al., 2016) pipeline v1.6.0.

For the relative microbiome matrix, each sample was downsized to 10,000 reads by random selection of reads. Samples with less than 10,000 reads were excluded (two samples). The taxonomy of reads was assigned using RDP classifier 2.12 (Wang et al., 2007) .

The quantitative microbiome profiling matrix was built as described by Vandeputte and colleagues (Vandeputte et al., 2017) . In short, samples were downsized to even sampling depth, defined as the ratio between sampling size (16S rRNA gene copy number corrected sequencing depth) and microbial load (average total cell count per gram of frozen faecal material). 16S rRNA gene copy numbers were retrieved from the ribosomal RNA operon copy numberdatabaserrnDB (Stoddard et al., 2015) . The copy number corrected sequencing depth of each sample was rarefied to the level necessary to equate the minimum observed sampling depth in the cohort. Samples with resulting rarefied read counts <150 were excluded from QMP analyses. Rarefied genus abundances were converted into numbers of cells per gram.

Characterization of the microbiota profiles below genus-level was performed using the DADA2 (Naftali et al., 2016) pipeline sequence variants, with taxonomy assignment by RDP classifier v 2.12 (Wang et al., 2007) . A phylogenetic tree for the sequence variants was reconstructed using the dada2 sequence variant ll OPEN ACCESS

The RAVEN toolbox and its use for generating a genome-scale metabolic model for Penicillium chrysogenum

Particle-counting immunoassay (PACIA) of pregnancy-specific beta 1-glycoprotein, a possible marker of various malignancies and Crohn's ileitis

Genetic factors in chronic inflammation: single nucleotide polymorphisms in the STAT-JAK pathway, susceptibility to DNA damage and Crohn's disease in a New Zealand population

Benchmark and integration of resources for the estimation of human transcription factor activities

The treatment-naive microbiome in new-onset Crohn's disease

Mu opioid signaling protects against acute murine intestinal injury in a manner involving Stat3 signaling

Impaired estrogen signaling underlies regulatory T cell loss-of-function in the chronically inflamed intestine

Comparative host-coronavirus protein interaction networks reveal pan-viral disease mechanisms

New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0

Structural host-microbiota interaction networks

Dynamics of the human gut microbiome in inflammatory bowel disease

NF-kappa B activation correlates with disease phenotype in Crohn's disease

Distinct histopathologic and molecular alterations in inflammatory bowel disease-associated intestinal adenocarcinoma: c-MYC amplification is common and associated with mucinous/signet ring cell differentiation

Toll-like receptors 2 and 4 are up-regulated during intestinal inflammation

Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0

Dirichlet multinomial mixtures: generative models for microbial metagenomics

R Package ''Coin (Comprehensive R Archive Network (CRAN))

Bacteriaautophagy interplay: a battle for survival

Small RNAs -big players in plantmicrobe interactions

Bacteria penetrate the normally impenetrable inner colon mucus layer in both murine colitis models and patients with ulcerative colitis

Host-microbe interactions have shaped the genetic architecture of inflammatory bowel disease

Network topology measures

Role of the enteric microbiota in intestinal homeostasis and inflammation. Free Radic

HPIDB-a unified resource for host-pathogen interactions

Immunohistochemical expression of CDX2, b-catenin, and TP53 in inflammatory bowel disease-associated colorectal cancer

NetCooperate: a networkbased tool for inferring host-microbe and microbe-microbe cooperation

Myeloid-derived cullin 3 promotes STAT3 phosphorylation by inhibiting OGT expression and protects against intestinal inflammation

ViRBase: a resource for virus-host ncRNAassociated interactions

Loss of STAT3 in CD4+ T cells prevents development of experimental autoimmune diseases

Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

FLASH: fast length adjustment of short reads to improve genome assemblies

Singlecell analysis of crohn's disease lesions identifies a pathogenic cellular module associated with resistance to anti-TNF therapy

IUPred2A: context-dependent prediction of protein disorder as a function of redox state and protein binding

Integrated microbiota and metabolite profiles link Crohn's disease to sulfur metabolism

Blocking Smad7 restores TGF-beta1 signaling in chronic inflammatory bowel disease

DirichletMultinomial: Dirichlet-Multinomial Mixture Model Machine Learning for Microbiome Data

Activation pattern of signal transducers and activators of transcription (STAT) factors in inflammatory bowel diseases

Signal transducers and activators of transcription 3 signaling pathway: an essential mediator of inflammatory bowel disease and other forms of intestinal inflammation

Distinct microbiotas are associated with ileum-restricted and colon-involving crohn's disease

Computational approaches for prediction of pathogen-host protein-protein interactions

CRAN -Package Vegan (Comprehensive R Archive Network (CRAN))

A microbial signature for Crohn's disease

Discovering causal pathways linking genomic events to transcriptional states using Tied Diffusion through Interacting Events (TieDIE)

Metagenomic analysis of crohn's disease patients identifies changes in the virome and microbiome related to disease status and therapy, and detects potential interactions and biomarkers

Intestinal non-canonical NFkB signaling shapes the local and systemic immune response

LocDB: experimental annotations of localization for Homo sapiens and Arabidopsis thaliana

Neutrophils can adhere via alpha4beta1-integrin under flow conditions

Pathways on demand: automated reconstruction of human signaling networks

Host-microbe protein interactions during bacterial infection

Butyrate inhibits inflammatory responses through NFkappaB inhibition: implications for Crohn's disease

Noncanonical NF-kB

Combinatorial regulation of transcription factors and microRNAs

Targeted interplay between bacterial pathogens and host autophagy

CIS3/ SOCS3/SSI3 plays a negative regulatory role in STAT3 activation and intestinal inflammation

Double-stranded RNA induces cathelicidin expression in the intestinal epithelial cells through phosphatidylinositol 3-kinaseprotein kinase Cz-Sp1 pathway and ameliorates shigellosis in mice

A gp130-Src-YAP module links inflammation to epithelial regeneration

UniProt: the universal protein knowledgebase

The human protein atlas: a spatial map of the human proteome

OmniPath: guidelines and gateway for literature-curated signaling pathway resources

Extracellular vesicles in inflammatory bowel disease: small particles, big players

Quantitative microbiome profiling links gut community variation to microbial load

Expression levels of 4 genes in colon tissue might be used to predict which patients will enter endoscopic remission after vedolizumab therapy for inflammatory bowel diseases

Proteopathogen2, a database and web tool to store and display proteomics identification results in the mzIdentML standard

IBD-associated dysplastic lesions show more chromosomal instability thansporadic adenomas

Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy

STAT3 genotypic variation and cellular STAT3 activation and colon leukocyte recruitment in pediatric Crohn disease

An integrative view of microbiome-host interactions in inflammatory bowel diseases

Genome-wide gene expression differences in Crohn's disease and ulcerative colitis from endoscopic pinch biopsies: insights into distinctive pathogenesis

Epithelial inducible nitric oxide synthase causes bacterial translocation by impairment of enterocytic tight junctions via intracellular signals of Rhoassociated kinase and protein kinase C zeta

Commensal bacterial endocytosis in epithelial cells is dependent on myosin light chain kinase-activated brush border fanning by interferon-g

Targeted disruption of SMAD3 results in impaired mucosal immunity and diminished T cell responsiveness to TGF-beta

Commensal bacterial internalization by epithelial cells: an alternative portal for gut leakiness

ReactomePA: an R/ Bioconductor package for reactome pathway analysis and visualization

clusterProfiler: an R package for comparing biological themes among gene clusters

Enteric dysbiosis promotes antibiotic-resistant bacterial infection: systemic dissemination of resistant and commensal bacteria through epithelial transcytosis

The contribution of long non-coding RNAs in Inflammatory Bowel Diseases

Strain-level dissection of the contribution of the gut microbiome to human metabolic disease

Metaproteomics reveals associations between microbiome and intestinal

bacterial species with alterations between CD (n = 67) and non-IBD subjects (n = 27) at both the levels of abundance and transcriptional activity, species with enhanced abundances displaying a t-statistic of R 2 were selected. From the same study, bacterial species contributing to alterations in the community composition as measured by the dysbiotic index between CD and non-IBD subjects (at both the levels of abundance and transcriptional activity) displaying enhanced abundances displaying a t-statistic of R 2 and Q-value % 0.05 were selected. The bacterial species identified above were subjected to a literature search to check if they have been observed in at least one other study to have increased abundance levels in CD patients. This was done to improve the reliability of our inferences about the bacterial species inferred as being important in the context of CD (Table S1). Genes with enhanced abundances/transcriptional activity in CD were determined from the metagenomic and metatranscriptomic profiles

De-novoprediction of K-numbers for the chosen proteomes were performed using BlastKOALA

All statistical tests used were two-sided. All p values were corrected for multiple testing using the Benjamini-Hochberg method (reported as FDR) when multiple tests were performed on lists (n > 1) of features (e.g. taxa-metadata or taxa-taxa associations), also when performing multiple pairwise group (n > 2) comparisons (e.g. Kruskal-Wallis test with post-hoc Dunn test)

Microbiome and physiological features associations

Taxa unclassified at genus level or present in less than 20% of samples were excluded from the statistical analyses. Pearson or spearman correlations were used respectively for linear or rank-order correlations between continuous variables. Mann-Whitney U was used to test median differences of continuous variables between two different groups. For more than two groups, Kruskal-Wallis test with post-hoc Dunn test were used

Counts were then normalized for library size using the DESeq2 package (Love et al., 2014), and only protein-coding genes (Ensemble hg 19 reference build) having an average of at least 10 normalised read counts were considered for further analysis (n = 14,263)

Only those proteins (# = 720) with at least two different sources of evidence were considered for further analysis. To identify the possible drugs or perturbations which can mimic or reverse the signatures of the DE-Gs (in the CD patients) in the microbiota-induced ileal and rectal networks, the DEG signatures were matched against the Library of Integrated Network-based Cellular Signatures (LINCS) L1000 dataset

This study did not generate new unique reagents.Data and code availability d The raw datasets used for the study have been deposited at ArrayExpress (accession number: E-MTAB-10395). d Additional information supporting the findings will be made available upon request from the lead contact Dr. Padhmanand Sudhakar (padhmanand.sudhakar@kuleuven.be).

For the validation cohort, patients were selected based on (a) their identification as having active CD as assessed by the endoscopic scores and (b) having paired fecal microbiome profiling datasets. Non-IBD controls were similarly chosen (paired ileal gene expression and fecal microbiome profiling) from the Flemish Gut Flora Project. The ethics committee of the University Hospitals Leuven approved the study (IRB approvals, B322201213950/S53684 and B322201627472/S57662). All individuals gave written informed consent.