We searched the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database for Crohn’s disease data with the following inclusion criteria: 1) key words“(Crohn OR CD OR IBD OR inflammatory bowel disease) AND (anti-TNF OR infliximab OR adalimumab OR certolizumab OR golimumab)”; 2) Homo sapiens; 3) expression profiling by array OR high throughput sequencing; 4) submitted date <01/10/2022; and 5) datasets or series. After that, we reviewed every data under the following exclusive criteria: 1) recruited Crohn patients<6; 2) no Crohn’s anti-TNFα therapy record; 3) no anti-TNF response record; 4) responder or nonresponders <3; and 5) no clear endoscopy evaluations. For discovery cohorts, we tried to choose homogeneous data. For example, all the anti-TNFα therapy response outcomes were based on the endoscopy measure after therapy. Since adequate discovery data of candidates were lacking, we chose mixed Crohn’s disease majority’s data without a clear subtype for validation cohort’s selection in a slightly relaxed range.

We performed differential expression analyses between responders and nonresponders to infliximab therapy using the limma R package (Team, 2013; Ritchie et al., 2015). Due to the limited recruited Crohn’s disease cohorts with the anti-TNFα response, for obtaining sufficient differential expressed genes, all significant differentially expressed genes (p < 0.05) were further analyzed using Gene Ontology (GO) enrichment in the Metascape database (https://metascape.org) (Zhou et al., 2019). The common differential genes in discovery cohorts were obtained by a Venn package in the UpsetR (Conway et al., 2017). Furthermore, the interactions among these genes were obtained from the STRING database (https://string-db.org/) (Szklarczyk et al., 2021). To identify the core network, we calculated the top differentially expressed genes using the degree algorithm in the CytoHubba package of Cytoscape software (version 3.8.2) (Shannon et al., 2003; Chin et al., 2014).

We further performed prediction analysis of microarrays (PAM), a type of classification based on nearest centroids. The PAM R package provides an accurate predictor that may outperform much more complicated methods (Dabney, 2005; Korkola et al., 2009; Arijs et al., 2010). Owing to the relatively large group of nonresponsive patients in GSE16879, we selected this dataset for PAM. Genes with a minimum classification error were further analyzed and combined with the top differentially expressed genes. Interactions of the aforementioned genes were further evaluated using the degree algorithm of the CytoHubba package. Immune cell scoring in the GSE16879 dataset was calculated based on the xCell website using a curve fitting approach for linear comparison of cell types and a novel spillover compensation technique for separating them (Aran et al., 2017). Different immune cell type profiles were displayed.

The top differentially expressed genes were further assessed in the discovery cohort, and the top five differentially expressed genes were further validated in the discovery cohort of responders and nonresponders with IBDs (validation cohort). Receiver operating characteristic (ROC) curve analysis was then performed to evaluate the predictive values of these top genes in response to anti-TNFα therapy in patients with CD.

We explored Crohn’s disease data entitled “PREDICT 2021 paper: CD” on the single-cell portal (https://singlecell.broadinstitute.org/single_cell) held by the Broad Institute of MIT and Harvard. The Crohn’s disease single-cell data include 27 volunteers’ 201,883 single-cell transcriptomes. Then we obtained the cell types’ tSNE (t-Distributed Stochastic Neighbor Embedding) map and TLR2 expression tSNE map on the interactive visualization web tools.

We then performed coculture of immune cells (THP1 cells, a human monocytic cell line) and colonic epithelial cells (Caco2, a human colonic epithelial cell line) to validate the bioinformatics results. Caco2 and THP1 cells were obtained from the Cancer Research Institute of Central South University, China. THP1 and Caco2 cells were cultured in RPMI-1640 or MEN medium (Gibco; Thermo Fisher United States) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher United States) at 37°C in an atmosphere containing 5% CO₂. THP1 cells were stimulated using different concentrations of lipopolysaccharide (LPS; Sigma-Aldrich, Merck KGaA) or infliximab (Remicade; Cilag AG, Sweden), and Cell Counting Kit-8 (CCK8) assays (Dojindo, Japan) were performed to determine the appropriate concentration to use in subsequent experiments. THP1 cells were then transfected with a Toll-like receptor (TLR) 2 overexpression vector or TLR2 small interfering RNA (siRNA) using Lipo 2000 (Invitrogen, Carlsbad, CA, United States). Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and Western blotting were further performed to validate the overexpression and knockdown of TLR2 in THP1 cells (Proteintech, United States). The TLR2 inhibitor C29 (TargetMol, United States) was added to block TLR2 overexpression in the coculture system. Caco2 cells were treated with 100 ng/ml LPS for 24 h and then with 100 μg/ml infliximab for another 24 h in the lower chambers of 6-well Transwell plates. The cells were then cocultured with TLR2-overexpressing THP1 cells cultured with 50 μM C29. Total RNA from the Caco2 cells was obtained after 24 h of coculture, and epithelial inflammation and tight junctions (representing epithelial permeability) were assessed using qPCR. RT-qPCR and Western blotting were carried out as described previously (Luo et al., 2019). All primers and short hairpin sequences are listed in the supplemental data.

Statistical analysis of bioinformatics data was performed using the R package. One-way analysis of variance (ANOVA) was used to compare multiple randomized sets of data. Unpaired t-tests were used to compare double randomized data. Results with p values less than 0.05 were considered significant.

A total of 30 data were included in the selection of further analysis. After carefully checking every data by the exclusion criteria, 25 data were excluded. Due to the limited eligible data for analysis, we choose total Crohn’s disease data (GSE111761, GSE107865, GSE52746, GSE16879) as discovery cohorts, which will help in revealing a reliable differential network between responders and nonresponders. For the validation cohort, we choose another mixed IBD data (GSE159034) dominated by Crohn’s disease patients without clear subtype information. The detailed information on included data is listed in Table 1. We gave priority to samples with expression profiles before the initial infliximab therapy in patients with CD (i.e., GSE16879). If there was no sampling information, we analyzed the results. First, our discovery analysis contains 47 responders and 30 nonresponders among patients with CD receiving anti-TNFα therapy. Second, differentially expressed gene networks and pathways were obtained under the p-value <0.05 for enough differential expressed genes. We performed PAM analysis and gene set enrichment analysis (GSEA) to validate the core predictor. Validation and ROC curve analyses were further performed in an independent Crohn’s disease cohort and a relatively large cohort. Coculture experimental validation was conducted. The study analysis flow is depicted in Figure 1. The four data series were downloaded from GEO and were normalized using the limma R package (Ritchie et al., 2015) (Supplementary Figure S1A).

From the analysis of differentially expressed genes between responders and nonresponders, due to the insufficient significant differential genes under the adj p-value, we identified 4690, 3422, 1020, and 2811 differentially expressed genes in GSE16879, GSE52746, GSE107865, and GSE111761 datasets, respectively, under p < 0.05 (Figure 2A, Supplementary Table S1). Some of the datasets shared common differentially expressed genes (purple) and interactions among genes (blue), as shown in Figure 2B. Overlaps among the four differentially expressed genes revealed that 32 common differentially expressed genes exhibited shared characteristics between responders and nonresponders of different cohorts (Figure 2C) by a Venn package in the UpsetR (Conway et al., 2017). C-X-C motif chemokine receptor (CXCR) 1, CXCR2, TLR2, and triggering receptor expressed on myeloid cells (TREM1) were common differentially expressed genes. Furthermore, GO enrichment (Zhou et al., 2019) of all differentially expressed genes showed that innate immunity was involved in the response of patients with CD to the anti-TNFα therapy. Specifically, lymphocyte activation, T-cell activation, leukocyte migration, and cell adhesion were the top GO enrichment pathways (Figure 2D, Supplementary Table S1). STRING is a database of known and predicted protein–protein interactions (Szklarczyk et al., 2021), including direct (physical) and indirect (functional) associations based on computational prediction, knowledge transfer between organisms, and interactions aggregated from other (primary) databases. Thus, interactions among common differentially expressed genes were obtained from the STRING database (Figure 2E). The degree algorithm based on this interaction network (Chin et al., 2014) revealed the core components between unresponsive and responsive subgroups. TLR2, TREM1, CXCR1, formyl peptide receptor 1 (FPR1), and FPR2 were the top five differential genes.

To further evaluate the predictive value of the newly built model, PAM was performed based on the nearest centroid classification to classify responders and nonresponders from the GSE16879 dataset. In the PAM model, the threshold corresponding to the lowest misclassification rate was 3.938, and 10 risk genes corresponding to this threshold were considered (Figures 3A,B, Supplementary Table S2). The larger range between the no-score (nonresponders) and yes-score (responders) indicated a better value of classification (Figure 3C, Supplementary Table S2). The good classification efficiency supported the roles of these 10 genes, which included matrix metalloproteinase (MMP) 1, MMP3, regulator of G protein signaling 2, and alpha-2-macroglobulin, in our predictive model (Figure 3D, Supplementary Table S2). A combination model of the top differentially expressed genes and lowest misclassification genes was built using interaction data from the STRING database (Supplementary Figure S1C). A set of core genes, including TLR2, tissue inhibitor of metalloproteinases 1, CXCR1, and TREM1, emerged after the degree calculation based on the combination model in Cytoscape (Figure 3F).

GSEA is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states (e.g., phenotypes) (Subramanian et al., 2005). Unlike differential analysis, which focuses on individual gene differences, GSEA derives its power by focusing on gene sets, that is, groups of genes that share common biological functions, chromosomal locations, or regulation (Subramanian et al., 2005). It includes three steps, namely, calculation of an enrichment score (ES), estimation of the significance level of the ES, and adjustment for multiple hypothesis testing (Subramanian et al., 2005). Thus, we next performed GSEA to detect the core genes between responders and nonresponders. When we treated the response to infliximab as a phenotype in the GSE16879 dataset, a Th1 immune response gene showed an enrichment score of 0.54 (p = 0.015), representing one of the top three enrichment gene sets (Figure 3E). Importantly, we confirmed that TLR2 was one of the genes in the core enrichment gene list for this Th1 gene set (Supplementary Figure S1B, Supplementary Table S2). However, no other previous core genes were found in the top three core enrichment genes.

Immune cell scoring of the GSE16879 dataset was then performed by uploading the expression data to the xCell website (Aran et al., 2017). A heatmap of immune cell scoring demonstrated divergent immune landscapes between responders and nonresponders (Figure 4A, Supplementary Table S3). The nonresponders had higher immune scores than responders and controls. Scores for macrophages, activated dendritic cells, natural killer cells, and neutrophils in all samples showed significant differences between responders and nonresponders (Figure 4B, Supplementary Table S3; one-way ANOVA, p < 0.05 or 0.001). Thus, we hypothesized that TLR2 in immune cells may have biological effects on the colonic response to infliximab. Moreover, after tSNE reanalysis of TLR2 expression and cell types in public single-cell sequencing data of Crohn’s disease patients. Data from the Single Cell Portal support the major contributor of TLR2 expression in colonic tissues is the mononuclear phagocyte system, and lead to our THP-1 selection in the coculture validation (Figure 4C). Based on these findings, we chose TLR2 as an experimental biomarker to validate the bioinformatics findings.

To validate and evaluate these newly built predictors’ value, we first displayed the basic expression profiles of top differential genes using a scatter plot for GSE16879 (Figure 5A). The expression signature in the discovery cohort was also validated in independent IBD cohort (GSE159034) including 9 Crohn’s disease patients. Importantly, we confirmed the significant differences in TLR2, and consistent difference in TREM1, CXCR1, FPR1, and FPR2 between responder and nonresponder groups (Figure 5B). ROC analysis using the GSE16879 dataset also demonstrated the significant predictive value of the top five differentially expressed genes [Figure 5C; TLR2, area under the curve (AUC): 0.829, p = 0.001, 95% confidence interval (CI): 0.680–0.979; TREM1. AUC: 0.844, p < 0.001, 95% CI: 0.716–0.873; CXCR1, AUC: 0.841, p < 0.001, 95% CI: 0.708–0.974; FPR1, AUC: 0.894, p < 0.001, 95% CI: 0.778–1.0; FPR2, AUC: 0.824, p < 0.001, 95% CI: 0.678–0.969].

A coculture system of THP1 and Caco2 cells was established to elucidate the effects of TLR2 expression on the colonic response to anti-TNFα therapy (Figure 6A). Efficient TLR2 knockdown was achieved in THP1 cells using siRNA (Supplementary Figure S1D), and we evaluated the effects of TLR2 overexpression and knockdown on TLR2 expression in THP1 cells using Western blotting (Figure 6B). CCK8 assays revealed the appropriate LPS and infliximab concentrations to use in coculture (Figure 6C). Importantly, we showed that infliximab could alleviates the amplified inflammation [i.e., IL-1β, IL-6, and monocyte chemotactic protein 1 (MCP1)] and improves the reduced epithelial permeability [i.e., occludin and zona occludens (ZO)-1] in the LPS-induced model. Furthermore, the overexpression of TLR2 in THP1 cells amplified colonic epithelial inflammation, as measured by IL-1β, IL-6, and MCP1 (Figure 5C; one-way ANOVA, p < 0.05 or 0.001), and reduced the expression of tight junction proteins, such as occludin and ZO-1, even in the context of infliximab rescue. The TLR2 inhibitor C29 alleviated TLR2 overexpression-induced amplification of epithelial inflammation and impaired permeability during the infliximab therapy. By contrast, TLR2 knockdown showed an effective response to infliximab rescue after the LPS induced inflammation. (Figure 6C; one-way ANOVA, p < 0.05 or 0.001).