Analysis of m6A-Related Signatures in the Tumor Immune Microenvironment and Identification of Clinical Prognostic Regulators in Adrenocortical Carcinoma

Adrenocortical carcinoma (ACC) is a rare endocrine malignancy with a high rate of mortality and recurrence. N6-methyladenosine methylation (m6A) is the most common modification to affect cancer development, but to date, the potential role of m6A regulators in ACC prognosis is not well understood. In this study, we systematically analyzed 21 m6A regulators in ACC samples from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database. We identified three m6A modification patterns with different clinical outcomes and discovered a significant relationship between diverse m6A clusters and the tumor immune microenvironment (immune cell types and ESTIMATE algorithm). Additionally, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) revealed that the m6A clusters were strongly associated with immune infiltration in the ACC. Next, to further explore the m6A prognostic signatures in ACC, we implemented Lasso (Least Absolute Shrinkage and Selection Operator) Cox regression to establish an eight-m6A-regulator prognostic model in the TCGA dataset, and the results showed that the model-based high-risk group was closely correlated with poor overall survival (OS) compared with the low-risk group. Subsequently, we validated the key modifications in the GEO datasets and found that high HNRNPA2B1 expression resulted in poor OS and event-free survival (EFS) in ACC. Moreover, to further decipher the molecular mechanisms, we constructed a competing endogenous RNA (ceRNA) network based on HNRNPA2B1, which consists of 12 long noncoding RNAs (lncRNAs) and 1 microRNA (miRNA). In conclusion, our findings indicate the potential role of m6A modification in ACC, providing novel insights into ACC prognosis and guiding effective immunotherapy.

Introduction

Adrenocortical carcinoma (ACC) is an uncommon endocrine malignancy with an annual incidence of 0.7–2.0 cases per million (1). Despite its rarity, the 5-year survival rate in most series is <35% (2). Currently, the only curative therapy for localized ACC is surgery. Even with complete excision, rates of local recurrence typically range from 19–34% (3, 4). Adjuvant treatments which aim to decrease recurrence, including chemotherapy and radiotherapy, show limited therapeutic effectiveness (5). Nevertheless, the most widely used tumor, lymph node, and metastasis (TNM) classification system remains unacceptable for heterogeneous outcomes and poor survival (6, 7). Therefore, unraveling the genomic properties underlying ACC is crucial for developing effective treatments and predicting individual survival and recurrence risk.

N6-methyladenosine (m6A), which was first discovered in the 1970s, is recognized as the most prominent and abundant form of internal modification that occurs in messenger RNAs (mRNAs) and long non-coding RNAs (lncRNAs) in many eukaryotic species (8, 9). m6A methylation is thought to affect every aspect of RNA metabolism, including RNA splicing, translocation, stability, and translation into protein (10). The m6A modification is dynamically deposited by three types of homologous factors: methyltransferases (“writers”), demethylases (“erasers”), and m6A binding proteins (“readers”) (11). Methyltransferases, with core members, METTL3, METTL14, WTAP, ZC3H13, RBM15, and RBM15B, catalyze the methyl group directly attached to the nitrogen on the sixth carbon of the aromatic ring of an adenosine residue (12). Demethylases, which mediate the m6A removal process, with core members FTO and ALKBH5, selectively remove the methyl code from specific mRNAs (11). m6A binding proteins, including the YTHDF family (YTHDF1/2/3), nuclear heterogeneous riboprotein family (HNRNPA2B1 and HNRNPC), and eukaryotic initiation factor (eIF or EIF1A), aim to decode RNA methylation and recognize the m6A motif (11).

To date, N6-methyladenosine, a potential biomarker, has been reported to actively participate in various important physiological processes such as stem cell differentiation, circadian periods, and DNA damage response in vivo (13–15). Aberrant expression and mutation of m6A were confirmed to result in the abnormal processes, including dysregulation of cell death and proliferation, developmental defects, and impaired self-renewal capacity (16, 17). Recent studies have demonstrated that abnormal m6A methylation modification is closely associated with a variety of human diseases, especially cancer, including bladder cancer, head and neck squamous cell carcinoma, gastric cancer, breast cancer, hepatocellular carcinoma, and colorectal cancer (18). For instance, in breast cancer, high FTO levels are significantly associated with poor survival rates. Furthermore, in a series of in vitro and in vivo assays, FTO dramatically alleviated and degraded BNIP3 (a tumor suppressor) via a YTHDF2-independent mechanism to induce cancer cell proliferation, colony formation, and metastasis (19). Additionally, Zewei Tu reportedly built a ceRNA network and established a 9 m6A-related lncRNA prognostic model in lower-grade glioma patients (20). Accumulating evidence has shown that m6A-related mRNAs and lncRNAs can serve as novel potential targets for predicting prognosis and developing personalized treatments for many types of cancer. However, little is known about the relationship between the effect of m6A methylation modification and ACC.

The tumor microenvironment (TME), which includes cancer cells, stromal cells, and distant recruited cells, such as infiltrating immune cells (myeloid cells and lymphocytes), bone marrow-derived cells (BMDCs), and secreted factors such as cytokines and chemokines, play a crucial role in tumor progression and affect the clinical benefit from novel strategies of immunological checkpoint blockade (ICB, PD-1/L1, and CTLA-4) (21, 22). In advanced ACC, no investigated therapy has offered long-term disease control, except for immune checkpoint blockade. A phase II study indicated that pembrolizumab (an anti-PD-1 monoclonal antibody) can provide meaningful and durable antitumor activity (23). Emerging studies have focused on research interests that enhance the in-depth understanding of the heterogeneity and complexity of the TME to improve immunotherapy strategies by comprehensive analysis of particular m6A regulators (24, 25). For instance, inhibiting ALKBH5 may enhance the efficacy of anti-PD-1 therapy in melanoma patients by mediating the level of Mct4/Slc16a3, which is involved in regulating suppressive lymphocyte Treg and myeloid-derived suppressor cell accumulations in the TME (24). Currently, none of these studies have extended their research into the frontiers of knowledge in ACC.

In this study, we performed a retrospective analysis based on The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases to estimate the effect of m6A-related genes on the prognostic value. We identified multiple m6A regulators and related lncRNAs or microRNAs (miRNAs) as potential biomarkers by shaping individual TME characterizations.

Materials and Methods

Data Processing of the ACC Dataset

The public RNA sequencing, mutation expression, and full clinical information of ACC were downloaded from TCGA and GEO. Patients without survival information were excluded from further evaluation. The RNA sequencing data (FPKM value) and somatic mutation data from TCGA-ACC (The Cancer Genome Atlas - Adrenocortical carcinoma) were downloaded from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/) and gathered as a training set for further analysis. In total, six eligible data from GEO (GSE10927, GSE19750, GSE33371, GSE76019, GSE76021, and GSE49280) were downloaded and an averaging method with the affy and simpleaffy packages was used to perform background adjustment and quantile normalization.

Consensus Clustering of m6A Regulators

We first selected 21 m6A RNA methylation regulators from previously published articles (26, 27). These 21 m6A regulators included 8 writers (METTL3, METTL14, RBM15, RBM15B, WTAP, KIAA1429, CBLL1, and ZC3H13), 2 erasers (ALKBH5, FTO), and 11 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, and ELAVL1). Based on the expression of the 21 m6A modulators, the patients were classified into three groups using the optimal k-means clustering (“kmeans” function in R). Cluster analysis was performed using the ConsensusClusterPlus R package with cycle computation 1,000 times to ensure stability and reliability (28). The overall survival (OS) between different clusters was calculated using the Kaplan-Meier method.

Identification of Differentially Expressed Genes (DEGs) Between m6A Patterns

To identify the DEGs between three clusters in the TCGA-ACC cohort, the empirical Bayesian approach of the limma R package was applied in the standard comparison mode. The significance criteria for determining DEGs was set to | logFC | > 1 and P-value < 0.05. To investigate the pathways enriched in the different subgroups, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) biological processes by applying a threshold P-value < 0.05, minimum count of 5, and enrichment factor > 0.15. Gene set enrichment analysis (GSEA) was used to evaluate all genes based on their log2 fold change and assess the functions associated with subtypes by implementing the clusterProfiler R package.

Comparison of Immune Cell Infiltration Among m6A Patterns

To explore the degree of immune cell infiltration among the three subgroups, we applied the ESTIMATE algorithm, in which R script was downloaded from the website (https://sourceforge.net/projects/estimateproject/) to calculate the estimate scores, immune scores, and stromal scores for further predicting tumor purity and analyzing the TME (29). To explore the differences in immune cell subtypes among multiple clusters, we utilized the CIBERSORT package to assess the proportions of 22 immune cell subtypes based on TCGA-ACC samples. The results with P < 0.05 in CIBERSORT analysis were used for further analysis. The Mann-Whitney U test was used to compare the differences among the three subgroups.

Least Absolute Shrinkage and Selection Operator Cox Regression and Validation of the Prognostic m6A Signatures

To enhance the prediction accuracy and interpretability of the statistical model, Lasso Cox regression analysis was carried out to examine the relationship between m6A prognosis signatures and ACC risk. Using the “glmnet” software package of R, 8 m6A-related genes were screened to construct the best prognostic model. The risk score was generated using the following formula: risk score = Expression_mRNA1 × Coefficient_mRNA1 + Expression_mRNA2 × Coefficient_mRNA2 +…Expression_mRNAn × Coefficient_mRNAn. According to the predictive model, the patients were divided into high-risk and low-risk groups using the median cutoff of risk score. The Cox proportional hazard regression model includes age and TNM stage. The hazard ratio (HR) from Cox regression analysis was used to distinguish the prognostic factors positively or negatively. A gene with HR > 1 was considered a risk gene, and a gene with HR < 1 was considered a protective gene. Subsequently, the Kaplan-Meier survival method was used to evaluate the availability of the prognostic model, and the sensitivity and specificity of the receiver operating characteristic (ROC) curve were used to evaluate the prognostic accuracy of the signature building. Similarly, the five validated GEO cohorts were further calculated to validate the prognostic value of selected m6A-related genes by OS and PFS (progression-free survival) analyses.

Gene Mutation Screening and Analysis

The SNP dataset was based on VarScan2 variant aggregation and masking data in TCGA. Here, we analyzed the SNP mutation and carried out visualization using the R maftools package (30). The SNP2APA database was designed to explore the effects of single nucleotide polymorphisms (SNPs) and provided OS across different cancer types (31). Based on this database, we searched for specific SNPs that have prognostic value in ACC, and assessed the relationship between these SNPs and key m6A-related genes.

Construction of the ceRNA Network

Using the TCGA cohort, differentially expressed genes between m6A patterns were identified with the standards of | log2(Fold change) | > 1 and P < 0.05 using the R package “limma.” Perl programming language was applied to target miRNAs-lncRNAs and miRNAs-mRNAs in the prediction analysis. Furthermore, miRcode was used to collect and target experimentally validated lncRNAs (32). StarBase v3.0 was used to predict miRNA-mRNA interactions (http://starbase.sysu.edu.cn/) (33). The ceRNA network was visualized using the “Cytoscape” software (34).

Statistical Analysis

Most analyses were performed using R software (version 3.6.1, http://www.R-project.org). Kaplan-Meier curves and the log-rank test were used to compare the OS between various subgroups based on the expression of m6A-related genes. Univariate and multivariate Cox proportional hazard regression analyses were used to evaluate the independent prognostic value of the clinical characteristics of OS. The prognostic ability of the predictive models for 1/3/5-year OS was evaluated by ROC curves (R package “timeROC”) and the area under the curve (AUC) values. In all analyses, all statistical P-values were bilateral, and P < 0.05 was considered statistically significant.

Clinical Specimens

Adrenocortical carcinoma tissues and normal adrenocortical tissues were obtained from five patients who received operation at the Third Xiangya Hospital (Changsha, Hunan, China) from January 2016 to December 2020. The patients were diagnosed by pathological analysis and were not subjected to chemotherapy or radiotherapy. The Institutional Review Board of the Ethics Committee of Third Xiangya Hospital approved the consent procedure, and written informed consent was provided by all patients in this research.

Immunohistochemistry

Paraffin-embedded ACC tissues and normal tissues were sliced, dewaxed, hydrated, and antigen-repaired, then endogenous peroxidase was blocked; anti-HNRNPA2B1 (1:100, 14813-1-AP, Proteintech Group, Wuhan, China), anti-LRPPRC (1:100, 21175-1-AP, Proteintech Group), and anti-ELAVL1 (1:100, 11910-1-AP, Proteintech Group) were added and incubated them together at 4°C overnight, respectively. Polymer enhancers were incubated for 30 min at room temperature, then biotinlabeled secondary antibodies were added and incubated for 30 min at room temperature. Next, the sections were stained by using diaminobenzidine staining solution, followed by counterstaining with hematoxylin, and then the sections were mounted in glycerol-vinyl-alcohol. Two independent professional pathologists were blinded to analysis the data and histopathological features of the patients, and also evaluated the IHC scores according to the scoring standards.

Results

Consensus Clustering of m6A Genes in Three Clusters With Different Clinical Outcomes of ACC

Here, the clinical data and corresponding gene expression profiles of ACC patients were downloaded from the TCGA and GEO databases. The workflow is shown in Figure 1. We first analyzed 21 m6A regulators and mapped the correlation between m6A patterns and ACC survival. As a result, the expression levels of m6A regulators and clinical characteristics were obtained from TCGA and GEO. Pearson correlation analysis was performed to determine the relationship between m6A regulators (Figure 2A). Some highly correlated (|correlation coefficient| ≥ 0.5, P < 0.05) m6A regulators were identified, such as YTHDF1, YTHDF2, HNRNPC, KIAA1429, ELAVL1, HNRNPA2B1, CBLL1, and YTHDF3. Based on the ConsensusClusterPlus R package, the TCGA-ACC cohort was clustered into different groups by consensus expression of m6A regulators. When the consensus matrix k value was equal to 3, there was the least crossover between the ACC samples (Figures 2B,C and Supplementary Figure 1). The OS difference between different clusters was calculated using the Kaplan-Meier method. Significantly better OS was found in patients specifically in cluster1 compared with other clusters (Figure 2D). Further, we plotted a boxplot (Figure 2E) and heatmap (Figure 4A) for visualizing the expression of the 21 m6A regulatory factors in clusters and found that the expression of CBLL1, ELAVL1, HNRNPA2B1, HNRNPC, KIAA1429, LRPPRC, RBM15, RBM15B, WTAP, YTHDF2, and YTHDF3 in cluster3 was higher than that in other clusters (P < 0.01), while the expression of ALKBH5, IGF2BP1, METTL3, and YTHDF1, was higher in cluster2 than in other clusters.

FIGURE 1

Figure 1. Study flow chart.

FIGURE 2

Figure 2. Consensus clustering of m6A genes. (A) The Pearson correlation analysis was used to search the relationship among m6A regulators. (B) Consensus clustering matrix for k = 3. (C) Consensus clustering cumulative distribution function (CDF) and relative change in area under CDF curve for k = 2 to 5. (D) Kaplan-Meier curves of OS for three clusters in ACC. (E) The expression of the 21 m6A regulatory factors in clusters (***P < 0.001; **P < 0.01; *P < 0.05).

The Interaction and Correlation Among the m6A Regulators in Three Patterns

To explore the potential biological differences among the three different m6A modification patterns, we identified DEGs by comparing clusters with the threshold of |logFC (fold change)| ≥ 1 and adj. P < 0.05. When comparing cluster1 and cluster2, there were 371 up-regulated and 292 down-regulated genes. GO analysis of biological processes showed that the DEGs were enriched in extracellular structure organization, matrix organization, and humoral immune response. Cellular component analysis indicated that DEGs were abundant in the extracellular matrix and collagen-containing extracellular matrix. Molecular function analysis indicated that DEGs were mainly located in receptor regulator and ligand activities. Furthermore, KEGG analysis showed that the DEGs were enriched in tyrosine metabolism, viral protein interaction with cytokine and cytokine receptor, and IL-17 signaling pathways (Figures 3A–C). Additionally, GSEA was performed for further signaling pathway enrichment analysis, and in a comparison between clusters1 and cluster2, the IL-17 signaling pathway, TNF, and NF-kappa B signaling pathway were enriched relative to cluster2 (Figure 3D). These signaling pathways are related to core biological carcinogenic processes, most of which are involved in the regulation of immune checkpoint expression, and have potential for further exploration of the effect of m6A modifications on immunotherapy (25, 35, 36).

FIGURE 3

Figure 3. The interaction and correlation among clusters. (A) The GO and KEGG analysis of DEGs between cluster1 and cluster2. (B) The GO and KEGG analysis of DEGs between cluster1 and cluster3. (C) The GO and KEGG analysis of DEGs between cluster2 and cluster3. (D) The GSEA analysis of DEGs between cluster1 and cluster2.

Immune Landscape in ACC Patients

The ESTIMATE algorithm provided stromal and immune scores, and tumor purity for all ACC samples. The heatmap of m6A-related gene expression and the stromal, immune, and ESTIMATE scores, and tumor purity are shown and clustered in Figure 4A. From the clustering, we found that such m6A genes showed a similar expression trend with tumor purity, and an opposite trend with stromal, immune, and ESTIMATE scores, indicating that the m6A pathway plays an important role in the tumor immune microenvironment and determines tumor progression and metastatic dissemination. Moreover, the density and location of immune cells can be quantified as a tangible indicator by an immune score. Here, there was a significant difference in the immune scores between the m6A clusters, and cluster1 showed the highest immune score. Furthermore, the ESTIMATE and stromal scores were also calculated, and the expression of cluster1 was higher than that of cluster2. Conversely, the distribution of tumor purity was different from the stromal, immune, and ESTIMATE scores, and cluster1 showed a lower tumor purity score than the others (Figure 4B).

FIGURE 4

Figure 4. Immune characteristics among three m6A patterns. (A) The heatmap of m6A regulators from 3 clusters and ESTIMATE algorithm. (B) Different expression of ESTIMATE score, immune score, stromal score, and tumor purity in three m6A Patterns. (C) Differences in the levels of infiltration of the 22 immune cells in three m6A patterns (***P < 0.001; **P < 0.01; *P < 0.05).

The 22 different immune cell types among different clusters were analyzed using the CIBERSORT algorithm. The results revealed that the macrophages M0, MI, M2 macrophages, dendritic cells activated, dendritic, cells resting, eosinophils mast cells resting, natural killer (NK) cells activated, and NK cells resting accounted for a large proportion of immune cell infiltration. Furthermore, cluster1, which had better survival, displayed a greater number of M1, M2 macrophages, and NK cells activated compared with other clusters with worse prognosis (Figure 4C). In addition, we found that the levels of dendritic cells, macrophages M0, and NK cells activated in cluster1 were significantly lower than those in cluster2 and cluster3. The outcome revealed that m6A-related patterns may remarkably suppress or strengthen the expression of specific immune cell types, thus potentially influencing the response to immunotherapy.

Prognostic Analysis of Risk Model and m6A Genes

To develop a signature for prognosis prediction of ACC, we performed Lasso Cox regression analysis on 21 m6A-related genes based on the TCGA database. Next, we obtained eight genes (METTL14, ZC3H13, FTO, YTHDF1, YTHDF3, HNRNPA2B1, LRPPRC, and ELAVL1) to build the risk model, and the coefficients of these genes were used to calculate the risk score (Figure 5A). The risk score = METTL14 × (−0.1750) + ZC3H13 × (−0.0212) + FTO × (−0.0984) + YTHDF1 × 0.0159 + YTHDF3 × (−0.0073) + HNRNPA2B1 × 0.0405 + LRPPRC × 0.0437 + ELAVL1 × 0.0376. Patients with ACC were separated into low-risk or high-risk groups with the median cutoff of risk score. As shown in Figures 5B,C, we found that as the risk score increased, high-risk patients had significantly worse OS than low-risk patients (P = 1.617e⁻⁰⁸). Univariate and multivariate analyses were used to evaluate the prognostic value of age, sex, M, N, T, and clinical stage (Supplementary Figure 2); however, only T stage was significantly correlated with OS. The ROC curve showed that the risk score had strong predictive ability, with an AUC of 0.844, 0.945, and 0.893 in 1, 3, and 5 years compared with other factors (Figure 5D). These results indicate that the risk model may serve as an important indicator for evaluating the prognosis of ACC.

FIGURE 5

Figure 5. Risk model from m6A-related genes. (A) Lasso Cox regression analysis of 21 m6A-related genes. (B) Overall survival analysis for patients in high/low risk. (C) The distributions of risk scores, alive/dead status, and expression of three m6A-related genes. (D) The ROC curve of risk score and clinical characteristics.

Validation of the m6A-Genes in the GEO Dataset

To verify and identify the key genes of m6A patterns, we evaluated the prognostic values of the eight genes and the risk model in the GEO datasets. Three datasets (GSE10927, GSE19750, and GSE33371) containing OS statistics and two datasets (GSE76019 and GSE76021) containing event-free survival (EFS) statistics were selected for further validation. Fortunately, the results indicated that the risk model have potential value of prediction in all GEO datasets (Supplementary Figure 3). Considering the similarities of identified genes in the TCGA and GEO data, it is believed that the overlapping m6A regulators might be significant, including HNRNPA2B1, LRPPRC, FTO, YTHDC1, and ELAVL1 (Supplementary Table 1). HNRNPA2B1 was validated as a significant indicator of poor OS and EFS based on all GEO datasets (Figure 6A). LRPPRC might be regarded as a crucial biomarker of poor survival, which was successfully validated on four GEO datasets (GSE10927, GSE33371, GSE76019, and GSE76021) (Figure 6B). Furthermore, ELAVL1 was re-verified as a potential biomarker based on two GEO containing EFS data (Figure 6C). However, other m6A-related genes received a failed verification on most GEO datasets (Supplementary Figures 3, 4). Subsequently, compared with the normal tissue, the expression of these three m6A regulators was upregulated in ACC by performing IHC (Supplementary Figure 5).

FIGURE 6

Figure 6. Identifying the key gene of m6A patterns. (A) Kaplan-Meier OS curves and EFS curves for patients in HNRNPA2B1. (B) Kaplan-Meier OS curves and EFS curves for patients in LRPPRC. (C) Kaplan-Meier OS curves and EFS curves for patients in ELAVL1.

SNP Analysis of m6A Genes Among Three Patterns

Tumor mutational burden (TMB) was considered a promising indicator and exhibited predictive utility in identifying responders and non-responders to immune checkpoint inhibitors (37, 38). Here, the most frequent variants were missense mutations, followed by nonsense mutations and splice sites. SNP was responsible for most variants, and single nucleotide variants (SNVs) mostly occurred as C > T and C > A. The top ten mutated genes in ACC samples were TTN, MUC16, PKHD1, TP53, CTNNB1, CNTNAP5, SVEP1, LRP1, HMCN1, and ASXL3 (Figure 7A). However, we failed to observe significant differences in these genes among the three clusters and other clinical characteristics (Figure 7B).

FIGURE 7

Figure 7. Landscape of TMB in ACC. (A) The summary plot of genetic alteration in ACC. (B) The oncoplot of genetic alteration in ACC.

To further analyze the relationship between m6A key genes and SNPs, we sorted out all the positive polymorphisms based on the SNP2APA database. Eight polymorphisms from four mRNAs (MID1IP1, CEBPZ, CNPY2, and TIPRL) were considered positive by the Kaplan-Meier survival method (Supplementary Table 2). The results of Pearson's correlation analysis indicated that these mRNAs were significantly correlated with the expression level of key m6A-related mRNAs, especially HNRNPA2B1, LRPPRC, and ELAVL1 (Supplementary Figure 6).

ceRNA Network Construction Based on the Key Gene

To reveal the function regulated by the m6A regulator pathway, we obtained 318 lncRNAs and 138 miRNAs in cluster2 compared with cluster1 (Supplementary Tables 3, 4). The heatmap of DEGs with |logFC (fold change)| ≥ 2 and FDR < 0.05 is shown in Figures 8A,B. To further understand how m6A-related lncRNAs mediate mRNA expression by sponging miRNAs, we constructed a ceRNA network based on m6A-related lncRNAs. Twenty-nine lncRNAs were selected from the miRcode database, which targeted 85 miRNAs and 19 miRNAs, were further identified after taking the intersection of DEGs of miRNAs. We then used StarBase to search for mRNAs, and a total of 455 mRNAs were selected based on three databases (miRTarBase, miRDB, and TargetScan) (Supplementary Table 5). Furthermore, taking the intersection of these mRNAs and 21 m6A-related genes as potential functional molecules, these data indicated that 12 lncRNAs (C17orf100, C8orf31, SNHG14, PLCL2-AS1, HCG11, HOTAIR, LINC00332, IDI2-AS1, ZFY-AS1, TTTY15, LINC00461, and DIO3OS) may dysregulate the behavior of hsa-mir-211 so that it promotes the expression of key m6A gene HNRNPA2B1 (Figures 8C,D).

FIGURE 8

Figure 8. ceRNA network construction of key m6A-related gene. (A) The heatmap of significant differential expression lncRNAs. (B) The heatmap of significant differential expression miRNAs. (C) The Venn plot of 455mRNAs and 21 m6A regulators. (D) The ceRNA network of the 12 m6A-related lncRNAs (blue) and target miRNAs (purple) and mRNAs (red).

Pan-Cancer Analysis of the Key Gene

To further confirm the key role of HNRNPA2B1 in the m6A modification process, we adopted a pan-cancer analysis and downloaded all the data from the UCSC Cancer Genomics Browser (https://genome-cancer.ucsc.edu), which offers interactive visualization and exploration of TCGA genomic and clinical data. HNRNPA2B1 expression was re-evaluated and was notably found to significantly impact prognosis in ACC. High HNRNPA2B1 expression was associated with decreased disease-specific survival (DSS) (P < 0.001), disease-free interval (DFI) (P = 0.034), and progression-free interval (PFI) (P < 0.001) (Figures 9A–C). In addition, as shown in Figure 9D, the expression of HNRNPA2B1 was significantly correlated with TNM stage. Furthermore, we analyzed the correlation among HNRNPA2B1 expression, ESTIMATE score, and infiltrating immune cells. The results showed that HNRNPA2B1 levels were significantly negatively correlated with immune score (r = −0.43, P = 0.00011), stromal score (r = −0.40, P = 0.00033), and resting mast cells (r = −0.56, P = 0.00012) (Figure 9E). Additionally, increased TMB has been linked to PD1/PD-L1 therapeutic response, and we found that HNRNPA2B1 mRNA levels were correlated to multiple types of cancers, including ACC, stomach adenocarcinoma (STAD), thymoma (THYM) etc. Microsatellite instability (MSI) is a pattern of hypermutation caused by defects in the mismatch repair system and has been known to be both predictive and prognostic to better profile responses to anti-PD-1 immunotherapy (39). Next, we found that HNRNPA2B1 was associated with diffuse large B-cell lymphoma (DBLC) and lung squamous cell carcinoma (LUSC), but not ACC, as shown in the radar plot (Figure 9F). Moreover, Spearman's correlation was used to show the relationship between HNRNPA2B1 and immune-related mRNAs in 33 cancer types (Figure 9G).

FIGURE 9

Figure 9. The pan-cancer analysis of HNRNPA2B1. (A) The relationship between HNRNPA2B1 expression and DSS. (B) The relationship between HNRNPA2B1 expression and DFI. (C) The relationship between HNRNPA2B1 expression and PFI. (D) The HNRNPA2B1 expression was significantly correlated with TNM stage. (E) The correlation among HNRNPA2B1 expression and immune landscape. (F) Pan-cancer analysis of the relationship between HNRNPA2B1 expression and TMB (above) or MSI (below). (G) The correlation among HNRNPA2B1 expression and immune-related mRNAs.

Discussion

With a heterogeneous clinical characteristic and poor OS, treating this complex malignant tumor, adrenocortical carcinoma, is a substantial clinical challenge (40). Further, the current TNM classification system remains inapplicable to predict prognosis. Here, we used univariate and multivariate analyses to evaluate the clinical prognostic value, and the results indicated that most clinical characteristics, including M stage and N stage, failed to guide treatment options. Previous studies have shown that the invasion and proliferation of ACC is regulated by genomic molecular characteristics (41). m6A is the most common and plentiful modification to affect cancer development through the regulation of m6A methyltransferases, demethylases, and binding proteins, but to date, the potential role of m6A regulators in ACC prognosis is not well-understood. Based on the TCGA dataset, we identified three clusters according to the optimal k-means clustering, and we observed a significant difference in OS among the three clusters, suggesting that the expression of m6A-related regulators is intimately related to the prognosis and malignancy of ACC.

Accumulating studies have focused on the tumor immune microenvironment and ACC is characterized as a highly immunogenic tumor, with 86.3% of ACC specimens showing high rates of tumor infiltrating lymphocytes (TILs) (42). Here, GO, KEGG, and GSEA analyses revealed that the DEGs among clusters were enriched in immune-related signaling pathways, such as IL-17 signaling pathway. Furthermore, the ESTIMATE score, especially the immune score, was significantly correlated with the expression of m6A patterns. Similarly, the infiltration of immune cells (M1, M2 macrophages, and NK cells activated) was significantly increased and the infiltration of immune cells (dendritic cells, macrophages M0, and NK cells activated) was significantly decreased in the low-risk score group (cluster1) compared with the high-risk score groups (cluster2 or cluster3). These results indicate a comprehensive evaluation of the m6A modification patterns that will facilitate understanding the characteristics of TME cell infiltration and promote individualized novel therapies by determining the response to immunotherapy.

Here, we systematically explored the effects of multiple m6A regulators on OS in ACC and attempted to construct a risk model for prediction. Next, we distinguished a prognostic risk signature with eight identified m6A regulators (METTL14, ZC3H13, FTO, YTHDF1, YTHDF3, HNRNPA2B1, LRPPRC, and ELAVL1), which divided the OS in ACC into high-risk subgroups with high mortality and low-risk subgroup with remarkably better survival. Notably, compared with the previous prognostic markers (T, N, M clinical stage), our prognostic risk signature can achieve higher accuracy, with AUC values >0.8. In summary, the risk signature we constructed might be viewed as a new potential and promising biomarker that can provide more precise clinical applications and an efficient guide for treatment. However, owing to the limited number of samples in most GEO databases, this model needs to be re-confirmed in other databases with large populations. Meanwhile, we applied survival analysis to further detect the eight identified genes in the TCGA and GEO datasets, and we served overlapping m6A regulators as significant biomarkers, particularly HNRNPA2B1, which was proven to have a positive association with poor OS, EFS, DSS, DFI, and PFI in TCGA and all GEO datasets.

The role of HNRNPA2B1 in cancer has recently garnered increasing attention. On the one hand, HNRNPA2B1 acted as a nuclear m6A reader that recruited the miRNA microprocessor complex protein DGCR8 to a subset of precursor miRNAs and mediated the mature miRNA processing (43). On the other hand, HNRNPA2B1 functioned as an adaptor and modulated the molecular changes to alternative splicing combined with METTL3 (44). Here, the results showed that lower immune scores, stromal scores, and ESTIMATE scores were significantly associated with higher HNRNPA2B1 expression, which demonstrated that this m6A gene signature played a non-negligible role in shaping diverse stromal and immune TME landscapes, implying that HNRNPA2B1 may affect the therapeutic efficacy of immune checkpoint blockade. Moreover, using the CIBERSORT algorithm, we found that HNRNPA2B1 mediated the TME infiltration patterns to accelerate ACC progression partly by regulating macrophages M0. Tumor mutational burden, defined as the total number of somatic coding mutations per million bases, has emerged as a notable biomarker of response to immunotherapy (45). Recent studies showed that HNRNPA2B1 with high frequency of mutation may have an influence on promoting tumorigenesis in melanoma (46). Similarly, by analyzing the mutation annotation files of the TCGA cohort, the results showed that HNRNPA2B1 has a close correlation with TMB in multiple types of cancer, particularly ACC. Based on the SNP2APA database, we identified four highly variant mutated genes, most of which were highly correlated with the HNRNPA2B1 expression level. These findings indicated that high HNRNPA2B1 expression is related to dysregulation of the TME landscape and a sharp accumulation of gene mutations, thus becoming a promising therapeutic target for ACC.

Furthermore, several selected lncRNAs were reported to be associated with cancer progression (47), but there have been few reports on lncRNAs regarding ACC progression, and how m6A-related genes act in an lncRNA-miRNA-dependent manner during ACC progression is still unknown. In our research, we constructed a ceRNA network to target m6A-related miRNA and lncRNAs, which consist of twelve lncRNAs and one miRNA. Therefore, we should pay more attention to filter key lncRNAs that could predict OS and EFS and have a close relationship with HNRNPA2B1. Finally, we identified two lncRNAs, HOTAIR and IDI2-AS1, according to the above-mentioned conditions. Located on chromosome.12q13.13, lncRNA HOTAIR (HOX Transcript Antisense Intergenic RNA) was regarded as a regulator of chromatin states (48). In immune cells, HOTAIR has the ability to interfere with the TME landscape by inducing IκBα degradation, with the consequent activation of NF-κB pathways and secretion of pro-inflammatory cytokines (49, 50). Moreover, HOTAIR also led to the downregulation of the tumor suppressor gene SETD2, promoting MSI and high TMB (50, 51). In ACC, HOTAIR is overexpressed in tumor tissues compared with normal tissues. For in vitro experiments, HOTAIR can prompt the progression of ACC by shortening the cell cycle and promoting the proliferation of ACC cells (52). Hence, we planned to take the next step to deeply explore the interaction of HOTAIR and HNRNPA2B1 in vitro. In our study, IDI2-AS1 was also found to play a significant role in ACC. However, there are few reports on how IDI2-AS1 acts as an oncogene. Thus, we carried out a pan-cancer analysis of IDI2-AS1 and hoped our results help to identify the prognostic value. We found that high IDI2-AS1 expression was obviously associated with poor DSS, DFI, PFI, and high TMB rate. Furthermore, the findings revealed that IDI2-AS1 had a stable association with the expression of PD1 (PDCD1) and PD-L1 (CD274) to serve as a potential prognostic marker or therapeutic targets of cancers integrated with m6A-related regulators (Supplementary Figure 7). Moreover, previous studies have reported the different functions of miR-211, which can facilitate platinum chemosensitivity by blocking the DNA damage response (53) or by stimulating the emergence of BRAF inhibitor resistance (54). Most speculations that miR-211 may regulate the m6A modification pathway need to be confirmed by further validation or experiments.

In conclusion, this study is the first to comprehensively identify and systematically profile the gene signatures of m6A-related regulators in ACC. The different m6A modification patterns played an important role in the heterogeneity and complexity of the TME. We also developed an eight-gene-signature prognostic model, in particular HNRNPA2B1, which might determine the clinical progression of ACC. Moreover, we constructed a ceRNA network to further decipher the molecular mechanisms based on HNRNPA2B1. Our results indicate that m6A genes are promising predictive indicators that may provide novel insights into ACC therapeutic strategies and guide effective immunotherapy.

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.

Ethics Statement

Written informed consent was obtained from the individual(s), and minor(s)' legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author Contributions

KC: designed the study. YJ: analyzed, interpreted the data, and wrote original draft. ZW, DH, YZ, and XH: wrote this manuscript. LG, MX, XC, and YC: edited and revised the manuscript. All authors have seen and approved the final version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (81874137), the Outstanding Youth Foundation of Hunan Province (2018JJ1047), the Hunan Province Science and Technology Talent Promotion Project (2019TJ-Q10), Hunan Cancer Hospital Climb Plan (No. QH201905), and the Special Foundation for Cancer Research from the National Cancer Center (NCC2017L01, and NCC2017A17), China.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.637933/full#supplementary-material

Supplementary Figure 1. Consensus clustering matrix for k = 2 to 5.

Supplementary Figure 2. (A) The univariate analysis of age, sex, M, N, T, and clinic stage. (B) The multivariate analysis of age, sex, M, N, T, and clinic stage.

Supplementary Figure 3. Prognostic value of risk model, YTHDF1 and YTHDF3. (A) Kaplan-Meier OS curves and EFS curves for patients in risk model. (B) Kaplan-Meier OS curves and EFS curves for patients in YTHDF1. (C) Kaplan-Meier OS curves and EFS curves for patients in YTHDF3.

Supplementary Figure 4. Prognostic value of FTO, METTL14 and ZC3H13. (A) Kaplan-Meier OS curves and EFS curves for patients in FTO. (B) Kaplan-Meier OS curves and EFS curves for patients in METTL14. (C) Kaplan-Meier OS curves and EFS curves for patients in ZC3H13.

Supplementary Figure 5. Verification of key m6A regulators expression in immunohistochemistry. (A) The expression of HNRNPA2B1 in tumor and normal tissue. (B) The expression of LRPPRC in tumor and normal tissue. (C) The expression of ELAVL1 in tumor and normal tissue. (***P < 0.001; **P < 0.01).

Supplementary Figure 6. Pearson's correlation analysis of positive polymorphisms and key m6A regulators.

Supplementary Figure 7. The pan-cancer analysis of IDI2-AS1. (A) The relationship between IDI2-AS1 expression and OS. (B) The relationship between IDI2-AS1 expression and DSS. (C) The relationship between IDI2-AS1 expression and DFI. (D) The relationship between IDI2-AS1 expression and PFI. (E) The correlation among IDI2-AS1 expression and immune landscape. (F) Pan-cancer analysis of the relationship between IDI2-AS1 expression and TMB. (G) The correlation among IDI2-AS1 expression and PD1/PD-L1 from starBase V3.0 project.

Supplementary Table 1. Cox analysis of 21 m6A genes in ACC patients.

Supplementary Table 2. Positive polymorphisms based on SNP2APA database.

Supplementary Table 3. All differential expression lncRNAs in cluster2 compared with cluster1.

Supplementary Table 4. All differential expression miRNAs in cluster2 compared with cluster1.

Supplementary Table 5. Targeted mRNAs base on miRTarBase, miRDB, and TargetScan databases.

Abbreviations

ACC, adrenocortical carcinoma; TNM, tumor, lymph node, and metastasis; m6A, N6-methyladenosine; TME, the tumor microenvironment; BMDCs, bone marrow-derived cells; TCGA, the cancer genome atlas; GEO, the gene expression omnibus database; GDC, genomic data commons; OS, the overall survival; DEG, differentially expressed gene; GSEA, gene set enrichment analysis; LASSO, the least absolute shrinkage and selection operator; HR, the hazard ratio; ROC, the receiver operating characteristic; SNP, single nucleotide polymorphism; AUC, the area under the curve; EFS, the event-free survival; TMB, the tumor burden; SNV, single nucleotide variants; DDS, the disease-specific survival; DFI, the disease-free interval; PFI, the progression-free interval; TILs, tumor infiltrating lymphocytes; HOTAIR, HOX transcript antisense intergenic RNA.

References

1. Kebebew E, Reiff E, Duh QY, Clark OH, McMillan A. Extent of disease at presentation and outcome for adrenocortical carcinoma: have we made progress. World J Surg. (2006) 30:872–8. doi: 10.1007/s00268-005-0329-x

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Allolio B, Fassnacht M. Clinical review: adrenocortical carcinoma: clinical update. J Clin Endocrinol Metab. (2006) 91:2027–37. doi: 10.1210/jc.2005-2639

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Gonzalez RJ, Tamm EP, Ng C, Phan AT, Vassilopoulou-Sellin R, Perrier ND, et al. Response to mitotane predicts outcome in patients with recurrent adrenal cortical carcinoma. Surgery. (2007) 142:867–75. doi: 10.1016/j.surg.2007.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Icard P, Chapuis Y, Andreassian B, Bernard A, Proye C. Adrenocortical carcinoma in surgically treated patients: a retrospective study on 156 cases by the French Association of Endocrine Surgery. Surgery. (1992) 112:972–9.

PubMed Abstract | Google Scholar

5. Else T, Kim AC, Sabolch A, Raymond VM, Kandathil A, Caoili EM, et al. Adrenocortical carcinoma. Endocr Rev. (2014) 35:282–326. doi: 10.1210/er.2013-1029

CrossRef Full Text | Google Scholar

6. Assié G, Jouinot A, Fassnacht M, Libé R, Garinet S, Jacob L, et al. Value of molecular classification for prognostic assessment of adrenocortical carcinoma. JAMA Oncol. (2019) 5. doi: 10.1001/jamaoncol.2019.1558

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Fassnacht M, Johanssen S, Quinkler M, Bucsky P, Willenberg HS, Beuschlein F, et al. European Network for the Study of Adrenal, Limited prognostic value of the 2004 International Union Against Cancer staging classification for adrenocortical carcinoma: proposal for a Revised TNM Classification. Cancer. (2009) 115:243–50. doi: 10.1002/cncr.24030

CrossRef Full Text | Google Scholar

8. Alarcón CR, Lee H, Goodarzi H, Halberg N, Tavazoie SF, Alarcón CR, et al. N6-methyladenosine marks primary microRNAs for processing. Nature. (2015) 519:482–5. doi: 10.1038/nature14281

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. (2017) 18:31–42. doi: 10.1038/nrm.2016.132

CrossRef Full Text | Google Scholar

10. Lan Q, Liu PY, Haase J, Bell JL, Hüttelmaier S, Liu T. The critical role of RNA m6A methylation in cancer. Cancer Res. (2019) 79:1285–92. doi: 10.1158/0008-5472.CAN-18-2965

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chen XY, Zhang J, Zhu JS. The role of m6A RNA methylation in human cancer. Mol Cancer. (2019) 18:103. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Meyer KD, Jaffrey SR. Rethinking m6A readers, writers, and erasers. Annu Rev Cell Dev Biol. (2017) 33:319–42. doi: 10.1146/annurev-cellbio-100616-060758

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M. m6A mRNA methylation facilitates resolution of naive pluripotency toward differentiation. Science. (2015) 347:1002–6. doi: 10.1126/science.1261417

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Fustin JM, Doi M, Yamaguchi Y, Hida H, Nishimura S, Yoshida M. RNA-methylation-dependent RNA processing controls the speed of the circadian clock. Cell. (2013) 55:793–806. doi: 10.1016/j.cell.2013.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xiang Y, Laurent B, Hsu CH, Nachtergaele S, Lu Z, Sheng W, et al. RNA m(6)A methylation regulates the ultraviolet-induced DNA damage response. Nature. (2017) 543:573–6. doi: 10.1038/nature21671

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Tong J, Cao G, Zhang T, Sefik E, Amezcua Vesely MC, Broughton JP, et al. m(6)A mRNA methylation sustains Treg suppressive functions. Cell Res. (2018) 28:253–6. doi: 10.1038/cr.2018.7

CrossRef Full Text | Google Scholar

17. Pinello N, Sun S, Wong JJ. Aberrant expression of enzymes regulating m6A mRNA methylation: implication in cancer. Cancer Biol Med. (2018) 15:323–34. doi: 10.20892/j.issn.2095-3941.2018.0365

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hong K. Emerging function of N6-methyladenosine in cancer. Oncol Lett. (2018) 16:5519–24. doi: 10.3892/ol.2018.9395

CrossRef Full Text | Google Scholar

19. Niu Y, Lin Z, Wan A, Chen H, Liang H, Sun L, et al. RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3. Mol Cancer. (2019) 18:46. doi: 10.1186/s12943-019-1004-4

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tu Z, Wu L, Wang P, Hu Q, Tao C, Li K, et al. N6-Methylandenosine-related lncRNAs are potential biomarkers for predicting the overall survival of lower-grade glioma patients. Front Cell Dev Biol. (2020) 8:642. doi: 10.3389/fcell.2020.00642

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, et al. Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell. (2017) 32:42–56.e6. doi: 10.1016/j.ccell.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. (2012) 366:2443–54. doi: 10.1056/NEJMoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Raj N, Zheng Y, Kelly V, Katz SS, Chou J, Do RKG, et al. PD-1 blockade in advanced adrenocortical carcinoma. J Clin Oncol. (2020) 38:71–80. doi: 10.1200/JCO.19.01586

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li N, Kang Y, Wang L, Huff S, Tang R, Hui H, et al. ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment. Proc Natl Acad Sci USA. (2020) 117:20159–70. doi: 10.1073/pnas.1918986117

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yang S, Wei J, Cui YH, Park G, Shah P, Deng Y, et al. m(6)A mRNA demethylase FTO regulates melanoma tumorigenicity and response to anti-PD-1 blockade. Nat Commun. (2019) 10:2782. doi: 10.1038/s41467-019-10669-0

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Shi H, Wei J, He C. Where, when, and how: context-dependent functions of RNA methylation writers, readers, and erasers. Mol Cell. (2019) 74:640–50. doi: 10.1016/j.molcel.2019.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. (2018) 28:616–24. doi: 10.1038/s41422-018-0040-8

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Method. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yang Y, Zhang Q, Miao YR, Yang J, Yang W, Yu F, et al. SNP2APA: a database for evaluating effects of genetic variants on alternative polyadenylation in human cancers. Nucleic Acids Res. (2020) 48:D226–32. doi: 10.1093/nar/gkz793

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Jeggari A, Marks DS, Larsson E. miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics. (2012) 28:2062–3. doi: 10.1093/bioinformatics/bts344

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic acids Res. (2014) 42:D92–7. doi: 10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kim HD, Song GW, Park S, Jung MK, Kim MH, Kang HJ, et al. Association between expression level of PD1 by tumor-infiltrating CD8(+) T cells and features of hepatocellular carcinoma. Gastroenterology. (2018) 155:1936–50.e17. doi: 10.1053/j.gastro.2018.08.030

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wang X, Yang L, Huang F, Zhang Q, Liu S, Ma L, et al. Inflammatory cytokines IL-17 and TNF-alpha up-regulate PD-L1 expression in human prostate and colon cancer cells. Immunol Lett. (2017) 184:7–14. doi: 10.1016/j.imlet.2017.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Cristescu R, Mogg R, Ayers M, Albright A, Murphy E, Yearley J, et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. Science. (2018) 362:eaar3593. doi: 10.1126/science.aar3593

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. (2017) 9:34. doi: 10.1186/s13073-017-0424-2

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Mandal R, Samstein RM, Lee KW, Havel JJ, Wang H, Krishna C, et al. Genetic diversity of tumors with mismatch repair deficiency influences anti-PD-1 immunotherapy response. Science. (2019) 364:485–91. doi: 10.1126/science.aau0447

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Assié G, Letouzé E, Fassnacht M, Jouinot A, Luscap W, Barreau O, et al. Integrated genomic characterization of adrenocortical carcinoma. Nat Genet. (2014) 46:607–12. doi: 10.1038/ng.2953

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wasserman JD, Novokmet A, Eichler-Jonsson C, Ribeiro RC, Rodriguez-Galindo C, Zambetti GP, et al. Prevalence and functional consequence of TP53 mutations in pediatric adrenocortical carcinoma: a children's oncology group study. J Clin Oncol. (2015) 33:602–9. doi: 10.1200/JCO.2013.52.6863

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Landwehr LS, Altieri B, Schreiner J, Sbiera I, Weigand I, Kroiss M, et al. Interplay between glucocorticoids and tumor-infiltrating lymphocytes on the prognosis of adrenocortical carcinoma. J Immunother Cancer. (2020) 8:e000469. doi: 10.1136/jitc-2019-000469

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Alarcón CR, Goodarzi H, Lee H, Liu X, Tavazoie S, Tavazoie SF. HNRNPA2B1 is a mediator of m(6)a-dependent nuclear RNA processing events. Cell. (2015) 162:1299–308. doi: 10.1016/j.cell.2015.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Spitale RC, Flynn RA, Zhang QC, Crisalli P, Lee B, Jung JW, et al. Structural imprints in vivo decode RNA regulatory mechanisms. Nature. (2015) 519:486–90. doi: 10.1038/nature14263

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Allgäuer M, Budczies J, Christopoulos P, Endris V, Lier A, Rempel E, et al. Implementing tumor mutational burden (TMB) analysis in routine diagnostics-a primer for molecular pathologists and clinicians. Transl Lung Cancer Res. (2018) 7:703–15. doi: 10.21037/tlcr.2018.08.14

CrossRef Full Text | Google Scholar

46. Li T, Gu M, Deng A, Qian C. Increased expression of YTHDF1 and HNRNPA2B1 as potent biomarkers for melanoma: a systematic analysis. Cancer Cell Int. (2020) 20:239. doi: 10.1186/s12935-020-01309-5

CrossRef Full Text | Google Scholar

47. Loewen G, Jayawickramarajah J, Zhuo Y, Shan B. Functions of lncRNA HOTAIR in lung cancer. J Hematol Oncol. (2014) 7:90. doi: 10.1186/s13045-014-0090-4

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Rinn JL, Kertesz M, Wang JK, Squazzo SL, Xu X, Brugmann SA, et al. Functional demarcation of active and silent chromatin domains in human HOX loci by noncoding RNAs. Cell. (2007) 29:1311–23. doi: 10.1016/j.cell.2007.05.022

PubMed Abstract | CrossRef Full Text

49. Xia Y, Shen S, Verma IM. NF-kappaB, an active player in human cancers. Cancer Immunol Res. (2014) 2:823–30. doi: 10.1158/2326-6066.CIR-14-0112

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Botti G, Scognamiglio G, Aquino G, Liguori G, Cantile M. LncRNA HOTAIR in tumor microenvironment: what role? Int J Mol Sci. (2019) 20:2279. doi: 10.3390/ijms20092279

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Li F, Mao G, Tong D, Huang J, Gu L, Yang W, et al. The histone mark H3K36me3 regulates human DNA mismatch repair through its interaction with MutSalpha. Cell. (2013) 153:590–600. doi: 10.1016/j.cell.2013.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yan ZC, He L, Qiu JH, Deng WT, Lu JR, Yuan Z, et al. LncRNA HOTAIR participates in the development and progression of adrenocortical carcinoma via regulating cell cycle. Eur Rev Med Pharmacol Sci. (2018) 22:6640–49. doi: 10.26355/eurrev_201810_16139

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wang T, Hao D, Yang S, Ma J, Yang W, Zhu Y, et al. miR-211 facilitates platinum chemosensitivity by blocking the DNA damage response (DDR) in ovarian cancer. Cell Death Dis. (2019) 10:495. doi: 10.1038/s41419-019-1715-x

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Díaz-Martínez M, Benito-Jardón L, Alonso L, Koetz-Ploch L, Hernando E, Teixidó J. miR-204-5p and miR-211-5p contribute to BRAF inhibitor resistance in melanoma. Cancer Res. (2018) 78:1017–30. doi: 10.1158/0008-5472.CAN-17-1318

PubMed Abstract | CrossRef Full Text | Google Scholar