Abstract Background Papillary thyroid carcinoma (PTC) is the most common endocrine malignancy and its incidence has increased over the last few decades. The molecular mechanisms underlying PTC tumorigenesis and progression are still unclear. Patients and methods The microRNA (miRNA) expression patterns of PTC were revealed by miRNA microarray analysis and validated with The Cancer Genome Atlas (TCGA) data. Promoter DNA methylation rates of miR-204 were analyzed by Agena Methylation MassAR-RAY analysis and validated with TCGA data. The underlying molecular mechanisms of miR-204 involved in PTC were studied by bioinformatics analyses. Results A total of 181 differentially expressed miRNAs (89 downregulated and 92 upregulated miRNAs) between PTC and normal tissues were detected in this study. We identified miR-204 as one of the most significantly downregulated miRNAs in PTC. Downregulation of miR-204 was related to PTC extrathyroidal extension, high T-stage, lymph metastasis, BRAF V600E mutation, and aggressive tall cell variant. The Agena MassARRAY results indicated that 12 CpG sites located at the promoter of miR-204 were hypermethylated in PTC tissues compared to normal tissues. The promoter methylation rates of miR-204 in PTC were negatively correlated with the expression levels of miR-204 and its host gene TRPM3. Downregulated miR-204 expression was related to several important pathways and mechanisms involved in tumorigenesis and progression. Conclusion Promoter DNA methylation-silenced miR-204 could serve as a potential diagnostic biomarker of PTC. Downregulation of miR-204 may play an important role in PTC via its involvement in many tumor-related pathways. Novel target genes and putative mechanisms of miR-204 in PTC need to be further validated. Keywords: papillary thyroid carcinoma, miRNA, promoter DNA methylation, bioinformatics Introduction Papillary thyroid carcinoma (PTC), a differentiated thyroid cancer, is the most common endocrine malignancy, and its incidence has increased over the last few decades. The prognosis of PTC is excellent. If treatment based on the current guidelines is given, cases of recurrence or persistence are rare. Surgery and radioiodine therapy followed by levothyroxine substitution remain the established therapeutic treatment.[33]^1 However, in patients with metastatic and/or radioiodine refractory tumors, additional treatment options are urgently needed. Over the last two decades, our understanding of the genetic mechanisms of the thyroid cancer has dramatically expanded. BRAF mutation is the first single gene mutational marker used for the diagnosis of PTC. Currently, several more advanced molecular tests are available for clinical use. These tests are based on multiple molecular markers and have an increasing impact on the clinical management of patients with thyroid nodules.[34]^2 The molecular mechanisms underlying PTC tumorigenesis and progression are still unclear. Therefore, a better understanding of these mechanisms will facilitate the development of novel diagnostic and even therapeutic targets for PTC. MicroRNAs (miRNAs) are noncoding RNAs that range in size from 19 to 25 nucleotides. They are important molecules with a vital role in several cellular functions, such as cell differentiation, cell cycle regulation, apoptosis, and angiogenesis in human diseases.[35]^3^,[36]^4 Typically, miRNAs basepair to target the 3′ UTR of the mRNAs within the RNA-induced silencing complex and block translation and stimulate mRNA transcript degradation.[37]^5 Numerous studies have reaffirmed the existence of mature miRNAs in the nucleus, and the nucleus-cytoplasm transport mechanism has also been illustrated. Moreover, active regulatory functions of nuclear miRNAs, including posttranscriptional gene silencing, transcriptional gene silencing, and transcriptional gene activation (TGA), in which miRNAs bind nascent RNA transcripts, gene promoter regions, or enhancer regions and exert further effects via epigenetic pathways, have been identified.[38]^6 MiRNAs have emerged as important regulators of human carcinogenesis by affecting the expression of key oncogenes and tumor suppressor genes.[39]^7 A new strategy of miRNA therapy has been developed in the last 5 years. There are two miRNA-based treatment approaches for cancer: miRNA reduction or inhibition therapy and miRNA replacement or restoration therapy.[40]^3 Accumulating evidence has revealed that miRNAs can serve as diagnostic and therapeutic targets for human cancers. Research on miRNA expression in thyroid tumors has provided new insights into the development of novel biomarkers that can be used to diagnose thyroid cancer and optimize its management.[41]^8 In the present study, we analyzed the expression patterns of miRNAs in five pairs of PTC vs normal tissues and compared our findings with the results derived from data deposited in The Cancer Genome Atlas (TCGA). We identified miR-204 as one of the most significantly downregulated miRNAs in PTC. The relationship between miR-204 and the clinicopathologic characteristics of PTC and the putative mechanism underlying the downregulation of miR-204 were explored. Promoter DNA methylation rates of miR-204 were analyzed by Agena Methylation Mas-sARRAY analysis and validated with TCGA data. Furthermore, we studied the underlying molecular mechanism of miR-204 involved in PTC by bioinformatics analyses. This study revealed that methylation-silenced miR-204 acts as an important biomarker in PTC and miR-204 may be involved in various pathway abnormalities in the tumorigenesis and progression of PTC. Patients and methods Patient samples Fifteen pairs of human PTC and matched adjacent normal tissues were collected from the Department of General Surgery of Xiangya Hospital of Central South University, People’s Republic of China. All tissues were immediately snap-frozen in liquid nitrogen. The cases selected were patients who had not received previous chemotherapy or radiotherapy before surgery and were based on clear pathological diagnosis. Among the 15 pairs of tissues used in our study, ten pairs of tissues were used for methylation array and quantitative real-time PCR (qRT-PCR) analyses, whereas others were used for miRNA expression analysis. This study was approved by the Institutional Ethics Committee of Xiangya Hospital. According to the Declaration of Helsinki, written informed consent permitting the use of their samples and publication of relevant data for the purpose of the present study was obtained from all patients. RNA extraction and miRNA microarrays Total RNA was extracted with TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) and purified using a mirVana™ miRNA Isolation Kit (AM1561). RNA dephosphorylation reactions were performed at 37°C for 30 minutes using a calf intestinal phosphatase master mix. Dephosphorylated RNA was labeled using the T4 RNA ligase labeling method with 4.5 µL ligation master mix (including 1.0 µL 10× T4 RNA ligase buffer, 3.0 µL Cyanine 3-pCP, and 0.5 µL T4 RNA ligase). The labeling reaction was performed at 16°C for 2 hours. The hybridization mixture was composed of 17.0 µL labeled miRNA sample, 1.0 µL Hyb Spike-In, 4.5 µL 10× gene expression (GE) (Agilent, Santa Clara, CA, USA) blocking agent, and 22.5 µL 2× Hi-RPM hybridization buffer. Hybridization was performed on a microarray chip (Agilent-070156 human miRNA V21.0 microarray), which contained the latest 2,549 miRNA probes spotted in triplicates. Hybridization reactions were performed at 20 rpm for 16 hours in an Agilent hybridization oven. The hybridized five pairs of chips were scanned with an Agilent G2565CA laser confocal scanner and visualized with Agilent feature extraction (v 10.7) software. Deregulated miRNAs were characterized by a fold change (FC) >2 and P-value <0.05. The microarray data have been deposited in the Gene Expression Omnibus (GEO) database (accession code [42]GSE113629). Cell culture and qRT-PCR Human PTC cell lines BCPAP and TPC-1 and the normal cell line Nthy-ori 3-1 were used in this study. Cells were cultured in RPMI 1640 supplemented with FBS (10%). Isolated total RNA from tissues and cells was reversed transcribed using the miRNA Reverse Transcription Kit (Thermo Fisher Scientific, Inc) according to the manufacturer’s protocol. cDNA was mixed with TaqMan MicroRNA Assay Kit and amplified using Applied Biosystems Instruments (ABI) 7900HT system. The primer sequences[43]^9 were as follows: U6 forward, 5′-GCTTCGGCAGCACATATACTAAAAT-3′, reverse 5′-CGCTTCACGAATTTGCGTGTCAT-3′; miR-204 forward, 5′-CTGTCACTCGAGCTGCTGGAATG-3′, and reverse, 5′-ACCGTGTCGTGGAGTCGGCAATT-3′. Identification of the miR-204 promoter and validation of promoter DNA methylation with TCGA data As intronic miRNAs are transcribed as a result of their host genes, the putative promoter region of miR-204 was predicted from the genomic sequence of the miR-204 host gene TRPM3 for the 2,000 upstream and 1,000 downstream sequences around its first exon based on the validated National Center for Biotechnology Information (NCBI) nucleotide sequences ([44]Figure S1). The promoter CpG islands were predicted using designing primers for methylation PCRs at Chinese Academy of Medical Sciences (MethPrimer, [45]http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi). DNA methylation data of TCGA samples were downloaded and processed with the database of DNA methylation and gene expression in human cancer (MethHC, [46]http://MethHC.mbc.nctu.edu.tw). The methylation levels of the identified CpG sites in tumor and normal tissues were analyzed on the basis of the gene chip data (methylation 450K) downloaded from the University of California Santa Cruz Xena Cancer Browser (UCSC) ([47]http://xena.ucsc.edu/). DNA extraction and methylation analyses Genomic DNA was extracted from fresh tissues using a QIAamp DNA Mini Kit (Qiagen NV, Venlo, the Netherlands) and modified by sulfite with an EZ DNA Methylation Kit (ZYMO, Irvine, CA, USA) according to the manufacturer’s protocols. Modified DNA was used for PCR in a total reaction volume of 10 µL as follows: ten cycles of melting (45 seconds at 94°C), annealing (48 seconds at 62°C), and extension (60 seconds at 72°C) and 35 cycles of melting (45 seconds at 94°C), annealing (48 seconds at 57°C), and extension (60 seconds at 72°C). PCR primers were designed using the online tools of Agena Bioscience, Inc. (San Diego, CA, USA) ([48]https://agenacx.com/online-tools/) to amplify miR-204 promoter regions with CpG islands. Two promoter regions (regions 9 and 15 that contain 50 and 18 CpG sites, respectively) were selected for the next amplification. The following primers were used: region 9, 5′-GTGGGTTTTGTATTTTTTAGAGAAG-3′ and 3′-AACCCCTACTTAAAACTTAACCTTTTCC-5′ and region 15, 5′-GTTTTTTTAAGGGTGAGAGTTAGGG-3′ and 3′-CAAACACCTAAAATATTCTTACTATTCTC-5′. PCR products (2 µL) modified by alkaline phosphatase were used as a template for in vitro transcription. RNase A cleavage (MassCleave Kit; Agena Bioscience, Inc.) was used for the reverse. The purified samples were spotted on a 384-pad SPectroCHIP (Agena Bioscience, Inc.) using a MassARRAY Nanodispenser RS 1000 (Agena Bioscience, Inc.). The spotted chips were placed in the MassARRAY Compact System (Agena Bioscience, Inc.) for detection, followed by spectral acquisition on a MALDI-TOF analyzer (Agena Bioscience, Inc.). The test data were analyzed and quantified by EpiTYPER software (Agena Bioscience, Inc.). To validate the promoter methylation regulation of miR-204, PTC cell lines were treated with the DNA methylation inhibitor 5-aza-2′-deoxycytidine (5-Aza). Cells were seeded in six-well plates, allowed to attach for 96 hours, and treat with different dose of 5-Aza. miR-204 expression levels were then detected by qRT-PCR. TCGA RNA sequencing and clinical data collection The TCGA thyroid carcinoma dataset (TCGA-THCA) consisted of 507 cases of PTC and 59 normal tissue samples. Merged level 3 mRNA and miRNA sequencing data and clinical information were downloaded from cBioPortal for Cancer Genomics ([49]http://www.cbioportal.org/index.do) and Firebrowse ([50]http://firebrowse.org). PTC cases and normal samples with both complete mRNA and miRNA sequencing and clinical data were identified. The RNAseq by expectation-maximization and reads per kilobase per million (RPKM) miRNA mapped values were used to quantify mRNA and miRNA expression levels, respectively. Screening differentially expressed miRNAs (DEMis) and genes The Limma package (version 3.32.5) of R (version 1.1.143) was used to identify DEMis between PTC and normal tissues from TCGA. According to the median expression values of miR-204, patients from TCGA were divided into the miR-204 high expression group and low expression group. The Limma package of R was used to identify differentially expressed genes (DEGs) between these two groups. Hierarchical analyses of DEMis identified from the miRNA microarray and TCGA data were performed by Cluster 3.0 and then visualized via Java TreeView 1.16r4 ([51]http://www.treeview.net). Gene functional eanalysis Bioinformatics methods were utilized to identify the pathways that are significantly associated with the downregulation of miR-204. We performed functional annotation clustering using the Database for Annotation, Visualization and Integrated Discovery ([DAVID], [52]https://david.ncifcrf.gov/) online system. Application of ggplot2 (version 2.2.1) in R studio was used to visualize the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. To further investigate the associations between downregulation of miR-204 and cancer-related pathways, we conducted gene set enrichment analysis (GSEA) using GSEA JAVA version 2.1 ([53]http://www.broadinstitute.org/gsea). GSEA calculates a gene set enrichment score that is estimated from predefined gene sets. The expression level of miR-204 was served as a phenotype, and KEGG pathway enrichment maps were used for visualization of the GSEA results. Statistical analysis All statistical analyses were performed using the SPSS 21.0 package and R studio. Application of Limma and ggplot2 in R studio was used as described above. Continuous variables are presented as the mean ± SD and were analyzed using the Student’s t-test or one-way ANOVA. Survival analysis was performed using Kaplan–Meier and log-rank tests. A P-value of <0.05 was considered statistically significant. Results Downregulation of miR-204 in PTC is associated with aggressive clinicopathologic characteristics Five pairs of PTC tumor and para-tumor tissues were analyzed for expression of miRNAs using the most updated miRNA array panel containing 2,549 human mature miR-NAs (release 21.0 from the miRBase database). With a fold change >2 and P-value <0.05 as the cutoff point, 181 DEMis (89 downregulated and 92 upregulated miRNAs; [54]Figure 1) were detected in this study. Among these DEMis, miR-146b, miR-551b, miR-375, miR-31, and miR-221 were the most significantly upregulated miRNAs, whereas miR-5703, miR-630, miR-718, miR-4646, and miR-6794 were the most significantly downregulated miRNAs (ranked by fold change values, for more details, see [55]GSE113629). Then, we screened DEMis with TCGA data. Differential analysis revealed 86 DEMis between PTC and normal tissues, of which 66 were downregulated and 20 were upregulated. We confirmed only six downregulated miRNAs in this independent and large TCGA cohort (miR-139, miR-486, miR-3652, miR-1275, miR-144, and miR-204). Consistent with the result from TCGA-THCA, qRT-PCR results indicated that miR-204 expression was markedly downregulated in thyroid cancer tissues ([56]Figure 2A). Notably, miR-204 was also lowly expressed in PTC cell lines compared to its expression level in normal cell line ([57]Figure 2B). Dysregulated miR-204 has been observed in various types of cancers, including thyroid cancer.[58]^10^,[59]^11 However, the underlying mechanisms of miR-204 in PTC remain to be further elucidated. To address the clinicopathological significance of the decreased expression of miR-204 in PTC, we used clinical data from TCGA PTC cohort. The data revealed that downregulation of miR-204 was frequently found in PTC tissues characterized by extra-thyroidal extension, high T stage, lymph metastasis, BRAF V600E mutation, and aggressive tall cell variant ([60]Table 1). By contrast, no significant differences were found based on patient gender, age, tumor multifocality, M stage, and TNM stage. Furthermore, we divided patients into two groups depending on the median expression value of miR-204 (miR-204 low expression group and high expression group). However, Kaplan–Meier survival analysis demonstrated that patients with lower miR-204 expression showed no differences in disease-free and overall survival rates compared to patients with higher miR-204 expression. Figure 1. [61]Figure 1 [62]Open in a new tab Heatmap of the 181 dysregulated miRNAs. Notes: Hierarchical analysis of the 181 dysregulated miRNAs based on the comparison of their expression values in PTC and normal tissues. Data have been deposited in GEO (accession code [63]GSE113629). Abbreviations: GEO, Gene Expression Omnibus; miRNAs, microRNAs; PTC, papillary thyroid carcinoma; T, tumor tissues; N, normal tissues. Figure 2. [64]Figure 2 [65]Open in a new tab Downregulated miR-204 expression in PTC tissues and cell lines. Notes: (A) miR-204 expression levels are downregulated in human PTC tissues compared to those in paired adjacent normal tissues. (B) miR-204 expression levels are downregulated in human PTC cell lines compare to that in normal cell line. *P<0.05. Abbreviation: PTC, papillary thyroid carcinoma. Table 1. The relationship between miR-204 expression and clinicopathological features in PTC patients from TCGA data set Variables Number of patients Mean expression of miR-204 (log2[RPKM]) P-values __________________________________________________________________ Gender  Male 136 4.28±2.15 0.763  Female 369 4.08±2.14 Age (years)  ≤45 238 4.21±2.13 0.617  >45 267 4.06±2.16 Histologic type  Follicular variant 102 5.83±2.16 <0.001[66]^*,[67]^a  Classic variant 357 3.78±1.95  Tall cell variant 37 2.66±1.03 Multifocality  Unifocal 269 4.06±2.14 0.932  Multifocal 228 4.21±2.13 Extrathyroidal extension  None 336 4.45±2.17 <0.001[68]^*  Minimal (T3) + moderate/advanced (T4a) 151 3.30±1.74 T stage  T1 + T2 311 4.42±2.20 0.005[69]^*  T3 + T4 192 3.66±1.94 N stage  N0 231 4.74±2.26 <0.001[70]^*  N1 224 3.38±1.70 M stage  M0 281 3.98±1.98 0.267  M1 9 4.04±2.52 AJCC TNM stage  Stage I + II 336 4.36±2.17 0.062  Stage III + IV 167 3.67±2.01 BRAF mutation  Without mutation 265 5.15±2.19 <0.001[71]^*  BRAF V600E 240 3.00±1.40 Disease-free status  Disease free 444 4.16±2.12 0.154[72]^b  Recurred/progressed 47 3.84±2.36 Overall survival status  Living 489 4.13±2.14 0.355[73]^b  Deceased 16 4.30±2.03 [74]Open in a new tab Notes: Mean expression of miR-204 is represented as the mean value ± SD. Patients with missed data were excluded. ^* P<0.05. ^a P-value by one-way ANOVA. ^b P-value by log-rank test. Abbreviations: AJCC, American Joint Committee on Cancer; PTC, papillary thyroid carcinoma; RPKM, per kilobase per million; TCGA, the Cancer Genome Atlas. Downregulated expression of miR-204 in PTC is correlated with promoter DNA methylation of its host gene TRPM3 Little is known about the genetic and regulatory factors contributing to the altered expression of miR-204 in PTC. Approximately, half of all miRNA genes are encoded in the introns of protein-encoding genes, and these miRNA genes are likely to be susceptible to transcriptional repression by aberrant DNA methylation of CpG islands located in the host genes.[75]^12 The transcriptional regulation of miR-204 may be controlled by the regulatory mechanism of its host gene. We found that miR-204 is located within the introns of TRPM3 (located at 9q21, [76]Figure 3A) using the UCSC genome browser. The chromatin patterns, genomic structure of the miR-204 and TRPM3 gene, and the location of neigh-boring genes are shown in [77]Figure S1. TCGA data indicated a coordinated downregulation of miR-204 expression with that of host gene TRPM3 (with log[2]FCs =−9.99 and −3.41 respectively, data not shown). Two CpG islands around the miR-204 promoter were identified ([78]Figure 3B). Therefore, we assessed the methylation status of promoter CpG islands of miR-204/TRPM3. Agena MassARRAY results indicated that 12 CpG sites located at the promoter of miR-204 were hypermethylated in primary tumor tissues compared to the matched normal tissues ([79]Table 2). Given the limited sample size of our study, data from TCGA database were used to confirm the DNA methylation status of the promoter of miR-204. Four CpG sites (cg00840310, cg02654358, cg14223966, and cg23553442) in the promoter region of TRPM3 were identified, which are also parts of our assessed CpG islands. Data from MethHC showed a significant negative correlation between miR-204 expression and methylation levels in its promoter region ([80]Figure 4A). Furthermore, three of the four CpG sites were strongly hypermethylated in PTC tissues compared to normal tissues ([81]Figure 4B). As shown in [82]Figure 4C and D, pretreatment with 5-Aza in PTC cells for 48 hours or 96 hours led to remarkable miR-204 expression, indicating a significantly negative correlation between niR-204 expression and methylation levels. Importantly, we observed higher methylation rates of the miR-204 promoter in tumor tissues than in normal tissues in 13 cancers of the total 18 cancers with the data deposited in MethHC ([83]Figure 5). Among these 13 types of cancers, downregulated miR-204 expression was also observed in nine cancers (data not shown). Thus, we revealed that promoter DNA methylation is a common mechanism of miR-204 downregulation in cancers, including PTC. Figure 3. [84]Figure 3 [85]Open in a new tab Promoter DNA methylation analysis of miR-204. Notes: (A) Schematic representation of the location of miR-204 in chromosome 9. miR-204 is located within the introns of TRPM3 (located at 9q21). (B) CpG island prediction and primer design of the miR-204 promoter region. Two CpG islands around the miR-204 promoter were identified using MethPrimer. Regions 9 and 15 that completely cover the two CpG islands were selected for the amplification. Each numbered region represents a region that can be amplified. Table 2. Comparison of methylation levels in CpG sites between PTC and adjacent normal tissues CpG sites Number of pairs Tissue types Mean methylation level P-value __________________________________________________________________ 9-CpG-1, 2, 3 10 PTC 0.292±0.097 <0.001 NT 0.128±0.168 9-CpG-5, 6 10 PTC 0.273±0.078 <0.001 NT 0.084±0.020 9-CpG-19, 20 10 PTC 0.272±0.068 <0.001 NT 0.127±0.031 9-CpG-21 10 PTC 0.281±0.060 <0.001 NT 0.121±0.059 9-CpG-35 10 PTC 0.303±0.220 0.005 NT 0.068±0.084 9-CpG-36, 37 10 PTC 0.287±0.035 <0.001 NT 0.127±0.032 9-CpG-38, 39 10 PTC 0.244±0.068 <0.001 NT 0.084±0.020 9-CpG-40, 41 10 PTC 0.392±0.056 <0.001 NT 0.195±0.106 9-CpG-44 10 PTC 0.292±0.055 <0.001 9-CpG-47 10 NT 0.121±0.059 <0.001 PTC 0.299±0.045 15-CpG-8 10 NT 0.121±0.056 0.009 PTC 0.463±0.096 15-CpG-13 10 NT 0.0279±0.172 <0.001 PTC 0.479±0.050 NT 0.299±0.090 [86]Open in a new tab Notes: CpG sites represent region number-CpG site numbers. Mean methylation level is represented as the mean value ± SD. Abbreviations: NT, normal tissue; PTC, papillary thyroid carcinoma. Figure 4. [87]Figure 4 [88]Open in a new tab Promoter DNA methylation analysis of miR-204 with TCGA data. Notes: (A) Correlation analysis between miR-204 expression levels and promoter methylation levels. The expression levels of miR-204 in tumor tissues were inversely correlated with the promoter DNA methylation rates. (B) Methylation levels of CpG sites from TCGA data. Three of the four CpG sites were hypermethylated in tumor tissues compared to normal tissues. One CpG site was hypomethylated (cg14223966). (C and D) qRT-PCR demonstrated increased expression of miR-204 in PTC cells after treatment with 5-Aza-dye for 48 and 98 hours compared to mock-treated cells. *P<0.05. Abbreviations: 5-Aza, 5-aza-2′-deoxycytidine; PTC, papillary thyroid carcinoma; qRT-PCR, quantitative real-time PCR. Figure 5. [89]Figure 5 [90]Open in a new tab miR-204 promoter methylation rates in multiple cancers. Notes: In 13 types of cancers, higher methylation rates were observed in cancer tissues than in normal tissues. *P<0.05, **P<0.005. Identification of related pathways involved in the downregulation of miR-204 Differential analysis revealed 1,273 DEGs (615 upregulated genes and 658 downregulated genes) between the miR-204 low expression group and high expression group. To gain a better understanding of gene function and signaling pathways related to the downregulation of miR-204 in PTC, we used the DAVID online tool to classify these DEGs using the GO and KEGG pathway references. For miR-204-related DEGs, the top enriched GO