Figure S1 displayed the flowchart for the study. Gene expression profiles of ulcerative colitis were retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). A total of 803 RNA-seq samples were downloaded from datasets GSE87473, GSE92415, and GSE206285, including untreated 743 UC samples and 60 healthy control samples. It should be noticed that though 364 patients with UC were treated with ustekinumab (data extracted from database GSE206285), RNA extraction and microarray analysis of biopsy samples collected when the data set was uploaded at baseline. Baseline is the period of time in a clinical study when a patient has been screened to for the study but has not yet started medication. To convert the probe expression matrix to the gene expression matrix, we used the platform annotation file. The average of multiple probes corresponding to the same gene is calculated. To remove the batch effect between datasets and merge them, we used the “ComBat” function in the SVA package. Equally, the CAC dataset was obtained from dataset GSE37283. For validation, GSE53306 and GSE16879 datasets were also downloaded and processed. The included datasets and baseline clinical characteristics of patients with UC were provided in Supplementary Table S1 .

Based on RNA-seq data obtained from GEO (743 UC samples and 60 normal samples), we used the LIMMA package to identify differentially expressed genes (DEGs). Volcano and heatmap visualizations of the results were created using the ggplot2 package in R. DEGs were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the ClusterProfiler package in R. GO analyses were conducted for cellular components (CC), biological processes (BP), and molecular functions (MF).

To calculate immune cell infiltration and examine the disease immune microenvironment, we used the algorithm CIBERSORT and iterated 1000 times (27). Using the ggplot2 package, we visualized the results of immune cell infiltration. For visualization of immune cell infiltration correlations, Corrplot was used to calculate Spearman analysis correlations between immune cells.

The WGCNA package was used to identify genes associated with immune cell infiltration. The gene expression matrix outliers were filtered using hierarchical clustering. Based on the expression data, Pearson correlation coefficients were calculated between two genes to construct a similarity matrix. Using a soft threshold of 15, the similarity matrix was transformed into an adjacency matrix and then into a topology matrix with a topological overlap measure (TOM). Modules are identified using a dynamic pruning tree and clustered based on 1-TOM distance. Setting the merge threshold to 0.25 resulted in the identification of 14 modules. By intersecting DEGs with genes associated with neutrophils (red modules), we identified neutrophil-associated differentially expressed genes (NADEGs). A search for hub genes and functional analysis of NADEGs in Cytoscape (3.9.1) using Cytohubba and Cluego.

We performed an unsupervised cluster analysis using ConsensuClusterPlus to identify different neutrophil infiltrating subtypes. Functional annotation and GSVA were used to identify biological differences between subtypes. GSVA was performed on gene sets from the MSigDB database (C2.Cp.ke.v7.2) with statistical significance set at adj.P<0.05. To assess immune cell infiltration, ssGSEA (single sample gene set enrichment analysis) was performed. In addition, we used the ESTIMATION algorithm to assess the stromal and immune score in each sample. Additionally, we compared clinical characteristics (age and disease extent) between subtypes.

Using the LIMMA package, the DEGs between various subtypes were found, and the criterion was | log2 FC | >0.585 with an adjusted p-value< 0.05. To test the predictive value of the identified biomarkers, we assessed the diagnostic validity of UC and control sample differentiation using the area under the ROC curve (AUC) values. Using Spearmen analysis, the relationship between the detected biomarkers and the quantity of invading immune cells was investigated.

The experiments were conducted with 8-10 weeks male mice under specific pathogen-free conditions (SPF). The AOM/DSS model was constructed using C57BL/6 mice, which is well-established for replicating the immune response in colitis (28). All mice were purchased from Vital River (Beijing, China). The Animal Health Research Committee of the First Hospital of Jilin University approved all animal experiments. The mice were divided into three groups, control, DSS and AOM/DSS, with three mice in each group. AOM/DSS models were constructed as follows: mice were intraperitoneally injected with 10 mg/kg of AOM on the first day. A week later, mice were fed 2% DSS solution (40-50 KDa DSS (Sigma Aldrich) dissolved in water at a final concentration of 2%) for 5 consecutive days, and so on for 3 cycles. The mice in the DSS group were fed with 2% DSS for 5 consecutive days and were regularly fed water for the remaining days. After 6 weeks, mice were euthanized and colonic tissues were collected for hematoxylin and eosin (H&E) staining, flow cytometric analysis, and RT-qPCR.

Colon tissues were taken from different groups of mouse colons, and we extracted total RNA from the tissues using the Animal Total RNA Isolation kit (FOREGENE, China) and determined its purity and concentration. We used the Uni RT&qPCR kit (transgen, China) to reverse convert RNA to cDNA and performed RT-qPCR on an instrument (ABI QuantStudio 3, USA). We used β-actin as an internal reference and 2-ΔΔCt to calculate the relative gene expression levels and Graphpad 9.0 to visualize the data. A detailed list of primer sequences used for RT-qPCR is shown in Table S2 .

Tissues from the colon were embedded in paraffin after being treated with 4 percent paraformaldehyde (PFA) overnight. To prepare a glass slide for H&E staining, 5 um-thick pieces of the colon were cut out and set flat on it.

Mice were sacrificed to obtain tissues and cells for detection. 5ml cold PBS was injected into the enterocele of mice, gently shaken the body, and collected PBS into 15ml centrifuge tubes to obtain peritoneal irrigation fluid (PIF). To obtain neutrophils in peripheral blood (PB), cells were separated by density gradient centrifugation (Mouse peripheral blood lymphocyte isolation kit, solarbio, China), and neutrophils were obtained after lysed red blood cells in the sedimentation. Peyer’s patches (PP) were carefully detached from the small intestine, and mesenteric lymph nodes (MLN) were carefully detached from the colon and placed in PBS. To collect lamina propria mononuclear cells (LPMC), the colon of mice was cut into small pieces (about 1mm3 in volume), mixed tissue pieces with digestion buffer (PBS containing 5 mM EDTA and penicillin/streptomycin) and shaken for 20 mins to remove intraepithelial lymphocytes. Then, the tissue pieces were digested with 1 mg/ml collagenase IV and 10 U/ml DNase I at 37°C for 1 h. All cells were obtained by passing the tissues through a 70 mm strainer.

5 * 105 cells were washed twice in FACS buffer and stained for CD45 Pacific Blue (S18009F), CD11b APC/Cyanine7 (M1/70), Ly 6g PerCP/Cyanine5.5 (1A8) for 30 mins. For intracellular stain, cells were treated with the transcription factor staining buffer set according to the manufacturer’s specification (eBioscience, USA). Then, cells were stained with the following antibodies: rabbit anti-mouse myeloperoxidase (Abcam, England), STAT3 Phospho (Tyr705) PE antibodies for 30 mins, then cells were washed with FACS buffer and stained with goat anti-rabbit IgG H&L Alexa Fluor^® 488 for 30 mins. Then, cells were washed twice for flow cytometric analysis. Antibodies for CD45, CD11b, Ly6g, and pSTAT3 were purchased from Biolegend. Antibodies for myeloperoxidase and goat anti-rabbit IgG H&L were purchased from Abcam. Flow cytometric analyses were conducted using an Arial flow cytometer (BD Biosciences). Data were evaluated using FlowJo software (Version 10; FlowJo).

The statistical analyses were conducted using R (version 4.2.0). Continuous variables were compared between groups, and normally distributed variables were compared using the t-test. Diagnostic biomarkers were evaluated using ROC curve analysis. Spearman correlation was used to analyze the relationship between infiltrating immune cells and gene biomarkers. All statistical analyses were two-sided, and P < 0.05 was considered statistically significant.

We collected 743 UC patients and 60 healthy controls from the GSE87473, GSE92415, and GSE206285 datasets. Using the battle algorithm in the SVA package, the combined dataset was obtained by removing batch effects from the GEO dataset. DEGs between UC and healthy controls were analyzed by the LIMMA package (|LogFC|>2, p<0.05). There were found to be 65 DEGs in total, of which 22 had down-regulation and 44 had upregulation ( Table S3 ). Heatmap and volcano plots illustrated these upregulated and downregulated genes ( Figures 1A, B ). The biological features of these DEGs were then ascertained by the GO analysis and KEGG analysis. GO functional enrichment analysis revealed that DEGs were significantly enriched in granulocyte chemotaxis, neutrophil chemotaxis, and neutrophil migration in the BP category. There was a significant enrichment in the cytoplasmic vesicle lumen and secretory granule lumen in the CC category. The MF terms were enriched with cytokine activity, chemokine receptor binding, and cytokine receptor binding ( Figure 1C and Table 1 ). KEGG pathway analysis revealed regulated terms (FDR<0.05), including IL-17 signaling pathway, Cytokine-cytokine receptor interaction, NF-κB signaling pathways, and TNF signaling pathway ( Figure 1D and Table 2 ). We also compared the respective significantly enriched pathways in the healthy control and UC cohorts. As shown in Figures 1E, F , immune-related pathways were enriched in the UC cohort compared to healthy controls. These results suggested that inflammation-related pathways are enriched in UC patients.

Based on CIBERSORT estimates, we determined the relative infiltration abundance of immune cell types in UC and healthy control samples. The majority of immune cells infiltrating the UC were more abundant than those in the control ( Figure 2A and Table S4 ). In UC samples, neutrophils, mast cells, activated CD4 T cells, T cell gamma delta, M0 macrophages, and M1 macrophages were the major infiltrating immune cells, suggesting the involvement of multiple immune cells. Figure 2B showed the correlation between immune cells. Next, we used the WGCNA package to construct co-expression networks and calculated average linkages and Pearson correlation values to cluster the samples. The logarithm of node connectivity K (K) was negatively correlated with the logarithm of node probabilities (P(K)) with a correlation coefficient above 0.85. Based on the gene expression matrix and the proportion of immune cell types, the optimal soft threshold power was 15 ( Figure S2 ). We created a hierarchical clustering tree using a dynamic hybrid cut method, where genes with similar expressions were tightly grouped to form a branch representing a gene module. Different branches of this tree represent different gene modules, with high co-expression levels of genes within the modules. As a result, 14 main modules were identified ( Figure 2C ). In order to select important modules to research the relationship between modules and immune cells, we found that the red module had the strongest correlation with neutrophils (R=0.83, P=3e-193) based on Pearson correlation coefficients between modules and sample characteristics. The infiltration characteristics of neutrophils in UC were further explored in combination with the enrichment results of DEGs. Therefore, 24 neutrophil-associated differentially expressed genes (NADEGs) were obtained by focusing on the red module and intersecting the genes in the red module with DEGs ( Figure 2D ). The core genes in the NADEGs were then identified using the CytoHubba plugin in Cytoscape. Using the Degree algorithm, we identified the top 5 hub genes as IL1B, CXCL8, CXCL1, TIMP1, and FCGR3B ( Figure 2E ). ClueGO then analyzed the NADEGs for pathway and process enrichment. RAGE receptor binding, cellular response to UV-A, IL-17 signaling pathway, and neutrophil chemotaxis were highly enriched ( Figure 2F ).

Using the ConsensusClusterPlus package in R software, different neutrophil infiltration subtypes were identified using consensus clustering based on 24 NADEGs, and two clusters (cluster A and cluster B) were identified ( Figure 3A ). There were 403 samples in cluster A and 340 samples in cluster B. PCA showed that NADEGs could completely distinguish between the two subtypes ( Figure 3B ). Then, the heatmap and boxplot plot illustrated the differential expression levels of the 24 NADEGs between the two clusters ( Figures 3C, D ). The expression levels in cluster A were all lower than those in cluster B. GSVA enrichment analysis was then conducted to understand the biological changes behind the different clusters. Figure 3E and Table S5 showed that subtype B is significantly enriched in immune-related pathways such as complement and coagulation cascades, cytokine-cytokine receptor interaction, chemokine signaling pathway, NOD-like receptor signaling pathway, natural killer cell-mediated cytotoxicity, TOLL like receptor signaling pathway, and JAK-STAT signaling pathway. Next, we compared the differences in immune cell infiltration between the two clusters using the ssGSEA method ( Table S6 ). In Figure 3F , there was a higher infiltration of immune cells, especially neutrophils and various types of macrophages, in cluster B compared to cluster A. The ESTIMATE package was also used to assess the microenvironment score (immune score, stromal score, and estimated score) for both subtypes ( Figure 3G ).

Age and disease extent are risk factors for colitis-association colorectal cancer, and patients with UC who are younger and have more extensive lesions are more likely to convert to CAC (29, 30). We observed a different distribution of age and extent of disease between the two subtypes. As shown in Figure 4A , a greater proportion of patients in cluster A were over 30 years of age (64.30%), while more than half of the patients under 30 years of age were from cluster B (55.30%). In terms of the extent of the disease, more patients with extensive lesions were from cluster B (62.50%), while 66.67% of patients in cluster A presented with limited lesions ( Figure 4B ). Following that, we examined the GSE37283 dataset for DEGs between UC and CAC. Using the LIMMA package, we set |LogFC>2|, p<0.05 to identify 14 genes associated with CAC ( Figure S3 ). We next compared the differences in the expression of the CAC-associated genes between different neutrophil subtypes. We found that most CAC-associated genes were highly expressed in cluster B, and therefore we concluded that cluster B with highly infiltrated neutrophils was more likely to transform from UC to CAC ( Figure 4C ). The major histocompatibility complex and T-cell stimulating factor were investigated concerning the two subtypes. Major histocompatibility complex and T-cell stimulating factor expression levels were higher in cluster B, but not in TNFRSF14 and TNFRSF25 ( Figures 4D, E ). Recent studies have shown that highly infiltrated neutrophils in UC could form neutrophil extracellular traps that contribute to the transformation of colitis into cancer (24). From the review, we collected mostly information on ligands and receptors that promote the NET formation, downstream signals, and molecules identified as NET adhesion frameworks (31). According to our results ( Figure 4F ), almost genes related to neutrophil extracellular traps were highly expressed in subtype B, suggesting that more neutrophil extracellular traps might be present in cluster B (32). Furthermore, we compared the differential expression of immune checkpoints between the two subtypes, and we found that most immune checkpoints, including PD-L1 and CTLA4, were highly expressed in subtype B ( Figure 4G ).

To better find genes that could diagnose UC or even CAC, we identified 5 DEGs (AQP9, FCGR3B, GPR109B, PROK2, S100A12) between neutrophil-associated subtypes using the LIMMA package. Firstly, we examined the correlation between these 5 genes, which was highly correlated as illustrated in Figure 5A . Next, we accessed the correlations between these genes and immune cells. These genes were also correlated with other immune cells besides neutrophils ( Figure 5B ). To account for the ability of each gene to distinguish between UC and normal samples, the area under the ROC curve (AUC) was calculated to assess the diagnostic value of these genes ( Figure 5C ). Next, we compared the expression of 5 genes in the GSE92415 and GSE16879 cohort receiving anti-TNF-α treatment to determine whether the effect of anti-TNF-α treatment could be assessed. As shown in Figure 5D and Figure S4 , the expression levels of all 5 genes were decreased in the response-yes group, suggesting that these genes could also be used as markers to assess the response of anti-TNF-α treatment. Furthermore, we investigated whether these five genes could be used to distinguish between patients with active UC and those with remission UC in the GSE53306 cohort. Apart from FCGR3B, the expression of the remaining genes was higher in active UC, as shown in Figure 5E .

DSS-induced colitis and AOM/DSS-induced CAC were used to further validate the expression of biomarker genes. The colonic length of mice with different models exhibited significant differences from control mice ( Figures 6A, B ). Histological analysis indicated the presence of damage to the intestinal epithelium in mice with enterocolitis, and mice in the AOM/DSS group had significant dysplasia ( Figure 6C ). We also examined the expression levels of the IL-6, TNF-α, and biomarker genes using RT-qPCR. Gene expression was increased in the DSS group compared to the Ctrl group and was even more significantly upregulated in the mice with AOM/DSS group ( Figure 6D ).

Then we proceeded to further investigate the variation in neutrophil infiltration with different mouse models (Control, DSS, and AOM/DSS). The initiation of colitis resulted in a mass recruitment of CD11b+ Ly 6g+ neutrophils in PB, MLN, and PIF ( Figures 7A, B, E, F , and Figure S5 ). The proportion of neutrophils in PP of mice with AOM/DSS was higher than that of DSS ( Figure 7C ) but declined in LP ( Figure 7D ). Next, we compared the expression of MPO in neutrophils by flow analysis to compare the differences between the different groups (33, 34). In addition, DSS mice showed an increased level of MPO in PB, PP, LP, MLN, and PIF ( Figures 7G–L ), and even higher in AOM/DSS mice ( Figures 7H–J, L ). JAK-STAT signaling pathway has been demonstrated to play a potential role in inflammation and tumor development. The elevated levels of pSTAT3 in PB, LP, MLN and PIF of AOM/DSS mice ( Figures 7M, O–Q ) were consistent with the results of GSVA ( Figure 3E ). However, for DSS mice, the pSTAT3 levels in PP did not change significantly ( Figure 7N ). These results established an increased percentage of neutrophils and MPO, pSTAT3 expression in DSS and AOM/DSS mouse models. These results also indicated that the neutrophils and the JAK-STAT pathway might be involved in the transformation of colitis to colitis-associated colorectal cancer.