Identification and Validation of a Prognostic Signature Based on Methylation Profiles and Methylation-Driven Gene DAB2 as a Prognostic Biomarker in Differentiated Thyroid Carcinoma

Recurrence is the major death cause of differentiated thyroid carcinoma (DTC), and a better understanding of recurrence risk at early stage may lead to make the optimal medical decision to improve patients' prognosis. The 2015 American Thyroid Association (ATA) risk stratification system primary based on clinic-pathologic features is the most commonly used to describe the initial risk of persistent/recurrent disease. Besides, multiple prognostics models based on multigenes expression profiles have been developed to predict the recurrence risk of DTC patients. Recent evidences indicated that aberrant DNA methylation is involved in the initiation and progression of DTC and can be useful biomarkers for clinical diagnosis and prognosis prediction of DTC. Therefore, there is a need for integrating gene methylation feature to assess the recurrence risk of DTC. Gene methylation profile from The Cancer Genome Atlas (TCGA) was used to construct a recurrence risk model of DTC by successively performed univariate Cox regression, LASSO regression, and multivariate Cox regression. Two Gene Expression Omnibus (GEO) methylation cohorts of DTC were utilized to validate the predictive value of the methylation profiles model as external cohort by receiver operating characteristic (ROC) curve and survival analysis. Besides, CCK-8, colony-formation assay, transwell, and scratch-wound assay were used to investigate the biological significance of critical gene in the model. In our study, we constructed and validated a prognostic signature based on methylation profiles of SPTA1, APCS, and DAB2 and constructed a nomogram based on the methylation-related model, age, and AJCC_T stage that could provide evidence for the long-term treatment and management of DTC patients. Besides, in vitro experiments showed that DAB2 inhibited proliferation, colony-formation, and migration of BCPAP cells and the gene set enrichment analysis and immune infiltration analysis showed that DAB2 may promote antitumor immunity in DTC. In conclusion, promoter hypermethylation and loss expression of DAB2 in DTC may be a biomarker of unfavorable prognosis and poor response to immune therapy.


Introduction
Thyroid cancer has been the fifth most common cancer expected to be diagnosed in female and is the most common endocrine cancer according to the global cancer statistics 2020 [1]. Differentiated thyroid carcinoma (DTC) is the most common thyroid cancer, accounting for more than 90% of cases [1]. Although the long-term prognosis is favorable in the majority of patients with DTC, especially in the most popular subtype of papillary thyroid cancer (PTC), up to 30% patients experience recurrent disease after initial therapy [2,3]. Accurate assessment of the recurrence risk would be conducive to the clinical decision making, so as to reduce the recurrence rate and improve the quality of life for DTC patients [4]. The 2015 American Thyroid Association (ATA) risk stratification system, mainly based on the pathological characteristics such as tumor size, lymph node metastasis, and genetic feature such as BRAFV600E mutation, is the most commonly used to describe the initial risk of persistent/recurrent disease [5]. The response to therapy further integrates the real-time biochemical and imaging results, which is helpful for the dynamic monitoring of recurrence. Of note, over the recent decades, the transcriptional profiles and clinical prognosis information of tumor samples integrated by The Cancer Genome Atlas (TCGA) are beneficial for researchers to construct gene signatures for diagnosis and predicting the prognosis of patients with cancer [6][7][8][9][10][11]. However, these currently used gene signatures that solely dependent on data from TCGA cohort are not integrated to the risk stratification system and still limited to predict the recurrence risk early and accurately of patients with PTC.
DNA methylation is involved in transcriptional regulation and genome stability [12], and aberrant DNA methylation can lead to the development of many types of cancer [13,14], also including PTC [15]. Moreover, DNA methylation signatures can be used as biomarkers for clinical diagnosis and prognosis prediction of cancer [16]. For instance, Langdon et al. identified a panel of DNA methylation biomarkers for early diagnosis of renal cell carcinoma. Besides, DNA methylation of SPEG located at Chr2:220325443-220326041 could potentially regulate oropharyngeal cancer survival [17]. Shen et al. reported that hypermethylation of SSTR2 promoter could be used to predict prognosis of laryngeal squamous cell carcinoma in males [18]. The methylation profiles of GSTP1, HIC1, and RASSF1A were associated with the recurrence of bladder cancer [19]. Accumulating studies have reported that aberrant DNA methylation can be useful biomarkers for clinical diagnosis of DTC. For example, Wei et al. reported that hypermethylation of both PTEN and DAPK is capable of discriminating malignant thyroid tissues from benign lesions [20]. Besides, a panel of 3 DNA methylation biomarkers discriminated thyroid cancers from benign thyroid lesions was identified by Park et al. [21]. Several studies have reported aberrant DNA methylation of tumor suppressor genes was associated with tumor aggressiveness in PTC [22,23]. However, a panel of gene methylation profiles for predicting the recurrence of DTC is still lacking. Hence, developing and validating a signature based on methylation profiles is of important clinical significance.
In our study, we established and validated a panel of methylation profiles of SPTA1, APCS, and DAB2 for predicting the recurrence risk of DTC patients with data from TCGA. The significant effectiveness of our methylationrelated panel in predicting recurrence of DTC patients was validated by GSE51090 and GSE97466 cohorts. According to previous literatures [24], promoter hypermethylation of DAB2 associated with the loss expression of DAB2 in many types of human cancers, such as nasopharyngeal carcinoma [25], breast cancer [26], and esophageal squamous cell carcinomas [27].
Our study first reported that DAB2 was hypermethylation and low expressed in DTC samples and associated with proliferation, survival, and migration of DTC cells.

Materials and Methods
2.1. Data Acquisition. The methylation profiles, transcriptional profiles, and clinical information of TCGA-THCA were downloaded from TCGA website (https://portal.gdc .cancer.gov/). The data of 58 normal thyroid tissues and 497 PTC patients with survival information were extracted for our study (Table S1). Two external methylation cohorts of DTC were downloaded from Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/). The survival information of the 2 GEO cohorts (GSE51090, GSE 97466) was extracted from previous articles (Tables S2  and S3) [28,29]. For methylation profiles of genes, mean values were taken for multiple probes with an identical gene symbol. DTC patients with local recurrence or distant metastasis or biochemical evidence of disease were considered as recurrence DTC patients. Besides, GSE51090 cohorts and GSE97466 cohorts were integrated as an independent GEO cohort for validating the results. ComBat function was utilized to remove the batch effect between GSE51090 cohorts and GSE97466 cohorts ( Figure S1).

Constructed and Validated a Model Based on
Methylation Profiles of 3 DMGs. The TCGA cohort was utilized to construct a DMGs methylation profiles-related recurrence risk model of DTC patients. Univariate Cox regression, lasso regression, and stepwise multivariate Cox regression were successively performed to identify DMGs closely related to the recurrence of PTC with "survminer," "survival," and "glmnet" package [33]. To validate the prognostic prediction value of the model, receiver operating characteristic (ROC) curve analysis and survival analysis were performed by R software with "survivalROC" package  3 Disease Markers [34], and data form both TCGA and GEO cohorts. The risk-Score_me of 497 PTC patients were calculated as follows: X-Title software was utilized to find the best cutoff risk-Score_me which discriminate PTC patients with High Score and Low Score.

Tumor
Infiltrating Cells Analysis. ESTIMATE [35] and xCell [36] (https://xcell.ucsf.edu/) analysis were performed by R software according to the instruction in the website.

Constructed and Validated a Nomogram Based on
Methylation Profiles and Clinical Parameters. RiskScore_ me and clinical parameters (AJCC_T, AJCC_N, AJCC_M, stage, age, histological_type, and sex) of 497 PTC patients were enrolled for further analysis. Univariate and stepwise multivariate Cox regression analysis were performed to construct a recurrence risk model combined riskScore_me, age, and AJCC_T. To validate the prognostic prediction value of the model, ROC curve analysis, predict function, calibration curve, and survival analysis were performed by R software. Finally, a nomogram based on methylation profiles, age, and AJCC_T was established for predicting recurrence probability of DTC patients by R software with "rms" [39] and "foreign" package. Decision curve analysis (DCA) was performed to compare the predicting value of our model with other 2 previous models by R software with "ggDCA" package [40].
2.7. GSEA Analysis. GSEA (Gene Set Enrichment Analysis) [41,42] was performed by R software with "clusterProfiler" package. The KEGG pathways with P < 0:05 and FDR < 0:05 were considered as pathways that are significantly correlated with expression of DAB2 in PTC samples.

Human Thyroid Carcinoma
Specimens. 34 thyroid carcinoma specimens along with paired adjacent normal thyroid specimens collected from Affiliated Hospital of Jining Medical University was approved by the Ethics Committee of Affiliated Hospital of Jining Medical University. The approval number was 2021-08-C015. All participants provided written informed consent.
2.11. Plasmids and Lentivirus Production. DAB2 coding sequence was cloned into the viral skeleton plasmid PCDH-3 flag-tagged vector for stable expression in BCPAP cells.
2.12. Western Blotting. Cell proteins were extracted by denatured buffer and then quantified by pierce BCA protein assay (Thermo Scientific). The protein lysate was separated on  and HRP-conjugated secondary antibody (Sigma), followed by being exposed with enhanced chemiluminescence (ECL).
2.13. Cell Growth Assay. Lentivirus infected stable BCPAP cells were seeded into 96-well plates and cultured in 10% FBS DMEM (2000 cells per well, five parallel wells). Then, the cells were collected at different points in time, and cell number in each well was counted by CCK-8 reagent. The absorbance at 450 nm of each well was measured to calculate the cell proliferation.
2.14. Colony-Forming Assay. Lentivirus infected stable BCPAP cells were seeded into 6-well plate at a density of 200 cells per well and cultured in DMEM with 10% FBS for 3 weeks. Then, colonies per well were stained by 0.1% crystal violet and counted by ImageJ software. 2.16. Scratch-Wound Assay. BAPCP cells were seeded in 6well plates and then incubated for 24 h to reach approximately 80% confluence. The cell monolayer was scratched using a sterile 100 μL pipette tip. Then, the cells were treated with serum-free medium. The wounds were photographed at 0 and 36 h, and migration distance of the cell was calculated by ImageJ (migratory ratio: 0-36 h width of wound/ 0 h width of wound).
2.17. Statistical Analysis. Two-tailed t test was utilized to analyze the difference between two groups. Log-rank test was utilized for survival analysis.

Constructed and Validated a Recurrence Risk Model of DTC Based on Methylation
Profiles of 3 DMGs. The methylation profiles of 10362 genes in the intersection of TCGA and 2 GEO cohorts were extracted for our research (Figure 1(a)). 61 DMGs were identified between DTC and normal thyroid samples (Figure 1(b)). By performing univariate Cox regression with data from TCGA cohort, DMGs   9 Disease Markers with P < 0:01 were selected for further analysis (Figure 1(c)). Then by successively performing lasso regression (Figures 1(d) and 1(e)) and stepwise multivariate regression analysis, methylation profiles of SPTA1, APCS, and DAB2 were utilized to construct a recurrence risk mode of PTC patients (Figure 1(f)). In the TCGA cohort, the area under the ROC curve (AUC) of 1-year, 3-year, and 5-year recurrence free survival were 0.758, 0.618, and 0.691, respectively (Figure 2(a)). The riskScore distribution, survival status of PTC patients, and the methylation profiles heat map of 3 genes between Low Score and High Score PTC samples are shown in Figures 2(b)-2(d). Kaplan-Meier result indicated that PTC patients in Low Score group had longer recurrence free survival (RFS) than those in High Score group (Figure 2(e)). As external validated cohorts, the AUC of 1year, 3-year, and 5-year recurrence free survival were 0.575, 0.659, and 0.703 in the GES51090 cohort, respectively (Figure 3(a)) and PTC patients in High Score group had worse RFS than those in Low Score group (Figure 3(b)). Similarity, the AUC of 1-year and 3-year recurrence free survival were 0.596 and 0.714 (Figure 3(c)) and High Score PTC patients had worse RFS (Figure 3(d)) in the GES97466 cohort. Finally, the data of 2 GEO cohorts were integrated as a single external validated GEO cohort for  Disease Markers further analysis. The batch effects between GSE51090 and GSE97466 were removed by "combat" function of R "sva" package ( Figure S1). The result showed that the AUC of 1year, 3-year, and 5-year recurrence free survival were 0.656, 0.724, and 0.762 (Figure 3(e)) and High Score PTC patients had shorter RFS (Figure 3(f)) in the single GEO cohort. Obviously, the results of validation were consistent between internal and external cohorts.

Constructed a Nomogram of DTC Based on Methylation
Profiles of 3 Genes, Age, and AJCC_T. By performing univariate Cox (Figure 4(a)), 3 clinical factors (age, AJCC_T, and stage) (P < 0:05) were selected to experience stepwise multivariate Cox regression analysis with riskScore_me and a recurrence risk model of DTC combing age, AJCC_T, and risk_me was constructed (Figure 4(b)). The calibration curve for 3-year and 5-year recurrence free survival shown a consistency between predicted value of model and actual observed value (Figures 4(c) and 4(d)). The AUC of 1-year, 3-year, and 5-year RFS were remarkably increased to 0.838, 0.734, and 0.775, respectively (Figure 4(e)). The Kaplan-Meier analysis indicated that the DTC patients in Low Score group had significantly longer RFS than those in High Score group (Figure 4(f)). Finally, we constructed a nomogram based on methylation-related recurrence risk model, age, and AJCC_T for predicting RFS probability of DTC patients ( Figure 5).  (Figures 6(a)-6(d)). However, by performing DCA, we found that the model based on methylation profiles had outstanding performance in predicting 5-year RFS of patients with DTC. Of note, the nomogram which combined risk_ me, age, and AJCC_T would dramatically improve the long-term treatment and management of DTC patients ( Figure 6(e)).

DAB2 Was Abnormally Expressed Between DTC Samples
and Normal Thyroid Tissues. For better understanding the biological significance of genes in the model, GEPIA was utilized to investigate the differences in mRNA expression between DTC and normal thyroid tissues. Boxplot showed that the mRNA expression of DAB2 was downregulated in DTC samples compare to normal thyroid tissue (P < 0:05) (Figure 7(a)) while the mRNA expressions of SPTA1 and APCS (Figures 7(b) and 7(c)) were no significant difference between cancer and para-cancer tissues. Then, we found that DAB2 expression was lower in DTC tissues than that in normal thyroid tissues at the protein level by IHC staining for DAB2 in our specimens (Figure 7(d)) and found that DTC patients with high expression of DAB2 had a longer RFS compare to that with low expression DAB2 (Figure 7(e)).
Overall, DAB2 associated with pathological feature and prognosis of DTC patients.

DAB2 Overexpression Inhibited Proliferation, Colony-Formation and Migration of BCPAP Cells.
To explore whether does DAB2 impacts DTC cells, we constructed and validated a DAB2 stable expressed BCPAP cells in our study (Figure 8(a) and Figure S2). In the results from CCK8 assay, colony-formation assay demonstrated that overexpression of DAB2 inhibited BCPAP cells proliferation (Figure 8(b)) and colony-formation (Figure 8(c)). Besides,

12
Disease Markers transwell assay (Figure 8(d)) and scratch-wound assay (Figure 8(e)) showed that DAB2 inhibited the migratory ability of BCPAP cells. For better understanding the potential biological significance of DAB2 in the progress of DTC, all 497 PTC patients were assigned into High DAB2 and Low DAB2 group based on the median mRNA expression of DAB2 using data from TCGA. The differential expression genes (DEGs) with FDR < 0:05 and jlog 2 fold changej > = 1 between High DAB2 and Low DAB2 group were identified with "limma" package ( Figure S3A), and the result of KEGG analysis using these DEGs as input showed that DAB2 regulated the migration and proliferation of DTC may through regulating signaling pathways such as Cell adhesion molecules, NF-kappa B signaling pathway, ECMreceptor interaction, and PI3K-Akt signaling pathway ( Figure S3B).

DAB2 may Promote Antitumor Immunity in PTC.
According to the result of KEGG analysis, DAB2 may participate in the T cell receptor signaling pathway ( Figure S3B). To further investigate the association between DAB2 expression and immune microenvironment in DTC, the GSEA, ESTIMATE, and xCell analysis were performed. We found that some antitumor immunity related pathways were enriched in High DAB2 group, such as T cell receptor signaling pathway, natural killer cell mediated cytotoxicity, antigen processing and presentation, and B cell receptor signaling pathway. (Figure 9(a)). Both the xCell and ESTIMATE analysis showed that High DAB2 group had higher abundance of immune infiltration score than Low DAB2 group (Figures 9(b) and 9(c)). The cell type infiltration abundance analysis showed that High DAB2 group had higher infiltration abundance of CD4+memory/ naïve T cells, CD8+T cells, B cells, and DC cells than Low DAB2 group (Figure 9(d)). Moreover, we found that the expression of DAB2 was positively correlated with the expression of classical biomarker of immune cells in PTC samples, such as CD3G (for T cells), MS4A1 (for B cells), KLRD1 (for NK cells), and immune checkpoint PDL1 (Figures 9(e)-9(h)). Based on the above results, we inferred that PTC patients with high expression of DAB2 would be benefit from immune-based therapy.

Discussion
Although there have been some studies on the relationship between abnormal DNA methylation and prognosis of DTC patients [21]. No reliable model considering the DNA methylation-driven gene signature to predict recurrence and metastasis risk accurately was reported before. In our study, we constructed a recurrence risk model with methylation profiles of SPTA1, APCS, and DAB2 by stepwise multivariable Cox regression with data from TCGA. SPTA1, which have been linked to hereditary elliptocytosis and hereditary spherocytosis [45], were reported as a possible tumor driver gene in prostate cancer [46] and papillary thyroid carcinoma [47]. APCS is a gene that codes serum amyloid P-component, which is one of the main acute-phase      17 Disease Markers reactants and has reported to a biomarker for survival in nonsmall cell lung cancer after thoracic radiotherapy [48]. Infrequent promoter hypermethylation of DAB2 have been found to play a critical role in tumorigenesis and progression according to previous studies [24,26]. Of note, we developed and validated a nomogram combined the methylationrelated model, age, and AJCC_T for predicting the recurrence probability of DTC, which has not been published before to the best of our knowledge. DCA indicated that our model showed more net benefit for 5-year RFS than age, AJCC_T, He et al.'s, and Zhang et al.'s model, and when we combined our model with clinical parameters, DTC patients would gain the greatest benefit for 1-year, 3-year, and 5-year RFS.
Analysis of the 3 genes in the model at the transcriptional level showed that DAB2 is a methylation-driven gene in DTC. The mRNA expression of DAB2 was downregulated significantly in DTC compared to normal thyroid tissue. However, SPTA1 and APCS both in DTC and normal thyroid tissues were in a state of hypermethylation and low expressed. Therefore, the function of DAB2 was selected to be further explored in the DTC progression. Numerous literatures reported that DAB2 involved in the migration, invasion, and proliferation of tumors [49][50][51][52]. However, the roles of DAB2 in DTC have rarely been explored. In our study, we found that DAB2 expression was lower in DTC tissues than that in normal thyroid tissues at the protein level, which was validated by tissue microarray staining for DAB2 in specimens collected from Affiliated Hospital of Jining Medical University. Moreover, we found that overexpression of DAB2 could inhibit proliferation and migration of BAPCP cells in vitro experiments. Enrichment of pathways using data from TCGA suggested that DAB2 was relevant to some pathways that played various essential roles in the tumor growth and aggressive phenotypes such as NFkappa B signaling pathway and PI3K-Akt signaling pathway. Activation of NF-kappa B is reported in DTC and correlates with tumor growth, drug-induced apoptosis [53,54], and aggressiveness of DTC [55]. The activation of PI3K-Akt signaling pathway 395 also involves in the initiation and progression of thyroid cancer [56,57]. Interestingly, we firstly focused on the relationship between DAB2 expression and immune microenvironment in DTC. Pervious study reported that, loss of DAB2 expression induced immune tolerance via accumulation of TGF-β in breast cancer [58]. In our study, we found that expression of DAB2 was correlated with immunoscore in DTC samples and may be a biomarker for predicting the response to immune therapy in patients with DTC.
Although we successfully developed and validated a nomogram that may predict the recurrence of DTC patients, there were still some limitations in this study. Firstly, our external validating cohorts only included 115 DTC patients that were not completely sufficient, moreover, due to the insufficient clinical characteristics in the current study cohorts, there is a need to build a better prognostic nomogram from more centers with complete clinical information and sequencing data. Secondly, further experimental verification on molecular mechanism of DAB2 in DTC is required to consolidate our results.

Conclusions
We successfully constructed and validated a nomogram based on methylation profiles of 3 genes, age, and AJCC_T for predicting recurrence probability of DTC patients. Besides, we found that promoter hypermethylation and loss expression of DAB2 may be a biomarker of unfavorable prognosis and poor response to immune therapy in DTC.

Data Availability
The main results of this study are based on public data from TCGA (https://portal.gdc.cancer.gov/) and GEO (https:// www.ncbi.nlm.nih.gov/geo/).The clinical information and IHC data for human specimens are shown in Table S4.

Ethical Approval
All experimental protocols and human specimens collected from Affiliated Hospital of Jining Medical University were approved by the Ethics Committee of Affiliated Hospital of Jining Medical University [number was 2021-08-C015]. All participants provided written informed consent.

Conflicts of Interest
There are no competing interests declared for all authors.