Characteristics of immune cell inltration landscape in colorectal cancer to Aid in Immunotherapy

The tumor microenvironment (TME) is a complex environment composed of a variety of stromal cells and immune cells that inltrate the tumor space. Recent clinical work has clearly shown such intratumoral immune cell inltration (ICI) to be closely related to colorectal cancer (CRC) patient survival, yet the specic landscape of inltrating immune cells associated with this cancer type remains to be claried. We utilized two computational algorithms to evaluate the ICI status of 712 CRC patients, stratifying these patients into two ICI status-based patterns and assigning ICI scores through the use of principal component analyses. We found that the overall survival (OS) of patients with higher ICI scores was signicantly longer than that of patients with low ICI scores. When ICI scores were combined with the results of tumor mutational burden (TMB) analyses, we determined that CRC patients with both high ICI scores and low TMB exhibited the best survival outcomes. High expression of MARCO in those patients with low ICI scores was correlated with reduced natural killer (NK) and effector T cell inltration and with increased regulatory T cell inltration, suggesting that these factors may be linked to poor patient prognosis. These results suggest ICI scores to be a valuable biomarker for the prognostic evaluation of CRC patients. Future efforts to analyze the ICI patterns of larger patient sample cohorts will help to extend these analyses, offering new insights into the role of the TME in cancer progression while highlighting novel immunotherapeutic approaches to treating this cancer type.


Introduction
Colorectal cancer (CRC) is among the most common and deadliest cancers in developed nations, with an estimated 400,000 CRC diagnoses and 212,000 deaths globally each year [1]. Immunotherapeutic treatments serve to augment the ability of natural host defense mechanisms to recognize and remove tumor cells, and such regimens have proven effective in synergistically enhancing the survival of patients with a range of tumor types [2][3][4]. However, immunotherapies generally only bene t a small percentage of patients, and there is thus a clear need to better understand which therapeutic markers can be analyzed to gauge CRC patient responses to immunotherapeutic treatment [2][3][4][5].
The tumor microenvironment (TME) is a complex setting that in uences cancer development and progression [6]. The TME contains a diverse array of stromal cells, lymphatic structures, and blood vessels, all of which can in uence and be in uenced by the oncogenic mutations driving the onset and evolution of a given cancer. The in ltration of various immune cell types into the TME has been shown to be an effective prognostic biomarker for tumor invasion status, offering insights into cancer grade, stage, and metastasis [7][8]. Tumor-associated macrophages (TAMs), for example, can produce cytokines including interleukin-10 (IL-10) and transforming growth factor-B (TGF-B) that can suppress immune cell activation, thereby promoting tumor proliferation and contributing to poorer patient prognosis [9][10][11]. In contrast, higher levels of tumor-in ltrating lymphocytes (TILs) including both CD4 + and CD8 + T cells are related to improved patient survival and immunotherapeutic responsiveness [12].
Immune checkpoint blockade therapies rely on the activation of extant TILs, enhancing their ability to recognize tumor cells and to thereby engage appropriate immunotherapeutic responses [13][14]. Merely identifying these TILs, however, is not su cient to reliably characterize the complexities of the TME, as many patients with high TIL levels nonetheless exhibit resistance to immunotherapy treatment [15][16].
Such effects may be in part attributable to TAM-derived immunosuppressive cytokine production, or to stromal in ltration of the tumor, as this can interfere with TIL accumulation [14][15]. These prior results suggest that dynamic interactions between cell types are more important in the TME than are any individual cell population from a prognostic perspective. Prior studies have not performed a high-level analysis of the prognostic relevance of different immune cell in ltration (ICI) populations in CRC.
Herein, we employed the CIBERSORT and ESTIMATE algorithms to characterize CRC patient tumor gene expression pro les in an effort to comprehensively characterize intratumor immune cell landscape characteristics within these tumors [17][18]. We then strati ed CRC patient tumors into two subtypes based upon the observed degree of ICI, and we established ICI scores that were successfully used to gauge patient prognosis and which may offer value for the design of novel immunotherapeutic regimens.
Tumor-Infiltrating Immune Cell Consensus Clustering The CIBERSORT R package 21 was used to estimate tumor infiltration by different immune cell populations based upon 1000 permutations and LM22 markers, with both immune cell and matrix contents being evaluated to yield specific immune and matrix scores. CRC patient samples were then clustered according to sample-specific ICI patterns. A PAM-based approach relying upon Euclidean and Ward's linkages was utilized for unsupervised clustering with the ConsensuClusterPlus R package, and clustering was conducted 1000 times to ensure stable classification.
Identification of ICI-related gene expression patterns ICI profiles were used to group CRC patients into different ICI clusters, after which differentially expressed genes (DEGs) between these clusters were compared with the following cutoff criteria using the limma R package: adjusted p < 0.05; absolute fold-change > 1.4.

Dimensional Reduction and ICI Score Calculation
An unsupervised clustering approach was initially used to classify patients according to DEG expression values, with DEGs that were positively and negatively correlated with key clustering-related features being respectively referred to as ICI gene features A and B. Dimensional reduction was then performed for these ICI-related genetic signatures using a Boruta algorithm, and principal component 1 in the resultant principal component analysis (PCA) was then extracted as a signature score. A method similar to that previously used to compute a gene expression grade index [42] was then employed to calculate individual patient ICI scores as follows: Somatic Alteration Data Compilation Mutational data corresponding to CRC patients in the TCGA database was downloaded and used to estimate tumor mutational burden (TMB). CRC driver genes were identified with the maftool R package, and somatic mutations therein were evaluated for those with high and low ICI scores. The top 20 most frequently mutated driver genes were then retained for further analysis.
Statistical Analyses GraphPad Prism 7.0 and SPSS 21.0 (IBM, NY, USA) were used for all statistical testing. Data between groups were compared via the Wilcoxon test. Patient classification into two ICI subtypes was achieved with the X-tile software via an iterative approach aimed at reducing computational batch effects [43].
Kaplan-Meier plotter was employed to generate all survival curves, which were analyzed with log-rank tests. Chi-squared tests were employed to assess correlations between ICI score subgroups and somatic mutation frequencies, while correlations were assessed via Spearman's correlation analyses. A two-tailed p < 0.05 was the significance threshold for this study.

Assessment of the CRC-related TME Immune Cell In ltration Landscape
We began by employing the ESTIMATE and CIBERSORT algorithms to gauge the relative enrichment of different immune cell populations in samples of tumors from 712 CRC patients (Fig. 1A). Data were obtained from the TCGA and GEO databases (Accession number: GSE17538), and analyzed tumor cells exhibited matched immune cell in ltration (ICI) pro les. The ConsesusClusterPlus R package was then used to classify these CRC patients into distinct ICI-related subtypes via an unsupervised clustering approach (Fig. 1B).
Signi cant differences in overall survival (OS) were evident between patients with the two independent ICI-related CRC subtypes identi ed in this study (p = 0.047; Fig. 1C). To better understand the biological basis for these prognostic relationships, we further assessed the immune cell makeup of the TME in these patient clusters (Clusters A and B). Patients in cluster A had a better prognosis, and exhibited tumors predicted to contain higher levels of naive B cells, M1 macrophages, memory CD4 T cells, CD8 T cells, NK cells, monocytes, and plasma cells. Patients in ICI cluster B had a median survival of 7 years, and exhibited tumors predicted to contain higher levels of regulatory T cells (Tregs), M0 macrophages, and activated mast cells ( Fig. 1D-E). Correlation coe cient heatmaps were also prepared to visualize the immune cell in ltration landscape of the TME in these CRC patient cohorts (Fig. 1F).
We further assessed the relative expression of the key immune-related gene MARCO in these two ICI patient clusters, revealings signi cantly increased MARCO expression among patients in cluster B relative to those in cluster A. Differences between immune cell types and MARCO expression levels in these patient cohorts were then evaluated (Fig. 1G).

Immune Gene-Related Subtype Analyses
To better understand the biological basis for the ICI subtypes of CRC patients identi ed above, we conducted a differential transcriptomic analysis of these patient cohorts with the R limma package. Through unsupervised clustering of 34 identi ed DEGs, these patients were strati ed into two genomic clusters (clusters A and B) ( Fig. 2A). In total, 23 of these DEGs were positively correlated with this gene cluster and were denoted as 'ICI gene signature A', whereas the remaining DEGs were referred to as 'ICI gene signature B'. The Boruta algorithm was then used to conduct a dimensionality reduction for these two ICI gene signatures in order to reduce computational noise associated with these analyses [19]. The clusterPro ler R package was then used to prepare a heatmap separating the transcriptomic pro les of these 34 DEGs across genomic clusters (Fig. 2B) as described previously [20], with signi cantly enriched biological processes associated with these genes being shown in Figs. 2C and 2D.
The prognostic relevance of these ICI gene clusters was next evaluated using the Kaplan-Meier plotter tool, revealing that patients in cluster A had a better prognosis relative to patients in cluster B (p = 0.002; Fig. 2E). Signi cant differences in MARCO gene expression were also detected between these two genomic clusters (p < 0.001; Fig. 2F), with the expression of this key immunological gene being signi cantly higher in samples from patients in cluster B relative to those from patients in cluster A.

Evaluating Correlations between ICI Scores and Somatic Variation
A PCA approach was next employed to compute aggregate ICI scores for the two ICI gene signatures detailed above in an effort to quantitatively evaluate the ICI landscape of CRC. The resultant scores (ISA and ISB for genes in signatures A and B, respectively) were then calculated for each patient by summing the relevant individual scores, yielding a prognostic signature which was de ned as an ICI score. The patients in our study cohort were then separated into two groups based upon whether they had high or low ICI scores, with the cutoff value used to differentiate these two groups being established with the Xtile software. Patient distributions in these two gene clusters are shown in Fig. 3A. We then evaluated immune activity in each patient group before assessing the prognostic relevance of our ICI scores, revealing that most immune checkpoint-and immune activity-related genes were signi cantly overexpressed in samples from the low ICI group with the exception of PIGR, IGLJ3, CLCA1, ITLN1, IGHM, PLA2G2A, CLCA4, ZG16, UGT2B17, REG4, and SPINK4 (Fig. 3B). TAMs expressing MARCO can interfere with the activation of NK and T cells, suppressing their cytokine production, proliferation, and cytotoxicity, while simultaneously aiding Treg proliferation (21). In line with such a model, we found that high MARCO expression in individuals with low ICI score subtypes was correlated with higher Treg in ltration and reduced NK and effector T cell levels, potentially explaining the poorer prognosis for these patients. A gene set enrichment analysis (GSEA) further showed the regulation and vascular signaling pathways to be signi cantly enriched among patients in the low ICI score group, while the butanoate and retinol signaling pathways were signi cantly enriched among those in the high ICI group (Figs. 3C).
Many prior studies have shown that tumors bearing a high mutational burden are more likely to exhibit increased CD8 + T cell in ltration, aiding in their elimination. As such, TMB may predict a given cancer patient's responsiveness to immunotherapy treatment [22][23]. We therefore evaluated potential relationships between TMB and ICI scores in CRC pateints, evaluating somatic variant distributions in CRC driver genes in the low and high ICI score patient subgroups using maftools [24]. We then retained the top 20 driver genes with the highest mutational frequency for further analysis (Fig. 3D-3E).

Assessment of the Prognostic Value of ICI Scores
Finally, we conducted an in-depth analysis of the prognostic relevance of the ICI scores computed above using the Kaplan-Meier plotter tool. The OS of patients in the high ICI score group was signi cantly better than that of those in the low ICI score group (p = 0.004; Fig. 4A). We then evaluated the synergistic value of ICI and TMB scores as predictors of CRC patient survival in a strati ed survival analysis which found TMB status to be unrelated to ICI score-based prognostic analyses such that higher ICI scores were linked to better patient OS regardless of TMB status (High TMB + High ICI score (HH) vs. High TMB + Low ICI score (HL), Low TMB + high ICI score (LH) vs. Low TMB + Low ICI score (LL), p = 0.006; Fig. 4B). Together these data suggested that ICI scores serve as prognostic biomarkers that are independent of TMB status. The mortality rate for patients with high ICI scores was also decreased relative to that of patients with low ICI scores (Fig. 4C), and the ICI scores of surviving patients were signi cantly higher than those of nonsurviving patients (p = 0.027; Fig. 4D) We additionally evaluated the relationship between the prognostic value of ICI scores and other clinically relevant parameters in CRC patients. This approach revealed high ICI scores to be associated with signi cantly better survival among individuals ≤ 65-years-old (p = 0.024; Fig. 4E) and those with advanced CRC (p = 0.01; Fig. 4F),. This also remained true regardless of patient gender (male p = 0.04; female p = 0.026; Fig. 4G-4H).

Discussion
Immunotherapeutic treatment represents a highly e cacious approach to controlling the growth of many tumors, and can improve the quality of life for advanced CRC patients. However, many patients fail to respond to immunotherapies, and the Association for Cancer Immunotherapy has emphasized the importance of identifying those patients most likely to bene t from these therapeutic regimens [25].
Herein, we developed a quantitative approach to evaluating the immunological TME associated with CRC patient tumors. The results of this study indicate that ICI scores represent a valuable prognostic biomarker that can be used to predict outcomes in these patients, and we found that higher ICI scores as well as increased intratumoral CD4+, CD8+, plasma cell, and macrophage in ltration were linkedto better patient outcomes, in line with prior research [26][27]. This underscores the potential for already extant immune responses to achieve anti-tumor effects and to thereby affect the way cancer patients respond to immunotherapy. CRC tumors are thought to exhibit among the highest levels of immune cell in ltration on average [28]. We thus hypothesized that comprehensively characterizing ICI pro les and related patterns of gene expression would represent a novel approach to developing patient-speci c evaluation and treatment strategies. We began by characterizing the TME associated with CRC tumors at a molecular level, enabling us to identify immune-related genes based upon ICI gene clusters. ICI gene cluster B was associated with lower immune scores, matrix scores, and with decreased immune cell in ltration consistent with what is often referred to as a 'cold' immune phenotype. In contrast, ICI gene cluster A was associated with higher immune scores and in ammatory cell in ltration, and these patients had a more favorable prognosis as well as increased CD8 + T cell, activated CD4 + T cell, and plasma cell in ltration [29][30]. Indeed, many prior studies have highlighted the effects of the TME on cancer patient OS [31]. The anti-tumor immune responses observed for patients in ICI gene cluster A suggest that they are likely to attain more bene ts from immunotherapeutic treatment, and these clusters may thus offer value for the development of more e cacious immunotherapies. However, given the signi cant heterogeneity associated with the immunological TME in a given patient, it is vital that ICI patterns be evaluated on a patient-by-patient basis. Consistent with such an approach, one research group has reported to use of an individually-tailored tumor-speci c biomarker model that can more reliably gauge breast cancer patient prognosis [32].
Using the Boruta algorithm, we established an ICI score that we were then able to analyze as a biomarker of CRC patient prognosis. GSEA analyses of these genes composing these ICI-related patterns revealed the regulation and vascular signaling pathways to be signi cantly enriched among samples from those with low ICI scores, while butanoate and retinol signaling pathways were enriched among those with high ICI scores. Recent wor has clari ed the relationship between genetic mutations and patient sensitivity to immunotherapy [33][34]. By leveraging this fact and evaluating key immunological parameters including consensus ICI scores, tumor driver mutation analyses [35], TMB measurements [36], LOH HLA (loss of heterozygosity at the HLA locus) assessments [37][38], PD-L1 expression analyses [39], and the detection of key immune gene expression-related signatures [40][41], it is possible to more reliably classify particular cancers. Integrated ICI scores offered new insights regarding differences in variant frequencies in many genes when comparing samples from the low and high ICI score groups, with some of these genes being directly linked to therapeutic sensitivity or resistance. When we performed a strati ed analysis, we determined the prognostic value of ICI scores to be independent of TMB status for CRC patients, suggesting that these two metrics offer distinct insights into the immunobiology of patient tumors, enabling the more robust assessment of patient outcomes. After employing two different algorithms to classify 712 CRC patients based upon their ICI pro les, we identi ed two ICI patterns and were able to employ a PCA approach to derive ICI scores therefrom. Higher ICI scores were associated with signi cantly higher patient OS relative to lower ICI scores, and this remained true in patients with advanced disease and in those ≤ 65-years-old. High MARCO expression in those with low ICI score subtype disease was correlated with increased Treg in ltration and reduced levels of NK and effector T cells, which may be linked to poorer patient prognosis.
Lastly, we explored ICI scores and associated patient characteristics, and found that simultaneous analyses of ICI scores and TMB status may improve our ability to reliably gauge CRC patient prognosis, offering a means of potentially identifying patients at a higher risk of tumor recurrence, although further research is required to test this possibility. ICI scores nonetheless represent a powerful tool for estimating tumor-speci c immune tness, and can be used to predict which patients are most likely to bene t from immunotherapy, thereby helping to improve CRC patient survival.

Conclusions
Overall, we found that ICI scores can be reliably leveraged to assess CRC patient prognosis, with higher ICI scores being associated with prolonged OS. Strati ed analyses revealed the prognostic value of ICI scores to be independent of TMB status for this cancer type, and this, together with observed GSEA and predictive outcomes, suggests that TMB and ICI status are independent tumor characteristics that may predict patient responses to immunotherapeutic treatment. We further found that that high MARCO expression in CRC patients with low ICI scores was correlated with reduced NK and effector T cell in ltration and increased Treg in ltration, potentially thereby contributing to poorer patient outcomes. Future studies of larger CRC patient cohorts will continue to offer new insights into the role of the TME as a driver of cancer progression, and may aid in the design of novel immunotherapeutic approaches to CRC treatment. All data generated or analyzed during this study are included in this published article.

Con icts of Interest
The authors declare there is no con ict of interest.     Evaluation of the prognostic value of ICI scores (A) The survival outcomes of CRC patients in the low and high ICI score groups in the TCGA dataset were evaluated via Kaplan-Meier curves and log-rank tests (p=0.004). (B) TCGA-CRC patients were strati ed according to TMB and ICI scores, after which differences in OS outcomes among cohorts were analyzed (p=0.006, log-rank test). (C) Survival outcomes and clinical responsiveness in those with low and high ICI scores. (D) Differences in ICI scores between surviving and non-surviving patients (Wilcoxon test, p = 0.027). (E) Survival outcomes for patients ≤ 65years-old with low and high ICI scores were evaluated (p=0.024; log-rank test). (F) Survival outcomes for CRC patients with advanced disease with low and high ICI scores were evaluated (p=0.01; log-rank test).
(G and H) The relationship between survival and ICI scores in male or female patients was evaluated (p=0.04; log-rank test).