Identification of prognostic and immunotherapy-related eRNA ID2-AS1 in bladder cancer

AbstractEnhancer RNAs (eRNAs) can participate in enhancer regulation and target gene transcription, thus affecting the occurrence and development of tumors. In this study, we identified eRNAs closely related to bladder cancer (BLCA). Gene expression profiles and clinical information from The Cancer Genome Atlas (TCGA) database were used in this study. The Atlas of Noncoding RNAs in Cancer (TANRIC) co-expression data was also studied to evaluate correlations between the inferred levels of eRNA and its predicted target genes. Moreover, we evaluated differences in tumor microenvironment between high and low ID2-AS1 expression groups, and predicted the response of high- and low-expression groups to immune checkpoint inhibitor (ICI) treatment. Finally, we analyzed the prognostic value of ID2-AS1 in different tumors. ID2-AS1 and ID2 were identified as eRNAs and target genes related to the prognosis of BLCA. Low ID2-AS1 levels were associated with advanced age, low overall survival, high histological grade, and late BLCA staging. ID2-AS1 appeared to regulate epithelial mesenchymal transition, mitotic spindle assembly, and angiogenesis, thereby affecting BLCA progression. The ID2-AS1 high-expression group had better ICI treatment response. In addition, ID2-AS1 also had prognostic value in other cancers. ID2-AS1 helps predict prognostic and immunotherapeutic effects in BLCA.


Introduction
Bladder carcinoma (BLCA) is the second-most common urological malignancy worldwide, with an estimated 83,730 new cases and 17,200 deaths in 2021. [1] Ninety percent of BLCAs are caused by malignant transformation of urothelial cells. Approximately 25% of patients with BLCA have muscle-invasive bladder cancer (MIBC) or metastasis, and 75% have nonmuscle invasive bladder cancer (NMIBC). [2] If not diagnosed and treated early, NMIBC cells are likely to exhibit muscle invasion and metastasis. As a result, approximately 2/3 of advanced BLCA patients present with MIBC at their first clinical visit. [3] Half of these patients are reported to develop distant metastases within 2 years. [4] Once a tumor has developed to the metastatic stage, surgery combined with general chemotherapy are inadequate for treatment, resulting in a median survival of 14 months. [3,5] According to research, the response rate for novel immunotherapy and immune checkpoint inhibitors (ICIs) is 30% or less. [6] Therefore, there remains an urgent clinical demand for the discovery of novel molecular biomarkers by molecular profiling of BLCA cases.
With the development of high-throughput gene detection technology and the establishment of large-scale gene expression data sets, greater attention is being paid to long noncoding RNAs linked to BLCA. [7][8][9] The expression of lncRNAs in tissues and cell types has distinct specificity, making them potential cancer biomarkers. [7] Among these lncRNAs, enhancer RNAs (eRNAs) have attracted considerable attention because of their functions in mediating enhancer and gene transcription, and frequent overlap with noncoding risk gene loci associated with disease. [10][11][12] It has been shown in BLCA that downregulation of P2RY2e diminishes cancer promotion, which thus represents a potential therapeutic target. [13] Estrogen might promote the development of BLCA by inducing SMAD7e production. [14] Medicine Furthermore, high expression of MARC1 eRNA is associated with malignant cell behavior. [15] To the best of our knowledge, only 3 studies have been conducted regarding BLCA-associated eRNAs.
In this study, we screened differentially expressed eRNAs (DEEs) in normal and tumor tissues, and explored prognostic eRNAs. ID2-AS1 was found to be significantly associated with survival of patients with BLCA. High expression of ID2-AS1 resulted in elevated purity of tumor cells in the tumor microenvironment (TME), which suggested that ID2-AS1 affects the cellular composition of the TME and thus affects the progression and treatment of BLCA. [16,17] 2. Materials and methods

TCGA data analysis
We obtained RNA expression data for various cancers from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). At the same time, we also obtained clinical data (age, gender, pathological stage, etc.) and case survival information from TCGA database. To explore possible target genes regulated by eRNAs, we identified 1585 eRNAs and 2303 predicted target genes using the PreSTIGE algorithm. We also studied TANRIC co-expression data (https://bioinformatics.mdanderson.org/public-software/tanric) to assess correlations between the predicted levels of eRNA and predicted target gene expression. Among these eRNAs, only those that were significantly correlated with the overall survival (OS) rate (P < .05) and target gene expression levels (R > .4, P < .001) were considered as candidate key eRNAs in BLCA. Subsequently, we analyzed the relationship between the expression of eRNA and clinicopathological features of BLCA. In addition, we validated the prognostic value of the selected eRNA at the pan-cancer level in all 32 types of malignancies. All datasets in the present study were downloaded from public databases, thus ethical approval was not required.

Gene bioinformatics enrichment analysis
Gene set variation analysis (GSVA) uses standardized RNA-seq data from TCGA database. As an unsupervised gene set enrichment method, GSVA quantifies the pathway level through gene level quantification to analyze the difference in pathway levels between the ID2-AS1 low-expression group and ID2-AS1 high-expression cohort. We downloaded the "Hallmark" gene sets from the Molecular Signatures Database (http://software. broadinstitute.org/gsea/index.jsp) for GSVA through the package "GSVA" in R. The "MCPcounter" package in R was used for analysis of microenvironment cell populations (MCPs) and quantification of immune cells from transcriptomic data. The results were delivered by the package "pheatmap" in R. Subsequently, the package "limma" in R was used to identify pathways with the most significant differences between patients in the signature group.

Evaluation of immune cell infiltration in TME and immune checkpoints Inhibitor responses
We evaluated the enrichment of immune and stromal cell gene markers. We used the ESTIMATE algorithm to calculate immunity and matrix scores. In addition, we chose a new TME scoring algorithm to explore differences between high-and low-risk groups based on Immune-related gene pairs. The CIBERSORT algorithm was used to infer the relative proportion of 22 infiltrating immune cells in each sample. In addition, the tumor immune dysfunction and rejection (TIDE) algorithm was used to evaluate ICI treatment responses.

Statistical analysis
All analyses in this study were performed using R software (v.3.6.3); P < .05 was considered statistically significant.

Identification of DEEs and prognostic eRNAs
A total of 1585 eRNAs and 2303 predicted target genes were identified using the PreSTIGE algorithm. [18] We extracted eRNA expression data from TCGA-BLCA dataset and found 379 DEEs. Next, according to TCGA-BLCA patient clinical prognostic information, we filtered out 53 of 379 DEEs that were significantly correlated with overall survival (Table 1, Kaplan-Meier log-rank test, P < .05). ID2-AS1 was closely related to the prognosis of bladder cancer (P = .003), and positively correlated with ID2 (R = .71). In addition, ID2 has been shown to play an important role in a variety of tumors. So we chose ID2-AS1 in this study.

LncRNA ID2-AS1 is a key eRNA in BLCA
In tumors, the expression level of lncRNA ID2-AS1 was lower than that in normal tissues (Fig. 1A). In addition, in patients with BLCA, the ID2-AS1 high-expression group showed better overall survival than the low-expression group (Fig. 1B). There is a strong correlation between ID2 level and predicted target gene ID2-AS1 expression (Fig. 1C). Moreover, the summary of clinical characteristics of TCGA-BLCA was shown in Table 2. ID2-AS1 expression analysis was performed for different clinical feature subclasses (age, gender, WHO grade, American Joint Committee on Cancer (AJCC) stage, and TNM stage). We observed that patients with advanced BLCA (high WHO grade, high AJCC stage, and MIBC) had lower expression levels ( Fig. 2A-G). Subsequently, OS curves using combinations of different clinical features and ID2-AS1 expression levels showed that lower expression levels were associated with poorer OS in patients aged > 60 years, AJCC-Stage I-II, T0-2 stage, and N0 stage subgroups ( Figure 3A-D).

Multiple aggressive, immunosuppressive, and angiogenic hallmarks related to low ID2-AS1 expression
To investigate differentially expressed features between ID2-AS1 high-and low-expression groups, we performed GSVA of the TCGA-BLCA cohort. Notably, gene sets related to tumor necrosis factor-α signaling via nuclear factor-κB, Myc targets, mTOR complex 1 signaling, epithelial mesenchymal transition, mitotic spindle, and angiogenesis signaling pathways were upregulated in the ID2-AS1 low-expression group (Fig. 4A-B).

Differences in TME between ID2-AS1 high-and lowexpression groups, and prediction of whether ID2-AS1 expression is correlated with ICI therapeutic responses
The ESTIMATE algorithm was used to evaluate tumor purity for each patient. There was a significant difference in tumor purity between the high-and low-expression groups (P <.001, Figure 5A). Furthermore, the TIDE algorithm, which was established to predict immunotherapy responses through transcriptomic data, was used to explore whether the expression level of ID2-AS1 could predict immunotherapeutic benefit in the TCGA-BLCA cohort. The results for each patient are shown in Figure 5B. These results revealed that the percentage of ICI therapy responders was significantly higher for the high-expression group (44%) compared to the low ID2-AS1 expression group (30%) (Chi-squared test, P =.00367; Figure 5B). ID2-AS1 expression was robustly positively correlated with ICI therapeutic response in BLCA patients.
To further characterize the levels of different immune cell subsets in the TME and chemotactic signals involved, we employed the microenvironment cell population-counter (MCP -counter) method, which allows for robust quantification of the absolute abundance of 8 immune and 2 stromal cell populations in heterogeneous tissues from transcriptomic data [19] (Figure 5C). As shown in Figure 5D, compared to the ID2-AS1 low-expression group, cytotoxic T-lymphocytes and T cells were upregulated in the high-expression group (P < .001 and P < .05). Strikingly, fibroblasts, a type of stromal cell, were downregulated in the high-expression group (P < .001).

Discussion
In recent years, studies have found that eRNAs regulate gene expression by promoting the enhancer promoter chain, and knockout of eRNAs is usually accompanied by a decrease in target gene expression. A large body of evidence shows that abnormal expression of eRNAs in various cancers is closely related to tumorigenesis. In various types of cancer, eRNAs may modulate the occurrence and progression of tumor cells by regulating certain signaling pathways. However, few researchers have examined BLCA prognosis-related eRNAs. Our study provides insights into the potential role of ID2-AS1 in tumor immunology and its use as a new BLCA-related molecular biomarker.
In our study, we searched for eRNAs related to BLCA prognosis. Through a series of comprehensive analyses, we found that ID2-AS1 and its corresponding target gene, ID2, had a significant impact on survival of BLCA patients. Therefore, we regarded ID2-AS1 as a primary eRNA in BLCA. The results showed that low expression of ID2-AS1 represented an independent prognostic factor in patients with BLCA. We further evaluated the clinical characteristics and relationships between ID2-AS1 expression and BLCA. Interestingly, low expression of ID2-AS1 predicted advanced tumor staging better, and the same results were also found in different clinical subgroup analyses. As an important member of the DNA-binding inhibitor family, ID2 has been shown in previous studies to play key roles in the growth and invasion of various tumor cells. [20,21] Upregulation of ID2 indicates poor prognosis of BRCA patients after breast-conserving surgery. [21] In head and neck squamous cell carcinoma, ID2 is a new molecule involved in the regulation of tumor cell invasion and lymph node metastasis. [20] These findings coincide with the results of our pan-cancer analysis.
Subsequently, the results of GSVA revealed that there were significant differences in biological behavior between the ID2-AS1 high-expression group and the ID2-AS1 low-expression group. Compared to the high-expression group, the low-expression cohort scored higher in epithelial mesenchymal transformation, mitotic spindle formation, and angiogenesis signaling pathways. There is evidence that these signals are closely related to the occurrence and development of tumors. In a study by Zhou et      . Gene expression differences for ID2-AS1 were analyzed in 24 different cancer types. ns: P ≥ .05, *P < .05, **P < .01, ***P < .001, ****P < .0001. www.md-journal.com al, ID2-AS1 and ID2 inhibited twist-induced epithelial mesenchymal transformation in hepatocellular carcinoma cells to inhibit tumor metastasis, [22] which is consistent with our results. Tumor tissues are mainly composed of tumor cells, stromal cells, and immune cells. Tumor purity represents the proportion of tumor cells in the tumor tissue. Interestingly, in studies of glioma, gastric cancer, and colon cancer, patients with low tumor purity tended to have higher malignancy and poorer prognosis. [23][24][25] In our study, the ID2-AS1 high expression group had a higher tumor density and better OS. TME immune cell infiltration is considered an important and unpredictable phenomenon to predict the prognosis of various cancers and immunotherapeutic responses. The difference analysis of immune cell infiltration showed that the infiltration abundance of CTLs and T cells in the high expression group was significantly higher, and the infiltration abundance of fibroblasts in the low expression group was higher, which is consistent with the results of GSVA to a certain extent. We used the TIDE algorithm to predict the response of patients with BLCA to immunotherapy. It was encouraging that patients in the high expression group are more likely to benefit from immunotherapy.
The significant advantage of our research is the comprehensive analysis conducted, which revealed that ID2-AS1 is an independent prognostic factor for BLCA. We also predicted the relationship between high expression of ID2-AS1 and positive ICI treatment, and found that it is associated with OS in many tumors. Our research has potential value in future tumor prognosis and immunotherapy. However, we did not verify the prognostic ability of ID2-AS1 in the external database; therefore, the predictive value of ID2-AS1 in terms of immunotherapeutic responses needs to be verified in future experiments. Other researchers can use this study as a cornerstone to further reveal the cellular mechanisms involving ID2-AS1 in the occurrence and development of various tumors.

Conclusions
In conclusion, ID2-AS1 is an independent protective factor in BLCA, and may become a potentially powerful tool to improve the diagnosis and prognosis of BLCA patients. Hence, ID2-AS1 is a potential promising immune-related therapeutic target for the precise treatment of a variety of malignant tumors in the future.