Progress in Microbes and Molecular Biology Original Research Article 1 In-depth characterization of miRNome in papillary thyroid cancer with BRAF V600E mutation Azliana Mohamad Yusof1, Francis Yew Fu Tieng2, Rohaizak Muhammad3, Shahrun Niza Abdullah Suhaimi3, Isa Mohamed Rose3, Sazuita Saidin1, Rahman Jamal1, Imilia Ismail5, Nurul-Syakima Ab Mutalib1* 1UKM Medical Molecular Biology Institute (UMBI), Universiti Kebangsaan Malaysia, Jalan Yaacob Latif, 56000 Cheras, Kuala Lumpur, Malaysia 2Department of Biomedical Science, Faculty of Medicine and Health Sciences, Universiti Putra Malaysia, 43400 Serdang, Selangor Darul Ehsan, Malaysia 3Department of Surgery, Faculty of Medicine, Universiti Kebangsaan Malaysia, Jalan Yaacob Latif, 56000 Cheras, Kuala Lumpur, Malaysia 4Department of Pathology, Faculty of Medicine, Universiti Kebangsaan Malaysia, Jalan Yaacob Latif, 56000 Cheras, Kuala Lumpur, Malaysia 5School of Biomedicine, Faculty of Health Sciences, Universiti Sultan Zainal Abidin, Gong Badak Campus, 21300 Kuala Terengganu, Terengganu, Malaysia Abstract: MicroRNAs (miRNAs) are small non-coding RNAs, which play a critical regulatory role in papillary thyroid carcinoma (PTC). BRAF V600E is a hotspot mutation occurring in a majority of PTC cases and is proposed to be associated with poor clinical outcomes. The relationship between BRAF V600E status and miRNA expression in PTC has not been comprehensively studied. In this study, we aimed to identify the differentially expressed miRNAs in PTCs with and without BRAF V600E in an unbiased manner. Five fresh frozen thyroid cancer tissues paired with their respective adjacent normal tissues from PTC patients were subjected to BRAF V600E genotyping using Sanger sequencing and small RNA deep sequencing (miRNAseq). MiRNAs differentially expressed between BRAF V600E-positive and BRAF V600E-negative PTC tissues were validated in silico using The Cancer Genome Atlas (TCGA) THCA datasets containing 420 samples. MiRNA target prediction and pathway enrichment analysis were performed to identify biological pathways altered in this cancer. We identified 174 differentially expressed miRNAs; 80 were significantly over-expressed, while 94 were under- expressed (adj. p-value < 0.1; log2 fold change ≤ -1 or ≥ 1). Fifteen miRNAs were significantly differentially expressed only in BRAF V600E-positive PTC, and eight of these were validated in TCGA THCA dataset (hsa-miR-212, -132, -135b-3p/5p, -200b, -200a-3p/5p, -27a-3p/5p, -29a and -1296). Subsequent analysis revealed significant enrichment of cancer-related pathways including proteoglycans in cancer, ECM-receptor interaction and MAPK pathways in BRAF V600E-positive PTC. Using the miRNAseq and in silico validation using TCGA THCA study, we identified eight miRNAs that were differentially expressed in PTC tissues with BRAF V600E. This study also complemented the existing knowledge about deregulated miRNAs in PTC development. Keywords: microRNA; papillary thyroid cancer; BRAF V600E, next-generation sequencing. Received: 9th October 2019 Accepted: 10th November 2019 Published Online: 18th March 2020 Citation: Mohamad Yusof A, Tieng FYF, Muhammad R, et al. In-depth characterization of miRNome in papillary thyroid cancer with BRAF V600E mutation. Prog Microbes Mol Biol 2019; 2(1): a0000043. https://doi.org/10.36877/pmmb. a0000043 Introduction The incidence rate of thyroid cancer had risen significantly in many countries worldwide in the past decade. In 2012, there were approximately 230,000 and 70,000 new cases of thyroid cancer among women and men, respectively [1,2]. According to the United States Cancer Statistics 2019, 52,070 new thyroid cancer cases with 2,170 deaths were estimated by 2019 in both sexes [3]. Most thyroid cancers originated from follicular epithelial cells, which were further divided into well-differentiated papillary thyroid carcinoma and follicular carcinoma, poorly differentiated carcinoma and anaplastic carcinoma [4–6]. Papillary thyroid carcinoma (PTC) is the most prevalent histological type that contributed to 85 to 90% of reported cases [7]. Several of its clinicopathological features are correlated to poor prognosis in PTC patients, which include the male gender, larger tumour size, older age, extrathyroidal extension, thyroid capsule invasion, Copyright 2019 by Mohamad Yusof A et al. and HH Publisher. This work under licensed under the Creative Commons Attribution- NonCommercial 4.0 International Lisence (CC-BY-NC 4.0) *Correspondence: Nurul-Syakima Ab Mutalib, UKM Medical Molecular Biology Institute (UMBI); syakima@ppukm.ukm.edu.my 2 lymph node metastasis as well as BRAF mutation [8]. Although PTC patients had a high survival rate [1,9], these prognostic factors had been shown to affect the overall survival among PTC patients [10]. MiRNAs are small non-coding RNAs that comprise of 19 to 22 nucleotides, which negatively regulate gene expression. Each miRNA can take part in many cellular pathways, and thus, miRNAs are involved in many different diseases [11] including thyroid cancers [12–15]. In addition, it was reported that different histopathological types of thyroid tumours have discrete miRNA profiles [16]. miR-146b [17,18], miR-221 and miR-222 are among the commonly upregulated miRNAs in papillary thyroid carcinoma [16,19–23]. These non-coding RNAs are promising biomarkers to identify aggressive PTC cases [17,24]. BRAF is a serine-threonine kinase that is activated by RAS binding and protein recruitment to the cell membrane [25,26]. Activation of MEK along with the MAPK signalling pathway is activated by BRAF phosphorylation [27]. The most frequent genetic changes in PTC are point mutations of BRAF which are observed in 35 to 70% of PTC cases [26,28]. More than 95% of BRAF mutations detected in thyroid cancers are thymine to adenine transversion at position 1799 (T1799A), resulting in the substitution of valine by glutamate at residue 600 (V600E) [29–31]. Various studies had shown that the BRAF V600E mutation was associated with lymph node metastasis, therefore, it had received the attention as a diagnostic and prognostic molecular marker in recent years to improve the diagnosis and identify individuals at increased risk of PTC recurrence [30,32]. In the past few years, the evolution of next-generation sequencing technologies has allowed global expression profiling of miRNAs [33] and the discovery of novel human miRNAs [14,34]. In this study, we aimed to identify the differentially expressed miRNAs in PTCs with BRAF V600E using small RNA sequencing. To achieve this objective, five fresh-frozen tumour tissues paired with their normal adjacent thyroid tissues were included in the study. The mutational analysis of BRAF V600E gene and differentially expressed miRNAs were performed using Sanger sequencing and small RNA sequencing, respectively. The results were validated with an in silico approach using the dataset obtained from The Cancer Genome Atlas (TCGA) THCA. Materials and Methods Clinical Specimens Five pairs of tumour-adjacent normal fresh frozen tissues were collected from patients diagnosed with PTC from the UKM Medical Centre (UKMMC). This study was approved by the Universiti Kebangsaan Malaysia Research Ethics Committee (UKMREC; UKM 1.5.3.5/244/UMBI-2015-002). Informed consent was obtained from all the study participants. The tissues were dissected, snap-frozen and stored in liquid nitrogen. All samples were cryosectioned and stained using haematoxylin and eosin and the percentage of tumour cells and normal cells contents were assessed by a pathologist. Only tumour samples with at least 80% cancerous cells and normal adjacent thyroid tissues with less than 20% necrosis were selected for further analysis. Nucleic Acid Isolation All fresh frozen tissues were subjected to nucleic acid extraction using Allprep DNA/RNA/miRNA Universal Kit (Qiagen, Hilden, Germany) according to the manufacturer’s recommendations. The integrity of RNA was assessed using Agilent BioAnalyzer 2100 (Agilent Technologies, CA, U.S.), while the quality of DNA was assessed using 1% agarose gel. The quantity and purity of RNA and DNA were assessed using Qubit 2.0 fluorometer (Thermo Scientific, MA, U.S.) and Nanodrop 2000c Spectrometer (Thermo Scientific, MA, U.S.), respectively. BRAF V600E Genotyping PCR amplification of genomic regions of interest was performed using BRAF V600E forward primer 5’-TGCTTGCTCTGATAGGAAAATG-3’ and BRAF V600E reverse primer 5’-AGCATCTCAGGGCCAAAAAT-3’ [35]. Amplification was performed in a reaction volume of 25 μl containing 50 ng DNA template, 10 μM each primer (Integrated DNA Technologies, IA, U.S.),10x PCR Gold Buffer without MgCl2, dNTP Mix (10 mM), MgCl2 solution (25 mM) , AmpliTaq Gold® (5U/μl) (Applied Biosystems, CA, U.S.) and nuclease-free water. PCR conditions were as follows; 95°C for 4 minutes; 35 cycles of 95°C for 45 seconds, 50°C for 30 seconds and 72°C for 1 minute; 72°C for 5 minutes; and held at 4°C. PCR products were visualized by electrophoresis on 1.5% agarose gel with an expected size of ~ 228 bp. PCR purification was conducted using QIAquick PCR Purification Kit (Qiagen, Hilden, Germany) as per manufacturer’s instruction. Subsequently, DNA sequencing was performed using ABI Prism 3130xl Genetic Analyzer (Applied Biosystem, CA, USA). Library Preparation and miRNA Sequencing RNA samples from tumour samples and their adjacent normal tissues were processed into libraries using TruSeq Small RNA Sample Prep Kit (Illumina, CA, USA). Briefly, 3′ and 5′ adapters were sequentially ligated to the ends of small RNAs fractionated from 2 μg of total RNA, and reverse transcribed to generate cDNA. The cDNA was amplified using a common primer complementary to the 3′ adapter, and a primer containing 1 of 48 index sequences. Samples were size-selected (140–160 bp fragments) on a 6% polyacrylamide gel, purified, quantified and pooled for multiplexed sequencing. The resulting pooled libraries were normalized to 2 nM and were hybridized to oligonucleotide-coated single-read flow cells for cluster generation using HiSeq® Rapid SR Cluster Kit v2 on Hiseq 2500 (Illumina, CA, USA). Subsequently, the clustered pooled miRNA libraries were sequenced on the HiSeq 2500 for 50 sequencing cycles using HiSeq® Rapid SBS Kit v2 (50 Cycle) (Illumina, CA, USA). In-depth characterization... 3 Bioinformatics and Statistical Analyses Pre-processing of data was executed in BaseSpace software (Illumina, CA, USA), and FASTQ files were generated. miRNA Analysis app version 1.0.0 was used for determination of differentially expressed miRNAs using the workflow described by Cordero et al, 2012 [25]. Briefly, the pipeline includes 3’ end adapter removal using Cutadapt, annotation to miRBase v21, mapping using SHRIMP aligner and differential analysis of miRNAs using DESeq2. Benjamini and Hochberg’s [36] correction was applied to ensure a false discovery rate (FDR) less than 0.1 and absolute log2 fold change ≤ -1 or ≥ 1 were considered for further analysis. Heatmaps were created using GeneE from the Broad Institute (http://www.broadinstitute.org/ cancer/software/GENE-E). MiRNA target prediction via DIANA-TarBase v7.0 [37] and pathway enrichment analysis was performed using DIANA-miRPath v3.0 [38]. Other statistical analyses were performed using GraphPad Prism 6 unless stated otherwise. In silico Validation We used the TCGA-generated level 3 miRNA sequencing data from THCA project [24]. These data were accessed from ‘https://tcga-data.nci.nih.gov/tcga/dataAccessMatrix. htm’ on April 13, 2016. The normalised expression (reads per million or RPM) of all miRNAs was log2-transformed and used for fold change calculation. We then performed the Students’ unpaired t-test with Benjamini Hochberg false discovery rate (FDR) multiple testing correction and log2 fold change calculation using Bioconductor version 3.1 (BiocInstaller 1.18.2) [39] in R version 3.2.0 (R Development Core Team, 2008). Results Demographic Data All patients in the small discovery set were women with a mean age of 51.6 years. Majority of the tumours were located in the right lobe of the thyroid and all tumours were larger than 1 cm. All of the patients had lymph node metastasis at diagnosis. Validation cohort from TCGA THCA studies included 229 BRAF V600E- positive PTC, 132 BRAF V600E-negative PTC and 59 unpaired normal thyroid samples. miRNAseq Analysis With the rapid run mode using HiSeq Rapid SBS Kit v2, we achieved an average of 5.6 M reads per sample with Q30 (3,779,968 to 8,308,285 reads in each sample). The percentage of mapped reads in normal samples were significantly lower than in tumour samples (Supplementary Table 1). Figure 1 illustrated the quality control statistics of miRNAseq experiment for BRAF V600E-positive versus adjacent normal samples. In Figure 1(A), approximately 45 to 75% of the reads were identified as isomiR (known precursor) and 25 to 55% were the already known mature miRNAs. Figure 1(B- D) illustrated the Principle Component Analysis (PCA) plots of miRNAs family, mature miRNAs and precursor miRNAs expression profiles. PCA analysis revealed that samples formed distinct clusters, withall tumour and normal samples clustering based on their respective groups. Differentially Expressed miRNAs There were 174 miRNAs significantly differentially expressed in PTC BRAF V600E-positive versus their normal adjacent thyroid tissues. Eighty miRNAs were upregulated, while 94 miRNAs were downregulated (log2 fold change ≤ -1 or ≥ 1; FDR p-value < 0.1). The volcano plot (Figure 2) illustrated the significantly differentially expressed miRNAs. In addition, the unsupervised hierarchical clustering analysis and heatmap illustrated in Figure 3, clearly demonstrated the difference of the miRNA expression profiles between tumour and normal samples. Mohamad Yusof A et al. Figure 1. Quality control statistics of miRNAseq experiment. (A) Distribution of miRNA sequences among the subcategories in BRAF V600E PTC and their adja- cent normal. (B-D) Principle Component Analysis (PCA) plots of miRNAs family, mature miRNA and precursor miRNAs. The plots showed that global miRNAs family, mature miRNAs and precursor miRNAs expression pattern clearly differentiate the samples according to respective group 4 In-depth characterization... Figure 2. MiRNAs with statistical significance after a t-test are shown in red and blue on a volcano plot above. Those with no significance are shown in black. The red dots represent significantly upregulated miRNAs while blue dots represented significantly downregulated miRNAs. Figure 3. Hierarchical clustering and heatmap representation of differentially expressed miRNAs in discovery set. The list of differentially expressed miRNAs were filtered using FDR-adjusted p value <0.1, absolute log2 fold change ≥1 or ≤-1. In BRAF V600E-positve PTC versus normal-adjacent tissues, miRNAseq revealed 174 differentially expressed miRNAs. Ninety-four (94) were under-expressed while 80 were over-expressed. 5 Mohamad Yusof A et al. Pathway Enrichment Analysis of the Deregulated miRNAs A single miRNA might play specific roles in the pathogenesis of PTC; however, the miRNA pathway as a whole might also be of importance. Therefore, DIANA- mirPath [37], a web-based computational tool, was used to identify pathways that were potentially altered by the expression of multiple miRNAs, and to incorporate miRNAs into molecular pathways. Based on the 80 upregulated miRNAs, there were 44 significantly enriched pathways identified by the meta-analysis algorithm using “Pathway union” option (p-value < 0.05) (Figure 4A). Thirty-eight out of the 80 upregulated miRNAs were involved in the proteoglycans in the cancer pathway, targeting 176 genes. Figure 4B illustrated the hierarchical clustering analysis and the heatmap derived from the 44 significantly enriched pathways. The number of miRNAs and experimentally validated predicted targets involved in each pathway were tabulated in Table 1. We also performed the same analysis on 94 significantly downregulated miRNAs. Only five significantly enriched pathways were identified: fatty acid biosynthesis, fatty acid metabolism, ECM-receptor interaction, fatty acid elongation and proteoglycans in cancer pathway with FDR adj. p-value < 1E-325, 2.57E-07, 7.70E-05 and 0.000352, respectively. TCGA THCA in silico Validation To validate the observed significance of the 174 miRNAs in BRAF V600E-positive PTC as compared to their normal adjacent thyroid tissues, we analysed the expression of these miRNAs in TCGA THCA dataset. Comparison between BRAF V600E-positive PTC versus normal thyroid tissues revealed 123 significantly upregulated and 319 downregulated miRNAs (Supplementary Table 2). To identify miRNAs which were deregulated only in BRAF V600E-positive PTC, the intersection of the differentially expressed miRNAs was performed (Figure 5). Since TCGA THCA dataset did not comprehensively annotate the miRNAs based on 3p and 5p arm, the intersection was performed without taking into account the arm annotation. From the intersection, 15 miRNAs were deregulated in BRAF V600E-positive PTC versus normal thyroid tissues in both and TCGA THCA datasets and the current study. From these 15 miRNAs, eight miRNAs were in concordance with TCGA THCA data expression levels which included hsa-miR-212, hsa-miR-132, hsa-miR- 135b-3p and 5p, hsa-miR-200b, hsa-miR-200a-3p and 5p, hsa-miR-27a-3p and 5p, hsa-miR-29a and hsa- miR-1296 (log2 fold change 1.59, 1.43, 2.36, 1.72, 2.12, 1.31, 1.52, 1.60, 1.69 and -1.12; adj. p-value < 0.1), respectively (Figure 6). Figure 4. Pathway enrichment analysis of 80 significantly upregulated miRNAs in BRAF V600E-positive PTC and miRNAs versus pathway heatmap (clustering based on sig- nificance levels) using the DIANA miRPath v3.0 interface. (A) A list of 44 significant KEGG targeted pathways with their log10 p-value. (B) Hierarchical cluster and heatmap. Dark red indicate lower significant values. The attached dendograms on both axes represent the hierarchical clustering results for miRNAs and pathways respectively. MiRNAs are clustered together by exhibiting similar pathway targeting patterns and pathways are clustered together by related miRNAs. 6 In-depth characterization... Figure 5. Venn diagram of the differentially expressed miRNAs in discovery and in silico validation dataset. From the intersection, there were 15 miRNAs that were differentially expressed in BRAF V600E-positive PTC versus normal. From these 15 miRNAs, eight deregulated miRNAs from discovery data were in concordance with TCGA THCA data (seven were upregulated and one was downregulated). Figure 6. Box plot of eight concordance miRNAs. All miRNAs were significantly expressed (adjusted p-value <0.1; log2 fold change ≤-1 or ≥1). 7 Mohamad Yusof A et al. Table 1. Number of miRNAs and experimentally validated targets involved in each of the 44 enriched pathways. KEGG pathway Adj p-value No. of genes No. of miRNAs Prion diseases <1E-325 18 8 MicroRNAs in cancer <1E-325 140 12 Oocyte meiosis <1E-325 83 19 TGF-beta signaling pathway <1E-325 62 19 Thyroid hormone signaling pathway <1E-325 95 20 ECM-receptor interaction <1E-325 50 21 Colorectal cancer <1E-325 59 21 Prostate cancer <1E-325 81 21 Fatty acid biosynthesis <1E-325 7 23 Pathways in cancer <1E-325 303 23 Fatty acid metabolism <1E-325 36 24 Cell cycle <1E-325 109 24 Protein processing in endoplasmic reticulum <1E-325 138 24 Hepatitis B <1E-325 118 25 Chronic myeloid leukemia <1E-325 70 27 Glioma <1E-325 58 28 p53 signaling pathway <1E-325 66 29 Viral carcinogenesis <1E-325 177 31 Lysine degradation <1E-325 39 33 Hippo signaling pathway <1E-325 119 35 Adherens junction <1E-325 66 37 Proteoglycans in cancer <1E-325 173 38 FoxO signaling pathway 1.75E-14 104 22 Bacterial invasion of epithelial cells 2.16E-14 65 21 Bladder cancer 3.91E-13 34 18 Endometrial cancer 4.77E-12 47 18 Small cell lung cancer 2.19E-10 72 17 Transcriptional misregulation in cancer 5.84E-10 133 17 Melanoma 9.28E-10 56 20 Signaling pathways regulating pluripotency of stem cells 7.09E-09 98 14 Thyroid cancer 1.00E-08 26 16 Endocytosis 1.43E-08 159 16 Pancreatic cancer 2.29E-07 62 13 Focal adhesion 6.68E-07 142 12 PI3K-Akt signaling pathway 1.00E-06 182 14 Ubiquitin mediated proteolysis 2.39E-06 105 14 Non-small cell lung cancer 5.32E-06 48 14 Renal cell carcinoma 1.47E-05 58 13 Estrogen signaling pathway 2.93E-05 60 8 Other types of O-glycan biosynthesis 0.000125 21 11 Central carbon metabolism in cancer 0.000955 50 9 Neurotrophin signaling pathway 0.001194 75 8 Shigellosis 0.002684 50 10 Steroid biosynthesis 0.009698 13 12 8 In-depth characterization... Discussion There was various type of BRAF mutations reported for malignant tumours including PTC such as BRAF V600E, BRAF V600D, BRAF V600Q, BRAF V600V and BRAF V600L [40,41]. For our study, we only focused on BRAF V600E mutation as it was the most common BRAF mutation in PTC, which comprised more than 90% of cases [42]. Many studies had demonstrated an association of this BRAF mutation with the aggressive clinicopathological characteristics of PTC such as extrathyroidal invasion, lymph node metastasis and recurrence of PTC [8,43]. It was suggested that the involvement of BRAF V600E mutation in the activation of RAS/RAF/MAPK pathway could result in higher deregulation of miRNA expression [44]. While previously published studies utilized microarray and real-time PCR, we used small RNA deep sequencing to determine the deregulation of miRNAs in BRAF V600E- positive PTC patients in an unbiased manner. There were more downregulated miRNAs in the tumour samples as compared to their adjacent normal thyroid tissues (80 upregulated versus 94 downregulated miRNAs). Deregulation of hsa-miR-146b-5p, hsa-miR-146b-3p, hsa- miR-222-3p, hsa-miR-221-3p, hsa-miR-204-5p and hsa- miR-7-5p were further reconfirmed in this study, signifying that dysregulation of these miRNA was common in PTC versus normal thyroid tissues. The association between BRAF V600E status and miRNA expression in PTC had been controversial. A large-scale analysis of TGCA data had demonstrated that the BRAF V600E mutation was one of the key drivers of PTC [24]. Other studies had shown that downregulation of hsa-miR-7-5p and hsa-miR-204-5p in PTC were associated with the BRAF V600E mutation [23,44]. There were however contradictory findings that showed that BRAF V600E mutation is not related with the aggressiveness of PTC and thus, cannot serve as prognosis marker for PTC [45–47]. MiR-200a and miR-200b were the members of miR- 200 family. Upregulation of these two miRNAs was observed in PTC [48], follicular thyroid carcinoma (FTC) and follicular adenoma [16], while being downregulated in anaplastic thyroid carcinoma [49]. It was suggested that miR-200b downregulates the tumour suppressor genes [48]. These miRNAs were shown to play a crucial role in tumour cell metastasis progression or epithelial-mesenchymal transition (EMT) [50] and inhibit angiogenesis [51]. Hsa-miR- 200a and hsa-miR-200b were also upregulated in thyroid cell-derived cell lines and tissues with BRAF V600E [16,48]. In this study, miR-200a and miR-200b were upregulated in BRAF V600E-positive cases in both discovery and TCGA THCA datasets, suggesting their relation to BRAF V600E mutation. Among the significantly enriched pathways in BRAF V600E-positive PTC were proteoglycans in cancer, cell cycle pathway and ECM-receptor interaction pathway. There were 173 genes with 38 miRNAs in proteoglycans in cancer, 109 genes with 24 miRNAs in cell cycle pathways and 50 genes with 21 miRNAs involved in the ECM- receptor interaction pathway. Proteoglycans (PGs) are key molecular constituents of the ECM and cell surfaces and play important roles in integrating signals from growth factors, chemokines and integrins, cell to cell interactions as well as matrix adhesion [52,53]. PGs act in a context- dependent manner; some have pro- and anti-angiogenic activities, while others can directly stimulate cancer growth by controlling key signalling pathways [52,54,55]. Aberrant accumulation of PGs in human thyroid cancer was first reported in 1984 [56]. Using the glycoproteomics approach, Arcinas and colleagues studied the expression of proteoglycans in various thyroid cancer cell lines. Two transmembrane heparan sulfate proteoglycans, syndecan-1 and syndecan-4, were uniquely expressed in FTC-133 and XTC-1 respectively. In addition, a GPI- anchored proteoglycan, glypican-1, was identified in three thyroid cancer cell lines (FTC-133, XTC-1 and DRO-1). The extracellular heparan sulfate proteoglycan basement membrane-specific core protein or perlecan was only detected in ARO, a dedifferentiated thyroid cancer cell line [57]. Numerous studies support the important role of proteoglycans as miRNA targets in cancer progression [58,59]. Deregulation of miRNAs results in atypical expression patterns of proteoglycans and their biosynthetic enzymes, thus leading to abnormal cell proliferation, apoptosis, adhesion, migration, invasiveness and epithelial-to-mesenchymal transition [60–62]. Therapeutic strategies targeting the microRNA– proteoglycan are emerging for cancers such as in melanoma and medulloblastoma [60,63,64]. While the relationship between miRNAs expression and ECM- receptor interaction pathway in regards to BRAF V600E has been reported [65,66], there is currently no published evidence linking miRNA regulation to proteoglycans in thyroid cancers, especially in the BRAF V600E-positive PTC. In order to validate our results in a larger dataset, we re-analysed TCGA THCA dataset containing 229 BRAF V600E-positive PTC, 132 BRAF V600E-negative PTC and 59 normal thyroid tissues. From the intersection of our discovery data with the TCGA data, 15 miRNAs were significantly deregulated in BRAF V600E-positive as compared to normal thyroid tissues and eight miRNAs were in concordance in term of expression levels. The remaining seven miRNAs showed the opposite trend of expression. The discrepancies could be due to the fact that TCGA THCA datasets used were from unmatched samples. Secondly, at the time of the data retrieval, TCGA THCA had not yet annotated the miRNAs according to their 3p and 5p arm and thus, the expression level between these two arms could not be differentiated. Nevertheless, we were able to reconfirm the expression of commonly deregulated miRNAs in PTCs and provide a new list of miRNAs related to BRAF V600E. Conclusion In conclusion, our study illustrated the interplay between BRAF V600E status and differentially expressed miRNAs in PTC. This information would add to the understanding of the molecular mechanisms of miRNAs in BRAF V600E-positive PTC. Although these findings were needed to be validated in larger sample size, they 9 could serve as a basis for the identification of a potential diagnostic or prognostic biomarker of PTC. In addition, functional studies to clarify further the mechanisms of miRNA regulation in BRAF V600E-positive PTC were warranted. Declaration of interest The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported. Funding This research was funded by the Fundamental Research Grant Scheme (FRGS) from the Ministry of Education Malaysia (FRGS/1/2014/SKK01/UKM/03/1). Author contributions A Mohamad Yusof, NS Ab Mutalib involved in the specimen collections, library preparation and sequencing, data analyses, acquisition of data and drafting the manuscript. FYF Tieng performed the TCGA analyses. S. Saidin performed the BRAF V600E genotyping. I. Mohamed Rose assessed tumour percentage of the tissues. SN Abdullah Suhaimi and R Muhammad were thyroid surgeons involved in specimen retrieval. I Ismail and R Jamal provided critical review on the manuscript. All authors read and approved the final manuscript. Acknowledgements The authors would like to thank the TCGA Thyroid Cancer Research Group for providing the data used in this study. We also thank Chin Minning, Field Application Specialist from Science Vision Sdn. Bhd. for the technical assistance during the next-generation sequencing experiment. References 1. Ferlay J, Soerjomataram I, Dikshit R, et al. Cancer incidence and mor- tality worldwide: Sources, methods and major patterns in GLOBO- CAN 2012: Globocan 2012. Int J Cancer, 2015; 136(5): E359–E386. doi:10.1002/ijc.29210 2. Torre LA, Bray F, Siegel RL, et al. Global cancer statistics, 2012. CA: Cancer J Clin, 2015; 65(2): 87–108. doi:10.3322/caac.21262 3. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA A Cancer J Clin, 2019; 69(1): 7–34. doi:10.3322/caac.21551 4. Nikiforov YE, Nikiforova MN. Molecular genetics and diagno- sis of thyroid cancer. Nat Rev Endocrinol, 2011; 7(10): 569–580. doi:10.1038/nrendo.2011.142 5. DeLellis R, Lioyd R, Heitz P, et al. WHO classification of tumors. In Pathology and genetics of tumours of endocrine organs. Lyon: IARC Press; 2004. Available from: https://publications.iarc.fr/Book-And- Report-Series/Who-Iarc-Classification-Of-Tumours/Pathology-And- Genetics-Of-Tumours-Of-Endocrine-Organs-2004. Accessed Novem- ber 7, 2019. 6. Nikiforov Y, Biddinger PW, Thompson LDR. Diagnostic pathology and molecular genetics of the thyroid. Philadelphia: Wolters Kluwer Health/Lippincott Williams & Wilkins; 2012. 7. Katoh H, Yamashita K, Enomoto T, et al. Classification and general considerations of thyroid cancer. Ann Clin Pathol, 2015; 3(1): 1045. 8. Haugen BR, Alexander EK, Bible KC, et al. 2015 American thyroid association management guidelines for adult patients with thyroid nod- ules and differentiated thyroid cancer: The American thyroid associa- tion guidelines task force on thyroid nodules and differentiated thyroid cancer. Thyroid, 2016; 26(1): 1–133. doi:10.1089/thy.2015.0020 9. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin, 2018; 68(6): 394–424. doi:10.3322/caac.21492 Mohamad Yusof A et al. 10. Chrisoulidou A, Boudina M, Tzemailas A, et al. Histological subtype is the most important determinant of survival in metastatic papillary thy- roid cancer. Thyroid Res, 2011; 4: 12. doi:10.1186/1756-6614-4-12 11. Mo Y-Y. MicroRNA regulatory networks and human disease. Cell Mol Life Sci, 2012; 69(21): 3529–3531. doi:10.1007/s00018-012-1123-1 12. Lu J, Getz G, Miska EA, et al. MicroRNA expression profiles classify human cancers. Nature, 2005; 435(7043): 834–838. doi:10.1038/na- ture03702 13. Pallante P, Visone R, Ferracin M, et al. MicroRNA deregulation in hu- man thyroid papillary carcinomas. Endocr Relat Cancer, 2006; 13(2): 497–508. doi:10.1677/erc.1.01209 14. Dettmer M, Perren A, Moch H, et al. Comprehensive MicroRNA expres- sion profiling identifies novel markers in follicular variant of papillary thyroid carcinoma. Thyroid, 2013; 23(11): 1383–1389. doi:10.1089/ thy.2012.0632 15. Pallante P, Battista S, Pierantoni GM, Fusco A. Deregulation of microR- NA expression in thyroid neoplasias. Nat Rev Endocrinol, 2014; 10(2): 88–101. doi:10.1038/nrendo.2013.223 16. Nikiforova MN, Tseng GC, Steward D, et al. MicroRNA expression pro- filing of thyroid tumors: Biological significance and diagnostic utility. J Clin Endocrinol Metab, 2008; 93(5): 1600–1608. doi:10.1210/jc.2007- 2696 17. Lee JC, Zhao JT, Clifton-Bligh RJ, et al. MicroRNA-222 and mi- croRNA-146b are tissue and circulating biomarkers of recurrent pap- illary thyroid cancer. Cancer, 2013; 119(24): 4358–4365. doi:10.1002/ cncr.28254 18. Chou C-K, Chen R-F, Chou F-F, et al. miR-146b is highly expressed in adult papillary thyroid carcinomas with high risk features including extrathyroidal invasion and the BRAF(V600E) mutation. Thyroid, 2010; 20(5): 489–494. doi:10.1089/thy.2009.0027 19. Swierniak M, Wojcicka A, Czetwertynska M, et al. In-depth character- ization of the microRNA transcriptome in normal thyroid and papillary thyroid carcinoma. J Clin Endocrinol Metab, 2013; 98(8): E1401–1409. doi:10.1210/jc.2013-1214 20. He H, Jazdzewski K, Li W, et al. The role of microRNA genes in papil- lary thyroid carcinoma. Proc Natl Acad Sci USA, 2005; 102(52): 19075– 19080. doi:10.1073/pnas.0509603102 21. Pallante P, Visone R, Croce CM, et al. Deregulation of microRNA ex- pression in follicular-cell-derived human thyroid carcinomas. Endocr Relat Cancer, 2010; 17(1): F91–104. doi:10.1677/ERC-09-0217 22. Tetzlaff MT, Liu A, Xu X, et al. Differential expression of miRNAs in papillary thyroid carcinoma compared to multinodular goiter using formalin fixed paraffin embedded tissues. Endocr Pathol, 2007; 18(3): 163–173. doi:10.1007/s12022-007-0023-7 23. Mancikova V, Castelblanco E, Pineiro-Yanez E, et al. MicroRNA deep-sequencing reveals master regulators of follicular and papillary thyroid tumors. Mod Pathol, 2015; 28(6): 748–757. doi:10.1038/mod- pathol.2015.44 24. Cancer Genome Atlas Research Network. Integrated genomic charac- terization of papillary thyroid carcinoma. Cell, 2014; 159(3): 676–690. doi:10.1016/j.cell.2014.09.050 25. Fakhruddin N, Jabbour M, Novy M, et al. BRAF and NRAS mutations in papillary thyroid carcinoma and concordance in BRAF mutations between primary and corresponding lymph node metastases. Sci Rep, 2017; 7(1): 4666. doi:10.1038/s41598-017-04948-3 26. Li X, Abdel-Mageed AB, Kandil E. BRAF mutation in papillary thyroid carcinoma. Int J Clin Exp Med, 2012; 5(4): 310–315. 27. Jiang L, Chu H, Zheng H. B-Raf mutation and papillary thyroid car- cinoma patients. Oncol Lett, 2016; 11(4): 2699–2705. doi:10.3892/ ol.2016.4298 28. Xing M. BRAF mutation in thyroid cancer. Endocr Relat Cancer, 2005; 12(2): 245–262. doi:10.1677/erc.1.0978 29. Davies H, Bignell GR, Cox C, et al. Mutations of the BRAF gene in human cancer. Nature, 2002; 417(6892): 949–954. doi:10.1038/na- ture00766 30. Xing M, Alzahrani AS, Carson KA, et al. Association between BRAF V600E mutation and recurrence of papillary thyroid cancer. J Clin On- col, 2015; 33(1): 42–50. doi:10.1200/JCO.2014.56.8253 31. Yarchoan M, LiVolsi VA, Brose MS. BRAF mutation and thyroid cancer recurrence. J Clin Oncol, 2015; 33(1): 7–8. doi:10.1200/ JCO.2014.59.3657 32. Lupi C, Giannini R, Ugolini C, et al. Association of BRAF V600E muta- tion with poor clinicopathological outcomes in 500 consecutive cases of papillary thyroid carcinoma. J Clin Endocrinol Metab, 2007; 92(11): 4085–4090. doi:10.1210/jc.2007-1179 33. Kozomara A, Griffiths-Jones S. miRBase: Integrating microRNA anno- tation and deep-sequencing data. Nucleic Acids Res, 2011; 39(Database issue): D152–157. doi:10.1093/nar/gkq1027 34. Kozomara A, Griffiths-Jones S. miRBase: Annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res, 2014; 42(Database issue): D68–73. doi:10.1093/nar/gkt1181 35. Kondo T, Nakazawa T, Murata S, et al. Enhanced B-Raf protein expres- sion is independent of V600E mutant status in thyroid carcinomas. Hum Pathol, 2007; 38(12): 1810–1818. doi:10.1016/j.humpath.2007.04.014 36. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practi- cal and powerful approach to multiple testing. J R Stat Soc Ser B (Meth- odol), 1995; 57(1): 289–300. 37. Vergoulis T, Kanellos I, Kostoulas N, et al. mirPub: A database for searching microRNA publications. Bioinf, 2015; 31(9): 1502–1504. doi:10.1093/bioinformatics/btu819 38. Vlachos IS, Zagganas K, Paraskevopoulou MD, et al. DIANA-miRPath v3.0: Deciphering microRNA function with experimental support. Nu- cleic Acids Res, 2015; 43(Web Server issue): W460–W466. doi:10.1093/ nar/gkv403 39. Gentleman RC, Carey VJ, Bates DM, et al. Bioconductor: Open soft- 10 In-depth characterization... ware development for computational biology and bioinformatics. Ge- nome Biol, 2004; 5(10): R80. doi:10.1186/gb-2004-5-10-r80 40. Mao C, Zhou J, Yang Z, et al. KRAS, BRAF and PIK3CA mutations and the loss of PTEN expression in Chinese patients with colorectal cancer. PLoS ONE, 2012; 7(5): e36653. doi:10.1371/journal.pone.0036653 41. Kansal R, Quintanilla-Martinez L, Datta V, et al. Identification of the V600D mutation in Exon 15 of the BRAF oncogene in congenital, be- nign langerhans cell histiocytosis. Genes Chromosomes Cancer, 2013; 52(1): 99–106. doi:10.1002/gcc.22010 42. Yu L, Ma L, Tu Q, et al. Clinical significance of BRAF V600E mu- tation in 154 patients with thyroid nodules. Oncol Lett, 2015; 9(6): 2633–2638. doi:10.3892/ol.2015.3119 43. Yip L, Kelly L, Shuai Y, et al. MicroRNA signature distinguishes the degree of aggressiveness of papillary thyroid carcinoma. Ann Surg On- col, 2011; 18(7): 2035–2041. doi:10.1245/s10434-011-1733-0 44. Saiselet M, Gacquer D, Spinette A, et al. New global analysis of the microRNA transcriptome of primary tumors and lymph node me- tastases of papillary thyroid cancer. BMC Genomics, 2015; 16: 828. doi:10.1186/s12864-015-2082-3 45. Ahn D, Park JS, Sohn JH, et al. BRAFV600E mutation does not serve as a prognostic factor in Korean patients with papillary thyroid car- cinoma. Auris Nasus Larynx, 2012; 39(2): 198–203. doi:10.1016/j. anl.2011.07.011 46. Gandolfi G, Sancisi V, Torricelli F, et al. Allele percentage of the BRAF V600E mutation in papillary thyroid carcinomas and corresponding lymph node metastases: No evidence for a role in tumor progression. J Clin Endocrinol Metab, 2013; 98(5): E934–942. doi:10.1210/jc.2012- 3930 47. Givens DJ, Buchmann LO, Agarwal AM, et al. BRAF V600E does not predict aggressive features of pediatric papillary thyroid carcinoma. Laryngoscope, 2014; 124(9): E389–393. doi:10.1002/lary.24668 48. Cahill S, Smyth P, Denning K, et al. Effect of BRAFV600E muta- tion on transcription and post-transcriptional regulation in a papillary thyroid carcinoma model. Mol Cancer, 2007; 6: 21. doi:10.1186/1476- 4598-6-21 49. Fuziwara CS, Kimura ET. MicroRNA deregulation in anaplastic thy- roid cancer biology. Int J Endocrinol. doi:10.1155/2014/743450 50. Gao Y, Wang C, Shan Z, et al. miRNA expression in a human papillary thyroid carcinoma cell line varies with invasiveness. Endocr J, 2010; 57(1): 81–86. doi:10.1507/endocrj.k09e-220 51. Pecot CV, Rupaimoole R, Yang D, et al. Tumour angiogenesis regula- tion by the miR-200 family. Nat Commun, 2013; 4: 2427. doi:10.1038/ ncomms3427 52. Iozzo RV, Sanderson RD. Proteoglycans in cancer biology, tumour microenvironment and angiogenesis. J Cell Mol Med, 2011; 15(5): 1013–1031. doi:10.1111/j.1582-4934.2010.01236.x 53. Martínez-Aguilar J, Clifton-Bligh R, Molloy MP. Proteomics of thyroid tumours provides new insights into their molecular composition and changes associated with malignancy. Sci Rep, 2016; 6: 23660. doi:10.1038/srep23660 54. Nagarajan A, Malvi P, Wajapeyee N. Heparan sulfate and heparan sulfate proteoglycans in cancer initiation and progression. Front Endocrinol (Lausanne), 2018; 9. doi:10.3389/fendo.2018.00483 55. Mitselou A, Skoufi U, Tsimogiannis KE, et al. Association of Syn- decan-1 with angiogenesis-related markers, extracellular matrix components, and clinicopathological features in colorectal carci- noma. Anticancer Res, 2012; 32(9): 3977–3985. 56. Shishiba Y, Yanagishita M, Tanaka T, et al. Abnormal accumulation of proteoglycan in human thyroid adenocarcinoma tissue. Endocri- nol Jpn, 1984; 31(4): 501–507. doi:10.1507/endocrj1954.31.501 57. Arcinas A, Yen T-Y, Kebebew E, et al. Cell surface and secreted protein profiles of human thyroid cancer cell lines reveal distinct glycoprotein patterns. J Proteome Res, 2009; 8(8): 3958–3968. doi:10.1021/pr900278c 58. Maurel M, Jalvy S, Ladeiro Y, et al. A functional screening iden- tifies five microRNAs controlling glypican-3: Role of miR-1271 down-regulation in hepatocellular carcinoma. Hepatol, 2013; 57(1): 195–204. doi:10.1002/hep.25994 59. Hannafon BN, Sebastiani P, de las Morenas A, et al. Expression of microRNA and their gene targets are dysregulated in preinvasive breast cancer. Breast Cancer Res, 2011; 13(2): R24. doi:10.1186/ bcr2839 60. Ibrahim SA, Hassan H, Götte M. MicroRNA regulation of pro- teoglycan function in cancer. FEBS J, 2014; 281(22): 5009–5022. doi:10.1111/febs.13026 61. Marini F, Luzi E, Brandi ML. MicroRNA role in thyroid cancer development. J Thyroid Res. doi:10.4061/2011/407123 62. Qiu J, Zhang W, Zang C, et al. Identification of key genes and miR- NAs markers of papillary thyroid cancer. Biol Res, 2018; 51(1): 45. doi:10.1186/s40659-018-0188-1 63. Liu X, Tang Q, Chen H, et al. Lentiviral miR30-based RNA inter- ference against heparanase suppresses melanoma metastasis with lower liver and lung toxicity. Int J Biol Sci, 2013; 9(6): 564–577. doi:10.7150/ijbs.5425 64. Asuthkar S, Velpula KK, Nalla AK, et al. Irradiation-induced an- giogenesis is associated with an MMP-9-miR-494-syndecan-1 regulatory loop in medulloblastoma cells. Oncogene, 2014; 33(15): 1922–1933. doi:10.1038/onc.2013.151 65. Nucera C, Lawler J, Hodin R, et al. The BRAFV600E mutation: What is it really orchestrating in thyroid cancer? Oncotarget, 2010; 1(8): 751–756. 66. Geraldo MV, Kimura ET. Integrated analysis of thyroid can- cer public datasets reveals role of post-transcriptional regulation on tumor progression by targeting of immune system media- tors. PLoS ONE, 2015; 10(11): e0141726. doi:10.1371/journal.