Three independent datasets were included in this study: one dataset from the TCGA (n = 451) and two single-cell RNA sequencing (scRNA-seq) datasets from GEO (GSE146771 and GSE132465). In addition, the patient cohort YJSHC (the Yijishan Hospital cohort, n = 108) was also involved. The characteristics of patients in the TCGA were downloaded from http://www.cbioportal.org on 18 April 2020. The clinical information of enrolled patients in the TCGA and YJSHC in this study are provided in Supplementary Table 1 and Supplementary Table 2. The YJSHC consisted of patients with COAD from Yijishan Hospital affiliated with Wannan Medical College (Wuhu, China), who underwent surgery between August 2011 and December 2014. All of the research was carried out in accordance with the Declaration of Helsinki and approved by the ethics committee of Wannan Medical College, Yijishan Hospital. Additionally, informed consent was obtained from every patient.

We collected formalin-fixed, paraffin-embedded surgical specimens, which were used for tissue microarray construction and subsequent immunohistochemistry studies. First, tissue microarray sections were rehydrated, and treated with hydrogen peroxide and then heat-mediated antigen retrieval by microwaving with sodium citrate. Next, the slides were incubated with the indicated antibodies obtained from Abcam (Cambridgeshire, United Kingdom) in a humidified box at 4°C overnight. The details for the four antibodies used in this study are provided in Supplementary Table 3. Next, the color reaction was realized by a DAB substrate kit and nucleus counterstaining was performed with hematoxylin. The level of protein expression was assessed based on the intensity of staining and the extent of staining at 200 × under a microscope. To determine the extent of staining, a quantity score (0–4) denoted 0, 0%; 1, 1–10%; 2, 11–50%; 3, 51–80% and 4, 81–100% of positive cells, was used. The staining intensity was divided into three grades: weak, moderate, and strong staining. It should be noted that the corresponding intensity scores ranged from 1 to 3. The final IHC score was calculated by multiplying the quantity and intensity scores.

The relative abundance of the tumor infiltrating lymphocytes (TILs) in COAD tissues with different CXCL11 mRNA expression statuses was calculated by CIBERSORT¹ (Newman et al., 2015) and TISIDB² (Ru et al., 2019). For each cancer type, GSVA (gene set variation analysis) was adopted to infer the relative proportion of various TILs based on the immune-related gene expression profile of 28 TIL types from Charoentong’s study (Charoentong et al., 2017). The correlation between CXCL11 and TILs was measured by Spearman’s test. In this study, we also investigated the differential expression of CXCL11 between tumor and adjacent normal tissues across all TCGA tumors in TIMER³ (Li et al., 2017).

The R package “limma” was used for differential expression analysis. Differentially expressed genes between CXCL11-high and CXCL11-low tumor samples were defined as the up-regulated or down-regulated genes with a fold change ≥ 2 and false discovery rate (FDR) P < 0.05.

The significantly enriched pathways in CXCL11-high tumor samples were identified by gene set enrichment analysis (GSEA) performed by “MSigDB” (Molecular Signatures Database, version 6.0). Additionally, the enrichment score was assessed by determining the expression of a gene set. When the majority of a gene set had elevated expression, we considered the gene set to be a “enriched” gene set.

External single-cell transcriptomic data of COAD patients were retrieved from GEO (GSE146771 and GSE132465). Single-cell data analysis was performed using scanpy v.1.4⁴. Since single-cell technologies currently capture only a portion of the transcripts in any cell, we generated “pseudobulks” as technical replicates by summing the raw UMI counts for each gene from the same patient.

In this study, R and corresponding R packages were utilized for statistical analysis. We used the Cutoff Finder platform⁵ to automatically determine the cutoff points (Budczies et al., 2012). The results were all displayed as the mean ± SD. Statistical significance was defined by a two-tailed P < 0.05. Student’s t-test was adopted for the analysis of continuous variables. Spearman’s correlation was performed to analyze correlations. We performed survival analysis with Kaplan-Meier method, and survival was compared among multiple groups with the log-rank test.

First, we evaluated the distribution of CXCL11 in tumor tissue across all TCGA cancers in TIMER and found that CXCL11 expression was significantly upregulated in the majority of tumors, including COAD and READ (Figure 1A). Next, TISIDB was utilized to examine the association between CXCL11 and clinical outcome across all TCGA tumors, and we found that COAD patients but not READ patients with high mRNA expression of CXCL11 had a better prognosis, instead of READ (Figure 1B). Therefore, we focused on the prognostic significance of CXCL11 in COAD, and we downloaded a dataset from the TCGA (n = 451). The results revealed that high CXCL11 mRNA expression was associated with prolonged survival (P = 0.0053, Figure 1C).

To understand the underlying mechanism of CXCL11 in COAD, we investigated the associations of CXCL11 with the tumor immune microenvironment. TISIDB was employed to investigate which kinds of TILs might be regulated by CXCL11 across all TCGA tumors (Supplementary Figure 1). We found that CXCL11 was positively correlated with Act CD8 (activated CD8⁺ T cells; r = 0.55, P < 0.001), NKT (natural killer T cells; r = 0.438, P < 0.001), and Act DC (activated dendritic cells; r = 0.488, P < 0.001) in COAD, as shown in Figures 2A–C.

Furthermore, the differential status of the tumor immune microenvironment in different groups of COAD patients determined by CXCL11 mRNA expression patterns was validated by CIBERSORT (Figure 2D). Consistently, the high CXCL11 mRNA expression group had a higher proportion of antitumor immune cells, such as: CD8T (CD8⁺ T cells, P < 0.001), NKa (activated natural killer cells, P < 0.05). Additionally, the high CXCL11 mRNA expression group had a lower proportion of protumor immune cells, such as: M0 (M0 macrophages, P < 0.05), NKr (resting natural killer cells, P < 0.001), and Mono (monocytes, P < 0.05). The results suggested that CXCL11 was associated with antitumor immunity in COAD, which partially explained the association of CXCL11 with a better prognosis.

As CXCL11 was correlated with antitumor immunity, we aimed to identify the differentially expressed genes and enriched biological pathways mediated by CXCL11 in the tumor immune microenvironment. First, we performed GSEA and identified that the biological pathways enriched in CXCL11-high tumor samples were the chemokine signaling pathway (NES = 2.3, P < 0.001), cytokine-cytokine receptor interaction (NES = 2.44, P < 0.001), antigen processing and presentation (NES = 2.36, P < 0.001), T cell receptor signaling pathway (NES = 2.07, P < 0.001), natural killer cell mediated cytotoxicity (NES = 2.31, p < 0.001), intestinal immune network for IgA production (NES = 2.25, P < 0.001), interferon alpha response (NES = 2.57, P < 0.001), and interferon gamma response (NES = 2.74, P < 0.001), which are known to contribute to the antitumor immunity, as shown in Figures 2E–L.

We then studied genetic alterations and found that the differentially expressed genes in tumors in the high CXCL11 mRNA expression group included various immune-activated genes. A volcano plot is shown in Figure 3A, and we observed that the CXCL11- high group was enriched in CXCL9, CXCL10, CCL4, CCL5, IFNG, CD8A, and PD-L1. As shown in Figures 3B,C, CXCL11 was positively correlated with a gene set associated with DC-, NK cell and T cell recruiting CCL4 (r = 0.543, P < 0.001), CCL5 (r = 0.59, P < 0.001), CXCL9 (r = 0.717, P < 0.001), and CXCL10 (r = 0.821, P < 0.001). In some ways, CXCL11 was found to contribute to the antitumor immune microenvironment in COAD.

Cytotoxic genes including IFNG, GZMA, GZMB, GZMK, GZMM, and PRF1 were positively correlated with CXCL11 mRNA expression (Spearman’s q = 0.49, 0.55, 0.22, 0.37, 0.21, and 0.37; all P-values are < 0.001), as shown in Figure 3D. Interestingly, several immunosuppressive molecules including PDCD1, PD-L1, and CTLA4 were positively correlated with CXCL11 mRNA expression (Spearman’s q = 0.35, 0.56, and 0.44; P < 0.001, P < 0.001, and P < 0.001, respectively), as shown in Figure 3E. Furthermore, we found upregulation of CXCL11 expression in tumor tissue compared with adjacent normal tissue (Figure 3F) in the COAD single-cell RNA-seq datasets GSE146771, and CXCL11 expression was positively correlated with PD-L1 (r = 0.66, P < 0.001; Figure 3G). In addition, the UMAP plots showing the expression levels of CXCL11 and PD-L1 in different clusters of GSE146771 were in Supplementary Figures 2A–C. We also found the upregulation of CXCL11 expression in tumor tissue compared with adjacent normal tissue (Figure 3H) in the COAD single-cell RNA-seq datasets GSE132465, and CXCL11 expression was positively correlated with PD-L1 (r = 0.78, P < 0.001; Figure 3I).

To further validate the effect of CXCL11, we enrolled the COAD patient cohort from Yijishan Hospital (YJSHC). First, we verified the upregulated protein expression of CXCL11 in tumor tissue compared to adjacent normal tissue (Figures 4A,B) by IHC staining. Next, the patients were divided into two groups based on their CXCL11 expression and the respective clinical outcomes were analyzed. We found that the patients within the YJSHC with a high abundance of intratumoral CXCL11⁺ cells had better OS (overall survival, P = 0.001) within the YJSHC, as shown in Figure 4C. Additionally, multiple Cox regression analysis was used to verify that the prognostic value of intratumoral CXCL11⁺ cell infiltration kept significant after adjusting for confounders (HR = 0.35; 95% CI 0.17–0.72; P = 0.004; Figure 4D). It was confirmed that CXCL11 could be an independent prognostic biomarker.

Further analysis of the YJSHC was conducted to validate the immunomodulatory effect of CXCL11⁺ cell infiltration. As shown in serial sections (Figure 5A), tumor with a high abundance of CXCL11⁺ cell infiltration tended to have a high abundance of intratumoral CD8⁺ T cell and CD56⁺ NK cell infiltration, which was associated with antitumor immunity. In addition, a high abundance of CXCL11⁺ cell infiltration was associated with a high abundance of intratumoral PD-L1⁺ cells. Accordingly, the expression of PD-L1, CD8A, and CD56 was higher in the CXCL11-high group than in the CXCL11-low group (Figures 5B–D), and the positive correlations between the expression of CXCL11 and PD-L1, CD8A, or CD56 (r = 0.62, P < 0.001; r = 0.34, P < 0.001; r = 0.32, P < 0.001, respectively; Figures 5E–G) were significant within the YJSHC.