Pathogenically high-risk variants among 31 driver genes and non-driver genes were detected in both cancer tissues and benign thyroid nodules

In total, 592 mutations were detected, including missense mutations (426/592, 71.96%), frame-shift insertions (27/592, 4.56%), frame-shift deletions (17/592, 2.87%), nonsense mutations (6/592, 1.01%) and synonymous variant (116/592, 19.59%). The pathogenic loci were further screened with ClinVar, and the pathogenicity of mutation sites was predicted to belong to variations of uncertain significance, among which probably harmful (115/476, 24.16%) and possibly harmful (57/476, 11.97%) loci were denoted pathogenically high-risk variants. A total of 172 pathogenically high-risk variants were detected, including somatic mutations and germline mutations. These variants were distributed among 31 genes in two major categories: driver genes (13 genes) and non-driver genes (18 genes). Driver genes included gain-of-function oncogenes (BRAF, RET, CTNNB1, NRAS, HRAS, KRAS, AKT1, PIK3CA, TERT and GNAS) and loss-of-function oncogenes (TP53, PTEN and APC). Non-driver genes included DNA damage repair genes (CHEK2, ATM and BRCA1), immune-associated genes (HLA-A, MUC6 and IGSF3), chromosomal modifier genes (KMT2C, PRAMEF2 and NEFH), membrane transporter-associated genes (OBSL1 and SYVN1), hormone-associated genes (TSHR), metabolism-associated genes (GGT1), keratin-associated genes (KRTAP4-8), ankyrin-associated genes (POTED), genes of the otopetrin domain protein family (OTOP1), an FSHD region gene family member (FRG2C) and a tubulin polyglutamylase complex subunit (TPGS2) (Fig. 1A, B).

To further explore the presence of germline variations, mutations were screened with the conditions that variant allele frequency was over 40% and mutations were present in both cancer tissues and para-carcinoma tissues. RET mutation (c. G2071A: p.G691S) was detected, which suggested the possibility of hereditary medullary thyroid carcinoma. Besides, germline mutations in POTED, TPGS2, IGSF3, KMT2C, FGR2C, MUC6, PRAMEF2, HLA-A, GNAS, OBSL1, OTOP1, NEFH, TSHR and CHEK2 were also detected (Table S5). Sequencing data for cancer tissues and benign nodules were filtered with reference to para-carcinoma tissues. And a total of 156 somatic mutations were detected, among which the top ten most frequently mutated genes were BRAF, ATM, GGT1, TSHR, KMT2C, CHEK2, PIK3CA, TERT, OTOP1 and RET.

Furthermore, we compared the mutation profiles between cancer tissues and benign nodules. BRAFV600E occurred more frequently in cancer tissues than benign nodules (92/124, 74.19% vs. 10/58, 17.24%, p < 0.0001), whereas BRAF mutations at other loci (such as the Raf-like Ras-binding domain) were detected primarily in benign nodules (10/58, 17.24% vs. 1/124, 0.81%, p < 0.0001). Mutations in TERT gene were dispersed over the entire sequence, but occurred primarily in the telomerase ribosomal nuclear protein complex-RNA binding domain (25/124, 20.16%). The frequency of TERT gene mutations was significantly higher in cancer tissues than benign nodules (27/124, 21.77% vs. 2/58, 3.45%, p = 0.094), thus suggesting the value of this gene in distinguishing benign from malignant thyroid nodules. Moreover, RET mutations accounted for a significantly higher proportion of mutations in cancer tissues than benign nodules (10/124, 8.06% vs. 1/58, 1.72%, p = 0.094). Furthermore, the mutation frequencies of DNA damage repair genes ATM and metabolism-associated genes GGT1 were significantly higher in cancer tissues than benign nodules (Fig. 1C). Furthermore, more mutation sites of genes in DNA damage repair signaling were observed in cancer tissues than benign nodules (Fig. 1D).

Gene fusion is another gene alteration type in thyroid carcinoma [28,29,30,31,32,33,34]. Gene fusions were detected at both the RNA and DNA levels. At the DNA level, 2 patients (2/124, 1.61%) had fusions in KRTAP4-8, KRTAP4-4 and NEFH. CCDC6-RET (2/124, 1.61%), NCOA4-RET (6/124, 4.84%), ETV6-NTRK3 (5/124, 4.03%) and TPM3-NTRK1 (1/124, 0.81%) were detected at the RNA level, and 13 were patients with cancer tissues (Fig. 1E). And fusions were the most common fusion types reported in previous studies (Fig. 1F).

Comparison of mutation profiles of PTC tissues between in-house and TCGA data

Mutated genes were compared between cancer tissues and the thyroid carcinoma cohort from The Cancer Genome Atlas (TCGA) database. The mutation frequencies in BRAF, NRAS, CHEK2, ATM, NRAS and HRAS were relatively high in both TCGA cohort and our cohort. In addition to the high frequency mutated genes in TCGA, the mutation frequencies of PIK3CA were high in cancer tissues from our cohort, thus suggesting possible activation of the PI3k/Akt/mTOR pathway. Moreover, the mutation frequencies of GGT1, TSHR, KMT2C, TERT, OTOP1 and RET were also relatively high in our cohort (Fig. 1G).

Fig. 1

Profiles of genetic alterations significantly differ between thyroid carcinomas and benign nodules. (A) Overall analysis of gene alteration profiles, indicating that most mutated genes have significantly higher mutation ratios in cancer tissues than benign nodules. (B) The number of each particular type of variations in detected genes. (C) Mutation frequencies of all genes sequenced among cancer tissues, benign nodules and para-carcinoma tissues. (D) Mutation sites of genes in DNA damage repair signaling pathways. (E) Significantly higher frequency of gene fusions detected at the RNA level than the DNA level. (F) Fusion sites of the four fusion types: CCDC6-RET (2/124, 1.61%) is the fusion of CCDC6 exon 1 and RET exon 12, NCOA4-RET (6/124, 4.84%) is the fusion of NCOA4 exon 8 and RET exon 12, ETV6-NTRK3 (5/124, 4.03%) is the fusion of ETV6 exon 4 and NTRK3 exon 14, and TPM3-NTRK1 (1/124, 0.81%) is the fusion of TPM3 exon 7 and NTRK1 exon 10. (G) Mutation proportion of mutated genes in our cohort and TCGA database

Concomitant pathogenically high-risk variants are detected more frequently in PTC tissues than benign nodules and correlate with lymph node metastasis

We analyzed the co-mutation status of gene mutations and fusions in PTC tissues. BRAFV600E was concomitant with ATM (37.90%, p = 0.052), TERT (20.97%, p = 0.003) and CHEK2 (23.39%, p = 0.004). TERT was co-mutated with ATM (14.52%, p = 0.025), TSHR (13.71%, p < 0.0001), CHEK2 (12.9%, p < 0.0001), RET (4.03%, p = 0.024), GGT1 (14.52%, p < 0.0001), PIK3CA (8.87%, p = 0.011) and BRCA1 (0.81%, p = 0.057). RET was co-mutated with PIK3CA (4.03%, p = 0.031), PTEN (2.42%, p = 0.004), AKT1 (0.81%, p = 0.028), CHEK2 (4.03%, p = 0.057) and BRCA1 (0.81%, p = 0.001). Furthermore, ATM (8.87%, p = 0.005), GNAS (0.81%, p = 0.003), HRAS (0.81%, p = 0.066) and KRTAP4-8 (0.81%, p = 0.003) were frequently co-mutated with gene fusions (Fig. 2A).

Correlation analysis of gene alterations and lymph node metastasis (LN) showed that RETMUT (p = 0.033), ATMMUT (p = 0.031) and TERTMUT (p = 0.005) were positively correlated with lymph node metastasis. Meanwhile, although the correlation between lymph node metastasis and BRAFV600E (p = 0.136), as well as PIK3CAMUT (p = 0.103), was not significant, they were also likely to be associated with lymph node metastasis (Table S6). We further analyzed the correlation between the genes above co-mutant with other genes and lymph node metastasis. BRAFV600E and ATMMUT or TERTMUT were highly correlated with lymph node metastasis (BRAFV600E and ATMMUT, p = 0.004; BRAFV600E and TERTMUT, p = 0.013). TERTMUT concomitant with ATMMUT, CHEK2MUT, GGT1MUT and PIK3CAMUT was highly associated with lymph node metastasis (TERTMUT and ATMMUT, p = 0.001; TERTMUT and CHEK2MUT, p = 0.011; TERTMUT and GGT1MUT, p = 0.019; TERTMUT and PIK3CAMUT, p = 0.003). PIK3CAMUT concomitant with GGT1MUT and CHEK2MUT was highly associated with lymph node metastasis (PIK3CAMUT and GGT1MUT, p = 0.03; PIK3CAMUT and CHEK2MUT, p = 0.017). Concomitant RETMUT and CHEK2MUT were highly correlated with lymph node metastasis (p = 0.052). Besides, concomitant ATMMUT and gene fusions were correlated with lymph node metastasis, with marginal significance (p = 0.143) (Fig. 2B). Therefore, co-mutations of driver genes had better performance than single mutations in predicting lymph node metastasis.

Fig. 2

Concomitant pathogenically high-risk variants are detected more frequently in cancer tissues than benign nodules and are correlated with lymph node metastasis. (A) Analysis of correlations among mutated genes in cancer tissues. (B) Frequency of concomitant gene alterations in cancer tissues with or without lymph node metastasis

ATM deficiency in PTC indicates high proliferation and aggressiveness of carcinoma cells

We classified 31 mutant genes in pathway enrichment analysis. A total of 11 pathways were enriched: thyroid hormone signaling, apoptosis, mTOR signaling, PI3K-Akt signaling, autophagy, VEGF signaling, MAPK signaling, p53 signaling, cAMP signaling, Ras signaling and homologous recombination (Figure S3A). In the analysis of the frequencies of mutated genes in each pathway, the mutation rates of genes in apoptosis (p < 0.0001), cAMP signaling (p < 0.0001), mTOR signaling (p < 0.0001), homologous recombination ( p < 0.0001), p53 signaling (p < 0.0001) and MAPK signaling (p < 0.0001) were significantly higher in cancer tissues than benign nodules (Figure S3B).

Furthermore, we simultaneously conducted pathway enrichment analysis between tumors with and without lymph node metastasis. Genes associated with apoptosis (p = 0.002), p53 signaling (p = 0.022) and homologous recombination (p = 0.021) were mutated more frequently in tumors with than without lymph node metastasis (Fig. 3A). The homologous recombination repair gene ATM existed in all three pathways and correlated with lymph node metastasis. In the cohort, we identified a new ATM gene mutation absent from TCGA database. The mutation was located in the FAT domain, a unique motif at the extreme C terminus in PIK-related family members (Fig. 3B), and was predicted to be pathogenic. The formalin fixed paraffin embedded samples of 48 cancer tissues were collected and stained to detect ATM protein expression. Low expression of ATM was observed frequently in samples with ATMMUT and was significantly associated with lymph node metastasis (Fig. 3C). And ATM mutation was highly correlated with lymph node metastasis in cancer tissues (Fig. 3D).

To explore the effects of ATM on the pro-tumoral effects of pathogenic mutations in driver genes, we established ATM knock-down thyroid cancer cell lines (B-CPAP, TPC-1, Cal-62 and Nthy-ori-3-1) (Figure S4). Compared with NC-transfected cells, B-CPAP, TPC-1, Cal-62 and Nthy-ori-3-1 cells had significantly increased viability after knock-down of ATM expression (Fig. 3E). After knock-down of ATM expression, the number of apoptotic cells significantly decreased (B-CPAP vs. B-CPAPsi−ATM = 10.69% ± 0.48% vs. 3.96% ± 0.06%, p < 0.0001; Cal-62 vs. Cal-62si − ATM = 9.56% ± 0.42% vs. 2.18% ± 0.13%, p < 0.0001; TPC-1 vs. TPC-1si − ATM = 12.15% ± 0.25% vs. 4.81% ± 0.62%, p = 0.0001; Nthy-ori-3-1 vs. Nthy-ori-3-1si − ATM = 11.54% ± 0.43% vs. 8.06% ± 0.36%, p = 0.0009; Fig. 3F). The invasiveness of cells was also explored. The invasive ability of B-CPAPsi−ATM, Cal-62si − ATM, TPC-1si − ATM and Nthy-ori-3-1si − ATM cells was significantly greater than that of NC-transfected cells (B-CPAP vs. B-CPAPsi−NC vs. B-CPAPsi−ATM = 56.33 ± 4.5 vs. 63.33 ± 8.73 vs. 105 ± 6.53, p = 0.0057; Cal-62 vs. Cal-62si − NC vs. Cal-62si − ATM = 104.67 ± 4.99 vs. 105.33 ± 4.64 vs. 212.33 ± 7.85, p < 0.0001; TPC-1 vs. TPC-1si − NC vs. TPC-1si − ATM = 163 ± 8.6 vs. 137.33 ± 8.65 vs. 212 ± 13.14, p = 0.0026; Nthy-ori-3-1 vs. Nthy-ori-3-1si − NC vs. Nthy-ori-3-1si − ATM = 137 ± 12.33 vs. 132 ± 6.68 vs. 236.67 ± 13.27, p = 0.0026; Fig. 3G). Thus, ATM deficient thyroid cancer cells bearing BRAFV600E, KRASG12R or CCDC6-RET displayed higher proliferative ability, lower apoptotic rate and greater aggressiveness than NC-transfected cells in vitro.

Fig. 3

Deficiency in the homologous recombination repair gene ATM predicts elevated proliferation and aggressiveness in thyroid carcinoma. (A) Mutations in genes associated with apoptosis, the p53 signaling pathway and the homologous recombination pathway are more frequent in tumors with lymph node metastasis. (B) Mutation sites of ATM are focused in the FAT domain. (C) Low expression of ATM is relatively more frequent in samples with ATMMUT and is significantly associated with lymph node metastasis. (D) ATM mutation was correlated with lymph node metastasis in cancer tissues. (E) Significantly greater relative viability of B-CPAP, TPC-1, Cal-62 and Nthy-ori-3-1 was shown after ATM knock-down than in NC-transfected cells. (F) Significantly decreased apoptosis was shown after decreasing ATM expression. (G) Significantly greater invasive capability of B-CPAPsi−ATM, Cal-62si − ATM, TPC-1si − ATM and Nthy-ori-3-1si − ATM was shown than NC-transfected cells

Combining gene alterations and clinical features with ultrasonographic radiomics significantly improves the predictive efficacy for lymph node metastasis in PTC

To explore the effect of radiomics on predicting lymph node metastasis in our cohort, we extracted 854 radiomic features from 108 PTC samples. The Glmnet package was used to filter the smallest regularization parameter lambda. The model was fitted by minimizing a combination of the loss function and the regularization term (Fig. 4A). Area under curve (AUC) with standard errors can be used to determine the optimum penalty lambda (λ) by cross-validation (Fig. 4B). Three features (wavelet-LHL_glrlm_Short Run Emphasis, wavelet-HLH_glszm_Size Zone Non-Uniformity and wavelet-HLH_glszm_Small Area Low Gray Level Emphasis) associated with lymph node metastasis were identified, and a radiomics signature model for predicting lymph node metastasis of PTC was constructed integrating these three features. Furthermore, a radiomics score (Radscore) was calculated for each sample based on the features selected. Cancer tissues with lymph node metastasis had significantly higher Radscores than those without lymph node metastasis in the training set, and cancer tissues with lymph node metastasis also tended to have higher Radscores in the test set (Fig. 4C, Table S7).

Fig. 4

Feature selection using the least absolute shrinkage and selection operator (LASSO) algorithm. (A) A coefficient path plot was generated showing how the coefficients of each variable changed at different regularization levels. (B) AUC (red dots) with standard errors (error bar) can be used to determine the optimum penalty lambda (λ). (C) Boxplot of Radscores between cancer tissues with or without lymph node metastasis from the training set (left) and the test set (right)

Next, we explored the efficacy of gene alterations, clinical features and radiomics in the prediction of lymph node metastasis. Mutated genes associated with lymph node metastasis in the univariate analysis, including RET, ATM, PIK3CA, TERT, GGT1 and fusions, were incorporated to construct a gene signature model. And a clinical feature model based on clinical characteristics (including age, gender, extrathyroidal extension, Hashimoto’s thyroiditis, tumor stage, family history and tumor size) was also constructed. The multi-feature integration nomogram model was then optimized integrating gene alterations, clinical features and radiomic features (Fig. 5A). The calibration plots of multi-feature integration nomogram model showed good calibration in the training set (Fig. 5B, C). We compared the lymph node metastasis prediction efficacy of the gene signature model, clinical feature model, radiomics signature model, their pairwise combination and the multi-feature integration nomogram model. The AUC of the multi-feature integration nomogram model (0.8695, 95% CI:0.7909–0.9482) was higher than that of the gene signature model (0.7649, 95% CI: 0.6614–0.8683), clinical feature model (0.7, 95% CI: 0.5826–0.8174), and radiomics signature model (0.7147, 95% CI: 0.6003–0.8291), as well as their pairwise combinations (Fig. 5D, Figure S5A, Table 1). These findings indicated that combining gene alterations and clinical features with radiomic signatures increased the efficacy of the model for predicting lymph node metastasis of thyroid carcinoma. The Decision Curve Analysis (DCAs) of the training data based on the seven models are presented in Fig. 5E (DCAs of the test data are presented in Figure S5B). For the prediction of lymph node metastasis of PTC, the multi-feature integration nomogram model had the highest overall net benefit.

Fig. 5

Construction of a multi-feature integration nomogram model and comparison of the seven models. (A) The multi-feature integration nomogram model, integrating radiomic features, mutated genes and clinical features, developed in the training set. (B) Calibration curves for the multi-feature integration nomogram model in the training set. (C) Calibration curves for the multi-feature integration nomogram model in the test set. (D) ROC curves of the seven models in the training set. (E) DCAs for the seven models in the training set

Table 1 The AUC of models predicting lymph node metastasis.