MicroRNA-Based Risk Score for Predicting Tumor Progression Following Radioactive Iodine Ablation in Well-Differentiated Thyroid Cancer Patients: A Propensity-Score Matched Analysis

Simple Summary The three-tiered American Thyroid Association (ATA) risk stratification helps clinicians tailor decisions regarding follow-up modalities and the need for postoperative radioactive iodine (RAI) ablation and radiotherapy. However, a significant number of well-differentiated thyroid cancers (DTC) progress after treatment. Current follow-up modalities have also been proposed to detect disease relapse and recurrence but have failed to be sufficiently sensitive or specific to detect, monitor, or determine progression. Therefore, we assessed the predictive accuracy of the microRNA-based risk score in DTC with and without postoperative RAI. We confirm the prognostic role of triad biomarkers (miR-2f04, miR-221, and miR-222) with higher sensitivity and specificity for predicting disease progression than the ATA risk score. Compared to indolent tumors, a higher risk score was found in progressive samples and was associated with shorter survival. Consequently, our prognostic microRNA signature and nomogram provide a clinically practical and reliable ancillary measure to determine the prognosis of DTC patients. Abstract To identify molecular markers that can accurately predict aggressive tumor behavior at the time of surgery, a propensity-matching score analysis of archived specimens yielded two similar datasets of DTC patients (with and without RAI). Bioinformatically selected microRNAs were quantified by qRT-PCR. The risk score was generated using Cox regression and assessed using ROC, C-statistic, and Brier-score. A predictive Bayesian nomogram was established. External validation was performed, and causal network analysis was generated. Within the eight-year follow-up period, progression was reported in 51.5% of cases; of these, 48.6% had the T1a/b stage. Analysis showed upregulation of miR-221-3p and miR-222-3p and downregulation of miR-204-5p in 68 paired cancer tissues (p < 0.001). These three miRNAs were not differentially expressed in RAI and non-RAI groups. The ATA risk score showed poor discriminative ability (AUC = 0.518, p = 0.80). In contrast, the microRNA-based risk score showed high accuracy in predicting tumor progression in the whole cohorts (median = 1.87 vs. 0.39, AUC = 0.944) and RAI group (2.23 vs. 0.37, AUC = 0.979) at the cutoff >0.86 (92.6% accuracy, 88.6% sensitivity, 97% specificity) in the whole cohorts (C-statistics = 0.943/Brier = 0.083) and RAI subgroup (C-statistic = 0.978/Brier = 0.049). The high-score group had a three-fold increased progression risk (hazard ratio = 2.71, 95%CI = 1.86–3.96, p < 0.001) and shorter survival times (17.3 vs. 70.79 months, p < 0.001). Our prognostic microRNA signature and nomogram showed excellent predictive accuracy for progression-free survival in DTC.


Introduction
The incidence of thyroid cancer (TC) in the United States (US) was estimated to be over 52,890 cases (12,720 men and 40,170 women) in 2020 and is projected to be the fourth most common cancer by 2030 [1]. The annual incidence has increased by 211% over the last two decades, mostly due to increased detection of papillary thyroid microcarcinoma (PTMC) [2][3][4]. Papillary thyroid cancer (PTC) and follicular thyroid cancer (FTC) account for most thyroid cancers [5]. These well-differentiated thyroid cancers (DTC) are generally treated with surgical resections followed by adjuvant radioactive iodine (RAI) therapy to ablate the remnant or residual thyroid tissue [6,7]. Following surgery, patients suffer from a significant risk of complications, need for hormone replacement therapy, and lengthy postoperative surveillance with unnecessarily higher health care costs and diminished quality of life.
Based on the revised 2015 American Thyroid Association (ATA) guidelines, DTC patients were categorized into low, intermediate, and high-risk groups according to the estimated risk of recurrence and cancer relapse [8]. This three-tiered risk stratification system helps to tailor decisions regarding the need for postoperative thyrotropin suppression, radioactive iodine ablation, or radiotherapy, as well as the frequency and modality of follow-up studies required [9]. Despite the good prognosis of DTC, up to 30% of patients experience recurrences in the thyroid bed or neck lymph nodes after initial treatment [10]. The modalities currently used for TC assessment, including ultrasound and fine-needle aspiration biopsy, offer only a snapshot of the disease in a single point of time and do not describe tumor behavior over time. The sonographic appearance and cytopathology findings may sometimes be non-informative and inconclusive [11]. Several clinical parameters have also been proposed but have failed to be sufficiently sensitive or specific to detect, monitor, or determine progression. Serum thyroglobulin (sTg) has been reported as a predictor for treatment efficacy during ablative radioiodine treatment [12]; however, some patients still exhibit elevated sTg levels even after receiving adjuvant RAI therapy [13,14]. As reported by Yim et al. [15], 11% of PTC patients who underwent bilateral thyroidectomy followed by RAI remnant ablation developed recurrence, with only 36.1% of these showing high sTg levels.
Similarly, in a study by Hirsch et al., 47% had a persistent disease despite re-treatment with RAI [16]. To date, there are no molecular markers that can predict tumor recurrence or persistence. Hence, there is a critical need to discover new biomarkers to biologically define which cancers have an aggressive form and optimize the selection criteria of management plans. In the absence of such a prognostic panel, prediction for progression will continue to be a challenging practice that will impede well-being and make patients and clinicians reluctant to choose active surveillance instead of surgical intervention.
As central regulators of gene expression, microRNAs (miRNAs) are attracting increasing attention because of their association with tumor development and progression. MiRNAs are short non-coding RNAs of around 18-23 nucleotides that regulate virtually all biological functions via post-transcriptional gene silencing [17]. MiRNAs can act as either oncogenes or tumor suppressor genes in thyroid cancer [18]. Altered expression levels of miRNAs influence apoptosis, migration and proliferation, angiogenesis, and local immune response. Distinct miRNA expression profiles are also associated with well-defined clinicopathological features of thyroid cancer [19] and prognosis/disease progression [20], as depicted in in vitro and in vivo studies.
There is increasing interest in the association of miRNA expression with chemo-and radiosensitivity for predicting or modulating resistance [21]. For example, upregulation of the miRNA-221/-222 and miRNA-17-92 cluster significantly increases the radioresistance of cancer cells through the downregulation of phosphatase and tension homolog (PTEN) and pAkt activity [22,23], while miR-145 treatment effectively increases the sensitivity of cells to radiation [24]. Furthermore, emerging evidence also demonstrates miRNAs as promising therapeutic targets in thyroid cancer. Upregulation of the miR-17-92 cluster could provide a promising therapeutic modality to counteract ATC progression [25,26]. Similarly, miR-204-5p upregulation plays a protective role by inhibiting PTC cell proliferation through regulating IGFBP5 expression [27,28]. In contrast, inhibiting oncomiRs inducing metastasis such as miR-146a and miR-146b in PTC cells through targeting IRAK1 or Wnt/β-catenin pathways might provide ancillary therapeutic strategies in conjunction with the current status quo regimens [29][30][31][32]. Therefore, reversing the altered miRNA signature may pave the road toward a cure.
No specific miRNA panel has overcome the hurdle of predicting recurrence and response to therapy, especially in the absence of lymph node metastasis and extrathyroidal extension. Therefore, we aimed to identify and validate a microRNome signature to predict recurrence at the time of surgery in well-differentiated thyroid cancer patients following radioactive iodine ablation. Analysis of TC datasets in The Cancer Genome Atlas (TCGA) database revealed several deregulated miRNAs in patients who developed recurrent/persistent disease postoperatively. Of these, miR-221, miR-222, and miR-204 consistently predicted recurrence before and after radioactive remnant ablation treatment. We validated our findings in genome-wide miRNA expression profiling studies and patient samples. Our results demonstrate the putative role of the triad biomarker as an effective prognostic signature that accurately predicts recurrence following RAI treatment in welldifferentiated thyroid cancer patients.

Bioinformatic Selection of MiRNAs
Transcriptomic signatures of 495 thyroid cancer patients were retrieved from the Genomic Data Commons (GDC) data portal for the Cancer Genome Atlas thyroid cancer dataset (TCGA-THCA) (https://www.cancer.gov/about-nci/organization/ccg/research/ structural-genomics/tcga) (accessed on 15 March 2021), and 1035 miRNAs from miRNAseq were included. Clinical, pathological, and molecular information was obtained from cBioPortal for Cancer Genomics (https://www.cbioportal.org) (accessed on 15 March 2021) and FireBrowse (http://www.firebrowse.org/) (accessed on 15 March 2021). Outcomes of interest were disease recurrence and/or progression. Patients with incomplete recurrence data or unmatched miRNA samples were excluded. Ultimately, 448 non-recurrent and 47 recurrent cancer patients were included. We classified TCGA-THCA cohorts according to the updated 2015 ATA risk stratification for structural disease recurrence into low (≤5%), intermediate (5-20%), or high-risk (≥20%) groups and reported the percentage of the expected risk of recurrence as a quantitative score [8]. A systematic search was performed in the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/) (accessed on 15 March 2021), an online public functional genomics data repository for high-throughput datasets.

Study Population and Propensity Score-Matched Groups
We retrospectively reviewed 788 patients recruited between January 2010 and December 2015 from Elbayan Pathology Laboratory, Port-Said, Egypt. Samples were collected during thyroidectomy due to thyroid cancer diagnosis. The study population comprised adult cohorts (aged 18 years or older) diagnosed with well-differentiated thyroid cancer (PTC or FTC) according to the International Classification of Oncological Diseases, 4th edition. Patients did not receive any treatment before operative resection. Exclusion criteria included Hürthle cell thyroid carcinoma, poorly differentiated thyroid carcinoma, anaplastic (undifferentiated) carcinoma, medullary thyroid cancer, thyroid lymphoma, thyroid cancer arising from a thyroglossal duct cyst, and thyroid cancer in malignant struma ovarii. Patients with incomplete follow-up or missing data were also excluded. As depicted in the selection process of samples in Figure 1, a total of 222 paired cancer and non-cancer adjacent tissues were eligible. Since the presence of confounders may favor the use of RAI ablation and lead to biased analysis, a confounder-elimination process was conducted using a 1:1 propensity score analysis. Matching yielded two similar datasets of 68 paired samples (cancer and non-cancer tissues) of cohorts who underwent surgical resection of tumor vs. those who had RAI treatment after thyroidectomy. Workflow for patient selection. ATA: American Thyroid Association. We reviewed 788 thyroid cancer samples for patients who underwent subtotal/total thyroidectomy for papillary or follicular thyroid carcinoma. Histopathological analysis was performed, and samples with insufficient tissue available for molecular work or lack of paired non-cancer tissues were excluded. After accounting for selection bias through propensity score analysis, two groups were established with similar general baseline features: (1) 34 paired cancer and non-cancer tissues for patients who underwent thyroidectomy and (2) another 34 paired tissues for those who received RAI following tumor resection.

Study Variables and Clinical Assessment
From patient records, we obtained demographic and clinicopathological characteristics, treatment strategies, response to therapy, recurrence, and mortality. Demographic variables included age at diagnosis, sex, and year of diagnosis. Disease characteristics included histologic subtype (PTC or FTC), TNM stage (8th edition), presence or absence of multifocal disease, and minor or gross extrathyroidal extension. Patients were classified according to surgery into none, lobectomy, subtotal/near-total thyroidectomy, or total thyroidectomy. Data regarding treatment consisted of the extent of surgery, use of radioactive iodine, and use of other treatment modalities (e.g., external beam radiotherapy). Adjuvant RAI therapy was defined as empirical RAI treatment performed after the first reoperation in patients with locoregionally recurrent PTC who initially underwent total thyroidectomy and remnant ablation. Clinical recurrence was defined as the reappearance of pathologically proven malignant tissue and/or the appearance of metastatic lesions such as lung, bone, and/or brain metastases. No clinical evidence of disease (NCED) was defined as the absence of disease, based on physical examination, neck ultrasonography or neck computed tomography (CT) scan, and any other imaging performed during clinical evaluation at the end of follow-up, regardless of serum thyroglobulin concentration.
Disease-free survival (or relapse-free survival) measures the number of people expected to be free from cancer for a particular amount of time to the time of death. Progressionfree survival (PFS) was defined as the period from initial treatment to new neoplasm, imaging evidence of disease or disease recurrence, or death. Overall survival refers to the time beginning at the start of treatment and up to the time of death and includes those who survive without any evidence of cancer and those who survive but still have cancer present in their bodies. For individuals who have no tumor relapse, we use the last follow-up time without a recurrence event.

Tissue Sample Preparation and Histopathological Assessment
Achieved formalin-fixed paraffin-embedded (FFPE) tissues of 68 thyroid cancer and 68 paired non-tumor thyroid samples. Histopathological diagnosis of well-differentiated thyroid cancer was confirmed, and samples were assessed for subtype variant, grading, and staging by two independent pathologists. Laser microdissection (Leica PBI Laser Microdissection model 7) was employed to identify regions of cancer and control tissues in FFPE specimens. Tissues were cut into 4 µm serial sections and stored at 4 • C until use. A 4 µm thick section was used for hematoxylin and Eosin (H&E) staining, and 3-4 sections in Eppendorf tubes were utilized for downstream qRT-PCR experiments.

RNA Extraction and MicroRNA Quantification
Total RNA, including small RNAs, was purified from FFPE samples by xylene deparaffinization and Qiagen miRNeasy FFPE Isolation kit (Qiagen, Hilden, Germany; Catalog number 217504) following the manufacturer's protocol. Nucleic acid concentration and purity were determined using the Nanodrop ND-1000 spectrophotometer (NanoDrop Tech. Inc., Wilmington, DE, USA) using the wavelength-dependent extinction coefficient of 33 [35,36]. Samples were stored in aliquots at −80 • C until analysis. RNA (10 ng) was converted to complementary DNA (cDNA) using TaqMan MiRNA Reverse Transcription (RT) kit (P/N 4366596); Thermo Fisher, CA, USA), and the 5× of specific stem-loop primers or endogenous control primers for normalization were used separately. The three miRNAs TaqMan assays are depicted in Table S1. Reverse transcription (RT) was carried out in a T-Professional Basic, Biometra PCR system (Biometra, Goettingen, Germany). Each studied miRNA was specifically converted to complementary DNA using the "TaqMan MiRNA RT kit (P/N 4366596; Applied Biosystems, Foster City, CA, USA)" with 5× miRNA-specific stem-loop primers at the following amplification conditions: 16 • C for 30 min, 42 • C for 30 min, and 85 • C for 5 min, then held at 4 • C. Successful removal of DNA contaminants was confirmed using no-reverse transcription controls of representative samples. qRT-PCR reactions were conducted following the "Minimum Information for publication of quantitative real-time PCR experiments (MIQE)" guidelines [37]. For each specified miRNA quantification, the final volume reaction of 20 µL included 1.33 µL RT product for the specified miR, 2× TaqMan Universal PCR master mix with UNG (Applied Biosystems, P/N 4440043) and 1 µL 20× TaqMan small RNA assay or small nuclear RNA U6 (RNU6B) (assay ID 001093). Three other endogenous control assays (i.e., RNU48, let-7a, and miR-16) were tested for data normalization based on the recent recommendation related to qRT-PCR for miRNAs and endogenous control assessment in PTC archived specimens [38]. Given the consistency and stability of RNU6B across samples, it was applied for data normalization. The PCR was performed in StepOne Real-time PCR system (Applied Biosystems) and incubated as follows: 95 • C for 10 min, followed by 45 cycles of 92 • C for 15 s and 60 • C for 1 min. Reactions were run in triplicate, and standard deviation >2.0 was set as an outlier. Appropriate controls were included in each run [39].

Statistical Analysis
The estimated power of the present study is 96% for a total of 68 matched paired TC samples, medium effect size = 0.5, and alpha error probability = 0.05, using G*Power version 3.1.9.2. With threshold cycle values acquired from the ABI 7900 HT SDS version 2.0.1 software (Applied Biosystems; adjusted settings at automatic baseline and threshold at 0.15), the relative miRNA expression levels were determined by the LIVAK method (2 -ddCq ), where ddC q (delta delta quantitative cycle) = (C q microRNA of interest − C q endogenous control ) cancer tissue − (C q microRNA of interest − C q endogenous control ) non-cancer tissue [40]. Wilcoxon matched-pairs signed-rank test was used for comparison of miRNA expression between cancer vs. normal tissues. Data are reported as the median and interquartile range (IQR) and plotted in box plots. Co-expression was estimated through Spearman's correlation test and plotted in a correlation matrix.
To test the prognostic value of miRNAs, overall and subgroup analysis of all patients (35 progressed vs. 33 indolent courses) and patients who received RAI treatment (20 progressed vs. 14 indolent courses) were performed. Mann-Whitney U or Kruskal-Wallis tests were performed for the comparison between progressed and indolent groups. Receiver operating characteristic (ROC) analysis was performed to test the predictive accuracy of the miRNAs, and the Youden statistic was used to select the best cutoff in discriminating patients' progression following RAI using the pROC R package [41]. MiRNAs with area under the curve (AUC) greater than 0.75 and p < 0.05 were set to be significant. DeLong test was applied using MedCalc (www.medcalc.org/) (accessed on 10 May 2021) to compare the AUC of multiple markers [42]. Accuracy measures including sensitivity, specificity, positive and negative predictive value, and likelihood ratios were calculated, and pooled weighted estimates of the significant miRNAs were determined using Meta-DiSc v.1.4 [43] software for meta-analysis of test accuracy data.
Data exploration by principal component analysis was performed using psych, factoextra, and FactoMineR R packages. Univariate and multivariate Cox regression revealed that the multi-miRNAs signature plays a functional role independent of clinicopathological characteristics. Hazard ratio and 95% confidence interval (CI) were estimated. Patients were then divided into low-and high-risk groups by the median expression of a specific risk score formula for predicting tumor progression. The risk score for each patient was calculated based on a linear combination of the miRNA expression level weighted by the regression coefficient derived from the multivariate Cox regression, as follows: In this formula, n represents the number of miRNAs, Coie indicates the coefficient of every miRNA in the result of multivariate Cox regression analysis, and EV i represents the expression level of the miRNA [44,45].
The ability of the miRNA risk score to accurately predict progression events was assessed using Harrell's concordance index (C-statistic) and Brier score [45]. The Wilcoxon rank-sum test compared the c-statistic and Brier score for each of the miRNAs. A C-statistic of 1.0 represents ideal discrimination, indicating the model is ideal for predicting with a greater probability a patient experiencing an event compared with a patient who does not. Alternatively, a C-statistic of 0.5 indicates that the model is no better at classifying outcomes than random chance. As a priori, we set a C-statistic > 0.8 to indicate excellent discrimination, between 0.7 and 0.8 to indicate reasonable or good discrimination, and between 0.5 and 0.7 to indicate poor or weak discrimination. Brier score measures the accuracy of probabilistic predictions to a set of mutually exclusive discrete outcomes. Across all items i ∈ 1... N in a set of N predictions, the Brier score (BS) measures the mean squared difference between (a) the predicted probability assigned to the possible outcomes for item i and (b) the actual outcome Oi. The score ranges from zero to one and represents the square of the largest possible difference between a predicted probability and the actual outcome. Therefore, the lower the Brier score is for a set of predictions, the better the predictions are calibrated. Next, the Brier skilled score (BSS) was calculated to compare the performance of forecasts with that of a reference probability, which is the ATA risk score. Values closer to 1.0 indicate a perfect forecast, while values closer to 0 indicate unreliable forecast probabilities. We set a score of <0.1 to indicate predictive precision >90% [46]. Calculations of C-statistic and BS were performed using the DescTools R package with Cstat, Brier score functions, and 1000 Bootstrap replicates; then, BSS was calculated manually using the following calculations: f t is the forecast probability, O t is the actual outcome of the event at instance t (0 if not happened, 1 if happened), and N is the number of forecasting instances or the number of items for which the Brier score is being calculated). Ref is the BS of reference gold standard test. Fagan's Bayesian nomogram was constructed to plot post-test probability (PP) and likelihood ratios (LR) [47]. It is a graphical tool that allows clinicians to use the results of a diagnostic test to estimate a patient's probability of having the disease. The calculation formula is as follows: Prior probability = probability/(1-probability). Positive likelihood ratio = sensitivity/(1-specificity). Negative likelihood ratio = (1-sensitivity)/specificity. Posterior probability = prior odds * likelihood ratio. We considered likelihood ratio for a positive test (LR+) of more than 10 to significantly increase the probability of disease ("rule in" disease), and for patients who have a negative result, a very low LR-below 0.1-virtually rules out the chance that a person has the disease.
Patients were categorized into high-risk and low-risk groups based on the cutoff value of miRNA risk score at 0.86. Kaplan-Meier plots were applied to illustrate the relationship between high-risk and low-risk groups and survival using the Survminer R package. Log-Rank (Mantel-Cox), Gehan-Breslow-Wilcoxon, and Tarone-Ware tests were applied to investigate the difference in the two curves at different time points. Univariate and multivariate Cox regression models were employed, and a Cox nomogram was constructed using regplot and survival R packages. Two-sided p-values < 0.05 were regarded as significant. Statistical analysis was carried out using IBM SPSS Statistics for Windows, Version 27.0. (IBM Corp. Armonk, NY), GraphPad prism v9.1.1 software (GraphPad Software, San Diego, CA, USA), and R software version 3.4.2 (R Foundation).

Target Gene Prediction, Functional Enrichment Analysis, and External Validation
For the three significant miRNAs, Qiagen Ingenuity Pathway Analysis (IPA) software was used to construct causal networks using complex omics data retrieved from thousands of published articles and curated publicly available datasets within the con-

TCGA and Microarray Thyroid Cancer Cohorts
A total of 495 thyroid cancer patients (448 non-recurrent and 47 recurrent) in the TCGA were screened. Their mean age was 47.2 ± 15.7 years, 73.1% (N = 362) were women, and 66.9% (N = 331) were white. Patients in the recurrence group were more likely to be older (p = 0.040) and white (p = 0.018). Higher prevalence of recurrence was found in patients with tumor stage T3/T4 (p = 0.009), distant metastasis (p < 0.001), and TERT mutation (p = 0.011). The median overall survival was 31.0 months (IQR = 17.4-51.9). Upon screening oncologic and clinicopathologic data of TC patients, we found that none of the recurrent patients received radioactive iodine therapy. A systematic search in the GEO database (up to 14 May 2021) revealed a lack of transcriptomic data following radioactive iodine in thyroid cancer patients. Therefore, we could not identify miRNAs with the predictive role of recurrence following RAI treatment using either TCGA or GEO datasets.

Discovery of Candidate Markers Associated with Progression
DIANA-miRPath v.3.0 revealed 469 miRNAs significantly enriched in the thyroid cancer KEGG pathway (KEGG ID: hsa05216) ( Table S2). The dbDEMC database, cancer vs. normal, metastasis vs. non-metastasis, and high-grade vs. low-grade tumors were compared (Table S3). The analysis yielded 367 significant deregulated miRNAs (193 upregulated and 174 downregulated) in cancer vs. normal comparison. Of these, 21 upregulated and 19 downregulated genes were differentially expressed at absolute fold change (FC) greater than 1. As prognostic biomarkers, 349 miRNAs (170 upregulated and 179 downregulated) were significant compared to the advanced tumor stage vs. lower stage. Filtration at FC > 1.0 showed 12 significant upregulated and eight downregulated miRNAs. Of these, six miRNAs were removed as they showed a paradoxical expression profile (Table S4).
The intersection between highly expressed miRNAs (FC > 1.0) with diagnostic and prognostic values and those identified in DIANA-miRPath v.3.0 depicted three common miRNAs (Figure 2A). Both miR-221-3p and miR-222-3p were upregulated in cancer compared to controls and in an advanced stage compared to lower disease stage ( Figure 2B). They were highly enriched in cancer-related pathways as pathways in cancer, proteoglycans in cancer, Hippo signaling pathway, p53 signaling pathway, cell cycle, adherens junction, and cancer-specific KEGG pathways.

Characteristics of Papillary Thyroid Cancer Patients
For the matched cohorts, the mean age at diagnosis was 38.5 years (range 18-80), and 72.1% (N = 49) were women. As demonstrated in Table 1, clinical and pathological characteristics of patients who received postoperative RAI matched those who did not.  Table 2.

MicroRNA Predictive Performance for Progression Following RAI Treatment
In comparison between indolent and progressive tumors, there was a remarkable downregulation of miR-204 in progressive cases (median = −1.7, IQR = −2. In contrast, the ATA risk score did not differ significantly between non-progressive and progressive groups ( Table 2).
In the PCA plot, exploratory data analysis showed clear discrimination between tumor specimens who remained indolent and those that progressed, with higher levels (long arrows) of miR-221 and miR-222 in progressive tumors, while the miR-204 level was higher in indolent samples. The miRNA discrimination ability performed slightly better in patients following radioactive iodine ( Figure 5A,B). Univariate and multivariate Cox regression analyses showed miR-204, miR-221, and miR-222 as independent risk factors for tumor progression (Table 5).     Log-rank test was used to compare high-risk and low-risk groups categorized based on the microRNA risk score above and below 0.86.

Prognostic Value of MicroRNA Risk Score and Nomogram Construction
The miRNA risk score was generated as follows: Fagan's Bayesian nomogram shows that posterior probabilities for tumor progression increased from 53% to 97% (95% CI = 82-100%) if the miRNA risk score exceeded 0.86 and decreased from 53% to 11% (95% CI = 5-25%) if the score was below 0.86 ( Figure 5F).
Compared to the low-risk group, the high-risk group had a threefold increased progression risk (HR = 2.71, 95% CI = 1.86-3.96, p < 0.001). Kaplan-Meier survival curves showed shorter survival times in the high-risk group of patients. In the overall analysis, patients with risk score >0.86 had disease-free survival (17.3 months, 95% CI = 10.06-24.55) compared to those with lower risk score (70.79 months, 95% CI = 59.12-82.45, p < 0.001). Similarly, subgroup analysis of patients who received RAI ablation treatment showed lower survival in the high-risk group (14.7 months, 95% CI = 7.82-21.7) compared to the low-risk group (49.0 months, p < 0.001) ( Figure 5D,E). The risk score performed best for predicting progression in the whole cohort (C-statistics = 0.943, Brier = 0.083) and RAI subgroup (C-statistic = 0.978, Brier = 0.049). However, the scores did not discriminate well for other pathological features and clinical outcomes (Table 6). IQR: interquartile range. A C-statistic of 1.0 represents ideal discrimination, indicating the model is ideal for predicting with a greater probability a patient experiencing an event compared with a patient who does not. Alternatively, a C-statistic of 0.5 indicates that the model is no better at classifying outcomes than random chance. As a priori, we set a C-statistic >0.8 to indicate excellent discrimination, between 0.7 and 0.8 indicates reasonable or good discrimination, and between 0.5 and 0.7 indicates poor or weak discrimination. The Brier score ranges from zero to one and represents the square of the largest possible difference between a predicted probability and the actual outcome. Therefore, the lower the Brier score for a set of predictions, the better the predictions are calibrated. Significant p-values < 0.05, Brier score, and/or C-Statistic results are shown in the table in bold.
For clinical implementation by physicians, we constructed a Cox nomogram to predict progression-free survival within two-and five-years following diagnosis using mi-croRNA risk score, radioactive ablation treatment, and demographic features ( Figure 6A). An example for the interpretation of nomogram and cox regression results is shown in Figure 6B. The concordance index was 0.805 ± 0.037, indicating that the model has good discrimination ability.  (Table S8). Literature screening results are depicted in Table 7. The current nomogram is derived from well-differentiated thyroid cancer cohorts who underwent surgery at a single center. The outcome measured was a post-operative progression. Cox proportional hazard model was applied. (B) Example for using the nomogram. Assumed having a 20-year-old female patient with a body mass index (BMI) of 20 Kg/m 2 , whose tissue microRNA risk score was high at 2.5, and received radioactive iodine (RAI) ablation. Each variable will be scored on its points scale. The scores for all variables are then added to obtain the total score, and a vertical line is drawn from the total-points row to estimate the probability of survival rates within 2 years and 5 years. *** indicates p < 0.001. Experiments were performed in qRT-PCR.

Discovery of the Regulatory Network
The predicted regulatory network showed that upregulated miR-221 and miR-222 and downregulated miR-204 leads to a subsequent cascade promoting thyroid cancer pathway (Figure 7). A systematic review demonstrated multiple deregulated signaling pathways and mechanisms leading to tumor development (Table S9)

Discussion
Despite evidence of the prominent diagnostic, prognostic, and predictive role of miRNAs in DTC, accurate prediction of progression is highly challenging in TC patients. Without a convenient marker to guide an informed decision between radical treatment and continued observation for PTC, there is an urgent need to identify the biological behavior of such tumors. This work serves as an initial proof-of-principle study that the triple miRNA biomarkers (miR-221, -222, and -204) can predict tumor progression following radioactive iodine ablation in well-differentiated TC patients.
The propensity score matching analysis approach was followed in the current study for specimen selection to allow data matching in general baseline factors and establish two similar datasets for investigating the expression of the selected miRNAs. There was no significant difference in microRNA expression in RAI and non-RAI groups, with consistent deregulation of our microRNA panel. Results showed miR-221 and miR-222 upregulation and miR-204 downregulation to exhibit good predictive accuracy for recurrence, even following RAI therapy. Radiation oncology-associated miRNAs were found to modulate cell death and proliferation after irradiation [101]. However, the molecular basis of gene regulation in cells exposed to radioactive iodine is not fully understood. Some molecular markers as the presence of BRAF V600E and TERT promoter mutations strongly indicate loss of iodine uptake rate and impairment of the iodide-metabolizing machinery [102]. However, around half of PTC tumors not harboring these mutations are non-RAI avid, highlighting a complicated mechanism underlying tumor recurrence/persistence following RAI ablation.
Using loss-of-function assays, Li et al. recently confirmed that the oncogenic long non-coding RNA LINC00514 could promote proliferation/migration and invasion and suppress apoptosis of PTC via miR-204-3p/CDC23 axis [118]. As part of the urothelial carcinoma-associated 1 (UCA1)/miR-204/bromodomain containing 4 (BRD4) axis, miR-204 plays an essential role in PTC cell proliferation/invasion and shows potential for therapeutic applications in PTC patients [119]. Interestingly, the latter investigators found that the long non-coding RNA UCA1 and miR-204 could negatively regulate each other, and the UCA1 may compete with BRD4 for miR-204 binding results in downregulating miR-204, promoting BRD4 expression and affecting PTC progression. Collectively, these findings support the prognostic value of miR-204 downregulation in the proposed model in our samples.
Although the present study has a modest sample size due to strict inclusion/exclusion criteria for selecting well-differentiated thyroid cancer patients subjected to postoperative radioactive iodine ablation, it was large enough to achieve significance in the predictive power of each cohort. In addition, using FFPE samples may seem disadvantageous from the point of view of some investigators, but miRNA expression signatures can be obtained with relative ease and stability using quantitative RT-PCR in tumor biopsy tissue [110,120]. Thus, the proposed prognostic risk score can be calculated as a part of the pathology workup.
Our novel three miRNAs panel and nomogram could accurately identify tumors that are likely to acquire more progressive behavior in PTC patients so that improved management strategies can be developed, avoiding unnecessary tissue biopsy and surgical intervention. Assessment of the prognostic value of this panel in a large-scale multicenter prospective setting is recommended to support the clinical utility and validity of this model. Evaluation of the identified panel with other treatment modalities such as irradiation and a combination of neck dissection and postoperative radiotherapy is warranted.

Conclusions
Our predictive panel/nomogram can define which cancers will have an aggressive phenotype, providing a new paradigm for managing patients diagnosed with localized low-risk DTC. Not only would this have an enormous positive impact on our ability to longitudinally monitor thyroid cancer for evidence of disease progression in a prospective clinical trial, but identifying the signature specific for tumor aggressiveness would unravel new pathophysiological mechanisms and open new horizons to tackle cancer with noninvasive diagnostics and innovative new miRNA-based therapeutics.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13184649/s1, Table S1: Tested microRNAs in the study cohorts, Table S2: Enriched microRNAs in thyroid cancer KEGG pathway (hsa05216), Table S3: MiRNA-seq experiments of thyroid carcinoma with different comparisons, Table S4: Meta-profiling of thyroid cancer diagnostic and prognostic microRNAs, Table S5: Pathway enrichment analysis of 104 deregulated microRNAs in thyroid cancer, Table S6: Downregulated miR-204-5p in various cancer types, Table S7: Upregulated miR-221-3p in various cancer types, Table S8: Upregulated miR-222-3p in various cancer types, Table S9: Summary of molecular pathways regulating the study microRNAs in cancer.  Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Tulane University, USA, (protocol code 2020-1636, 12 August 2020) and the Ethical Committee of Suez Canal University, Egypt (protocol code 4344, 11 March 2020). Patient data were obtained from the hospital medical records, which were anonymized and de-identified before analysis. Informed Consent Statement: Patient consent was waived since archived formalin-fixed paraffinembedded specimens were enrolled in this study. Data Availability Statement: Data are available from the corresponding author upon reasonable request and obtaining the approval of the "Office of Technology Transfer and Intellectual Property Development, Tulane University, USA".