1.
Introduction
Papillary thyroid carcinoma (PTC) is the most common malignant thyroid tumor, accounting for 80–90% of thyroid cancers. PTC is often characterized to have an early metastasis to the lymph node of the neck [1,2]. Although the prognosis of most patients with PTC is good, special histological subtypes of papillary carcinoma, iodium-resistant PTC, thyroid invasion outside the membrane, vascular invasion, and distant metastasis are often considered risk factors for postoperative PTC recurrence, leading to relatively poor prognosis [3,4]. At present, the molecular mechanism of the occurrence and development of PTC has not been clarified; therefore, an increasing number of studies are being performed to identify a molecular marker that can predict the occurrence, metastasis, and recurrence of PTC.
Circular RNA (circRNA) is a type of non-coding RNA that was first discovered in viroids in plants in 1976. circRNAs are single-stranded and covalently closed RNAs that lack a free cap end and poly (A) tail [5,6] and are widely distributed in mammalian cells [7]. circRNAs are stabilized in cells [8] and act as miRNA sponges to competitively bind to mRNA [9]. In recent years, the correlation between circRNAs and the occurrence and development of malignant tumors has become a focus of academic research. For example, the correlation between circRNA and the occurrence and development of bladder cancer [10], correlation between circRNA and the early diagnosis of breast adenocarcinoma [11], correlation between circRNA and the occurrence and development of colon cancer [12], and correlation between circRNA and the proliferation of liver cancer cells [13] have been discovered. Hence, the relationship between circRNAs and human malignancies has been evaluated and established; however, the diagnostic value and biological function of circRNAs in PTC have not been extensively studied.
In this study, microarray datasets for PTC and normal samples were downloaded from a public database. Differentially expressed mRNAs (dif-mRNAs), differentially expressed miRNAs (dif-miRNAs), and differentially expressed circRNAs (dif-circRNAs) between PTC and control tissues were investigated, followed by an enrichment analysis of the dif-mRNAs. To further explore the mechanism of key mRNAs, circRNAs, and miRNAs in PTC, protein–protein interaction (PPI) network and regulatory network analyses were performed. Therefore, this study provides novel therapeutic targets for PTC.
2.
Materials and methods
2.1. Data sources
The microarray datasets GSE93522 (platform: Agilent-069978 Arraystar Human CircRNA microarray V1), GSE73182 (platform: Agilent-035758 Human miRBASE 16.0 plus 031181), and GSE66783 (GPL19850 Agilent-060228 Human LncRNA v5 4X180K [Probe Name Version]) from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database were downloaded. The dataset GSE93522 includes circRNA expression data of 6 PTC and 6 paired adjacent noncancerous thyroid tissues, GSE73182 includes miRNA expression data of 19 primary PTC and 5 normal thyroids, and GSE66783 includes mRNA expression data of 5 PTC and 5 paired adjacent noncancerous thyroid tissues.
In addition, three microarray datasets were downloaded for external validation, including GSE173299 (circRNA expression data of three PTC and three normal samples), GSE104006 (miRNAs and mRNA expression data of 24 PTC and 5 normal samples), and GSE151181 (miRNA and mRNA expression data of 39 PTC and 13 normal samples).
2.2. Differential expression analysis
The original data were read using the limma package (Version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) [14] and preprocessed using the robust multi-array average (RMA). Data preprocessing included background correction, normalization, and expression calculations. Gene expression values were calculated by mapping probes to symbols using the microarray dataset and annotation files from the chip platform. The probes were removed when no genes were matched. The mean value was selected as the expression value when multiple probes matched one gene symbol. The dif-mRNAs, dif-miRNAs, and dif-circRNAs between samples from patients with PTC and control tissues were analyzed using a typical Bayesian method. The dif-miRNAs and dif-circRNAs were defined as miRNAs or circRNAs with |log fold change (FC)| greater than 1 and a P-value less than 0.05. The dif-mRNAs were defined as mRNAs with |log2 (FC)| greater than 2 and a P-value less than 0.05.
2.3. Pathway enrichment analysis of dif-mRNAs
The biological process (BP) terms of the Gene Ontology (GO) [15] annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [16] pathways were enriched for the upregulated and downregulated overlapping dif-mRNAs using the enrichment analysis tool DAVID. The number of enriched genes was set as count ≥ 2. Statistical significance was set at P < 0.05. The REViGO online tool was utilized to remove redundancy among the GO-BP terms, with the following parameter settings: allowed similarity, small (0.5); species, Homo sapiens; and semantic similarity measure, Lin.
2.4. PPI network and module analyses for dif-mRNAs
Using the STRING database (version: 10.0, http://www.string-db.org/) [17], the interactions between proteins encoded by dif-mRNAs were retrieved. The gene set consisted of all differentially expressed genes, and humans were used as the species. A PPI score of 0.7 was defined as having high confidence. Cytoscape software [18] (version: 3.2.0, http://www.cytoscape.org) was used to construct the PPI network based on the obtained interactions. Significant network modules were analyzed using the Molecular Complex Detection (MCODE) plug-in (version: 1.4.2, http://apps.cytoscape.org/apps/MCODE) [19] in the Cytoscape software. Modules with scores ≥ 10 were considered significant. The BP terms of the GO annotation and KEGG pathways were enriched for the genes in the significant clustering modules. The number of enriched genes was set as count ≥ 2. Statistical significance was set at P < 0.05.
2.5. Screening of key dif-circRNAs
The support vector machine (SVM)-recursive feature elimination (RFE) algorithm was used to screen the key dif-circRNAs. In brief, five-fold cross-validation was conducted for all the dif-circRNAs using the Caret package (version 6.0-86). Thereafter, SVM-RFE was performed for dimensionality reduction screening of dif-circRNAs, and the circRNAs contained in the model with the lowest RMSE (standard error) were selected as the key dif-circRNAs.
2.6. Construction of a competing endogenous RNA (ceRNA) network
The circRNA-miRNA interaction pairs and miRNA-mRNA interaction pairs were predicted for all dif-miRNAs by utilizing the Encyclopedia of RNA Interactomes (ENCORI, previous version: starBase v2.0, http://starbase.sysu.edu.cn/) with the following parameters: CLIP-Data ≥ 1, pan-Cancer ≥ 0 and Degradome-Data ≥ 0. The circRNA–miRNA and miRNA–mRNA interaction pairs were first screened using the key dif-circRNAs (obtained from SVM-RFE analysis) and mRNAs in PPI modules; subsequently, the circRNA-miRNA and miRNA-mRNA interaction pairs that involved the same miRNAs were retained. Based on the ceRNA hypothesis, two RNA molecules regulated by the same miRNA are co-expressed, and an miRNA usually inhibits the expression of its target gene (mRNA). Therefore, positive circRNA-mRNA interaction pairs and negative miRNA-mRNA interaction pairs were selected to construct ceRNA regulatory pairs.
2.7. Associations of key mRNAs with survival and clinical features
The expression data and clinical phenotypes of thyroid cancer were downloaded from TCGA database, and only the TPC samples with clinical phenotypes were selected. A total of 336 samples were obtained. The expression of key mRNAs (mRNAs in the ceRNA network) was extracted, and samples were assigned into high- and low-expression groups based on their median expression value. Kaplan-Meier survival analysis was performed to screen for prognostic mRNAs. In addition, associations between the expression of key mRNAs with clinical phenotype were analyzed, including age, sex, pathologic T/N/M stages, primary thyroid gland neoplasm location anatomic site, and tumor stage. The t-test was used for comparison between groups.
2.8. Validation of the expression of key RNAs
The expression of key circRNAs, miRNAs, and mRNAs between PTC and normal tissues was validated using external validation datasets and the UALCAN online tool. The t-test was used for comparison between groups.
3.
Results
3.1. Differential expression analysis
According to the screening criteria, there were 137 dif-circRNAs (115 upregulated and 22 downregulated), 18 dif-miRNAs (9 upregulated and 9 downregulated), and 1039 dif-mRNAs (454 upregulated and 585 downregulated) between PTC and normal tissues (Figure 1A, B). Principal component analysis showed that samples in the different groups could be clearly distinguished (Figure 1C).
3.2. Enrichment analysis of dif-mRNAs
The functions of dif-mRNAs were explored using enrichment analysis. The upregulated and downregulated mRNAs were enriched in 81 and 114 GO-BP terms, respectively. The significantly enriched GO-BP terms obtained after removing redundancy among the GO-BP terms using REViGO are shown in Figures 2A, B. It could be seen that the upregulated mRNAs were mainly involved in extracellular matrix disassembly, collagen catabolic process, leukocyte migration, regulation of cell adhesion mediated by integrin, and others (Figure 2A). The downregulated mRNAs were mainly implicated in negative regulation of growth, cell adhesion, cellular response to cadmium ions, among others (Figure 2B). Additionally, the significantly enriched KEGG pathways are shown in Figure 2C, D. The upregulated mRNAs were mainly involved in transcriptional misregulation in cancer, p53 signaling pathway, and cytokine−cytokine receptor interaction (Figure 2C) while the downregulated mRNAs were mainly involved in tyrosine/retinol metabolism, drug metabolism-cytochrome P450, TGF−beta signaling pathway, and the Hippo signaling pathway (Figure 2D).
3.3. PPI network and module analyses
The PPI network of dif-mRNAs comprised 361 nodes and 805 interactions (Figure 3A). FN1 (degree = 28), APOA1 (degree = 23), SERPINA1 (degree = 23), TIMP1 (degree = 22), APOE (degree = 22), GPC3 (degree = 21), APOA2 (degree = 20), NMU (degree = 20), and LPAR5 (degree = 20) were hub nodes with the highest degrees in the network. Three significant modules were identified separately from the PPI network (Figure 3B). Module A had 15 nodes and 105 interactions, of which FN1, APOA1, SERPINA1, APOE, TIMP1, GPC3 and APOA2 were also hub nodes. The mRNAs in module A were mainly enriched in GO-BP terms such as high-density lipoprotein particle assembly/clearance, positive regulation of cholesterol esterification, and lipoprotein biosynthetic process and KEGG pathways such as ECM–receptor interaction and focal adhesion. Module B had 21 nodes and 122 interactions, and the mRNAs in this module were mainly involved in GO-BP terms such as G-protein coupled receptor signaling pathway, positive regulation of ERK1 and ERK2 cascade, and cell–cell signaling and in KEGG pathways such as calcium signaling pathway and cytokine–cytokine receptor interaction. Module C had 10 nodes and 45 interactions, and the mRNAs in this module were mainly implicated in 4 GO-BP terms, namely negative regulation of viral genome replication, antibacterial humoral response, response to mechanical stimulus, and cellular response to interleukin-1 (Table 1).
3.4. Key dif-circRNAs screened via SVM-RFE
The expression of 137 dif-circRNAs was extracted, and SVM-RFE was utilized for the dimensionality reduction screening of dif-circRNAs. As described in the Methods section, the circRNAs with the lowest RMSE (standard error) contained in the model were selected as the key dif-circRNAs, and consequently, 37 dif-circRNAs were obtained (Figure S1, Table S1).
3.5. Construction of a ceRNA regulatory network
A total of 11 circRNA-miRNA interaction pairs were predicted between the key circRNAs and dif-miRNAs. Overall, 40 miRNA-mRNA interaction pairs were predicted between the mRNAs in modules and dif-miRNAs. Finally, 46 circRNA-miRNA-mRNA regulatory networks were constructed, which contained 7 circRNAs (6 upregulated and 1 downregulated), 11 miRNAs (6 upregulated and 5 downregulated), and 18 mRNAs (12 upregulated and 6 downregulated) (Figure 4A).
Furthermore, the ceRNA regulatory pairs were screened. Based on the ceRNA hypothesis, two RNA molecules regulated by the same miRNA are usually co-expressed, and miRNAs usually inhibit the expression of its target gene (mRNA). Therefore, positive circRNA-mRNA interaction pairs and negative miRNA-mRNA interaction pairs were selected to construct ceRNA regulatory pairs. A total of eight ceRNA regulatory pairs were obtained, including two circRNAs (circ_0004053 and circ_0028198), three miRNAs (miR-199a-5p, miR-199b-5p, and miR-7-5p), and five mRNAs, including APOA2, CCL20, LPAR5, MFGE8 and TIMP1 (Figure 4B). To investigate the correlation between miRNA and mRNA expression, two external datasets (GSE104006 and GSE151181) containing both miRNA and mRNA expression data were used. There were significant negative correlations between the expression of miR-7-5p and its targets LPAR5, MFGE8, and TIMP1 in both the GSE104006 and GSE151181 datasets (Figure S2, Table S2).
3.6. Association between the expression of key mRNAs with survival and clinical features
The five mRNAs in the ceRNA regulatory pairs were considered key mRNAs. We further investigated the association of these five key mRNAs with survival and clinical features (Figures 5 and S3). Survival analysis showed that only LPAR5 expression was associated with patient survival; in particular, high LPAR5 expression was associated with favorable patient survival (Figure 5A). In addition, APOA2 expression showed significant differences between metastatic and non-metastatic tumors, while the remaining four mRNAs, CCL20, LPAR5, MFGE8 and TIMP1, showed significant differences between metastatic and non-metastatic lymph nodes (Figure 5B).
3.7. Validation of the expression of key mRNAs
The expression of key circRNAs (circ_0004053 and circ_0028198), miRNAs (miR-199a-5p, miR-199b-5p, and miR-7-5p) and mRNAs (APOA2, CCL20, LPAR5, MFGE8 and TIMP1) in PTC and normal tissues was validated using external validation datasets. As shown in Figure 6A and 6B, the expression of the three miRNAs was confirmed to be downregulated in PTC compared to normal samples, and all the five mRNAs were upregulated in PTC compared to normal samples. The expression of miR-199a-5p in the GSE151181 dataset and expression of APOA2 in the GSE104006 dataset showed no significant differences. The expression of circ_0004053 and circ_0028198 was upregulated in PTC compared to that in normal samples (Figure 6C). In addition, the expression of key miRNAs and mRNAs in PTC and normal tissues was validated using the UALCAN online tool. As shown in Figure 7, the expression of the three miRNAs was downregulated in tumor samples compared to that in normal samples, especially in classical thyroid papillary carcinoma and follicular thyroid papillary carcinoma (Figure 7A). The five mRNAs showed significantly higher expression in tumor samples than in normal samples; in particular, the expression of these five mRNAs was significantly higher in classical thyroid papillary carcinoma and tall thyroid papillary carcinoma (Figure 7B).
Box plots showing the expression of key miRNAs (upper panel) and mRNAs (lower panel) between PTC and normal samples using the external GSE151181 dataset (A) and the GSE104006 dataset (B). Box plots showing the expression of circRNAs in PTC and normal samples using the external GSE173299 dataset (C). The lines in the middle of each box represent the median expression values of the corresponding mRNA. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.
Box plots showing the expression differences of key miRNAs (A) and mRNAs (B) between normal and different types of tumor tissue samples, including classical thyroid papillary carcinoma, follicular thyroid papillary carcinoma, Tall thyroid papillary carcinoma and others. Blue blocks represent normal samples, and the other blocks represent tumor histology, including classical, follicular, and tall thyroid papillary carcinoma. The black line in each box represents the median expression value of the corresponding RNA.
3.
Discussion
circRNAs are non-coding RNAs usually expressed in a specific tissue or developmental stage and have been considered promising alternative biomarkers owing to their high stability imparted by their covalently closed cyclic structures [20,21]. Various circRNAs have been identified in tumor tissues, functioning as ceRNAs by sponging miRNAs and regulating the transcription and translation of target genes. Through these mechanisms, circRNAs significantly influence several diseases and malignant tumors [22,23,24]. In this study, 137 circRNAs were found to be differentially expressed in PTC tissues, of which 37 key circRNAs were obtained using the SVM-RFE algorithm. In addition, 1039 dif-mRNAs were screened in PTC, of which 46 key mRNAs in PPI modules were identified. Among the 37 key dif-circRNAs, 46 key mRNAs, and 18 dif-miRNAs, 46 circRNA-miRNA-mRNA regulatory pairs were predicted. Finally, eight regulatory pairs in line with the ceRNA hypothesis were obtained, including two circRNAs (circ_0004053 and circ_0028198), three miRNAs (miR-199a-5p, miR-199b-5p, and miR-7-5p), and five mRNAs (APOA2, CCL20, LPAR5, MFGE8, and TIMP1).
circ_0028198 has been reported to promote the pathophysiology of atherosclerosis and lipid disorders through regulation of lipoprotein metabolism and cholesterol homeostasis by decreasing cholesterol efflux from macrophages [25,26,27]; however, the role of circ_0028198 in the development and progression of tumors has not yet been demonstrated. circ_0004053 was identified as a novel circRNA. In our study, circ_0004053 was found to regulate the expression of APOA2 and CCL20 by targeting miR-119a-5p and miR-119b-5p, while circ_0028198 targeted miR-7-5p to regulate the expression of APOA2, LPAR5, MFGE8, and TIMP1. miR-119a-5p and miR-119b-5p are all members of the miR-199 family, which has been found to be involved in different cancers. For example, Zhang et al. suggested that miR-199 could repress epithelial–mesenchymal transition (EMT) and the invasion of tumor cells by inhibiting the expression of Snail in hepatocellular carcinoma [28]. Koshizuka et al. revealed that the miR-199 family showed antitumor activity in head and neck squamous cell carcinoma, inhibiting the migration of tumor cells by decreasing the expression of integrin α3 [29]. miR-7 has been reported to regulate the expression of multiple oncogenes [30]. For example, Xiao et al. showed that miR-7-5p could inhibit tumor growth and metastasis by targeting neuro-oncological ventral antigen 2 in non-small cell lung cancer both in vitro and in vivo [31]. In addition, miR-7 has been reported to be involved in colorectal cancer growth and metastasis [32], reversing chemotherapy resistance in breast cancer [33] and inhibiting tumor progression in pancreatic cancer [34] and esophageal cancer [35]. miR-7-5p may be a diagnostic biomarker of PTC [36]; additionally, it has been reported to be involved in regulating the proliferation and migration of PTC cells [37]. These studies suggest the important role of miR-199 and miR-7-5p in tumors; moreover, circ_0004053 and circ_0028198 target these miRNAs and thus may play key roles in PTC.
The ceRNA regulatory network analysis suggested that circ_0004053 and circ_0028198 regulate the expression of APOA2, CCL20, LPAR5, MFGE8, and TIMP1. Further analysis revealed that the expression of these five mRNAs showed significant differences between metastasis/lymph node metastasis and non-metastasis/non-lymph node metastasis. In addition, survival analysis showed that LPAR5 expression was associated with patient survival. The APOA2 gene encodes apolipoprotein A2, which has been suggested as a biomarker for pancreatic cancers [38,39]. CCL20 encodes a C-C motif chemokine ligand for the receptor CCR6; thus, it promotes cell migration and invasion in thyroid cancer by activating CCR6 [40]. LPAR5 encodes a lysophosphatidic acid receptor, which contributes to the metastasis and tumorigenesis of PTC by activating the PI3K/AKT pathway and EMT [41,42]. MFGE8 is a milk fat globule-epidermal growth factor that has been found to play crucial roles in regulating transforming growth factor-β-induced EMT in endometrial epithelial cells [43] and cardiac fibrosis [44]. Tissue inhibitors of metalloproteinases (TIMPs) have been considered to contribute to cancer hallmarks. TIMPs are always aberrantly expressed in tumors, and a high expression of TIMP1 is correlated with tumor progression and worse survival [45]. Zhang et al. suggested that TIMP1 was implicated in the lymphatic metastasis of PTC by regulating the invasive and migratory activities of PTC cells [46]. These data indicate that the five aforementioned mRNAs play significant roles in mediating EMT, metastasis, cell invasion, and migration, consistent with the findings that the expression of these five mRNAs exhibit significant differences between metastasis/lymph node metastasis and non-metastasis/non-lymph node metastasis.
Although we have validated the expression of key circRNAs, miRNAs, and mRNAs using external datasets, their expression should be further confirmed through biological experiments using clinical samples, and their clinical value should be further investigated based on clinical data with a large sample size. In addition, the ceRNA regulatory pairs predicted were based on positive expression trends between circRNAs and mRNAs and negative expression trends between miRNAs and mRNAs because expression correlation analysis could not be performed among the different datasets. We performed Pearson correlation analyses using external datasets to investigate the expression correlation between miRNAs and mRNAs, and no significant correlations were found between miR-199a/b-5p and their targets. Therefore, in the future, a dual luciferase assay should be used to verify the binding of an miRNA to its target, and the ceRNA regulatory pairs should be confirmed by conducting functional experiments.
4.
Conclusions
In conclusion, we identified eight ceRNA regulatory pairs, involving two circRNAs (circ_0004053 and circ_0028198), three miRNAs (miR-119a-5p, miR-119b-5p, and miR-7-5p), and five mRNAs (APOA2, CCL20, LPAR5, MFGE8, and TIMP1), in PTC. These five mRNAs were correlated with metastatic and non-metastatic lymph nodes in PTC. Our study reveals several potential targets and the underlying regulatory mechanisms in PTC.
Acknowledgments
This work was supported by the Hainan Provincial Department of Science and Technology, under Grant No. HN01ENT07, and the Clinical Medicine Project of the Affiliated Hospital of Qingdao University, under Grant No. 3400.
Conflict of interest
The authors declare that they have no competing interests.
DownLoad:
