A Hypoxia-Related lncRNA Signature Correlates with Survival and Tumor Microenvironment in Colorectal Cancer

Zhong, Xinyang; He, Xuefeng; Wang, Yaxian; Hu, Zijuan; Yu, Deshui; Huang, Huixia; Zhao, Senlin; Wei, Ping; Li, Dawei

doi:https://doi.org/10.1155/2022/9935705

Journal of Immunology Research

On this page

Abstract Introduction Methods and Materials Results Discussion Conclusion Data Availability Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

The Tumor Microenvironment and Immunity 2021

View this Special Issue

Research Article | Open Access

Volume 2022 | Article ID 9935705 | https://doi.org/10.1155/2022/9935705

A Hypoxia-Related lncRNA Signature Correlates with Survival and Tumor Microenvironment in Colorectal Cancer

Xinyang Zhong,^1,2Xuefeng He,^1,2Yaxian Wang,^1,2Zijuan Hu,^2,3,4,5Deshui Yu,⁶Huixia Huang,^2,3,4,5Senlin Zhao,^1,2Ping Wei ,^2,3,4,5and Dawei Li^1,2

Academic Editor: Xue-Li Zhang

Received21 Mar 2022

Revised25 Apr 2022

Accepted03 May 2022

Published08 Jul 2022

Abstract

The hypoxic tumor microenvironment and long noncoding RNAs (lncRNAs) are pivotal in cancer progression and correlate with the survival outcome of patients. However, the role of hypoxia-related lncRNAs (HRLs) in colorectal cancer (CRC) development remains largely unknown. Herein, we developed a hypoxia-related lncRNA signature to predict patients’ survival and immune infiltration. The RNA-sequencing data of 500 CRC patients were obtained from The Cancer Genome Atlas (TCGA) dataset, and HRLs were selected using Pearson’s analysis. Next, the Cox regression analysis was applied to construct a risk signature consisting of 9 HRLs. This signature could predict the overall survival (OS) of CRC patients with high accuracy in training, validation, and entire cohort. This signature was an independent risk factor and exerted predictive ability in different subgroups. Functional analysis revealed different molecular features between high- and low-risk groups. A series of drugs including cisplatin showed different sensitivities between the two groups. The expression pattern of immune checkpoints was also distinct between the two clusters in this model. Furthermore, the high-risk group had higher immune, stromal, and ESTIMATE score and a more repressive immune microenvironment than the low-risk group. Moreover, MYOSLID, one of the lncRNAs in this signature, could significantly regulate the proliferation, invasion, and metastasis of CRC.

1. Introduction

According to the latest epidemiological statistics, colorectal cancer (CRC) is an important factor that threatens the survival of mankind [1]. In 2022, it has been predicted that CRC will still be the third most common tumor and the second leading cause of cancer-related death in the USA [1]. Although a high incidence and mortality rate exist in CRC, effective treatments are still limited due to the complexity of tumor biology and the tumor heterogeneity [2]. Tumor microenvironment (TME) is a complex and heterogenous population that is composed of cancer cells, stromal cells such as cancer-associated fibroblasts, immune cells, cytokine networks, and oxygen [3]. Hypoxia is a common event in large amounts of solid tumors due to high consumption of oxygen and disordered intratumor blood vessels [4, 5]. Thus, cancer cells undergo metabolic reprogramming towards a more glycolytic phenotype and adjust to such a hypoxic environment, also known as the Warburg effect [6]. Enhanced glycolytic activity in cancer cells helps them obtain enough intermediate to synthesize biomacromolecules and promotes lactate secretion, which leads to dysregulation of extracellular pH [4, 6]. Hypoxia-induced pathophysiological processes could help cancer cells gain a survival advantage over other cells which leads to an immunosuppressive tumor microenvironment, an unsatisfactory chemotherapeutic effect, and poor patient survival [4, 7].

Long noncoding RNA (lncRNA) is widely involved in cancer progression but has no potential to code for protein [8]. Emerging studies have proved that extracellular oxygen levels regulate the expression of lncRNA [9]. Moreover, the expression pattern of lncRNAs varies under hypoxia and normoxia [9]. Huan et al. observed that the expression of lncRNA LUCAT1 in CRC was upregulated under hypoxia, which can promote the growth and drug resistance of CRC [10]. In addition, it has been demonstrated that dysregulated lncRNA levels in cancer can regulate some key regulators related to hypoxia such as HIF1A [11]. For example, lncRNA SNHG11 is upregulated and correlated with poor prognosis in CRC patients. Mechanically, SNHG11 can bind to pVHL and protect HIF1A from the degradation, accumulating HIF1A in the cancer cells [12]. Moreover, the accumulation of HIF1A is essential for tumor metastasis because it is involved in a lot of oncogenic processes such as epithelial-mesenchymal transition (EMT) and metabolic reprogramming [13]. Therefore, it has been established that hypoxia and lncRNA could regulate each other and form a complex network which has a strong impact on a series of cancer hallmarks [11].

Both lncRNAs and hypoxic TME contribute to tumor malignancy. Hence, discerning a hypoxia-related lncRNA (HRL) panel and establishing a relevant model could evaluate patient’s clinical prognosis [14]. In addition, it has been revealed that targeting certain HRL can be an effective method in cancer treatment [10, 12]. Recently, several studies have explored the possibility of developing a HRL signature to predict patients’ survival and tumor microenvironment in multiple cancer types [14–16]. Nevertheless, application of the HRL signature in CRC is still unknown. Therefore, this research constructed a HRL signature in CRC using The Cancer Genome Atlas (TCGA) database. This signature could not only offer prognostic value but also describe the TME features and estimate clinical treatments of each patient.

2. Methods and Materials

2.1. Data Acquisition and Processing

The expression matrix (FPKM) of CRC patients was downloaded from TCGA database (https://portal.gdc.cancer.gov/). Then, the FPKM value was transformed into the TPM value. The study included a total of 500 CRC patients with the following criteria: (1) overall survival (OS) time is more than 30 days. (2) Clinical data includes T, M, N, and AJCC stages, age, and gender. R package “caret” was applied to divide 500 CRC patients into training and validation group with a ratio of 7 : 3. Subsequently, 31 hypoxia-related genes were collected based on previous studies and listed in the Supplementary Table S1 [17, 18].

2.2. Identification of Hypoxia-Related lncRNAs

A total of 14048 RNAs were considered as lncRNA according to the annotation file from GENCODE (https://www.gencodegenes.org). Moreover, package “limma” was used to discover differently expressed lncRNAs ( and ) between tumorous and nontumorous samples. Further, differently expressed lncRNAs whose expression was correlated with the expression of at least 1 hypoxia-related gene ( and ) were regarded as HRLs.

2.3. Construction of the Hypoxia-Related lncRNA Model

Univariate Cox regression analysis was applied to select HRLs related to the overall survival (OS) of CRC patients (). Next, the associated lncRNAs were further selected by least absolute shrinkage and selection operator (LASSO) Cox regression analysis using the “glmnet” package. Finally, several lncRNAs and their relevant coefficients were confirmed which were applied to establish a HRL model. The risk .

2.4. Prognostic Evaluation

The Kaplan-Meier survival analysis was applied to find out the survival difference between different subgroups. The sensitivity and specificity of HRL signature were calculated by ROC analysis. Further, univariate and multivariate Cox regression analyses were used to validate whether the HRL signature was an independent risk factor.

2.5. Nomogram

The 1-, 3-, 5-, and 7-year survival probabilities of CRC patients in the training and total cohort were predicted by nomogram. Nomogram construction included risk score, AJCC stage, gender, and age for accurate prediction. Moreover, the discrimination and accuracy of the nomograms were evaluated by concordance index (C-index) and calibration curves.

2.6. Functional Enrichment Analysis

The R package “limma” was used to select differently expressed genes between high- and low-score groups for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Subsequently, the gene symbols of these genes were converted into ENTREZ ID. The gene set “h.all.v7.4.entrez.gmt” was downloaded from the msigDB database (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#H) for Gene Set Enrichment Analysis (GSEA) studies. Similarly, the gene symbols were converted into ENTREZ ID and sorted by descending log fold changes (logFC).

2.7. Drug Sensitivity and Immunotherapy Response Prediction

Using package “pRRophetic,” predicted drug sensitivity of CRC patients to certain drug was calculated [19]. By entering patient’s RNA-sequencing data, this package could calculate relevant IC₅₀ (half-maximal inhibitory concentration) value based on its built-in dataset [19]. To predict immunotherapy response, TIDE (Tumor Immune Dysfunction and Exclusion) algorithm was applied [20]. The predicted efficiencies of the immune checkpoint blockade (ICB) treatment of each CRC patient could be obtained after feeding their RNA-Seq data on the official website (http://tide.dfci.harvard.edu/).

2.8. Mutation Analysis

The copy number variation (CNV) data of the TCGA-COAD and TCGA-READ cohort was downloaded from TCGA database (https://portal.gdc.cancer.gov/). Further, to analyze the CNV data of high- and low-score groups, R package “maftools” was used [21]. The mutation landscape of both groups was visualized through waterfall plot. Tumor mutation burden (TMB) and differently mutated genes were further explored to reveal the intrinsic mechanism that leads to the biological difference between the two groups.

2.9. Immune Infiltration Evaluation

Immune, stromal, and ESTIMATE scores of CRC patients were calculated using ESTIMATE algorithm [22]. Further, the CIBSERSORT algorithm analyzed the 22 immune cells proportion of each sample by simply inputting TPM values of the sequencing data [23]. TIMER2.0 is a user-friendly web server that provides a rounded analysis of tumor-infiltrating immune cells in the TCGA cohorts [24]. The immune infiltration data including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, and EPIC was downloaded from TIMER2.0 (http://timer.cistrome.org/). Also, the expression of immune checkpoints and cytokines was analyzed in two groups.

2.10. ssGESA

Single-sample gene set enrichment (ssGSEA) can calculate the score of certain cell or pathway at the level of a single sample [25]. A series of immune gene sets were retrieved from a previous study and are listed in Supplementary Table S2 [26]. Further, the scores of these gene sets were calculated using “GSVA” package.

2.11. In Vitro Experiments

To further explore the potential therapeutic targets in CRC, a series of functional experiments were performed and the detailed methods are introduced in Supplementary methods.

2.12. Statistical Analysis

R software performed the majority of the bioinformatic analysis. Student’s -test was used to compare the difference between the two groups of normally distributed variables. The Wilcoxon rank-sum test was applied for variables that were not normally distributed. The chi-square test could compare the differences between the two groups of categorical variables. The Pearson regression analysis revealed the correlation between the two variables. The two-sided value of less than 0.05 was considered statistically significant.

3. Results

3.1. Construction of a Hypoxia-Related lncRNA Signature

Figure 1(a) presented the overall design of the research. We obtained 159 differently expressed hypoxia-related lncRNAs between normal and tumorous samples. Among the 159 lncRNAs, 27 lncRNAs were correlated with the OS of CRC patients (Figure 1(b)). Package “caret” randomly separated CRC samples into a training () and a validation () cohort with a ratio of 7 : 3. There was no significant survival difference between the two cohorts (Figure S1A). We used training group to construct a HRL signature and test and entire cohort () to verify the reliability of this model. LASSO regression analysis further removed 18 genes according to the value of lambda.min, and the remaining 9 genes with relevant regression coefficients were obtained to construct the prognostic model (Figures 1(c)–1(e)). As is shown in Figure 1(f), all of the 9 lncRNAs had positive correlations () with at least one hypoxia-related gene. The expression of LINC02257, LINC02188, MYOSLID, C6orf223, MYG1-AS1, and Lnc-SKA2-1 was upregulated in tumor samples, while the expression of LINC00702, LOC100129434, and LINC01915 was downregulated in tumor samples (Figure S1B). The TPM values and coefficients of the 9 lncRNAs were applied to generate the risk score: risk . According to the median value of each group, the patients were divided into high- and low-score group.

(a)

(b)

(c)

(d)

(e)

(f)

3.2. HRL Signature Correlated with Prognosis in CRC

The Kaplan-Meier survival analysis showed that the high-score group exhibited a much shorter OS time than the low-score group in training, testing, and entire group (Figures 2(a)–2(c)). Then, the calculated the AUC value of 1-, 3-, 5-, and 7-year survival in all 3 cohorts confirmed the outstanding predictive ability of HRL model (Figures 2(d)–2(f)). The 5-year survival rate is commonly regarded as an important indicator to evaluate the outcomes of tumor patients. The AUCs for 5-year OS prediction in training, validation, and entire cohort were 0.8, 0.63, and 0.76, respectively, providing an effective basis for clinicians to judge the prognosis of patients (Figures 2(d)–2(f)). The distribution of survival time, survival status, and risk score were visualized by risk curves and scatter plots (Figures 2(g)–2(i)). In addition, the expression pattern of the 9 HRLs between high- and low-score groups was also different (Figures 2(j)–2(l)). Overall, the results demonstrated that the HRL signature exerts outstanding prognostic ability in CRC patients.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

3.3. Clinical Performance of the lncRNA Signature

Univariate Cox regression analysis revealed the values of lncRNA signature, age, and AJCC stage to be less than 0.05, indicating their association with patients’ OS (Figures 3(a)–3(c)). Multivariate Cox regression analysis discovered that the values of HRL signature were all less than 0.01 and the values of HR were all more than 2, demonstrating that the HRL signature is an independent risk factor in all three cohorts () (Figures 3(d)–3(f)). The clinicopathological characteristics of both the training and validation sets were summarized in Table 1. We aimed to find out whether the HRL signature showed predictive abilities in different subgroups. Interestingly, the lncRNA signature could still exert predictive function in the majority of the subgroups (11 out of 12) (Figures 3(g)–3(r)). The relationship among clinical factors, gene expression, and risk level were shown using a heat map (Figure S2A). Moreover, the risk score in T3-4, N1-2, M1, and III-IV stage was significantly higher than in T1-2, N0, M0, and stages I-II (Figure S2B).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

(r)

3.4. Construction of a Nomogram

Risk score, AJCC stage, gender, and age were used to construct the nomograms for predicting the 1-, 3-, 5-, and 7-year OS in both training cohort and entire cohort (Figures 4(a) and 4(b)). In the study, the patients were given a nomogram-based risk score and were estimated for survival probability. In addition, 1-, 3-, 5-, and 7-year calibration curves were drawn for both cohorts, which indicated that the nomogram showed high accuracy and stability (Figures 4(c)–4(j)). The concordance index of the nomogram in training and entire cohort was 0.78 and 0.76, respectively.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

3.5. Gene Enrichment Analysis

A series of functional analyses including GO, KEGG, and GESA were performed to determine the potential pathway or process that leads to the different clinical outcomes between the two groups. The R package “limma” was applied to select differently expressed genes between high- and low-risk subgroups, and then, these genes were subjected to GO and KEGG analyses. Representative GO terms included biological process (BP), molecular function (MF), and cellular component (CC). The BP included positive regulation of cytosolic calcium ion concentration, extracellular structure organization, extracellular matrix organization, calcium ion transport, and BMP signaling pathway (Figure 5(a)); the MF included signaling receptor activator activity, receptor ligand activity, ion channel activity, extracellular matrix structural constituent, and cell adhesion molecule binding (Figure 5(b)); and the CC included transporter complex, ion channel complex, glutamatergic synapse, collagen trimer, and collagen-containing extracellular matrix (Figure 5(c)). Moreover, the KEGG results indicated that different signaling such as ECM-receptor interaction, pentose and glucuronate interconversions, chemical carcinogenesis-DNA adducts, focal adhesion, and bile secretion might be involved in the different prognosis between the two groups (Figure 5(d)). Gene sets from “h.all.v7.4.entrez.gmt” were used to analyze their activities in each group. Next, the 5 most enriched pathways in each group were selected after setting , values <0.25, and abs as the criteria to judge the enriched pathways. As is shown in Figures 5(e) and 5(f), pathways such as HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION and HALLMARK_HYPOXIA were enriched in the high-risk group, while HALLMARK_OXIDATIVE_PHOSPHORYLATION and HALLMARK_MYC_TARGETS were enriched in the low-risk group.

(a)

(b)

(c)

(d)

(e)

(f)

3.6. HRL Signature Correlated with the Therapeutic Effects of Chemotherapy and Immunotherapy

The IC₅₀ of several common drugs was calculated. We found the IC₅₀ of masitinib, phenformin, and ruxolitinib was higher in the high-risk group, indicating that the high-risk group patients tend to be more resistant to these drugs (Figure 6(a)). On the contrary, the high-risk group might benefit from pazopanib, foretinib, and cisplatin, a widely used drug in CRC treatment (Figure 6(a)). Further, the mutation rate of the top 20 driver genes between two groups was visualized through waterfall plots (Figure S3A–B). A higher mutation rate of MUC16, NIPBL, PTPRK, MYO9A, LPA, TRRAP, APOB, and COL7A1 was found in the high-risk group (Figure S3C). Interestingly, the mutation rate of MUC16 in the high-risk group was 35% while in the low-risk group was 21% (Figure S3A–B). Moreover, the mutation of MUC16 in solid tumors has been reported to be correlated with a higher tumor mutational burden [27]. Immunotherapy is emerging as a promising method in cancer treatment. We discovered a significantly higher expression of CTLA4, PDCD1/PD1, LAG3, HAVCR2/TIM3, and PDCD1LG2/PD-L2 in the high-risk group (Figure 6(b)). Additionally, patients in the low-risk group showed a much lower TIDE score than its counterpart (), suggesting that immunotherapy had better therapeutic effects in these patients (Figure 6(c)). MSI status and TMB are two important hallmarks in predicting the efficiency of ICB therapy [28]. The results demonstrated that the high-risk group occupied a higher proportion of MSI-H patients () (Figure 6(d)). What is more, the patients with MSI-L/MSS and low scores had better OS compared with those with MSI-L/MSS and high scores (Figure 6(e)). Similarly, the low-score MSI-H patients had better OS than the high score MSI-H patients (Figure 6(e)). However, the TMB analysis showed no significant difference between the two groups (Figure 6(f)), although the classification based on both TMB and risk score demonstrated that the risk score could predict patients’ OS regardless of their TMB status (Figure 6(g)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 6

Treatment prediction. (a) Box plots that showed the IC₅₀ of 6 drugs in high- and low-score group. (b) Box plots that showed the expression of immune checkpoints in high- and low-score group. (c) TIDE prediction score in high- and low-score group. (d) Bar chart that showed the proportion of MSS/MSI-L and MSI-H patients in high- and low-score group. (e) The Kaplan-Meier plots that compared the OS of MSS/MSI-L+high-risk, MSS/MSI-L+low-risk, MSI-H+high-risk, and MSI-H+low-risk group. (f) TMB levels of high- and low-score group. (e) The Kaplan-Meier plots that compared the OS of TMB-H+high-risk, TMB-H+low-risk, TMB-L+high-risk, and TMB-L+low-risk group.

3.7. Hypoxia-Related lncRNA Signature Correlated with Immune Infiltration

Exploring the infiltration pattern of immune cells could offer prognostic value and guide immunotherapy in CRC [28, 29]. Thus, a series of immune-related algorithms were applied to analyze the immune infiltration status of 500 CRC patients in the study. ESTIMATE algorithm revealed that the high-risk group had a higher stromal, immune, and ESTIMATE score (Figures 7(a)–7(c)). CIBERSORT results showed that the proportion of macrophages M0, mast cells activated, and T cells regulatory (Tregs) was higher, while the proportion of dendritic cells activated, dendritic cells resting, mast cells resting, plasma cells, T cells CD4 memory resting, and T cells CD4 memory activated was lower in the high-risk group (Figure 7(d)). Pearson analysis showed that the risk score had strong positive correlation with M0 macrophage and a strong negative correlation with T cells CD4 memory resting (Figures 7(e) and 7(f)). Further, the relationship between immune/stromal cells and HRL signature was explored using algorithms including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, and EPIC. We observed that the relationship between macrophage, dendritic cell, CD4 T cell, and risk score in CIBERSORT was consistent with the majority of these algorithms (Figure 7(g)). Surprisingly, the risk score had strong positive correlation with cancer-associated fibroblasts (CAFs) ( in XCELL, in MCPCOUNTER, and in EPIC) (Figure 7(g)). Next, the expression patterns of cytokines were analyzed as they play an important role in CRC immunology. It was observed that most of the protumor cytokines had positive correlation with the risk score (Figure 7(h)). Scores of 18 immunotherapy response-related pathways were calculated by ssGESA. Surprisingly, an inverse correlation was observed between risk score and the scores of the pathways whose value was less than 0.05 (Figure 7(i)).

Figure 7

Tumor microenvironment estimation. (a–c) Beeswarm plots that compare the stromal score, immune score, and ESTIMATE score between high- and low-risk groups. (d) Box plots that showed the proportion of 22 immune cells in high- and low-score group. (e) Scatter plot that showed the correlation between T cell CD4 memory resting and risk score. (f) Scatter plot that showed the correlation between macrophage M0 and risk score. (g) Relationship between risk score and the proportion of stromal cells. (h) Correlation heat map that showed the relationship between risk score and the expression of cytokines. (i) Correlation heat map that showed the relationship between risk score and the activities of immunotherapy response-related pathways.

3.8. MYOSLID Promotes the Progression of CRC

There have been limited studies that discuss the oncogenic roles of MYOSLID. The previous research has indicated that MYOSLID was upregulated in the head and neck squamous cell carcinoma and contributed to tumor metastasis [30]. Similarly, MYOSLID could enhance the degree of malignancy by promoting the expression of MCL1 in gastric cancer [31]. However, no study has been reported that discusses the role of MYOSLID in CRC. We explored the relationship between patients’ OS and expression of 9 HRLs and found that high expression of MYOSLID indicated a poor prognosis in TCGA CRC cohort (Figure 8(a) and Figure S4). Also, higher expression of MYOSLID can be found in the advanced T stage, N1-2 stage, and AJCC stage in TCGA CRC cohort (Figure 8(b)). In addition, our study found a strong positive correlation between MYOSLID and EMT-related genes including VIM, SNAI2, CDH2, and TGFB1 (Figure 8(c)). Similar correlation was also found between HIF1A and MYOSLID (Figure 8(c)). Therefore, we speculated that MYOSLID might also promote the progression of CRC. We examined the expression of MYOSLID in HCT116, HCT15, HCT8, DLD1, SW480, RKO, LOVO, CaCO2, and NCM460 and discovered that the expression of MYOSLID was highest in HCT15 (Figure 8(d)). Additionally, hypoxic environment can induce MYOSLID expression (Figure 8(e)). Two siRNAs were transfected into HCT15 to knock down MYOSLID (Figure 8(f)). It was observed that the knockdown of MYOSLID significantly reduced the proliferative ability of HCT15 (Figure 8(g)). Furthermore, results from the transwell assay and wound healing assay indicated that repressing MYOSLID expression could also inhibit invasion and migration in CRC (Figures 8(h)–8(k)). Therefore, we concluded that MYOSLID could function as an oncogene and could become a potential therapeutic target in CRC.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 8

The role of MYOSLID in CRC. (a) The Kaplan-Meier plot that compared the OS between high- and low-MYOSLID groups. (b) Box plots that compared the expression of MYOSLID in subgroups with different clinical features. (c) Relationship between risk score and VIM, SNAI2, CDH2, TGFB1, and HIF1A. (d) The expression of MYOSLID (RNA level) in different cell lines. (e) The expression of MYOSLID under hypoxia. (f) The expression of MYOSLID (RNA level) before and after knocking down MYOSLID. (g) The OD450 value of HCT15 before and after knocking down MYOSLID. (h–k) Relative cell number in transwell assay. (i and j) Relative migration distance in wound healing assay.

4. Discussion

Hypoxia has been regarded as a hallmark of the tumor microenvironment and exists in most of the solid tumors [4]. It has been studied that lack of oxygen could remodel tumor microenvironment by promoting abnormal angiogenesis and reducing extracellular pH levels and nutrient availability, resulting in immune suppression, drug resistance, and tumor progression [4, 32]. Herein, the 9 hypoxia-related lncRNAs were applied to construct a HRL signature. The HRL signature exhibited excellent predictive ability in predicting patients’ OS. Additionally, the signature divided patients into two clusters, which showed distinctively different immune checkpoint expression and immune cell infiltration patterns. Finally, we experimentally validated the oncogenic function of MYOSLID in CRC.

The past decade has widely studied the function of lncRNA in CRC development, and the application of lncRNA as a prognostic model has been reported in a number of studies [33]. It has been demonstrated that hypoxia could induce the aberrant expression of lncRNAs, such as H19, HOTAIR, and NEAT1 to regulate tumor biology [34]. Moreover, lncRNA could upregulate the expression of HIF by increasing HIF-1α expression at both transcription and posttranscription levels [35]. Recently, a series of studies revealed that several hypoxia-related lncRNA signatures could provide predictive value [14, 16, 36]. In bladder cancer, a lncRNA signature consisting of 4 hypoxia-related lncRNAs was constructed and was found to exhibit excellent predictive performance [16]. Hence, it is pivotal to establish a HRL signature in CRC. Using bioinformatic analyses, we obtained 9 HRLs which included LINC02257, LINC02188, MYOSLID, C6orf223, MYG1-AS1, Lnc-SKA2-1, LINC00702, LOC100129434, and LINC01915 for further analysis. A HRL signature was established to calculate the risk scores of each patient. The Kaplan-Meier analysis and ROC curve revealed that the HRL signature showed potent accuracy in survival prediction. The signature was also an independent risk factor and could predict the OS in different subgroups. Hence, this novel hypoxia-related lncRNA signature could offer predictive value in CRC.

In CRC, chemotherapy is still the mainstay of treatment other than surgery [2]. Surprisingly, patients in the high-risk group were predicted to be more sensitive to cisplatin, a widely applied chemotherapy drug in CRC. Immune checkpoint therapy can offer lasting therapeutic effects and successfully improves the OS of patients [37, 38]. Nevertheless, the vast majority of CRC patients showed resistance to immune therapy, which greatly restricts their application in the clinic [39]. It has been demonstrated that TMB, mismatch repair deficiency, expression of immune checkpoint, and tumor-infiltrating immune cells are predictive biomarkers for therapeutic efficiency of ICB [28]. We found that patients in the high-score group had high mRNA level of CTLA4, LAG3, PDCD1, HAVCR2, and PDCD1LG2. Further, a larger population of MSI-H patients can be found in the high-score group. However, the TIDE algorithm predicted that patients in the low-score group can receive more benefits from ICB, which might be the result of a higher proportion of antitumor immune cell population in the low-score group.

Emerging evidence has demonstrated that the development and prognosis of CRC are tightly closed to TME. The proportion of a series of antitumor immune cells in CRC was significantly lower than in normal colorectal mucosa [29]. In addition, it has been observed that the suppressive myeloid cells were dominant in immunosuppressive CRC TME [29]. Studies have revealed that hypoxia could regulate the biological function of nearby stromal cells to shape an immunosuppressive microenvironment [4]. HIF-1 signaling has been reported to help cancer cells escape immune attack by inducing the expression of PD-L1 [40]. Hypoxia could also change the function of myeloid-derived suppressor cells (MDSCs) and promote their differentiation toward tumor-associated macrophages in the TME [41]. CAFs also facilitate the formation of an immunosuppressive TME [42]. The expression of several immunosuppressive cytokines including IL6, IL10, and TGF-β that are secreted by CAFs was also seen to be upregulated in hypoxia, inhibiting T cell-mediated cytotoxicity to suppress immune activation [43]. Our results revealed that a higher stromal score existed in the high-risk group, which could be a result of the enrichment of CAFs. Additionally, the HRL score was positively correlated with the proportion of CAF in XCELL, MCPCOUNTER, and EPIC. In our study, the HRL signature divided CRC patients into two clusters with distinct immune infiltration patterns. The result that the high-risk score group had a high infiltration of immune cells seemed to be perplexing. In fact, upregulation of immunosuppressive immune cells could contribute to the progression of cancer and indicates poor survival outcomes. Our study found a higher proportion of macrophages M0 and T cells regulatory (Tregs) in the high-risk group and a higher proportion of dendritic cells activated, dendritic cells resting, mast cells resting, plasma cells, T cells CD4 memory resting, and T cells CD4 memory activated in the low-risk group. Furthermore, the HRL score was negatively correlated with immunotherapy response-related pathways and was positively correlated with several protumor cytokines, which partly explained the poor response of the high-risk group to ICB therapy in TIDE algorithm.

It has been observed that lncRNAs selected using bioinformatic analysis could be more likely to be potential therapeutic targets for cancer treatment [14]. LINC02257, an enhancer RNA with DUSP10 as the target gene, has been reported to be upregulated and associated with the survival outcomes in multiple kinds of cancers [44]. MYG1-AS1 was regarded as a hypoxia-related lncRNA in hepatocellular carcinoma [45]. The role of LINC00702 in tumor development has been studied in several researches. LINC00702 functions as a tumor-suppressor gene by promoting PTEN expression in colorectal cancer and nonsmall cell lung cancer (NSCLC) while as an oncogene in malignant meningioma and ovarian cancer [46–49]. LINC01915 has been reported to inhibit the formation of CAFs through miR-92a-3p/KLF4/CH25H axis, which was consistent with our analysis that LINC01915 was downregulated in CRC compared with normal tissues [50]. MYOSLID has been reported to be a hypoxia- and autophagy-related lncRNA in HNSCC [51, 52]. Furthermore, the role of MYOSLID in HNSCC and gastric cancer has been validated using a series of in vitro and in vivo experiments [30, 31]. Our research focused on exploring the function of MYOSLID in CRC. The knockdown effects of MYOSLID were studied by applying a series of functional assays. The results revealed that MYOSLID could promote the progression of CRC by enhancing proliferation, invasion, and migration of CRC cells. Hence, MYOSLID could be a potential target in CRC therapy.

5. Conclusion

In summary, we developed a HRL signature that could predict the survival outcomes of CRC patients. Additionally, the HRL signature was correlated with drug sensitivity, immune checkpoint expression, and immune infiltration of CRC. Finally, we observed that MYOSLID could affect the biological function of CRC. These indicated that hypoxia-related lncRNA signature could be a novel biomarker in CRC diagnosis and promising therapeutic target in cancer treatment.

Data Availability

In this study, all the data can be found in the following: (1) https://www.cancergenome.nih.gov/, (2) http://www.gsea-msigdb.org/gsea/msigdb/, (3) http://timer.comp-genomics.org/, (4) https://www.gencodegenes.org, and (5) http://tide.dfci.harvard.edu/.

Conflicts of Interest

The authors declare that they have no competing interests.

Authors’ Contributions

Xinyang Zhong, Xuefeng He, and Yaxian Wang contributed equally to this work. Dawei Li and Ping Wei designed the whole study. Xinyang Zhong wrote the whole draft. Xinyang Zhong, Xuefeng He, and Yaxian Wang performed the experiments and bioinformatic analyses. Zijuan Hu, Deshui Yu, and Huixia Huang made the figures in this study. Senlin Zhao, Dawei Li, and Ping Wei checked and revised the manuscript. All authors approved this manuscript.

Acknowledgments

This research was supported by grants from the National Science Foundation of China (81772583, 81972293, and 81972185) and Shanghai Rising-Star Program (19QA1402200).

Supplementary Materials

Supplementary Figure S1: (A) the Kaplan-Meier analysis of training cohort and validation cohort. (B) Heat map of expression matrix of 9 lncRNAs in tumor and normal samples. Supplementary Figure S2: (A) heat map that described the correlation among risk level, clinical features, and hypoxia-related lncRNAs. (B) Risk score comparison in different T, M, N, and AJCC stage. Supplementary Figure S3: (A and B) waterfall plots of the top 20 mutated genes in two groups. (C) Differentially mutated genes between high- and low-risk groups. Supplementary Figure S4: the Kaplan-Meier analysis of LINC02257, LINC02188, LINC00702, LOC100129434, LINC01915, C6orf223, MYG1-AS1, and Lnc-SKA2-1. Supplementary Table S1: 31 hypoxia-related genes. Supplementary Table S2: gene lists used in ssGSEA. Supplementary Table S3: sequences of siRNA and primers used in this study. Supplementary Methods: in vitro experiments in this study. (Supplementary Materials)

References

R. Siegel, K. D. Miller, H. E. Fuchs, and A. Jemal, “Cancer statistics, 2022,” CA: a Cancer Journal for Clinicians, vol. 72, no. 1, pp. 7–33, 2022.
View at: Publisher Site | Google Scholar
E. Dekker, P. J. Tanis, J. L. A. Vleugels, P. M. Kasi, and M. B. Wallace, “Colorectal cancer,” Lancet (London, England), vol. 394, article 10207, pp. 1467–1480, 2019.
View at: Publisher Site | Google Scholar
D. Hanahan and R. Weinberg, “Hallmarks of cancer: the next generation,” Cell, vol. 144, no. 5, pp. 646–674, 2011.
View at: Publisher Site | Google Scholar
X. Jing, F. Yang, C. Shao et al., “Role of hypoxia in cancer therapy by regulating the tumor microenvironment,” Molecular Cancer, vol. 18, no. 1, p. 157, 2019.
View at: Publisher Site | Google Scholar
A. Patel and S. Sant, “Hypoxic tumor microenvironment: opportunities to develop targeted therapies,” Biotechnology Advances, vol. 34, no. 5, pp. 803–812, 2016.
View at: Publisher Site | Google Scholar
M. Liberti and J. Locasale, “The Warburg effect: how does it benefit cancer cells?” Trends in Biochemical Sciences, vol. 41, no. 3, pp. 211–218, 2016.
View at: Publisher Site | Google Scholar
J. Walsh, A. Lebedev, E. Aten, K. Madsen, L. Marciano, and H. C. Kolb, “The clinical importance of assessing tumor hypoxia: relationship of tumor hypoxia to prognosis and therapeutic opportunities,” Antioxidants & Redox Signaling, vol. 21, no. 10, pp. 1516–1554, 2014.
View at: Publisher Site | Google Scholar
S. Liu, H. X. Dang, D. A. Lim, F. Y. Feng, and C. A. Maher, “Long noncoding RNAs in cancer metastasis,” Nature Reviews. Cancer, vol. 21, no. 7, pp. 446–460, 2021.
View at: Publisher Site | Google Scholar
H. Choudhry, A. Harris, and A. McIntyre, “The tumour hypoxia induced non-coding transcriptome,” Molecular Aspects of Medicine, vol. 47-48, pp. 35–53, 2016.
View at: Publisher Site | Google Scholar
L. Huan, T. Guo, Y. Wu et al., “Hypoxia induced LUCAT1/PTBP1 axis modulates cancer cell viability and chemotherapy response,” Molecular Cancer, vol. 19, no. 1, p. 11, 2020.
View at: Publisher Site | Google Scholar
X. Wang, D. Zhao, H. Xie, and Y. Hu, “Interplay of long non-coding RNAs and HIF-1α: a new dimension to understanding hypoxia-regulated tumor growth and metastasis,” Cancer Letters, vol. 499, pp. 49–59, 2021.
View at: Publisher Site | Google Scholar
L. Xu, L. Huan, T. Guo et al., “LncRNA SNHG11 facilitates tumor metastasis by interacting with and stabilizing HIF-1α,” Oncogene, vol. 39, no. 46, pp. 7005–7018, 2020.
View at: Publisher Site | Google Scholar
S. Unwith, H. Zhao, L. Hennah, and D. Ma, “The potential role of HIF on tumour progression and dissemination,” International Journal of Cancer, vol. 136, no. 11, pp. 2491–2503, 2015.
View at: Publisher Site | Google Scholar
S. Xu, L. Tang, Z. Liu, C. Luo, and Q. Cheng, “Hypoxia-related lncRNA correlates with prognosis and immune microenvironment in lower-grade glioma,” Frontiers in Immunology, vol. 12, article 731048, 2021.
View at: Publisher Site | Google Scholar
P. Tang, W. Qu, T. Wang et al., “Identifying a hypoxia-related long non-coding RNAs signature to improve the prediction of prognosis and immunotherapy response in hepatocellular carcinoma,” Frontiers in Genetics, vol. 12, article 785185, 2021.
View at: Publisher Site | Google Scholar
X. Chen, Y. Zhang, F. Wang et al., “A novel assessment model based on molecular subtypes of hypoxia-related LncRNAs for prognosis of bladder cancer,” Frontiers in Cell and Developmental Biology, vol. 9, article 718991, 2021.
View at: Publisher Site | Google Scholar
A. Eustace, N. Mani, P. N. Span et al., “A 26-gene hypoxia signature predicts benefit from hypoxia-modifying therapy in laryngeal cancer but not bladder cancer,” Clinical Cancer Research : An Official Journal of the American Association for Cancer Research, vol. 19, no. 17, pp. 4879–4888, 2013.
View at: Publisher Site | Google Scholar
Y. Ye, Q. Hu, H. Chen et al., “Characterization of hypoxia-associated molecular features to aid hypoxia-targeted therapy,” Nature Metabolism, vol. 1, no. 4, pp. 431–444, 2019.
View at: Publisher Site | Google Scholar
P. Geeleher, N. Cox, and R. Huang, “pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels,” PLoS One, vol. 9, no. 9, article e107468, 2014.
View at: Publisher Site | Google Scholar
P. Jiang, S. Gu, D. Pan et al., “Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response,” Nature Medicine, vol. 24, no. 10, pp. 1550–1558, 2018.
View at: Publisher Site | Google Scholar
A. Mayakonda, D. C. Lin, Y. Assenov, C. Plass, and H. P. Koeffler, “Maftools: efficient and comprehensive analysis of somatic variants in cancer,” Genome Research, vol. 28, no. 11, pp. 1747–1756, 2018.
View at: Publisher Site | Google Scholar
K. Yoshihara, M. Shahmoradgoli, E. Martínez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, no. 1, p. 2612, 2013.
View at: Publisher Site | Google Scholar
A. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015.
View at: Publisher Site | Google Scholar
T. Li, J. Fu, Z. Zeng et al., “TIMER2.0 for analysis of tumor-infiltrating immune cells,” Nucleic Acids Research, vol. 48, no. W1, pp. W509–W514, 2020.
View at: Publisher Site | Google Scholar
S. Shen, G. Wang, R. Zhang et al., “Development and validation of an immune gene-set based prognostic signature in ovarian cancer,” eBioMedicine, vol. 40, pp. 318–326, 2019.
View at: Publisher Site | Google Scholar
Y. Guo, R. Wang, J. Li et al., “Comprehensive analysis of m6A RNA methylation regulators and the immune microenvironment to aid immunotherapy in pancreatic cancer,” Frontiers in Immunology, vol. 12, article 769425, 2021.
View at: Publisher Site | Google Scholar
L. Zhang, X. Han, and Y. Shi, “Association of MUC16 mutation with response to immune checkpoint inhibitors in solid tumors,” JAMA Network Open, vol. 3, no. 8, article e2013201, 2020.
View at: Publisher Site | Google Scholar
B. Shum, J. Larkin, and S. Turajlic, “Predictive biomarkers for response to immune checkpoint inhibition,” Seminars in Cancer Biology, vol. 79, pp. 4–17, 2022.
View at: Publisher Site | Google Scholar
J. Li, X. Ma, D. Chakravarti, S. Shalapour, and R. A. DePinho, “Genetic and biological hallmarks of colorectal cancer,” Genes & Development, vol. 35, no. 11-12, pp. 787–820, 2021.
View at: Publisher Site | Google Scholar
H. Xiong, H. Li, Y. Xiao et al., “Long noncoding RNA MYOSLID promotes invasion and metastasis by modulating the partial epithelial-mesenchymal transition program in head and neck squamous cell carcinoma,” Journal of Experimental & Clinical Cancer Research : CR, vol. 38, no. 1, p. 278, 2019.
View at: Publisher Site | Google Scholar
Y. Han, N. Wu, M. Jiang et al., “Long non-coding RNA MYOSLID functions as a competing endogenous RNA to regulate MCL-1 expression by sponging miR-29c-3p in gastric cancer,” Cell Proliferation, vol. 52, no. 6, article e12678, 2019.
View at: Publisher Site | Google Scholar
R. Abou Khouzam, R. F. Zaarour, K. Brodaczewska et al., “The effect of hypoxia and hypoxia-associated pathways in the regulation of antitumor response: friends or foes?” Frontiers in Immunology, vol. 13, article 828875, 2022.
View at: Publisher Site | Google Scholar
Z. Jia, J. An, Z. Liu, and F. Zhang, “Non-coding RNAs in colorectal cancer: their functions and mechanisms,” Frontiers in Oncology, vol. 12, article 783079, 2022.
View at: Publisher Site | Google Scholar
P. Peng, K. W. Hsu, J. Chieh-Yu Lai, and K. J. Wu, “The role of hypoxia-induced long noncoding RNAs (lncRNAs) in tumorigenesis and metastasis,” Biomedical Journal, vol. 44, no. 5, pp. 521–533, 2021.
View at: Publisher Site | Google Scholar
D. Barth, F. Prinz, J. Teppan, K. Jonas, C. Klec, and M. Pichler, “Long-noncoding RNA (lncRNA) in the regulation of hypoxia-inducible factor (HIF) in cancer,” RNA, vol. 6, no. 3, p. 27, 2020.
View at: Publisher Site | Google Scholar
P. Gu, L. Zhang, R. Wang et al., “Development and validation of a novel hypoxia-related long noncoding RNA model with regard to prognosis and immune features in breast cancer,” Frontiers in Cell and Developmental Biology, vol. 9, article 796729, 2021.
View at: Publisher Site | Google Scholar
P. Sharma, B. A. Siddiqui, S. Anandhan et al., “The next decade of immune checkpoint therapy,” Cancer Discovery, vol. 11, no. 4, pp. 838–857, 2021.
View at: Publisher Site | Google Scholar
S. Topalian, C. Drake, and D. Pardoll, “Immune checkpoint blockade: a common denominator approach to cancer therapy,” Cancer Cell, vol. 27, no. 4, pp. 450–461, 2015.
View at: Publisher Site | Google Scholar
Y. Xie, Y. Chen, and J. Fang, “Comprehensive review of targeted therapy for colorectal cancer,” Signal Transduction and Targeted Therapy, vol. 5, no. 1, p. 22, 2020.
View at: Publisher Site | Google Scholar
L. You, W. Wu, X. Wang et al., “The role of hypoxia-inducible factor 1 in tumor immune evasion,” Medicinal Research Reviews, vol. 41, no. 3, pp. 1622–1643, 2021.
View at: Publisher Site | Google Scholar
C. Corzo, T. Condamine, L. Lu et al., “HIF-1α regulates function and differentiation of myeloid-derived suppressor cells in the tumor microenvironment,” The Journal of Experimental Medicine, vol. 207, no. 11, pp. 2439–2453, 2010.
View at: Publisher Site | Google Scholar
R. Mhaidly and F. Mechta-Grigoriou, “Role of cancer-associated fibroblast subpopulations in immune infiltration, as a new means of treatment in cancer,” Immunological Reviews, vol. 302, no. 1, pp. 259–272, 2021.
View at: Publisher Site | Google Scholar
L. Ziani, S. Buart, S. Chouaib, and J. Thiery, “Hypoxia increases melanoma-associated fibroblasts immunosuppressive potential and inhibitory effect on T cell-mediated cytotoxicity,” Oncoimmunology, vol. 10, no. 1, article 1950953, 2021.
View at: Publisher Site | Google Scholar
J. Xiao, Y. Liu, J. Yi, and X. Liu, “LINC02257, an enhancer RNA of prognostic value in colon adenocarcinoma, correlates with multi-omics immunotherapy-related analysis in 33 cancers,” Frontiers in Molecular Biosciences, vol. 8, article 646786, 2021.
View at: Publisher Site | Google Scholar
M. Cheng, J. Zhang, P. B. Cao, and G. Q. Zhou, “Prognostic and predictive value of the hypoxia-associated long non-coding RNA signature in hepatocellular carcinoma,” Hereditas, vol. 44, no. 2, pp. 153–167, 2022.
View at: Publisher Site | Google Scholar
D. Yu, X. Wang, and Z. Jin, “Linc00702 inhibits cell growth and metastasis through regulating PTEN in colorectal cancer,” European Review for Medical and Pharmacological Sciences, vol. 24, no. 7, pp. 3624–3632, 2020.
View at: Publisher Site | Google Scholar
W. Yu, D. Li, X. Ding et al., “LINC00702 suppresses proliferation and invasion in non-small cell lung cancer through regulating miR-510/PTEN axis,” Aging, vol. 11, no. 5, pp. 1471–1485, 2019.
View at: Publisher Site | Google Scholar
T. Li, J. Ren, J. Ma et al., “LINC00702/miR-4652-3p/ZEB1 axis promotes the progression of malignant meningioma through activating Wnt/β-catenin pathway,” Biomedicine & Pharmacotherapy = Biomedecine & Pharmacotherapie, vol. 113, article 108718, 2019.
View at: Publisher Site | Google Scholar
L. Wang, T. Y. Ye, H. Wu, S. Y. Chen, J. R. Weng, and X. W. Xi, “LINC00702 accelerates the progression of ovarian cancer through interacting with EZH2 to inhibit the transcription of KLF2,” European Review for Medical and Pharmacological Sciences, vol. 23, Supplement_3, pp. 201–208, 2019.
View at: Publisher Site | Google Scholar
M. Zhou, S. Wang, D. Liu, and J. Zhou, “LINC01915 facilitates the conversion of Normal fibroblasts into cancer-associated fibroblasts induced by colorectal cancer-derived extracellular vesicles through the miR-92a-3p/KLF4/CH25H Axis,” ACS Biomaterials Science & Engineering, vol. 7, no. 11, pp. 5255–5268, 2021.
View at: Publisher Site | Google Scholar
C. Yang and X. Zheng, “Identification of a hypoxia-related lncRNA biomarker signature for head and neck squamous cell carcinoma,” Journal of Oncology, vol. 2022, Article ID 6775496, 2022.
View at: Publisher Site | Google Scholar
L. Shen, N. Li, Q. Zhou, Z. Li, and L. Shen, “Development and validation of an autophagy-related LncRNA prognostic signature in head and neck squamous cell carcinoma,” Frontiers in Oncology, vol. 11, article 743611, 2021.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Xinyang Zhong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

349

Downloads

594

Citations