Next Article in Journal
Using Y-Chromosomal Haplogroups in Genetic Association Studies and Suggested Implications
Next Article in Special Issue
Identification of Key Pathways and Genes in the Dynamic Progression of HCC Based on WGCNA
Previous Article in Journal
Isoform Sequencing and State-of-Art Applications for Unravelling Complexity of Plant Transcriptomes
Previous Article in Special Issue
SNCA Is a Functionally Low-Expressed Gene in Lung Adenocarcinoma
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Molecular Network-Based Identification of Competing Endogenous RNAs in Thyroid Carcinoma

1
College of Life Sciences, Zhejiang Sci-Tech University, Hangzhou 310018, China
2
School of Mathematics and Statistics, Hainan Normal University, Haikou 570100, China
3
Institute of Basic Medical Sciences, Wannan Medical College, Hefei 241000, China
4
Department of Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai, New York City, NY 10029, USA
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this study.
Genes 2018, 9(1), 44; https://doi.org/10.3390/genes9010044
Submission received: 12 December 2017 / Revised: 10 January 2018 / Accepted: 11 January 2018 / Published: 19 January 2018
(This article belongs to the Special Issue Computational Approaches for Disease Gene Identification)

Abstract

:
RNAs may act as competing endogenous RNAs (ceRNAs), a critical mechanism in determining gene expression regulations in many cancers. However, the roles of ceRNAs in thyroid carcinoma remains elusive. In this study, we have developed a novel pipeline called Molecular Network-based Identification of ceRNA (MNIceRNA) to identify ceRNAs in thyroid carcinoma. MNIceRNA first constructs micro RNA (miRNA)–messenger RNA (mRNA)long non-coding RNA (lncRNA) networks from miRcode database and weighted correlation network analysis (WGCNA), based on which to identify key drivers of differentially expressed RNAs between normal and tumor samples. It then infers ceRNAs of the identified key drivers using the long non-coding competing endogenous database (lnCeDB). We applied the pipeline into The Cancer Genome Atlas (TCGA) thyroid carcinoma data. As a result, 598 lncRNAs, 1025 mRNAs, and 90 microRNA (miRNAs) were inferred to be differentially expressed between normal and thyroid cancer samples. We then obtained eight key driver miRNAs, among which hsa-mir-221 and hsa-mir-222 were key driver RNAs identified by both miRNA–mRNA–lncRNA and WGCNA network. In addition, hsa-mir-375 was inferred to be significant for patients’ survival with 34 associated ceRNAs, among which RUNX2, DUSP6 and SEMA3D are known oncogenes regulating cellular proliferation and differentiation in thyroid cancer. These ceRNAs are critical in revealing the secrets behind thyroid cancer progression and may serve as future therapeutic biomarkers.

1. Introduction

Thyroid cancer is the most common malignancy of endocrine organs with enormous heterogeneity in terms of morphological features and prognosis [1]. Although most thyroid carcinomas tend to be biologically indolent and have a good prognosis, there are a few associated with more aggressive clinical manifestations [2]. Thyroid carcinoma is popularly classified into four classes including anaplastic thyroid carcinoma (ATC), follicular thyroid carcinoma (FTC), medullary thyroid carcinoma (MTC), and papillary thyroid carcinoma (PTC) [3], among which PTC is most common [4,5]. In 2017, there are approximately 56,870 new thyroid cancer incidences, representing about 3.4% of all new cancer cases worldwide [6]. Similarly, in the United States, the incidences of thyroid carcinoma increased steadily by around 6.6% each year from 2002 to 2009, which is the highest among all cancers [7]. In fact, thyroid cancer has become the fifth most common cancer in women [7,8].
Even though thyroid cancer harbors several highly universal genetic alterations, some of which are unique to this cancer [3], it is still very challenging to predict this cancer due to the complex disease progression process and complicated molecular interactions involved in it. Over the past decades, numerous studies have been performed to predict this cancer based on molecular, morphological, and immunological features [9], most of which focused on the detection of cancer-related protein-coding genes. However, it is known that protein-coding RNAs only cover approximately 2% of the total transcripts in mammalian [10], which urges the need to study the functions of non-coding RNAs (ncRNAs), especially long ncRNAs (lncRNAs) [11].
Recently, significant progresses have been achieved in exploring lncRNA biology [12]. For example, the lncRNA PVT1 was shown to be up-regulated in PTC, FTC, and ATC; MALAT1 was inferred to be involved in the regulation of cell cycle and migration [13]. However, most studies only focus on a subset of lncRNAs with specific regulatory mechanisms, while less is known on a transcriptome wide scale [14]. Moreover, the interaction mechanism between many kinds of RNAs are still elusive. Recent studies suggest that RNAs regulate each other’s expression levels by competing for a limited pool of microRNAs (miRNAs) in some circumstances [15,16]. In particular, Poliseno and his colleagues proposed a hypothesis: there exists an intricate post-transcriptional regulatory network mediated by miRNAs, in which non-coding RNAs and protein-coding RNAs compete for binding to miRNAs and regulate each other’s expression via sharing one or more miRNA response elements (MREs) [15]. The ncRNAs and protein-coding RNAs are called competing endogenous RNAs (ceRNAs) in the hypothesis.
ceRNA hypothesis demonstrates a new level of post-transcriptional regulation [17]. Given that even complete relief from repression by a miRNA usually has only mild effects on an individual mRNA, this theory highlights the importance of sharing binding sites for different miRNAs to yield substantial crosstalk [18,19,20]. ceRNAs may play a major role in certain dis-regulated or transient cellular states. For instance, it has been shown that the expression of tumor-suppressor gene PTEN can be regulated by its miRNA-mediated competitors VAPA, CNOT6L, SERINC1 or ZNF460 [21]. Especially, such mechanisms seem to be of particular relevance in cancer. For example, the lncRNA linc-MD1 has been shown to regulate the skeletal muscle cell differentiation clock by sponging miRNAs from its competitors, thereby enacting a ceRNA mechanism. In this ceRNA mechanism, MAML1 and MEF2C compete with linc-MD1 for miR-133 and miR-135, respectively [22,23]. Nevertheless, it is very difficult to build the exact ceRNA network and use it to understand RNA competing mechanisms. Fortunately, there are many well-established RNA databases such as long-non-coding RNA-associated diseases (LncRNADisease) database [24], the Human miRNA Disease Database (HMDD) [25], and database of Differentially Expressed MiRNAs in human Cancers dbDEMC [26]. In addition, miRNA-target interaction databases including miRcode [27] and miRanda [28,29,30,31,32,33,34,35], and ceRNA databases such as long non-coding competing endogenous database (lnCeDB) [36] have been developed, which provide much useful information.
In this study, we develop a novel pipeline called Molecular Network-based Identification of ceRNA (MNIceRNA) to identify ceRNAs in thyroid carcinoma. MNIceRNA first performs differential RNA analyses using edgeR [37], and then constructs gene co-expression and regulatory networks using machine learning based methods such as weighted correlation network analysis WGCNA [38] and known interaction data downloaded from databases such as miRcode. Finally, MNIceRNA focuses on thyroid carcinoma associated key driver genes (KDGs) and constructs a ceRNA network according to the lnCeDB [36]. The functions of the identified KDGs are also explored.

2. Materials and Methods

2.1. Data Collection and Pre-Processing

We downloaded RNA expression profiles of thyroid cancer and control samples from the Genomic Data Commons (GDC) data portal [39,40] and patient’s clinical information (see Table 1) from The Cancer Genome Atlas (TCGA) database [39,40]. Specifically, there are 559 samples used in this study, including 501 primary tumor samples and 58 solid tissue normal samples. The Genome research project of ENCyclopedia of DNA Elements (GENCODE) (GRCh38) (v25) catalogue (http://www.gencodegenes.org/) was used as a reference to quantify lncRNAs and mRNAs. In summary, 15,540 lncRNAs and 19,848 mRNAs from RNA-Sequencing (RNA-Seq) and 1881 miRNAs from miRNA-Seq were retrieved.

2.2. Differential Gene Expression Analysis

We applied edgeR to identify differentially expressed RNAs [37]. Specifically, the gene read counts were first processed with one scaling normalized factor from Trimmed Mean of M values (TMM) [41]. The negative binomial (NB) model was then applied to calculate the significance of RNAs, followed by an adjustment of p-values using the Benjamini–Hochberg method [42]. The cut-off values for significantly expressed RNAs were: (1) the false discovery rate (FDR, the adjusted p value) < 0.001; and (2) |log2 fold change (FC)| > 2 [5].

2.3. Construction of Gene Regulatory Network

We reconstructed the regulatory network using data combining lncRNAs, mRNAs, and miRNAs. The lncRNA–miRNA interactions and miRNA–mRNA interactions were downloaded from miRcode [27]. We then adopted a software called key driver analysis (KDA) [43] to identity key drivers in the regulatory network. Specifically, KDA takes a set of genes G and a directed gene network N as inputs. In our study, G is the differentially expressed RNAs and N is the regulatory network [44].

2.4. Construction of Gene Co-Expression Network

We inferred the co-expression network for a set of 32,209 RNAs including lncRNA, miRNA and mRNA using weighted gene co-expression network (WGCNA) algorithm [45], which was then visualized by Cytoscape 3.4.0 [12,46].

2.5. Survival Analysis

Survival analysis was performed using Cox proportional hazards regression models, with RNA expression in samples established as a binary variable. Because patient’s age and tumor stage have been interpreted to deeply influence molecular traits and clinical effect in thyroid cancer, we limited our initial cohort to primary tumor >Stage I patients. Specifically, we performed survival analyses for 9 (8 miRNA from miRcode, and 1 miRNA from WGCNA) miRNA, 80 lncRNA, and 190 mRNA (Table S1). p-values generated under this model were corrected for multiple-hypothesis testing using the Benjamini–Hochberg correction, with a significance threshold of FDR < 0.05 [47].

2.6. Function Enrichment

For clustering analysis of the significantly expressed RNAs, a pairwise complete-linkage hierarchical clustering method was employed to calculate the Euclidean distance. The results were shown using a heat map generated from the software packages cluster 3.0 [48] and TreeView [49]. In addition, differentially expressed genes (DEGs) were annotated by the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool (V6.8) [50,51], and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was also used to discover the potential pathways involved [12].

3. Results

3.1. Differentially Expressed RNAs between Primary Tumor and Control Samples

A total of 60,488 genes from RNA-Seq including 15,540 (26%) lncRNAs and 19,848 (33%) mRNAs were detected according to GENCODE (GRCh38) (v25) annotation. We also obtained 1881 miRNAs from miRNA-Seq. Similar to the genotype-tissue expression (GTEx) study [52], we performed a few data processing steps to require that genes to have at least 0.1 fragments per kilobase million (FPKM) in 2 or more individuals followed by quantile normalization across genes. In total, 14,848 lncRNAs, 16,575 mRNAs, and 786 miRNAs were left after the filtering. Based on the differential analyses, we inferred 463 up-regulated differentially expressed lncRNAs and 135 down-regulated differential lncRNAs, 812 up-regulated differential mRNAs and 213 down-regulated differential mRNAs, and 82 up-regulated differential miRNAs and eight down-regulated miRNAs respectively (Table S2). We selected a few top differential RNAs in each category and plotted their expression heat-maps in Figure 1.

3.2. Enriched Functions of Differentially Expressed RNAs

mRNAs: The functional enrichment analysis revealed that up-regulated DEGs were significantly enriched in the synthesis and degradation of extracellular matrix (ECM) related terms such as ECM organization (GO:0030198, FDR = 7.76 × 10 7 ) and ECM disassembly (GO:0022617, FDR = 0.05), collagen catabolic process (GO:0030574, FDR = 2.85 × 10 5 ), and cell adhesion (GO:0007155, FDR = 4.3 × 10 7 ) (Table S3). Based on the pathway enrichment analysis, the over-represented pathways of up-regulated DEGs was Neuroactive ligand-receptor interaction (hsa04080, FDR = 7.05 × 10 5 ) (Table S4). We also plotted the top enriched gene ontology (GO) and KEGG terms in Figure 2 for a better view. Our results are basically in line with Qiu et al. [5], especially ECM–receptor interaction, activation of MAPK activity and positive regulation of MAPK cascades [53].
lncRNAs: There are only 3 differentially expressed lnRNAs annotated in GO, namely RP11-161M6.2, ELFN2 and LINC00473. ELFN2 and LINC00473 were assigned to GO term “negative regulation of phosphatase activity (GO:0010923)” and “transcription, DNA-templated (GO:0006351)”, respectively. RP11-161M6.2, also called lipase maturation factor 1 (LMF1), played significant roles in many GO terms including triglyceride metabolic process (GO:0006641), endoplasmic reticulum (ER) to Golgi vesicle-mediated transport (GO:0006888), protein secretion (GO:0009306), protein glycosylation in Golgi (GO:0033578), chylomicron remnant clearance (GO:0034382), lipid digestion (GO:0044241), positive regulation of lipoprotein lipase activity (GO:0051006), protein maturation (GO:0051604), regulation of cholesterol metabolic process (GO:0090181), and regulation of triglyceride metabolic process (GO:0090207). In addition, there is only one lncRNA MIR205HG annotated in KEGG pathway “MiRNAs in cancer (hsa05206)”.

3.3. Key Driver Analysis

We have miRNA targets on 79,940 lncRNAs and 375,324 mRNAs and a total of 455,264 interaction pairs were collected from miRcode. We called this regulatory network the Database Network. We then mapped DEGs onto the Database Network, and performed key driver analysis to infer key genes driving the DEGs. As a result, we identified eight key driver miRNAs: hsa-mir-507, hsa-mir-375, hsa-mir-31, hsa-mir-144, hsa-mir-221, hsa-mir-222, hsa-mir-184, and hsa-mir-187. As annotated in DAVID database, hsa-mir-375, hsa-mir-31, hsa-mir-221, hsa-mir-222 and hsa-mir-184 were miRNAs in cancer, among which hsa-mir-184 was also relevant to endothelial dystrophy-iris hypoplasia-congenital cataract-stromal thinning (EDICT) syndrome. Interestingly, Zhang et al. [13] found that hsa-mir-375 [55] was up-regulated in PTC and MTC; hsa-mir-144 was up-regulated in PTC; hsa-mir-187 was up-regulated in PTC and FTC; and hsa-mir-221 [56,57] and hsa-mir-222 [56,58] were up-regulated in PTC, FTC, and ATC [55], which confirm the intimate association of our key drivers to thyroid cancer.
To study the influence of network structure to MNIceRNA, we reconstructed the co-expression network for the set of 32,209 RNAs using WGCNA. Sixty-five modules were identified, among which blue module is contains the largest number of DEGs. These genes of module were enriched in oxidation-reduction process and thyroid hormone synthesis. We then used this network to infer key driver genes and acquired 273 key drivers, including three miRNAs, 190 mRNAs and 80 lncRNAs, among which two (out of three) miRNAs, hsa-mir-221 and hsa-mir-222, overlapped with previous key driver miRNAs. The miR221-222 cluster is located on the X chromosome. Many previous studies have shown that the miR221-222 cluster was in the downstream of the MAPK pathway and involved in the regulation of cell cycle and apoptosis [13]. They were well known for deregulation in various malignancies and were among the first group of miRNAs shown to be deregulated in thyroid carcinoma [58,59]. In terms of mechanism, it was shown that miR221-222 cluster plays functions in PTC through negatively regulating p27 [60] and in ATC through their interaction with p27, RECK, and PTEN [59]. Similar other malignancies, up-regulation of miR221-222 cluster was associated with increased treatment resistance and recurrence rate, worse prognosis, and more invasive disease course [59,61,62]. It was likely that their tumorigenic property in thyroid was associated with their functions in tumor invasion and epithelial-mesenchymal transition (EMT) as shown in other malignancies [61,63]. They also act as potential biomarkers of thyroid malignancy with worse prognosis.

3.4. Competing Endogenous RNA Network Reveals Competing Endogenous Mechanisms of Long Non-Coding RNAs and Messenger RNAs

We focused on KDGs and constructed a ceRNA network according to lnCeDB database. As a result, 97 ceRNAs were obtained, most of which were also KDGs identified using the WGCNA network. Specifically, there were 34 ceRNAs related to hsa-mir-375, among which 11 were significant in the survival analysis (Figure 3). As examples, RUNX2 and SEMA3D are known oncogenes [5] and DUSP6 is critical for PTC and the MAPK pathway [64,65].

3.5. Survival Analysis of Key Driver Genes

We studied the association of key driver RNAs with patients’ survival, which can be used to evaluate their prognostic potential, and plotted in Figure 4 the Kaplan–Meier overall survival curves for miR-375 and a few clinical traits including cancer stage, gender and age. As a result, Figure 4A shows that cancer stage is significantly associated with overall survival with a p-value of 0.002, while gender and age (59 and above) were not significantly associated. Interestingly, Figure 4E,F shows that miR-375 is significantly associated to survival (p-value 0.03), indicating that it might be a prognostic maker and drug target for thyroid cancer. In addition, we also identified 44 significant mRNAs and 13 significant lncRNAs (p value < 0.05) in the survival analyses. Among them, lncRNAs RP5-1024C24.1 and LINC01539 were verified to be differentially expressed in PTC and exhibit specific topological characteristics in the lncRNA–mRNA co-expression network [66]. LncRNA RP5-1024C24.1 is an antisense transcript of metallophosphoesterase domain containing 2 (MPPED2). A previous study has shown that MPPED2 functions as a tumor suppresser in neuroblastoma tumorigenesis, and thus the low expression of RP5-1024C24.1 might promote tumor progression [66]. In addition, LINC01539 are deregulated in PTC [13].

4. Discussion

MNIceRNA combines information of RNA co-expression, regulation, and competing endogenous mechanisms to identify potential thyroid cancer biomarkers. It infers highly meaningful biomarkers as validated by literature mining and survival analyses. Specifically, MNIceRNA identified a cohort of 463 up-regulated differential lncRNAs, 812 differential mRNAs, and 82 differential miRNAs, which were mainly enriched in the ECM organization and degradation pathway. Similarly, there were 135 down-regulated differential lncRNAs, 213 down-regulated differential mRNAs, and eight down-regulated differential miRNAs, mainly enriched in cancer (hsa05206). Five up-regulated mRNAs, namely ELF3, HMGA2, LCN2, MET, and RUNX2, are known oncogenes [5]. In addition, as the up-regulated differential RNA exhibiting the most TF activity, PLAU encodes a serine protease that acts as an activator in the ECM degradation in tumor development [67]. PLAU also plays crucial roles in tumor invasion and metastasis [68]. Its overexpression has been found in cancer associated fibroblasts, which contribute to the tumor growth and progression [69]. Recently, a few studies suggested that PLAU was induced by PTC [70]. In contrast, a set of down-regulated mRNAs such as EGR2, GPC3, IGFBPL1, LRP1B, NR4A3, and PROX1 were identified as TSGs [5]. As a vital TF exhibiting higher transcriptional activity in PTC samples than in control, EGR2 functions as a tumor suppressor that is generally decreased in numerous cancer types such as ovarian [71] and gastric cancer [72]. These collectively suggest that PLAU might server as a significant TF promoting PTC progression, while EGR2 might be a tumor suppressor. However, more experimental validations are needed [5]. As for lncRNA, AC007255.8 is an antisense transcript of proline rich 15 (PRR15), which is overexpressed in advanced stage human colorectal cancer [73] and its expression was correlated with patient age [66]. AC079630.2 is mostly enriched in the pathway “Transcriptional misregulation in cancer”. Luo et al. showed that AC079630.2 exhibited high diagnostic ability to distinguish normal tissue and PTC tissue [74]. HOX transcript antisense RNA (HOTAIR) has been shown to be deregulated in a great number of human cancers such as oral cancer, nasopharynx, breast, esophagus, lung, liver, pancreas, colon, endometrium and cervix [75]. A study based on the TCGA and Gene Expression Omnibus (GEO) showed that HOTAIR was associated with poor survival of thyroid cancer patients [76].
In addition, the KDA analysis revealed eight key miRNAs, hsa-mir-375/507/222/221/187/184/144/31, among which hsa-mir-375 was significantly associated with patients’ survival. Previous studies suggested that overexpression of miR-375 was frequently observed in thyroid cancer patients and thus it might play a critical oncogenic role and be a potential diagnostic biomarker in thyroid cancer [77,78,79,80,81]. For example, the expression level of miR-375 was comparable between primary tumors and matched lymph node metastases, suggesting that its expression pattern in nodal metastases could prominently reflect that of the primary tumor in MTC [81]. Lassalle et al. found that the up-regulation of miR-375 was accompanied by reduced cell growth and synergistically improved sensitivity to vandetanib [80]. They also observed that miR-375 could enhance PARP cleavage and decline AKT phosphorylation, which play vital roles in cell growth. In addition, Wang et al. indicated that overexpression of miR-375 could inhibit the PTC cells proliferation and this inhibition was caused by the induction of cell apoptosis [82]. Moreover, it was found that miR-375 was upregulated by more than 35-fold in Follicular Variant Papillary Thyroid Carcinoma (FVPTC) and 20-fold in classic PTC, and not upregulated in normal thyroid tissue, hyperplastic nodules, and follicular carcinomas [83]. Thus, hsa-mir-375 is associated with early- and late-stage malignant progression and might be a novel clinical biomarker for thyroid cancer.
Moreover, using the lnCeDB database, MNIceRNA identified 34 ceRNAs (including six lncRNAs and 28 mRNAs) associated with hsa-mir-375, among which 11 were significant in the survival analysis. As examples, microarray analysis revealed that hypothyroidism induces significant reductions in KCNK2 transcripts [84]. CLDN1 is an essential molecule in tight junctions, CLDN1 was highly expressed in normal thyroid epithelium, but reduced in Hashimoto's Thyroiditis (HT) injured thyroid epithelial cells [85]. In addition, SEMA3D had superior diagnostic accuracy independently of the cytology in six datasets including The Cancer Genome Atlas (TCGA) thyroid dataset. This gene exhibited differences in the correlation coefficients between benign and malignant samples and could be an effective clinical biomarker for diagnosis of thyroid cancer [86]. Similarly, Steffen et al. [87] assembled a miRNA–protein target network for 677 human miRNAs and 18,880 targets which are listed in the TargetScan (http://www.targetscan.org). In this list, we found that hsa-mir-507 has a corresponding protein complex, corum-1422, and hsa-mir-221-222 cluster regulates nine protein complexes, of which ITGB1 [88], ITGB3 [88] and CD47 [89] were associated with the thyroid cancer. However, none of these protein complexes were putative ceRNAs we predicted.
It is worth mentioning that, although MNIceRNA identified many literature- and experiment-validated miRNAs, lncRNAs and mRNAs associated with thyroid cancer, the extent and relevance of ceRNA effect in vivo is still poorly understood. Recent experimental studies have suggested that the miRNA-mediated competition between ceRNAs could constitute an additional level of posttranscriptional regulation, playing important roles in many biological contexts. The sensitivity analysis shows that binding free energy and repression mechanisms are key ingredients for cross-talk between ceRNAs to arise. Interactions emerging in specific ranges of repression values, can be symmetrical (one ceRNA influences another and vice versa) or asymmetrical (one ceRNA influences another but not the reverse), and may be highly selective, while possibly be limited by noise [90]. On the other hand, the reporter assays demonstrated that [91] only active miRNA families with low total miRNA:target ratios are susceptible to ceRNA inductions even up to approximately 10,000 additional target copies per cell. In summary, there were many criteria for validating ceRNA, such as quantitative measurements of miRNA and target abundance [92], miRNA concentration and the size and affinities of the competing target pool [91], miRNA:target ratio, the absolute concentration of the effective target pool [91], timescales, steady-state, kinetic parameters [23,93], cellular concentrations of RBPs and miRNAs [20], and so on. In addition, PTR model was introduced in Figliuzzi et al. [90,93] by characterizing the transient response of the system to perturbations to validate ceRNA. The method was divided into two parts. The first part focuses on small perturbations by analyzing (in Fourier space) the linearized dynamics of a system of N ceRNAs jointly targeted by a single miRNA species. The second part focuses instead on large perturbations by using numerical analysis and analytical estimations to characterize the emergence of nonlinear response. The PTR model can be used to analyze the stoichiometric relationship of miR-375 and its target sites by manipulating TA through controlled expression of a validated target of miR-375 in thyroid [92]. Finally, the titration mechanism [93], argonaute individual-nucleotide resolution cross-linking and immune-precipitation (iCLIP) [91], Hermes systematically [94], stochastic model [95], Gillespie algorithm [95,96] and dynamical method [93] could also be used to validate ceRNAs. In the future, we will adopt these methods and perform experiments to validate the ceRNAs we predicted.
Taking together, our study identified many literature-validated RNAs critical to thyroid cancer progression and proposed a few novel RNAs to function as competing endogenous RNAs for thyroid cancer. However, we are fully aware that the limited sample size and information on miRNA–lncRNA–mRNA interactions might restrict the power of our conclusions. More experimental validations are suggested to confirm the contribution of our proposed RNAs in thyroid cancer.

5. Conclusions

In summary, we proposed a more multifaceted approach to construct ceRNAs network, and identified a set of crucial genes that could be used as biomarkers for thyroid carcinoma therapy, such as hsa-mir-375, AC012668.2, and SEMA3D. In the future, more attention should be paid to the construction of ceRNA networks and the validation of biomarkers or RNA competing endogenous interactions.

Supplementary Materials

The following are available online at www.mdpi.com/2073-4425/9/1/44/s1. Table S1: Survival analysis of key driver RNAs, Table S2: Up- and down-regulated differential RNAs, Table S3: Function enrichment of differential expressed RNAs on Gene ontology biological processes, Table S4: Function enrichment of differential expressed RNAs on KEGG pathways.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (61762035 to Y.Y.), and research grants (LY18F020027 to Y.Y.) from the Zhejiang Provincial Natural Science Foundation of China.

Author Contributions

J.Y. and Y.Y. conceived the concept of the work. M.L. and X.X. performed the analysis. M.L., J.Y. and Y.Y. wrote the manuscript. All authors helped to discuss and improve the work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Nikiforov, Y.E.; Nikiforova, M.N. Molecular genetics and diagnosis of thyroid cancer. Nat. Rev. Endocrinol. 2011, 7, 569–580. [Google Scholar] [CrossRef] [PubMed]
  2. Livolsi, V.A. Papillary thyroid carcinoma: An update. Mod. Pathol. 2011, 24 (Suppl. 2), S1–S9. [Google Scholar] [CrossRef] [PubMed]
  3. Xing, M. BRAF mutation in thyroid cancer. Endocr. Relat. Cancer 2005, 12, 245–262. [Google Scholar] [CrossRef] [PubMed]
  4. Agrawal, N.; Akbani, R.; Aksoy, B.A.; Ally, A.; Arachchi, H.; Asa, S.L.; Auman, J.T.; Balasundaram, M.; Balu, S.; Baylin, S.B. Integrated genomic characterization of papillary thyroid carcinoma. Cell 2014, 159, 676–690. [Google Scholar] [CrossRef] [PubMed]
  5. Qiu, J.; Zhang, W.; Xia, Q.; Liu, F.; Li, L.; Zhao, S.; Gao, X.; Zang, C.; Ge, R.; Sun, Y. RNA sequencing identifies crucial genes in papillary thyroid carcinoma (PTC) progression. Exp. Mol. Pathol. 2016, 100, 151–159. [Google Scholar] [CrossRef] [PubMed]
  6. Yapa, S.; Mulla, O.; Green, V.; England, J.; Greenman, J. The Role of Chemokines in Thyroid Carcinoma. Thyroid 2017, 27, 1347–1359. [Google Scholar] [CrossRef] [PubMed]
  7. National Cancer Institute. SEER Cancer Statistics Review 1975–2009. National Cancer Institute: Bethesda, MD, USA. Available online: http://seer.cancer.gov/archive/csr/1975_2009_pops09 (accessed on 12 December 2017).
  8. Udelsman, R.; Zhang, Y. The epidemic of thyroid cancer in the United States: The role of endocrinologists and ultrasounds. Thyroid 2014, 24, 472–479. [Google Scholar] [CrossRef] [PubMed]
  9. Matson, D.R.; Hardin, H.; Buehler, D.; Lloyd, R.V. AKT activity is elevated in aggressive thyroid neoplasms where it promotes proliferation and invasion. Exp. Mol. Pathol. 2017, 103, 288–293. [Google Scholar] [CrossRef] [PubMed]
  10. Birney, E.; Stamatoyannopoulos, J.A.; Dutta, A.; Guigó, R.; Gingeras, T.R.; Margulies, E.H.; Weng, Z.; Snyder, M.; Dermitzakis, E.T.; Thurman, R.E. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature 2007, 447, 799–816. [Google Scholar] [CrossRef] [PubMed]
  11. Zhang, Y.; Tao, Y.; Liao, Q. Long non-coding RNA: A crosslink in biological regulatory network. Brief. Bioinform. 2017. [Google Scholar] [CrossRef]
  12. Wang, Q.; Yang, H.; Wu, L.; Yao, J.; Meng, X.; Jiang, H.; Xiao, C.; Wu, F. Identification of Specific Long Non-Coding RNA Expression: Profile and Analysis of Association with Clinicopathologic Characteristics and BRAF Mutation in Papillary Thyroid Cancer. Thyroid 2016, 26, 1719–1732. [Google Scholar] [CrossRef] [PubMed]
  13. Zhang, R.; Hardin, H.; Chen, J.; Guo, Z.; Lloyd, R.V. Non-Coding RNAs in Thyroid Cancer. Endocr. Pathol. 2016, 27, 12–20. [Google Scholar] [CrossRef] [PubMed]
  14. Nagano, T.; Fraser, P. No-nonsense functions for long non-coding RNAs. Cell 2011, 145, 178–181. [Google Scholar] [CrossRef] [PubMed]
  15. Poliseno, L.; Salmena, L.; Zhang, J.; Carver, B.; Haveman, W.J.; Pandolfi, P.P. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature 2010, 465, 1033–1038. [Google Scholar] [CrossRef] [PubMed]
  16. Seitz, H. Redefining MicroRNA Targets. Curr. Biol. 2009, 19, 870–873. [Google Scholar] [CrossRef] [PubMed]
  17. Xia, T.; Liao, Q.; Jiang, X.; Shao, Y.; Xiao, B.; Xi, Y.; Guo, J. Long non-coding RNA associated-competing endogenous RNAs in gastric cancer. Sci. Rep. 2014, 4, 6088. [Google Scholar] [CrossRef] [PubMed]
  18. Salmena, L.; Poliseno, L.; Tay, Y.; Kats, L.; Pandolfi, P.P. A ceRNA hypothesis: The rosetta stone of a hidden RNA language? Cell 2011, 146, 353–358. [Google Scholar] [CrossRef] [PubMed]
  19. Tay, Y.; Rinn, J.; Pandolfi, P.P. The multilayered complexity of ceRNA crosstalk and competition. Nature 2014, 505, 344–352. [Google Scholar] [CrossRef] [PubMed]
  20. Jens, M.; Rajewsky, N. Competition between target sites of regulators shapes post-transcriptional gene regulation. Nat. Rev. Genet. 2015, 16, 113–126. [Google Scholar] [CrossRef] [PubMed]
  21. Tay, Y.; Kats, L.; Salmena, L.; Weiss, D.; Shen, M.T.; Ala, U.; Karreth, F.; Poliseno, L.; Provero, P.; Cunto, F.D. Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. Cell 2011, 147, 344–357. [Google Scholar] [CrossRef] [PubMed]
  22. Cesana, M.; Cacchiarelli, D.; Legnini, I.; Santini, T.; Sthandier, O.; Chinappi, M.; Tramontano, A.; Bozzoni, I. A long non-coding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 2011, 147, 358–369. [Google Scholar] [CrossRef] [PubMed]
  23. Martirosyan, A.; Figliuzzi, M.; Marinari, E.; Martino, A.D. Probing the limits to microRNA-mediated control of gene expression. PLoS Comput. Biol. 2016, 12, e1004715. [Google Scholar] [CrossRef] [PubMed]
  24. Geng, C.; Wang, Z.; Wang, D.; Qiu, C.; Liu, M.; Xing, C.; Zhang, Q.; Yan, G.; Cui, Q. LncRNADisease: A database for long-non-coding RNA-associated diseases. Nucleic Acids Res. 2013, 41, D983. [Google Scholar]
  25. Li, Y.; Qiu, C.; Tu, J.; Geng, B.; Yang, J.; Jiang, T.; Cui, Q. HMDD v2.0: A database for experimentally supported human microRNA and disease associations. Nucleic Acids Res. 2014, 42, 1070–1074. [Google Scholar] [CrossRef] [PubMed]
  26. Yang, Z.; Ren, F.; Liu, C.; He, S.; Sun, G.; Gao, Q.; Yao, L.; Zhang, Y.; Miao, R.; Cao, Y. dbDEMC: A database of differentially expressed miRNAs in human cancers. BMC Genom. 2010, 11 (Suppl. 4), S5. [Google Scholar] [CrossRef] [PubMed]
  27. Ashwini, J.; Marks, D.S.; Erik, L. miRcode: A map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics 2012, 28, 2062–2063. [Google Scholar]
  28. Betel, D.; Koppal, A.; Agius, P.; Sander, C.; Leslie, C. Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol. 2010, 11, R90. [Google Scholar] [CrossRef] [PubMed]
  29. Betel, D.; Wilson, M.; Gabow, A.; Marks, D.S.; Sander, C. The microRNA.org resource: Targets and expression. Nucleic Acids Res. 2008, 36, D149. [Google Scholar] [CrossRef] [PubMed]
  30. Enright, A.J.; John, B.; Gaul, U.; Tuschl, T.; Sander, C.; Marks, D.S. MicroRNA targets in Drosophila. Genome Biol. 2003, 5, R1. [Google Scholar] [CrossRef] [PubMed]
  31. Hofacker, I.L.; Fontana, W.; Stadler, P.F.; Bonhoeffer, L.S.; Tacker, M.; Schuster, P. Fast folding and comparison of RNA secondary structures. Monatshefte für Chemie 1994, 125, 167–188. [Google Scholar] [CrossRef]
  32. John, B.; Enright, A.J.; Aravin, A.; Tuschl, T.; Sander, C.; Marks, D.S. Human MicroRNA targets. PLoS Biol. 2004, 2, e363. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Landgraf, P.; Rusu, M.; Sheridan, R.; Sewer, A.; Iovino, N.; Aravin, A.; Pfeffer, S.; Rice, A.; Kamphorst, A.O.; Landthaler, M. A Mammalian microRNA Expression Atlas Based on Small RNA Library Sequencing. Cell 2007, 129, 1401–1414. [Google Scholar] [CrossRef] [PubMed]
  34. Mccaskill, J.S. The equilibrium partition function and base pair binding probabilities for RNA secondary structure. Biopolymers 1990, 29, 1105–1119. [Google Scholar] [CrossRef] [PubMed]
  35. Zuker, M.; Stiegler, P. Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information. Nucleic Acids Res. 1981, 9, 133–148. [Google Scholar] [CrossRef] [PubMed]
  36. Das, S.; Ghosal, S.; Sen, R.; Chakrabarti, J. lnCeDB: Database of human long non-coding RNA acting as competing endogenous RNA. PLoS ONE 2014, 9, e98965. [Google Scholar] [CrossRef] [PubMed]
  37. Robinson, M.D.; Mccarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed]
  38. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted gene co-expression network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef] [PubMed]
  39. GDC Data Portal. Available online: http://portal.gdc.cancer.gov/projects/TCGA-THCA (accessed on 21 December 2016).
  40. The Cancer Genome Atlas. Available online: http://cancergenome.nih.gov/ (accessed on 12 December 2017).
  41. Fang, S.M.; Hu, B.L.; Zhou, Q.Z.; Yu, Q.Y.; Zhang, Z. Comparative analysis of the silk gland transcriptomes between the domestic and wild silkworms. BMC Genom. 2015, 16, 60. [Google Scholar] [CrossRef] [PubMed]
  42. Benjamini, Y.; Hochberg, Y. Controlling The False Discovery Rate—A Practical And Powerful Approach to Multiple Testing. J. R. Stat. Soc. 1995, 57, 289–300. [Google Scholar]
  43. Zhang, B.; Zhu, J. Identification of key causal regulators in gene networks. Lect. Notes Eng. Comput. Sci. 2013, 2205, 1309–1312. [Google Scholar]
  44. Yang, J.; Qiu, J.; Wang, K.; Zhu, L.; Fan, J.; Zheng, D.; Meng, X.; Yang, J.; Peng, L.; Fu, Y. Using molecular functional networks to manifest connections between obesity and obesity-related diseases. Oncotarget 2017, 8, 85136–85149. [Google Scholar] [CrossRef] [PubMed]
  45. Zhang, B.; Horvath, S. A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 2004, 4, 17. [Google Scholar] [CrossRef] [PubMed]
  46. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef] [PubMed]
  47. Zou, A.E.; Zheng, H.; Saad, M.A.; Rahimy, M.; Ku, J.; Kuo, S.Z.; Honda, T.K.; Wangrodriguez, J.; Xuan, Y.; Korrapati, A. The non-coding landscape of head and neck squamous cell carcinoma. Oncotarget 2016, 7, 51211–51222. [Google Scholar] [CrossRef] [PubMed]
  48. De Hoon, M.J.; Imoto, S.; Nolan, J.; Miyano, S. Open Source Clustering Software. Bioinformatics 2004, 20, 1453–1454. [Google Scholar] [CrossRef] [PubMed]
  49. Page, R.D. TREEVIEW: An application to display phylogenetic trees on personal computers. Comput. Appl. Biosci. 1996, 12, 357–358. [Google Scholar] [PubMed]
  50. Huang, D.W.; Sherman, B.T.; Lempicki, R.A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 2009, 4, 44–57. [Google Scholar] [CrossRef] [PubMed]
  51. Huang, D.W.; Sherman, B.T.; Lempicki, R.A. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009, 37, 1–13. [Google Scholar] [CrossRef] [PubMed]
  52. Consortium, G. The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans. Science 2015, 348, 648–660. [Google Scholar] [CrossRef] [PubMed]
  53. Murugan, A.K.; Munirajan, A.K.; Alzahrani, A.S. Long non-coding RNAs: Emerging players in thyroid cancer pathogenesis. Endocr. Relat. Cancer 2017. [Google Scholar] [CrossRef]
  54. Wu, D.; Wang, B.; Shang, J.; Song, J.; Zhang, H. miR-31 Reduces Cell Growth of Papillary Thyroid Carcinoma by RNA-Binding Protein HuR. Clin. Lab. 2015, 61, 1625–1634. [Google Scholar] [CrossRef] [PubMed]
  55. Amaral, P.P.; Dinger, M.E.; Mercer, T.R.; Mattick, J.S. The Eukaryotic genome as an RNA machine. Science 2008, 319, 1787–1789. [Google Scholar] [CrossRef] [PubMed]
  56. Venkatesh, T.; Suresh, P.S.; Tsutsumi, R. Non-coding RNAs: Functions and applications in endocrine-related cancer. Mol. Cell. Endocrinol. 2015, 416, 88–96. [Google Scholar] [CrossRef] [PubMed]
  57. Vester, B.; Wengel, J. LNA (locked nucleic acid): High-affinity targeting of complementary RNA and DNA. Biochemistry 2004, 43, 13233–13244. [Google Scholar] [CrossRef] [PubMed]
  58. Guo, Z.; Hardin, H.; Montemayor-Garcia, C.; Asioli, S.; Righi, A.; Maletta, F.; Sapino, A.; Lloyd, R.V. In Situ Hybridization Analysis of miR-146b-5p and miR-21 in Thyroid Nodules: Diagnostic Implications. Endocr. Pathol. 2015, 26, 157–163. [Google Scholar] [CrossRef] [PubMed]
  59. Fuziwara, C.S.; Kimura, E.T. MicroRNA deregulation in anaplastic thyroid cancer biology. Int. J. Endocrinol. 2014, 2014, 743450. [Google Scholar] [CrossRef] [PubMed]
  60. Visone, R.; Russo, L.; Pallante, P.; De, M.I.; Ferraro, A.; Leone, V.; Borbone, E.; Petrocca, F.; Alder, H.; Croce, C.M. MicroRNAs (miR)-221 and miR-222, both overexpressed in human thyroid papillary carcinomas, regulate p27Kip1 protein levels and cell cycle. Endocr. Relat. Cancer 2007, 14, 791–798. [Google Scholar] [CrossRef] [PubMed]
  61. Pallante, P.; Battista, S.; Pierantoni, G.M.; Fusco, A. Deregulation of microRNA expression in thyroid neoplasias. Nat. Rev. Endocrinol. 2014, 10, 88–101. [Google Scholar] [CrossRef] [PubMed]
  62. Chruścik, A.; Lam, A.K. Clinical pathological impacts of microRNAs in papillary thyroid carcinoma: A crucial review. Exp. Mol. Pathol. 2015, 99, 393–398. [Google Scholar] [CrossRef] [PubMed]
  63. Ogata, H.; Goto, S.; Sato, K.; Fujibuchi, W.; Bono, H.; Kanehisa, M. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 27, 29–34. [Google Scholar] [CrossRef]
  64. Zhang, S.; Wang, Y.; Chen, M.; Sun, L.; Han, J.; Elena, V.K.; Qiao, H. CXCL12 methylation-mediated epigenetic regulation of gene expression in papillary thyroid carcinoma. Sci. Rep. 2017, 7, 44033. [Google Scholar] [CrossRef] [PubMed]
  65. Lee, E.K.; Chung, K.W.; Yang, S.K.; Park, M.J.; Min, H.S.; Kim, S.W.; Kang, H.S. DNA methylation of MAPK signal-inhibiting genes in papillary thyroid carcinoma. Anticancer Res. 2013, 33, 4833–4839. [Google Scholar] [PubMed]
  66. Du, Y.; Xia, W.; Zhang, J.; Wan, D.; Yang, Z.; Li, X. Comprehensive analysis of long non-coding RNA-mRNA co-expression patterns in thyroid cancer. Mol. Biosyst. 2017, 13, 2107–2115. [Google Scholar] [CrossRef] [PubMed]
  67. Choi, H.S.; Song, M.; Song, M.K.; Kim, Y.J.; Ryu, J.C. Analysis of differentially expressed genes by Mirex ‘persistent organic pollutant’ in HepG2 cells. Toxicol. Environ. Health Sci. 2011, 3, 245–253. [Google Scholar] [CrossRef]
  68. Blasi, F.; Sidenius, N. The urokinase receptor: Focused cell surface proteolysis, cell adhesion and signaling. FEBS Lett. 2010, 584, 1923–1930. [Google Scholar] [CrossRef] [PubMed]
  69. Casey, T.M.E.J.; Crocker, A.; White, J.; Tessitore, J.; Stanley, M.; Harlow, S.; Bunn, J.Y.; Weaver, D.; Muss, H.; Plaut, K. Cancer associated fibroblasts stimulated by transforming growth factor beta1 (TGF-β1) increase invasion rate of tumor cells: A population study. Breast Cancer Res. Treat. 2008, 110, 39–49. [Google Scholar] [CrossRef] [PubMed]
  70. Delys, L.; Detours, V.; Franc, B.; Thomas, G.; Bogdanova, T.; Tronko, M.; Libert, F.; Dumont, J.E.; Maenhaut, C. Gene expression and the biological phenotype of papillary thyroid carcinomas. Oncogene 2007, 26, 7894–7903. [Google Scholar] [CrossRef] [PubMed]
  71. Unoki, M.; Nakamura, Y. Growth-suppressive effects of BPOZ and EGR2, two genes involved in the PTEN signaling pathway. Oncogene 2001, 20, 4457–4465. [Google Scholar] [CrossRef] [PubMed]
  72. Wu, Q.; Jin, H.; Yang, Z.; Luo, G.; Lu, Y.; Li, K.; Ren, G.; Su, T.; Pan, Y.; Feng, B. MiR-150 promotes gastric cancer proliferation by negatively regulating the pro-apoptotic gene EGR2. Biochem. Biophys. Res. Commun. 2010, 392, 340–345. [Google Scholar] [CrossRef] [PubMed]
  73. Meunier, D.; Patra, K.; Smits, R.; Hägebarth, A.; Lüttges, A.; Jaussi, R.; Wieduwilt, M.J.; Quintanilla-Fend, L.; Himmelbauer, H.; Fodde, R. Expression analysis of proline rich 15 (Prr15) in mouse and human gastrointestinal tumors. Mol. Carcinog. 2011, 50, 8–15. [Google Scholar] [CrossRef] [PubMed]
  74. Luo, Y.L.L.; He, R.; Wen, D.; Deng, G.; Yang, H.; He, Y.; Ma, W.; Cai, X.; Chen, J.; Chen, G. RNA-sequencing investigation identifies an effective risk score generated by three novel lncRNAs for the survival of papillary thyroid cancer patients. Oncotarget 2017, 8, 74139–74158. [Google Scholar] [CrossRef] [PubMed]
  75. Huarte, M. The emerging role of lncRNAs in cancer. Nat. Med. 2015, 21, 1253–1261. [Google Scholar] [CrossRef] [PubMed]
  76. Li, H.M.; Yang, H.; Wen, D.Y.; Luo, Y.H.; Liang, C.Y.; Pan, D.H.; Ma, W.; Chen, G.; He, Y.; Chen, J.Q. Overexpression of LncRNA HOTAIR is Associated with Poor Prognosis in Thyroid Carcinoma: A Study Based on TCGA and GEO Data. Horm. Metab. Res. 2017, 49, 388–399. [Google Scholar] [CrossRef] [PubMed]
  77. Hudson, J.; Duncavage, E.; Tamburrino, A.; Salerno, P.; Xi, L.; Raffeld, M.; Moley, J.; Chernock, R.D. Overexpression of mir-10a and mir-375 and downregulation of yap1 in medullary thyroid carcinoma. Exp. Mol. Pathol. 2013, 95, 62–67. [Google Scholar] [CrossRef] [PubMed]
  78. Mian, C.; Pennelli, G.; Fassan, M.; Balistreri, M.; Barollo, S.; Cavedon, E.; Galuppini, F.; Pizzi, M.; Vianello, F.; Pelizzo, M.R. MicroRNA profiles in familial and sporadic medullary thyroid carcinoma: Preliminary relationships with ret status and outcome. Thyroid 2012, 22, 890–896. [Google Scholar] [CrossRef] [PubMed]
  79. Gundara, J.S.; Zhao, J.T.; Gill, A.J.; Clifton-Bligh, R.; Robinson, B.G.; Delbridge, L.; Sidhu, S.B. Nodal metastasis microRNA expression correlates with the primary tumour in MTC. ANZ J. Surg. 2012, 84, 235–239. [Google Scholar] [CrossRef] [PubMed]
  80. Lassalle, S.; Zangari, J.; Popa, A.; Ilie, M.; Hofman, V.; Long, E.; Patey, M.; Tissier, F.; Belléannée, G.; Trouette, H. MicroRNA-375/SEC23A as biomarkers of the in vitro efficacy of vandetanib. Oncotarget 2016, 7, 30461–30478. [Google Scholar] [CrossRef] [PubMed]
  81. Shi, L.; Zhao, S.M.; Luo, Y.; Zhang, A.W.; Wei, L.H.; Xie, Z.Y.; Li, Y.Y.; Ma, W. Mir-375: A prospective regulator in medullary thyroid cancer based on microarray data and bioinformatics analyses. Pathol. Res. Pract. 2017, 213, 1344–1354. [Google Scholar] [CrossRef] [PubMed]
  82. Wang, X.Z.; Hang, Y.K.; Liu, J.B.; Hou, Y.Q.; Wang, N.; Wang, M.J. Over-expression of microRNA-375 inhibits papillary thyroid carcinoma cell proliferation and induces cell apoptosis by targeting ERBB2. J. Pharmacol. Sci. 2016, 130, 78–84. [Google Scholar] [CrossRef] [PubMed]
  83. Dettmer, M.; Perren, A.; Moch, H.; Komminoth, P.; Nikiforov, Y.E.; Nikiforova, M.N. Comprehensive microRNA expression profiling identifies novel markers in follicular variant of papillary thyroid carcinoma. Thyroid 2013, 23, 1383–1389. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  84. Sabrina, L.B.; Sophie, D.; Arnaud, C.; Chloé, B.; Franck, A.; Gilles, T.; Gilles, L.; Sepideh, S.; Isabelle, B.; Pond, A.L. Microarray analysis reveals complex remodeling of cardiac ion channel expression with altered thyroid status: Relation to cellular and integrated electrophysiology. Circ. Res. 2003, 92, 234–242. [Google Scholar]
  85. Zhu, J.; Zhang, Y.; Zhang, W.; Zhang, W.; Fan, L.; Wang, L.; Liu, Y.; Liu, S.; Guo, Y.; Wang, Y. MicroRNA-142-5p contributes to Hashimoto’s thyroiditis by targeting CLDN1. J. Transl. Med. 2016, 14, 166. [Google Scholar] [CrossRef] [PubMed]
  86. Gomezrueda, H.; Palacioscorona, R.; Gutiérrezhermosillo, H.; Trevino, V. A robust biomarker of differential correlations improves the diagnosis of cytologically indeterminate thyroid cancers. Int. J. Mol. Med. 2016, 37, 1355–1362. [Google Scholar] [CrossRef] [PubMed]
  87. Sass, S.; Dietmann, S.; Burk, U.; Brabletz, S.; Lutter, D.; Kowarsch, A.; Mayer, K.F.; Brabletz, T.; Ruepp, A.; Theis, F. MicroRNAs coordinately regulate protein complexes. BMC Syst. Biol. 2011, 5, 136. [Google Scholar] [CrossRef] [PubMed]
  88. Cockburn, J.G.; Richardson, D.S.; Gujral, T.S.; Mulligan, L.M. Ret-mediated cell adhesion and migration require multiple integrin subunits. J. Clin. Endocrinol. Metab. 2010, 95, 342–346. [Google Scholar] [CrossRef] [PubMed]
  89. Rath, G.M.; Schneider, C.; Dedieu, S.; Rothhut, B.; Soula-Rothhut, M.; Ghoneim, C.; Sid, B.; Morjani, H.; Btaouri, H.E.; Martiny, L. The C-terminal CD47/IAP-binding domain of thrombospondin-1 prevents camptothecin- and doxorubicin-induced apoptosis in human thyroid carcinoma cells. Biochim. Biophys. Acta Mol. Cell Res. 2006, 1763, 1125–1134. [Google Scholar] [CrossRef] [PubMed]
  90. Figliuzzi, M.; Marinari, E.; De, M.A. MicroRNAs as a selective channel of communication between competing RNAs: A steady-state theory. Biophys. J. 2013, 104, 1203–1213. [Google Scholar] [CrossRef] [PubMed]
  91. Bosson, A.D.; Zamudio, J.R.; Sharp, P.A. Endogenous miRNA and target concentrations determine susceptibility to potential ceRNA competition. Mol. Cell 2014, 56, 347–359. [Google Scholar] [CrossRef] [PubMed]
  92. Denzler, R.; Agarwal, V.; Stefano, J.; Bartel, D.P.; Stoffel, M. Assessing the ceRNA hypothesis with quantitative measurements of miRNA and target abundance. Mol. Cell 2014, 54, 766–776. [Google Scholar] [CrossRef] [PubMed]
  93. Figliuzzi, M.; De, M.A.; Marinari, E. RNA-based regulation: Dynamics and response to perturbations of competing RNAs. Biophys. J. 2014, 107, 1011–1022. [Google Scholar] [CrossRef] [PubMed]
  94. Sumazin, P.; Yang, X.; Chiu, H.S.; Chung, W.J.; Iyer, A.; Llobet-Navas, D.; Rajbhandari, P.; Bansal, M.; Guarnieri, P.; Silva, J.; et al. An extensive microRNA-mediated network of RNA-RNA interactions regulates established oncogenic pathways in glioblastoma. Cell 2011, 147, 370–381. [Google Scholar] [CrossRef] [PubMed]
  95. Bosia, C.; Pagnani, A.; Zecchina, R. Modelling competing endogenous RNA networks. PLoS ONE 2013, 8, e66609. [Google Scholar] [CrossRef] [PubMed]
  96. Martirosyan, A.; Marsili, M.; De, M.A. Translating ceRNA susceptibilities into correlation functions. Biophys. J. 2017, 113, 206–213. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Heat-maps on expressions of top differential RNAs for thyroid cancer: (A) differential long non-coding RNAs (lncRNAs); (B) differential micro RNAs (miRNAs); and (C) differential messenger RNAs (mRNAs). X-axis represents samples, while Y-axis represents the biological elements studied.
Figure 1. Heat-maps on expressions of top differential RNAs for thyroid cancer: (A) differential long non-coding RNAs (lncRNAs); (B) differential micro RNAs (miRNAs); and (C) differential messenger RNAs (mRNAs). X-axis represents samples, while Y-axis represents the biological elements studied.
Genes 09 00044 g001
Figure 2. Plot of the differentially expressed genes enriched GO and KEGG: (A) The plot of the enriched GO (biological process) with DEGs of the mRNAs. X-axis represents the percentage of enriched genes. (B) The plot of the enriched KEGG with DEGs of the mRNAs. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DEGs: Differentially expressed genes; FDR: False discovery rate. miRNAs: The differentially expressed miRNAs were significantly enriched in miRNAs in cancer (hsa05206, FDR = 4.44 × 10 14 ), which include hsa-mir-221, hsa-mir-222, hsa-mir-31, hsa-mir-34a, hsa-mir-373, hsa-mir-375, hsa-mir-451a, hsa-mir-483, hsa-mir-520a, hsa-mir-520c, hsa-mir-520g, and hsa-mir-520h. Among them, Wu et al. [54] found that miR-31 was significantly down-regulated in papillary thyroid carcinoma patients. Furthermore, down regulation of miR-31 increased the proliferation, migration, and invasion of ovarian carcinoma cells. In addition, they revealed that the human antigen R (HuR) was a target for miR-31 and knock down of HuR resulted in enhanced cell viability and decreased cell migration rate.
Figure 2. Plot of the differentially expressed genes enriched GO and KEGG: (A) The plot of the enriched GO (biological process) with DEGs of the mRNAs. X-axis represents the percentage of enriched genes. (B) The plot of the enriched KEGG with DEGs of the mRNAs. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DEGs: Differentially expressed genes; FDR: False discovery rate. miRNAs: The differentially expressed miRNAs were significantly enriched in miRNAs in cancer (hsa05206, FDR = 4.44 × 10 14 ), which include hsa-mir-221, hsa-mir-222, hsa-mir-31, hsa-mir-34a, hsa-mir-373, hsa-mir-375, hsa-mir-451a, hsa-mir-483, hsa-mir-520a, hsa-mir-520c, hsa-mir-520g, and hsa-mir-520h. Among them, Wu et al. [54] found that miR-31 was significantly down-regulated in papillary thyroid carcinoma patients. Furthermore, down regulation of miR-31 increased the proliferation, migration, and invasion of ovarian carcinoma cells. In addition, they revealed that the human antigen R (HuR) was a target for miR-31 and knock down of HuR resulted in enhanced cell viability and decreased cell migration rate.
Genes 09 00044 g002
Figure 3. Key drivers for Thyroid carcinoma in the long non-coding competing endogenous database (lnCeDB): (A) yellow nodes represent miRNAs (diamond) and lncRNAs (ellipse), while the purple represent mRNAs; and (B) yellow nodes represent lncRNAs, purple represent mRNAs, and diamond represent significant RNAs in the survival analysis from key driver analysis (KDA) of the weighted correlation network analysis (WGCNA). The above relationship are competing endogenous (ceRNAs) from the lnCeDB database.
Figure 3. Key drivers for Thyroid carcinoma in the long non-coding competing endogenous database (lnCeDB): (A) yellow nodes represent miRNAs (diamond) and lncRNAs (ellipse), while the purple represent mRNAs; and (B) yellow nodes represent lncRNAs, purple represent mRNAs, and diamond represent significant RNAs in the survival analysis from key driver analysis (KDA) of the weighted correlation network analysis (WGCNA). The above relationship are competing endogenous (ceRNAs) from the lnCeDB database.
Genes 09 00044 g003
Figure 4. Thyroid patients’ survival analyses on miR-375 and a few clinical traits: (A,B) survival outcomes of different stages and gender; (C,D) different age groups have different effects on over-all survival; and (E,F) survival outcomes according to relatively high and low expression, of which (E) represents survival outcomes of has-mir-375 in all samples, and (F) represents survival outcomes of has-mir-375 in higher than stage I tumor.
Figure 4. Thyroid patients’ survival analyses on miR-375 and a few clinical traits: (A,B) survival outcomes of different stages and gender; (C,D) different age groups have different effects on over-all survival; and (E,F) survival outcomes according to relatively high and low expression, of which (E) represents survival outcomes of has-mir-375 in all samples, and (F) represents survival outcomes of has-mir-375 in higher than stage I tumor.
Genes 09 00044 g004
Table 1. Clinical information of the 559 samples used in this study.
Table 1. Clinical information of the 559 samples used in this study.
CharacteristicNumbers
Sample typePrimary tumor501
Solid tissue normal58
AgeMedian47
Range [years]15~89
SexMale152
Female407
Vital statusAlive539
Dead20
StageI315
II59
III124
IV2

Share and Cite

MDPI and ACS Style

Lu, M.; Xu, X.; Xi, B.; Dai, Q.; Li, C.; Su, L.; Zhou, X.; Tang, M.; Yao, Y.; Yang, J. Molecular Network-Based Identification of Competing Endogenous RNAs in Thyroid Carcinoma. Genes 2018, 9, 44. https://doi.org/10.3390/genes9010044

AMA Style

Lu M, Xu X, Xi B, Dai Q, Li C, Su L, Zhou X, Tang M, Yao Y, Yang J. Molecular Network-Based Identification of Competing Endogenous RNAs in Thyroid Carcinoma. Genes. 2018; 9(1):44. https://doi.org/10.3390/genes9010044

Chicago/Turabian Style

Lu, Minjia, Xingyu Xu, Baohang Xi, Qi Dai, Chenli Li, Li Su, Xiaonan Zhou, Min Tang, Yuhua Yao, and Jialiang Yang. 2018. "Molecular Network-Based Identification of Competing Endogenous RNAs in Thyroid Carcinoma" Genes 9, no. 1: 44. https://doi.org/10.3390/genes9010044

APA Style

Lu, M., Xu, X., Xi, B., Dai, Q., Li, C., Su, L., Zhou, X., Tang, M., Yao, Y., & Yang, J. (2018). Molecular Network-Based Identification of Competing Endogenous RNAs in Thyroid Carcinoma. Genes, 9(1), 44. https://doi.org/10.3390/genes9010044

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop