Novel Hypoxia-Associated Gene Signature Depicts Tumor Immune Microenvironment and Predicts Prognosis of Colon Cancer Patients

Abstract

Hypoxia, a typical hallmark of numerous tumors, indicates poor infiltration of antitumor lymphocytes, as well as facilitates the development, progression, and drug resistance of malignant cells. Here, the present research was performed to identify novel hypoxia-related molecular markers and their correlation to the tumor immune microenvironment (TIME) in colon cancer. The expression of hypoxia-related gene signature was extracted from The Cancer Genome Atlas (TCGA) COAD cohort. Based on this signature, a risk score model was constructed using the Lasso regression model. Its discrimination ability and stability were validated in another independent cohort (GSE17536) from Gene Expression Omnibus (GEO) database. Moreover, molecular biology experiments (quantitative real-time PCR and multiple immunohistochemistry) were performed to validate the results of bioinformatics analyses. Three hub genes, including PPFIA4, SERPINE1, and STC2, were chosen to build the risk score model. All of these genes were increasingly expressed in the hypoxia subgroup (HS). Compared with the normoxia subgroup (NS), HS had worse pathological features (T, N, M, and stage) and overall survival (OS), more expression of immune checkpoint molecules, poorer infiltration of some pro-inflammation immune cells (CD4⁺ T cells and CD8⁺ T cells), and enriched infiltration of M0/M2 macrophages. After the risk model was proven to be valuable and stable, a nomogram was built based on this model and some clinicopathological factors. Moreover, it had been identified that three hub genes were all increasingly expressed in hypoxic conditions by quantitative real-time PCR (qPCR). The results of multiple immunohistochemistry (mIHC) also showed that higher expression of hub genes was associated with poorer infiltration of pro-inflammation immune cells (CD8⁺ T cells and M1 macrophages) and richer infiltration of anti-inflammation immune cells (Treg cells and M2 macrophages). In conclusion, the present study uncovered the relations among hypoxia, TIME, and clinicopathological features of colon cancer. It might provide new insight and a potential therapeutic target for immunotherapy.

Introduction

Colon cancer (CC) is one of the most common malignancies worldwide and responsible for more than 0.5 million deaths in 2020 (Sung et al., 2021). Compared with 2018, the numbers of new cases and deaths of CC have rapidly increased (Bray et al., 2018; Sung et al., 2021). Despite the advance in medical technology, patients with CC still have a relatively high mortality rate, 13.1% in the transitioning and 4.7% in the transitioned countries (Sung et al., 2021). Owing to the heavily medical and financial burden caused by CC, it is urgent to develop novel methods to improve the diagnostic and therapeutic efficiency for these patients.

Recently, immunotherapy has attached much attention from the public for its promising therapeutic efficiency. It is well established that the killing effect induced by immunotherapy relies on some immune cells that are recognized as tumor suppressors (Hiam-Galvez et al., 2021). While in the complex tumor microenvironment (TME), besides these tumor suppressors, there are many other components that promote the development and progression of cancer (Anderson and Simon, 2020). These promoters and the tumor itself deprive of the oxygen and nutrient and subsequently produce a hypoxic and acidic TME, which significantly restrain the function of those antitumor immune cells (Kaymak et al., 2021). Previous studies have well established that hypoxia is a hallmark of tumor growth, survival, and metastasis of CC and confers to resistance to immunotherapy (Hsu et al., 2020; Singhal et al., 2021). Therefore, the establishment of a hypoxia-related gene signature may help to comprehend the immunogenomic profile of CC and provide a useful prognostic tool for CC patients.

Here, based on The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases, we developed a hypoxia-related gene signature to depict the tumor immune microenvironment (TIME) and predict the overall survival (OS) of CC patients. Moreover, we also did quantitative real-time PCR (qPCR) and multiple immunohistochemistry (mIHC) to verify the results of bioinformatics analyses.

Materials and Methods

Training and Validation Cohort

The FPKM RNA-seq data (398 tumors and 39 normal tissue samples) and related clinical information of 385 CC patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/) using the GDC API tools on 7 July 2021. Then, 282 patients with complete overall survival (OS) data were included in the training set. Meanwhile, the GSE17536 cohort was applied as the independent validation set. It included 177 colon tumor tissue samples and was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse17536). The FPKM RNA-seq data from TCGA were transformed to log₂(x+1). Then, they were normalized to eliminate the batch effect using limma (version 3.48.1) and sva (version 3.40.0) packages.

Gene Signature List

200 hypoxia-related genes were accessed from the HALLMARK_HYPOXIA gene set of Gene set Enrichment Analysis (GSEA) database (http://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_HYPOXIA.html). The complete gene list was contained in Supplementary Table S1.

The immune-related genes were obtained from the Tracking Tumor Immunophenotype database (http://biocc.hrbmu.edu.cn/TIP/index.jsp) (Xu et al., 2018). This gene list contained negative regulatory, positive regulatory, T cell, CD8⁺ T cell, CD4⁺ T cell, dendritic cell, eosinophil, macrophage, monocyte, neutrophil, nature kill (NK) cell, Th1 cell, Th17 cell, Th 2 cell, Th22 cell, and Treg cell-correlated genes.

Procedure of Developing Risk Score Model

First, the fold change (FC) of the 200 hypoxia-related genes between tumor and normal tissue samples was calculated using the limma package. Genes with log₂|FC|>1 & adjusted p-value < 0.05 were identified as the differentially expressed genes (DEGs). Meanwhile, the statistically prognostic genes were identified using univariate Cox regression analysis. Then, the DEGs (Supplementary Table S2) and prognostic genes (Supplementary Table S3) were intersected to identify hub genes.

Based on the least absolute shrinkage and selection operator (LASSO) regression analysis, the formula of the risk score model was built as follows:The index represents a significantly prognostic gene of the Lasso regression analysis and stands for the beta coefficients of these genes.

Differences in RNA Expression and Clinical Characteristics Between Subgroups

To compare the differences between subgroups in RNA expression/gene function/clinical characteristics, we used Rtsne (version 0.15) and pheatmap (version 1.0.12) packages/GSEA analysis/stats (version 4.1.0) package.

Evaluation of Tumor Immune Microenvironment and Drug Response

To investigate the association between tumor immune microenvironment (TIME) and the risk score model, we used CIBERSORT and microenvironment cell populations-counter (MCP-counter) to estimate the infiltration of different immune cells (Newman et al., 2019; Becht et al., 2016). While for single-gene analysis, we used the TIMER webserver to evaluate the relation between six different types of immune cells and the target gene (https://cistrome.shinyapps.io/timer/) (Li et al., 2016; Li et al., 2017).

To evaluate the response of immune checkpoint blockade (ICB), we used the website tool ImmuCellAI (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/), which was based on ssGSEA analysis (Miao et al., 2020; Miao et al., 2021). But the pRRophetic package (version 0.5) was applied to compare the differences in drug response for cytotoxic and targeted medicine (Geeleher et al., 2014).

Model Visualization, Validation, and Comparison

To visualize the risk score model, we used the rms package (version 6.1-0) to create a nomogram that could predict the 1-, 3-, and 5-year OS of CC patients. It contained clinical factors (age, pathological M, and stage) and the risk score.

The discrimination ability of the risk score model was assessed using operating characteristic curve (ROC) analysis. Based on a series of different binary classification methods (critical or cutoff value), it could calculate the true positive (TP) and false positive (FP). The curve was drawn with TP or sensitivity as the ordinate, and with FP or 1-specificity as the abscissa. The area under the curve (AUC) was used for quantitative analysis in ROC analysis. Generally (Sung et al., 2021) AUC between 0.5 and 0.7 would be considered lower accuracy (Bray et al., 2018); AUC between 0.7 and 0.9 would be considered to be valuable (Hiam-Galvez et al., 2021); AUC above 0.9 would be considered high accuracy; however, AUC = 0.5 meant that the model had no diagnostic value. In this study, ROC analysis was performed using the timeROC package (version 0.4).

The calibration plot was used to assess the 1-, 3-, and 5-year OS, based on the Cox proportional hazard model. Meanwhile, the Kaplan-Meier (K-M) survival analysis was performed according to different subgroups (NS vs. HS) in both training and validation sets. The calibration and K-M analysis were performed using the rms and survival (version 3.2-7) packages.

In this study, Decisive Curve Analysis (DCA) was used to estimate prediction ability among different models. The abscissa of the DCA plot was threshold probability, and the ordinate was the net benefit (NB) after therapeutic advantages minus the disadvantage. In general, the farther the curve was from the extreme curves the better its prediction power represented. The DCA analysis was performed using the ggDCA (version 1.2).

Cell Culture and Culture Conditions

The human colon epithelial cell line (FHC) and colon cancer cell lines (HCT-8, RKO, SW480, and SW620) used in vitro experiment was purchased from the Cell Bank of the Chinese Academy of Science and authenticated by the supplier. These cell lines were grown in DMEM (Gibco) supplemented with 10% fetal calf serum (FCS). We have identified the source of cell lines by STR profiling. Meanwhile, the cells were routinely tested for mycoplasma contamination (MycoAlert PLUS Mycoplasma Detection Kit, Lonza).

In normoxic conditions, cell lines were maintained at 37°C in the humidified incubator with 5% CO₂ (Thermo Scientific). Cell lines would be plated at the desired density (60%–70%) 24 h before the placement into a hypoxia incubator (BioSpherix). The condition of hypoxia treatment was set to 1% O₂ and 5% CO₂ for 24 h.

RNA Extraction and qRT-PCR

RNA was isolated using TRI reagent solution (Sigma) followed by the on-column RNeasy mini kit and DNase treatment (Qiagen, Germany). cDNA synthesis was performed using the Transcription First Strand cDNA Synthesis Kit (Roche). qRT-PCR was performed using ABI 7900T PCR System (Applied Biosystems). Gene expression using SYBR Magic was normalized to the expression of β-actin. The primers used in the present study were supplied in Supplementary Table S4.

Multiple Immunohistochemistry Staining

Four-micron slices, cut from the paraffin block of tissues, were mounted onto charged slides and baked at 60°C for 1 h as the first step. Then, they were dewaxed with xylene for 10 min and stained with 100%, 90%, and 70% ethanol for 10 min per concentration. After being washed with deionized water for 2 min, these slides would be soaked in neutral buffered formalin for 30 min. Next, Opal manual kit (PerkinElmer) was used to stain the slides according to the manufacturer’s instructions. After nonspecific antigen sites were blocked, slides were incubated with antigen-specific primary antibody overnight at 4°C; secondary antibody incubation was performed for 1 h at room temperature. Then, dyes contained in the kit (Opal TSA) would be applied for immunofluorescence staining. We found slide stained with 3 markers plus 4 colors was the optimal choice. AR9 buffer would be used for antigen retrieval after three steps (incubation of primary antibody, secondary antibody, and dye staining) were finished. Finally, the slides were incubated with DAPI for nuclear DNA staining.

All of the primary antibodies used in mIHC are listed in Supplementary Table S5.

Statistical Analysis

Continuous and categorical (frequencies and percentages) variables were analyzed using independent t, chi-square, or 2-tailed Fisher exact tests, respectively. Meanwhile, ranked data were analyzed using the Mann-Whitney U test. The discrimination of the prediction model was assessed using ROC analysis. The OS was defined as the period from the date of surgery to the date of death due to any cause. OS between different groups was measured using the Log-rank method of K-M analysis. Cox regression analysis was used to assess time-event-dependent OS status of CC patients. The correlations of RNA expression among different hub genes were measured using spearman analysis. A p-value less than 0.05 was considered statistically significant. All statistical analyses were carried out using R (version 4.0.3; https://www.r-project.org/) and R studio (version 1.3.1093; https://www.rstudio.com/) software.

Results

Baseline Characteristics of Training and Validation Sets

All detailed information on baseline characteristics of both training and validation sets are listed in Table 1.

In the training set, 271 CC patients had complete clinical and pathological data and the remaining 11 patients only had follow-up information. Among the 271 patients, 124 patients were female (43.97%) and 147 were male (52.13%). The average age was 65.06 ± 12.70 years. Meanwhile, 224 patients (79.43%) had advanced disease (Stages II–IV), among them 39 patients (13.83%) with distant metastasis.

Also in the validation set, data of average age, the constituent ratio of gender, and pathological stage were provided. The average age and constituent ratio of pathological were comparable between the two sets. The information on pathological grade was only available in the validation set. There were 16 patients with grade 1 (9.04%), 134 with grade 2 (75.71%), and 27 with grade 3 (15.25%) disease.

Searching Procedure of Hub Hypoxia-Related Genes and Development of Risk Score Model

The complete pipeline of this study is shown in Figure 1.

FIGURE 1

We defined hub hypoxia-related genes as those that were differentially expressed between tumor and normal tissue samples and were statistically associated with the prognosis of CC patients.

According to these criteria, we first extracted the expression of 200 hypoxia-related genes from the TCGA-COAD cohort. Second, 64 DEGs (Figure 2A) and 11 prognostic genes (Figure 2B) were identified. Then, three hub genes (PPFIA4, SERPINE1, and STC2) were established as the intersection of DEGs and prognostic genes (Figure 2C).

FIGURE 2

After the confirmation of hub genes, Lasso regression analysis was used to construct the risk score model (Figures 2D,E). The formula was built as follows:

Risk score = 1.582 × PPFIA4 + 0.249×SERPINE1 + 0.279 × STC2

After the risk score of each patient was calculated, these patients were divided into two subgroups, according to the median of the risk score. The normoxia subgroup (NS) represented the one having a lower expression level of hub genes, while the hypoxia subgroup (HS) was defined as the one with a higher expression level of hub genes.

To evaluate the distinguishing ability of the risk score, principal component analysis (PCA), t-SNE (t-distributed stochastic neighbor embedding), and heatmap were used. The results of the PCA and t-SNE methods suggested that two subgroups could be separated clearly and stably according to the risk score in the training (Figure 2F) and validation (Figure 2G) sets. Figures 2H,I shows that the expression of all three hub genes was also higher in the HS compared with that in the NS in both sets.

After the grouping method was proven to be acceptable, we then investigated the differences in clinicopathological features, survival, and immune status between the two subgroups.

Different Clinicopathological and Prognostic Characteristics

First, we compared clinicopathological characteristics between two subgroups in both training and validation sets. In the training set, we found that there were more patients with metastatic lymph nodes, venous invasion, lymphatic invasion, perineural invasion, and proficient mismatch repair (pMMR) in the HS. Meanwhile, these patients had the worse pathological stage, T (invasion depth), N (lymph node metastasis), and M (distant metastasis). In the validation set, patients in the HS had poorer pathological stage and grade, and younger age (Table 1).

Second, compared with NS, HS had a higher mortality rate (19% vs. 10% in the training set; 51% vs. 31% in the validation set) (Figures 3A,B). The results of the K-M analysis suggested that patients in HS had significantly poorer OS in both training (p = 0.004) and validation (p = 0.029) sets (Figures 3C,D).

FIGURE 3

The Validation of Discrimination Ability and Stability of the Risk Score Model

In the training set, the results of the ROC analysis showed that the AUCs of our risk score model were larger than 0.6 (1-year: 0.645; 3-year: 0.700; 5-year: 0.669) (Figure 4A). Then, we used Cox regression analysis to determine whether the risk score was an independent prognostic factor. The univariate Cox analysis suggested that age, pathological stage, distant metastasis, and the risk score were associated with the prognosis (Figure 4B). The multivariate Cox analysis identified that age, distant metastasis, and the risk score were the independent prognostic factors (Figure 4C).

FIGURE 4

In the validation set, the AUCs of 1, 3, and 5 years were all above 0.5 (1-year: 0.644; 3-year: 0.596; 5-year: 0.583) (Figure 4D). Through univariate and multivariate Cox analysis, we found that the pathological stage and the risk score were the independent prognostic factors (Figures 4E,F).

Taken together, the results confirmed that the risk score model based on the hub hypoxia-related genes had valuable discrimination ability and stability. It could predict the prognosis of patients with colon cancer accurately in the different study populations.

Visualization and Calibration of the Risk Score Model

We decide to use the clinicopathological features (age, M, and pathological stage) in the TCGA-COAD cohort, along with risk score, to build a nomogram for clinical application.

First, in Figure 5A, the result of calibration analysis showed that the blue, red, and green lines, which represented the performance of 1-, 3-, and 5-year prognostic prediction, were just close to the diagonal. It meant that the risk score was a stable prognostic model. Next, we compared the prognostic ability of three models, including the risk score, clinical characteristics (age, M, and pathological stage), and risk score plus clinical characteristics using Decisive Curve Analysis (DCA). The results showed that the performance of the multiple-factor model (risk score plus clinical characteristics) was slightly better than that of the single-factor model (the risk score model alone); however, both of their performance was significantly better than that of the sole clinical characteristics model (Figure 5B).

FIGURE 5

Since the diagnostic ability of the risk score model was validated to be valuable and stable, we then visualized it using a “nomogram” (Figure 5C). Patients with colon cancer could predict their 1-, 3-, and 5-year OS according to the information of age, pathological stage, status of distant metastasis, and our risk score model. It also could help doctors to predict the prognosis of patients with colon cancer accurately and easily. After the construction of the nomogram, we also evaluated the prediction power of the nomogram; the AUCs of 1-, 3-, and 5-year OS were 0.871, 0.788, and 0.807, respectively (Figure 5D).

Meanwhile, we also did calibration analysis, DCA, and ROC analysis, and built a nomogram in the validation set. The results showed that the risk model had stable prediction performance and the combined usage of the risk model with clinicopathological features had better prediction power than that of mono-marker (Supplementary Figure S1). The AUCs of 1-, 3-, and 5-year OS of nomogram in the validation set were 0.768, 0.724, and 0.676, respectively (Supplementary Figure S1D).

Analyses of Differentially Expressed Genes and Functional Enrichment

After the confirmation of the correlation between risk score and the prognosis of CC patients, we then investigated different bio-functions between NS and HS. First, we used the limma package to identify DEGs between NS and HS. There were 72 DEGs (log₂|FC|>1 & adjusted p < 0.05) between different subgroups. Because bio-functional analysis should be built on enough DEGs, we reset the standard of DEGs to |FC|>1.5. As a result, there were 163 DEGs (4 down and 159 up in the HS) (Figure 6A).

FIGURE 6

Through GO and KEGG analyses, we found that the bio-functions of upregulated genes in the HS were enriched in proliferation, differentiation, and tumorigenesis-related signaling pathways, including PI3K-Akt, Hedgehog, and Wnt signaling pathways. Meanwhile, these DEGs were also involved in Growth factor and Cytokine binding and extracellular matrix reorganization and reconstruction (Figure 6B). Moreover, GSEA analysis revealed that compared with NS, HS was enriched in hypoxia, epithelial–mesenchymal transition (EMT), angiogenesis, and KRAS, IL2/STAT5, and Hedgehog signaling pathways (Figure 6C).

Taken together, HS had enriched bio-functions in tumorigenesis, proliferation, and differentiation. Meanwhile, some signaling pathways that would induce the resistance of antitumor drugs (EMT and KRAS) were also found to be involved in HS. These might partially explain why the patients in HS had worse pathological features and prognosis.

Tumor Immune Microenvironment and the Expression of Different Molecules

Previous studies have found that the infiltration of types of immune cells was significantly associated with the clinical outcomes of CRC patients (Galon et al., 2006; Craig et al., 2020; Picard et al., 2020). Based on the previous analyses, we found that the OS between NS and HS were statistically different. We wondered whether it was also associated with different TIME between the two subgroups. Then, we analyzed the immune-cell infiltration of each CC sample in both the training and validation sets using CIBERSORT.

The result showed that the HS had a significantly smaller number of CD8⁺ T and resting memory CD4⁺ T cells and a larger number of M0 macrophages in the training (Figures 7A,B) and validation (Supplementary Figure S2A) sets. While in the training set, we could also find that the number of plasma and resting dendritic cells was smaller in the HS than that in the NS. Compared with that, HS had a higher infiltration level of M2 macrophages and neutrophils (Figure 7B). Moreover, to ensure the stability of the results, we then used another algorithm (MCP-counter) to evaluate the infiltration of immune cells in different subgroups. The results showed that the infiltration of T cells, CD8⁺ T cells, and cytotoxic lymphocytes was significantly fewer in HS, whereas the infiltration of monocytic cells, endothelial cells, and fibroblasts was enriched in HS (Figure 7C). These results revealed that in the hypoxic microenvironment, the infiltration of pro-inflammation immune cells was significantly impeded, whereas the anti-inflammation immune cells and fibroblast were enriched in the hypoxic conditions that would also hinder the proliferation and migration of antitumor immune cells (Vitale et al., 2019; Davidson et al., 2021).

FIGURE 7

We then investigated the differential expression of negative regulatory immune-related genes, immune checkpoint molecules, human leukocyte antigen, CXCL, and CCL chemokines. First, we found that numerous negative regulator immune-related genes, including some immune checkpoint genes, were highly expressed in HS (Figure 7D). Second, we also found that almost all of the immune checkpoint molecules were highly expressed in HS, except LAG3 and ICOSLG (Figure 7E). Multiple studies also found that the highly expressed immune checkpoint molecules would inhibit the antitumor immunity, resulting in a poor prognosis for colorectal cancer patients (Neupane et al., 2021; Kudo-Saito et al., 2021; Gordon et al., 2017). Besides, the expression of some human leukocyte antigens (HLAs) was lowly expressed in HS, which was also consistent with the results of fewer infiltration of pro-inflammation immune cells in HS (Figure 7F). A plenty of studies have confirmed that impediments to processing HLA would hinder the identification, migration, and infiltration ability of tumor-infiltrating lymphocytes (TILs), thus facilitating the proliferation and invasion of malignancy (Dong et al., 2021; Maggs et al., 2021; Kawazu et al., 2022).

Next, we found that most chemokines, including CXCL5, 6, 8, 9, 10, 11, and 12, CCL1, 2, 3, 4, 7, 8, 11, 12, 13, 19, 21, 23, 26, and 28, were significantly highly expressed in HS (Figures 7G,H). Other studies revealed that the upregulated expression of chemokines, including CXC and CC chemokine families, could impede the TILs infiltrating into a tumor, support the growth of malignant cells, and facilitate the migration of myeloid-derived suppressor cells (MDSCs), which would cause the drug resistance to the chemotherapy and immunotherapy (Korbecki et al., 2020; Bullock and Richmond, 2021; Matsuo et al., 2021).

All in all, these findings revealed that the TIME and the expression of immune checkpoint molecules, immune-related genes, and chemokines were quite different between NS and HS, which might be associated with different pathological features and prognosis between the two subgroups. Next, based on the findings mentioned above, we would evaluate the differences in therapeutic responses between NS and HS.

Drug Response of Immune Checkpoint Blockade, Cytotoxic, and Targeted Medicine

We used several methods to evaluate the therapeutic responses between NS and HS, including the website tool (ImmuCellAI: http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) for ICB and pRRophetic package for targeted and cytotoxic medicine.

First, in NS, there were 38 CC (26.9% of all CC patients in NS) patients who responded to ICB therapy. Compared with that, only 17 CC (12.1% of all CC patients in HS) patients would benefit from ICB therapy. The difference between the two subgroups was statistically significant (p = 0.002) (Figure 8A). Second, we evaluated three cytotoxic drugs that were most widely used in clinical application, including 5-fluorouracil (5-FU), cisplatin, and gemcitabine. The result showed that the half inhibitory concentration (IC50) of cisplatin was statistically lower in NS (Figure 8B). For targeted medicine, we found that the IC50s of bosutinib, imatinib, sorafenib, and sunitinib were all lower in NS (Figure 8C).

FIGURE 8

Although most targeted drugs were not approved to be used as the first-line therapy for CRC by FDA, some studies revealed that tyrosine kinase inhibitor (TKI), including bosutinib, imatinib, and sunitinib, could enhance the infiltration of cytotoxic and effector T cells, which would directly affect the efficacy of immunotherapy (Roulleaux Dugage et al., 2021; Tazzari et al., 2021; Hirata et al., 2022). These studies revealed that the application of TKI might be positively related to the infiltration of pro-inflammation immune cells, which meant that the combination therapy of TKI and ICB might receive better efficacy than monotherapy. Recently, in a CRC mouse model, researchers found that the combination therapy of TKI and ICB could reduce tumor-stromal volume and increase the infiltration of CD8⁺ T cells and the activation of immune-related pathways (Yorita et al., 2021).

Analyses of Correlated Expression and Quantitative Real-Time PCR for Hub Genes

Although our previous work revealed that the expression of three hub genes (PPFIA4, SERPINE1, and STC2) might be positively related, we still wondered whether there were direct correlations among these genes. The results suggested that all three genes were indeed positively expressed (Figure 9A). Meanwhile, we found that the correlation between the expression of PPFIA4 and SERPINE1 was statistically highest (Figures 9B–D). To ensure the accuracy of the results, we also used the website tool TIMER to evaluate the correlations among three hub genes. The results confirmed that these genes were statistically positively related (Figures 9E–G).

FIGURE 9

All 200 genes included in our analysis were in the hallmark hypoxia gene set of the GSEA database. Next, we used qPCR to confirm whether three hub genes were highly expressed in the hypoxic conditions. We used five types of colon cell lines, including 1 colon epithelial cell line (FHC) and 4 CRC cell lines (HCT-8, RKO, SW480, and SW620). Among them, HCT-8 and RKO were dMMR/MSI cell lines, while SW480 and SW620 were pMMR/MSS. For the hypoxic conditions, cell lines were cultivated in a hypoxia incubator with 1% O₂ and 5% CO₂ for 24 h. The result showed that all three hub genes were statistically highly expressed in the hypoxic conditions (Figure 9H).

Analysis of Each Hub Gene

All previous analyses were about the risk score model that was constructed with hub genes and related coefficients. We then analyzed the correlation of the expression of single hub gene with clinicopathological features and immune infiltration (Supplementary Figures S3–5).

The results suggested that a higher level of expression of all three hub genes was found in the tumor samples and was associated with a poorer prognosis for CC patients. Meanwhile, all three genes were positively related to the infiltration of macrophages, neutrophils, and NK cells. However, PPFIA4 and STC2 were negatively associated with CD8⁺ T cells, T cells, and cytotoxic cells, which were following the previous results of TIME in different subgroups that CC patients with higher risk scores had a lower level of infiltration of CD8⁺ T cells (Figure 7).

Multiple Immunohistochemistry Staining

Finally, we wanted to validate the results of the CIBERSORT analysis, which suggested that the expression of hub hypoxia-related genes (PPFIA4, SERPINE1, and STC2) was associated with the infiltration of some immune cells. We collected 35 colon cancer tissue samples from the Department of General Surgery, Changzhou Wujin Hospital. All the samples were fixed by formalin and embedded in paraffin.

The results showed that the expression of three hub genes was positively related (Figure 10A). NS had a relatively lower expression of hypoxia-related genes. We divided 35 samples into a lower expression subgroup (NS: 18) and a higher expression subgroup (HS: 17), according to the median of the number of positively expressed cells. Then, we checked their immune status, including CD8⁺ T cells, Treg cells, the expression of PD-L1, M1, and M2 macrophages. In Figures 10B,C, the results suggested that the NS had larger colonization of CD8⁺ T cells and M1 macrophages (CD80) than HS, which was insistence with our previous results in Figures 7A,B. Meanwhile, similar to the results shown in Figures 7B,E, M2 macrophages (CD163) and PD-L1 were statistically less in NS, compared with those in HS. Moreover, we also compared the pathological feature between NS and HS of the 35 colon cancer patients. The results showed that the OS and pathological variates, including pathological stage, invasion depth, and lymphocytic metastasis, were significantly better than those of the HS (Figures 10D,E).

FIGURE 10

Combined with these results, we identified that NS, which represented normoxia TIME, had relatively “hot” TIME. It had enriched colonization of pro-inflammation immune cells (CD8⁺ T cells, M1 macrophages), less population of anti-inflammation immune cells (M2 macrophages), and higher expression of immune checkpoint molecules (PD-L1). Compared with HS, the TIME of NS would be more suitable for TILs to survive and exhibit their tumor-killing function. On the other hand, the results of mIHC were similar to those analyzed by CIBERSORT, indicating its stable performance.

Discussion

In the present study, we have identified a hypoxia-related gene signature, which contained three hub genes, to classify CC patients into NS and HS, and subsequently predicted the different TIME and prognosis between different subgroups. Based on this signature, a risk score model was constructed and proven to be a valuable and stable prognostic tool for CC patients. Moreover, we identified that the TIME and the gene expression were quite different, which might induce distinct drug responses between NS and HS. Based on ImmuCellAI and pRRophetic package, we found that patients in NS were more sensitive to ICB and targeted therapies. Finally, three hub genes (PPFIA4, SERPINE1, and STC2) were confirmed to be highly expressed in the hypoxic conditions by qPCR. Also, the infiltrations of CD8⁺ T cells and M1 macrophages were proven to be negatively related to the expression of these genes by mIHC.

Our present work identified that hypoxia-related genes and the risk score model built on them were significantly correlated with the prognosis of CC patients. Similar results were found in ovarian cancer (Chen et al., 2021), triple-negative breast cancer (Yang et al., 2021a), osteosarcoma in children (Jiang et al., 2021a), acute myeloid leukemia (Jiang et al., 2021b), and so on. It indicated that along with other factors, hypoxia might play an important role in the development and progression of cancer.

In our risk score model, three hub genes, including PPFIA4, SERPINE1, and STC2, were all identified to be differentially expressed between normal and tumor tissues and be closely associated with the prognosis and infiltration of immune cells in CC patients. PPFIA4, belonging to the PPFIA family of kinesin-cargo linkers, was first identified and characterized in silico in 2003 (Katoh and Katoh, 2003). Recently, it has been proved to be related to CC cell proliferation and migration by enhancing tumor glycolysis (Huang et al., 2021). Besides, previous studies suggested that PPFIA4 was also the key prognostic gene in thyroid and prostate cancer (Xu et al., 2021a; Xu et al., 2021b). Compared with PPFIA4, the role of SERPINE1 in CC had been more deeply investigated. Some studies found that the expression of SERPINE1 was negatively associated with tumor grade and response to adjuvant therapy of CC patients (Halamkova et al., 2011; Cheng et al., 2018). Meanwhile, SERPINE1 has also been proven to play an important role in remodeling TME and enhancing tumor progression in CC (Wang et al., 2021). For STC2, it was recognized as a regulator in CC cell biological processes, and silencing STC2 could effectively suppress cancer cell proliferation, survival, and migration (Li et al., 2019). Moreover, higher expression of STC2 mRNA in tumor tissues was correlated with larger tumor size, presence of venous invasion, lymphatic invasion, distant metastasis, and poorer prognosis of CC patients (Watanabe et al., 2021). However, many studies we mentioned above were only based on the bioinformatics analyses, which were not stable and persuasive enough. Therefore, besides bioinformatics analyses, we also used cell lines and tumor tissues to perform qPCR and mIHC for validation.

The limit of nutrients and oxygen, which also restrains the proliferation of tumor cells, stimulates tumors to enhance the growth of new vasculatures. However, these newly formed vessels are leaky for their discontinuous endothelium, which will induce high permeability and permeation (Majidpoor and Mortezaee, 2021). The disorganized vasculatures, along with the high level of metabolic rate and the low efficiency of ATP producing method of tumor cells, cause a severe hypoxic condition in TME (Bertout et al., 2008). In addition, the hypoxia TME simultaneously promotes the famous “Warburg effect,” enhancing glycolysis and lactic acid production catalyzed by the lactate dehydrogenase A (LDH-A) (Harris, 2002). Subsequently, it will result in acidic pH, which impairs cytotoxicity and proliferation of types of immune cells by reducing their chemotaxis, respiratory activity, and bactericidal ability (Jing et al., 2019). Taken together, the hypoxia and acidic TME greatly suppress the antitumor immune function and thus induce tumor survival and metastasis. Among numerous immune cells infiltrating in TME, macrophages are the principal component, which can differentiate into the tumor-associated macrophages (TAMs) that have been identified to be preferentially located in almost all tumor hypoxic regions (Hegde et al., 2021). Different from M1 TAMs, M2 TAMs play an anti-inflammatory role in the TME by secreting immunosuppressive molecules, including IL-10, human leukocyte antigen G (HLA-G), TGF-β (Komohara et al., 2016). They directly with MDSCs restrain the infiltration of antitumor T cells and their secretion of IFN-γ. In the present study, less infiltration of CD8⁺ T cells and more infiltration of M2 macrophages were found in the high-risk group in both bioinformatics and mIHC analyses. Meanwhile, M2 TAMs were proven to express increased programmed cell death ligand 1 (PD-L1). It has been established that increased expression of PD-L1 was related to poor prognosis in numerous malignancies (Pérez-Ruiz et al., 2020). Accordant with the previous study, we also found that the expression of immune checkpoint molecules was higher in the high-risk group, which indicated that they might be associated with hypoxia condition. A recent study proved that tumor cells might escape immune attacks from both innate and adaptive immune systems by secreting hypoxia-inducible factor 1 (HIF-1) (You et al., 2021). Since suppressive TAMs, immune checkpoint molecules, and HIF-1 were all important negative factors, blocking them has predictably received promising results in enhancing the infiltration of tumor-infiltrating lymphocytes (TILs), thus improving their tumor-killing effect (Yang et al., 2021b; Caushi et al., 2021; Yap et al., 2021; You et al., 2021).

Based on the knowledge that hypoxia was a key barrier to antitumor immunity, studies have focused on how to target the hypoxic metabolic production or reverse the hypoxic condition in TME. First, enhancing tumor oxygenation is an option, which could be applied by carbogen breathing and intervention to reduce O₂ consumption by the tumor (Kheir et al., 2012; Zou et al., 2018). Second, hypoxia-activated prodrugs (HAPs) were designed to specifically target those hypoxic tumor cells. These HAPs could bring both genotoxic agents and non-genotoxic effectors (Penketh et al., 2012; Skwarska et al., 2021). Third, targeting acidosis or hypoxia-acidosis-related pathways in TME is another choice (Singleton et al., 2021). Moreover, the combination therapy could also enhance the antitumor immunity and reduce drug toxicity and resistance. For instance, the objective response rates (ORR) of monotherapy of PD-1 and CTLA-4 blockers were 10%–16% for ipilimumab and 10%–40% for nivolumab and pembrolizumab (Robert et al., 2011; Ribas et al., 2015; Robert et al., 2015). While the combination of them could significantly increase the ORR to 61% and reduce the incidence of grade 3–4 adverse events to 46% (Long et al., 2017). The combination of ICBs with chemotherapy, targeted therapy, radiation, or intratumoral therapy also showed gratifying therapeutic results in treating different types of malignancies (Meric-Bernstam et al., 2021). The reason is these therapeutic methods could enhance the activation and infiltration of TILs by inhibiting angiogenesis, normalizing vasculature, reconstructing immunosupportive TME, increasing antigen presentation, co-stimulating molecules, and so on.

Recently, thanks to the advance in bioinformatics, we have deeper insight into the genomics of human beings. The function of more gene signatures has been identified. Using the existing gene signature, scientists can filter out hub genes and build risk models that have prognostic power and can depict different tumor microenvironments of different cancer patients.

These risk models were built using different statistical methods. For instance, Guan et al. (2020) established a prediction model that could separate gastric cancer patients into two different subgroups using the immune-related gene signature (Guan et al., 2020). They used the ssGSEA score and the hierarchical clustering algorithm. One subgroup had a higher expression of the immune-related score and better prognosis and another had a lower score and poorer prognosis. They successfully associated the tumor immune microenvironment with the prognosis in the population of gastric cancer patients. Based on the immune checkpoint-related gene signature, another study built a risk model to predict the prognosis of hepatocellular carcinoma patients using Lasso and Cox regression analyses (Zhao et al., 2020). Their model could also divide liver cancer patients into 2 subgroups, which had different prognosis and tumor immune microenvironments. The third example is that Bagaev et al. (2021) built a model, which had 29 knowledge-based functional gene signatures, to separate cancer patients into 4 subtypes, including immune-enriches (fibrotic), immune-enriched (non-fibrotic), fibrotic, and depleted (Bagaev et al., 2021). They used the ssGSEA score and the Louvain clustering method. These subtypes had a quite different infiltration of immune cells and prognosis. The most important thing is that this model is pan-cancer conserved, which means that it can be used in most types of cancer and has great clinical application potential. Like the previous studies, we also used Lasso and Cox regression analyses to identify hub genes and built a risk model subsequently. We believe that along with the advance in bioinformatics and statistics, more powerful prediction models will be developed.

There are several limitations in the present study. First, all analyses were based on retrospective data from public databases, which will induce recall and selection biases. Second, although related work to eliminate the batch effect has been done, the potential risk still exists when analyses are based on data from two different databases. Finally, our biomolecular experiments were just to verify the results of bioinformatics analyses. They have not deeply uncovered the bio- and molecular mechanism of hypoxia involved in the development and progression of CC. However, our further work will continuously focus on this field.

Conclusion

In summary, based on the hypoxia-related genes, we constructed a risk score model to predict the prognosis of colon cancer patients. Moreover, we deeply analyzed the differences, including functional enrichment, infiltration of immune cells, expression of different genes (immune checkpoint genes, human leukocyte antigen, CXCL, and CCL chemokines), and the therapeutic responses, between high-risk and low-risk subgroups. Finally, we performed qPCR and multiple immunohistochemistry (mIHC) for validation. This stud might provide new insights into the association among hypoxia, clinical prognosis, TIME, and therapy.

References

• 1

  AndersonN. M.SimonM. C. (2020). The Tumor Microenvironment. Curr. Biol.30 (16), R921–R925. Epub 2020/08/19PubMed PMID: 32810447. 10.1016/j.cub.2020.06.081PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 2

  BagaevA.KotlovN.NomieK.SvekolkinV.GafurovA.IsaevaO.et al (2021). Conserved Pan-Cancer Microenvironment Subtypes Predict Response to Immunotherapy. Cancer Cell39 (6), 845–865. Epub 2021/05/22PubMed PMID: 34019806. 10.1016/j.ccell.2021.04.014PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 3

  BechtE.GiraldoN. A.LacroixL.ButtardB.ElarouciN.PetitprezF.et al (2016). Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol.17 (1), 218. Epub 2016/10/22PubMed PMID: 27765066; PubMed Central PMCID: PMCPMC5073889. 10.1186/s13059-016-1070-5PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 4

  BertoutJ. A.PatelS. A.SimonM. C. (2008). The Impact of O2 Availability on Human Cancer. Nat. Rev. Cancer8 (12), 967–975. Epub 2008/11/07PubMed PMID: 18987634; PubMed Central PMCID: PMCPMC3140692. 10.1038/nrc2540PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 5

  BrayF.FerlayJ.SoerjomataramI.SiegelR. L.TorreL. A.JemalA. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A Cancer J. Clin.68 (6), 394–424. Epub 2018/09/13PubMed PMID: 30207593. 10.3322/caac.21492PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 6

  BullockK.RichmondA. (2021). Suppressing MDSC Recruitment to the Tumor Microenvironment by Antagonizing CXCR2 to Enhance the Efficacy of Immunotherapy. Cancers13 (24), 6293. Epub 2021/12/25PubMed PMID: 34944914; PubMed Central PMCID: PMCPMC8699249. 10.3390/cancers13246293PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 7

  CaushiJ. X.ZhangJ.JiZ.VaghasiaA.ZhangB.HsiueE. H.-C.et al (2021). Transcriptional Programs of Neoantigen-specific TIL in Anti-PD-1-treated Lung Cancers. Nature596 (7870), 126–132. Epub 2021/07/23PubMed PMID: 34290408. 10.1038/s41586-021-03752-4PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 8

  ChenX.LanH.HeD.XuR.ZhangY.ChengY.et al (2021). Multi-Omics Profiling Identifies Risk Hypoxia-Related Signatures for Ovarian Cancer Prognosis. Front. Immunol.12, 645839. Epub 2021/08/06PubMed PMID: 34349753; PubMed Central PMCID: PMCPMC8327177. 10.3389/fimmu.2021.645839PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 9

  ChengX.HuM.ChenC.HouD. (2018). Computational Analysis of mRNA Expression Profiles Identifies a Novel Triple-Biomarker Model as Prognostic Predictor of Stage II and III Colorectal Adenocarcinoma Patients. Cmar10, 2945–2952. Epub 2018/09/15PubMed PMID: 30214289; PubMed Central PMCID: PMCPMC6118290. 10.2147/cmar.S170502PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 10

  CraigS. G.HumphriesM. P.AlderdiceM.BinghamV.RichmanS. D.LoughreyM. B.et al (2020). Immune Status Is Prognostic for Poor Survival in Colorectal Cancer Patients and Is Associated with Tumour Hypoxia. Br. J. Cancer123 (8), 1280–1288. Epub 2020/07/21PubMed PMID: 32684627. 10.1038/s41416-020-0985-5PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 11

  DavidsonS.ColesM.ThomasT.KolliasG.LudewigB.TurleyS.et al (2021). Fibroblasts as Immune Regulators in Infection, Inflammation and Cancer. Nat. Rev. Immunol.21, 704–717. Epub 2021/04/30PubMed PMID: 33911232. 10.1038/s41577-021-00540-zPubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 12

  DongL.YangX.WangY.JinY.ZhouQ.ChenG.et al (2021). Key Markers Involved in the Anticolon Cancer Response of CD8+ T Cells through the Regulation of Cholesterol Metabolism. J. Oncol.2021, 1–11. Epub 2021/12/04PubMed PMID: 34858500; PubMed Central PMCID: PMCPMC8632400. 10.1155/2021/9398661CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 13

  GalonJ.CostesA.Sanchez-CaboF.KirilovskyA.MlecnikB.Lagorce-PagesC.et al (2006). Type, Density, and Location of Immune Cells within Human Colorectal Tumors Predict Clinical Outcome, Science313, 1960–1964. Epub 2006/09/30PubMed PMID: 17008531. 10.1126/science.1129139PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 14

  GeeleherP.CoxN.HuangR. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One9 (9), e107468. Epub 2014/09/18PubMed PMID: 25229481; PubMed Central PMCID: PMCPMC4167990. 10.1371/journal.pone.0107468PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 15

  GordonS. R.MauteR. L.DulkenB. W.HutterG.GeorgeB. M.McCrackenM. N.et al (2017). PD-1 Expression by Tumour-Associated Macrophages Inhibits Phagocytosis and Tumour Immunity. Nature545 (7655), 495–499. Epub 2017/05/18PubMed PMID: 28514441; PubMed Central PMCID: PMCPMC5931375. 10.1038/nature22396PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 16

  GuanX.XuZ.-Y.ChenR.QinJ.-J.ChengX.-D. (2020). Identification of an Immune Gene-Associated Prognostic Signature and its Association with a Poor Prognosis in Gastric Cancer Patients. Front. Oncol.10, 629909. Epub 2021/02/26PubMed PMID: 33628738; PubMed Central PMCID: PMCPMC7898907. 10.3389/fonc.2020.629909PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 17

  HalamkovaJ.KissI.PavlovskyZ.TomasekJ.JarkovskyJ.CechZ.et al (2011). Clinical Significance of the Plasminogen Activator System in Relation to Grade of Tumor and Treatment Response in Colorectal Carcinoma Patients. neo58 (5), 377–385. Epub 2011/07/13PubMed PMID: 21744990. 10.4149/neo_2011_05_377PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 18

  HarrisA. L. (2002). Hypoxia - a Key Regulatory Factor in Tumour Growth. Nat. Rev. Cancer2 (1), 38–47. Epub 2002/03/21PubMed PMID: 11902584. 10.1038/nrc704PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 19

  HegdeS.LeaderA. M.MeradM. (2021). MDSC: Markers, Development, States, and Unaddressed Complexity. Immunity54 (5), 875–884. Epub 2021/05/13PubMed PMID: 33979585. 10.1016/j.immuni.2021.04.004PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 20

  Hiam-GalvezK. J.AllenB. M.SpitzerM. H. (2021). Systemic Immunity in Cancer. Nat. Rev. Cancer21, 345–359. Epub 2021/04/11PubMed PMID: 33837297; PubMed Central PMCID: PMCPMC8034277. 10.1038/s41568-021-00347-zPubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 21

  HirataA.SawaiE.HenmiM.ShibasakiC.MizoguchiY.NarumiK.et al (2022). Imatinib Mesylate Exerted Antitumor Effect by Promoting Infiltration of Effector T Cells in Tumor. Biol. Pharm. Bull.45 (1), 34–41. Epub 2022/01/05PubMed PMID: 34980779. 10.1248/bpb.b21-00493PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 22

  HsuT.-S.LinY.-L.WangY.-A.MoS.-T.ChiP.-Y.LaiA. C.-Y.et al (2020). HIF-2α Is Indispensable for Regulatory T Cell Function. Nat. Commun.11 (1), 5005. Epub 2020/10/08PubMed PMID: 33024109; PubMed Central PMCID: PMCPMC7538433. 10.1038/s41467-020-18731-yPubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 23

  HuangJ.YangM.LiuZ.LiX.WangJ.FuN.et al (2021). PPFIA4 Promotes Colon Cancer Cell Proliferation and Migration by Enhancing Tumor Glycolysis. Front. Oncol.11, 653200. Epub 2021/06/08PubMed PMID: 34094943; PubMed Central PMCID: PMCPMC8173052. 10.3389/fonc.2021.653200PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 24

  JiangF.MaoY.LuB.ZhouG.WangJ. (2021). A Hypoxia Risk Signature for the Tumor Immune Microenvironment Evaluation and Prognosis Prediction in Acute Myeloid Leukemia. Sci. Rep.11 (1), 14657. Epub 2021/07/21PubMed PMID: 34282207; PubMed Central PMCID: PMCPMC8289869. 10.1038/s41598-021-94128-1PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 25

  JiangF.MiaoX.-L.ZhangX.-T.YanF.MaoY.WuC.-Y.et al (2021). A Hypoxia Gene-Based Signature to Predict the Survival and Affect the Tumor Immune Microenvironment of Osteosarcoma in Children. J. Immunol. Res.2021, 1–13. Epub 2021/08/03PubMed PMID: 34337075; PubMed Central PMCID: PMCPMC8299210. 10.1155/2021/5523832CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 26

  JingX.YangF.ShaoC.WeiK.XieM.ShenH.et al (2019). Role of Hypoxia in Cancer Therapy by Regulating the Tumor Microenvironment. Mol. Cancer18 (1), 157. Epub 2019/11/13PubMed PMID: 31711497; PubMed Central PMCID: PMCPMC6844052. 10.1186/s12943-019-1089-9PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 27

  KatohM.KatohM. (2003). Identification and Characterization of Human PPFIA4 Gene In Silico. Int. J. Mol. Med.12 (6), 1009–1014. Epub 2003/11/13. PubMed PMID: 14612982. 10.3892/ijmm.12.6.1009PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 28

  KawazuM.UenoT.SaekiK.SaxN.TogashiY.KanasekiT.et al (2022). HLA Class I Analysis Provides Insight into the Genetic and Epigenetic Background of Immune Evasion in Colorectal Cancer with High Microsatellite Instability. Gastroenterology162162 (3), 799799–812812. Epub 2021/10/24PubMed PMID: 34687740. 10.1053/j.gastro.2021.10.010CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 29

  KaymakI.WilliamsK. S.CantorJ. R.JonesR. G. (2021). Immunometabolic Interplay in the Tumor Microenvironment. Cancer Cell39 (1), 28–37. Epub 2020/10/31PubMed PMID: 33125860; PubMed Central PMCID: PMCPMC7837268. 10.1016/j.ccell.2020.09.004PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 30

  KheirJ. N.ScharpL. A.BordenM. A.SwansonE. J.LoxleyA.ReeseJ. H.et al (2012). Oxygen Gas-Filled Microparticles Provide Intravenous Oxygen Delivery. Sci. Transl. Med.4 (140), 140ra88. Epub 2012/06/30PubMed PMID: 22745438. 10.1126/scitranslmed.3003679PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 31

  KomoharaY.FujiwaraY.OhnishiK.TakeyaM. (2016). Tumor-associated Macrophages: Potential Therapeutic Targets for Anti-cancer Therapy. Adv. Drug Deliv. Rev.99 (Pt B), 180–185. Epub 2015/12/02PubMed PMID: 26621196. 10.1016/j.addr.2015.11.009PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 32

  KorbeckiJ.KojderK.SimińskaD.BohatyrewiczR.GutowskaI.ChlubekD.et al (2020). CC Chemokines in a Tumor: A Review of Pro-cancer and Anti-cancer Properties of the Ligands of Receptors CCR1, CCR2, CCR3, and CCR4. Ijms21 (21), 8412. Epub 2020/11/14PubMed PMID: 33182504; PubMed Central PMCID: PMCPMC7665155. 10.3390/ijms21218412PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 33

  Kudo-SaitoC.OgiwaraY.ImazekiH.BokuN.UemuraY.ZhangR.et al (2021). CD11b(+)DIP2A(+)LAG3(+) Cells Facilitate Immune Dysfunction in Colorectal Cancer. Am. J. Cancer Res.11 (11), 5428–5439. Epub 2021/12/08. PubMed PMID: 34873470; PubMed Central PMCID: PMCPMC8640801. PubMed Abstract | Google Scholar

  • Google Scholar

• 34

  LiB.SeversonE.PignonJ.-C.ZhaoH.LiT.NovakJ.et al (2016). Comprehensive Analyses of Tumor Immunity: Implications for Cancer Immunotherapy. Genome Biol.17 (1), 174. Epub 2016/08/24PubMed PMID: 27549193; PubMed Central PMCID: PMCPMC4993001. 10.1186/s13059-016-1028-7PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 35

  LiQ.ZhouX.FangZ.PanZ. (2019). Effect of STC2 Gene Silencing on Colorectal Cancer Cells. Mol. Med. Rep.20 (2), 977–984. Epub 2019/06/08PubMed PMID: 31173256; PubMed Central PMCID: PMCPMC6625197. 10.3892/mmr.2019.10332PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 36

  LiT.FanJ.WangB.TraughN.ChenQ.LiuJ. S.et al (2017). TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res.77 (21), e108–e110. Epub 2017/11/03PubMed PMID: 29092952; PubMed Central PMCID: PMCPMC6042652. 10.1158/0008-5472.Can-17-0307PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 37

  LongG. V.AtkinsonV.CebonJ. S.JamesonM. B.FitzharrisB. M.McNeilC. M.et al (2017). Standard-dose Pembrolizumab in Combination with Reduced-Dose Ipilimumab for Patients with Advanced Melanoma (KEYNOTE-029): an Open-Label, Phase 1b Trial. Lancet Oncol.18 (9), 1202–1210. Epub 2017/07/22PubMed PMID: 28729151. 10.1016/s1470-2045(17)30428-xPubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 38

  MaggsL.SadagopanA.MoghaddamA. S.FerroneS. (2021). HLA Class I Antigen Processing Machinery Defects in Antitumor Immunity and Immunotherapy. Trends Cancer7 (12), 1089–1101. Epub 2021/09/08PubMed PMID: 34489208; PubMed Central PMCID: PMCPMC8651070. 10.1016/j.trecan.2021.07.006PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 39

  MajidpoorJ.MortezaeeK. (2021). Angiogenesis as a Hallmark of Solid Tumors - Clinical Perspectives. Cell Oncol.44 (4), 715–737. Epub 2021/04/10PubMed PMID: 33835425. 10.1007/s13402-021-00602-3CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 40

  MatsuoK.YoshieO.NakayamaT. (2021). Multifaceted Roles of Chemokines and Chemokine Receptors in Tumor Immunity. Cancers13 (23), 6132. Epub 2021/12/11PubMed PMID: 34885241; PubMed Central PMCID: PMCPMC8656932. 10.3390/cancers13236132PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 41

  Meric-BernstamF.LarkinJ.TaberneroJ.BoniniC. (2021). Enhancing Anti-tumour Efficacy with Immunotherapy Combinations. Lancet397 (10278), 1010–1022. Epub 2020/12/08PubMed PMID: 33285141. 10.1016/s0140-6736(20)32598-8PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 42

  MiaoY.-R.XiaM.LuoM.LuoT.YangM.GuoA.-Y. (2021). ImmuCellAI-mouse: a Tool for Comprehensive Prediction of Mouse Immune Cell Abundance and Immune Microenvironment Depiction. Bioinformatics38, 785–791. Epub 2021/10/13PubMed PMID: 34636837. 10.1093/bioinformatics/btab711CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 43

  MiaoY. R.ZhangQ.LeiQ.LuoM.XieG. Y.WangH.et al (2020). ImmuCellAI: A Unique Method for Comprehensive T‐Cell Subsets Abundance Prediction and its Application in Cancer Immunotherapy. Adv. Sci.7 (7), 1902880. Epub 2020/04/11PubMed PMID: 32274301; PubMed Central PMCID: PMCPMC7141005. 10.1002/advs.201902880PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 44

  NeupaneP.MimuraK.NakajimaS.OkayamaH.ItoM.Thar MinA. K.et al (2021). The Expression of Immune Checkpoint Receptors and Ligands in the Colorectal Cancer Tumor Microenvironment. Anticancer Res.41 (10), 4895–4905. Epub 2021/10/02PubMed PMID: 34593437. 10.21873/anticanres.15303PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 45

  NewmanA. M.SteenC. B.LiuC. L.GentlesA. J.ChaudhuriA. A.SchererF.et al (2019). Determining Cell Type Abundance and Expression from Bulk Tissues with Digital Cytometry. Nat. Biotechnol.37 (7), 773–782. Epub 2019/05/08PubMed PMID: 31061481; PubMed Central PMCID: PMCPMC6610714. 10.1038/s41587-019-0114-2PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 46

  PenkethP. G.ShyamK.BaumannR. P.IshiguroK.PatridgeE. V.ZhuR.et al (2012). A Strategy for Selective O6-Alkylguanine-DNA Alkyltransferase Depletion under Hypoxic Conditions. Chem. Biol. Drug Des.80 (2), 279–290. Epub 2012/05/05PubMed PMID: 22553921; PubMed Central PMCID: PMCPMC3399964. 10.1111/j.1747-0285.2012.01401.xPubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 47

  Pérez-RuizE.MeleroI.KopeckaJ.Sarmento-RibeiroA. B.García-ArandaM.De Las RivasJ. (2020). Cancer Immunotherapy Resistance Based on Immune Checkpoints Inhibitors: Targets, Biomarkers, and Remedies. Drug Resist. Updat.53, 100718. Epub 2020/08/01PubMed PMID: 32736034. 10.1016/j.drup.2020.100718PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 48

  PicardE.VerschoorC. P.MaG. W.PawelecG. (2020). Relationships between Immune Landscapes, Genetic Subtypes and Responses to Immunotherapy in Colorectal Cancer. Front. Immunol.11, 369. Epub 2020/03/27PubMed PMID: 32210966; PubMed Central PMCID: PMCPMC7068608. 10.3389/fimmu.2020.00369PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 49

  RibasA.PuzanovI.DummerR.SchadendorfD.HamidO.RobertC.et al (2015). Pembrolizumab versus Investigator-Choice Chemotherapy for Ipilimumab-Refractory Melanoma (KEYNOTE-002): a Randomised, Controlled, Phase 2 Trial. Lancet Oncol.16 (8), 908–918. Epub 2015/06/28PubMed PMID: 26115796. 10.1016/s1470-2045(15)00083-2PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 50

  RobertC.LongG. V.BradyB.DutriauxC.MaioM.MortierL.et al (2015). Nivolumab in Previously Untreated Melanoma withoutBRAFMutation. N. Engl. J. Med.372 (4), 320–330. Epub 2014/11/18PubMed PMID: 25399552. 10.1056/NEJMoa1412082PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 51

  RobertC.ThomasL.BondarenkoI.O'DayS.WeberJ.GarbeC.et al (2011). Ipilimumab Plus Dacarbazine for Previously Untreated Metastatic Melanoma. N. Engl. J. Med.364 (26), 2517–2526. Epub 2011/06/07PubMed PMID: 21639810. 10.1056/NEJMoa1104621PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 52

  Roulleaux DugageM.JonesR. L.TrentJ.ChampiatS.DumontS. (2021). Beyond the Driver Mutation: Immunotherapies in Gastrointestinal Stromal Tumors. Front. Immunol.12, 715727. Epub 2021/09/08PubMed PMID: 34489967; PubMed Central PMCID: PMCPMC8417712. 10.3389/fimmu.2021.715727PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 53

  SinghalR.MittaS. R.DasN. K.KerkS. A.SajjakulnukitP.SolankiS.et al (2021). HIF-2α Activation Potentiates Oxidative Cell Death in Colorectal Cancers by Increasing Cellular Iron. J. Clin. Invest.131 (12), 143691. Epub 2021/04/30PubMed PMID: 33914705; PubMed Central PMCID: PMCPMC8203462. 10.1172/jci143691PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 54

  SingletonD. C.MacannA.WilsonW. R. (2021). Therapeutic Targeting of the Hypoxic Tumour Microenvironment. Nat. Rev. Clin. Oncol.18 (12), 751–772. Epub 2021/07/31PubMed PMID: 34326502. 10.1038/s41571-021-00539-4PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 55

  SkwarskaA.CalderE. D. D.SneddonD.BollandH.OdyniecM. L.MistryI. N.et al (2021). Development and Pre-clinical Testing of a Novel Hypoxia-Activated KDAC Inhibitor. Cell Chem. Biol.28 (9), 1258–1270. e13. Epub 2021/04/29PubMed PMID: 33910023; PubMed Central PMCID: PMCPMC8460716. 10.1016/j.chembiol.2021.04.004PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 56

  SungH.FerlayJ.SiegelR. L.LaversanneM.SoerjomataramI.JemalA.et al (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A Cancer J. Clin.71 (3), 209–249. PubMed PMID: 33538338. 10.3322/caac.21660CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 57

  TazzariM.BergamaschiL.De VitaA.ColliniP.BarisellaM.BertolottiA.et al (2021). Molecular Determinants of Soft Tissue Sarcoma Immunity: Targets for Immune Intervention. Ijms22 (14), 7518. Epub 2021/07/25PubMed PMID: 34299136; PubMed Central PMCID: PMCPMC8303572. 10.3390/ijms22147518PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 58

  VitaleI.ManicG.CoussensL. M.KroemerG.GalluzziL. (2019). Macrophages and Metabolism in the Tumor Microenvironment. Cell Metab.30 (1), 36–50. Epub 2019/07/04PubMed PMID: 31269428. 10.1016/j.cmet.2019.06.001PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 59

  WangS.PangL.LiuZ.MengX. (2021). SERPINE1 Associated with Remodeling of the Tumor Microenvironment in Colon Cancer Progression: a Novel Therapeutic Target. BMC Cancer21 (1), 767. Epub 2021/07/04PubMed PMID: 34215248; PubMed Central PMCID: PMCPMC8254339. 10.1186/s12885-021-08536-7PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 60

  WatanabeT.ShiozawaM.KimuraY.HiroshimaY.HashimotoI.KomoriK.et al (2021). Clinical Significance of Stanniocalcin2 mRNA Expression in Patients with Colorectal Cancer. Anticancer Res.41 (4), 2117–2122. Epub 2021/04/05PubMed PMID: 33813422. 10.21873/anticanres.14983PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 61

  XuF.XuH.LiZ.HuangY.HuangX.LiY.et al (2021). Glycolysis-Based Genes Are Potential Biomarkers in Thyroid Cancer. Front. Oncol.11, 534838. Epub 2021/05/14PubMed PMID: 33981593; PubMed Central PMCID: PMCPMC8107473. 10.3389/fonc.2021.534838PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 62

  XuL.DengC.PangB.ZhangX.LiuW.LiaoG.et al (2018). TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res.78 (23), 6575–6580. Epub 2018/08/30PubMed PMID: 30154154. 10.1158/0008-5472.Can-18-0689PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 63

  XuZ.XuL.LiuL.LiH.JinJ.PengM.et al (2021). A Glycolysis-Related Five-Gene Signature Predicts Biochemical Recurrence-free Survival in Patients with Prostate Adenocarcinoma. Front. Oncol.11, 625452. Epub 2021/05/07PubMed PMID: 33954109; PubMed Central PMCID: PMCPMC8092437. 10.3389/fonc.2021.625452PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 64

  YangF.HeZ.DuanH.ZhangD.LiJ.YangH.et al (2021). Synergistic Immunotherapy of Glioblastoma by Dual Targeting of IL-6 and CD40. Nat. Commun.12 (1), 3424. Epub 2021/06/10PubMed PMID: 34103524. 10.1038/s41467-021-23832-3PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 65

  YangX.WengX.YangY.ZhangM.XiuY.PengW.et al (2021). A Combined Hypoxia and Immune Gene Signature for Predicting Survival and Risk Stratification in Triple-Negative Breast Cancer, Aging13. 19486–19509. (undefined)Epub 2021/08/04PubMed PMID: 34341184. 10.18632/aging.203360PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 66

  YapT. A.ParkesE. E.PengW.MoyersJ. T.CurranM. A.TawbiH. A. (2021). Development of Immunotherapy Combination Strategies in Cancer. Cancer Discov.11 (6), 1368–1397. Epub 2021/04/04PubMed PMID: 33811048; PubMed Central PMCID: PMCPMC8178168. 10.1158/2159-8290.Cd-20-1209PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 67

  YoritaN.YugeR.TakigawaH.OnoA.KuwaiT.KuraokaK.et al (2021). Stromal Reaction Inhibitor and Immune-Checkpoint Inhibitor Combination Therapy Attenuates Excluded-type Colorectal Cancer in a Mouse Model. Cancer Lett.498, 111–120. Epub 2020/11/02PubMed PMID: 33129954. 10.1016/j.canlet.2020.10.041PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 68

  YouL.WuW.WangX.FangL.AdamV.NepovimovaE.et al (2021). The Role of Hypoxia‐inducible Factor 1 in Tumor Immune Evasion. Med. Res. Rev.41 (3), 1622–1643. Epub 2020/12/12PubMed PMID: 33305856. 10.1002/med.21771PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 69

  ZhaoE.ChenS.DangY. (2020). Development and External Validation of a Novel Immune Checkpoint-Related Gene Signature for Prediction of Overall Survival in Hepatocellular Carcinoma. Front. Mol. Biosci.7, 620765. Epub 2021/02/09PubMed PMID: 33553243; PubMed Central PMCID: PMCPMC7859359. 10.3389/fmolb.2020.620765PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar

• 70

  ZouM.-Z.LiuW.-L.LiC.-X.ZhengD.-W.ZengJ.-Y.GaoF.et al (2018). A Multifunctional Biomimetic Nanoplatform for Relieving Hypoxia to Enhance Chemotherapy and Inhibit the PD-1/pd-L1 Axis. Small14 (28), 1801120. Epub 2018/06/09PubMed PMID: 29882235. 10.1002/smll.201801120PubMed Abstract | CrossRef Full Text | Google Scholar

  • CrossRef
  • Google Scholar