Comprehensive analysis of the prognosis and immune infiltration for CXC chemokines in colorectal cancer

Ethics statement

The collection of colorectal tissues was approved by the Ethics Committee of the Taizhou People’s Hospital, Suzhou, China (project no. KY 202005401), and was conducted according to the relevant guidelines. The datasets were obtained from the relevant published literature.

Oncomine analysis

The transcription levels of CXC chemokines in different cancers were downloaded from a cancer microarray database (Oncomine, https://www.oncomine.org). The differences in expression levels of CXC chemokines at the mRNA level in CRC and corresponding normal tissues were evaluated by Student’s t-test. The cutoffs for the identification of differentially expressed genes were set as P < 0.0001 and a fold change in expression ≥2.

Gene expression profiling interactive analysis (GEPIA) dataset

GEPIA, which is an interactive web server, was used to assess the RNA sequencing expression, co-expression, and survival data of 8, 587 normal samples and 9, 736 tumors from the Genotype Tissue Expression projects and The Cancer Genome Atlas (TCGA) according to the standard processing pipeline (http://gepia.cancer-pku.cn). GEPIA was used to compare the mRNA expression levels of CXC chemokines between CRC and normal samples. The correlation between CXC chemokines expression and patient prognosis in colon adenocarcinoma (COAD), rectal adenocarcinoma (READ), and colon and rectal adenocarcinoma (COAD + READ) were analyzed in terms of disease-free survival (DFS), overall survival (OS), and the tumor stage. Moreover, the relationships between CXC chemokine expression and EGFR, KRAS and NRAS of COAD, READ, and COAD + READ were evaluated.

Immunohistochemistry (IHC) analysis

All samples were obtained from CRC patients who underwent initial surgical resection in Taizhou People’s Hospital. CRC samples and adjacent noncancerous colorectal tissues were collected for IHC analysis using the following primary detection antibodies: anti-CXCL2 (Bioss, Cat#bs200208R, China), anti-CXCL3 (Bioss, Cat#bs2547R), anti-CXCL8 (PTG, Cat#27095-1-AP, China), anti-CXCL9 (PTG, Cat#22355-1-AP, China), anti-CXCL10 (PTG, Cat#10937-1-AP, China), anti-CXCL11 (PTG, Cat#10707-1-AP, China), and anti-CXCL14 (PTG, Cat#10468-1-AP, China). Immunohistochemistry analysis was used following the manufacturer’s protocol.

cBioPortal for cancer genomics analysis

The CRC datasets (TCGA, Firehose legacy) involving data from 640 cases with pathology reports were selected for further analysis of CXC chemokines using cBioPortal (https://www.cbioportal.org/). The genomic profiles included protein expression z-scores (RPPA), mutations, mRNA expression z-scores (RNA Seq V2 RSEM), and putative copy-number alterations from genomic identification of significant targets in cancer (GISTIC). Moreover, co-expression was assessed according to cBioPortal and using the R package “ggcorrplot”.

GeneMANIA and STRING

GeneMANIA, which is a web-based dataset, was used to analyze the internal associations of gene sets (https://genemania.org/). STRING is a database of predicted protein-protein interactions, including direct (physical) and indirect (functional) associations (https://string-db.org/). Therefore, the interactions of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 were analyzed at the gene and protein expression levels by GeneMANIA and STRING.

Database for annotation, visualization and integrated discovery (DAVID) analysis

As a bioinformatics resource, DAVID (version 6.8) is an open database focusing on gene functional annotations and enrichment pathway analysis (https://david.ncifcrf.gov/home.jsp). In our study, DAVID was used to annotate the cellular component (CC), gene ontology (GO) terms for biological process, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and molecular function categories (MF) of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17.

KEGG analysis

The KEGG (https://www.kegg.jp/) is a database resource for understanding the high-level functions of biological systems. The signaling pathway enrichments of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 were identified and adjusted according to KEGG and Cytoscape.

Tumor IMmune Estimation Resource (TIMER)

TIMER2.0 is a comprehensive resource for the systematic analysis of immune infiltrates across different cancer types from TCGA (http://timer.cistrome.org/) [13, 14]. We identified the association between CXCL2/3/8/9/10/11/14 expression in COAD and READ and the abundances of six immune infiltrates, including CD4⁺ T cells, CD8⁺ T cells, B cells, macrophages, dendritic cells (DCs), and neutrophils. TIMER2.0 generates a heatmap table of the Spearman’s correlations between the expression of CXCL2/3/8/9/10/11/14 and the abundance of immune cell type as well as its subtypes across COAD and READ. Moreover, correlations between genetic markers of tumor-infiltrating immune cells and CXCL2/3/8/9/10/11/14 expression were evaluated using correlation modules (Gene_Corr). According to TIMER2.0, the positive correlation is defined as P < 0.05, rho > 0; the negative correlation is defined as P < 0.05, rho < 0; not significant is defined as P ≥ 0.05.

Transcription levels of CXC chemokines in patients with CRC

The transcription levels of CXC chemokine family members (excluding CXCL15) were assessed in human cancers. As indicated as Table 1 and Figure 1, we used Oncomine to compare the mRNA expression levels of CXC chemokines in CRC with those in normal samples. This analysis showed that the transcription levels of CXCL1/2/3/4/5/6/7/8/9/10/11/16/17 were significantly upregulated in CRC patients. However, the mRNA expression levels of CXCL12/13/14 were significantly downregulated in CRC patients.

Figure 1. The mRNA expression levels of CXC chemokines in different types of cancers (Oncomine). The figure presents the numbers of datasets with statistically significant mRNA over-expression (red) or downregulated expression (blue) of CXC chemokines.

The transcription levels of CXCL1 in colorectal carcinoma were higher than those in colorectal tissues in all 10 datasets. In the Notterman colon dataset, the mRNA expression levels of CXCL1 in colon adenocarcinoma were increased (fold change = 12.328, P = 7.62E-9). In the Skrzypczak colorectal 2 dataset, CXCL2 was overexpressed at the mRNA level in colon adenoma (fold change = 12.924, P = 7.11E-10), colon epithelial carcinoma (fold change = 7.451, P = 2.29E-10), colon epithelial adenoma (fold change = 10.404, P = 3.08E-7), and colon carcinoma (fold change = 8.527, P = 6.09E-9). In the Gaspar colon dataset, CXCL3 was overexpressed in colon epithelial adenoma (fold change = 2.117, P = 5.34E-8) compared with that in colon tissues. In the Skrzypczak colorectal 2 dataset, CXCL4 was overexpressed in colon carcinoma (fold change = 5.228, P = 5.84E-10) and colon epithelial adenoma (fold change = 3.165, P = 8.43E-7). The mRNA levels of CXCL5 in rectosigmoid adenocarcinoma (fold change = 11.693, P = 8.34E-14), colon adenocarcinoma (fold change = 2.739, P = 3.60E-22), colon mucinous adenocarcinoma (fold change = 23.08, P = 9.97E-11), and rectal adenocarcinoma (fold change = 2.538, P = 3.20E-13) were significantly higher in TCGA dataset. In addition, analysis of the Skrzypczak colorectal dataset indicated the mRNA expression levels of CXCL6 were significantly unregulated in patients with colorectal carcinoma (fold change = 4.156, P = 3.39E-10).

A similar trend was found in the Hong and Skrzypczak colorectal datasets, with mRNA expression levels of CXCL7 found to be significantly higher in CRC (fold change = 14.815, P = 4.41E-10 and fold change = 6.016, P = 1.49E-6, respectively). The results from nine datasets indicated that the mRNA expression levels of CXCL8 were significantly upregulated in patients with colorectal carcinoma. In the Skrzypczak colorectal 2 dataset, CXCL9 was overexpressed in colon epithelial carcinoma (fold change = 7.004, P = 2.28E-7) compared with that in non-cancerous colon tissues. Furthermore, in the Skrzypczak colorectal 2 dataset, the mRNA expression levels of CXCL10 were upregulated in colon carcinoma (fold change = 6.673, P = 3.40E-10). In the Skrzypczak colorectal 2 dataset, the mRNA levels of CXCL11 were increased in colon epithelial carcinoma (fold change = 29.164, P = 9.03E-13), colon carcinoma (fold change = 32.137, P = 3.09E-9), colon adenoma (fold change = 4.306, P = 6.60E-7), and colon epithelial adenoma (fold change = 5.543, P = 2.83E-5). In the Hong colorectal dataset, CXCL16 was upregulated in colon carcinoma (fold change = 2.636, P = 2.08E-18). In addition, the results from TCGA colorectal dataset indicated that the mRNA expression levels of CXCL17 were significantly higher in rectosigmoid adenocarcinoma (fold change = 2.664, P = 1.25E-5), rectal mucinous adenocarcinoma (fold change = 3.492, P = 1.60E-5), and colon mucinous adenocarcinoma (fold change = 4.178, P = 2.10E-7).

Expression levels and co-expression analysis of CXC chemokines in patients with CRC

The comparative analysis of the CXC mRNA expression in normal vs. COAD and READ tissues was conducted according to the information from the GEPIA database. As shown as Figure 2, the mRNA expression levels of CXCL1/2/3/4/5/8/9/10/11/13/14/16 were significantly upregulated in COAD and READ tissues vs. normal tissue, while CXCL12 expression was downregulated. In addition, the results shown in Supplementary Table 1 revealed significant associations in COAD + READ between tumor stage and expression of CXCL1 (P = 0.0345), CXCL2 (P = 0.0156), CXCL3 (P = 0.0344), CXCL9 (P = 0.0122), CXCL10 (P = 0.0108), and CXCL11 (P = 0.0227) expression.

Figure 2. The mRNA expression pattern of CXC chemokines (CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17) from GEPIA between CRC tissues (red) and normal tissues (blue). The P-value was set at 0.05, ^* indicate that the results are statistically significant.

We also analyzed potential correlations of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 co-expression with EGFR, KRAS and NRAS according to the mRNA expression levels using the GEPIA database for CRC. The results shown in Table 2 indicated that EGFR expression was significantly correlated with the expression of CXCL3, CXCL12 and CXCL16 in COAD + READ. In addition, a significant association was identified between the expression of KRAS and CXCL16 in COAD + READ. Furthermore, NRAS expression was significantly correlated with the expression of CXCL14, CXCL16 and CXCL17 in COAD + READ.

Association of CXC chemokine expression with the prognosis of CRC

Using GEPIA, we also analyzed the correlation between CXC chemokine expression and the OS and DFS of patients with CRC (Supplementary Table 2). The OS curves shown in Figure 3A demonstrated that high mRNA expression levels of CXCL2 (P = 0.052), CXCL3 (P = 0.015), CXCL8 (P = 0.032), and CXCL14 (P = 0.039) were correlated with longer OS in 181 CRC patients (COAD + READ). As illustrated as Figure 3B, increased mRNA expression levels of CXCL2 (P = 0.042), CXCL3 (P = 0.047), CXCL8 (P = 0.05), and CXCL14 (P = 0.045) were significantly correlated with longer OS in 145 COAD patients. These results indicated that CXCL2/3/8/14 are independent risk factor for OS among CRC patients.

Figure 3. The prognostic value of different expressed CXC chemokines in CRC patients in OS (GEPIA). The OS curve of CXCL2/3/8/14 in (A) COAD + READ and (B) COAD.

The correlations of CXC chemokine expression with the DFS of CRC patients were also investigated (Figure 4.) As shown in in Figure 4A, increased mRNA expression levels of CXCL9 (P = 0.059), CXCL10 (P = 0.0019), and CXCL11 (P = 0.0096) were associated with longer DFS in COAD + READ patients. Moreover, high mRNA levels of CXCL9 (P = 0.026), CXCL10 (P = 0.0011), and CXCL11 (P = 0.0045) were significantly associated with longer DFS in COAD patients (Figure 4B).

Figure 4. The prognostic value of different expressed CXC chemokines in CRC patients in DFS (GEPIA). The DFS curve of CXCL9/10/11 in (A) COAD + READ and (B) COAD.

IHC analysis

As shown in Figure 5, IHC analysis of CXCL2/3/8/9/10/11/14 proteins revealed overexpression of CXCL2/3/8/9/10/11/14 proteins in CRC tissues compared with the adjacent noncancerous tissues.

Figure 5. The expression protein of CXCL2/3/8/9/10/11/14 in CRC tissues and noncancerous tissues (IHC).

Alteration, co-expression and interaction analyses of CXC chemokines at the gene and protein levels in patients with CRC

We investigated the alterations and correlations of CXC chemokine expression at the gene and protein levels using the cBioPortal database for CRC (TCGA, Firehose legacy). CXC chemokine expression was were altered in 131 of 640 (21%) of CRC samples. As shown as in Figure 6A and 6B, we observed alterations in the expression of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 in 1.2%, 3%, 2.3%, 2.5%, 2.7%, 2.1%, 0.4%, 1.6%, 2.7%, 1.6%, 2%, 2.7%, 1.8%, 2.1%, 2.3%, and 1.6% of the queried CRC samples, respectively. We also analyzed the co-expression patterns among the chemokines (CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17) at the mRNA level using the cBioPortal database for CRC. The results shown in Figure 6C indicated a low to moderate correlation in the expression patterns among CXCL4/5/6/7/8/11/13/14/15/16/17; a moderate to high correlation in the expression of CXCL9 with CXCL11, and CXCL9 with CXCL12; and a high correlation in the expression of CXCL1 with CXCL2, CXCL1 with CXCL3, CXCL2 with CXCL3, CXCL9 with CXCL10, and CXCL10 with CXCL11.

Figure 6. CXC chemokines gene expression and mutation analysis in CRC (cBioPortal). (A) Summary of alterations in CXC chemokines (cBioPortal). (B) OncoPrint visual summary of alteration on a query of CXC chemokines (cBioPortal). (C) Correction between different CXC chemokines in CRC (cBioPortal). (D) Gene-gene interaction network (GeneMANIA) and protein-protein interaction network (STRING) among CXC chemokines.

Using GeneMANIA, we identified the relationships among CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 at the gene expression level. As shown as in Figure 6D, several relationships were identified among CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17, including co-localization, co-expression, and genetic interactions. In addition, we analyzed the interactions among CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 at the protein expression level using STRING. In total, 111 edges and 16 nodes were identified in resulting the network.

Cellular functions and pathways of CXC chemokines in patients with CRC

As shown in Supplementary Figure 1, the functions of CXCL1/2/3/4/5/6/7/8/9/10/11/12/13/14/16/17 were predicted by analyzing the GO terms and KEGG pathways using DAVID. The GO enrichment terms were classified into three aspects of function: CC (3 items), MF (5 items), BP (10 items), and KEGG pathway (9 items). The CC analysis demonstrated a significant enrichment of genes related to the extracellular space, external side of plasma membrane, and extracellular region (Supplementary Figure 1A). The MF analysis revealed enrichment in chemokine activity, CXCR3 chemokine receptor binding, heparin binding, CXCR chemokine receptor binding, and receptor binding (Supplementary Figure 1B). In addition, BP was associated with cell-cell signaling, chemotaxis, regulation of cell proliferation, positive regulation of leukocyte chemotaxis, immune response, chemokine-mediated signaling pathway, inflammatory response, response to lipopolysaccharide, cell chemotaxis, and G-protein coupled receptor signaling pathway (Supplementary Figure 1C). As shown as Supplementary Figure 1D, we also found significant enrichment in KEGG pathways, including chemokine signaling pathway (Supplementary Figure 2A), cytokine-cytokine receptor interaction (Supplementary Figure 2B), legionellosis, Salmonella infection, toll-like receptor signaling pathway, NOD-like receptor signaling pathway, rheumatoid arthritis, pertussis, and TNF signaling pathway.

Correlation analysis between CXC chemokine expression and immune infiltration levels in CRC

Tumor-infiltrating lymphocytes are used to analyze the sentinel lymph node status and survival in cancers. According to TIMER2.0, we observed the correlation between the CXCL2/3/8/9/10/11/14 expression and the abundance of immune cell type as well as its subtypes estimated by all six algorithms across COAD and READ. As shown in Figure 7, the heatmap table represented the purity-adjusted partial spearman’s rho values as the degree of correlation. Moreover, a significant association was observed between the CXCL8/9/10/11 expression and immune infiltration in COAD and READ.

Figure 7. The association between the CXC chemokines expression and immune infiltration level of multiple immune cells types estimated by all six algorithms in ad heatmap table across COAD and READ. The red indicates a statistically significant positive association (P ≤ 0.05, rho > 0), and the blue indicates a statistically significant negative association (P ≤ 0.05, rho < 0). White denotes a non-significant result (P > 0.05).

Correlation analysis between CXC chemokine expression and immune markers

To further investigate the association between the expression of CXCL2/3/8/9/10/11/14 and the diverse immune cell infiltration in CRC patients, we analyzed the relationship between CXC chemokines and sets of markers of various immune cells based on TIMER2.0, including tumor-associated macrophages (TAMs), neutrophils, DCs, and T helper 1 (Th1), Th2 and regulatory T (Treg) cells (Figure 8). Specifically, the CXCL 8/9/10/11 expression in COAD and READ was correlated with high levels of infiltration by TAMs, DCs, neutrophils, Th1, Th2 and Tregs. In addition, the results shown in Figure 8 confirmed the role of CXCL2/3/8/9/10/11/14 in the colorectal adenocarcinoma microenvironment.

Figure 8. The association between the CXC chemokines expression and immune markers in ad heatmap table across COAD and READ. The red indicates a statistically significant positive association (P ≤ 0.05, rho > 0), and the blue indicates a statistically significant negative association (P ≤ 0.05, rho < 0). White denotes a non-significant result (P > 0.05).