Firstly, to identify the transcriptomic data of CRC samples, we extensively searched the datasets from the GEO database and downloaded three datasets through the geoquery package (22) in the R software. The GSE110223 dataset included 13 pairs of patients with confirmed colorectal adenocarcinomas and normal colon tissues. The GSE110224 and GSE113513 datasets contained 17 and 14 cases of CRC and matched normal samples, respectively. We also obtained the clinical information and mRNA data (including 480 colon cancer tissues and 41 normal tissues) from the TCGA data portal. The total TCGA queue was equally classified into two groups, namely, TCGA validation set 1 (n = 219) and TCGA validation set 2 (n = 219). Furthermore, 496 CFRGs and 2,484 IRGs were acquired respectively from the GeneCards and Analysis Portal (ImmPort) database (23).

Before screening the differentially expressed genes (DEGs), all datasets were first normalized using the “limma” package and PCA diagrams were drawn to evaluate the differences between groups. We finally identified the DEGs with the cutoff conditions of |log₂FC| >1 and p < 0.05. The volcano plot of DEIRGs and Venn diagram were plotted by the R “ggplot2” package for visualization.

To elucidate the latent functions of DECFRGs, we worked with the org.Hs.eg.db package to transfer gene ID and used the clusterProfiler package (24) (version: 3.14.3) to clarify the potential functional enrichment analysis. Enrichment analysis includes KEGG pathway analyses and GO analyses, including BP, CC, and MF. p < 0.05 denotes statistically significant differences. To visualize the above enrichment analysis, the ggplot2 package was used to represent the top items per ontology with bubble diagram and KEGG pathways with a network diagram.

Univariate Cox regression analysis was performed to calculate the association between predicted gene expression levels and OS of CRC patients. Genes were suggested to provide a meaningful prognostic latent capacity at a p-value <0.05. Then, we built a prognostic model for OS based on the LASSO regression via the R package “glmnet” (25, 26). Finally, we established a survival risk model that divided the patients into low-risk and high-risk groups after delimiting prognostic genes and their LASSO coefficients. Kaplan–Meier plots were used to evaluate the prognosis of CRC patients with a cutoff value of p < 0.05. Additionally, the time-dependent ROC curve at specified time points (1 year, 3 years, and 5 years) was conducted to access the accuracy and sensitivity of survival prediction. The risk score = confident * expression of gene1 + confident * expression of gene2 + confident * expression of gene3 + confident * expression of gene4 +… + confident *expression of gene N.

TIMER is a tool to analyze and visualize the correlation between levels of immune infiltrate and diverse tumor types. The correlation between CXCL8 expression and tumor-infiltrating immune cells was investigated through TIMER and ssGSEA in the GSVA package. Next, we used the R package “Estimate” to obtain the three scores of CXCL8, including immune score stromal score, stromal score, and estimate score. The TISIDB database provided us the associations for CXCL8 with immunomodulators, lymphocytes, and drugs. Furthermore, immune checkpoint molecules play a crucial role in cancer therapy. Thus, we explored the relationship between CXCL8 and the CTLA-4 inhibitor ipilimumab, PD-1 inhibitors, PD-L1 inhibitors, and so on.

Human colorectal tumor specimens were provided by the Department of Oncology Surgery, The Second Clinical Medical College of Lanzhou University. All experiments involving human tissues complied with the principles of the Declaration of Helsinki and have been approved by the Institutional Review Board of the Second Hospital of Lanzhou University. TRIzol was used to obtain total RNA extraction. SuperScript™ III First-Strand Synthesis SuperMix was used to synthesize for qRT-PCR. Then, we perform real-time PCR with Power SYBR^® Green PCR Master Mix kit. The primers used are shown in Supplementary Table 2 .

Extraction of total protein from CRC cells and tissue samples is carried out by RIPA and a protease inhibitor 99:1 mixture. After centrifugation, the collected supernatant was further analyzed. Protein concentration was determined by BCA assay, analysis, and denaturation. Then, total protein was loaded onto 12% PAGE gels and transferred onto PVDF membranes. Afterwards, the membrane was blocked with 5% nonfat milk for 2 h at 37°C on a shaker and incubated with the primary CXCL8/IL-8 antibodies and GAPDH antibodies at 4°C overnight on a shaker. Then, the secondary anti-mouse or anti-rabbit antibodies were used to incubate the membranes at room temperature for 1 h, followed by visual analysis using the ECL developer on the Tianneng exposure imaging system (TANO, Shanghai, China).

Whole tissues were paraffin embedded and serially cut into 4-μm-thick sections, The sections were baked in an oven, dewaxed in xylene, soaked in alcohol, and washed with distilled water. After antigen retrieval, sample peroxidase activity was blocked with 3% H₂O₂. Afterwards, primary antibodies (anti-CXCL8/IL-8, CD3, CD4, CD8, CD20, CD56, CD68, and CD163) were used to incubate sections overnight at 4°C. The next day, sections were incubated with anti-rabbit antibody for 30 min. Last, sections were stained with 3,3’-diaminobenzidine.

Statistical analyses were made with the GraphPad Prism 8.0.1 software. Statistical significance was determined using unpaired two-tailed Student’s t-test: *p < 0.05.

As shown in Figure 1A , 700 DEGs, namely 291 upregulated genes and 409 downregulated genes, were identified in the GSE110224 dataset. Meanwhile, we screened the DEGs in the GSE110223 dataset with a threshold at |logFC| > 1 and adjusted p at < 0.05, for a total of 480 genes, 182 of which were upregulated and 298 were downregulated ( Figure 1B ). In the GSE113513 dataset, there are 1,539 DEGs (namely 585 upregulated genes and 954 downregulated genes) ( Figure 1C ). The threshold was set at |logFC| > 1 and adjusted p at < 0.05. After intersecting the DEGs with CFRGs, we finally came to 17 genes ( Figure 1D ). To explore the biological function of these 17 DECFRGs, we performed GO and KEGG enrichment analyses using the clusterProfiler package in R.

KEGG enrichment analysis indicated that the signal pathways most related to the DECFRGs were rheumatoid arthritis, cytokine–cytokine receptor interaction, and IL-17 signaling pathway ( Figure 1E ). Furthermore, the top 4 enriched terms in the aspect of BP were extracellular matrix disassembly, response to nutrient levels, collagen catabolic process, and humoral immune response. The enriched GO items included endoplasmic reticulum lumen, platelet alpha granule lumen, and platelet alpha granule in the CC group. DECFRGs were mostly enriched in serine-type endopeptidase activity, metallopeptidase activity, and cytokine activity in MF ( Figure 1F ). Detailed GO and KEGG information is provided in Supplementary Table 1 . Finally, nine DEIRGs were identified through the Venn diagram ( Figure 1G ).

In the TCGA dataset containing 521 patients, duplicate samples and non-clinical samples were excluded. Univariate analysis was conducted and LASSO regression was adopted to conduct IRGPI. Lambda.min was selected to produce the model with a higher accuracy rate ( Figure 2A ); thus, we ended up with five prognostic genes: CXCL8, MMP12, GDF15, SPP1, and NR3C2 ( Figure 2B ). Univariate analysis showed that CXCL8 is remarkably related to OS of patients with COAD ( Figure 2C ). After multiplying gene expression levels with LASSO coefficients, we obtained the risk score = (−0.0843825) * expression of CXCL8 + (−0.043955) * expression of MMP12 + (−0.127046) * expression of GDF15 + 0.09601238 * expression of SPP1 + (−0.0656746) * expression of NR3C2.

The expression levels of these selected four genes were significantly increased in tumor tissue than normal tissue, and another gene had a lower expression compared to normal samples ( Figure S1 ). To confirm the expression of hub genes, we also collected tumor and normal tissues from Lanzhou Second Hospital and performed RT-qPCR. The results were consistent with online data. The results suggested that CXCL8, MMP12, GDF15, and SPP1 had a higher expression, while NR3C2 had a lower expression compared with the normal group ( Figure 3 ).

First, according to the median value of IRGPI risk score, we calculated each patient’s risk score and classified the critical values of the high- and low-risk group in the TCGA training set ( Figures 4A, D ), in TCGA validation set 1 ( Figures 4B, E ), and in TCGA validation set 2 ( Figures 4C, F ), and the IRGPI thresholds are −1.354126396, −1.460113703, and −1.248139089, respectively. In the training set, patients in the high-risk group had a significantly poorer OS, and those included in the low-risk group had a higher OS ( Figure 4J , p < 0.05). As shown in Figure 4G , the time-dependent ROC curve respectively indicated an AUC at 1 year, 3 years, and 5 years of 0.618, 0.613, and 0.587, respectively. When dividing the entire TGCA dataset into validation sets, we discovered that IRGPI also presented a high correctness of the predictive OS in TCGA validation set 1 and TCGA validation set 2. The validation datasets were also applied to the training set, and the AUC for 1-year, 3-year, and 5-year survival was 0.636, 0.587, and 0.653 in TCGA validation set 1, respectively ( Figures 4H, K ). Finally, validation set 2 had a p-value of <0.02 and an AUC of 0.535, 0.550, and 0.536 for 1-year, 3-year, and 5-year survival, correspondingly ( Figures 4I, L ). The outcome suggested that IRGPI could be a good dependable target and better than other features.

In Table 1 , we explored the relationship between clinicopathologic factors and OS via univariate analyses. The results demonstrated that there was a significant relationship between clinicopathologic factors (including age, T/N/M stage, CEA levels, and OS; p < 0.05). Additionally, we built a prognostic nomogram to forecast the survival risk of each CRC patient ( Figure 5A ). Although there was a good consistency between predictive and actual 1-year survival in the entire TCGA set, there was a bad consistency between predictive and actual 3- and 5-year survival ( Figures 5B–D ).

CXCL8, as one of the genes in IRGPI, is one of the first and most intensively studied chemokines and is secreted by different cell types (27, 28). We have previously predicted that CXCL8 associated with immunity. Therefore, the correlation of CXCL8 with immune cell infiltration was evaluated. Firstly, the stromal score and immune score of CXCL8 showed a higher infiltration of stromal and immune cells when CXCL8 was highly expressed. Furthermore, the estimated score showed higher CXCL8 had a lower purity of tumor cells ( Figure 6A , p < 0.0001). Then, we researched the relation between these immune cells and CXCL8. Among the 24 types of immune cells, the relative proportion of B cells, T cells, CD8 T cells, macrophages, and Treg had a significantly positive correlation with CXCL8 ( Figure 6B ). We further used TIMER to assess the immune infiltration levels of CXCL8. The outcomes suggested that there is a notable correlation between CXCL8 and tumor purity, CD8+ T cells, macrophages, neutrophils, and DCs ( Figure 6C ). Together, those results reveal the strong role of CXCL8 to reflect tumor microenvironment.

In order to further understand the correlation between CXCL8 and immune infiltration, we studied the relationship between CXCL8 expression and various immune characteristics, including TILs, immunomodulators, chemokines, and drugs. We obtained the correlation analysis between CXCL8 expression and several immune features from the TISIDB database. Figure 7A indicates the positive associations between CXCL8 and TILs, including Imm_B_abundance, Act_B_abundance, Act_CT4_abundance, Act_CT8_abundance, Tem_CD4_abundance, and Tcm_CD4_abundance. Figure 7B shows positive relationships between immune inhibitors including BTLA, IDO1, CSF1R, CD96, CD244, and CD274 with CXCL8. Figure 7C shows positive correlations between CXCL8 expression and chemokines, including CCL2_exp, CCL17_exp, CCL11_exp, CXCL11_exp, CXCL12_exp, and CCL14_exp. Figure 7D shows correlations between CXCL8 expression and drugs, including ABT510 targeting FGF2, HGF, VEGFA, CXCL8, Tapinaroftargeting IL2, IL6, IL12B, CXCL8, Rivanicline targeting CHRNB2,CXCL8 and MDX-018 targeting CXCL8, respectively. Thus, it could be confirmed that CXCL8 widely affected immune infiltration in the tumor microenvironment.

To further investigate the response of immune checkpoints, we conducted the relationship between CXCL8 and some immune checkpoints using the TCGA database. The major immune checkpoints include CTLA-4, PD-1, PD-L1, and PD-L2. After investigating the expression of some key immune checkpoints, the results indicated that the expression of CTLA-4 ( Figure 8A ), TIGIT ( Figure 8B ), BTLA ( Figure 8C ), PD-1 ( Figure 8D ), PD-L2 ( Figure 8E ), TIM-3 ( Figure 8F ), CD274 ( Figure 8G ), CD27 ( Figure 8H ) and LAG3 ( Figure 8I ) was positively related with expression of CXCL8 (p<0.05). Among these immune checkpoints, CXCL8 was highly correlated with PD-L2, TIM-3 and CD274. Taken together, these results suggested that CXCL8 had a strong response to immune response and immunotherapy.

To identify the expression of CXCL8, we determined the CXCL8 protein expression level through Western blot both in vivo and in vitro. Consistent with the previous experiments, Western blot suggested that the protein expression of CXCL8 is higher in tumor tissues than normal tissues ( Figure 9A , p < 0.05). Furthermore, we identified CXCL8 expression in several types of COAD cells. As shown in Figure 9B , CXCL8 levels were higher in all COAD cells than in control cells. To evaluate the role of CXCL8 in immune infiltration, we evaluated CXCL8 and various immune cell markers (including CD3, CD4, CD8, CD20, CD56, CD68, and CD163) via immunohistochemistry (IHC). As shown in Figure 9C . CXCL8 had a positive correlation with CD3, CD4, CD8, CD20 (MS4A1), and CD163 ( Figure 9D , p < 0.05) in COAD; however, there is no significant correlation with CD56 and CD68. The results verified that CXCL8 was closely related with immune cell infiltration in COAD.