lncRNA-miRNA-mRNA ceRNA Network in Thyroid Carcinoma from The Cancer Genome Atlas

Increasing evidence indicates that the competitive endogenous RNA (ceRNA) hypothesis, that is, long non-coding RNA (LncRNA) can competitively bind microRNA (miRNA) through miRNA response elements to affect the expression of target RNA, and dysregulation of LncRNA expression plays a key role in tumor progression. The papillary thyroid carcinoma that we studied is the most signicant pathological type of thyroid cancer, but its ceRNA network has not been extensively evaluated. We analyzed level-3 data from RNA-Seq of 58 para-carcinoma tissues and 501 patients with primary papillary thyroid carcinoma (PTC) using the DEseq software package and downloaded clinical information from The Cancer Genome Atlas (TCGA) to nd potential biomarkers or therapeutic targets. As a result, 149 differential miRNAs were selected, including 117 up -regulated, 32 down-regulated, and 3099 differential mRNAs, including 1976 up-regulated, 1123 down-regulated, and 434 differential lncRNAs, including 331 up-regulated and 103 down-regulated (Fold Change > 2, P < 0.05). The interactions between these differentially expressed RNA groups constitute the ceRNA network of PTC. Moreover, we used the microde database to predict the miRNAs that may be acted by the above screened differential lncRNAs and intersected with the selected miRNAs, and further predicted the target genes of the intersecting miRNAs by TargetScan, miRTarBase and miRDB, and intersected with the selected mRNAs. From the constructed ceRNA network we can see that Linc00460 may cause the invasion and metastasis of PTC by competitively inhibiting hsa-mir-150 and upregulating the expression of its downstream target gene EREG. Our study identied a series of lncRNAs associated with PTC progression and prognosis, and this complex ceRNA interaction network provides guidance for better understanding of the molecular mechanisms of PTC and can be used as an effective diagnostic tool for PTC invasion, metastasis and prognosis. Kaplan-Meier analysis of the differentially expressed RNAs associated with PTC pathogenesis conrmed that the lncRNAs

groups constitute the ceRNA network of PTC. Moreover, we used the microde database to predict the miRNAs that may be acted by the above screened differential lncRNAs and intersected with the selected miRNAs, and further predicted the target genes of the intersecting miRNAs by TargetScan, miRTarBase and miRDB, and intersected with the selected mRNAs. From the constructed ceRNA network we can see that Linc00460 may cause the invasion and metastasis of PTC by competitively inhibiting hsa-mir-150 and upregulating the expression of its downstream target gene EREG. Our study identi ed a series of lncRNAs associated with PTC progression and prognosis, and this complex ceRNA interaction network provides guidance for better understanding of the molecular mechanisms of PTC and can be used as an effective diagnostic tool for PTC invasion, metastasis and prognosis. Kaplan-Meier analysis of the differentially expressed RNAs associated with PTC pathogenesis con rmed that the lncRNAs AC097717.1, C20orf203, EMX2OS were potentially associated with the prognosis of PTC (P<0.05).

Introduction:
In recent years, the incidence of thyroid cancer has been increasing constantly, and thyroid cancer has become the most common and fastest growing malignant tumor of the endocrine system 1 , while papillary thyroid carcinoma (PTC) is the main type of thyroid cancer (about 80%) 2 .Although effective by thyroidectomy and hormonal therapy, a large number of patients still experience progressive disease and metastasis 3 .Therefore, exploring the molecular mechanisms of PTC and developing effective diagnostic and therapeutic targets are urgently needed.
As people's understanding of the genome deepped, it has been found that dysregulation of LncRNA expression plays a vital role in the progress of tumor 4 and is more complex and diverse in the regulation of biological processes (5). It is an RNA with a transcription length of more than 200 base pairs, no signi cant open reading frame, and no protein-coding function, but has been shown to be involved in many steps of gene transcription, including epigenetic and genetic levels (4), which mainly regulate cell growth, differentiation, and apoptosis (6), with both cancer-promoting and tumor-inhibiting effects. Therefore, LncRNA can be used as a recognized biomarker for the prognosis of these diseases.
Moreover, miRNA, a small ncRNA belongings to the same non-coding RNA (ncRNA), can regulate gene expression by attaching to miRNA-responsive elements of target transcription and at the posttranscriptional level (7), which is considered as an important regulator of cellular gene expression networks. In 2011, the advent of competitive endogenous RNAs (ceRNAs) introduced novel regulatory mechanisms driven by different non-coding RNAs, which mainly involve competitive cross-control of LncRNAs, mRNAs, miRNAs (8). They sponge miRNAs through common miRNA response elements (MREs), which in turn upregulate the expression of target RNAs, thereby weakening miRNA activity (9).
Currently, lncRNAs-miRNA-mRNA networks that may lead to tumor progression have been widely constructed in gastric cancer, breast cancer, liver cancer, and glioblastoma.
In recent years, an increasing number of LncRNAs have been used in the study of PTC, such as LncRNA H19 can upregulate the target gene YES1 by sponge miR-17-5p, which has been con rmed to play a ceRNA role (10). Another example is LncRNA Gas5 which can regulate PTEN expression in PTC by sponging miR-222-3p (11). However, genome-wide analysis of LncRNA-miRNA-mRNA regulatory networks for PTCs with large-scale samples is still lacking.
In summary, the aim of our study was to identify novel, speci c LncRNAs associated with PTC and to gain an insight into the characteristics of LncRNAs as biomarkers for the diagnosis and treatment of patients with papillary thyroid carcinoma, thereby further improving prognosis. After we obtained aberrantly expressed LncRNA, miRNA, and mRNA, and established an LncRNA-related PTC ceRNA network based on ceRNA theory, we found genes that express sequentially upregulated in the tissue, classical PTC tissue, and highly cellular PTC tissue, and concluded the possible mechanism of tumor malignant phenotype evolution of PTC.

Materials And Methods: Data Collection
We downloaded and extracted clinical data and level-3 data of RNA-Seq for PTC from the TCGA (The Cancer Genome Atlas) database. In addition, we extracted LncRNAs from RNA-Seq and exclude thyroid cancers of particular pathological types and preserve only PTC. The inclusion criteria were as follows: histological samples with complete analysis data, histological diagnosis of PTC without other malignancies; the exclusion criteria were as follows: tissue samples without complete data available for analysis, histological diagnosis isn't PTC, patients with other malignancies. For PTC samples containing detailed RNA expression data, we will also obtain their corresponding clinical details, such as survival time, survival status, gender, stage, etc.
Identi cation of DEmRNAs, DEmiRNAs, and DElncRNAs from the TCGA Database We Used the DEseq package in R to identify DEmRNAs, DEmiRNAs, and DElncRNAs that differentially expressed between PTC and nontumor tissues, Fold Change > 2 and P < 0.05 were considered signi cantly different and identi ed 501 PTC and 58 adjacent noncancerous tissue samples from the TCGA database for comprehensive bioinformatics analysis using online tools.

Construction of lncRNAs and mRNA co-expression networks
We designed lncRNA-mRNA networks as the basis for lncRNA-mediated regulation of mRNA abundance, and therefore, selected DElncRNAs and DEmRNAs by comparing PTC with normal tissues. To determine the respective RNA content, we set the selection criteria to Fold Change > 2 and P < 0.05, and lncRNAs -mRNA co-expression networks were constructed from the associated DElncRNAs and their respective target mRNAs.

Construction of ceRNA Network
TargetScan, miRTarBase and miRDB databases were used to predict target mRNAs that bound DEmiRNAs, and miRcode software was used to predict lncRNA binding to DEmiRNAs. To strengthen the correlation between mRNAs and DEmiRNAs, we cross-processed the dataset, and the lncRNA-miRNA-mRNA ceRNA network was formed on the basis of DEmiRNA-DElncRNA and DEmiRNA-DEmRNA interactions, and further visualized using Cytoscape software.

Survival analysis
Using the median expression of the differentially expressed RNAs as the cutoff value, patients with PTC were divided into high and low expression groups. The difference in survival time between the two groups was evaluated by Kaplan-Meier survival curve and log-rank test analysis.To screen prognostic signature DElncRNAs, mRNA, miRNA , the R package "survival" was used to analyze the correlation between lncRNAs and the survival time in the clinical data of PTC obtained from TCGA. P<.05 were considered as the cut-off value.

Result:
Differential expression of mRNA, lncRNA, miRNA in PTC Based on the cutoff threshold, a total of 3099 DEmRNAs (1976 upregulated and 1123 downregulated RNAs) was selected from both PTC and tissues adjacent to carcinoma. A total of 434 DElncRNAs (331 upregulated and 103 downregulated RNAs) and 149 DEmiRNAs (117 upregulated and 32 downregulated RNAs) were differentially expressed in PTC compared with tissues adjacent to carcinoma genes with an absolute fold change >2 and FDR value <0.05 were considered as discriminatively expressed. Heatmaps and volcanos were constructed to display the distinctive pattern of RNAs expressed in PTC compared with tissues adjacent to carcinoma (Figure 1). Thereafter, all differentially expressed RNAs were used for further analyses. Also, we used the online tool Venny for a comprehensive analysis of PTC and the concordant expression of DElncRNAs and DEmRNAs in TCGA (Figure 2).

ceRNA network in PTC
In the next step, we further understand the function of DElncRNAs by predicting the potential interactions between the above differential genes based on the ceRNA hypothesis through computational analysis ( Figure 3).

Discussion
Due to the lack of understanding in the pathogenesis of PTC and the lack of effective therapeutic targets, it is still di cult to improve the survival rate even after aggressive treatment such as maximum resection and radioiodine, and a considerable number of patients will experience poor prognosis such as recurrence or metastasis. A number of studies have reported that lncRNA has complex roles in the regulation of multiple signaling pathways and play an important role in the pathological process of tumors (12), which may become potential therapeutic targets (12)(13). However, there are few reports on the expression pro le of lncRNAs in PTC. To understand and identify potential therapeutic targets for PTC, we discussed the pathogenesis of PTC and use the interaction between differentially expressed RNA groups to constitute the ceRNA network of PTC. The ceRNA network hypothesis explains the relationship between lncRNAs, miRNAs, and mRNAs (14), nds genes that are sequentially up-regulated in tissues, classical PTC tissues, and highly cellular PTC tissues, and derives the possible mechanism of tumor malignant phenotype evolution of PTC.
LncRNAs can affect the stability of coding genes by regulating the activity of miRNAs, thereby regulating a variety of biological activities of tumor cells, including transcriptional regulation, post-transcriptional regulation, and phenotypic regulation (15).In previous studies, Zhou et al found that the expression of lncRNA CASC2 was associated with lymph node metastasis in PTC patients, which inhibited the proliferation, apoptosis and migration of cells (16).Also, Song et al con rmed that lncRNA ENST0000053963 could act as an oncogenic factor in PTC through MAPK signaling (17), suggesting that lncRNAs could function both as tumor suppressors and tumor promoters. In this study, according to the TCGA we newly discovered many lncRNAs aberrantly expressed in PTC tissues, of which Linc00460 may promote the invasion and metastasis of cancer, Linc00460 has been con rmed to be closely related to the progression of various cancers, such as its promotion of cell migration and invasion by regulating hnRNPK and induction of epithelial-mesenchymal transition in cell lung cancer. It promotes the occurrence of NPC by regulating the miR-149-5p/IL6 signaling pathway (18), which is compatible with our ndings. However, the speci c function and mechanism of Linc00460 remain unclear, and previous studies have shown that the up-regulation of Linc00460 expression is closely related to TNM stage and lymph node metastasis (19), which indicates that similar to other lncRNAs, Linc00460 may act by regulating cell proliferation and migration in PTC tissues, and involves the regulation of cell cycle in the mechanism of regulating proliferation.
On this basis, we further explored the mechanism by which Linc00460 promotes PTC tumorigenesis. In recent years, a large number of reports have shown that ceRNA networks are closely related to tumor occurrence, such as hepatocellular carcinoma, prostate cancer, and glioblastoma, and most lncRNAs and miRNAs can be used as members of ceRNA regulatory networks, which regulate the expression of their downstream targets by competitively binding to miRNAs. In this paper, we found by bioinformatics analysis that Linc00460 upregulates the expression of its downstream target gene EREG by competitively inhibiting hsa-mir-150 and may be used as a biomarker for PTC diagnosis and poor prognosis. Notably the present study is the rst to demonstrate that high expression of AC097717.1, C20orf203, EMX2OS is negatively associated with the survival time of patients with PTC and could be used as biomarkers for predicting the prognosis of the disease. The construction of lncRNA-miRNA-mRNA ceRNA network reveals the relationship between novel speci c lncRNAs and other RNAs in PTC. In addition, the close link between differentially expressed RNAs and the prognosis of PTC patients is being widely explored, so the results of this study have signi cant clinical implications. Figure 1 The volcano plot of the DEmRNA (A), DElncRNA (B), and DEmiRNA (C) and the heatmap of the DEmRNA (D), DElncRNA (E), and DEmiRNA (F).