Exploring Somatic Alteration Associating With Aggressive Behaviors of Papillary Thyroid Carcinomas by Targeted Sequencing

Wisely differentiating high-risk papillary thyroid carcinoma (PTC) patients from low-risk PTC patients preoperatively is necessary when comes to making a personalized treatment plan. It is not easy to stratify the risk of patients according to sonography or lab results before surgery. This study aims to seek out potential mutation gene markers that may be helpful in stratifying the risk of PTC. A custom panel of 439 PTC relevant and classic tumor metabolic pathway relevant genes was designed. Targeted capture sequencing was performed on 35 pairs of samples from 35 PTC tumors and 35 para-tumor thyroid tissues obtained during surgery. Variant calling and detection of cancer gene mutations were identified by bio-information analysis. Ingenuity Pathway Analysis (IPA) was performed to do functional enrichment analysis of high-frequency mutant genes. Immunohistochemistry (IHC) was performed on 6 PTC patients to explore the expression of protein associated with interested genes. Event-free survival (EFS) was calculated to determine which genes might affect the prognosis of patients. We have identified 32 high-frequency mutant genes in PTC including BRAF. RBL2 was found to be significantly correlated to event-free survival, FOXO1, MUC6, PCDHB9, NOTCH1, FIZ1, and RTN1 were significantly associated with EFS, while BRAF mutant was not correlated to any of the prognosis indicators. Our findings in this study might open more choices when designing thyroid gene panels used in FNA samples to diagnose PTC and predict the potentially aggressive behavior of PTC.


INTRODUCTION
Thyroid carcinoma is the most common malignant tumor in the endocrine system which accounts for 96% malignant tumors of the head and neck. Most occurrences of thyroid carcinoma are papillary thyroid carcinoma (PTC) (1)(2)(3)(4). Although overall prognosis of PTC is good and the mortality is stable at the rate of approximately 0.5 deaths per 100,000 (5,6), clinical evidence has shown the performances of PTC patients are quite different (3,7). On average, 50% of the patients are presented with lymph node metastases (LNM) at diagnosis, and in some cases, patients with micro papillary thyroid carcinoma could have macronodular lung metastasis (8,9).
Considering the different nature of PTC and potential complication of thyroidectomy, doctors have to make a personalized treatment plan by wisely differentiating high-risk PTC patients from low-risk ones (10). Surgery is the most effective way to treat patients with PTC (11). Meanwhile, some other techniques are also applied gradually in PTC treatment, such as ultrasound-guided laser, radiofrequency, microwave ablation and percutaneous ethanol injection therapy (12)(13)(14). A strategy of active surveillance (AS) instead of immediate surgery is also proposed for confirmed low-risk PTC patients (15).
But how should high-risk PTC patients be wisely differentiated from low-risk PTC patients preoperatively? Over the years, multiple staging systems have been developed to predict the risk of mortality in patients with PTC. Each of the systems uses particular combinations of some clinical factors including age at diagnosis, size of the primary tumor, specific tumor histology, and extrathyroidal spread of the tumor (direct extension of the tumor outside the thyroid gland, loco-regional metastases, and/or distant metastases), which tend to stratify patients into one of several categories with different risks of death from thyroid cancer (11,16). Pre-operating TSH lever, TMN stage when diagnosed, and histology subtypes are also factors associated with PTC progression (5,7). However, preoperative evaluation still has difficulties with demonstrating all of the clinical information because some of factors can only be confirmed by post-operation pathological results. Fine needle aspiration (FNA) is an efficient method for evaluating thyroid nodules that has high sensitivity and specificity, and gene testing based on FNA samples could be a potential tool to stratify the risk of PTC. The key problem for this possibility is which genes and types of the alteration should be chosen to identify for screening the high-risk patients with PTC. Figure 1 demonstrates the workflow of FNA of the thyroid nodule and how gene testing participates in the management of the thyroid nodule.
Presently, next generation sequencing (NGS) technology provides the simultaneous analysis of hundreds of genes of interest, using targeted sequencing panels and NGS-based molecular tests performed on thyroid cancer to seek more valuable information to differentiate III/IV nodules from malignant ones (17)(18)(19). Currently, there is less evidence showing specific gene mutations reported for the recurrence or aggressive behaviors of PTC. Thus, in this study, a custom panel was designed and targeted sequencing would be performed to seek potential mutation gene markers helpful for stratifying the risk of PTC.

General Information and Groups of Patients
The study included 41 PTC patients operated on in hospital between June 2016 and January 2017. 70 tissue samples including tumor tissue (n=35) and para-carcinoma tissue (n=35) samples were obtained during surgery and preserved in liquid nitrogen immediately and 6 paraffin-embedded samples obtained from other patients were sliced. All patients provided written informed consent for sample collection and analysis. The study protocol was approved by Shanghai 6 th people's hospital institutional review board and adhered to the Declaration of Helsinki.
(i) Candidate samples for targeted NGS The cohort of 35 patients in the study group was categorized into two groups: recurrence group (n=8) and no-recurrence group (n=27). The standard of recurrence was as below: new lymph nodes appearing confirmed as metastasized by cytological exam; New mass appearing around thyroid bed confirmed as metastasized by cytological exam; New lymph nodes with obvious malignant characteristics such as liquefaction or calcification with increased TG level; Distant metastasis; Obvious increase in TG without positive imaging findings but undergoing significant reduction after iodine therapy. Meeting any of these rules would be considered as recurrence. If the patient did not meet any of the rules during the following up period, it would be considered as no-recurrence.
Besides being grouped as recurrence or not, the patients were also regrouped into aggressive PTC group (n=18) and mild PTC group (n=17). The aggressive PTC group contained patients with extra-thyroidal invasion or severe lymph nodes involvement (≥5) at surgery or with recurrence while the mild group did not meet any of the aggressive behaviors mentioned above.
(ii) Candidates for Immunohistochemistry (IHC) Samples of tumor tissue and normal thyroid tissue from 6 PTC patients in the assessment group which were prepared for validating the NGS findings were randomly selected among 318 patients who underwent thyroid surgery in our hospital. The information of the 6 patients was listed in Table S1. One patient died 3 years after their first surgery because of hemoptysis due to PTC recurrence. One patient was confirmed to have multiple lymph nodes metastasis months after surgery. One patient was confirmed to have bone metastasis during surgery. These three patients belonged to the high-risk assessment group, and the other three patients belonged to low-risk assessment group, which were considered relatively mild because of unilateral nodules, no metastatic lymph nodes, and no recurrence in 5 years of follow-up.

Following Up
All of the patients were followed up and the last dates of followup were recorded on August 21 st 2020. The events of death, recurrence, or metastasis were recorded by follow-up secretary.

Design of SureSelect Custom Panel of PTC Relevant Genes
We performed targeted capture sequencing of prior-collected PTC genes in 70 tissue samples including 35 tumor samples and 35 para-carcinoma tissues, as described previously.
We designed a SureSelect custom panel of PTC relevant genes (Agilent Technologies, Santa Clara, CA) ( Table S2). The gene list in the panel includes genes recurrently mutated in PTC; genes mutated in stage IV PTC analyzed and extracted from TCGA and COSMIC database; genes closely related to the relevant biological pathways including cell cycle, cell growth and proliferation, survival/cell death regulation, TP53 signaling, Notch signaling, DNA Damage response, telomere maintenance, RTK signaling family, PI3K-AKT-mTOR signaling, RAS-RAF-MEK-ERK/JNK, regulation of ribosomal protein synthesis and cell growth angiogenesis folate transport invasion and metastasis; and genes mentioned in relative literature and commercial panel. There is 3.342 Mbp in each panel including 439 genes, 8541 regions and 109254 probes. Extracted genomic DNA was fragmented and bait-captured according to manufacturer protocols. Captured DNA libraries were then sequenced using a HiSeq 2000 sequencer (Illumina, San Diego, CA). Raw image files were processed by Illumina base-calling Software 1.7 for base-calling with default parameters and the sequences of each individual were generated as 90/100bp pair-end reads.

Variant Calling and Detection of Genes With High-Frequency Somatic Mutation
Raw FASTQ data were generated from the Illumina pipeline, the adapter sequence in the raw data and low-quality reads were removed in the QC step, and the remaining "clean reads" were used for downstream bioinformatics analysis. Standard bioinformatics analysis was firstly performed on the clean data by Burrows-Wheeler Aligner (BWA) to do the sequence alignment. The human genome build37 (hg19) was used as the reference genome for this project which is available at http:// hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/. Picard (http://picard.sourceforge.net/) was used to fix mate information of the alignment, add read group information and remove duplicate reads caused by polymerase chain reaction (PCR). After these processes, the obtained BAM files were used to do the variant calling. Single Nucleotide Polymorphisms (SNPs) was obtained by GATK (http://www.broadinstitute.org/gsa/wiki), and annotated on the basis of known protein-coding and non-proteincoding genes taken from the NCBI RNA reference sequences collection (RefSeq). To remove false-positive mutations, these SNPs were validation and compared in databases like dbSNP (The Single Nucleotide Polymorphism Database) and 1000 Genomes (The 1000 Genomes Project), and their functionality and conservation annotation were carried out according to SIFT, PolyPhen-2, PhyloP, etc. Meantime, the InDels were also called and annotated as the similar pipeline. Somatic Single Nucleotide Variants (SNVs) and Small Insertion/Deletion (InDels) were detected by Varscan (20) and GATK, which identified the tumor specific somatic substitutions by comparing PTC tumor and normal para-carcinoma in a pair. Then they were given functional annotation and classification by ANNOVAR in the above similar filtering pipeline. QualByDepth (QD), FisherStrand (FS), RMSMappingQuality (MQ), MappingQualityRankSum (MQRankSum) and ReadPosRankSum were applied to filter out variants with low quality. Only variants which would cause an obvious protein-coding change were kept, where specifically, the variants with an ANNOVAR annotation of missense_variant, splice_acceptor_variant, splice_donor_variant, splice_ region_variant, start_lost,stop_gained, stop_lost, initiator_ codon_variant, were counted. Mutation plots were generated using R package complexHeatmap (21).

Functional Enrichment Analysis of Genes With High-Frequency Somatic Mutation by Ingenuity Pathway Analysis
We captured the high-frequency mutant genes of PTC through counting the mutation frequency of the genes with proteincoding change SNVs in all patients. These genes are generally thought to be related to the pathogenicity of thyroid cancer. To further uncover potential disease-associated biological functions regulated by these candidate driver genes, we analyzed these genes by IPA database with diverse functional annotation and enrichment analysis (22).

Mutant Genes With High-Frequency Somatic Mutation Associated to Aggressive Behaviors of PTC
Gene mutation frequency was presented as frequencies or percentages, Fisher's exact test was carried out between "mild PTC" and "aggressive PTC" to search the mutant genes associated with high-risk progression of PTC patients, where statistical significance was defined as one-sided values of p<0.05.
To guarantee clinically detectable/observable protein expression of these mutant genes in thyroid tissue, these candidates were all screened in THE HUMAN PROTEIN ATLAS database (23), which can help provide efficient targets for Immunohistochemistry (IHC) validation of potential protein expression/function change involved in PTC (24).

Immunohistochemistry
Samples of 6 PTC patients were taken to detect the potential protein expression/function alterations of final gene candidates we found, where the 6 patients of IHC assessment group were independent from 35 patients of study group mentioned before ( Table S1). The catalogs of the antibodies for the 5 proteins are (BS71182 from Bioworld, 101107-T10 from Yiqiao Shenzhou, Ab189647 from Abcam and CY5319 from Abcam).
The samples were fixed in 4% buffered paraformaldehyde, decalcified in 50 mM EDTA, embedded in paraffin, and sectioned at 5-mm thickness. Samples were deparaffinized and rehydrated. Twelve high-power fields were randomly selected and photographed in each sample, and the percentage of positive area was used to evaluate coding proteins of target genes. Immunochemistry was used to evaluate protein expressions of candidate mutant genes. The average optical density also called IOD (integrated optical density/area) of positive reactions was analyzed with Image Pro Plus 6.0 software.

Exploration of the Potential Prognostic Mutant Gene
For each PTC patient, the event-free survival (EFS) was defined as the time between surgical treatment and recurrence. Survival curves were plotted with the Kaplan-Meier method, and differences between mutant and non-mutant groups were compared by the log-rank test. Statistical significance was defined as two-sided values of P < 0.05.

Statistics
Clinical data was presented either as frequencies and percentages, or as means and standard deviations, or medians and interquartile range. Categorical variables were compared using either Pearson chi-square test or Fisher's exact test. Continuous variables were compared using either the independent t test or one-way analysis of variance. Statistical significance was defined as two-sided values of P < 0.05. Of note, as limited number of samples and targeted sequencing were used in this study, the (high) mutation frequency instead of multitest correction was applied to further control the potential false-positive.
Statistical analyses were conducted with SPSS version 23.0 (SPSS Inc). Figure 2 demonstrates the method and flow of the study.

Code and Data
The raw data generated for this study can be found in The National Omics Data Encyclopedia (NODE), project ID OEP002494. The related analysis and visualization scripts and supported data can be accessed at https://github.com/ graceyu1985/TargetedPTC.

Patient Characteristics
Our study group cohort underwent NGS consisting of 35 patients diagnosed with classic PTC. The median age was 43.9 years, ranging from 15 to 80 years, and these patients were predominantly female (60%, 21/35). The median PTC tumor size was 1.8 cm, ranging from 0.3 to 5.0 cm. 62.9%(22/35) of the patients underwent total thyroidectomy, 54.3%(19/35) of the patients underwent unilateral or bilateral central part lymph node dissection while 45.7% (16/35) of the patients got central & lateral cervical lymph node dissection. After surgery, 28.6% (10/ 35) of the patients got iodine RAI therapy of which 5 patients did this for prophylactic treatment. Patients were staged using the eighth edition of the American Joint Cancer Committee/Union for International Cancer Control (AJCC/UICC) staging system. Also, patients' recurrence risk was stratified by ATA guidelines (11). One patient died of hemoptysis in the second year of following up. Another 7 patients were confirmed to have recurrence during the 4-year following up. 27 patients did not show any ultrasonographic evidence or lab results indicating recurrence. The patients were divided into two groups according to whether they recurred, no less than 5 metastatic lymph nodes confirmed after the first surgery, the preoperative ATA risk classification and RAI therapy was related to recurrence(p=0.031, 0.012 and 0.048). Tumor size, age, gender, extra thyroidal invasion, BRAFv600E mutation, multifocality, and TMN Stage did not show significant relation to recurrence ( Table 1).

Mutation Profiling of Patients With PTC
We analyzed 439 genes including 8541 regions and 109254 probes using the panel on matched tissue samples from 35 PTC patients. The statistical results are shown in Figure 3a and Table S3. The average sequencing depth of the target region is about 367.77 X. on average, 99.82% of the target region is covered by at least one read, and 99.12% of the target region is covered by at least 10 reads. Mutations were identified in 35 (100%) PTCs, of which BRAF (62.86%) and NCOR2 (51.43%)  mutations were exceeded among 50% patients. Multiple highfrequency mutant genes (with mutations in >20% patients) were found in PTCs by the NGS assay and further analysis is required to reveal PTC-associated biological functions regulated by these mutant genes. We analyzed 32 high-frequency mutant genes (mutant in >20% patients) by IPA (Ingenuity Pathway Analysis) database with diverse functional annotation and enrichment analysis. The results from IPA analysis (Table S4) also supplied many other functional evidences for our selected genes. For example, in the pathway enrichment analysis results (Figure 3b), "Notch Signaling" pathway is significantly enriched and mutant genes can participate in many novel pathways such as: "Regulation of the Epithelial-Mesenchymal Transition Pathway", and "Tumor Microenvironment Pathway".
Of note, we have also summarized the mutation frequency of another PTC driver gene RAS as shown in Figure S1. Actually, RAS tends to have a remarkably high mutation frequency, and especially have more mutations in "aggressive PTC" than those in "mild PTC". However, the dominant mutation type annotated of RAS was "intron variant" in our data. Thus, other PTC driver genes listed in Figure 3c were considered in the following analysis and discussion.

Genetic Comparison of "Mild PTC" and "Aggressive PTC"
We further examined the gene mutations in the subgroups "mild PTC" (i.e.17 patients from study group) and "aggressive PTC" (i.e. 18 patients from the study group). We carried out fisher's exact test analysis between "mild PTC" and "aggressive PTC" to search the differentially mutated genes associated with PTC aggressive behaviors. The genes with one-sided values of P < 0.05 were final candidates. We found different altered genes through different group methods, and the genes are listed in Table 2. After being screened in THE HUMAN PROTEIN ATLAS database, we finally confirmed several candidates including "MDC1, NCOR2, RBL2" as well as BRAF for the following IHC experiment to further confirm their protein expression state and change during PTC occurrence and progression. Of note, we have also tried to build a logistic regression (LR) model based on these mutant genes for predicting risk classification. The LR model was built on 70% of samples, and evaluated to predict the remaining 30% of samples for belonging to "mild PTC" or "aggressive PTC", which achieved AUC=0.64 and ci= [0.33,0.95] ( Figure S2).

Protein Expression of Mutant Gene Candidates in PTC
Immunochemistry was used to evaluate protein expression of MDC1, NCOR2, RBL2, and BRAF, which can provide additional complementary clinical assessment of the functional protein loss of these coding genes (25). According to the IHC results, all of these candidates have higher protein expressions in tumors no matter their categorization as mild cases or aggressive cases than those in normal thyroid tissue (Figure 3c and Table S5). Particularly, the average expression of MDC1 protein was higher in aggressive cases than those in mild cases (0.434 ± 0.057 vs 0.177 ± 0.021). And the average expression of RBL2 protein was lower in aggressive cases than those in mild cases (0.165 ± 0.016 vs 0.371 ± 0.054).

Prognostic Mutant Gene for PTC
In addition, the prognostic potential of those four mutant gene candidates were also evaluated by Kaplan-Meier method. RBL2 was found to be remarkably related to event-free survival (P=0.016). Considering many potential functional contributions except for protein expression would also be caused by mutation, other differentially mutated genes associated with PTC disease progression were re-screened again, and it was actually found that FOXO1, MUC6, PCDHB9, NOTCH1, FIZ1, RTN1 could be significantly associated with EFS ( Figure 3d).

DISCUSSION AND CONCLUSION
Several studies have applied NGS on PTCs (3,9,26,27). The Cancer Genome Atlas (TCGA) Network described the genomic characterization of 496 PTCs, and generated data using whole genome sequencing. The TCGA project identified the TERT promoter mutation, which accounts for approximately 1% of PTCs, but shows association with a high risk of recurrence (19,27). Nikiforova et al. applied the ThyroSeq series NGS panel targeting 12 genes or 14 genes to differentiate FNA III/IV nodules (18,28). Picarsic et al. analyzed 17 pediatric PTCs by using ThyroSeq v2 panel and found ETV6-NTRK3 fusion was associated with aggressive histologic features. In this study, we focused on core genes related to aggressive behavior of PTC. We not only use recurrence as an indicator to select genes, but also considered severe lymph nodes involvement and extrathyroidal invasion as high risk indicators based on the evidence in the 2015 American Thyroid Association (ATA) management guidelines for adult patients with thyroid nodules and differentiated thyroid cancer (11). In this guideline, it was said that neck lymph node (LN) metastases are found in up to 70% of cases of PTC, and initial LN metastases as well as extra-thyroid infiltration have also been identified as an independent risk factor for distant metastases (29,30). So, we considered all of the high-risk factors and found that RBL2 was significantly correlated to event-free survival, FOXO1, MUC6, PCDHB9, NOTCH1, FIZ1, RTN1 could be significantly associated with EFS, while BRAF mutant was not correlated to any of the prognosis indicators. BRAFV600E mutation is considered as the most common mutation in PTC, but its role in predicting prognosis is unclear (31,32). Recent data suggest that specific molecular profiles, such as the coexistence of BRAF with other oncogenic mutations (such as PIK3CA, AKT1), TERT promoter, or TP53 mutations, may serve as more specific markers of less favorable outcomes of PTC (11). In our study, BRAFV600E showed a high mutation rate (62.86%) in the 35 patients, but it did not show any differentiating capability in picking cases of "aggressive PTC". Other oncogenic mutations including TERT promoter and TP53 mutations also did not show differentiating capability either. Another interesting finding is that RAS tend to have a remarkably high mutation frequency, and especially have more mutations in "aggressive PTC" than those in "mild PTC". However, the dominant mutation type annotated of RAS was "intron variant" in our data. PTCs were separated into two groups of distinct downstream signaling pathways: the BRAFV600E-like cohort and the RAS-like cohort and most of the RAS-like PTCs were follicular variant PTCs (FVPTCs) (19,27). In our cohort, there are all classic PTCs, it might help explain why meaningful RAS mutation was rare in our results. So, we did not do further bio-informative analysis in the RAS mutation pathway.
In this study, we have detected some mutant gene related protein candidates to explore the potential influence of these gene signatures, and successfully found that the expression of MDC1 and RBL2 proteins were different in aggressive cases than mild cases.
Rbl2(RB transcriptional corepressor like 2) belongs to the retinoblastoma (Rb) family which are considered tumor suppressors (33). RB is closely related to two genes in mice and humans (p107 and 130, which are not commonly mutated in tumors) (34). The retinoblastoma protein (RB) has been implicated in many cellular processes, such as regulation of the cell cycle, DNA-damage responses, DNA repair, DNA replication, protection against apoptosis, and differentiation. They could contribute to the function of RB as a tumor suppressor. The Rb can inhibit cell cycle progression through disabling the E2F family of cell-cycle-promoting transcription factors (35). In our study, we found RBL2 mutated in tumor samples of the group confirmed as recurrence and its protein also called Rbl2 expressed in normal thyroid samples and PTC samples. However, the expression level was lower in patients with recurrence than in tumor samples belonging to patients without recurrence. There were few articles discussing the role of RB in the PTC driving process, and our study was the first to show a potential relation between RB and the prognostic results of thyroid cancer, which deserves further study.
Mediator of DNA damage checkpoint 1 (MDC1) is a scaffold protein involved in the early steps of DDR functioning as a tumor suppressor (36,37). The DNA damage response (DDR) is necessary to maintain genome integrity and prevent the accumulation of oncogenic mutations (38,39). Consequently, proteins involved in the DDR often serve as tumor suppressors, carrying out the crucial task of keeping DNA fidelity intact (40). Studies have found a correlation between MDC1 loss and tumor aggression (41). In breast cancer and pancreatic cancer, low MDC1 expression correlates with more invasive, later stage tumors, indicating that MDC1 expression may serve to limit the progression of disease (39,(42)(43)(44). By contrast, it was found that MDC1 is upregulated in most cervical cancer samples (45), laryngeal squamous cell carcinoma samples (SCC) (46), and nasopharyngeal carcinoma  (39,48). In this study, we observed mutation of MDC1 in the aggressive group, which indicates that the mutation might relate to more aggressive behavior. There are several studies that have reported a correlation between HBV and thyroid cancer (49,50). And we might also boldly make a hypothesis that virus and viral carcinogenesis could play an important role in the drive and progression of PTC mediated by MDC1 with inducing epithelialmesenchymal-transition. The study of MDC1 and thyroid cancer is few, and our results just unveiled the tip of the iceberg of its relations. This association between MDC1 expression and its potentially oncogenic and tumor suppressive (dual) roles indicate an opportunity for further research such as a 3D tumor microenvironment on thyroid tumor development to understand its complex role in thyroid cancer (39,51). Of course, there are still several limits in the study. The numbers of samples were only 35 pairs which was relatively small to capture more mutant gene candidates. It is a single institute investigation using one classic type of PTC rather than all types of thyroid cancers so the sample size was relative small. The resulting data would still be insufficient for making decisions on either patient diagnosis or treatment in clinical practice. To overcome such limitations, an expansion of the study cohort with significantly increased patient sample size in our future work should help investigate and validate more and more new clinical findings on the basis of next-generation sequencing and immunohistochemistry technologies, such as building a more accurate and robust LR model for predicting risk classification.
As we discussed above, scholars have applied multiple panels in helping characterize III/IV nodules in FNA but not specified on aggressive PTCs. Future studies are expected to establish the impact of molecular profiling involving multiple mutations on clinical management of patients with a PTC smaller than 1 cm. Furthermore, characterization of tumor mutation profiles (e.g. exon sequencing) combined with US scan and laboratory results (e.g. immunohistochemistry) might also show a more comprehensive picture of a malignant thyroid nodule and it will also open more choices for therapeutic plans for PTC patients.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Omics Data Encyclopedia (NODE) [accession: OEP002494].

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Shanghai 6th people's hospital institutional review board. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
YL, WG, and XC contributed equally to this work. YCZ and YC obtained tumor and para-carcinoma tissue samples and completed follow-up. AJ performed clinic data statistical analysis. XC performed targeted capture sequencing of PTC genes. YL edited the manuscript. WG modified this article. YYZ and YW designed experiments. BH was the consultant, and supervised the project. XY and TZ performed Bioinformatics Analysis. All authors contributed to the article and approved the submitted version.