Progression Risk Assessment of Post-surgical Papillary Thyroid Carcinoma Based on Circular RNA-Associated Competing Endogenous RNA Mechanisms

Background: Accurate risk assessment of post-surgical progression in papillary thyroid carcinoma (PTC) patients is critical. Exploring key differentially expressed mRNAs (DE-mRNAs) regulated by differentially expressed circular RNAs (circRNAs) via the ceRNA mechanism could help establish a novel assessment tool. Methods: ceRNA network was established based on differentially expressed RNAs and correlation analysis. DE-mRNAs within the ceRNA network associated with progression-free interval (PFI) of PTC were identified to construct a prognostic ceRNA regulatory subnetwork. least absolute shrinkage and selection operator (LASSO)–Cox regression was applied to identify hub DE-mRNAs and establish a novel DE-mRNA signature in predicting PFI of PTC. Results: Six hub DE-mRNAs, namely, CLCNKB, FXBO27, FXYD6, RIMS2, SPC24, and CDKN2A, were identified to be most significantly related to the PFI of PTC, and a prognostic DE-mRNA signature was proposed. A nomogram incorporating the DE-mRNA signature and clinical parameters was established to improve the progression risk assessment in post-surgical PTC, which was superior to the American Thyroid Association risk stratification system and distant Metastasis, patient Age, Completeness of resection, local Invasion, and tumor Size (MACIS) score American Joint Committee on Cancer staging system. Conclusions: Based on the circRNA-associated ceRNA RNA mechanism, a DE-mRNA signature and prognostic nomogram was established, which may improve the progression risk assessment in post-surgical PTC.


INTRODUCTION
Thyroid cancer, originating from thyroid follicular epithelial cells, is the most common malignancy of the endocrine system (La Vecchia et al., 2015). Thyroid cancer is currently the eighth most common malignancy around the world (Antonelli and La Motta, 2017), and its incidence is rapidly increasing worldwide. There are four main types of thyroid cancer, including papillary thyroid carcinoma (PTC), follicular thyroid carcinoma, anaplastic thyroid carcinoma, and medullary thyroid carcinoma. Of these, PTC accounts for more than 80% of all thyroid cancer cases. Among all types of thyroid cancer, PTC has the highest incidence with a higher proportion of young patients. Genetic mutations and environmental exposures are the major risk factors for PTC in which BRAF V600E is the most common mutation. The primary treatment for PTC is thyroidectomy. Most PTCs grow slowly, and 90% of patients have a relatively good prognosis after treatment, with a low post-operative recurrence rate <10% (Haugen et al., 2016). However, a small percentage of patients experience a high risk of recurrence and metastasis after surgery, resulting in a poor prognosis.
The molecular mechanisms leading to PTC recurrence remain unclear. Moreover, there is a lack of satisfactory assessment tools to precisely evaluate the post-surgical risk of recurrence in PTC. The American Thyroid Association (ATA) currently recommends the use of American Joint Committee on Cancer (AJCC) staging system and distant Metastasis, patient Age, Completeness of resection, local Invasion, and tumor Size (MACIS) score to predict mortality after PTC and proposed the risk stratification system for risk assessment of postsurgical progression in PTC (Haugen et al., 2016). However, the accuracy of currently available assessment tools is inadequate to meet the requirements of clinical work. The challenge of treating PTC lies in balancing side effects and treatment benefits (Cabanillas et al., 2016). Patients with low-and high-risk PTC should be differentiated to adopt different treatment strategies. Patients with high risk of recurrence require more extensive surgical resection and adjuvant radioactive iodine therapy to improve post-surgical prognosis and prevent recurrence. Secondary surgery for recurrence can result in additional surgical trauma with a higher risk of recurrent laryngeal nerve injury. Alternatively, PTC patients with low risk of recurrence should receive a more conservative treatment to reduce surgical trauma and improve quality of life. Unnecessary post-surgical radioactive iodine therapy and thyroid-stimulating hormone (TSH) suppression therapy should also be avoided. Longterm subclinical hyperthyroidism caused by a high dose of TSH suppression can lead to multiple potential side effects, including osteoporosis, atrial fibrillation, heart discomfort, and an increased risk of fractures and heart disease in elderly patients (Biondi and Cooper, 2018). Unnecessary radioactive iodine therapy also increases the risk of developing other cancers. Therefore, more accurate assessment tools to assess the risk of progression in post-surgical PTC patients are urgently required.
Non-coding RNAs such as microRNAs (miRNAs), long non-coding RNAs, and circular RNAs (circRNAs) are major components of the human transcriptome. circRNAs were first discovered in human cells in 1986 and are critical in the pathophysiology of several human diseases including cancer (Kos et al., 1986). circRNAs are products of a reverse splicing mechanism and are resistant to RNase R (Chen and Yang, 2015). circRNAs contain a single-stranded covalent closedloop structure that have neither a 5 ′ -3 ′ polarity nor a polyadenylated tail. circRNAs are primarily located in the cytoplasm and regulate gene expression by acting as an miRNA sponge at the transcriptional or post-transcriptional stage, affecting protein translation and thus the regulation of cellular processes. According to the competitive endogenous RNA (ceRNA) hypothesis, non-coding RNAs such as circRNAs can positively regulate mRNA expression by competitive binding of its shared miRNAs (Salmena et al., 2011). Recent studies further reveal that several circRNAs differentially expressed in cancers can exert a tumor-promoting or tumorsuppressing role via ceRNA mechanism, thus regulating apoptosis, proliferation, invasion, and metastasis (Zhong et al., 2018;Liu et al., 2020). mRNAs regulated by differentially expressed circRNAs (DE-circRNAs) may be more likely to play a functional role in cancer. Particularly in PTC, the upregulated SOX2 by circ_0005273/hsa-miR-1183 axis promotes the proliferation, invasion, and metastasis of PTC . Transforming growth factor α, upregulated by hsa-circ-0060060/hsa-miR-144-3p axis in PTC, was identified to promote the proliferation and autophagy of PTC and was associated with poor prognosis . The establishment of circRNA-miRNA-mRNA regulatory networks in PTC may provide key targets and novel regulatory pathways in understanding the progression of PTC.
The primary objective of this study was to explore key differentially expressed mRNAs (DE-mRNAs) regulated by DE-circRNA-associated ceRNA regulatory mechanism in PTC to elucidate the post-surgical progression of PTC and establish a novel assessment tool.

Identification of DE-circRNAs, DE-miRNAs, and DE-mRNAs in PTC
Expression data of circRNAs and related clinical data for PTC were downloaded from the GSE93522 dataset based on GPL19978 platform (Peng et al., 2017). An annotation file provided by the manufacturers was used to match the probes with corresponding circRNA IDs. The median ranking value was used to determine the expression value if multiple probes matched a single circRNA ID. The 'LIMMA' R package was used to identify DE-circRNAs between cancer and normal tissues (Ritchie et al., 2015). The cutoff value was set at |Log2FC| > 1 and P < 0.05. The basic characteristics of DE-circRNAs were obtained from circBase (https://www.circbase. org/), and the corresponding structures were analyzed using the Cancer-Specific circRNA Database (CSCD) (https://gb.whu.edu. cn/CSCD/#/).
Normalized RNA sequencing data in the form of millions of transcripts (TPM) for miRNA and mRNA and related clinical data of post-surgical PTC were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/, up to July 1, 2020). TCGA-Thyroid Cancer (THCA) dataset originally included 507 cases, 510 tumor samples, and 58 normal tissue samples. A total of 495 cases of PTC with matching tumor tissue and clinical information, and 58 normal samples of thyroid tissue were eventually included in the current analysis following the removal of five cases without matching tumor samples, two cases with poorly differentiated oncolytic carcinoma or follicular carcinoma, five cases with history of neoadjuvant therapy, and eight samples of metastases. The 'LIMMA' R package was used to identify DE-miRNAs and DE-mRNAs between cancer and normal tissues with a cutoff value set at | Log2FC | > 1 and a false discovery rate (FDR) < 0.05 (Ritchie et al., 2015). GEPIA (http://gepia.cancer-pku.cn) incorporating RNA sequencing expression data of tumors and normal samples from TCGA and the Genotype-Tissue Expression projects was used to validate the differential expression level of a specific DE-mRNA (Tang et al., 2017). Data of copy number alterations and mutation were downloaded from Cbioportal (http://www. cbioportal.org/).

Construction of the circRNA-miRNA-mRNA Network
CIRCinteractome is a computational tool that enables the prediction and mapping of binding sites for RNA-binding proteins and miRNAs on reported circRNAs (Dudekula et al., 2016). circRNA-miRNA relationships were predicted using CIRCinteractome. miRWalk is an open-source platform that provides predicted and validated miRNA-binding sites with a machine learning algorithm, incorporating data from TargetScan (conserved site context scores, version 7.1), miRDB (release 5.0), and validated information from miRTarBase (version 7.0) (Sticht et al., 2018). miRNA-mRNA relationships were predicted using miRWalk 3.0 with score ≥0.95. Based on differential expression data, only relationship pairs with negative correlation were retained. The correlation between DE-miRNA and DE-mRNA expression with potential ceRNA regulatory relationships was further analyzed using Pearson correlation coefficient in TCGA-THCA dataset. Only miRNA-mRNA relationships with r < −0.2 and P < 0.05 were retained in the final circRNA-miRNA-mRNA network. To construct a ceRNA regulatory subnetwork associated with progression-free interval (PFI) of PTC, univariate Cox regression analysis was performed based on the DE-mRNAs involved. Only DE-mRNAs associated with PFI of PTC and the corresponding circRNA-miRNA and miRNA-mRNA pairs were retained in the PFI-related subnetwork. The ceRNA network was visualized using Cytoscape v.3.8.0 (https://www. cytoscape.org/). The Sankey diagram that represented the ceRNA regulatory relationship was drafted using Origin 2020 (https:// www.originlab.com/).

Bioinformatics Analysis of the circRNA-miRNA-mRNA Network
To investigate the underlying biological function of ceRNA regulatory relationships in PTC, DAVID (https://david.ncifcrf. gov/) was employed to explore enriched biological processes, cellular components, molecular functions, and significantly relevant signal pathways of DE-mRNAs involved, using default parameters (Huang da et al., 2009). P < 0.05 was regarded statistically significant. The results were visualized using the OmicShare tools, a free online platform for data analysis (http:// www.omicshare.com/tools).

Identification of PFI-Related Hub DE-mRNAs and Establishment of the Prognostic DE-mRNA Signature
In this study, PFI was selected as the main endpoint (https://wiki. nci.nih.gov/plugins/servlet/mobile#content/view/24279961). DE-mRNAs associated with PFI were identified based on TCGA-THCA dataset using a univariate Cox proportional hazards regression model. Normalized gene expression data were transformed on the base-2 logarithm for further survival analysis. DE-mRNAs with P < 0.05 were considered statistically significant for further analysis. Only cases with follow-up >30 days were included for survival analysis. The 492 eligible TCGA cases were subsequently randomly divided into a training dataset and a validation dataset in a 7:3 ratio. Least absolute shrinkage and selection operator (LASSO) penalized Cox regression analysis was performed to select prognostic hub DE-mRNAs related to PFI and constructed a prognostic DE-mRNA signature in patients with PTC based on a linear combination of the regression coefficient derived from LASSO-Cox regression model coefficients (β) multiplied with its normalized mRNA expression level in the training dataset. Patients were divided into the high-risk and the low-risk group based on the optimal cutoff value determined by X-Tile (Camp et al., 2004). Performance of the prognostic DE-mRNA signature was assessed using Harrell's concordance index (C-index), area under the curve (AUC) of the receiver operating characteristic (ROC) curve, and Kaplan-Meier analysis. The validation dataset and TCGA-THCA dataset were further utilized for validation.

Gene Set Enrichment Analysis
Gene Set Enrichment Analysis (GSEA) was performed to compare the molecular alteration in the high-risk group and the low-risk group against previously annotated gene sets (Subramanian et al., 2005). Samples from TCGA-THCA dataset were divided into the high-risk group and the low-risk group according to the optimal cutoff value determined by X-Tile. Thereafter, GSEA was run on javaGSEA v4.0.3 based on the Molecular Signatures Database v7.1. H: Hallmark gene sets, C2: curated gene sets and C5: gene otology gene sets were searched to identify enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, biological processes, cellular components, molecular functions, and specific well-defined biological states or processes associated with poorer survival of the high-risk group. Gene sets with |NES| >1 and FDR <0.25 were considered statistically significant.

Validation of the Independent Prognostic Role of the DE-mRNA Signature
To validate the independent prognostic role of the DE-mRNA signature, univariate and multivariate Cox regression analyses were performed on the prognostic DE-mRNA signature and clinical parameters, including age; PFI status; mutational status of BRAF V600E, RAS, RET, NTRK1, and TERT; sex; histological subtype; T stage; N stage; M stage; AJCC stage; residual tumor status; extrathyroidal extension; multifocality; and anatomic site in TCGA dataset. Parameters with P < 0.25 in univariate analysis were further included in the multivariate Cox regression analysis. P < 0.05 was considered statistically significant.

Building and Validation of a Prognostic Nomogram
Following a test for collinearity, independent prognostic parameters and relevant clinical parameters were included to construct a prognostic nomogram to predict 1-, 2-, 3-, 4-, and 5-year PFI of PTC patients in the entire TCGA dataset with the stepwise Cox regression model. Discrimination of the nomogram was assessed with the C-index calculated by a bootstrap method with 1,000 resamples. The predicted progression-free survival of nomogram against observed survival rates was plotted using the calibration curve. The performance of the prognostic nomogram was further assessed with ROC curves and Kaplan-Meier analysis. The ATA risk stratification system, MACIS score, and AJCC staging system were used as references for comparison of AUCs. Information of ATA risk stratification and MACIS score for each case was obtained from the official TCGA publication (Cancer Genome Atlas Research, 2014).

Statistical Analysis
Statistical analysis was conducted using R software v3.4.3 (https:// www.r-project.org/), and GraphPad Prism v8.01 (https://www. graphpad.com/). Categorical variables were analyzed with χ 2test or Fisher exact test. Continuous variables of two groups were analyzed with Student t-test. Multiple groups of continuous variables were analyzed with a one-way analysis of variance test. DE-mRNAs associated with progression-free survival were identified based on univariate and multivariate Cox regression analysis. Hazard ratio (HR) and 95% confidence interval (CI) were calculated. Correlation between two variables was analyzed with Pearson correlation coefficient. ROC curves were analyzed using the "timeROC" R package, and the AUCs were compared with the method described by Blanche et al. (2013). A two-tailed P < 0.05 was considered statistically significant if there was no special statement.

Identification of DE-circRNAs, DE-miRNAs, and DE-mRNAs in PTC
Workflow of the current study is shown in Figure 1A. DE-circRNAs, miRNAs, and  Prediction of Potential Circular RNA-Associated ceRNA Regulatory Relationships and Functional Enrichment Analysis circRNA-miRNA relationships were predicted using Circular RNA Interactome. miRNA-mRNA relationships were predicted using mirWalk 3.0. Based on the targeted relationship combined with differential expression data, 93 potential circRNA-miRNA pairs and 1,105 miRNA-mRNA pairs were identified, which involved 54 DE-circRNAs, 14 DE-miRNAs, and 698 DE-mRNAs ( Figures 1B-D and Supplementary Tables 4, 5).
To investigate the underlying biological function of ceRNA regulatory relationships in PTC, functional enrichment analysis was performed based on the 698 DE-mRNAs involved (Figures 2A-D and Supplementary Table 6). Biological processes associated with malignant phenotype were significantly enriched, which included regulation of cell adhesion, extracellular matrix disassembly and organization, positive regulation of cell proliferation, and angiogenesis. Processes associated with intracellular signal transduction were also significantly enriched, particularly the positive regulation of the mitogen-activated protein kinase cascade. In terms of KEGG pathways, DE-mRNAs were significantly associated with transcriptional misregulation and miRNAs in cancer. Multiple tumor-associated signaling pathways were also enriched, including p53 and PI3K-Akt signaling pathways.

Construction of the circRNA-miRNA-mRNA Network
The correlation between DE-miRNA and DE-mRNA expression with potential ceRNA regulatory relationships in PTC was further analyzed using Pearson correlation coefficient. Only miRNA-mRNA relationships with negatively correlated expression were retained. Correlation analysis of miRNAs and potential target mRNAs in TCGA-THCA dataset are shown in Supplementary Table 7. A ceRNA network was constructed, which included 25 DE-circRNAs, 6 DE-miRNAs, and 150 DE-mRNAs based on 32 circRNA-miRNA pairs and 153 miRNA-mRNA pairs (Figure 3).
To construct a ceRNA regulatory subnetwork associated with PFI in PTC, univariate Cox regression analysis was performed based on the 150 DE-miRNAs involved. A total of 27 DE-mRNAs were identified to be associated with PFI of PTC. Following retention of the corresponding regulatory circRNA-miRNA and miRNA-mRNA pairs, a ceRNA regulatory subnetwork associated with PFI in PTC was constructed. This included 11 DE-circRNAs, 3 DE-miRNAs and 27 DE-mRNAs based on 11 circRNA-miRNA pairs, and 27 miRNA-mRNA pairs ( Figure 4A). hsa_circ_0003645, hsa_circ_0089153, hsa_circ_0005699, hsa_circ_0007146, hsa_circ_0038718, hsa_circ_0001658, hsa_circ_0008784, hsa_circ_0000965, hsa_circ_0001917, hsa_circ_0008354, and hsa_circ_0049271 were DE-circRNAs within the PFI-related subnetwork. The basic characteristics of the 11 DE-circRNAs are shown in Table 1. The corresponding structures of DE-circRNAs were analyzed using the CSCD and shown in Figure 4B.  To identify hub PFI-related DE-mRNAs within the ceRNA subnetwork, 492 TCGA cases with follow-up >30 days were randomly divided into the training dataset and the validation dataset in a 7:3 ratio ( Table 2). The training dataset was employed to identify hub DE-mRNAs related to PFI using LASSO-Cox regression. Subsequently, six DE-mRNAs, including CLCNKB, FXBO27, FXYD6, RIMS2, SPC24, and CDKN2A, were identified as hub PFI-related DE-mRNAs (Figures 5A,B). ceRNA regulatory relationships of the six hub DE-mRNAs are shown in Figure 5C. Downregulated CLCNKB, FBXO27, and FXYD6 with HR <1 were negatively regulated by hsa-miR-146b-3p. Upregulated RIMS2 with HR >1 was negatively regulated by hsa-miR-139-5p. Moreover, upregulated SPC24 and CDKN2A with HR  Figure 6A). Thereafter, patients were divided into a high-risk and a low-risk group according to the optimal cutoff value determined by X-tile. Kaplan-Meier analysis revealed that high-risk patients with significantly poorer prognosis could be discriminated by the DE-mRNA signature   ( Figure 6D). As the incremental increase of risk score and the deterioration of prognosis occurred, the expression of FXYD6, CLCNKB, and FBXO27 gradually decreased, and the expression of RIMS2, SPC24, and CDKN2A gradually increased ( Figure 6G).
In general, the prognostic DE-mRNA signature based on hub PFI-related DE-mRNAs within the ceRNA subnetwork was accurate in predicting PFI of PTC in the training dataset.  Figure 6C). Thereafter, patients in each dataset were divided into a high-risk and a low-risk group according to the optimal cutoff values determined by X-tile. Kaplan-Meier analysis revealed that high-risk patients with significantly poorer prognosis could be discriminated by the DE-mRNA signature in both datasets (Figures 6E,F). Moreover, heatmaps revealed that the expression of the six DE-mRNAs gradually changed along with risk score increase and prognosis deterioration (Figures 6H,I). The prognostic DE-mRNA signature based on the hub PFI-related DE-mRNAs within the ceRNA subnetwork was confirmed to be robust in predicting PFI of PTC in the validation dataset and the entire TCGA-THCA dataset.

Analysis of Molecular, Clinical, and Mutational Relevance of the DE-mRNA Signature
To explore the underlying molecular mechanism of the DE-mRNA signature, the entire TCGA-THCA dataset was divided into a high-risk and a low-risk group based the optimal cutoff value of the DE-mRNA signature determined by X-tile. GSEAs were utilized to compare molecular alteration in the high-risk group and the low-risk group against previously annotated gene sets in the Molecular Signatures Database v7.1. Regarding KEGG pathways, molecular alteration in the high-risk group was significantly enriched in the p53 signaling pathway, cell cycle, etc. (Figures 7A,B). Analysis of hallmark gene sets further revealed that molecular alteration in the high-risk group was significantly associated with genes encoding cell cycle-related targets of E2F transcription factors ( Figure 7C). HALLMARK_MYC_TARGETS_V1 and HALLMARK_MYC_TARGETS_V2 was also significantly enriched indicating that the molecular alteration in the high-risk group was potentially regulated by MYC (Figures 7D,E). Full results of the GSEA are shown in Supplementary Table 8.
Analysis of clinical relevance revealed that the risk score increased significantly along with the escalation of T stage (P < 0.05; Figure 7F). Patients with lymph node and distal metastasis had a significantly higher risk score (P < 0.05; Figures 7G,H). Regarding AJCC stage, advanced stages (III and IV) were  associated with a higher risk score than early stages (I and II) (P < 0.05; Figure 7I). Patients with microscopic and macroscopic residual tumors also had a significantly higher risk score, further indicating the relevance of the DE-mRNA signature to the aggressiveness of PTC (P < 0.05; Figure 7J). Analysis of histological subtypes revealed that aggressive subtypes of PTC, including tall cell, hobnail variant, and columnar cell carcinoma, were associated with a significantly higher risk score (P < 0.05; Figure 7K). Analysis of mutational relevance revealed that the BRAF V600E mutation and RAS mutation was associated with a significantly higher and lower risk score, respectively (P < 0.05; Figures 7L,M). The mutational pattern of PTC further revealed that CLCNKB, FBXO27, and FXYD6 had a relatively lower expression level in samples with the BRAF V600E mutation and were negatively associated with hsa-miR-146b-3p expression ( Figure 7N). Alternatively, SPC24 and CDKN2A negatively regulated by hsa-miR-139-3p had a relatively higher expression level in samples with the BRAF V600E mutation. RIMS2 also had a relatively higher expression level in samples with the BRAF V600E mutation and was negatively associated with hsa-miR-139-5p expression. Regarding the RAS mutation, the six DE-mRNAs and the corresponding miRNA regulators had opposite expression patterns to the BRAF V600E mutation.

Validation of the Independent Prognostic Role of the DE-mRNA Signature in PTC
A total of 389 cases from the entire TCGA-THCA dataset with complete clinical information were utilized to identify prognostic factors associated with PFI in PTC. Clinical information for analysis included DE-mRNA signature, age, PFI status, mutational status of BRAF V600E, RAS, RET, NTRK1, and TERT, sex, histological subtype, T stage, N stage, M stage, AJCC stage, residual tumor status, extrathyroidal extension, multifocality, and anatomic site. Reasons for exclusion are listed in Supplementary Table 9. Baseline characteristics of patients who were included for prognostic factor evaluation are shown in Table 3. Patients in the high-risk group had significantly worse prognosis; higher T stage, N stage, M stage, and AJCC stage; and a larger portion with the BRAF V600E mutation and extrathyroidal extension. Tumors of the high-risk group were also more likely to be unifocal. Univariate Cox regression analysis revealed that DE-mRNA signature, age, histological type, aggressive subtype, T stage, N stage, AJCC stage, and neoplasm largest dimension were prognostic factors associated with PFI of PTC (P < 0.05; Table 4). Adjustment for parameters with P < 0.25 based on univariate analysis and multivariate Cox regression analysis further validated that the DE-mRNA signature was an independent prognostic factor associated with PFI in PTC ( Table 5). The   ) and corresponding DE-miRNAs (hsa-miR-146b-3p, hsa-miR-139-3p, and hsa-miR-139-5p) and the mutational profiles (BRAF, RAS, RET, NTRK1, and TERT) of PTC. Data of genetic alteration were obtained from the cBioPortal for Cancer Genomics (https://www.cbioportal.org). *P < 0.05, **P < 0.01, and ****P < 0.0001.
RAS mutation was also identified to be an independent prognostic factor.

Building and Validation of a DE-mRNA Signature-Based Prognostic Nomogram
Based on cases from the entire TCGA-THCA dataset with complete clinical information, a prognostic nomogram was established with the DE-mRNA signature and prognostic clinical factors using stepwise Cox regression ( Figure 8A). Factors involved in the nomogram included the DE-mRNA signature, RAS mutation status, aggressive subtype, N stage, and tumor size. The C-index of the nomogram in predicting PFI of PTC was 0.790 (95% CI, 0.652-0.927). The calibration curve revealed that the predicted risk of progression by the nomogram was close to the observed risk of progression ( Figure 8C). Patients were subsequently divided into the high-risk and the low-risk group according to the optimal cutoff value determined by Xtile. Kaplan-Meier analysis revealed that high-risk patients had significantly poorer prognosis than low-risk patients ( Figure 8B). The ATA risk stratification system is recommended to estimate the risk of post-surgical recurrent disease in PTC. MACIS score and AJCC staging system are employed to predict post-surgical mortality of PTC. The accuracy of nomogram in predicting PFI of PTC was analyzed with ROC curves using the ATA risk stratification system, MACIS score, and AJCC staging system as reference (Figures 8D-H

DISCUSSION
Accumulating evidence has shown that circRNAs and their mediated ceRNA regulatory networks play a vital role in pathogenesis and progression of various human cancers, including PTC. hsa-circ-0058124 was previously identified to be upregulated in PTC and associated with poor prognosis. By sponging hsa-miR-218-5p, hsa-circ-0058124 promoted proliferation, invasion, and metastasis of PTC via the NOTCH3/GATAD2A signaling axis (Yao et al., 2019). Through regulation of ABCA9 and MTA1 via sponge hsa-miR-1179 and hsa-miR-1205, upregulated hsa_circ_0039411 promoted proliferation, migration, and invasion of PTC cells and inhibited cell apoptosis (Yang et al., 2020). CircRNAs are resistant to exonase and RNase R, enabling them to be more stable than other types of ncRNAs, as well as play an important regulatory role in cancer cell. The current understanding of the circRNA-related ceRNA network in PTC is limited, requiring further study. In this study, DE-circRNAs, DE-miRNAs, and DE-mRNAs with potential ceRNA regulatory relationships in PTC were identified through comprehensive analysis of targeting relationship prediction and differential expression data. Based on the further correlation analysis of expression between DE-miRNAs and DE-mRNAs with potential ceRNA regulatory relationships, a circRNA-miRNA-mRNA regulatory network was constructed in PTC. DE-mRNAs within the ceRNA network associated with PFI of PTC were further identified to construct a ceRNA regulatory subnetwork associated with PFI. Six hub DE-mRNAs, namely, CLCNKB, FXBO27, FXYD6, RIMS2, SPC24, and CDKN2A, were identified to be most significantly related to the PFI of PTC, which were regulated by three DE-miRNAs, including hsa-miR-146b-3p, hsa-miR-139-5p, and hsa-miR-139-3p. They were subsequently regulated by 11 DE-circRNAs via the ceRNA mechanism, including hsa_circ_0003645, hsa_circ_0089153, hsa_circ_0005699, hsa_circ_0007146, hsa_circ_0038718, hsa_circ_0001658, hsa_circ_0008784, hsa_circ_0000965, hsa_circ_0001917, hsa_circ_0008354, and hsa_circ_0049271. These newly identified key circRNA-miRNA-mRNA regulatory relationships may help elucidate the molecular mechanism of progression in post-surgical PTC. Such insights may provide novel therapeutic targets for treatment of recurrence in PTC. Hub DE-mRNAs identified in the ceRNA network were potential predictors of PFI in PTC.
Although most patients with PTC have a relatively good prognosis, a portion of patients will eventually develop postsurgical recurrence. These high-risk PTC patients should be distinguished early enough to adopt more aggressive treatment options, additional adjuvant radioactive iodine therapy, thyroid  hormone suppression therapy with a higher dosage, and timely intervention to prevent post-surgical progression if necessary. In contrast, for the remaining low-risk PTC patients, aggressive treatment options should be avoided to improve quality of life. Therefore, it is critical to accurately predict the post-surgical risk of progression in patients with PTC. Currently, the ATA guidelines recommend the use of AJCC staging system and MACIS score to predict the risk of post-operative mortality. ATA risk stratification system was recommended to assess postsurgical risk of recurrence. The existing risk assessment tools are not sufficiently accurate, and a novel prediction system should be established to more accurately predict the prognosis of PTC patients. In the current study, we established a circRNAassociated ceRNA network in the PTC. Based on this, six hub PFI-related DE-mRNAs of the PTC were identified. A six-DE-mRNA signature was subsequently established, which was able to distinguish high-risk PTC patients from the low-risk ones and could accurately predict the PFI. Nomogram has been widely applied in oncology to assess the prognosis of cancer patients (Balachandran et al., 2015). Nomogram can integrate various prognostic factors, including molecular and clinical parameters, and provide a visual graphical interface for personalized prediction of clinical events. In this study, to establish a more accurate assessment tool for prognosis evaluation in post-surgical PTC, a prognostic nomogram was established incorporating the DE-mRNA signature and clinicopathological parameters. This  nomogram was able to accurately predict PFI of PTC and was better than the ATA risk stratification system, MACIS score, and AJCC staging system recommended by ATA guidelines. Currently, 11 DE-circRNAs were identified to regulate the six hub DE-mRNAs via the ceRNA mechanism in the study. The role of hsa_circ_0089153 in PTC has previously been reported. hsa_circ_0089153 was upregulated in clinical specimens of PTC . In vitro experiments indicated that downregulation of hsa_circ_0089153 expression could inhibit proliferation, migration, and invasion of PTC cells. The luciferase reporter assay confirmed that hsa_circ_0089153 sponged hsa-miR-145 to mediate upregulation of ZEB2 expression, thereby playing a carcinogenic role in PTC. hsa_circ_0089153 was also reported to be upregulated in gastric cancer and bladder cancer (Wei et al., 2019;Zhuang et al., 2020). Our study indicated that upregulated hsa_circ_0089153 may sponge hsa-miR-139-3p and upregulated SPC24 and CDKN2A expression to promote PTC progression. hsa_circ_0003645 was previously identified to be upregulated in non-small cell lung cancer tissue . hsa_circ_0038718 was reported to be upregulated in hepatocellular carcinoma (Qiu et al., 2019). Our result suggested that hsa_circ_0003645 and hsa_circ_0038718 may also play an oncogenic role as a sponge of hsa-miR-139-3p. The function of hsa_circ_0001917 in PTC has not been reported. However, hsa_circ_0001917 was also identified to be upregulated in hepatocellular carcinoma (Qiu et al., 2019). Our result suggested that the upregulated  Six hub DE-mRNAs markedly associated with the PFI of PTC were identified in the current study. Upregulated SPC24 and CDKN2A were identified as targets of hsa-miR-139-3p. The tumor-suppressing role of hsa-miR-139-3p has been identified in multiple tumors (Sannigrahi et al., 2017). SPC24 is an important component of kinetochore-associated NDC80 complex. This complex mediates chromosome segregation and spindle checkpoint activity. SPC24 plays an important role in maintaining the integrity of kinetochore (McCleland et al., 2004). SPC24 was previously identified to be highly expressed in anaplastic thyroid cancer. Knockdown of SPC24 expression inhibited cell growth and invasion and promoted tumor cell apoptosis (Yin et al., 2017). SPC24 was also reported to be highly expressed in hepatocellular carcinoma and was an independent predictor of survival (Zhu et al., 2015). Downregulation of SPC24 inhibited growth and invasion of tumor cells and promoted apoptosis. The oncogenic role of SPC24 was also identified in breast cancer and lung cancer (Zhou et al., , 2018. Here, SPC24 was identified to be negatively regulated by hsa-miR-139-3p. This regulatory role between SPC24 and hsa-miR-139-3p was previously observed in bladder cancer (Yonemori et al., 2016). CDKN2A is traditionally known as a tumor-suppressor gene coding for two proteins, including the p16INK4a and p14arf (Bockstaele et al., 2006). CDKN2A is involved in cell cycle regulation. However, mounting data suggest that CDKN2A may play a dual role in multiple tumors. The p14arf that it codes plays an important role in invasion and metastasis and is associated with a poor prognosis (Fontana et al., 2019). Upregulation of p14arf has been identified in multiple hematological malignancies, aggressive types of Bcell lymphomas, and bladder cancers (Sánchez-Aguilera et al., 2002;Berggren et al., 2003;Lee et al., 2003). The oncogenic function of p14arf is associated with the autophagy regulation (Fontana et al., 2019). Downregulation of p14arf was shown to inhibit progression of lymphomas with the MYC mutation (Humbey et al., 2008). Particularly in PTC, p16INK4a, and p14arf coded by CDKN2A were both identified to be upregulated in thyroid tumorigenesis (Ferru et al., 2006). Wild-type p14arf has The calibration plot for internal validation of the nomogram. The X axis represents the predicted risk of progression, whereas the Y axis represents the observed risk of progression. (D) The time-dependent ROC curves for 1-year PFI predicted by the nomogram in comparison with the ATA risk stratification system, the MACIS score, and the AJCC staging system. (E) The time-dependent ROC curves for 2-year PFI predicted by the nomogram in comparison with the ATA risk stratification system, the MACIS score, and the AJCC staging system. (F) The time-dependent ROC curves for 3-year PFI predicted by the nomogram in comparison with the ATA risk stratification system, the MACIS score, and the AJCC staging system. (G) The time-dependent ROC curves for 4-year PFI predicted by the nomogram in comparison with the ATA risk stratification system, the MACIS score, and the AJCC staging system. (H) The time-dependent ROC curves for 5-year PFI predicted by the nomogram in comparison with the ATA risk stratification system, the MACIS score, and the AJCC staging system. AUCs were compared with the method described by Subramanian et al. (2005).
been observed to delocalize into the cytoplasm in aggressive PTC. Here, we further identified that CDKN2A was upregulated in PTC and was associated with shorter PFI, and hsa-miR-139-3p was the potential negative regulator of CDKN2A. The function of CDKN2A in PTC progression requires further study. RIMS2 was also identified to be upregulated in PTC and was potentially regulated by hsa-miR-139-5p. hsa-miR-139-5p was identified to be downregulated in the primary tumor and further in PTC metastasis (Montero-Conde et al., 2020). hsa-miR-139-5p was also able to be sponged by circBACH2 and relieved the suppression of the target gene LMO4 in PTC (Cai et al., 2019). RIMS2 was identified as a novel target of hsa-miR-139-5p in this study. RIMS2 is a Rab effector and scaffold protein associated with exocytosis (Yoo et al., 2013). It was recently reported to have a high mutation rate in melanoma and mantle cell lymphoma (Hill et al., 2020;Zhang and Xia, 2020). In this study, RIMS2 was identified as a hub DE-mRNA in the ceRNA network and was associated with PFI of PTC. The function of RIMS2 in PTC requires further experimental validation. In this study, downregulated CLCNKB, FBXO27, and FXYD6 were identified to be novel targets of hsa-miR-146b-3p in PTC. hsa-miR-146b-3p was previously reported to be upregulated in PTC and positively associated with central lymph node metastases (Han et al., 2016). hsa-miR-146b-3p may promote invasion and metastasis of PTC by targeting NF2 (Yu et al., 2018). hsa-miR-146b-3p also targeted PAX8 in modulating the differentiated phenotype of PTC (Riesco-Eizaguirre et al., 2015). FXYD6 is known as a specific modulator of Na and K-ATPase and is expressed in multiple epithelial cells of the inner ear (Delprat et al., 2007). In accordance with the current study, FXYD6 was previously identified to be downregulated in PTC and was associated with poor prognosis . CLCNKB is a voltage-gated chloride channel participating in the regulation of cell volume, membrane potential stabilization, signal transduction, and transepithelial transport (Estévez et al., 2001). CLCNKB was previously reported to be downregulated in renal carcinoma (Murakami et al., 2007). Hypermethylation in deletions of CLCNKB in renal carcinoma further indicated its tumor-suppressing role in cancer (Girgis et al., 2012). FBXO27 is a component for substrate recognition in the SCF-type E3 ubiquitin ligase complex (Glenn et al., 2008). Its role in cancer is currently unknown. In this study, CLCNKB, FBXO27, and FXYD6 were identified to be negatively regulated by hsa-miR-146b-3p. Together, hsa_circ_0008354 and hsa_circ_0049271 formed an important part of the PFI-related ceRNA subnetwork we established. The potential tumorsuppressing role of CLCNKB, FBXO27, and FXYD6 in PTC requires further validation.
To the best of our knowledge, the ceRNA network we established and the prognostic DE-mRNA signature we proposed have not been reported previously. The nomogram incorporating the DE-mRNA signature and clinical parameters was robust in predicting PFI of PTC. DE-circRNAs, DE-miRNAs, and DE-mRNAs were potential therapeutic targets for prevention and treatment of recurrent PTC. We acknowledge that our study inevitably has some limitations. First, sequencing data and follow-up data of PTCs in our study were based on TCGA-THCA dataset. Most patients were from North America. Therefore, caution should be exercised in extrapolating the conclusions from findings to other populations. Second, the targeting relationship of ceRNA is based on bioinformatics speculation and requires further experimental validation. Finally, the biological function of certain DE-circRNAs, DE-miRNAs, and DE-mRNAs requires investigation with experiments in PTC.

CONCLUSIONS
This study revealed a circRNA-associated ceRNA network in PTC. Based on survival analysis, a ceRNA subnetwork associated with PFI in PTC was identified. The molecules within the PFI-related subnetwork may represent a promising target for the treatment of patients with recurrent PTC. With the hub DE-mRNAs identified within the subnetwork, a prognostic six-DE-mRNA signature was established. The nomogram incorporating the DE-mRNA signature and clinical parameters is robust in predicting the PFI of post-surgical PTC.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: the datasets analyzed during the current study are available in the Gene Expression Omnibus (https://www.ncbi. nlm.nih.gov/geo/) and the Cancer Genome Atlas (https://portal. gdc.cancer.gov/).

AUTHOR CONTRIBUTIONS
ZL, XL, and XX: conception and design. MW, SL, JH, and RL: development of methodology. MW and HY: analysis and interpretation of data. MW: writing of the manuscript. XX, SL, JH, and XL: review of the manuscript. ZL: study supervision. All authors contributed to the article and approved the submitted version.