Peng Gu1†
Lei Zhang2†
Ruitao Wang1†Wentao Ding1
Wei Wang1
Yuan Liu1Wenhao Wang3Zuyin Li4
Bin Yan1*Xing Sun1*- 1Department of General Surgery, Shanghai General Hospital, School of Medicine, Shanghai Jiao Tong University, Shanghai, China
- 2Department of Vascular Surgery, Intervention Center, Shanghai General Hospital, School of Medicine, Shanghai Jiao Tong University, Shanghai, China
- 3Department of Urology, Shanghai General Hospital, School of Medicine, Shanghai Jiao Tong University, Shanghai, China
- 4Department of Hepatobiliary Surgery, Peking University Organ Transplantation Institute, Peking University People’s Hospital, Beijing, China
Background: Female breast cancer is currently the most frequently diagnosed cancer in the world. This study aimed to develop and validate a novel hypoxia-related long noncoding RNA (HRL) prognostic model for predicting the overall survival (OS) of patients with breast cancer.
Methods: The gene expression profiles were downloaded from The Cancer Genome Atlas (TCGA) database. A total of 200 hypoxia-related mRNAs were obtained from the Molecular Signatures Database. The co-expression analysis between differentially expressed hypoxia-related mRNAs and lncRNAs based on Spearman’s rank correlation was performed to screen out 166 HRLs. Based on univariate Cox regression and least absolute shrinkage and selection operator Cox regression analysis in the training set, we filtered out 12 optimal prognostic hypoxia-related lncRNAs (PHRLs) to develop a prognostic model. Kaplan–Meier survival analysis, receiver operating characteristic curves, area under the curve, and univariate and multivariate Cox regression analyses were used to test the predictive ability of the risk model in the training, testing, and total sets.
Results: A 12-HRL prognostic model was developed to predict the survival outcome of patients with breast cancer. Patients in the high-risk group had significantly shorter median OS, DFS (disease-free survival), and predicted lower chemosensitivity (paclitaxel, docetaxel) compared with those in the low-risk group. Also, the risk score based on the expression of the 12 HRLs acted as an independent prognostic factor. The immune cell infiltration analysis revealed that the immune scores of patients in the high-risk group were lower than those of the patients in the low-risk group. RT-qPCR assays were conducted to verify the expression of the 12 PHRLs in breast cancer tissues and cell lines.
Conclusion: Our study uncovered dozens of potential prognostic biomarkers and therapeutic targets related to the hypoxia signaling pathway in breast cancer.
Background
Female breast cancer surpassed lung cancer as the most frequently diagnosed cancer (representing 11.7% of total cases and 24.5% of female cancers) worldwide in 2020 (Sung et al., 2021). Despite recent advances in high-quality prevention, early detection, and treatment services, breast cancer is still the leading cause of cancer-related death for female patients (15.5%) (Loibl et al., 2021a). Thus, it is of great importance to identify novel prognosticators and develop a more precise prognostic model to help optimize the individual treatment for patients.
Tumor cells grow in and continuously interact with an incredibly complex and dynamic network called tumor microenvironment (TME). TME consists of noncellular components, such as extracellular matrix, growth factors, cytokines, enzymes, and hormones, as well as diverse types of cells, including endothelial cells, pericytes, immune inflammatory cells, and fibroblasts (Hanahan and Weinberg, 2011). Hypoxia, as one of the most important abnormal microenvironments that tumor cells are continuously exposed to, has a significant impact on tumor biology, leading to a higher phenotypic heterogeneity (Marusyk et al., 2012). The hypoxic tissue areas of many breast cancers are heterogeneously distributed within the tumor mass (Vaupel et al., 2004). Oxygen tension (pO2) measured in normal breast tissue exhibited a mean (and median) of 65 mm Hg (Vaupel et al., 2007), whereas for breast cancers in T1b–T4 stages, the median pO2 was 28 mm Hg (Vaupel et al., 1991). The hypoxic response is mainly ascribed to HIF. The HIF is a heterodimer comprising an inducible α subunit (HIF-α) and a constitutively expressed β subunit (HIF-β), which functions as a transcriptional factor that regulates the expression of genes involved in diverse biological characteristics of tumors (Wang et al., 1995). The levels of HIF-1α have been implicated as an independent prognostic factor for patients with breast cancer (Bos et al., 2003). The hallmarks of cancer, such as angiogenesis (Niu et al., 2020), invasion and metastasis (Bristow and Hill, 2008), evading immune destruction (Hsing et al., 2012) and reprogramming of energy metabolism (Samanta and Semenza, 2018), have been validated to be tumor supportive and associated with the hypoxic tumor environment of breast cancer.
LncRNAs are defined as RNA transcripts longer than 200 nucleotides in length, with no potential to encode proteins (Kopp and Mendell, 2018). Accumulating evidence suggests that lncRNAs play crucial roles in the occurrence and development of a large variety of cancers including breast cancer (Gao et al., 2020; Liu et al., 2021). During tumor hypoxia, hypoxia-inducible factors (HIFs) are upregulated and either positively or negatively regulate lncRNAs through hypoxia response elements within their promoters. Or vice versa, some lncRNAs can regulate the HIF signaling pathway directly or indirectly (Choudhry et al., 2016). Moreover, increasing evidence reveals the potential of lncRNAs as diagnostic or prognostic biomarkers (Bin et al., 2018; Qian et al., 2020; Chen Y. et al., 2021).
Though there was a similar study proposing a prognostic signature based on HRLs (Zhao et al., 2021), which served to stratify patients with early-stage breast cancer, we aimed to develop and validate a more comprehensive and reliable prognostic model for predicting the survival of all breast cancer patients. By means of both bioinformatic analysis and experimental validation, we were also meant to mine some potential tumor-supportive or tumor-suppressive lncRNAs to indicate the further research direction on their functions in breast cancer.
Materials and Methods
Data Acquisition
We downloaded RNA sequencing data (1,109 breast cancer tissues and 113 matched normal tissues) and corresponding clinical and pathological information of patients with breast cancer from The Cancer Genome Atlas (TCGA) (https://gdc.cancer.gov/) database in March 2021. The clinical characteristics of the patients are listed in Table 1. A total of 200 hypoxia-related mRNAs were obtained from the Molecular Signatures Database V7.3 (Subramanian et al., 2005) (https://www.gsea-msigdb.org/gsea/msigdb/, M5891).
TABLE 1. Baseline patients characteristics (n = 1,097).
Screening of Hypoxia-Related lncRNAs
Differentially expressed hypoxia-related mRNAs (DEHmRNAs) between the breast cancer and normal tissues were identified via the differential expression analysis using the R package “DEseq” with the threshold value set as |log2 fold change (FC)| >1 and adjusted p value of <0.05. Heatmap and volcano plots were drawn using the R packages “pheatmap” and “EnhancedVolcano,” respectively. The PPI (protein–protein interaction) network was constructed using the STRING version 11.0 Program (Szklarczyk et al., 2017). Further, we performed Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis via the R package “clusterProfiler”. Finally, we identified the HRLs by Spearman correlation analysis of DEHmRNAs according to the criteria of |Correlation Coefficient| > 0.4 and p < 0.001. The mRNAs–lncRNAs network were constructed via Cytoscape 3.8.2.
Construction of a Hypoxia-Related lncRNA Prognostic Model
Univariate Cox regression analysis was used to identify prognostic hypoxia-related lncRNAs via the R package “survival” with p < 0.05 as the criteria. Heatmap and box plot were drawn using the R packages “pheatmap” and “ggplot2,” respectively. Spearman’s rank correlation was performed to assess the correlation among these genes using the R package “corrplot.” Then, we randomly divided patients into a training set and a testing set in the ratio of 1:1 using the R package “caret.” The R package “glmnet” was used to perform least absolute shrinkage and selection operator (LASSO) Cox regression analysis in the training set to generate a coefficient for each hypoxia survival-related lncRNA. Then, the risk score for each sample was calculated based on the formula: Risk Score (RS) =
Validation of the Model
Kaplan–Meier survival analysis was conducted to investigate the differences in OS and DFS between high- and low-risk groups using the R packages “survival” and “survminer,” no matter whether we classified the patients into different groups based on their clinicopathological characteristics. Receiver operating characteristic (ROC) curves and area under the curve (AUC) were used to evaluate the sensitivity and specificity of this prognostic model via the R package “survivalROC.” Additionally, we used “survival” R package to perform univariate and multivariate Cox regression analyses so as to estimate the prognostic value of the model constructed. Finally, we formulated a nomogram using the R package “rms,” which could assign points for independent prognostic factors to predict 1-, 3-, 5-, and 10-year OS of individual patients with breast cancer (Sui et al., 2020). Meanwhile, we calculated the concordance index (C-index) using “survcomp” and constructed calibration curves to evaluate the predictive power of the nomogram.
Drug Sensitivity Analysis
We used the R package “pRRophetic” to predict the differences in chemosensitivity from tumor gene expression levels for some chemotherapy drugs, which was decided by the half-maximal inhibitory concentration (IC50), between high- and low-risk groups based on our prognostic model (Geeleher et al., 2014).
Gene Set Enrichment Analysis
We performed gene set enrichment analysis (GSEA) (version 4.1.0, https://www.gsea-msigdb.org/gsea/index.jsp) to identify the related differential biological function and pathways between high- and low-risk groups as separated by the prognostic model using the following gene sets: Hallmark, GO, KEGG and BioCarta. Each analysis included 1,000 random permutations. The statistical significance level was set to be p < 0.05 and the false discovery rate as < 0.25. The plots were drawn using the R package “ggplot2.”
Immune Cell Infiltration Analysis
The CIBERSORT(Floberg et al., 2021), CIBERSORT-ABS(Yuan et al., 2021), TIMER (tumor immune estimation resource) (Jiang et al., 2021), xCELL (Sui et al., 2020), quanTIseq (Subramanian et al., 2005), EPIC (extended polydimensional immunome characterization) (Yeo et al., 2020), MCPcounter (Dienstmann et al., 2019) algorithms were used to further analyze the differences in immune cell infiltration between high- and low-risk groups. The immune score for each single sample was calculated using the R package “estimate.” Additionally, we calculated the correlation between infiltrating immune cells and risk scores using Spearman’s rank correlation with the criteria of p < 0.05. Subsequently, 11 immune checkpoints were selected from the previous literature (Marin-Acevedo et al., 2018).
RNA Extraction and Real-Time Quantitative PCR
Total RNA was extracted from 11 breast cancer tissues, paired adjacent normal tissues, and breast cell lines including MCF10A and MDA-MB-231 (under normoxia or hypoxia conditions) cell lines. RT-qPCR experiments were performed as previously described. The primers used in this study were shown in Supplementary Table S1. The mRNA quantification of the 12 PHRLs was based on the 2–∆∆Ct method and the expression levels were plotted by using ACTB as the reference gene.
Statistical Analysis
Data were processed using Perl language (version 5.26.1) and analyzed using R software (version 1.3.1073). The Wilcoxon test was applied to compare gene expression, risk scores, drug sensitivity (IC50), and immune cell infiltration scores between two independent groups. Univariate Cox regression and LASSO Cox regression analyses were applied to identify the most useful lncRNAs to build a prognostic model. Survival analysis was conducted using the Kaplan-Meier method and log-rank tests. The chi-square test was used to compare the relationship between risk and immune score level and other categorical clinicopathological factors. Correlation analysis was assessed by Spearman’s rank correlation. Statistical significance was set at p < 0.05 or 0.001 for diverse analysis.
Results
Identification of Differentially Expressed Hypoxia-Related mRNAs in Breast Cancer
We obtained 200 hypoxia-related mRNAs from the Molecular Signatures Database V7.3 (HALLMARK_HYPOXIA, M5891), which were previously validated to be hypoxia-regulated in multiple kinds of cancer including breast cancer (Wang et al., 2009; Xia and Kung, 2009; Leithner et al., 2014). By comparing their expression levels between 1,109 breast cancer tissues and 113 normal tissues using the R package “DEseq,” we identified 46 DEHmRNAs (Supplementary Table S2). The heatmap and volcano plot are shown in Figures 1A,B. A PPI network showed that almost all proteins encoded by the aforementioned genes interacted with each other except SRPX and TPBG (Figure 1C). Further, we performed KEGG (Figure 1D) and GO (Figures 1E–G) analyses to investigate the biological function of 46 DEHmRNAs. The top 3 most enriched KEGG pathways were “HIF-1 signaling pathway”, “Glycolysis/Gluconeogenesis”, and “Carbon metabolism”.
FIGURE 1. Differential expression and functional annotation of DEHmRNAs in breast cancer. (A) Heatmap and (B) volcano plot of DEHmRNAs between breast cancer tissues and normal tissues (HTSeq-Counts data from TCGA, “DESeq”). (C) A PPI network of DEHmRNAs. (D) KEGG pathway and (E–G) GO function enrichment analysis of DEHmRNAs.
Screening of Hypoxia-Related lncRNAs and Construction of a Co-Expression Network
Co-expression analysis between DEHmRNAs and lncRNAs based on Spearman’s rank correlation was performed to screen out 166 HRLs according to the criteria of |Correlation Coefficient| >0.4 and p < 0.001 (Supplementary Table S3), of which 46 were significantly upregulated and 41 were downregulated (Supplementary Figure S1 and Supplementary Table S4). Based on the mRNA–lncRNA co-expression pattern, we constructed a network to show the relationships among them (Supplementary Figure S2).
Development of Hypoxia-Related lncRNA Prognostic Model
By matching the sample IDs in the expression matrix and clinical information profile, we collected 1,090 samples for subsequent analysis after excluding 19 samples. We performed univariate Cox regression analysis to identify 20 PHRLs significantly associated with survival (p < 0.05) (Figure 2A). The differential expression of these 20 PHRLs is shown in Figures 2B,C. Most of them were notably correlated with each other (Figure 2D). Then, we randomly divided these 1,090 patients with breast cancer into two groups in the ratio of 1:1: a training set (n = 546) and a testing set (n = 544). LASSO-penalized Cox regression was performed in the training set, and finally a prognostic model was developed, which consisted of seven risk lncRNAs and five protective lncRNAs on the basis of the coefficient value (Figures 2E–G). We calculated the risk scores for all samples using the following formula: Risk score = (0.309031851 × Exp of TDRKH-AS1) + (0.213481355 × Exp of AC011978.2) + (0.170772984 × Exp of AC110995.1) + (0.051734285 × Exp of OTUD6B-AS1) + (0.046646469 ×Exp of YTHDF3-AS1) + (0.019875191 × Exp of AL512380.1) + (0.019087737 × Exp of MIR4435-2HG) + (–0.037067764 × Exp of HSD11B1-AS1) + (–0.063799579 × Exp of LINC02084) + (–0.162152639 × Exp of TRG-AS1) + (–0.399031951 × Exp of AL451085.3) + (–0.886350422 × Exp of AL109955.1). All the patients were further separated into high- and low-risk groups according to the median risk score in the training set.
FIGURE 2. Development of a hypoxia-related lncRNA prognostic model. (A) Forest plot showed the hazard ratio (HR, 95% CI) and the p value of 20 PHRLs by univariate Cox regression analysis. Heatmap (B) and box plot (C) showed the expression level of the 20 PHRLs in breast cancer tissues compared with normal tissues (Wilcoxon test). (D) Spearman rank’s correlation analysis based on the expression of the 20 PHRLs. (E and F) LASSO-penalized Cox regression was performed to filter out 12 optimal PHRLs. (G) 12 PHRLs and their coefficients. *p < 0.05, **p < 0.01, and ***p < 0.001.
Validation of the Prognostic Model
We performed Kaplan–Meier survival analysis to investigate whether there any differences in survival outcomes existed between the high and low-risk groups. The results shown in Figure 3A indicated that in the training set, patients in the high-risk group had shorter median OS compared with those in the low-risk group (log-rank test, p = 4.385e-08). This significant difference can also be seen in both the testing set and the total set (p = 0.0026 in the former, p = 8.537e-10 in the latter) (Figures 3B,C). Then we built time-dependent ROC curves for the three sets. The AUC for 1-, 3-, 5-, and 10-year OS was 0.734, 0.727, 0.741, and 0.786 in the training set and 0.681, 0.637, 0.612, and 0.749 in the testing set, respectively (Figure 3D). In the total set, it was 0.707, 0.691, 0.684, and 0.769 (Figures 3E,F). Additionally, we compared the DFS of 984 patients with breast cancer. As expected, the patients in the high-risk group had significantly shorter median DFS than those in the low-risk group when the best cut point of the risk score was set as 1.55 (Figures 3G–I).
FIGURE 3. Validation of the predictive ability of the hypoxia-related lncRNA prognostic model. Kaplan–Meier survival analysis revealed the differences in survival outcomes between high- and low-risk groups in the training (A), testing (B), and total sets (C). Time-dependent ROC analysis for evaluating the sensitivity and specificity of this prognostic model in the training (D), testing (E), and total sets (F). (G) Survival curves of DFS between high- and low-risk groups. (H and I) The best cutpoint of the risk score for (G).
We next visualized risk score distribution, survival status, and PHRLs expression in all three sets. The plots showed that, whether in the training (Figures 4A,D), testing (Figures 4B,E), or total Figures 4C,F) set, patients in the high-risk group had poorer survival and higher risk scores compared with those in the low-risk group. The heatmaps manifested the differential expression pattern of PHRLs between the high- and low-risk groups (Figures 4G–I).
FIGURE 4. Risk score distribution, survival status, and 12 PHRLs expression. (A,D,G) Training set. (B,E,H) Testing set. (C,F,I) Total set.
We excluded 305 patients with incomplete clinical and pathological information to further test the independent predictive ability of the risk model. In order to ensure the accuracy of the following analyses, we compared the difference in clinicopathological factors of breast cancer patients between training (n = 383) and testing (n = 402) sets (Supplementary Table S5). As expected, there was no significant difference between the two sets randomly sampled. Then, we performed univariate and multivariate Cox regression analyses. The results showed that the risk score was significantly associated with OS (p ≤ 0.001) after adjustment for age, stage, TNM stage, ER/PR status, and HER2 status in the training (Figures 5A,D), testing (Figures 5B,E), and total sets (Figures 5C,F), suggesting that the risk score could act as an independent prognostic factor.
FIGURE 5. Identification of the risk score as an independent prognostic factor by Cox regression analysis. (A,D) Training set. (B,E) Testing set. (C,F) Total set.
Moreover, we divided patients in the total set into multiple groups according to their clinicopathological characteristics and then used Kaplan–Meier survival analysis to verify the prognostic value of the risk model for them. The results shown in Figures 6A–M indicated that the prognosis of patients in the high-risk group was significantly worse than that of the patients in the low-risk group. It is worth noting that this model was applicable to patients with breast cancer of three different molecular subtypes [HR-positive/luminal, HER2-positive, or triple-negative breast cancer (TNBC)]. By comparing the risk scores of patients in different groups, we found that patients over 65 years and those with higher clinical stage had higher risk scores while patients with TNBC had lower risk scores (Supplementary Figure S3).
FIGURE 6. Survival curves of OS between high- and low-risk groups for patients classified in different ways. (A–C) Patients with HR-positive/luminal breast cancer, HER2-positive breast cancer, or TNBC. (D,E) Patients aged ≤65 or >65 years. (F,G) Patients with clinical stage I–II or III–IV. (H,I) Patients with T1–2 or T3–4. (J,K) Patients with M0 or M1. (L,M) Patients with N0 or N1–3.
Establishment and Validation of a Nomogram and Drug Sensitivity Analysis
The independent prognosticators, including age, stage, ER, PR, HER2 status, and risk score, were used to establish a nomogram for the total set, which could assign a point for each subgroup of these prognosticators. Then, we could calculate the total points to predict the 1-, 3-, 5-, and 10-year OS. As shown in Figure 7A, a 54-year-old patient with clinical stage II, risk score of 10.02, and TNBC had a total point of 236 and predicted 1-, 3-, 5-, and 10-year OS of 93.8, 70.3, 49.7, and 7.2%, respectively. The C-index was 0.816 (95% CI 0.760–0.873) in the training set (n = 383), 0.800 (95% CI 0.725–0.874) in the testing set (n = 402), and 0.797 (95% CI 0.746–0.848) in the total set (n = 785). Moreover, we constructed calibration curves for three sets to predict 5-year OS, which validated the great repeatability and reliability of the established nomogram (Figures 7B–D). Intrinsic and acquired resistance to chemotherapy remains a major challenge in effective breast cancer treatment. Hence, we calculated the differences in drug sensitivity, which was decided by IC50 between high- and low-risk groups. As predicted, patients in the high-risk group had lower chemosensitivity to paclitaxel (p = 2.3e-13), docetaxel (p = 3.8e-07), and doxorubicin (not significant) compared with those in the low-risk group (Figures 7E–G).
FIGURE 7. Establishment and validation of a nomogram for predicting the OS of patients with breast cancer along with the drug sensitivity analysis. (A) A nomogram considering multiple independent prognosticators, including age, stage, ER, PR, HER2 status, and risk score for predicting the 1-, 3-, 5-, and 10-year OS of patients in the total set. Calibration curves of the nomogram predicting 5-year OS in the training set (B), testing set (C), and total set (D). (We did not show the nomograms of the training set and the testing set.) Drug sensitivity of patients in high- and low-risk groups to paclitaxel (E), docetaxel (F), and doxorubicin (G) in the total set. ***p < 0.001.
Functional Annotation in High- and Low-Risk Groups by GSEA
We performed GSEA using the gene sets, including Hallmark (Figure 8A), GO (Figure 8B), KEGG (Figure 8C), and BioCarta (Figure 8D), to identify the related differential biological function and pathways between high- and low-risk groups. Remarkably, multiple immune-related signaling pathways were enriched in the low-risk group, such as inflammatory response, interferon gamma response, TNFA signaling via NFκB, T cell activation, regulation of natural killer cell mediated immunity, T cell differentiation, T cell receptor signaling pathway, B cell receptor signaling pathway, IL12 pathway, Th1Th2 pathway, NKT pathway, and so on. Based on the GSEA results, we would like to further explore the difference in immune infiltration between the groups.
FIGURE 8. Functional annotation in high- and low-risk groups by GSEA. The GSEA results based on the gene sets, including Hallmark (A), GO_BP (B), KEGG (C), and BioCarta (D). BP, Biological process.
Immune Cell Infiltration Analysis
With the CIBERSORT algorithm, we compared the differential infiltration levels of 22 kinds of immune cells for each sample between the high- and low-risk groups. As shown in Figure 9A, plasma cells (p < 0.001), resting NK cells (p = 0.049), M0 macrophages (p < 0.001), M2 macrophages (p < 0.001), and neutrophils (p < 0.001) were significantly enriched in the high-risk group while the proportions of naïve B cells (p < 0.001), CD8+T cells (p < 0.001), regulatory T cells (Tregs) (p < 0.001), and activated NK cells (p = 0.004) were significantly higher in the low-risk group. To validate the observed differences, we divided the patients into high- and low-immune score groups according to the median immune score. We then compared the risk scores between them (Wilcoxon test). Obviously, patients with high immune scores had lower risk scores (p < 0.001) (Figure 9B). Consistently, the immune scores of patients in the high-risk group were lower than those of patients in the low-risk group (p < 0.001) (Figures 9C,D). Additionally, we calculated the correlation between infiltrating immune cells and risk scores and exhibited it in a bubble chart (p < 0.05) (Figure 9E). Immune checkpoints are immunomodulators of both stimulatory and inhibitory pathways (Pardoll, 2012). We obtained 11 immune checkpoints from the previous literature (Marin-Acevedo et al., 2018) and then compared their expression level between the high- and low-risk groups. Most of them (GITR, OX40, CD137, CD40LG, CD28, CD278, CTLA4, VSIR, and CD223) were dramatically downregulated in the high-risk group (Figures 9F,G).
FIGURE 9. Differential immune features between risk groups. (A) CIBERSORT algorithm results demonstrated the differential infiltration levels of 22 kinds of immune cells for each sample in TCGA between the high- and low-risk groups (Wilcoxon test). (B) Comparison of immune scores of patients between the high- and low-risk groups (Wilcoxon test). (C) Comparison of risk scores of patients between the high- and low-immune score groups (Wilcoxon test). (D) A comprehensive heatmap integrating the expression of 12 PHRLs and the distribution of multiple clinicopathological factors (HER2, PR, ER, TNM, stage, age, and immune score) between the high- and low-risk groups (chi-square test). (E) A bubble chart exhibited the correlation between infiltrating immune cells and risk scores using CIBERSORT, CIBERSORT-ABS, TIMER, xCELL, quanTIseq, EPIC, and MCPcounter algorithms (Spearman’s rank correlation). (F,G) Heatmap and box plot showed the expression of 11 immune checkpoints in the high- and low-risk groups (Wilcoxon test). *p < 0.05, **p < 0.01, and ***p < 0.001.
Validation of the Expression of the 12 PHRLs in Breast Cancer Tissues and Breast Cell Lines
To verify the differential expression of the 12 PHRLs, 11 breast cancer tissues and the corresponding adjacent normal tissues were collected for RT-qPCR assay. Two lncRNAs were significantly overexpressed in tumor tissues while six were significantly underexpressed (Figure 10). We further explored the expression of the 12 PHRLs in MCF10A (normoxia, 21%O2), MDA-MB-231 (normoxia, 21%O2), and MDA-MB-231 (hypoxia, 1%O2). The results were shown in Figure 11. We next carried out the Kaplan–Meier survival analysis for the 12 PHRLs using GEPIA2 database (http://gepia2.cancer-pku.cn/) except AL512380.1 for its median expression being zero (Figure 12). Based on the transcriptomic data from TCGA and the Genotype-Tissue Expression (GTEx) database, we analyzed the differential expression of the 12 PHRLs between breast cancer tissues or TNBC tissues and normal breast tissues (Supplementary Figures S4 and S5). Taken together, we could also draw a conclusion that TDRKH-AS1, MIR4435-2HG, HSD11B1-AS1, TRG-AS1, and AL451085.3 were valuable prognostic indicators for patients with breast cancer.
FIGURE 10. The expression of 12 PHRLs in 11 breast cancer tissues and paired normal tissues. (A)TDRKH-AS1, (B) AC011978.2, (C) AC110995.1, (D) OTUD6B-AS1, (E) YTHDF3-AS1, (F) AL512380.1, (G) MIR4435-2HG, (H) HSD11B1-AS1, (I) LINC02084, (J) TRG-AS1, (K) AL451085.3 and (L) AL109955.1.
FIGURE 11. The expression of 12 PHRLs in breast cell lines. (A)TDRKH-AS1, (B) AC011978.2, (C) AC110995.1, (D) OTUD6B-AS1, (E) YTHDF3-AS1, (F) AL512380.1, (G) MIR4435-2HG, (H) HSD11B1-AS1, (I) LINC02084, (J) TRG-AS1, (K) AL451085.3 and (L) AL109955.1. hyp: hypoxia (1% O2) treatment for 24 h.
FIGURE 12. Kaplan-Meier survival analysis of breast cancer patients based on the expression of PHRLs using GEPIA2 database. (A) TDRKH-AS1, (B) AC011978.2, (C) AC110995.1, (D) OTUD6B-AS1, (E) YTHDF3-AS1, (F) MIR4435-2HG, (G) HSD11B1-AS1, (H) LINC02084, (I) TRG-AS1, (J) AL451085.3 and (K) AL109955.1.
Discussion
Breast cancer remains the most common cancer worldwide, with continuously increasing incidence and various factors affecting its development (Bray et al., 2015; Loibl et al., 2021b). Hypoxic areas are distributed heterogeneously throughout the tumor mass and surrounding environment. The crosstalk between tumor cells and the local microenvironment contributes to carcinogenesis, metastasis, and chemoresistance and even determines clinical outcomes. Over the past decade, an increasing number of lncRNAs have been verified to play crucial roles under hypoxia TME in breast cancer. For example, lncRNA BCRT1, which is transcriptionally regulated by HIF-1α under hypoxic conditions, promoted breast cancer cell proliferation and progression (Byrne et al., 2020). Interestingly, another study showed that hypoxia-induced lincRNA-P21 promoted the Warburg effect to increase ATP generation by regulating HIF-1α transcriptional activity in breast cancer (Yang et al., 2014). Most of these studies focused on the role of lncRNAs in the occurrence and development of breast cancer.
Although several articles concerning prognostic models in breast cancer have been published (Li et al., 2020; Pei et al., 2020; Zhao et al., 2021), in this study, we selected a well-recognized gene set “HALLMARK_HYPOXIA” which consists of 200 canonical hypoxia-related mRNAs to perform a differential expression analysis between breast cancer tissues and normal tissues and obtained 46 DEHmRNAs. The co-expression analysis was performed to screen out 166 HRLs. Using these HRLs, we finally constructed a novel prognostic model for patients with breast cancer, which consisted of 12 PHRLs: TDRKH-AS1, AC011978.2, AC110995.1, OTUD6B-AS1, YTHDF3-AS1, AL512380.1, MIR4435-2HG, HSD11B1-AS1, LINC02084, TRG-AS1, AL451085.3, and AL109955.1 (RBM38-AS1). The PHRLs above were totally different from those four of Zhao et al.’s: AL031316.1, AC004585.1, LINC01235, and ACTA2-AS1. The principle reason for this difference is the different hypoxia-related gene sets we chose, which led to only 13 overlapped DEHmRNAs and 34 overlapped HRLs between our study and theirs in the co-expression network (Supplementary Table S6). Nevertheless, the KEGG analysis for the 46 DEHmRNAs we filtered out verified the accurate and strong correlation with hypoxia. Moreover, our prognostic model could divide the patients with breast cancer into high- and low-risk groups more effectively in both training and testing sets and was found to have greater prognostic value by multiple verification methods. The AUCs for 1-, 3-, 5-, and 10-year OS in three sets were calculated to evaluate the predictive accuracy. Beside OS, the difference in DFS between high- and low-risk patients was also significant. It’s also worth noting that our prognostic model was applicable to all breast cancer patients with different clinical stages or molecular subtypes (HR-positive/luminal, HER2-positve, or TNBC) while Zhao et al.’s model applied only to early-stage breast cancer patients. When we compared the risk scores of patients in different groups, we found that patients over 65 years or those with later clinical stage had higher risk scores while patients with TNBC had lower risk scores. There may be two reasons accounting for this result: For one thing, among the 208 patients over 65 years in TCGA breast cancer database, only 21 of them were with TNBC; for another, there were merely 23 TNBC patients among a total of 178 patients with late stage (stage III or stage IV). This seems to indicate that the roles of age and clinical stage as clinicopathological factors outweigh the role of molecular subtypes to predict the survival of patients in the current study, which has also been proved by previous studies (Ferguson et al., 2013; O'Brien et al., 2018). The nomogram established in the present study, which took the clinicopathological factors (HER2, PR, ER, stage, and age) and risk scores into consideration, could precisely predict the OS of patients with breast cancer. The drug sensitivity analysis showed that patients in the high-risk group may be more likely to be resistant to chemotherapy drugs including paclitaxel, docetaxel, and doxorubicin. Though we failed to retrieve ideal GEO datasets including all the 12 PHRLs for independent validation, the validation process above was still sufficient to underscore the utility of our model in predicting the prognosis of breast cancer patients.
Among the 12 PHRLs, AC011978.2, AC110995.1, YTHDF3-AS1, AL512380.1, HSD11B1-AS1, AL451085.3, and AL109955.1 have not been reported to date. TDRKH-AS1 was reported to promote colorectal cancer (CRC) progression through the Wnt/β-catenin signaling pathway (Jiao et al., 2020). In addition, its differential expression with copy number alteration was positively associated with longer OS in lung adenocarcinoma (Wang L. et al., 2019). OTUD6B-AS1 was demonstrated to indicate poor prognosis in ovarian cancer, clear cell renal cell carcinoma (ccRCC), and breast cancer (as an immune-related lncRNA) while its overexpression inhibited ccRCC proliferation (Wang G. et al., 2019; Li and Zhan, 2019; Ma et al., 2020). Similarly, OTUD6B-AS1 played a tumor-suppressive role in thyroid carcinoma (Wang Z. et al., 2020), bladder carcinoma (Wang Y. et al., 2020), and CRC (Cai et al., 2021; Wang et al., 2021). In contrast, in hepatocellular carcinoma (HCC), OTUD6B-AS1 could enhance cell proliferation and invasion ability via the GSKIP/Wnt/β-catenin signaling pathway (Kong et al., 2020). MIR4435-2HG has been verified to be tumor supportive in multiple forms of cancer including breast cancer (Chen D. et al., 2021). Linc02084, as a low-risk immune-related lncRNA in ccRCC(Sun Z. et al., 2020), could also be used to construct a classifier for predicting early recurrence in HCC after curative resection (Lv et al., 2018). TRG-AS1 acted as a molecular sponge to stimulate tongue squamous cell carcinoma (He et al., 2020), HCC (Sun X. et al., 2020), and glioblastoma (Xie et al., 2019) progression. By means of both bioinformatic analysis and experimental validation, our research also indicated that TDRKH-AS1, MIR4435-2HG, HSD11B1-AS1, TRG-AS1, and AL451085.3 were likely to play tumor-supportive or tumor-suppressive roles in the development and progression of breast cancer and may be valuable independent prognostic indicators for patients with breast cancer.
Immune escape has emerged as a key mechanism for breast cancer progression and a crucial step in the preinvasive-to-invasive transition (Gil Del Alcazar et al., 2020). The hypoxic areas in solid tumors are highly infiltrated with immunosuppressive cells, such as tumor-associated macrophages (TAMs), myeloid-derived suppressor cells, and Tregs (Noman et al., 2014). Liang et al. showed that hypoxia-induced exosomal lncRNA BCRT1 contributed to M2 phenotype polarization of TAMs and enhanced its tumor-promoting function (Liang et al., 2020). Ben-Shoshan et al. revealed that hypoxia induced the differentiation of nonspecific CD4+ T cells into functionally active Foxp3 + CD4+CD25+ Treg cells to initiate an anti-inflammatory program via HIF-1α(Ben-Shoshan et al., 2008). Neutrophils with HIF2A gain of function displayed a reduction of apoptosis both ex vivo and in vivo (Thompson et al., 2014). The presence of CD8+ T cells in breast cancer is a reliable predictor of clinical outcome and treatment response (Ali et al., 2014; Byrne et al., 2020). By univariable analysis, Denkert et al. found that high tumor-infiltrating lymphocytes predicted longer disease-free survival in patients with HER2-positive breast cancer and TNBC treated with neoadjuvant therapy (Denkert et al., 2018). Our GSEA results showed that multiple immune-related signaling pathways were significantly enriched in the low-risk group. To validate these findings, we performed immune cell infiltration analysis. Obviously, patients with high immune scores had lower risk scores, indicating better prognosis, in line with the findings of Bruni et al. (2020), Jiang et al. (2019), Mlecnik et al. (2016), and Savas et al. (2016). Naïve B cells, CD8+T cells, activated NK cells, and Tregs were significantly enriched in the low-risk group, which was also partially consistent with the aforementioned previous findings. Considering the importance of immune checkpoint inhibitor–based immunotherapies, we further investigated the differences in the expression of 11 immune checkpoints between the high- and low-risk groups. We found that the expression GITR, OX40, CD137, CD40LG, CD28, CD278, CTLA4, VSIR, and CD233 were downregulated in the high-risk group, which corroborated the results of the study of Hu et al. that upregulated immune checkpoint genes were positively associated with high immune infiltration and favorable prognosis in patients with invasive breast carcinoma (Hu et al., 2020). The immune cell infiltration analysis based on our prognostic model suggests that patients in the low-risk group (with high prevalence of tumor-infiltrating cells lymphocytes, elevated immune-related signaling, and high mRNA expression levels of immune checkpoints) may benefit from immune checkpoint inhibitors while patients in the high-risk group (with low prevalence of tumor-infiltrating cells lymphocytes, downregulated immune-related signaling, and low mRNA expression levels of immune checkpoints) may not. Whether and how the PHRLs influence the immune microenvironment of breast cancer still remains to be explored.
Conclusion
In conclusion, we developed a novel prognostic model consisting of 12 hypoxia-related lncRNAs and an integrative nomogram that could predict the OS accurately and effectively for patients with breast cancer. Furthermore, we analyzed the immune cell infiltration conditions and drug sensitivity between high- and low-risk breast cancer classified based on the prognostic model. Our study uncovered dozens of potential prognostic biomarkers and therapeutic targets concerning the hypoxia signaling pathway in breast cancer.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
PG, BY, and XS conceived and designed the study. PG, WWH, YL, and LZ acquired the data and analyzed the data. PG drafted the manuscript. PG, LZ, and WD prepared the tables and figures. PG, LZ, RW, and WW revised the manuscript. All authors read and approved the final manuscript.
Funding
Natural Science Foundation of Shanghai (Grant No.19ZR1441300). National Natural Science Foundation of China (Grant No.81270556). Shanghai Sailing Program (Grant No.20YF1438800).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We sincerely acknowledge TCGA, GTEx, and GEPIA2 databases for providing their platforms and contributors for uploading their meaningful datasets.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.796729/full#supplementary-material
Abbreviations
LncRNA, Long noncoding RNA; TCGA, The Cancer Genome Atlas; GSEA, Gene Set Enrichment Analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; OS, Overall survival; DFS (disease-free survival); ER, Estrogen receptor; PR, Progesterone receptor; HR, Hormone receptor; HER2, Human epidermal growth factor receptor 2; TNBC, Triple-negative breast cancer; LASSO, Least absolute shrinkage and selection operator; ROC, Receiver operating characteristic; AUC, Area under the ROC; HRL, Hypoxia-related lncRNA; PHRLs, Prognostic hypoxia-related lncRNAs; DEHmRNAs, Differentially expressed hypoxia-related mRNAs.
References
Ali, H. R., Provenzano, E., Dawson, S.-J., Blows, F. M., Liu, B., Shah, M., et al. (2014). Association between CD8+ T-Cell Infiltration and Breast Cancer Survival in 12 439 Patients. Ann. Oncol. 25 (8), 1536–1543. doi:10.1093/annonc/mdu191
Ben-Shoshan, J., Maysel-Auslender, S., Mor, A., Keren, G., and George, J. (2008). Hypoxia Controls CD4+CD25+ Regulatory T-Cell Homeostasis via Hypoxia-Inducible Factor-1α. Eur. J. Immunol. 38 (9), 2412–2418. doi:10.1002/eji.200838318
Bin, X., Hongjian, Y., Xiping, Z., Bo, C., Shifeng, Y., and Binbin, T. (2018). Research Progresses in Roles of LncRNA and its Relationships with Breast Cancer. Cancer Cel Int 18, 179. doi:10.1186/s12935-018-0674-0
Bos, R., van der Groep, P., Greijer, A. E., Shvarts, A., Meijer, S., Pinedo, H. M., et al. (2003). Levels of Hypoxia-Inducible Factor-1? Independently Predict Prognosis in Patients with Lymph Node Negative Breast Carcinoma. Cancer 97 (6), 1573–1581. doi:10.1002/cncr.11246
Bray, F., Ferlay, J., Laversanne, M., Brewster, D. H., Gombe Mbalawa, C., Kohler, B., et al. (2015). Cancer Incidence in Five Continents: Inclusion Criteria, Highlights from Volume X and the Global Status of Cancer Registration. Int. J. Cancer 137 (9), 2060–2071. doi:10.1002/ijc.29670
Bristow, R. G., and Hill, R. P. (2008). Hypoxia, DNA Repair and Genetic Instability. Nat. Rev. Cancer 8 (3), 180–192. doi:10.1038/nrc2344
Bruni, D., Angell, H. K., and Galon, J. (2020). The Immune Contexture and Immunoscore in Cancer Prognosis and Therapeutic Efficacy. Nat. Rev. Cancer 20 (11), 662–680. doi:10.1038/s41568-020-0285-7
Byrne, A., Savas, P., Sant, S., Li, R., Virassamy, B., Luen, S. J., et al. (2020). Tissue-resident Memory T Cells in Breast Cancer Control and Immunotherapy Responses. Nat. Rev. Clin. Oncol. 17 (6), 341–348. doi:10.1038/s41571-020-0333-y
Cai, Y., Li, Y., Shi, C., Zhang, Z., Xu, J., and Sun, B. (2021). LncRNA OTUD6B-AS1 Inhibits many Cellular Processes in Colorectal Cancer by Sponging miR-21-5p and Regulating PNRC2. Hum. Exp. Toxicol. 40, 1463–1473. doi:10.1177/0960327121997976
Chen, D., Tang, P., Wang, Y., Wan, F., Long, J., Zhou, J., et al. (2021a). Downregulation of Long Noncoding RNA MR44352HG Suppresses Beast Cancer Progression via the Wnt/βcatenin Signaling Pathway. Oncol. Lett. 21 (5), 373. doi:10.3892/ol.2021.12634
Chen, Y., Zitello, E., Guo, R., and Deng, Y. (2021b). The Function of LncRNAs and Their Role in the Prediction, Diagnosis, and Prognosis of Lung Cancer. Clin. Translational Med. 11 (4), e367. doi:10.1002/ctm2.367
Choudhry, H., Harris, A. L., and McIntyre, A. (2016). The Tumour Hypoxia Induced Non-coding Transcriptome. Mol. Aspects Med. 47-48, 35–53. doi:10.1016/j.mam.2016.01.003
Denkert, C., von Minckwitz, G., Darb-Esfahani, S., Lederer, B., Heppner, B. I., Weber, K. E., et al. (2018). Tumour-infiltrating Lymphocytes and Prognosis in Different Subtypes of Breast Cancer: a Pooled Analysis of 3771 Patients Treated with Neoadjuvant Therapy. Lancet Oncol. 19 (1), 40–50. doi:10.1016/s1470-2045(17)30904-x
Dienstmann, R., Villacampa, G., Sveen, A., Mason, M. J., Niedzwiecki, D., Nesbakken, A., et al. (2019). Relative Contribution of Clinicopathological Variables, Genomic Markers, Transcriptomic Subtyping and Microenvironment Features for Outcome Prediction in Stage II/III Colorectal Cancer. Ann. Oncol. 30 (10), 1622–1629. doi:10.1093/annonc/mdz287
Ferguson, N. L., Bell, J., Heidel, R., Lee, S., Vanmeter, S., Duncan, L., et al. (2013). Prognostic Value of Breast Cancer Subtypes, Ki-67 Proliferation index, Age, and Pathologic Tumor Characteristics on Breast Cancer Survival in Caucasian Women. Breast J. 19 (1), 22–30. doi:10.1111/tbj.12059
Floberg, J. M., Zhang, J., Muhammad, N., DeWees, T. A., Inkman, M., Chen, K., et al. (2021). Standardized Uptake Value for 18F-Fluorodeoxyglucose Is a Marker of Inflammatory State and Immune Infiltrate in Cervical Cancer. Clin. Cancer Res. 27, 4245–4255. doi:10.1158/1078-0432.Ccr-20-4450
Gao, N., Li, Y., Li, J., Gao, Z., Yang, Z., Li, Y., et al. (2020). Long Non-coding RNAs: The Regulatory Mechanisms, Research Strategies, and Future Directions in Cancers. Front. Oncol. 10, 598817. doi:10.3389/fonc.2020.598817
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468
Gil Del Alcazar, C. R., Alečković, M., and Polyak, K. (2020). Immune Escape during Breast Tumor Progression. Cancer Immunol. Res. 8 (4), 422–427. doi:10.1158/2326-6066.Cir-19-0786
Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: the Next Generation. Cell 144 (5), 646–674. doi:10.1016/j.cell.2011.02.013
He, S., Wang, X., Zhang, J., Zhou, F., Li, L., and Han, X. (2020). TRG-AS1 Is a Potent Driver of Oncogenicity of Tongue Squamous Cell Carcinoma through microRNA-543/Yes-Associated Protein 1 axis Regulation. Cell Cycle 19 (15), 1969–1982. doi:10.1080/15384101.2020.1786622
Hsing, C.-H., Cheng, H.-C., Hsu, Y.-H., Chan, C.-H., Yeh, C.-H., Li, C.-F., et al. (2012). Upregulated IL-19 in Breast Cancer Promotes Tumor Progression and Affects Clinical Outcome. Clin. Cancer Res. 18 (3), 713–725. doi:10.1158/1078-0432.Ccr-11-1532
Hu, F.-F., Liu, C.-J., Liu, L.-L., Zhang, Q., and Guo, A.-Y. (2020). Expression Profile of Immune Checkpoint Genes and Their Roles in Predicting Immunotherapy Response. Brief Bioinform 22, bbaa176. doi:10.1093/bib/bbaa176
Jiang, Y.-Z., Ma, D., Suo, C., Shi, J., Xue, M., Hu, X., et al. (2019). Genomic and Transcriptomic Landscape of Triple-Negative Breast Cancers: Subtypes and Treatment Strategies. Cancer Cell 35 (3), 428–440. doi:10.1016/j.ccell.2019.02.001
Jiang, Y., Chen, S., Li, Q., Liang, J., Lin, W., Li, J., et al. (2021). TANK-binding Kinase 1 (TBK1) Serves as a Potential Target for Hepatocellular Carcinoma by Enhancing Tumor Immune Infiltration. Front. Immunol. 12, 612139. doi:10.3389/fimmu.2021.612139
Jiao, Y., Zhou, J., Jin, Y., Yang, Y., Song, M., Zhang, L., et al. (2020). Long Non-coding RNA TDRKH-AS1 Promotes Colorectal Cancer Cell Proliferation and Invasion through the β-Catenin Activated Wnt Signaling Pathway. Front. Oncol. 10, 639. doi:10.3389/fonc.2020.00639
Kong, S., Xue, H., Li, Y., Li, P., Ma, F., Liu, M., et al. (2020). The Long Noncoding RNA OTUD6B-AS1 Enhances Cell Proliferation and the Invasion of Hepatocellular Carcinoma Cells through Modulating GSKIP/Wnt/β-catenin Signalling via the Sequestration of miR-664b-3p. Exp. Cel Res. 395 (1), 112180. doi:10.1016/j.yexcr.2020.112180
Kopp, F., and Mendell, J. T. (2018). Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell 172 (3), 393–407. doi:10.1016/j.cell.2018.01.011
Leithner, K., Wohlkoenig, C., Stacher, E., Lindenmann, J., Hofmann, N. A., Gallé, B., et al. (2014). Hypoxia Increases Membrane Metallo-Endopeptidase Expression in a Novel Lung Cancer Ex Vivo Model - Role of Tumor Stroma Cells. BMC Cancer 14, 40. doi:10.1186/1471-2407-14-40
Li, N., and Zhan, X. (2019). Identification of Clinical Trait-Related lncRNA and mRNA Biomarkers with Weighted Gene Co-expression Network Analysis as Useful Tool for Personalized Medicine in Ovarian Cancer. EPMA J. 10 (3), 273–290. doi:10.1007/s13167-019-00175-0
Li, X., Li, Y., Yu, X., and Jin, F. (2020). Identification and Validation of Stemness-Related lncRNA Prognostic Signature for Breast Cancer. J. Transl Med. 18 (1), 331. doi:10.1186/s12967-020-02497-4
Liang, Y., Song, X., Li, Y., Chen, B., Zhao, W., Wang, L., et al. (2020). LncRNA BCRT1 Promotes Breast Cancer Progression by Targeting miR-1303/PTBP3 axis. Mol. Cancer 19 (1), 85. doi:10.1186/s12943-020-01206-5
Liu, S. J., Dang, H. X., Lim, D. A., Feng, F. Y., and Maher, C. A. (2021). Long Noncoding RNAs in Cancer Metastasis. Nat. Rev. Cancer 21, 446–460. doi:10.1038/s41568-021-00353-1
Loibl, S., Poortmans, P., Morrow, M., Denkert, C., and Curigliano, G. (2021a). Breast Cancer. Lancet. doi:10.1016/s0140-6736(20)32381-3
Loibl, S., Poortmans, P., Morrow, M., Denkert, C., and Curigliano, G. (2021b). Breast Cancer. The Lancet 397 (10286), 1750–1769. doi:10.1016/s0140-6736(20)32381-3
Lv, Y., Wei, W., Huang, Z., Chen, Z., Fang, Y., Pan, L., et al. (2018). Long Non-coding RNA Expression Profile Can Predict Early Recurrence in Hepatocellular Carcinoma after Curative Resection. Hepatol. Res. 48 (13), 1140–1148. doi:10.1111/hepr.13220
Ma, W., Zhao, F., Yu, X., Guan, S., Suo, H., Tao, Z., et al. (2020). Immune-related lncRNAs as Predictors of Survival in Breast Cancer: a Prognostic Signature. J. Transl Med. 18 (1), 442. doi:10.1186/s12967-020-02522-6
Marin-Acevedo, J. A., Dholaria, B., Soyano, A. E., Knutson, K. L., Chumsri, S., and Lou, Y. (2018). Next Generation of Immune Checkpoint Therapy in Cancer: New Developments and Challenges. J. Hematol. Oncol. 11 (1), 39. doi:10.1186/s13045-018-0582-8
Marusyk, A., Almendro, V., and Polyak, K. (2012). Intra-tumour Heterogeneity: a Looking Glass for Cancer? Nat. Rev. Cancer 12 (5), 323–334. doi:10.1038/nrc3261
Mlecnik, B., Bindea, G., Kirilovsky, A., Angell, H. K., Obenauf, A. C., Tosolini, M., et al. (2016). The Tumor Microenvironment and Immunoscore Are Critical Determinants of Dissemination to Distant Metastasis. Sci. Transl. Med. 8 (327). 327ra326. doi:10.1126/scitranslmed.aad6352
Niu, Y., Bao, L., Chen, Y., Wang, C., Luo, M., Zhang, B., et al. (2020). HIF2-Induced Long Noncoding RNA RAB11B-AS1 Promotes Hypoxia-Mediated Angiogenesis and Breast Cancer Metastasis. Cancer Res. 80 (5), 964–975. doi:10.1158/0008-5472.Can-19-1532
Noman, M. Z., Desantis, G., Janji, B., Hasmim, M., Karray, S., Dessen, P., et al. (2014). PD-L1 Is a Novel Direct Target of HIF-1α, and its Blockade under Hypoxia Enhanced MDSC-Mediated T Cell Activation. J. Exp. Med. 211 (5), 781–790. doi:10.1084/jem.20131916
O’Brien, K. M., Mooney, T., Fitzpatrick, P., and Sharp, L. (2018). Screening Status, Tumour Subtype, and Breast Cancer Survival: a National Population-Based Analysis. Breast Cancer Res. Treat. 172 (1), 133–142. doi:10.1007/s10549-018-4877-9
Pardoll, D. M. (2012). The Blockade of Immune Checkpoints in Cancer Immunotherapy. Nat. Rev. Cancer 12 (4), 252–264. doi:10.1038/nrc3239
Pei, J., Li, Y., Su, T., Zhang, Q., He, X., Tao, D., et al. (2020). Identification and Validation of an Immunological Expression-Based Prognostic Signature in Breast Cancer. Front. Genet. 11, 912. doi:10.3389/fgene.2020.00912
Qian, Y., Shi, L., and Luo, Z. (2020). Long Non-coding RNAs in Cancer: Implications for Diagnosis, Prognosis, and Therapy. Front. Med. 7, 612393. doi:10.3389/fmed.2020.612393
Samanta, D., and Semenza, G. L. (2018). Metabolic Adaptation of Cancer and Immune Cells Mediated by Hypoxia-Inducible Factors. Biochim. Biophys. Acta (Bba) - Rev. Cancer 1870 (1), 15–22. doi:10.1016/j.bbcan.2018.07.002
Savas, P., Salgado, R., Denkert, C., Sotiriou, C., Darcy, P. K., Smyth, M. J., et al. (2016). Clinical Relevance of Host Immunity in Breast Cancer: from TILs to the Clinic. Nat. Rev. Clin. Oncol. 13 (4), 228–241. doi:10.1038/nrclinonc.2015.215
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Sui, S., An, X., Xu, C., Li, Z., Hua, Y., Huang, G., et al. (2020). An Immune Cell Infiltration-Based Immune Score Model Predicts Prognosis and Chemotherapy Effects in Breast Cancer. Theranostics 10 (26), 11938–11949. doi:10.7150/thno.49451
Sun, X., Qian, Y., Wang, X., Cao, R., Zhang, J., Chen, W., et al. (2020a). LncRNA TRG-AS1 Stimulates Hepatocellular Carcinoma Progression by Sponging miR-4500 to Modulate BACH1. Cancer Cel Int 20, 367. doi:10.1186/s12935-020-01440-3
Sun, Z., Jing, C., Xiao, C., and Li, T. (2020b). Long Non-coding RNA Profile Study Identifies an Immune-Related lncRNA Prognostic Signature for Kidney Renal Clear Cell Carcinoma. Front. Oncol. 10, 1430. doi:10.3389/fonc.2020.01430
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., 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, 209–249. doi:10.3322/caac.21660
Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING Database in 2017: Quality-Controlled Protein-Protein Association Networks, Made Broadly Accessible. Nucleic Acids Res. 45 (D1), D362–d368. doi:10.1093/nar/gkw937
Thompson, A. A. R., Elks, P. M., Marriott, H. M., Eamsamarng, S., Higgins, K. R., Lewis, A., et al. (2014). Hypoxia-inducible Factor 2α Regulates Key Neutrophil Functions in Humans, Mice, and Zebrafish. Blood 123 (3), 366–376. doi:10.1182/blood-2013-05-500207
Vaupel, P., Schlenger, K., Knoop, C., and Höckel, M. (1991). Oxygenation of Human Tumors: Evaluation of Tissue Oxygen Distribution in Breast Cancers by Computerized O2 Tension Measurements. Cancer Res. 51 (12), 3316–3322.
Vaupel, P., Höckel, M., and Mayer, A. (2007). Detection and Characterization of Tumor Hypoxia Using pO2 Histography. Antioxid. Redox Signaling 9 (8), 1221–1236. doi:10.1089/ars.2007.1628
Vaupel, P., Mayer, A., and Höckel, M. (2004). Tumor Hypoxia and Malignant Progression. Methods Enzymol. 381, 335–354. doi:10.1016/s0076-6879(04)81023-1
Wang, G. L., Jiang, B. H., Rue, E. A., and Semenza, G. L. (1995). Hypoxia-inducible Factor 1 Is a basic-helix-loop-helix-PAS Heterodimer Regulated by Cellular O2 Tension. Proc. Natl. Acad. Sci. 92 (12), 5510–5514. doi:10.1073/pnas.92.12.5510
Wang, G., Zhang, Z.-j., Jian, W.-g., Liu, P.-h., Xue, W., Wang, T.-d., et al. (2019a). Novel Long Noncoding RNA OTUD6B-AS1 Indicates Poor Prognosis and Inhibits clear Cell Renal Cell Carcinoma Proliferation via the Wnt/β-Catenin Signaling Pathway. Mol. Cancer 18 (1), 15. doi:10.1186/s12943-019-0942-1
Wang, L., Zhao, H., Xu, Y., Li, J., Deng, C., Deng, Y., et al. (2019b). Systematic Identification of lincRNA‐based Prognostic Biomarkers by Integrating lincRNA Expression and Copy Number Variation in Lung Adenocarcinoma. Int. J. Cancer 144 (7), 1723–1734. doi:10.1002/ijc.31865
Wang, W., Cheng, X., and Zhu, J. (2021). Long Non-coding RNA OTUD6B-AS1 Overexpression Inhibits the Proliferation, Invasion and Migration of Colorectal Cancer Cells via Downregulation of microRNA-3171. Oncol. Lett. 21 (3), 193. doi:10.3892/ol.2021.12454
Wang, Y., Roche, O., Yan, M. S., Finak, G., Evans, A. J., Metcalf, J. L., et al. (2009). Regulation of Endocytosis via the Oxygen-Sensing Pathway. Nat. Med. 15 (3), 319–324. doi:10.1038/nm.1922
Wang, Y., Yang, T., Han, Y., Ren, Z., Zou, J., Liu, J., et al. (2020a). lncRNA OTUD6B-AS1 Exacerbates As2O3-Induced Oxidative Damage in Bladder Cancer via miR-6734-5p-Mediated Functional Inhibition of IDH2. Oxidative Med. Cell Longevity 2020, 1–22. doi:10.1155/2020/3035624
Wang, Z., Xia, F., Feng, T., Jiang, B., Wang, W., and Li, X. (2020b). OTUD6B-AS1 Inhibits Viability, Migration, and Invasion of Thyroid Carcinoma by Targeting miR-183-5p and miR-21. Front. Endocrinol. 11, 136. doi:10.3389/fendo.2020.00136
Xia, X., and Kung, A. L. (2009). Preferential Binding of HIF-1 to Transcriptionally Active Loci Determines Cell-type Specific Response to Hypoxia. Genome Biol. 10 (10), R113. doi:10.1186/gb-2009-10-10-r113
Xie, H., Shi, S., Chen, Q., and Chen, Z. (2019). LncRNA TRG-AS1 Promotes Glioblastoma Cell Proliferation by Competitively Binding with miR-877-5p to Regulate SUZ12 Expression. Pathol. - Res. Pract. 215 (8), 152476. doi:10.1016/j.prp.2019.152476
Yang, F., Zhang, H., Mei, Y., and Wu, M. (2014). Reciprocal Regulation of HIF-1α and LincRNA-P21 Modulates the Warburg Effect. Mol. Cel 53 (1), 88–100. doi:10.1016/j.molcel.2013.11.004
Yeo, J. G., Wasser, M., Kumar, P., Pan, L., Poh, S. L., Ally, F., et al. (2020). The Extended Polydimensional Immunome Characterization (EPIC) Web-Based Reference and Discovery Tool for Cytometry Data. Nat. Biotechnol. 38 (6), 679–684. doi:10.1038/s41587-020-0532-1
Yuan, M., Yu, C., Chen, X., and Wu, Y. (2021). Investigation on Potential Correlation between Small Nuclear Ribonucleoprotein Polypeptide A and Lung Cancer. Front. Genet. 11, 1842. doi:10.3389/fgene.2020.610704
Keywords: breast cancer, hypoxia, long noncoding RNA, prognosis, immune infiltration, nomogram
Citation: Gu P, Zhang L, Wang R, Ding W, Wang W, Liu Y, Wang W, Li Z, Yan B and Sun X (2021) Development and Validation of a Novel Hypoxia-Related Long Noncoding RNA Model With Regard to Prognosis and Immune Features in Breast Cancer. Front. Cell Dev. Biol. 9:796729. doi: 10.3389/fcell.2021.796729
Received: 17 October 2021; Accepted: 30 November 2021;
Published: 16 December 2021.
Edited by:
Zong Sheng Guo, Roswell Park Comprehensive Cancer Center, United StatesReviewed by:
Liangyou Gu, People’s Liberation Army General Hospital, ChinaGuo-Jiang Zhao, Capital Medical University, China
Copyright © 2021 Gu, Zhang, Wang, Ding, Wang, Liu, Wang, Li, Yan and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Bin Yan, yanbin8988@sjtu.edu.cn; Xing Sun, xingsun_sjtu@aliyun.com
†These authors have contributed equally to this work