MicroRNA Profile for Diagnostic and Prognostic Biomarkers in Thyroid Cancer

Simple Summary Thyroid nodules are frequently detected, but the majority of these nodules are benign. Molecular markers have become useful for diagnosing the minority of thyroid cancer. We performed high-throughput small RNA sequencing in a discovery cohort and identified three miRNAs (miR-136, miR-21, and miR-127) as potential biomarkers for thyroid cancer. We validated the diagnostic and prognostic utilities of three miRNAs in patients with thyroid cancer in an independent cohort. High expression of the three miRNAs could be used to differentiate thyroid cancers from benign tumors and tumors with extremely low malignant potential. In patients with thyroid cancers, a high expression of three miRNAs was associated with poor clinicopathological features and recurrent or persistent disease following surgery. Therefore, testing for high expression levels of these three miRNAs in thyroid nodules may be useful for diagnosing and assessing the recurrence of the thyroid cancer’s risk stratification. Abstract The challenge in managing thyroid nodules is to accurately diagnose the minority of those with malignancy. We aimed to identify diagnostic and prognostic miRNA markers for thyroid nodules. In a discovery cohort, we identified 20 candidate miRNAs to differentiate between noninvasive follicular thyroid neoplasms with papillary-like nuclear features (NIFTP) and papillary thyroid carcinomas (PTC) by using the high-throughput small RNA sequencing method. We then selected three miRNAs (miR-136, miR-21, and miR-127) that were differentially expressed between the PTC follicular variant and other variants in The Cancer Genome Atlas data. High expression of three miRNAs differentiated thyroid cancer from nonmalignant tumors, with an area under curve (AUC) of 0.76–0.81 in an independent cohort. In patients with differentiated thyroid cancer, the high-level expression of the three miRNAs was an independent indicator for both distant metastases and recurrent or persistent disease. In patients with PTC, a high expression of miRNAs was associated with an aggressive histologic variant, extrathyroidal extension, distant metastasis, or recurrent or persistent disease. Three miRNAs may be used as diagnostic markers for differentiating thyroid cancers from benign tumors and tumors with extremely low malignant potential (NIFTP), as well as prognostic markers for predicting the risk of recurrent/persistent disease for differentiated thyroid cancer.


Introduction
The prevalence of thyroid nodules varies with the detection methods. Thyroid nodules are detected in 5% of the adults on palpation [1] and 33% to 68% on ultrasound [2,3]. Although most thyroid nodules are benign, thyroid cancer is the most common endocrine malignancy. The age-standardized (world population) incidence rate of thyroid cancer available online at the Global Cancer Observatory (https://gco.iarc.fr/) is 10.2 per 100,000 females and 3.1 per 100,000 males in 2018 [4]. Papillary thyroid carcinoma (PTC) is the most common subtype, accounting for more than 80% of all thyroid cancers [5]. The incidence of thyroid cancer has been increasing over the past three decades in many countries, mainly because of an increase in the diagnosis of PTC and improvements in the detection and diagnosis of small cancers [6][7][8][9][10]. However, the mortality rate from thyroid cancers has remained stable over the same period despite the increasing incidence of thyroid cancer [4,10,11]. The age-standardized mortality rate was 0.42 per 100,000 people in 2018 [4]. The estimated overdiagnosis rates for thyroid cancers range from 50% to 90% of newly diagnosed cases in women, depending upon the regions and health care environment [4]. Therefore, to reduce overtreatment of nonfatal thyroid cancers, there is a clinical need for biomarkers to enable a more accurate diagnosis and risk assessment for thyroid cancers.
Several studies have evaluated the diagnostic and prognostic role of miRNA expression levels in patients with thyroid nodules [13,14,17]. Traditionally, the expression profiling of miRNAs has been studied by hybridization-based methods including microarrays. Recently, next-generation sequencing (NGS)-based RNA-seq has emerged as an alternative gene expression profiling technology for miRNAs. The NGS technology can cover all coding and noncoding RNAs and more accurately measure their expression level changes than can microarray platforms. Despite advances in molecular thyroid cancer diagnostics, there is a limited number of studies of NGS-based miRNA diagnostics for thyroid cancer in the literature.
In the present study, we aimed to obtain high-throughput miRNA-expression profiles and subsequently identify miRNA markers to differentiate thyroid cancers from noninvasive follicular thyroid neoplasm with papillary-like nuclear features (NIFTP) and benign thyroid nodules. Furthermore, the role of miRNA markers in the risk stratification to predict recurrence was further revealed by analyzing the association of their expression levels with clinicopathologic parameters in an independent cohort.
By applying unsupervised hierarchical clustering with a total of 712 miRNAs showing varying expression changes across thyroid samples (standard deviation (SD) > 0.7), we found that thyroid tumor samples were well differentiated from normal thyroid tissues. When the histological data was compared among the tumor samples, those with NIFTP or IEFVPTC were found to be differently clustered from those with classic PTC or TCVPTC ( Figure 1A).
By exploring the expression patterns, we identified a distinct subset of 62 miRNAs that showed significantly different expression patterns between patients with NIFTP/IEFV and classic PTC/TCVPTC ( Figure 1A). These miRNA expression levels were higher in patients with classic PTC/TCVPTC than in those with other histologic types. Among these 62 miRNAs, many of them that were involved in focal adhesion, cell proliferation, epidermal growth factor receptor (EGFR) signaling pathways, or histone deacetylase binding were significantly enriched ( Figure 1A).
To identify subsets of miRNAs exclusive to NIFTP, we analyzed differentially expressed miRNAs between the classic PTC and NIFTP, TCVPTC and NIFTP, or IEFVPTC and NIFTP patient subgroups by using Venn diagrams ( Figure 1B). Seven miRNAs were significantly upregulated in the NIFTP, whereas 13 miRNAs were significantly downregulated in the NIFTP ( Figure 1C; Table 2).

Validation of Potential miRNA Markers in TCGA Data
Using The Cancer Genome Atlas (TCGA) data set, we sought to verify these twenty miRNAs that were differentially expressed in the NIFTP of the original cohort. Among the candidates, three miRNAs (i.e., miR-136, miR-21, and miR-127) showed significant differences in expression levels between the follicular variant and the others (i.e., classic PTC, TCVPTC, or other variants) as shown in Figure 2, consistent with the results in the discovery cohort. These three miRNAs were applied into the next further downstream experimental validations.

Validation of Potential miRNA Markers in TCGA Data
Using The Cancer Genome Atlas (TCGA) data set, we sought to verify these twenty miRNAs that were differentially expressed in the NIFTP of the original cohort. Among the candidates, three miRNAs (i.e., miR-136, miR-21, and miR-127) showed significant differences in expression levels between the follicular variant and the others (i.e., classic PTC, TCVPTC, or other variants) as shown in Figure 2, consistent with the results in the discovery cohort. These three miRNAs were applied into the next further downstream experimental validations. P-values were obtained by 2-sample t-tests. CPM: counts per million mapped reads. *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Validation of Potential miRNA Markers in an Independent Cohort
The integrity of total RNA isolated from FFPE samples showed an RNA integrity number (RIN) ranging 2.0 to 3.0. Because expression of miRNA is highly stable, even in degraded human tissue RNA and cell samples [18][19][20], we used the RNA samples extracted from 233 FFPE tissue samples for the miRNA marker validation study by using a TaqMan-based qRT-PCR assay.
We observed that expression of three candidate miRNAs was more highly expressed in PTC than in other thyroid tumors ( Figure 3A-C). The three selected miRNA markers were further evaluated for their capacity to differentiate nonmalignant tumors, including follicular adenoma (FA) and NIFTP from malignant tumors, including PTC, follicular thyroid carcinoma (FTC), and Hürthle cell carcinoma ( Figure 3D-F). P-values were obtained by 2-sample t-tests. CPM: counts per million mapped reads. *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Validation of Potential miRNA Markers in an Independent Cohort
The integrity of total RNA isolated from FFPE samples showed an RNA integrity number (RIN) ranging 2.0 to 3.0. Because expression of miRNA is highly stable, even in degraded human tissue RNA and cell samples [18][19][20], we used the RNA samples extracted from 233 FFPE tissue samples for the miRNA marker validation study by using a TaqMan-based qRT-PCR assay.
We observed that expression of three candidate miRNAs was more highly expressed in PTC than in other thyroid tumors ( Figure 3A-C). The three selected miRNA markers were further evaluated for their capacity to differentiate nonmalignant tumors, including follicular adenoma (FA) and NIFTP from malignant tumors, including PTC, follicular thyroid carcinoma (FTC), and Hürthle cell carcinoma ( Figure 3D-F).
To estimate the area under the curve (AUC) value, we performed receiver operating characteristic (ROC) analysis. The AUC values of miR-136, miR-21, and miR-127 were 0.76, 0.83, and 0.83, respectively ( Figure 4A-C). Their optimal cutoff values and sensitivity and specificity for each marker of the diagnosis of malignancy are described in Figure 3.

Clinicopathologic Utility of the Three miRNA Markers in Thyroid Cancer
The expression levels of three miRNAs were classified into two sub-groups of high-or low-expression based on their median expression values. We then categorized all patients into four groups based on the number of markers showing high miRNA expression levels: all low (group 1, n = 88), one high (group 2, n = 24), two high (group 3, n = 28), or all high (group 4, n = 93).
Benign tumors (FA) and tumors with extremely low malignant potential (NIFTP) were mostly found in group 1 whereas malignant tumors (PTC, Hürthle cell carcinoma, and FTC) were mostly found in group 4 ( Figure 5A). Patients classified as high-risk according to the American Thyroid Association (ATA) recurrence risk stratification system were most frequently observed in group 4 ( Figure 5B). Recurrent or persistent diseases were mostly found in group 4 ( Figure 5C). Stage III and IV were only found in group 4 ( Figure 5D). Hürthle cell carcinoma. The area under the curve (AUC) indicates the probability that the classifier ranks a randomly chosen positive instance higher than a randomly chosen negative instance. The criteria for high miRNA expression levels of miR-136 (A), miR-21 (B), and miR-127 (C) were defined using the AUC.

Clinicopathologic Utility of the Three miRNA Markers in Thyroid Cancer
The expression levels of three miRNAs were classified into two sub-groups of highor low-expression based on their median expression values. We then categorized all patients into four groups based on the number of markers showing high miRNA expression levels: all low (group 1, n = 88), one high (group 2, n = 24), two high (group 3, n = 28), or all high (group 4, n = 93).
Benign tumors (FA) and tumors with extremely low malignant potential (NIFTP) were mostly found in group 1 whereas malignant tumors (PTC, Hürthle cell carcinoma, and FTC) were mostly found in group 4 ( Figure 5A). Patients classified as high-risk according to the American Thyroid Association (ATA) recurrence risk stratification system were most frequently observed in group 4 ( Figure 5B). Recurrent or persistent diseases were mostly found in group 4 ( Figure 5C). Stage III and IV were only found in group 4 ( Figure 5D). Multivariate generalized linear model analysis showed that high expression levels of three miRNA markers were significant predictors for recurrent or persistent disease (odds Multivariate generalized linear model analysis showed that high expression levels of three miRNA markers were significant predictors for recurrent or persistent disease (odds ratio, 3.56; 95% confidence interval [CI], 1.57-8.07) and distant metastases (odds ratio, 4.52; 95% CI, 1.71-11.93) in 133 patients with differentiated thyroid cancer ( Figure 6). ratio, 3.56; 95% confidence interval [CI], 1.57-8.07) and distant metastases (odds ratio, 4.52; 95% CI, 1.71-11.93) in 133 patients with differentiated thyroid cancer ( Figure 6).

Clinicopathologic Utility of the Three miRNA Markers in Patients with PTC
A subgroup analysis was conducted in 119 patients with PTC, as shown in Table 3. High expression levels of miRNAs (miR-136, miR-21, and miR-127) were associated with histologic variants (P < 0.001), extrathyroidal extension (P = 0.003), distant metastases (P = 0.038), recurrent or persistent disease (P = 0.013), and increased ATA recurrence risk (P < 0.001). In addition, expression levels of miR-21 were associated with lymph node metastasis and BRAF V600E mutation.

Clinicopathologic Utility of the Three miRNA Markers in Patients with PTC
A subgroup analysis was conducted in 119 patients with PTC, as shown in Table 3. High expression levels of miRNAs (miR-136, miR-21, and miR-127) were associated with histologic variants (P < 0.001), extrathyroidal extension (P = 0.003), distant metastases (P = 0.038), recurrent or persistent disease (P = 0.013), and increased ATA recurrence risk (P < 0.001). In addition, expression levels of miR-21 were associated with lymph node metastasis and BRAF V600E mutation.

Discussion
Many miRNAs have been investigated as biomarkers for thyroid cancer, and their number is continuously increasing. The main role of these miRNAs is the regulation of cancer-related signaling pathways involved in cell proliferation, survival, adhesion, invasion, and migration [14]. Generally, oncomiRs and tumor suppressor miRNAs are overexpressed and underexpressed in thyroid cancer, respectively. Group profiling strategies have been used to identify the dysregulation of miRNAs in thyroid cancer. However, there has been inconsistent miRNA profiling data in different studies, depending on the measurement methods and the algorithms necessary to associate the miRNAs. We performed small RNA-seq analysis to detect both novel and known miRNAs to differentiate benign tumors (FA) and tumors with extremely low malignant potential (NIFTP) from malignant thyroid tumors. Three candidate miRNAs were validated for diagnostic and prognostic utility in an independent cohort by TaqMan-based qRT-PCR assay.
NIFTP is a borderline thyroid tumor that has an indolent clinical outcome and overlapping histologic and molecular features with FA and with IEFVPTC [5]. It is difficult to distinguish these follicular-patterned tumors by clinical and cytologic examination before surgery. Molecular analysis cannot accurately differentiate these tumors because NIFTPs have molecular profiles similar to those of FA and IEFVPTC, with frequent RAS-like mutations [21][22][23]. However, recent studies have shown the miRNA expression profiles of NIFTPs to be different from those of other thyroid tumors [24,25].
We showed that the expression of three miRNA markers was higher in differentiated thyroid cancers than nonmalignant tumors (FA and NIFTP). A high-expression level of these miRNAs was an independent prognostic indicator for both distant metastases and recurrent or persistent disease in patients with differentiated thyroid cancers. Most cases of FA and NIFTP had low expression levels of three miRNA markers. In patients with PTC, a high-expression of miRNAs was associated with aggressive histologic variants, extrathyroidal extensions, distant metastases, recurrent or persistent disease, and a high recurrence risk as defined by the ATA.
Although this is an extremely limited study, the expression of miR-136-5p was higher in thyroid neoplasms, especially in invasive PTC with lymph node metastasis than it was in nontumor tissues [31]. In the TCGA data set, a high expression of miR-136-5p was associated with older age (>50 years), high pathologic T category, lymph node metastasis, and advanced cancer stage [32]. However, contradictory studies exist concerning whether miR-136 is an oncogene or a tumor suppressor. The expression of miR-136 acts as an oncogene in lung cancer [33] and gastric cancer [34]. A previous study suggested that the expression of miR-136 downregulates the PTEN expression and subsequent activation of AKT expression in gastric cancer cells [34]. However, other studies reported miR-136 plays a tumor suppressor role in glioma [35], triple-negative breast cancer [36], colon cancer [37], and gallbladder cancer [38]. A previous study reported that miR-136 suppressed the epithelial-to-mesenchymal transition of cancer cells by regulating the RAS protein activator like 2 (RASAL2) gene in triple-negative breast cancer [36]. Another study done on colon cancer reported that miR-136 suppressed the epithelial-to-mesenchymal transition of cancer cells by targeting the migration and invasion enhancer 1 (MIEN1) gene. Its expression level in colon cancer tissue was inversely correlated with tumor size, lymph node metastasis, and cancer stage [37]. In the gallbladder, overexpression of miR-136 suppressed cancer cell proliferation and enhanced apoptosis of cancer cells via inhibition of mitogen-activated protein kinase 4 (MAP2K4) gene expression [38].
As with miR-136, the miR-127 plays a dual role as an oncogene and a tumor suppressor, depending upon the cancer tissue type involved. In the TCGA data set, the miR-127 expression level was higher in stage III PTC than in stage II tumors [39]. In medullary thyroid cancer, miR-127 was more frequently upregulated in RET wild-type cancers [40]. In our study, the high expression of miR-127 was associated with thyroid cancers and a predictive marker for tumor recurrence. In contrast, miR-127 has also been reported to act as a tumor suppressor in various tumor types. The target genes downregulated by miR-127 have been shown as follows: B-cell lymphoma 6 (BCL6) in breast cancer [41], replication initiator 1 (REPIN1) in glioma [42], formin-like 3 (FMNL3) in esophageal cancer [43], cytochrome C oxidase assembly factor 1 homologue (COA1) in giant cell tumor of bone [44], and delta like non-canonical notch ligand 1 (DLK1) in melanoma [45].
Therefore, miR-136 and miR-127 play different roles depending on the cancer types and microenvironment. Further studies are required to reveal the oncogenic role of miR-136 and miR-127 in thyroid cancers.
Intracellular and extracellular miRNAs are relatively stable in clinical tissues, blood, and body fluids under various conditions [46][47][48]. These characteristics enable reliable detection of common miRNAs in different types of cytology and pathology specimens. Currently, two miRNA-based molecular tests are commercially available for indeterminate thyroid nodules on fine needle aspiration cytology [49]: ThyraMIR (Interpace Diagnostics, Parsippany, NJ, USA) using a panel of 10 miRNAs and RosettaGX Reveal (Rosetta Genomics, Philadelphia, PA, USA) using a panel of 24 miRNAs. The miRNA-base tests classify cytologically indeterminate thyroid nodules into benign or suspicious for malignancy by miRNA profiling. This practice highlights the clinical relevance of miRNA tests in the diagnosis of thyroid nodules. As the miRNA panels in both tests do not include our three miRNA markers, the diagnostic performance of miRNA testing may increase by incorporating our three miRNAs. However, further studies are needed to validate the clinical utility of miRNA markers on cytology samples.  Table 1. The present study was performed on the same cohort and tissue samples used in our previous study investigating DNA methylation biomarkers [50].
In addition, we used 233 FFPE tissues samples to validate our selected miRNA markers by qRT-PCR. The validation cohort consisted of FA (n = 43), NIFTP (n = 57), invasive EFVPTC (n = 22), classic PTC (n = 49), TCV PTC (n = 45), other variants of PTC (n = 3), follicular thyroid carcinoma (FTC, n = 12), and Hürthle cell carcinoma (n = 2). Hürthle cell adenoma was not included within the FA group. Thyroid tumors were pathologically classified according to the 2017 World Health Organization classification of tumors of endocrine organs [5]. Cancer stages were categorized according to the 8th edition of the American Joint Committee on Cancer (AJCC) staging manual [51]. Recurrence risk was evaluated using the ATA classification for risk of recurrence [52].

Total RNA Isolation and Small RNA Sequencing
Total RNA was isolated from the fresh-frozen tissues and 10-µm thick paraffinembedded tissue sections using the RecoverAll™ Total Nucleic Acid Isolation Kit (Life Technologies, Carlsbad, CA, USA), according to the manufacturer's instructions. The quality and quantity of the extracted total RNA were analyzed with an ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). In addition, the 2100 Agilent Bioanalyzer (Agilent Technologies, Waldbronn, Germany) was used for the estimation of the RNA integrity number (RIN) score.
Small RNAs (20-30 nt) were purified from 15% Novex TBE-Urea Gel (Invitrogen) and ligated first with the 5' RNA adaptor and then with the 3' RNA adaptor provided by Illumina TruSeq small RNA sample preparation protocol. In each step, the ligated product was PAGE-gel purified. After first-strand synthesis and 11 cycles of PCR amplification, the product was PAGE-gel purified and submitted for sequencing on an Illumina NextSeq500 at LAS (Seoul, Korea; https://www.las.kr/).

Small RNA Sequencing Data Analysis
After TruSeq small RNA adapters were eliminated among the sequenced reads by Trimmomatic software (v. 0.38), the remained sequence data was mapped to the human genome (GRCh38) using bowtei2 (v. 2.3.4). The quantification was carried out using the HTSeq software package. The dataset generated by small RNA-seq is available in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) public database under data series accession number GSE159330

Public miRNA-Sequencing Data Collection
Public miRNA-sequencing data of normal thyroid and PTC samples was obtained from the TCGA dataset (https://portal.gdc.cancer.gov/) to confirm the miRNA patterns of candidate miRNA markers.

MiRNA Expression Levels by qRT-PCR
To evaluate the miRNA expression level, we used the TaqMan MicroRNA assay (Thermo Fisher Scientific) according to manufacturer's instructions. Briefly, 10 ng of total RNA was reverse transcribed into cDNA with each miRNA-specific primer, and then 1.44 uL of cDNA was used for the CFX96 real-time PCR system (Bio-Rad, Hercules, CA, USA). The PCR cycles were as follows: initial denaturation at 95 • C for 10 min, followed by 40 cycles at 94 • C for 15 s, and 60 • C for one minute. The expression level of each miRNA was normalized to a noncoding small nuclear U6 RNA expression that was commonly used as a reference in previous studies [53][54][55]. The relative miRNA expression level was calculated using the following formula: 2ˆ(U6 ct − target miRNA ct ) + 15.

BRAF Mutational Analysis
Genomic DNA was isolated from 10-µm thick paraffin-embedded tissue sections using the RecoverAll™ Total Nucleic Acid Isolation Kit (Life Technologies). After PCR amplification of the extracted DNA, the sequences of BRAF exon 15 were analyzed by direct sequencing of amplicons, as described previously [50,56,57].

Statistical Analysis
To perform miRNA expression profiling for thyroid tissue samples, a hierarchical clustering algorithm with the centered correlation coefficient as the measure of similarity and complete linkage clustering was applied. For cluster analysis, the read counts per million fragments mapped (CPM) of each sample were used to estimate the expression level of each miRNA. The CPM data was normalized by the quantile method, log 2 transformed, and median centered across genes and samples. To estimate the significance of differences in gene expression between sample subgroups, the edgeR package, which uses a negative binomial model, was employed to detect differentially expressed miRNAs from the count data [58]. Expression differences in miRNAs were considered to be statistically significant if the P-value was < 0.05 and the fold difference in expression between two sample groups was ≥ 1.5. Function enrichment analysis was carried out to identify the most significant miRNA sets with miEAA software [6]. The significance of overrepresented gene sets was estimated by Fisher's exact test. Statistical analysis was mainly carried out in the R language environment (ver. 3.5.3).
The ROC and the respective area under the ROC curve (AUC) were calculated for each miRNA marker, using the ROCR package of the R software (version 3.4.0). ROC curve analysis estimated the optimal cutoff values maximizing sensitivity and specificity between low and high levels of miRNA expression. The relationship between clinicopathologic features and the expression levels of miRNAs was analyzed using parametric (chi-squared test) and non-parametric (Fisher's exact) assessments, where appropriate. Generalized liner model was performed to assess the associations of clinicopathologic variables and miRNA expression levels with the adverse clinical outcomes.

Conclusions
In conclusion, we report that upregulation of three miRNA markers, miR-136, miR-21 and miR-127, may serve as diagnostic and prognostic biomarker for differentiated thyroid cancer. These miRNA markers may be clinically useful for efficiently stratifying thyroid tumors.