Competing Endogenous RNA Network Analysis Reveals Pivotal ceRNAs in Adrenocortical Carcinoma

Objective: To construct ceRNA network and identify pivotal competing endogenous RNAs (ceRNAs) in adrenocortical carcinoma (ACC) using ceRNA network analysis. Methods: The RNA sequencing expression data of 77 ACCs in TCGA were obtained from GEPIA. Cancer specific ceRNAs, cancer specific microRNAs (miRNAs), and cancer specific messenger RNAs (mRNAs) were identified. The interaction of cancer specific miRNAs with cancer specific ceRNAs and cancer specific mRNAs were predicted. CeRNA network was constructed and visualized by Cytoscape 3.7.0 software. The genes in ceRNA network regulated GO terms and regulated pathways were performed by function analysis. Survival analysis of pivotal ceRNAs was performed for the pivotal lncRNAs. Result: Twenty-eight cancer specific ceRNAs, 149 cancer specific miRNAs, and 104 mRNAs were identified. CeRNA network was constructed including 10 ceRNAs, 35 miRNAs, and 34 mRNAs. The genes in ceRNA network regulated GO terms and were classified into three groups: cellular component (CC), molecular function (MF), and biological process (BP). The genes in ceRNA network regulated the following pathways: leukocyte transendothelial migration, and proteoglycans in cancer. Survival analysis showed that CTB-63M22.1 and RP1-241P17.4 were significantly associated with ACC patient disease free survival and overall survival. Conclusion: This study has constructed ceRNA networks in ACC. The study provides a set of pivotal ceRNAs for future investigation into the molecular mechanisms.


INTRODUCTION
Cancers are often associated with aberrant transcriptomes. Recent analyses of the human transcriptome have revealed that 2% of the genome encode protein-coding transcripts even though over three-quarters of the genome are transcribed (1,2). As genes do not function in isolation, they can be grouped into "networks" based on their interactions. A novel hypothesis, known as the ceRNA hypothesis presented by Salmena et al. (3), has emerged in recent years. CeRNAs, acting as miRNA sponges, can bind miRNAs competitively to indirectly regulate mRNAs using shared microRNA response elements. In theory, ceRNAs refer to all transcripts that may become the targets of miRNA such as long non-coding RNA (lncRNA), pseudogene RNA, and circular RNA. The function of miRNAs is relatively well-understood. MiRNAs can interact with target mRNAs to regulate the expression of the gene (4). Studies have shown that miRNAs are involved in the initiation and progression of cancers (4). While studies on lncRNAs and transcribed pseudogenes are still in their infancy, some studies showed that lncRNAs and pseudogenes played key regulatory roles in cancer biology (5). Recently, investigations about the ceRNA networks provide a better understanding of tumorigenesis (3,6). This hypothesis laid the foundations for an important discovery regarding cross talk between both coding and non-coding RNAs through miRNAs.
Though an increasing number of miRNAs (7,8) and lncRNAs (9) have been discovered to play a role in endocrine cancer. Our knowledge about the ceRNA network in endocrine cancer is still limited. The regulation landscape of ceRNA crosstalk in cancers is still largely unknown.
In this study, we focused on adrenocortical carcinoma (ACC). ACC is a rare endocrine malignancy, often with an unfavorable prognosis (10, 11); 5-years overall survival is below 40% (11). According to a report based on the Surveillance, Epidemiology, and End Results (SEER) database, the incidence of ACC is ∼0.72 per million cases per year leading to 0.2% of all cancer deaths in the United States (12). The recent identification of the genomic alterations had provided some molecular mechanisms (13), but ceRNA crosstalk in ACC is still unknown. In this study, we have taken full advantage of the rich data from the TCGA consortium. The aberrant expression profiles of RNA in ACC were filtered. The aberrant ceRNA network was constructed. This paper represented a significant leap forward in understanding the biological functions of ceRNAs in endocrine cancers.

Data Sources
The RNA sequencing expression data of 77 ACC patients in TCGA were obtained from GEPIA, a web server for cancer, normal gene expression profiling, and interactive analyses (http://gepia.cancer-pku.cn/index.html) (14). The dataset was strict to ACC, with [log 2 FC] > 3 and p-value < 0.01. ANOVA was used as differential methods. Both over-expressed and underexpressed expression data were downloaded.

Construction of ceRNA Network of ACC
The RNA expression data downloaded from GEPIA were annotated by R package RefSeq mRNA and RefSeq ncRNA. We then used R package transcript biotype to double check.
First, we identified cancer specific ceRNA, cancer specific miRNA, and cancer specific mRNA. The traditional view was that lncRNA functioned as a ceRNA. Recent articles suggested that pseudogenes also have ceRNA function (15). Therefore, in this paper lncRNA and pseudogenes were both identified as potential  Frontiers in Endocrinology | www.frontiersin.org cancer specific ceRNAs. Cancer specific miRNAs were retained from Oncomir (http://www.oncomir.org) (16). A list of miRNAs closely associated with patient survival or tumor development were retrieved with p < 0.05. mRNAs with absolute [log 2 FC] > 3 and p-value < 0.01 retained from GEPIA were identified as cancer specific mRNAs.

Functional Analysis
The biological function and pathway of the genes involved in the ceRNA network would demonstrate instructive information. Thus, we used Gene ontology analysis (GO) to identify characteristic biological attributes. Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) enrichment analysis was performed to identify functional attributes. p < 0.05 was set as the cut-off criterion.

Survival Analysis of Pivotal lncRNAs
The ceRNAs in ACC ceRNA network were identified as pivotal ceRNAs. CeRNAs' correlations with patient survival were featured in GEPIA. Available TCGA patient survival data were used for Kaplan-Meier survival analysis and to generate overall and disease-free survival plots. Group-cutoff was set as the median, the hazards ratio was calculated based on Cox PH Model and the 95% confidence interval was added as a dotted line. p < 0.05 was set as the cut-off criterion.

Identification of Cancer Specific ceRNAs, miRNAs, and mRNAs
Seventy-seven ACC patient data in TCGA were analyzed. Twenty-eight cancer specific ceRNAs (10 down regulated in ACC and 18 up regulated in ACC, Supplementary Table 1), 149 cancer specific miRNAs and 104 cancer specific mRNAs were identified (Supplementary Tables 2,3).

CeRNA Network Construction
To construct the ceRNA network, we next predicted the interaction of cancer specific miRNAs with cancer specific ceRNAs and cancer specific mRNAs. We overlapped the predicted targets of cancer specific ceRNAs with cancer specific miRNAs and overlapped the predicted targets of cancer specific miRNAs with cancer specific mRNAs. Finally, the ceRNA regulatory network was constructed (Figure 1A), including 10 ceRNAs, 35 miRNAs, and 34 mRNAs. In the down-regulated ceRNA regulatory network, there were 4 ceRNAs, 31 miRNAs, and 33 mRNAs involved ( Figure 1B). And in the up-regulated ceRNA regulatory network, there were six ceRNAs, 24 miRNAs, and 25 mRNAs involved ( Figure 1C).

Gene Ontology and Pathway Enrichment Analyses
We performed gene ontology analysis of the genes in the ceRNA network (Figures 2-4, Tables 1-3). The result showed that related to CC, the genes in the ceRNA network were mainly enriched in extracellular region. Related to MF, the genes in the ceRNA network were mainly enriched in fibronectin binding. Related to BP, the genes in the ceRNA network were mainly enriched in inflammatory response and immune response. We then performed pathway enrichment analysis of the genes in the ceRNA network and identified the ceRNAs regulated pathway. Functional analysis showed that the ceRNAs regulated pathways included leukocyte transendothelial migration, and proteoglycans in cancer.

Survival Analysis of Pivotal lncRNAs
We used GEPIA to perform survival analysis of pivotal ceRNAs (Figures 5, 6). p < 0.05 was set as the cut-off criterion. The results showed that CTB-63M22.1 and RP1-241P17.4 were significantly associated with ACC patient disease free survival and overall survival.

Salmena et al. (3) firstly proposed a hypothesis on how mRNAs
and lncRNAs "talk" to each other using miRNA as letters of a new language. The last few years have nominated ceRNA hypothesis as a novel layer of gene regulation. This ceRNA regulatory network greatly expands the functional genetic information and will aid in deciphering the underlying pathogenesis of human polygenic disease at the post-transcriptional level. Previously, Xiao et al. (20) have analyzed the gene expression data of ACC in The Cancer Genome Atlas (TCGA) and identified the potential prognostic mRNA biomarkers. To our knowledge, this is the first study to investigate the ceRNA network in ACC.
In this study, we extracted RNA expression profile from cases of ACC from TCGA. Twenty-eight cancer specific ceRNAs, 149 cancer specific miRNAs, and 104 mRNAs were identified. A ceRNA network was then constructed, consisting of 10 ceRNAs, 35 miRNAs, and 34 miRNAs. The Kyoto Encyclopedia of Genes and Genomes (KEGG), and the Gene Ontology (GO) databases were used to identify biological and functional attributes of the genes in the ceRNA network. The ceRNA network-regulated GO    terms were classified into three groups. The ceRNA network was mainly enriched in the following pathways: leukocyte transendothelial migration and proteoglycans in cancer. Survival analysis showed two potential prognostic associated ceRNAs: CTB-63M22.1 and RP1-241P17.4.
Of the genes that were involved in the ceRNA network, GO functional enrichment analysis showed that these genes were mainly enriched in the extracellular region, fibronectin binding, inflammatory response, and immune response. Perge's (21) previous study showed that miRNAs in exosome (an extracellular substance) might have potential diagnostic value in ACC. MiR-31 in serum exosomes possibly correlated with an aggressive profile in endocrine cancer metastasis (8). Interestingly, miR-31 was a miRNA identified in our ceRNA network. Griggs et al. (22) reported that fibronectin, a protein in the extracellular matrix, served as a growth factor delivery system in epithelialmesenchymal transition of cancers. Host immune response status played an important role in the clinical outcome of patients with cancer. Systemic inflammatory response was related with cancer progression (23,24).
Of the genes that were involved in ceRNA network, KEGG analysis showed that these genes were mainly enriched in leukocyte transendothelial migration and proteoglycans in cancer. Leukocyte transendothelial migration is generally activated in cancer progression (25). Proteoglycans (PGs), comprising structurally diverse constituents of the extracellular matrix and cell surfaces, have emerged as novel biomarkers and molecular players both within tumor cells and within their microenvironment. Recent studies have demonstrated that the expression of proteoglycans is regulated by miRNAs (26,27). The ceRNA network identified in our study provided useful clues of further study.
By analyzing TCGA data, 35 cancer specific miRNA were identified in the ceRNA network. Several miRNAs in our list have been investigated in ACC: Doghman et al. (28) identified that miR-99a was differentially expressed in childhood adrenocortical tumors. Mir-99a might involve in tumor cell growth via regulating mTOR signaling. Bimpak et al. (29) found miR-200 was differentially expressed in massive macronodular adrenocortical disease and miR-200b directly targeted matrin 3 expression in an adrenocortical cancer cell line. The understanding of other miRNAs identified in this study is limited, therefore, further molecular studies are require. These ACC-specific miRNAs might become potential biomarkers with specificity in the diagnosis and progression of ACC in future.
In our study, we identified 10 pivotal ceRNAs. Liu's study showed that H19 expression is maintained at low levels in hormonally active adrenocortical carcinomas, which is consistent with our study (30). There is some evidence indicating that the ceRNAs identified in our study are associated with other cancer types. Xiong et al. (31) reported that NEAT1, a ceRNA, accelerated lung adenocarcinoma deterioration. Teng et al.'s (32) recent study showed that lncRNA ARAP1-AS1, which functioned as a ceRNA, promoted the progression of bladder cancer by regulating the miR-4735-3p/NOTCH2 axis. Studies of other ceRNAs are limited, but they provide a good direction for future research. Of the 10 pivotal ceRNAs, CTB-63M22.1, and RP1-241P17.4 were associated with tumor prognosis. There are currently no studies on these two ceRNAs, so they are worthy of future research.
Although the findings of our study have important clinical implications, the limitations must also be noted. In addition to lncRNA and pseudogenes, cirRNA and tRNA may also function as ceRN, and their role in ceRNA network needs further research. Further studies with much larger sample sizes and longer follow-up are required to verify our results. The biological roles of ceRNAs in ACC need to be further investigated. In the future, molecular biology methods including qPCR, luciferase reporter systems and co-immunoprecipitation assays are helpful to validate our findings, thus unraveling the molecular mechanisms of ceRNA networks. CeRNA activity is associated with a multitude of parameters such as the number of miRNA binding types residing on coding and non-coding transcripts, and these parameters need further investigation.
In conclusion, this study has constructed ceRNA network in ACC. The study provides a set of pivotal ceRNAs for future investigation into the molecular mechanisms. This study provided novel insights to explore the underlying mechanism of endocrine cancers.

AUTHOR CONTRIBUTIONS
WL drew the figures and wrote the paper. FS analyzed the data.