Bioinformatic analysis and identification of potential prognostic microRNAs and mRNAs in thyroid cancer

Introduction

Thyroid cancer is a common endocrine malignancy which has increased rapidly worldwide in the past decades. The analysis of the Surveillance, Epidemiology, and End Result cancer registry data between 1975 and 2013 revealed that the incidence of thyroid cancer increased by 211% (Lim et al., 2017). The improvements in diagnosis and treatment for thyroid cancer substantially improve prognosis (Grigsby et al., 2006; Lang, Wong & Wan, 2013). According to the histopathological features, thyroid cancer is classified into three major categories: well-differentiated thyroid cancer (WDTC) which includes papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC); poorly-differentiated thyroid cancer (PDTC) and anaplastic thyroid cancer (ATC) (Zarkesh et al., 2017). While WDTC accounts for the majority part of all cases (Yu et al., 2013), and its mortality rates were controlled to less than 10% for 10 years, the mortality rates of PDTC and ATC were reported to be 38–57% and close to 100% respectively (Smallridge et al., 2012; Xu & Ghossein, 2016).

It has been widely accepted that different expression levels of specific genes are associated with cancer initiation. During the past few years, molecular cytogenetics studies have been used to investigate the molecular mechanisms of thyroid cancer. Multiple genes and cellular pathways were reported to participate in the occurrence and development of thyroid cancer. Alterations in the Ras-Raf-mitogen-activated protein kinase (MAPK) pathway are usually observed in WDTC (Hou et al., 2007). BRAF mutations were found in 32.4% of PTC cases and RAS mutations in 20–50% of FTC and 13 % of PTC cases (Nikiforova et al., 2003; Sahpaz et al., 2015). In the ATC, BRAF and RAS alterations were detected in 29% and 23% cases (Xu & Ghossein, 2016). Telomerase reverse transcriptase (TERT) promoter mutations were detected in 10% of PTC and related to clinically aggressive behaviors. In addition, mutations of both TERT promoter and BRAF/RAS have a tendency for co-occurrence (Landa et al., 2016; Song et al., 2016). TP53 mutation is considered to be a genetic event distinguishing WDTC from ATC. Several studies revealed that TP53 mutations were detected in 59% of ATC cases, compared with in 10% WDTC cases (Kunstman et al., 2015; Sykorova et al., 2015).

With the development of high-throughput technologies, expression profiling of multiple genes is a useful way to find different expression levels of specific genes between normal and tumor tissues. Altered expression levels usually indicate pathological conditions, and proteins coded by these differentially expressed genes may involve in different molecular pathways, biological process, and cellular behaviors during tumor progression. MicroRNAs are small noncoding RNAs which participate in the post-transcriptional regulation of gene expression (Chen & Kang, 2015). They function as negative regulators by binding to the 3′-untranslated region of candidate mRNAs, and repress the gene expression by inhibiting protein translation or degrading mRNAs (Wang et al., 2010). Accumulating evidence demonstrated that microRNAs could function as either oncogenes or tumor suppressors in various types of malignancies and regulate different carcinogenic processes (Liang, Li & Wang, 2017). In previous bioinformatic studies, a number of microRNAs and mRNAs were identified as predictors for the prognosis of thyroid cancer (Table S1). Due to the small sample sizes and different detection platforms, some results were controversial. In the current study, we analyzed the public microRNA and mRNA sequencing data from The Cancer Genome Atlas Project (TCGA, https://cancergenome.nih.gov/) and identified the differentially expressed microRNAs and mRNAs between thyroid cancer tissues and matched normal thyroid tissues. We then investigated the prognostic value of these differentially expressed microRNAs. In addition, combined with microRNA-mRNA interaction analysis, we analyzed the functions and pathways and constructed protein-protein interaction (PPI) network of differentially expressed mRNAs to investigate the underlying mechanisms of thyroid cancer occurrence and development.

Materials and Methods

Data processing

We obtained the clinical information of thyroid and RNA sequencing data from TCGA database (https://cancergenome.nih.gov/). A workflow of the data selection was shown in Fig. 1. Using edgeR package in R language, we normalized the microRNA and mRNA sequencing data and screened out differentially expressed microRNAs and mRNAs between thyroid cancer and normal tissues (Robinson, McCarthy & Smyth, 2010). Cut-off criteria of adjusted p-value < 0.01 and |log2 fold-change (FC)| > 2 were considered to be statistically significant.

Clinical information of thyroid cancer patients included sex (female and male), age at diagnosis (<45 or ≥45 years), race (White, Black, Asian and other), histological type (classical papillary thyroid cancer, follicular variant of papillary thyroid cancer, tall cell papillary thyroid cancer and other) and tumor-node-metastasis (TNM) stage. Using Kaplan–Meier survival curves by log-rank test, we evaluated the prognostic value of each differentially expressed microRNA with a threshold of p-value < 0.05. The microRNAs significantly associated with the overall survival were considered as prognostic microRNAs. A Cox proportion hazards model was used to evaluate the relative risk of these prognostic microRNAs on OS. Hazard ratios (HR) with 95% confidence intervals were obtained, any HR > 1.0 showed an increased risk of death. A P value < 0.05 was considered statistically significant and all tests were two-sided.

Results

Identification of differentially expressed microRNAs and mRNAs in thyroid cancer

A total of 573 samples with microRNA sequencing data and 568 samples with mRNA sequencing data were analyzed. The microRNA sequencing data included 514 thyroid cancer samples and 59 matched normal samples; and the RNA sequencing data included 510 thyroid cancer samples and 58 matched normal samples. Based the cut-off criteria (P < 0.01 and |log2FC| > 2.0), a total of 72 differentially expressed microRNAs and 1766 differentially expressed mRNAs were identified between thyroid cancer tissues and matched normal tissues, including 67 up-regulated and five down-regulated microRNAs; 1,370 up-regulated and 396 down-regulated genes. The results were presented as Volcano plot (Fig. 2).

Association between differentially expressed microRNAs and clinical features

To identify the potential prognostic microRNAs correlated with the overall survival of thyroid cancer patients, we investigated the associations between miRNAs expression and patient survival using Kaplan–Meier curve and Log-rank test. Seven microRNAs were related to the overall survival. While five of them (miR-146b, miR-184, miR-767, miR-6730 and miR-6860) were positively correlated with the overall survival, the other two microRNAs (miR-196a-2 and miR-509-3) were negatively related to the overall survival (Fig. 3). In the multivariate Cox analysis, only miR-184, miR-146b and miR-509-3 were significantly associated with overall survival when controlling for the remaining variables: age at diagnosis, sex, stage, tumor size, lymph nodes metastases, distant metastases and other diagnostic genes revealed by Kaplan–Meier curve (Table 1). Interestingly, all of the seven microRNAs were up-regulated in the tumor samples. MiR-146b, miR-509-3 and miR-6730 were significantly associated with tumor stage and lymph node metastases; only miR-509-3 was related to tumor size (Fig. 4).

Table 1:

Cox proportional hazards regression model analysis of factors associated with overall survival.

Variable	Multivariate hazard ration	95% CI	*P-value
mir-767 expression High vs. Low	0.316	0.071–1.412	0.131
mir-6730 expression High vs. Low	0.370	0.081–1.685	0.199
mir-184 expression High vs. Low	0.201	0.045–0.901	0.036
mir-6860 expression High vs. Low	0.729	0.142–2.749	0.706
mir-146b expression High vs. Low	0.251	0.065–0.970	0.035
mir-196a-2 expression High vs. Low	2.864	0.065–4.881	0.147
mir-509-3 expression High vs. Low	4.534	2.281–9.923	0.002
LPAR5 expression High vs. Low	0.043	0.012–0.453	0.005

DOI: 10.7717/peerj.4674/table-1

Notes:

*P values calculated by multivariate Cox analysis (adjusted for age at diagnosis, sex, stage, tumor size, lymph nodes metastases and distant metastases).

CL

confidence interval

Gene set enrichment analysis

Based on the cut-off criteria FDR < 0.05, only three functional gene sets were enriched: “Notch signaling pathway” “Glycosaminoglycan degradation” and “P53 signaling pathway” (Fig. 5).

Protein-protein interaction and mircoRNA-target network analysis

To predict protein interactions, the 1,766 differentially expressed genes were submitted to the Search Tool for the Retrieval of Interacting Gene (STRING) database. The PPI network consisted of 483 nodes and 1,437 edges. The target genes of prognostic microRNAs were compared with these differentially expressed genes, and six target genes regulated by three microRNAs were identified in the network (miR-767 was predicted to target COL10A1, PLAG1 and PPP1R1C; miR-146b was predicted to target MMP16; miR-196a-2 was predicted to target SYT9) (Fig. S1). Moreover, the network complex was further analyzed, and the most significant module was screened out using MCODE, which contained 26 nodes and 325 edges (Fig. 6). Then we performed KEGG pathway analysis of these 26 genes, our results demonstrated that they were mainly involved in neuroactive ligand–receptor interaction, chemokine signaling pathway and cAMP signaling pathway. We screened out the top 10 hub nodes with higher degrees using the plug-in CytoHubba in Cytoscape. These hub genes included neuropeptide Y (NPY), neuromedin U (NMU), kininogen 1 (KNG1), lysophosphatidic acid receptor 5 (LPAR5), C–C motif chemokine receptor 3 (CCR3), somatostatin (SST), pancreatic polypeptide (PPY), gamma-aminobutyric acid type B receptor subunit 2 (GABBR2), adenylate cyclase 8 (ADCY8) and serum amyloid A1 (SAA1). Survival analysis of the 10 hub genes demonstrated that only LPAR5 was significantly correlated with the overall survival (Fig. 3H). Multivariate Cox analysis demonstrated that LPAR5 was an independent risk factors for overall survival (Table 1). And the expression of LPAR5 was associated with tumor stage, size and lymph node metastases (Fig. 4).

Discussion

As one of the most common endocrine malignancies, it is important to investigate the molecular mechanisms of thyroid cancer occurrence and development. In the current study, we investigated potential microRNAs and mRNAs correlated with thyroid tumorigenesis using bioinformatics analysis. A total of 72 differentially expressed microRNAs and 1,766 differentially expressed genes were identified from TCGA database, including 67 up-regulated and five down-regulated microRNAs;1370 up-regulated and 396 down-regulated genes. MiR-146b, miR-184, miR-767, miR-6730, miR-6860, miR-196a-2, miR-509-3 and LPAR5 were correlated with the overall survival of thyroid cancer patients. GSEA analysis demonstrated that the gene sets “Notch signaling pathway” “Glycosaminoglycan degradation” and “P53 signaling pathway” were significantly enriched in thyroid cancer samples.

In the last decade, microRNAs have been revealed to modulate multiple processes of cancer development, including cancer cell proliferation, differentiation, apoptosis, migration and invasion. However, the studies of microRNAs’ effects on cancer prognosis was limited due to the small sample size, different detection platforms and the lack of long-term follow-up. To investigate the potential prognostic microRNAs of thyroid cancer, we analyzed the high-throughput data from TCGA database and identified seven prognostic microRNAs associated with the clinical outcome of thyroid cancer patients. Geraldo et al. stated that miR-146b-5p was highly expressed in papillary thyroid cancer and considered as a relevant diagnostic marker for this type of thyroid cancer. MiR-146b-5p was reported to promote the migration and invasion of papillary thyroid cancer cells via downregulating ZNRF3 and upregulating Wnt/catenin signaling pathway, and to promote thyroid follicular cell growth via downregulating TGF- β pathway by binding to the 3′-untranslated region of SMAD4 (Chou, Liu & Kang, 2017; Deng et al., 2015; Geraldo, Yamashita & Kimura, 2012). Inversely, the expression levels of miR-146b was positively related to the overall survival of patients with thyroid cancer. The expression levels of miR-509-3 were lower in multiple cancers, and miR-509-3-5P downregulation promoted the migration and invasion of gastric cancer cells by targeting PODXL. Overexpression of miR-509-5p markedly inhibited the proliferation, migration and invasion of non-small lung cancer cells via targeting YWHAG (Peng, Yu & Fu, 2016; Zhang et al., 2017). While miR-509-3 was overexpressed in thyroid cancer and negatively related to the overall survival, miR-184 was considered as a potential oncogenic microRNA of squamous cell carcinoma via promoting cancer cell proliferation. Inversely, miR-184 was down-regulated in aggressive human glioma and breast cancer cells, and inhibited cancer cell proliferation and invasion (Emdad et al., 2015; Feng & Dong, 2015). In the present study, miR-184 was positively related to the overall survival of patients with thyroid cancer. MiR-767 was reported to represses TET1/3 (two tumor suppressor genes) and identified as a hallmark of cancer (Loriot et al., 2014). The molecular mechanisms are still to be investigated. Previous studies have identified a number of miRNAs involved in thyroid carcinogenesis (Table S2). MiR-146b, miR-221 and miR-222 were the most three frequent miRNAs reported in thyroid cancer, they appeared to associated with high-risk features such as extrathyroidal extension, lymph node metastasis, distant metastasis and BRAF^V600E mutation. While in our present study, the results were partly different from the previous: miR-221, miR-222 and some other miRNAs were not identified. This might because of the following reasons: different study types (RT-PCR, array, high throughput sequencing), different platforms, different ways for normalization, different ways for DEG screening, and the cut-off criteria (in our study P < 0.01 and |log2FC| > 2.0). In our validation datasets GSE3467 and GSE73182, 130 differentially expressed mRNAs and five differentially expressed miRNAs were identified, among them, 95 differentially expressed mRNAs and 4 differentially expressed miRNAs were overlapping with TCGA dataset.

Through integrated bioinformatics analysis, we identified the most significant module which contained 26 nodes and 325 edges, and these genes were mainly enriched in neuroactive ligand–receptor interaction, chemokine signaling pathway and cAMP signaling pathway. The top 10 hub genes with higher degrees were NPY, NMU, KNG1, LPAR5, CCR3, SST, PPY, GABBR2, ADCY8 and SAA1. NPY encodes a neuropeptide which influences multiple physiological processes through G protein-coupled receptors (GPCRs) and MAPK. NPY was reported to promote inflammation-induced tumorigenesis via PI3-K/Akt pathway and miR-375-dependent apoptosis (Jeppsson, Srinivasan & Chandrasekharan, 2016). NMU was associated with increased breast cancer aggression (Martinez et al., 2017) and its overexpression induced regional metastasis of head and neck squamous cell carcinoma (Wang et al., 2016). KNG1 was identified as a potential marker of early colorectal cancer stages (Quesadacalvo et al., 2017). LPAR5 is a member of the rhodopsin class GPCRs. In our study, it was positively corrected with the overall survival of thyroid cancer patients. LPAR5 was down-regulated in primary undifferentiated nasopharyngeal carcinoma and promoted the LPA-induced migration of nasopharyngeal carcinoma cell lines (Yap et al., 2015). Additionally, LPAR5 negatively regulated cell motile and invasive activities of human sarcoma cell lines (Dong et al., 2014). CCR3 belongs to family 1 of the GPCRs, which enhance cellular proliferation, invasion, and migration through ERK and JNK signaling pathway (Dong et al., 2014). CCRs was reported to correlated with improved distant relapse-free survival in breast cancer (Gong et al., 2016). SST is a neuropeptide which affects neurotransmission, secretion and cell proliferation. The receptor of SST was reported to be a predictor of better response to therapy in medullary thyroid carcinoma (Kendler et al., 2017). In radiation-induced papillary thyroid cancer from chernobyl pediatric patients, GABBR2 was highly expressed (Stein et al., 2010), and related to tumor occurrence. ADCY8 is a membrane bound enzyme which is differentially expressed in endometrial cancer (Orchel et al., 2012). SAA1 is highly expressed in response to tissue injury and inflammation, and its highly expression is associated with chronic inflammation, lipid metabolism and tumor pathogenesis (Sun & Ye, 2016). SAA1 has been used as a non-invasive biomarker for the prognosis of many cancers, including stomach, breast, liver and lung neoplasms (Knebel et al., 2017; Upur et al., 2015).

Conclusion

In conclusion, our study identified the crucial microRNAs and mRNAs in thyroid cancer, and constructed the regulatory network between microRNAs and mRNAs through bioinformatics analysis. A total of 72 differentially expressed microRNAs and 1,766 differentially expressed genes were screened out. Among them, seven microRNAs were correlated with the overall survival. Among the hub genes identified from PPI network, LPAR5 may play important roles in thyroid cancer. Multivariate analysis demonstrated that miR-184, miR-146b, miR-509-3 and LPAR5 were an independent risk factors for prognosis, and they have the potential to be the targets for treatment of thyroid cancer. However, further experimental research is still required to confirm the functions of identified genes.

Supplemental Information