Integrated Characterization of MicroRNA and mRNA Transcriptome in Papillary Thyroid Carcinoma

The incidence rate of papillary thyroid carcinoma (PTC) has rapidly increased in the recent decades, and the microRNA (miRNA) is one of the potential biomarkers in this cancer. Despite good prognosis, certain features such as lymph node metastasis (LNM) and BRAF V600E mutation are associated with a poor outcome. More than 50% of PTC patients present with LNM and BRAF V600E is the most common mutation identified in this cancer. The molecular mechanisms underlying these features are yet to be elucidated. This study aims to elucidate miRNA–genes interaction networks in PTC with or without LNM and to determine the association of BRAF V600E mutation with miRNAs and genes expression profiles. Next generation sequencing was performed to characterize miRNA and gene expression profiles in 20 fresh frozen tumor and the normal adjacent tissues of PTC with LNM positive (PTC LNM-P) and PTC without LNM (PTC LNN). BRAF V600E was genotyped using Sanger sequencing. Bioinformatics integration and pathway analysis were performed to determine the regulatory networks involved. Based on network analysis, we then investigated the association between miRNA and gene biomarkers, and pathway enrichment analysis was performed to study the role of candidate biomarkers. We identified 138 and 43 significantly deregulated miRNAs (adjusted p value < 0.05; log2 fold change ≤ −1.0 or ≥1.0) in PTC LNM-P and PTC LNN compared to adjacent normal tissues, respectively. Ninety-six miRNAs had significant expression ratios of 3p-to-5p in PTC LNM-P as compared to PTC LNN. In addition, ribosomal RNA-reduced RNA sequencing analysis revealed 699 significantly deregulated genes in PTC LNM-P versus normal adjacent tissues, 1,362 genes in PTC LNN versus normal adjacent tissue, and 1,576 genes in PTC LNM-P versus PTC LNN. We provide the evidence of miRNA and gene interactions, which are involved in LNM of papillary thyroid cancer. These findings may lead to better understanding of carcinogenesis and metastasis processes. This study also complements the existing knowledge about deregulated miRNAs in papillary thyroid carcinoma development.

Introduction

Papillary thyroid cancer (PTC) is the most common type of thyroid malignancy and it accounts for 80% of thyroid cancers; with approximately 77% diagnosed in women (1). Women are three to four times more likely to develop thyroid cancer than men, an opposite trend compared to other non-gender-specific sites (2). In the past few decades, its occurrence has increased significantly in many countries (3). Despite a high survival rate, patients with certain clinicopathologic features have been associated with poorer prognosis, such as older age at diagnosis (4), gender (5), extrathyroidal invasion (6), BRAF V600E mutation status (7), and lymph node metastasis (LNM) (8). In addition, the presence of LNM is associated with locoregional recurrence and increases the risk of mortality even in young patients (8). LNM is one of the factor of PTC propensity to metastasize to distant organs (9). Due to that, lymph node resection in PTC patients is considered vital (10).

The molecular mechanisms underlying LNM in PTC remain poorly understood. One of the molecules of interest in terms of elucidating the molecular landscape of cancer is microRNA (miRNA). miRNAs are small non-coding RNAs that comprise of 19 to 22 nucleotides, which negatively regulate gene expression. The relationship between miRNA with tumor growth, disease progression, and response to treatment has been demonstrated through miRNA expression profiling studies, which indicate the ability of these molecules to be used as biomarkers for cancer diagnosis and prognosis (11). miRNA expression is not only able to distinguish between cancer samples and normal samples from various types of cancer, but it can also distinguish whether the cancer cases have metastasized to the lymph nodes or not (12). In addition, the involvement of miRNA in the process of tumorigenesis and metastasis in PTC has also been demonstrated (13). Furthermore, the integration of miRNA and mRNA analysis can assist in detecting potential targets and may show the evidence to explain the importance of biological association of LNM in PTC (14).

The knowledge of how miRNAs interact with their target mRNAs and the consequences of this interaction will bring us one step closer to understanding the biological impact of these regulation on cell fate and function. Therefore, this study aims to elucidate miRNA–genes interaction networks in PTC with and without LNM and to associate BRAF V600E mutation with miRNAs and genes expression profiles.

Materials and Methods

Clinical Specimens

Ten pairs of tumor and adjacent normal fresh frozen tissues were collected from patients diagnosed with PTC with or without LNM (PTC LNM-P and PTC LNN, respectively) from the Universiti Kebangsaan Malaysia Medical Centre (UKMMC). This study was carried out in accordance with the recommendations of UKM Research Ethics Committee (UKMREC; UKM 1.5.3.5/244/UMBI-2015-002). Written informed consent was obtained from all the study participants. The tissues were dissected, snap-frozen, and stored in liquid nitrogen. All tissues were cryosectioned and stained with hematoxylin and eosin and the percentage of tumor and normal cells were assessed by a pathologist (Isa Mohamed Rose). Only cancer samples with at least 80% cancerous cells and normal adjacent thyroid tissues with less than 20% necrosis were selected for further analysis. The 10 paired fresh frozen specimens were subjected to miRNA and mRNA sequencing while additional 11 formalin-fixed, embedded tissues (FFPE) samples were used for qPCR validation.

Nucleic Acid Isolation

All fresh frozen tissues and FFPE specimens were subjected to nucleic acid extraction using Allprep DNA/RNA/miRNA Universal Kit or Allprep DNA/RNA FFPE Kit (Qiagen, Germany), respectively according to manufacturer’s recommendations. The integrity of RNA was assessed using Agilent RNA 6000 Nano Kit and only total RNA with RNA Integrity Number (RIN) of at least six were used for subsequent steps. The miRNA content was quantified using Small RNA Kit (both from Agilent Technologies, Palo Alto, CA, USA) 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 (Invitrogen, USA) and Nanodrop 2000c Spectrometer (Thermo Scientific/Nanodrop, DE), respectively. The acceptable purity as indicated by A260/280 is 2.0–2.2 for the RNA and 1.8–2.0 for the DNA.

Library Preparation

microRNA libraries were prepared using the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, 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, normalized to 2 Nm, and pooled for multiplexed sequencing.

Meanwhile, whole transcriptome libraries were prepared using Ion Total RNAseq v2 kit (Life Technologies, USA) according to manufacturer’s protocol. Two micrograms of total RNA were used for ribosomal removal using RiboMinus Eukaryote kit (Life Technologies, USA). ERCC RNA Spike-In Control Mixes (Thermo Scientific, USA) were added into rRNA depleted-RNA before RNA fragmentation using RNase III. rRNA-depleted and fragmented RNA were quantified using Qubit fluorometer assay and Agilent RNA 6000 Nano Kit. Subsequently, the fragmented RNA was hybridized and ligated to adaptors, reverse-transcribed, purified, size-selected, and amplified. The final libraries were quantified using Agilent High Sensitivity DNA chip (Agilent Technologies, USA) and were normalized to 100–130 pM for the next step.

miRNA and mRNA Sequencing

For miRNA sequencing, the resulting pooled libraries were hybridized to oligonucleotide-coated single-read flow cells for cluster generation using HiSeq^® Rapid SR Cluster Kit v2 on Hiseq 2500 (Illumina, USA). The clustered pooled miRNA libraries were then sequenced on the HiSeq 2500 for 50 sequencing cycles using HiSeq^® Rapid SBS Kit v2 (50 Cycle) (Illumina, USA). As for mRNA sequencing, each library were subjected to clonal amplification on Ion Chef System using Ion PI IC 200 kit followed by sequencing using PI BC v2 chip on Ion Proton system (all from Life Technologies, USA). Two to three mRNA libraries were pooled into a single chip.

Bioinformatics, Integrative, and Statistical Analyses

Preprocessing of miRNA sequencing data was executed in BaseSpace software (Illumina, USA), and FASTQ files were generated. MiRNA Analysis app version 1.0.0 was used for determination of differentially expressed (DE) miRNAs using the workflow described by Cordero et al. (15). Briefly, the pipeline includes 3′ end adapter removal using cutadapt (16), annotation to miRBase v21 (17), mapping using SHRIMP aligner (18), and differential analysis of miRNAs using DESeq2 (19). In addition, Rank Product statistics were used to detect 3p/5p ratio changes on the same miRNA in two different conditions (20). Unsupervised hierarchical clustering and heatmaps were performed and created using Morpheus from Broad Institute (https://software.broadinstitute.org/morpheus/).Pathway enrichment analysis were executed using DIANA-miRPath v3.0 (21) and significance was determined using Fisher’s Exact Test.

For mRNA sequencing, preprocessing of data were conducted on Torrent Server using Torrent Suite v4.4.2 software (Life Technologies, USA) and FASTQ files were generated and then exported for data analysis using CLCBio Genomics Workbench v8.3 (CLC Bio, Denmark). The manufacturer’s analysis pipeline for transcriptome sequencing was used for differential expression of mRNA. Briefly, the pipeline includes removal of 3′ adapter, mapping to reference transcript (hg19) and differential expression analysis using the Empirical analysis of DGE (EDGE) (22). Only DE miRNAs or genes with p-value <0.05 and log2-fold change ≤−1.0 or ≥1.0 were subjected to further analysis.

Integrated miRNA and mRNA analysis were performed using MAGIA2 (23). Log2 normalized reads from the significant DE miRNAs and mRNAs were uploaded into MAGIA2 and positive and negative correlation analyses between miRNAs, transcription factors, and target mRNA were performed. Skip variability filter option was selected because the data uploaded only contain the significant miRNAs or mRNAs. miRNA targets were predicted using DIANA microT (24) (score threshold at 2.0, top 75% of the predictions distribution) and Spearman correlation coefficient was used as correlations measure (threshold <0.05). In addition, prediction of miRNA-transcription factor was also performed based on mirGen v2.0 (25) and TransmiR (26) databases.

BRAF V600E Genotyping

PCR amplification of genomic regions of interest was performed using BRAF V600E forward primer 5′-TGCTTGCTCTGATAGGAAAATG-3′ and reverse primer 5′-AGCATCTCAGGGCCAAAAAT-3′ (27). Amplification was performed in a reaction volume of 25 μl containing 50 ng DNA template, 10 μM each primer (Integrated DNA Technologies, USA), 10× PCR Gold Buffer without MgCl₂, dNTP Mix (10 mM), MgCl₂ solution (25 mM), AmpliTaq Gold^® (5 U/μl) (Applied Biosystems, USA) and nuclease-free water. PCR conditions were as follows; 95°C for 4 min; 35 cycles of 95°C for 45 s, 50°C for 30 s, and 72°C for 1 min; 72°C for 5 min; 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, Germany) as per manufacturer’s instruction. Subsequently, DNA sequencing was performed using ABI Prism 3130xl Genetic Analyzer (Applied Biosystem, USA).

Real-Time PCR Validation

Based on 3p/5p miRNA expression ratio analysis, six most significant DE miRNAs were chosen for validation. The Pick-&-Mix panel (Exiqon, Denmark) consisting of hsa-miR-193a-3p, hsa-miR-193a-5p, hsa-miR-376a-3p, hsa-miR-376a-5p, hsa-miR-205-3p, and hsa-miR-205-5p primers was used on the set of samples with hsa-miR-10b-5p and hsa-miR-191-5p as reference miRNA. This assay was carried out on CFX96 Touch Real-Time PCR Detection (BioRad, USA). Fold change was calculated using 2^−ΔΔCq (28) and statistical correlation between miRNA sequencing and qPCR was analyzed using R² on GraphPad Prism v7.0 (Graphpad Software, USA).

Results

Demographic Data

As summarized in Table 1, all subjects in our discovery phase were female patients and 70% of them were Malays. Mean age at diagnosis was 51.6 years in the PTC LNM-P group and 47.4 years in the PTC LNN group. Majority of the cancers were located at the right lobe. All subjects with PTC LNM-P had BRAF V600E mutation.

TABLE 1

Table 1. Clinicopathological features of patients in discovery and validation phase.

DE miRNAs

We achieved an average of 5.6 M reads per sample with Q30 (3,779, 968–8,308, 285 reads). miRNA expression analysis was divided into three comparisons, which were PTC LNM-P versus adjacent normal, PTC LNN versus adjacent normal, and PTC LNM-P versus PTC LNN. A total of 138 miRNAs were significantly DE (66 upregulated and 72 downregulated) in PTC LNM-P versus their normal adjacent group. Whereas for PTC LNN versus adjacent normal, 43 miRNAs were significantly DE (31 upregulated and 12 downregulated). However, no miRNA was DE when we compared PTC LNM-P versus PTC LNN. Unsupervised hierarchical clustering revealed that all samples were clustered according to their respective groups (Figures 1 and 2), proving that the miRNAs expression signatures were able to differentiate tumors from normal thyroid tissues even when the tissues originated from the same patients. The list of significant DE miRNAs was provided in Tables S1 and S2 in Supplementary Material. The miRNA sequencing data in fastq. format was deposited at the NCBI Sequence Read Archive (SRA) at http://www.ncbi.nlm.nih.gov/sra with submission IDs SUB2458655 and SUB2753230.

FIGURE 1

Figure 1. Hierarchical clustering and heatmap for microRNA (miRNA) expression profile in papillary thyroid carcinoma LNM-P versus adjacent normal comparison group. The list of miRNA were filtered with p-value <0.05 and log2-fold change ≤−1.0 or ≥1.0. A total of 138 miRNA have significant differentially expression. The sample labeled in yellow bar represents tumor sample while green bar represents normal samples. Color key illustrate the relative expression of miRNA in all samples. Red indicates high expression level while blue indicate low expression level. Dendrogram shows all samples were clustered well according to their respective groups.

FIGURE 2

Figure 2. Hierarchical clustering and heatmap for microRNA (miRNA) expression profile in papillary thyroid carcinoma LNN versus adjacent normal comparison group. The list of miRNA were filtered with p value <0.05 and log2-fold change ≤−1.0 or ≥1.0. A total of 43 miRNA have significant differentially expression. The sample labeled in yellow bar represents tumor sample while green bar represents normal samples. Color key illustrate the relative expression of miRNA in all samples. Red indicates high expression level while blue indicate low expression level. Dendrogram showed that all samples were clustered well according to their respective groups.

Pathway Enrichment Analysis of the Deregulated miRNAs

A single miRNA may play a specific role in the pathogenesis of PTC; however, the miRNA pathway as a whole may also be of importance. Therefore, the DIANA-mirPath v3.0 (21) was used to identify pathways that are potentially regulated by the expression of multiple miRNAs and to incorporate miRNAs into molecular pathways using experimentally validated miRNA target according to TarBase issue v7 (29). Using the upregulated miRNAs in PTC LNM-P, 47 pathways were enriched significantly (Figure 3). On the other hand, using downregulated miRNAs in PTC LNM-P, four pathways were enriched significantly, which include fatty acid biosynthesis (p < 1E−325), fatty acid metabolism (p < 1E−325), fatty acid elongation (p = 9.30E−10), and fatty acid degradation (p = 3.60E−02).

FIGURE 3

Figure 3. Pathway enrichment analysis on 66 upregulated microRNA (miRNA) in papillary thyroid carcinoma LNM-P and heatmap for miRNA versus pathway (clustered based on their significance value) using DIANA miRPath v3.0 software. Dendrogram on both axis represent hierarchical clustering resulting from miRNA and pathway. The miRNA were clustered together with similar pathway pattern and the pathways were clustered together based on related miRNAs.

For miRNAs, which were upregulated in PTC LNN, 30 pathways were significantly enriched (Figure 4) while eight pathways were enriched significantly using downregulated miRNAs, which include fatty acid biosynthesis (p < 1E−325), fatty acid metabolism (p = 1.31E−11), lysine degradation (p = 1.35E−11), proteoglycans in cancer (p = 9.75E−04), chronic myeloid leukemia (p = 3.55E−03), viral carcinogenesis (p = 8.51E−03), glioma (p = 2.13E−02), and pathways in cancer (p = 2.48E−02).

FIGURE 4

Figure 4. Pathway enrichment analysis on 31 upregulated microRNA (miRNA) in papillary thyroid carcinoma LNN and heatmap for miRNA versus pathway (clustered based on their significance value) using DIANA miRPath v3.0 software. Dendrogram on both axis represent hierarchical clustering resulting from miRNA and pathway. The miRNA were clustered together with similar pathway pattern and the pathways were clustered together based on related miRNAs.

Ratio of 3p/5p miRNA

In order to analyze the 3p and 5p expression ratio on the same miRNA in two different groups, PTC LNM-P and PTC LNN tumor sequencing data were analyzed using RankProd (20). Ninety-six miRNA had significant differential expression ratio and present on both group. Table 2 and Figure 5 show the list of top 20 miRNAs with significant 3p/5p ratio.

TABLE 2

Table 2. The list of top rank 20 microRNA (miRNA) with highest expression ratio.

FIGURE 5

Figure 5. The top 20 3p/5p expression ratio in both comparison groups. Several microRNA (miRNA) have significant expression ratio shows that arm selection preference in miRNA are different significantly between papillary thyroid carcinoma (PTC) LNM-P and PTC LNN.

DE Genes

Analysis of gene expression was also divided into three groups, which are PTC LNM-P versus adjacent normal, PTC LNN versus adjacent normal, and PTC LNM-P versus PTC LNN. In PTC LNM-P versus adjacent normal, 348 genes were upregulated and 351 were downregulated. As for PTC LNN versus adjacent normal group, 396 genes were upregulated and 966 were downregulated. For PTC LNM-P versus PTC LNN, 733 genes were upregulated and 843 were downregulated. Tables S3–S5 in Supplementary Material summarized the list of significant DE genes. The mRNA sequencing data in fastq. format was deposited at the NCBI SRA at http://www.ncbi.nlm.nih.gov/sra with submission ID SUB2756288.

Integrated Analysis of miRNA and Gene Expression

In order to obtain a comprehensive understanding of regulatory molecules in PTC LNM-P and PTC LNN, integrated analysis of miRNA-gene target was performed using MAGIA2 (23). In PTC LNM-P versus adjacent normal, integrated analysis revealed 47 miRNAs in 295 significant interactions while PTC LNN demonstrated 19 miRNAs in 70 significant interactions. Figures 6 and 7 illustrated 200 top significant miRNA–gene interactions in PTC LNM-P and PTC LNN, respectively. Interestingly, interactions that involved hsa-let-7i-5p and hsa-let-7f-5p only presented in PTC LNM-P.

FIGURE 6

Figure 6. Grand view for 200 main significant microRNA (miRNA)–gene interaction that was generated using miRNA and gene that have significant expression in papillary thyroid carcinoma LNM-P versus their adjacent normal group.

FIGURE 7

Figure 7. Grand view for 200 main significant microRNA (miRNA)–gene interaction that was generated using miRNA and gene that have significant expression in papillary thyroid carcinoma LNN versus their adjacent normal group.

Validation of 3p/5p miRNA Expression Ratio

Expression ratio of selected miRNAs were validated in 21 patients and all were in concordance with miRNA sequencing (R² = 0.9998) (Figure 8).

FIGURE 8

Figure 8. Correlation coefficient analyses for comparison between microRNA (miRNA) sequencing (miRNAseq) and qPCR validation of hsa-miR-193a-3p/5p, hsa-miR-376a-3p/5p and hsa-miR-205-3p/5p. The two technologies produced concordance result with R² = 0.9998, indicating the positive correlation of fold change of six selected miRNAs qPCR and miRNAseq.

Discussion

We have identified 138 significant DE miRNAs in PTC LNM-P compared to normal adjacent thyroid tissues; 66 miRNAs were upregulated while 72 miRNAs were downregulated. On the other hand, comparison between PTC LNN and normal adjacent thyroid revealed 43 significant DE miRNAs (31 upregulated and 12 downregulated). A few DE miRNAs from our findings such as hsa-miR-146b, hsa-miR-222, hsa-miR-221, hsa-miR-204, and hsa-miR-7 were commonly dysregulated in both comparisons (PTC LNM-P versus normal adjacent thyroid and PTC LNN versus normal adjacent thyroid) and have consistent expression levels with several other studies (30–33). Hsa-miR-146-5p and 3p in particular are the top upregulated miRNA in PTC compared to normal thyroid and were also highly expressed in our PTC LNM-P and PTC LNN versus respective normal thyroid, suggesting their roles in thyroid carcinogenesis in general. Yet, these miRNAs expression were not significant when we compared PTC LNM-P to PTC LNN. Even though there were no DE miRNAs in PTC LNM-P versus PTC LNN group, we did identified 733 upregulated and 843 downregulated genes. Tables S3–S5 in Supplementary Material summarized the list of significant DE genes. Pathway enrichment analyses of PTC LNM-P revealed alteration in mitogen-activated protein kinase (MAPK) signaling pathway (Figure 3), which could be due to presence of BRAF V600E mutation. On the contrary, this pathway was not enriched in PTC LNN analysis (Figure 4), which does not have BRAF V600E mutation. Our finding is concordant with Wang et al. (34) but is in disagreement with several others (35, 36). On closer inspection, the upregulation was only significant with a larger sample size with a modest fold change. Through functional assays, hsa-miR-146b-5p was shown to promote migration, invasion, and induced epithelial-to-mesenchymal (EMT) by downregulating ZNRF3 and modulating Wnt/β-catenin signaling. However, the function of the passenger strand of miR-146b, hsa-miR-146b-3p, is yet to be elucidated.

Several studies have shown that the 3p/5p miRNA arm selection were variable among human cancers (36, 37). To date, the expression ratio between 3p and 5p miRNA is not well-characterized and there has been only a few studies that investigated the arm selection preference in relation to cancer. Even though 5p and 3p miRNA were derived from the same pre-miRNA, both arms might have different expression in the same tissues where one arm is highly expressed in tumor tissues than normal tissues, or vice versa (36). Our study did not detect any significant DE miRNAs in PTC LNM-P as compared to PTC LNN. However, when we examined 3p and 5p miRNA expression ratio in PTC LNM-P versus PTC LNN, 96 miRNAs have significant expression ratio changes. Among the miRNAs that have significant expression ratios in this study include hsa-miR-205-3p/5p, hsa-miR-193a-3p/5p, and hsa-miR-367a-3p/5p. For example, hsa-miR-205-5p was more highly expressed than hsa-miR-205-3p in both PTC LNM-P and PTC LNN; however, the 3p/5p expression ratio increased from −9.14 in PTC LNN to −6.95 in PTC LNM-P, which is a 2.18 log-fold increment. This suggest that hsa-miR-205-5p was preferentially downregulated in PTC LNM-P while hsa-miR-205-3p expression remained the same. This preference during miRNA biogenesis process could potentially play a role in LNM. To the best of our knowledge, there is only one study previously conducted on 3p and 5p miRNA expression ratio preference in PTC and unfortunately, is inconsistent with our findings (13). Saiselet et al. (13) performed small RNA deep sequencing of only three PTC and compared the miRNA profiles with the matching normal tissues and metastatic lymph node from the same patients. The authors then suggested the lack of involvement of 3p/5p miRNA arm selection in PTC tumorigenesis and LNM progression (13). On the other hands, our analyses involved PTC LNM-P compared to PTC LNN; which were from different patient groups. The inconsistency with study by Saiselet and colleagues (13) could be attributed by different experimental design and this aspect of looking at miRNAs from a different angle certainly warrant further investigations.

Hsa-miR-205 is a regulatory miRNA where its dysregulation can lead to cancer formation, EMT, and metastasis (38). Up until now, there is no study conducted on hsa-miR-205-3p/5p expression ratio in PTC or thyroid cancer in general. However, this miRNA plays two different roles, namely as an oncogene in lung cancer, bladder cancer, cervical cancer and head and neck cancer, and as a tumor suppressor in breast cancer, colorectal cancer and glioma (39, 40). Through our miRNA profiling, hsa-miR-205 was upregulated in PTC LNM-P and PTC LNN when compared to their adjacent normal tissues, respectively. However, only hsa-miR-205-5p reached statistical significance in PTC LNN versus adjacent normal tissues. Our findings were consistent with the study conducted by Nikiforova et al. in 2008 where hsa-miR-205 was upregulated in PTC and anaplastic thyroid cancer as compared to hyperplastic nodule and normal thyroid tissues (41). However, downregulation of hsa-miR-205 was identified to be significant in metastatic PTC as compared to PTC without metastasis. This downregulated expression was associated with upregulation of VEGF in thyroid tumor tissue, hence it proved the role of miR-205 as tumor suppressor in a variety of thyroid cancer cell lines (42). These discrepancies and the limited references currently available highlight the needs to further investigate the exact role of miR-205 in PTC and LNM.

In order to identify the role of miRNAs and the pathways involved, pathways enrichment analyses were carried out. Among significant pathways resulted from these analyses were proteoglycans in cancer, ECM–receptor interaction and MAPK signaling pathway. There are 50 genes and 21 miRNAs in PTC LNM-P and 22 genes and 4 miRNAs in PTC LNN that are involved in ECM–receptor interaction pathway. Cellular matrix (ECM) plays an important role in homeostasis development and maintenance as well as tissue and organ architectures. The synthesis and degradation of ECM components such as collagen type I, collagen type IV, and fibronectin can cause modification, proliferation activation, migration, and activated adhesive endothelium cell to ECM and encourage angiogenesis process to occur (43). In a study conducted by Zhang et al. (44), expression of fibronectin-1 was significantly higher in the PTC cell line that was transfected with CXCR7 and ECM-receptor interaction was one of the enriched pathway resulted. As an important ECM component, FN1 plays a role in maintaining the survival, proliferation, adhesion, migration, and angiogenesis of epithelial cells (45). In addition, several studies reported that expression of FN1 was high in PTC and it can act as a potential biomarker in detecting PTC (46, 47). Many other studies have also proven the involvement of ECM in PTC progression (47, 48).

To obtain an understanding of the target mRNA expression that have significant expression, we integrated the miRNA and mRNA sequencing data. When the miRNA acts through the degradation of target mRNA, the miRNA and mRNA target expression profiles are expected to be inversely correlated (49). However, a positive correlation (miRNA highly expressed/mRNA highly expressed or miRNA lowly expressed/mRNA lowly expressed) is formed when the activity of miRNA is part of a series of complex regulatory and gene expression profiles resulting from different levels of regulation (49). There are two miRNAs that are involved in both comparison groups, namely, hsa-miR-7-5p and hsa-miR-181b-5p, hence suggesting that both miRNAs are involved in the pathogenesis of PTC in general and are not specific to LNM. Various studies have reported that hsa-miR-7-5p has a low expression level in PTC samples (13, 31, 50). Hsa-miR-7 plays a role in regulating cell growth, migration, and invasion by targeting PAK1 in thyroid cancer (51). The low expression level of hsa-miR-181b can cause cell growth inhibition and apoptosis by targeting CYLD, and this can be harnessed as a potential therapeutic target for the treatment of patients with PTC (52).

Let-7 (lethal-7) is an important gene involved in development in C. elegans, as well as being one of the first miRNA that was discovered (53). Let 7i-5p and let-7f-5p formed a unique miRNA-mRNA network that was only detected in PTC LNM-P versus normal adjacent tissues when the miRNA and mRNA sequencing data were integrated. However, the expression profiling in this study was inconsistent with one previous study (54). The expression of let-7f was not significant in cancer tissues with BRAF V600E. However, let-7f was highly expressed in tumor tissues with BRAF-wildtype (54). In a study involving breast cancer by Liu et al. (55), let-7f regulated the expression of β2-AR in breast cancer cells. In human breast cancer cell lines (MCF-7, SKBR3, and BT474), let-7f caused an increasing expression of β2-AR and was also shown to be associated with LNM (55). In addition, let-7i was downregulated in PTC with LNM-P compared to PTC LNN (56), and it is upregulated in primary PTC when compared with benign thyroid nodules and normal tissues (39). However, there were no significant changes in the expression when the results were validated using qPCR technique (22).

The T1799A nucleotide transversion in BRAF gene is an oncogenic mutation that is widely reported in PTC, and it occurs in an average of 45% of the cases (57). This mutation causes the change from glutamic acid to valine at codon 600 BRAF protein, resulting in the BRAF V600E variant. BRAF V600E constitutively increases the activity of serine/threonine protein and activates MAPK signaling pathways in human cancers (57). PTC patients with BRAF V600E are often linked with the aggressive clinicopathological characteristics such as LNM (58, 59). In this study, 11 out of 12 PTC LNM-P patients showed the presence of BRAF V600E, while it was detected in only two of the nine patients in the PTC LNN samples. From our study, we postulated that the BRAF V600E mutation is associated with LNM because almost all (11/12) of PTC LNM-P patients have this mutation. The study conducted by The Cancer Genome Atlas regarding the classification of PTC showed that miRNA expression is associated with the PTC phenotypes and plays an important role in the prognosis of PTC (30). Although definitive studies are still limited, miRNA clearly has the potential as an indicator of prognosis in PTC patients (60). In another study conducted by Cahill et al., which compared a PTC cell line with BRAF mutation and a normal thyroid cell line, the analysis revealed 15 highly expressed and 23 downregulated miRNAs (61). In addition, human thyroid cancers that are BRAF mutated or have other mutations have different miRNA expression profiles (41). Although the results of our study are consistent with most of other studies (58, 62), there are also other studies, which indicate that BRAF V600E mutation do not contribute to the aggressive clinicopathological characteristics in PTC patients (13, 63). These discrepancies could be due to the unique variation in each PTC patient.

Conclusion

To the best of our knowledge, this is the first study that involved the integrated analysis of miRNA and gene expression profiles from thyroid cancer and the corresponding normal adjacent thyroid tissues in the same patients and associate the expression signatures with LNM and BRAF V600E status. PTC LNM-P exhibited a higher number of dysregulated miRNAs and affected pathways compared to PTC LNN. Even though there was no significant DE miRNAs in PTC LNM-P versus PTC LNN, there were significant changes in the 3p and 5p expression ratios. This study proposed the miRNA-gene expression patterns and miRNA 3p/5p ratios that may become potentially useful to identify PTC patients with LNM. Further investigations, such as in vitro functional studies are warranted to fully understand the role of these miRNAs in LNM in this cancer.

Ethics Statement

This study carried out in accordance with the recommendations of UKM Research Ethics Committee (UKMREC; UKM 1.5.3.5/244/UMBI-2015-002). Written informed consent was obtained from all the study participants.

Author Contributions

AY and N-SM involved in the specimen collections, libraries preparations and sequencings, data analyses, acquisition of data, and drafting the manuscript. RM and SS are thyroid surgeons involved in specimen retrieval. IR assessed tumor percentage of the tissues. SS performed the BRAF V600E genotyping. RJ provides critical review on the manuscript. All authors read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors acknowledge the contributions of PPUKM-Biobank Team, UMBI Cancer Genomics Team, and Universiti Kebangsaan Malaysia for the supports.

Funding

This research was funded by Fundamental Research Grant Scheme (FRGS) from the Ministry of Higher Education Malaysia (FRGS/1/2014/SKK01/UKM/03/1).

Supplementary Material

The Supplementary Material for this article can be found online at https://www.frontiersin.org/articles/10.3389/fendo.2018.00158/full#supplementary-material.

References

1. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer (2015) 136:E359–86. doi:10.1002/ijc.29210

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rahbari R, Zhang L, Kebebew E. Thyroid cancer gender disparity. Future Oncol (2010) 6:1771–9. doi:10.2217/fon.10.127

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kitahara CM, Sosa JA. The changing incidence of thyroid cancer. Nat Rev Endocrinol (2016) 12:646–53. doi:10.1038/nrendo.2016.110

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Shi R, Qu N, Liao T, Wei WJ, Wang YL, Ji QH. The trend of age-group effect on prognosis in differentiated thyroid cancer. Sci Rep (2016) 6:27086. doi:10.1038/srep27086

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Jonklaas J, Nogueras-Gonzalez G, Munsell M, Litofsky D, Ain KB, Bigos ST, et al. The impact of age and gender on papillary thyroid cancer survival. J Clin Endocrinol Metab (2012) 97:E878–87. doi:10.1210/jc.2011-2864

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Yin DT, Yu K, Lu RQ, Li X, Xu J, Lei M. Prognostic impact of minimal extrathyroidal extension in papillary thyroid carcinoma. Medicine (2016) 95:e5794. doi:10.1097/MD.0000000000005794

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Czarniecka A, Kowal M, Rusinek D, Krajewska J, Jarzab M, Stobiecka E, et al. The risk of relapse in papillary thyroid cancer (PTC) in the context of BRAF V600E mutation status and other prognostic factors. PLoS One (2015) 10:e0132821. doi:10.1371/journal.pone.0132821

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Adam MA, Pura J, Goffredo P, Dinan MA, Reed SD, Scheri RP, et al. Presence and number of lymph node metastases are associated with compromised survival for patients younger than age 45 years with papillary thyroid cancer. J Clin Oncol (2015) 33:2370–5. doi:10.1200/JCO.2014.59.8391

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Schneider DF, Chen H, Sippel RS. Impact of lymph node ratio on survival in papillary thyroid cancer. Ann Surg Oncol (2013) 20:1906–11. doi:10.1245/s10434-012-2802-8

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Moo TA, Fahey TJ. Lymph node dissection in papillary thyroid carcinoma. Semin Nucl Med (2011) 41:84–8. doi:10.1053/j.semnuclmed.2010.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Iorio MV, Croce CM. microRNA involvement in human cancer. Carcinogenesis (2012) 33:1126–33. doi:10.1093/carcin/bgs140

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Yan LX, Huang XF, Shao Q, Huang MY, Deng L, Wu QL, et al. MicroRNA miR-21 overexpression in human breast cancer is associated with advanced clinical stage, lymph node metastasis and patient poor prognosis. RNA (2008) 14:2348–60. doi:10.1261/rna.1034808

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Saiselet M, Gacquer D, Spinette A, Craciun L, Decaussin-Petrucci M, Andry G, et al. New global analysis of the microRNA transcriptome of primary tumors and lymph node metastases of papillary thyroid cancer. BMC Genomics (2015) 16:828. doi:10.1186/s12864-015-2082-3

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Akhtar MM, Micolucci L, Islam MS, Olivieri F, Procopio AD. Bioinformatic tools for microRNA dissection. Nucleic Acids Res (2016) 44:24–44. doi:10.1093/nar/gkv1221

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Cordero F, Beccuti M, Arigoni M, Donatelli S, Calogero RA. Optimizing a massive parallel sequencing workflow for quantitative miRNA expression analysis. PLoS One (2012) 7:e31630. doi:10.1371/journal.pone.0031630

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J (2011) 17:10–2. doi:10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

17. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res (2014) 42:D68–73. doi:10.1093/nar/gkt1181

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Rumble SM, Lacroute P, Dalca AV, Fiume M, Sidow A, Brudno M. SHRiMP: accurate mapping of short color-space reads. PLoS Comput Biol (2009) 5:e1000386. doi:10.1371/journal.pcbi.1000386

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol (2010) 11:R106. doi:10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Breitling R, Armengaud P, Amtmann A, Herzyk P. Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments. FEBS Lett (2004) 573:83–92. doi:10.1016/j.febslet.2004.07.055

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Vlachos IS, Zagganas K, Paraskevopoulou MD, Georgakilas G, Karagkouni D, Vergoulis T, et al. DIANA-miRPath v3. 0: deciphering microRNA function with experimental support. Nucleic Acids Res (2015) 43:W460–6. doi:10.1093/nar/gkv403

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Robinson MD, Smyth GK. Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics (2007) 9:321–32. doi:10.1093/biostatistics/kxm030

CrossRef Full Text | Google Scholar

23. Bisognin A, Sales G, Coppe A, Bortoluzzi S, Romualdi C. MAGIA(2): from miRNA and genes expression data integrative analysis to microRNA-transcription factor mixed regulatory circuits (2012 update). Nucleic Acids Res (2012) 40:W13–21. doi:10.1093/nar/gks460

CrossRef Full Text | Google Scholar

24. Maragkakis M, Alexiou P, Papadopoulos GL, Reczko M, Dalamagas T, Giannopoulos G, et al. Accurate microRNA target prediction correlates with protein repression levels. BMC Bioinformatics (2009) 10:295. doi:10.1186/1471-2105-10-295

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Georgakilas G, Vlachos IS, Zagganas K, Vergoulis T, Paraskevopoulou MD, Kanellos I, et al. DIANA-miRGen v3.0: accurate characterization of microRNA promoters and their regulators. Nucleic Acids Res (2016) 44(D1):D190–5. doi:10.1093/nar/gkv1254

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wang J, Lu M, Qiu C, Cui Q. TransmiR: a transcription factor–microRNA regulation database. Nucleic Acids Res (2010) 38:D119–22. doi:10.1093/nar/gkp803

CrossRef Full Text | Google Scholar

27. Monti E, Bovero M, Mortara L, Pera G, Zupo S, Gugiatti E, et al. BRAF mutations in an Italian regional population: implications for the therapy of thyroid cancer. Int J Endocrinol (2015) 2015:138734. doi:10.1155/2015/138734

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2^{− ΔΔCT} method. Methods (2012) 25:402–8. doi:10.1006/meth.2001.1262

CrossRef Full Text | Google Scholar

29. Vlachos IS, Paraskevopoulou MD, Karagkouni D, Georgakilas G, Vergoulis T, Kanellos, et al. DIANA-TarBase v7. 0: indexing more than half a million experimentally supported miRNA: mRNA interactions. Nucleic Acids Res (2014) 43:D153–9. doi:10.1093/nar/gku1215

CrossRef Full Text | Google Scholar

30. Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell (2014) 159:676–90. doi:10.1016/j.cell.2014.09.050

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Mancikova V, Castelblanco E, Pineiro-Yanez E, Perales-Paton J, de Cubas AA, Inglada-Perez L, et al. MicroRNA deep-sequencing reveals master regulators of follicular and papillary thyroid tumors. Mod Pathol (2015) 28:748–57. doi:10.1038/modpathol.2015.44

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yoruker EE, Terzioglu D, Teksoz S, Uslu FE, Gezer U, Dalay N. MicroRNA expression profiles in papillary thyroid carcinoma, benign thyroid nodules and healthy controls. J Cancer (2016) 7:803. doi:10.7150/jca.13898

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Ab Mutalib NS, Yusof AM, Mokhtar NM, Harun R, Muhammad R, Jamal R. MicroRNAs and lymph node metastasis in papillary thyroid cancers. Asian Pac J Cancer Prev (2016) 17:25–35. doi:10.7314/APJCP.2016.17.1.25

CrossRef Full Text | Google Scholar

34. Deng X, Wu B, Xiao K, Kang J, Xie J, Zhang X, et al. MiR-146b-5p promotes metastasis and induces epithelial-mesenchymal transition in thyroid cancer by targeting ZNRF3. Cell Physiol Biochem (2015) 35:71–82. doi:10.1159/000369676

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Yang Z, Yuan Z, Fan Y, Deng X, Zheng Q. Integrated analyses of microRNA and mRNA expression profiles in aggressive papillary thyroid carcinoma. Mol Med Rep (2013) 8:1353–8. doi:10.3892/mmr.2013.1699

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Li SC, Liao YL, Ho MR, Tsai KW, Lai CH, Lin WC. miRNA arm selection and isomiR distribution in gastric cancer. BMC Genomics (2012) 13(Suppl 1):S13. doi:10.1186/1471-2164-13-S1-S13

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Tsai KW, Leung CM, Lo YH, Chen TW, Chan WC, Yu SY, et al. Arm selection preference of microRNA-193a varies in breast cancer. Sci Rep (2016) 6:28176. doi:10.1038/srep28176

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Vosgha H, Salajegheh A, Smith RA, Lam AK. The important roles of miR-205 in normal physiology, cancers and as a potential therapeutic target. Curr Cancer Drug Targets (2014) 14:621–37. doi:10.2174/156800961407140926105634

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li MH, Fu SB, Xiao HS. Genome-wide analysis of microRNA and mRNA expression signatures in cancer. Acta Pharmacol Sin (2015) 36:1200–11. doi:10.1038/aps.2015.67

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Yue X, Wang P, Xu J, Zhu Y, Sun G, Pang Q, et al. MicroRNA-205 functions as a tumor suppressor in human glioblastoma cells by targeting VEGF-A. Oncol Rep (2012) 27:1200–6. doi:10.3892/or.2011.1588

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Nikiforova MN, Tseng GC, Steward D, Diorio D, Nikiforov YE. MicroRNA expression profiling of thyroid tumors: biological significance and diagnostic utility. J Clin Endocrinol Metab (2008) 93:1600–8. doi:10.1210/jc.2007-2696

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Salajegheh A, Vosgha H, Rahman AM, Amin M, Smith RA, Lam AK. Modulatory role of miR-205 in angiogenesis and progression of thyroid cancer. J Mol Endocrinol (2015) 55:183–96. doi:10.1530/JME-15-0182

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Wang G, Wang L, Sun S, Wu J, Wang Q. Quantitative measurement of serum microRNA-21 expression in relation to breast cancer metastasis in Chinese females. Ann Lab Med (2015) 35:226–32. doi:10.3343/alm.2015.35.2.226

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang X, Dong H, Tian Y. MicroRNA Detection and Pathological Functions. Berlin: Springer-Verlag (2015).

Google Scholar

45. Paik JY, Ko BH, Jung KH, Lee KH. Fibronectin stimulates endothelial cell 18F-FDG uptake through focal adhesion kinase–mediated phosphatidylinositol 3-kinase/akt signaling. J Nucl Med (2009) 50:618–24. doi:10.2967/jnumed.108.059386

CrossRef Full Text | Google Scholar

46. Qu T, Li YP, Li XH, Chen Y. Identification of potential biomarkers and drugs for papillary thyroid cancer based on gene expression profile analysis. Mol Med Rep (2016) 14:5041–8. doi:10.3892/mmr.2016.5855

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhang S, Wang Y, Chen M, Sun L, Han J, Elena VK, et al. CXCL12 methylation-mediated epigenetic regulation of gene expression in papillary thyroid carcinoma. Sci Rep (2017) 7:44033. doi:10.1038/srep44033

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Geraldo MV, Kimura ET. Integrated analysis of thyroid cancer public datasets reveals role of post-transcriptional regulation on tumor progression by targeting of immune system mediators. PLoS One (2015) 10:e0141726. doi:10.1371/journal.pone.0141726

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Girardi C, De Pittà C, Casara S, Calura E, Romualdi C, Celotti L, et al. Integration analysis of microRNA and mRNA expression profiles in human peripheral blood lymphocytes cultured in modeled microgravity. Biomed Res Int (2014) 2014:296747. doi:10.1155/2014/296747

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Ab Mutalib NS, Othman SN, Yusof AM, Suhaimi SNA, Muhammad R, Jamal R. Integrated microRNA, gene expression and transcription factors signature in papillary thyroid cancer with lymph node metastasis. PeerJ (2016) 4:e2119. doi:10.7717/peerj.2119

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Yue K, Wang X, Wu Y, Zhou X, He Q, Duan Y. microRNA-7 regulates cell growth, migration and invasion via direct targeting of PAK1 in thyroid cancer. Mol Med Rep (2016) 14:2127–34. doi:10.3892/mmr.2016.5477

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Li D, Jian W, Wei C, Song H, Gu Y, Luo Y, et al. Down-regulation of miR-181b promotes apoptosis by targeting CYLD in thyroid papillary cancer. Int J f Clin Exp Pathol (2014) 7:7672–80.

PubMed Abstract | Google Scholar

53. Reinhart BJ, Slack FJ, Basson M, Pasquinelli AE, Bettinger JC, Rougvie AE, et al. The 21-nucleotide let-7 RNA regulates developmental timing in Caenorhabditis elegans. Nature (2000) 403:901–6. doi:10.1038/35002607

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Geraldo M, Yamashita A, Kimura E. MicroRNA miR-146b-5p regulates signal transduction of TGF-β by repressing SMAD4 in thyroid cancer. Oncogene (2012) 31:1910–22. doi:10.1038/onc.2011.381

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Liu D, Deng Q, Sun L, Wang T, Yang Z, Chen H, et al. A Her2-let-7-β2-AR circuit affects prognosis in patients with Her2-positive breast cancer. BMC Cancer (2015) 15:832. doi:10.1186/s12885-015-1869-6

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Wang Z, Zhang H, Zhang P, Li J, Shan Z, Teng W. Upregulation of miR-2861 and miR-451 expression in papillary thyroid carcinoma with lymph node metastasis. Med Oncol (2013) 30:577. doi:10.1007/s12032-013-0577-9

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Xing M. Molecular pathogenesis and mechanisms of thyroid cancer. Nat Rev Cancer (2013) 13:184–99. doi:10.1038/nrc3431

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Xing M, Alzahrani AS, Carson KA, Shong YK, Kim TY, Viola D, et al. Association between BRAF V600E mutation and recurrence of papillary thyroid cancer. J Clin Oncol (2015) 33:42–50. doi:10.1200/JCO.2014.56.8253

CrossRef Full Text | Google Scholar

59. Yarchoan M, LiVolsi VA, Brose MS. BRAF mutation and thyroid cancer recurrence. J Clin Oncol (2015) 33:7–8. doi:10.1200/JCO.2014.59.3657

CrossRef Full Text | Google Scholar

60. Han PA, Kim HS, Cho S, Fazeli R, Najafian A, Khawaja H, et al. Association of BRAF V600E mutation and microRNA expression with central lymph node metastases in papillary thyroid cancer: a prospective study from four endocrine surgery centers. Thyroid (2016) 26:532–42. doi:10.1089/thy.2015.0378

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Cahill S, Smyth P, Denning K, Flavin R, Li J, Potratz A, et al. Effect of BRAF V600E mutation on transcription and post-transcriptional regulation in a papillary thyroid carcinoma model. Mol Cancer (2007) 6:21. doi:10.1186/1476-4598-6-21

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Lupi C, Giannini R, Ugolini C, Proietti A, Berti P, Minuto M, et al. Association of BRAF V600E mutation with poor clinicopathological outcomes in 500 consecutive cases of papillary thyroid carcinoma. J Clin Endocrinol Metab (2007) 92:4085–90. doi:10.1210/jc.2007-1179

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Ito Y, Yoshida H, Maruo R, Morita S, Takano T, Hirokawa M, et al. BRAF mutation in papillary thyroid carcinoma in a Japanese population: its lack of correlation with high-risk clinicopathological features and disease-free survival of patients. Endocr J (2009) 56:89–97. doi:10.1507/endocrj.K08E-208

PubMed Abstract | CrossRef Full Text | Google Scholar