Hyperactive neutrophil chemotaxis contributes to anti-tumor necrosis factor-α treatment resistance in inflammatory bowel disease.

BACKGROUND AND AIM: Anti-tumor necrosis factor-α (anti-TNF-α) agents have been used for inflammatory bowel disease; however, it has up to 30% nonresponse rate. Identifying molecular pathways and finding reliable diagnostic biomarkers for patient response to anti-TNF-α treatment are needed.
METHODS: Publicly available transcriptomic data from inflammatory bowel disease patients receiving anti-TNF-α therapy were systemically collected and integrated. In silico flow cytometry approaches and Metascape were applied to evaluate immune cell populations and to perform gene enrichment analysis, respectively. Genes identified within enrichment pathways validated in neutrophils were tracked in an anti-TNF-α-treated animal model (with lipopolysaccharide-induced inflammation). The receiver operating characteristic curve was applied to all genes to identify the best prediction biomarkers.
RESULTS: A total of 449 samples were retrieved from control, baseline, and after primary anti-TNF-α therapy or placebo. No statistically significant differences were observed between anti-TNF-α treatment responders and nonresponders at baseline in immune microenvironment scores. Neutrophil, endothelial cell, and B-cell populations were higher in baseline nonresponders, and chemotaxis pathways may contribute to the treatment resistance. Genes related to chemotaxis pathways were significantly upregulated in lipopolysaccharide-induced neutrophils, but no statistically significant changes were observed in neutrophils treated with anti-TNF-α. Interleukin 13 receptor subunit alpha 2 (IL13RA2) is the best predictor (receiver operating characteristic curve: 80.7%, 95% confidence interval: 73.8-87.5%), with a sensitivity of 68.13% and specificity of 84.93%, and significantly higher in nonresponders compared with responders (P < 0.0001).
CONCLUSIONS: Hyperactive neutrophil chemotaxis influences responses to anti-TNF-α treatment, and IL13RA2 is a potential biomarker to predict anti-TNF-α treatment response.

Introduction

Inflammatory bowel disease (IBD) is an idiopathic and relapsing–remitting chronic inflammatory disorder characterized by a susceptible genetic background, causing immunological dysfunction and intestinal microbiome dysbiosis. The long‐standing mucosa inflammation destructs tight junctions, induces intestinal barrier injury and permeability, and increases the incidence of colonic neoplasia. It is estimated that the prevalence of IBD exceeds 0.3% in North America, Oceania, and many countries in Europe. With the incidence rising in the newly industrialized countries, including Brazil and Taiwan, thus, IBD places a large burden on public health services and healthcare economies.

Tumor necrosis factor‐α (TNF‐α) is a pleiotropic cytokine that participates in several pathological processes in IBD and is recognized as a pro‐inflammatory cytokine. A soluble, biologically active homotrimer TNF‐α originally from the monomeric TNF‐α claved by TNF‐α‐converting enzyme (TACE) via proteolysis. TNF‐α activity is mediated through binding to the TNF receptors I and II (TNFRI and TNFRII). This binding activates immune cells response and pro‐inflammatory cytokine and chemokine productions, such as IL‐1, IL‐6, IL‐8, and RANTES. It also increases the expression of adhesion molecules, production of matrix metalloproteinase, and induction of apoptosis. The use of anti‐TNF‐α compounds such as full monoclonal IgG1 antibodies (infliximab and adalimumab), pegylated anti‐TNF‐α F[ab′]2 fragment (certolizumab), and IgG1κ monoclonal antibody derived from immunizing genetically engineered mice with human TNF‐α (golimumab) has been approved for IBD patients, including Crohn's disease (CD) and ulcerative colitis (UC) and IBD unclassified (IBD‐U).

Although CD and UC are the distinct subtypes of IBD, these diseases present a certain level of similarities, including symptoms, pathological features, immune response, risk factors, and the biological pathways producing TNF‐α. In addition, studies found that up to 3% of CD patients will be reclassified as UC and vice versa after their primary diagnosis, 5–15% of IBD patients classified as IBD‐U, and a small portion of UC patients are later changed to CD or IBD‐U. More importantly, up to 30% of patients do not respond to anti‐TNF‐α blockers, , and the use of vedolizumab (anti‐IL‐12/23) and ustekinumab (anti‐integrin) may be efficacious in many patients that failed anti‐TNF‐α therapy. Thus, there is a clear need to identify potential anti‐TNF‐α treatment pathways in overall IBD patients with a view to better targeting anti‐TNF‐α treatment to more responsive cohorts and to minimize the adverse anti‐TNF‐α treatment effects.

Methods

Search strategy and data collection and integration

A searching strategy for publicly available datasets related to IBD patients who received anti‐TNF‐α therapy was designed for the NCBI Gene Expression Omnibus database dated December 31, 2020, using the keywords, “TNF,” “Tumor Necrosis Factor,” “anti‐TNF,” “anti‐Tumor Necrosis Factor,” “Infliximab,” “Adalimumab,” “Golimumab,” “inflammatory bowel disease,” “IBD,” “ulcerative colitis,” “UC,” “Crohn Disease,” and “CD.” The included datasets have to meet the following inclusion criteria: (i) colonic sample from IBD patients, (ii) transcriptomic data, (iii) raw data are available, (iv) anti‐TNF‐α treatment response status, (v) publicly accessible, and (vi) each of the original studies obtained approval from their local ethics committee and had written informed patient consent. Sample exclusion criteria were as follows: (i) subjects receive therapy other than anti‐TNF‐α, (ii) overlapped subjects, (iii) colonic samples other than large intestine, (iv) transcriptomic data other than Affymetrix, and (v) the posttreatment time point being over 3 months.

The eligible raw microarray datasets were collected and subjected to background correction, normalization, and summarization using the Robust Multichip Average (RMA) algorithm using Affy package version 1.66.0 individually. Mean value of multiple probe sets representing the same gene was calculated. Next, the ComBat function from sva package version 3.36.0 was implemented on the datasets to eliminate the study‐specific batch effects. ,

Composition of immune cells and immune‐related score evaluation

The evaluation of immune microenvironment scores and immune–stroma cell population are calculated using xCELL and ESTIMATE (Estimation of Stromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithms. The immune cell types were evaluated from gene expression profile using five different algorithms, including CIBERSORT, EPIC, MCP‐counter, xCELL, and Deconvolution‐To‐Estimate‐Immune‐Cells (DTEIC). Each of the algorithms was developed using their in‐house or publicly immune cells expression data and different statistical learning approaches. For instance, DTEIC utilized ε‐support vector regression and CIBERSORT applied linear support vector regression , ; MCP‐counter is a single‐sample scoring system, while xCELL requires heterogeneous dataset , ; ESTIMATE utilizes single‐sample Gene Set Enrichment Analysis (ssGSEA) to rank samples on the expression of two different 141 gene sets, and xCELL is based on the sets of cells values calculated from its algorithm. ,

Value for each of the immune cell type estimated from EPIC version 1.1, MCP‐counter version 1.2.0, and ESTIMATE version 1.0.13 were performed in r programming version 4.0.0, DTEIC was operated in Python 3.7 and CIBERSORT and xCELL were calculated by using their corresponding online tools. The source code can be found on the corresponding authors' GitHub page.

Functional enrichment analysis

Identification of differentially expressed genes (DEGs) between responders and nonresponders was calculated using limma package version 3.22.3, and the threshold for the DEGs has a Benjamini–Hochberg adjusted P‐value < 0.05 with absolute log2 fold change ≥ 0.75. EnhancedVolcano package version 1.6.0 was applied for volcano plot. Heatmap was generated by using pheatmap package version 1.0.12. All the packages are applied within the r programming environment. The differentially overexpressed genes were utilized for the pathway enrichment analysis using Metascape (http://metascape.org), a gene enrichment tool for understanding from previously pre‐defined gene sets in different enriched biological themes, including GO terms, Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, BioCarta, and MSigDB. For each gene inputted into the server, the enrichment score was calculated and clustered to match biological signaling pathways. Visualization of the selected pathways utilized Cytoscape version 3.8.0.

Lipopolysaccharide‐induced inflammation in neutrophils

To further confirm the outcomes from the functional enrichment analysis, experimental neutrophil data from Macaca mulatta were applied. Briefly, neutrophils were collected from the target site (at approximately 130 days of gestation). Inflammation was subsequently induced at this site via lipopolysaccharide (LPS) treatment. Subjects were either treated or not treated with adalimumab at 3 and 1 h before LPS, with samples taken at 16 h after LPS. The original study obtained approval from their local ethics committees. The pre‐processed raw RNA‐sequencing count data were retrieved from GSE145918. Then, the normalized log2 + 1 values were calculated using per million reads mapped (CPM) from the count matrix using edgeR version 3.32.0 under r programming environment.

Statistics

Statistical analysis was performed using r version 4.0.0. The pROC package version 1.16.2 in r programming environment was applied to conduct receiver operating characteristic curve analysis to evaluate diagnostic accuracy. The statistical significance was evaluated using a Mann–Whitney U‐test, and Benjamini and Hochberg adjustment was applied for the IBD treatment data. Differences were considered statistically significant at a P‐value of < 0.05, and < 0.05, < 0.01, < 0.001, and < 0.0001 are indicated with one, two, three, and four asterisks, respectively.

Results

Characteristics of studies included in the analysis

After the keyword searching, removal of ineligible and overlapped datasets from the total of 182 records, 5 transcriptomic data, including GSE16879: from the University Hospital of Gasthuisberg, Belgium, with ClinicalTrials.gov number NCT00639821 ; GSE23597: the multicenter, randomized, double‐blind, placebo‐controlled ACT‐1 study between March 2002 and March 2005 with ClinicalTrials.gov number NCT00036439 ; GSE52746: the colonic samples collected between November 2010 and November 2013 from the Department of Gastroenterology, Hospital Clinic of Barcelona, Spain ; GSE73661, UC samples collected from two phase III clinical trials of Vedolizumab (VDZ)—GEMINI I and GEMINI LTS at Leuven University Hospitals, Belgium (the dataset included patients received anti‐TNF‐α blockers) ; and GSE92415: the PURSUIT golimumab study conducted in multicenters, with ClinicalTrials.gov number NCT01988961. The five eligible microarray datasets were normalized and combined, and batch effects were corrected (Fig. S1). Eventually, a total of 449 samples, with 17 771 common gene symbols, were included in this study (Table 1).

Table 1

Summary of the included transcriptomic studies from large intestinal tissues in IBD patients

GSE no. (ref)	Affymetrix platform	Anti‐TNF‐α drug	Study location	IBD	Pretreatment				Posttreatment				Control	Time point (weeks)
					R	NR	Placebo		R	NR	Placebo
							R	NR			R	NR
1687932	Human Genome U133 Plus 2.0	Infliximab	Belgium	UC	8	16	—	—	8	16	—	—	6	4–6
				CD	12	7	—	—	11	7	—	—
2359733	Human Genome U133 Plus 2.0	Infliximab	USA	UC	25	7	5	8	20	7	3	6	—	8
5274634	Human Genome U133 Plus 2.0	Infliximab/adalimumab	Spain	CD	6	1	—	—	7	5	—	—	17	12
7366135	Human Gene 1.0 ST	Infliximab	Belgium	UC	8	15	—	—	8	15	—	—	12	6/12
9241536	HT HG‐U133+ PM	Golimumab	USA	UC	32	27	11	17	29	21	10	15	21	8
Summary Overall = 449				UC	73	65	16	25	65	59	13	21
				CD	18	8	—	—	18	12	—	—
				Total	91	73	16	25	83	71	13	21	56

CD, Crohn's disease; GSE, Gene Set Enrichment; IBD, inflammatory bowel disease; NR, nonresponder; R, responder; TNF‐α, tumor necrosis factor‐α; UC, ulcerative colitis.

Immune microenvironment cell population is significantly higher in nonresponders

The immune microenvironment scores from both ESTIMATE and xCELL identified that the baseline anti‐TNF‐α treatment nonresponders are significantly higher than the responders (ESTIMATE: P < 0.0001; xCELL: P = 0.0003) (Fig. 1a,b). The TNF‐α treatment responders showed a significant drop after their treatments (ESTIMATE: P < 0.0001; xCELL: P = 0.0004) while no significant changes in the nonresponders (ESTIMATE: P = 0.0650; xCELL: P = 0.11) (Fig. 1a,b). The immune and stroma scores (the two calculation factors for the immune microenvironment) are also significantly higher in baseline anti‐TNF‐α treatment nonresponders compared with the responders (Fig. S2A–D).

Figure 1

Microenvironment scores are significantly higher on baseline nonresponders compared with the responders. Immune microenvironment scores evaluated via (a) ESTIMATE and (b) xCELL algorithms. NR, nonresponder; R, responder. The y‐axes are the relative immune microenvironment scores from the corresponding algorithms. P‐value determined by Mann–Whitney U‐test with Benjamini and Hochberg adjustment. Asterisks denote statistically significant differences (*** P < 0.001 and **** P < 0.0001).

Neutrophils, endothelial cells, and B cells are significantly higher in nonresponders

To further our understanding of the immune cell‐type composition between the treatment responders and nonresponders, five different in silico flow cytometry approaches, including CIBERSORT, xCELL, EPIC, MCP‐counter, and DTEIC, were applied (Data S1). Across the algorithms, neutrophils (MCP‐counter: P < 0.0001; xCELL: P = 0.0084; and CIBERSORT: P = 0.0021) (Fig. 2a–c), endothelial cells (MCP‐counter: P = 0.0009; xCELL: P = 0.0183; and EPIC: P = 0.0337) (Fig. 2d–f), and B cells/B linage (MCP‐counter: P = 0.0042; xCELL: P = 0.0251; and EPIC: P = 0.0042) (Fig. 2g–i) are significantly higher in baseline treatment nonresponders compared with the responders. The three‐dimensional plots illustrated that neutrophils, endothelial cells, and B cells from both MCP‐counter and xCELL have positive Pearson correlations with each other using all the eligible data (Figs 2j,k and S3A,B).

Figure 2

Neutrophils, endothelial cells, and B cells are significantly higher on the baseline anti‐TNF‐α treatment nonresponders compared with responders. Immune cell population evaluated in five in silico flow cytometry, and (a–c) neutrophils, (d–f) endothelial cells, and (g–i) B cells can be recognized in three out of five algorithms. B‐cell, endothelial cell, and neutrophil populations are higher on baseline anti‐TNF‐α nonresponders compared with responders. The three‐dimensional plots illustrated that (j, k) neutrophil, endothelial cell, and B‐cell populations from MCP‐counter and xCELL algorithms have positive correlations with each other. The y‐axes are the relative immune cell population abundance from the corresponding algorithms. P‐value determined by Mann–Whitney U‐test. Asterisks denote statistically significant differences (* P < 0.05, ** P < 0.01, *** P < 0.001, and **** P < 0.0001).

Hyperactive chemotaxis contributes to anti‐tumor necrosis factor‐α treatment resistance in inflammatory bowel disease

The pretreatment anti‐TNF‐α subjects (responder: n = 91 and nonresponder: n = 73) were utilized for DEGs analysis and identified a total of 77 DEGs (upregulated genes = 64 and downregulated genes = 13) (Fig. 3a,b and Data S2). Principal component analysis does not have a clear separation between responder and nonresponder subjects in the DEGs (Fig. 3c). The differently upregulated genes compose of several gene families, including cytokines (CCL2, CCL3, CCL4, CXCL13, CXCL5, CXCL6, and CXCL8), chemokines (IL1B, IL6, IL11, and IL24), S100 protein family (S100A8, S100A9, and S100A12), selectin (SELE and SELL), matrix metalloproteinases (MMP1, MMP3, and MMP10), and formyl peptide receptors (FPR1 and FPR2). Metascape pathway enrichment analysis on the 64 highly expressed genes revealed that GO terms with chemotaxis are commonly found from the outcomes and may have a critical role affecting anti‐TNF‐α treatment (GO:0030595: leukocyte chemotaxis, GO:0002688: regulation of leukocyte chemotaxis, GO:1901623: regulation of lymphocyte chemotaxis, and GO:0050918: positive chemotaxis) (Fig. 3d,e and Data S3).

Figure 3

Hyperactive chemotaxis may be involved in anti‐TNF‐α treatment resistance in inflammatory bowel disease. To identify the molecular mechanisms between anti‐TNF‐α blocker responders (n = 91) and nonresponders (n = 73), global gene expression analysis was applied from five combined and normalized microarray datasets. The differentially expressed genes were identified and presented using (a) heatmap (the relative expression values were z‐score transformed) (response: , nonresponder; , responder; IBD: , CD; , UC; dataset: , GSE16879; , GSE23957; , GSE52746; , GSE73661; , GSE92415), (b) volcano plot, and (c) principal component analysis from a total of 64 upregulated and 13 downregulated genes based on the adjusted P‐value < 0.05 with absolute log2 fold change ≥ 0.75. TNF‐α treatment: , nonresponder; , responder. (d) The significantly upregulated genes were utilized for Metascape pathway enrichment analysis from the previously pre‐defined gene set. Enriched terms related to the chemotaxis‐related pathways are underlined. The y‐axis represents the top 20 gene set category, the x‐axis represents −log10 P‐value, and the color intensity of the bar represents the number of genes identified in each hallmark category. (e) The four subsets of enriched terms under the chemotaxis‐related pathways were selected and visualized using Cytoscape.

Tumor necrosis factor‐α blocker does not reduce chemotaxis in lipopolysaccharide‐induced inflammation in neutrophils

To demonstrate our finding in chemotaxis, we used the RNA‐sequencing data from an animal study. Briefly, neutrophils isolated from choriodecidua cells with LPS‐induced inflammation were treated with or without adalimumab. The list of genes from the four chemotaxis enrichment terms was matched with the neutrophil data to obtain the mean expression values of each sample (GO:0030595: leukocyte chemotaxis [31 out of 44 genes], GO:0002688: regulation of leukocyte chemotaxis [18 out of 22 genes], GO:1901623: regulation of lymphocyte chemotaxis [13 out of 16 genes], and GO:0050918: positive chemotaxis [18 out of 21 genes]) (Data S4). The data process workflow is in Figure S4. All the enrichment terms related to chemotaxis were significantly higher in the LPS‐exposed neutrophils; three out of the four enrichment terms do not have significant reduction in the anti‐TNF‐α‐treated group (Fig. 4).

Figure 4

Anti‐tumor necrosis factor‐α (anti‐TNF‐α) blocker does not reduce chemotaxis in lipopolysaccharide (LPS)‐induced inflammation in neutrophils. Mean expression level (across samples) of genes matching indicated GO term. P‐value determined by Mann–Whitney U‐test. Asterisks denote statistically significant differences (* P < 0.05).

Interleukin 13 receptor subunit alpha 2 is a diagnostic biomarker to predict tumor necrosis factor‐α treatment response

In order to find the best potential biomarker to predict anti‐TNF‐α respond IBD patients, receiver operating characteristic curve analysis was applied to all the genes using a for‐loop with pROC package under the r programming environment. Among them, IL13RA2 has the area under the curve of 80.7% (95% confidence interval: 73.8–87.5%) with the best sensitivity of 68.13% and specificity of 84.93% (Table 2, Fig. 5b, and Data S5). IL13RA2 was stand‐alone from the volcano plot (log2FC:1.678, adjusted P < 0.0001) (Fig. 2b), and the expression of IL13RA2 is significantly higher in pretreatment nonresponders compared with the pretreatment responders (P < 0.0001). The responders showed a bigger drop after the treatment compared with the nonresponders (responder: P < 0.0001; nonresponder: P = 0.0037), and the responders restored the expression level to normal control after the treatment (mean ± standard deviation: control: 4.721 ± 1.039; posttreatment, responder: 4.868 ± 1.364; P = 0.4969) (Fig. 5a).

Table 2

The top 15 genes to predict TNF‐α treatment response

Rank	Gene name	AUC (95% CI)	Threshold	Best specificity	Best sensitivity	LR+	LR−
1	IL13RA2	0.807 (0.738–0.875)	6.263	0.849	0.681	4.510	0.222
2	ADGRE2	0.786 (0.716–0.857)	6.079	0.795	0.703	3.429	0.292
3	ADGRL2	0.771 (0.698–0.843)	6.013	0.753	0.736	2.980	0.336
4	HGF	0.768 (0.695–0.841)	5.598	0.753	0.725	2.935	0.341
5	TLR1	0.764 (0.69–0.838)	5.95	0.836	0.637	3.884	0.257
6	NCF2	0.757 (0.683–0.831)	7.696	0.671	0.769	2.337	0.428
7	RGS5	0.757 (0.682–0.832)	7.105	0.767	0.681	2.923	0.342
8	FPR1	0.756 (0.682–0.83)	6.947	0.836	0.626	3.817	0.262
9	TMTC1	0.754 (0.68–0.828)	6.085	0.849	0.593	3.927	0.255
10	CCL4	0.754 (0.679–0.829)	7.955	0.781	0.659	2.304	0.434
11	PDE4B	0.754 (0.679–0.829)	7.681	0.671	0.758	3.009	0.332
12	IGFBP5	0.754 (0.678–0.829)	7.185	0.712	0.725	2.517	0.397
13	SRGN	0.753 (0.678–0.828)	9.849	0.808	0.615	3.203	0.312
14	PAPPA	0.751 (0.675–0.827)	5.946	0.603	0.835	2.103	0.475
15	TMEM71	0.749 (0.674–0.824)	5.672	0.849	0.615	4.073	0.246

AUC, area under the curve; CI, confidence interval; LR+, positive likelihood ratio; LR−, negative likelihood ratio; TNF‐α, tumor necrosis factor‐α.

Figure 5

Interleukin 13 receptor subunit alpha 2 (IL13RA2) can be a diagnostic biomarker to predict TNF‐α treatment response. (a) The expression of IL13RA2 is significantly higher in the pre‐TNF‐α nonresponders compared with the pre‐TNF‐α responders. (b) The expression of IL13RA2 has an area under the curve (AUC) of 80.7% (73.8–87.5%) with a sensitivity (Se) of 68.13% and specificity (Sp) of 84.93%. NR, nonresponder; R, responder. Asterisks denote statistically significant differences (** P < 0.01, *** P < 0.001, and **** P < 0.0001). Statistical significance was determined by two‐tailed Student's t‐test.

Discussion

Treatment resistance of anti‐TNF‐α is a critical issue in IBD patients. By integrating the existing raw data and increasing the statistical power, we revealed that immune microenvironment scores are higher in treatment resistance patients on baseline level (Fig. 1a,b), indicating a higher inflammatory burden in anti‐TNF‐α treatment nonresponders. Further in‐depth analysis uncovered that neutrophils, endothelial cells, and B cells contribute to the change in inflammatory burden (Fig. 3). Next, a total of 64 upregulated genes were identified (Fig. 4a,b), and neutrophil chemotaxis (four out of the top 12 enrichment terms) may contribute to anti‐TNF‐α treatment resistance in IBD patients (Fig. 4d,e). Utilizing an animal study model, mean expression level (across samples) of genes matching the four chemotaxis GO terms is upregulated in LPS‐induced neutrophils but no statistical changes in the three out of the four enrichment terms in the adalimumab‐treated group (Fig. 4).

In a typical inflammatory response, immune cells such as macrophages, dendritic cells, natural killer cells, and T lymphocytes release TNF‐α pro‐inflammatory cytokines, leading to the activation of endothelial cells and neutrophils. The activation of endothelial cells in colonic mucosa enhances vascular permeability and induces the recruitment of immune cells, leading to the activation of chemotaxis. The activation of neutrophils follows the tethering, rolling, crawling, and transmigration process from the blood vessel into the inflamed colonic tissues. When neutrophils engulf invasive gut microbiome, they release granule proteins and chromatin to form neutrophil extracellular traps and secrete antimicrobial peptides to mediate extracellular killing of microbial pathogens. However, hyperactive neutrophils trigger an unrestrained activity of the positive feedback amplification loops, leading to endothelial cells and the surrounding tissues damage, inducing resolution delay (IL‐6, TNF‐α, and IFN‐γ) and chemokines (IL‐8, CCL3, and CCL4), which further the recruitment of neutrophils, monocytes, and macrophages to the inflamed sites. The use of anti‐TNF‐α blockers significantly suppresses the infiltration of neutrophil and B‐cell population in the inflamed mucosa, and suppresses pro‐inflammatory mediators, such as calprotectin (S100A8/A9), IL‐8, IL‐6, and TNF‐α production, , and matched with our finding only in responders (Fig. 2a–c,g–i). The unwanted immunogenicity, however, has a high level of B cells due to the presence of antidrug antibodies. The presence of antidrug antibodies neutralizes, interferes, and/or alters the binding efficacy, as well as pharmacodynamic/pharmacokinetic properties of anti‐TNF‐α monotonical antibodies.

Several S100 calcium‐binding protein family genes are highly expressed and previously studied in anti‐TNF‐α treatment (Fig. 3a,b). Calprotectin is a calcium‐binding protein from the S100A8 and S100A9 monomers, representing up to 40% of neutrophil cytosolic proteins and constantly released from the inflamed region(s). S100A12, also known as calgranulin C, is released from neutrophils and participates in pro‐inflammatory process via the activation of the NF‐κB. A small‐scale study reported that the fecal calprotectin test (commonly used to distinguish between irritable bowel syndrome and IBD) and S100A12 may predict relapse after 1 year of infliximab treatment, while another fecal calprotectin study did not find the difference.

IL13RA2 is stand‐alone in the volcano plot with the highest fold change and the lowest P‐value (Fig. 3b and Data S4) and has the best area under the curve (80.7%, 95% confidence interval: 73.8–87.5%) outcome (Table 2, Fig. 5b, and Data S5). Early studies uncovered that IL13RA2 is active in mucosal biopsies on the UC or CD anti‐TNF‐α treatment nonresponders compared with the responders. , A small‐scale study demonstrated that soluble IL13RA2 protein cannot be detected in serum, and tissue expression of IL13RA2 could predict anti‐TNF‐α treatment in CD patients. IL13RA2 is a decoy receptor enable to bind IL‐13 cytokine, diminishes its JAK1/STAT6‐mediated effector functions, and activates activator protein 1 (AP‐1) to induce the secretion of TGF‐β. , The IL‐13 pathway is also dependent on the production of TNF‐α. Several IL‐13 targeting drugs have been tested to inhibit hyperactive immune response on Th2‐driven inflammatory diseases. However, insufficient protection was demonstrated by the phase IIa anrukinzumab (an IL‐13 monoclonal antibody) clinical trial on UC patients. Thus, blocking the IL‐13 pathway via IL13RA2 could be a new approach for treating IBD patients. IL13RA2 knockout mice in DSS‐induced acute colitis model showed a better recovery rate compared with the wild‐type mice and negatively regulate epithelial/mucosal healing. By neutralizing IL13RA2 in DSS‐induced IBD murine model using a monoclonal antibody, it presented a speedy recovery compared with the control group.

The study here identified many strengths but should be considered in the context of shortcomings. Firstly, we only focused on data from large intestine and eliminated ileum data from the Arijs et al. study due to the low number of ileum samples that can be integrated, and also reduce the gene expression variation between two different organ sites for downstream analysis. Secondly, our comparison does not include studies from vedolizumab and ustekinumab as it has limited datasets available online. Thirdly, the anti‐TNF‐α response criteria and the determination time points are slightly different between studies. As we can only rely on the information provided by the authors, and thus, our study has to accept the potential bias. Fourthly, anti‐TNF‐α is broadly used in colitis‐based diseases with a high percentage of treatment failure, and the diagnosis criteria of UC/CD/IBD‐U on each of the included studies may be slightly different with a certain percentage of misclassification. Therefore, our priority is to find the common patterns to minimize the treatment resistance rate in this study. Last but not least, the animal study in neutrophils is not from the colonic tissue sites, and some of the chemotactic factor markers such as IL‐8/CXCL8 and CSF3 have a significant reduction after the adalimumab treatment. We believed that some single markers may not represent a whole picture of chemotaxis. Thus, in the early future, IBD subtype analysis and more in‐depth study in the relation to hyperactive chemotaxis are needed.

In conclusion, nonresponders presented higher populations of neutrophils, endothelial cells, and B cells compared to responders at baseline level. IL13RA2 is a potential biomarker to predict anti‐TNF‐α treatment response.

Availability of data and materials

All the transcriptomics data were retrieved from publicly available datasets from the NCBI Gene Expression Omnibus (GEO) database. All the in silico cell sorting algorithms used in this study were based on the default settings recommended by the corresponding authors, either from web‐based or retrieved from the corresponding authors' GitHub page for academic research purposes.

Figure S1. (A) Before and (B) after the batch effects correction of the included samples.

Click here for additional data file.

Figure S2. Immune‐Stroma scores are significantly higher on baseline non‐responders compared to the responders. Scores evaluated using both xCELL and ESTIMATE algorithms for (A, B) immune and (C, D) stroma/stromal scores are significantly higher in the pre‐anti‐TNFα treatment nonresponders compared to responders. NR: non‐responder; R: responder. The y‐axes are the relative immune microenvironment scores from the corresponding algorithms. P‐value determines by Mann–Whitney U‐test with Benjamini and Hochberg adjustment. Asterisks denote statistically significant differences (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

Click here for additional data file.

Figure S3. Pearson correlation matrix of neutrophils, endothelial cells and B cells from (A) xCELL and (B) MCP‐Counter algorithms. Blue indicates positive correlation, and red indicate negative correlation. Darker colours are associated with stronger correlation coefficients.

Click here for additional data file.

Figure S4. The RNA‐sequencing data process workflow in LPS‐induced inflammation in neutrophils.

Click here for additional data file.

Data S1. Supporting Information

Click here for additional data file.

Data S2. Supporting Information

Click here for additional data file.

Data S3. Supporting Information

Click here for additional data file.

Data S4. Supporting Information

Click here for additional data file.

Data S5. Supporting Information

Click here for additional data file.