Identification hypoxia-related lncRNA in EC
To better understand the important role of hypoxia-related lncRNA in oncogenesis and progression, we first downloaded from the hallmark gene sets of hypoxia including 200 genes from Molecular Signatures Database. RNA-seq data from 550 tumor tissue samples and 35 normal samples were downloaded from TCGA. According to the correlation efficiency and probability cut-off value (Filter: |r| > 0.5 and P < 0.001), 677 lncRNAs were screened out as hypoxia-related lncRNA. In order to investigate the prognostic value of these hypoxia-related lncRNA in EC, univariate Cox regression analysis was performed based on the expression levels of these lncRNA in TCGA database. As a result, we found that 17 out of the 677 lncRNAs were significantly associated with overall survival (p < 0.05) (Figure 1A). The expression values of 17 hypoxia-related lncRNAs were extracted from EC patients. The expression profile of 17 prognostic associated hypoxia-related lncRNA was showed in the heatmap and box plot (Figure 1B-C). As shown in Figure 1C, 17 hypoxia-related lncRNAs were significantly abnormally expressed in EC tissues samples. Compared to normal tissue samples, lncRNAs AL021578.1, AP000593.3, LINC01843, SOS1-IT1 and NRAV were up-regulated, while AC079921.1, AC119674.1, AC114947.2, AC011466.1, CFAP58-DT, LINC01614, AC244517.7, TPM1-AS, AL035530.2, AC007786.2, AP003096.1 and AC020910.5 were down-regulated in EC tumor tissues. In addition, based on the co-expression relationship, we also built an interaction network based on hypoxia-related lncRNAs and the corresponding hypoxia associated genes (Figure 1D).
Classification of EC based on hypoxia-related lncRNA
In order to analysis the consensus cluster of hypoxia-related lncRNA in EC, we used the common clusterplus package to identify the different groups of hypoxia-related lncRNA based on their co-expression patterns in EC tissues from the TCGA database. Due to the grouping was suboptimal when they were divided into more than 2 clusters, we divided the hypoxia-related lncRNA into two groups based on their expression indices using k = 2 as the optimal value (Figure 2A-D). Next, we analyzed the relationship between these two clusters and the clinicopathological characteristics of EC patients. We found that consensus clustering could make significant differences in the clinical and molecular characteristics of the two EC clusters (Figure 2E). Cluster 1 patients were strikingly correlated with stage, grade, fustat and disease type by Chi-square test. Besides, compared with patients in cluster 2, EC patients in cluster 1 showed a shorter overall survival time (Figure 2F).
Prognostic value of hypoxia-related lncRNA and construction of a risk signature predicting prognosis
To identified the most powerful prognostic hypoxia-related lncRNAs, the last absolute shrinkage and selection operator (LASSO) Cox regression analysis to the 17 prognosis-related lncRNAs was conducted. We constructed prognostic models using the multivariate Cox proportional hazards regression analysis and the coefficient of each independent prognostic gene were calculated. The risk score was estimated based on the coefficients from the LASSO results. According to the median risk score, EC patients were assigned into low-risk and high-risk groups. The distribution of the hypoxia-related risk signature in the TCGA dataset and survival status of EC patients in different groups were shown in Figure 3A-B. After adjusting for clinicopathological features such as age, grade, stage, radiation therapy, surgical approach, and disease type, we found that age, grade, stage, disease type and risk score were correlated with the OS of EC patients in univariate analysis, while multivariate COX regression analysis showed that age, grade, stage, radiation therapy, surgical approach and risk score were independent risk factors for the prognosis of EC patients (Figure 3C-D).
The heatmap of these most significantly six lncRNAs were shown in Figure 3E, and we observed a strong correlation between the risk score and the clinicopathological characteristics such as stage, age and grade. The result of survival analysis showed that the high-risk group had significantly shorter survival time compared to low-risk group (Figure 3F). The AUC value of the ROC curve in risk score for 5-year survival is 0.691, which is obviously higher than that of ROC in grade (0.662), radiation therapy (0.371), surgical approach (0.461) and disease type (0.569), but lower than stage (0.702) (Figure 3G). In addition, the calibration plot revealed ideal agreements of these genes between the actual observations and 1-, 3- and 5-year OS predicted by the nomogram (Figure 3H). These results indicated that the prognostic index based on hypoxia-related lncRNA has the potential to predict the survival of EC patients.
Relationships between the risk score and clinicopathologic factors
Next, we built a complete prognostic model based on the entire set. The stratification analysis was done according to the age, grade, stage, radiation therapy and surgical approach. EC patients were stratified into age ≤60 and > 60 subgroup, grade G1-G2 and G3-G4 subgroup, stage I-II and III-IV subgroup, surgical approach minimally invasive and open subgroup, radiation therapy No and Yes subgroup. For the EC patients in age > 60 subgroup, grade G3-G4 subgroup, and stage III-IV subgroup, the survival time of patients was significantly shorter than that of patients in another group (Figure 4A), and the average risk score was much higher (Figure 4D-E). However, there is no difference between surgical approach and radiation therapy subgroup. Besides, we found the survival time of patients in high-risk group was significantly shorter in age > 60 subgroup, grade G3-G4 subgroup, stage III-IV subgroup, surgical approach minimally invasive subgroup, and radiation therapy No subgroup (Figure 4B-C).
Biological characteristics and pathway analysis of hypoxia-related risky lncRNA
The EC patients in TCGA were divided into high-risk and low-risk groups based on the median risk score. Then, the lncRNA significantly upregulated (fold change> 1 and p < 0.05) or downregulated (fold change < −1 and p < 0.05) were selected for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEEG) pathway, and Gene Set Enrichment Analysis (GESA). To elucidate the biological functions and pathways that were associated with the risk score, the hypoxia-related lncRNA between the high-risk and low-risk groups were used to perform GO enrichment and KEGG pathway analyses (Figure 5A-D). As expected, hypoxia-related lncRNAs were enriched in several cancer-related biological pathways, such as membrane potential regulation, cell-cell adhesion, and protein digestion signaling pathway, et al. The GSEA analysis results showed that different risk group was involved in multiple significantly enriched pathways, including cell cycle, wnt signaling, spliceosome, et al (Figure 5E). These results showed the two-risk group identified based on the six hypoxia-related risky lncRNAs were closely associated with the malignancy of EC.
The correlation of risk score with immune cell environment characterization in EC
To clarify the potential role of the two clusters (Cluster 1 and 2) in immune cell environment characterization of EC, we investigated the expression levels of the 22 immune cell types infiltration between the two clusters and immune signature by CIBERSORT analysis. As shown in Figure 6A, Cluster 2 patients exhibited higher levels of proportions in CD4 memory resting T cells, activated NK cells and M0 Macrophages. While a higher proportion of CD8 T cells, follicular helper T cells, M1 Macrophages were enriched in Cluster 1 patients. After further uncovering the role of risk score in immune cell environment characterization of EC, we found the risk score was positively correlated with activated dendritic cells, naive B cells, and gamma delta T cell, which meant that hypoxia-related lncRNAs in the signature might influence the tumor immune microenvironment by promoting activated dendritic cells, naive B cells, and gamma delta T cell to infiltrate into the tumor tissue. However, the neutrophils, resting dendritic cells, regulatory T cell (Tregs), CD8 T cells, and activated NK cells were negatively correlated with risk score (Figure 6B-I).
The risk score of hypoxia-related lncRNA was associated with immunotherapy, microsatellite instability and tumor mutation burden in EC
The combination of anti-CTLA-4 and anti-PD-1 treatment could increase the proportion of activated CD8+ cells and natural killer cells in the tumor microenvironment, and decrease the proportion of inhibitory immune cells, resulting in changes in the immune landscape, which achieved therapeutic effect in the mouse model and prolonged the tumor free survival. After analyzing the difference in response to immunotherapy between different risk score group. We found that low-risk group tended to respond effectively to immunotherapy such as anti-PD-1 and anti-CTLA-4 therapy (Figure 7A-D). High microsatellite instability (MSI-H) is associated with the response to immunotherapy treatment. Interestingly, we found that the hypoxia-related lncRNA risk score of EC patients with MSI-H was lower than that of EC patients with low MSI or microsatellite stability (MSS) (Figure 7E-F). Tumor genome mutation leads to the production of new antigens that is beneficial to immunotherapy. Tumor mutation burden (TMB) is an important biomarker, which can be used to predict the response of a variety of tumors to PD-1/PD-L1 targeted immunotherapy. The box plot showed that there was a difference in mutation frequency between different risk score groups (Figure 7G). We observed that the risk score was negatively correlated with TMB in correlation analysis (R=-0.23, P=1.2e-07) (Figure 7H). The EC patients were classified into the high and low TMB group. Kaplan-Meier analysis showed High TMB significantly correlated with better prognosis (Figure 7I). Then we evaluated the synergistic effect of risk score on prognosis stratification. The survival analysis showed that there was significant difference in survival rate according to risk score between different TMB subgroups (Figure 7J).
SOS1-IT1 is clinically relevant in EC and promotes EC cell growth
SOS1-IT1 was the most correlated prognostic lncRNAs in this model. Therefore, we will further evaluate the role of SOS1-IT1 in EC to verify the hypoxia-related lncRNA model. To investigate the clinical significance of SOS1-IT1 in EC, we analyzed its expression and clinical relevance in TCGA EC database. As shown in Figure 8A-B, SOS1-IT1 was overexpressed in tumor tissues. Analysis from clinical investigations suggested that the aberrant level of SOS1-IT1 was closely correlated with clinicopathological parameters of EC, including the age, disease type, fustat and grade (Figure 8C-F). Kaplan-Meier analysis revealed that high expression of SOS1-IT1 was significantly associated with a poor prognosis in EC patients (Figure 8G). To directly investigate the biological functions of SOS1-IT1 on EC cells, SOS1-IT1 was knocked down or overexpressed respectively in EC cell line Ishikawa cells. The real-time PCR analysis results indicated that the expression level of SOS1-IT1 was markedly decreased or increased in knocking down or overexpressing cells respectively (Figure 8H-I). The downregulation of SOS1-IT1 significantly inhibited EC cell growth. In contrast, overexpressing SOS1-IT1 yielded the opposite results (Figure 8J-K). These results indicate that SOS1-IT1 may lead to increased EC aggressiveness.
SOS1-IT1 is upregulated under hypoxia and directly transactivated by HIF-1α
To further confirm whether SOS1-IT1 was a functional effector of hypoxia in EC progression, we treated Ishikawa cells with hypoxia or its chemical inducer CoCl2 for 24 h, and found the expression level of SOS1-IT1 was significantly increased (Figure 9A). Besides, the expression level of SOS1-IT1 was significantly decreased after knocking down HIF-1α, which is the main signaling pathway response component to hypoxia, both in normoxia and hypoxia condition (Figure 9 B-C). After detecting of the upstream region (~1KB upstream) of SOS1-IT1 gene by promoter sequence analysis, we found a putative HIF-1α binding site in the promoter region of SOS1-IT1 gene (Figure 9D). We generated the mutant binding site of the reporter constructs, following with luciferase assays after transfection of different reporter constructs in EC cells. Under normoxia conditions, there was no much difference of luciferase activity in SOS1-IT1 wild type (WT) or mutant, as well as empty vector group. However, under hypoxic conditions, there was a nearly six-fold induction of luciferase activity in SOS1-IT1 WT construct. Mutation in SOS1-IT1 promoter decreased the luciferase activity to nearly basal levels caused by hypoxia (Figure 9E). Besides, the Chromatin immunoprecipitation (ChIP) assays confirmed that HIF-1α directly bind to the promoter region promoter of SOS1-IT1 (Figure 9F). Next, we knocked down HIF-1α, and found HIF-1α suppression remarkably repressed the luciferase density in cells with WT promoter, but not in mutant group (Figure 9G). Moreover, SOS1-IT1 silencing partially reduced EC cell proliferation ability by hypoxia condition or HIF-1α overexpression (Figure 9H-J). Totally, these results intensively indicated that SOS1-IT1 was upregulated under hypoxia and a direct transcriptional target of HIF-1α.