Identification and functional enrichment analysis of DEGs involved in THCA

The “Limma” package was used to analyze the DEGs between THCA and normal tissues. A total of 557 up-regulated and 746 down-regulated DEGs were obtained, with ZCCHC12 identified as the most highly expressed gene (Table S2). These results were presented in the volcano and heat maps (Fig. 1A, B). And then, we conducted function analysis based on up-and down-regulated genes respectively. GO functions annotated by the up-regulated genes included synapse organization, skin development, regulation of cell morphogenesis involved in differentiation, regulation of cell morphogenesis, regulation of axonogenesis, etc. (Fig. 1D). While top 20 KEGG metabolic pathways enriched by these genes included p53 signaling pathway, Transcriptional misregulation in cancer, Small cell lung cancer, Pyrimidine metabolism, Proteoglycans in cancer, etc. (Fig. 1C). GO terms annotated by down-regulated genes were urogenital system development, striated muscle tissue development, renal system development, regulation of ossification, regulation of epithelial cell proliferation, etc. (Fig. 1F). And the top 20 KEGG metabolic pathways enriched by down-regulated genes included Wnt signaling pathway, Tyrosine metabolism, Thyroid hormone synthesis, Retinol metabolism, Pyruvate metabolism, etc. (Fig. 1E).

Fig. 1

Screening for differentially expressed genes (DEGs) in THCA and functional enrichment analysis of these DEGs. A Volcano map for DEG screening with volcanoes plotted using fold change and corrected p-values. Red dots in the schematic indicate genes that are significantly differentially up-regulated, blue dots indicate genes that are significantly differentially down-regulated, and grey dots indicate genes that are not significantly regulated. B Heat map for clustering analysis of DEGs with different colors representing the expression trends in different tissues. Due to the large number of DEGs, only 50 up-regulated and 50 down-regulated genes with the greatest differential change are shown here. C KEGG analysis of highly expressed DEGs. D GO analysis of highly expressed DEGs. E KEGG analysis of DEGs with low expression. F GO analysis of DEGs with low expression. The different colors represent the significance of the differential enrichment results with the larger values representing smaller FDR values. The size of the circles represents the number of enriched genes with the larger circles representing the larger number of genes

Expression of ZCCHC12 and its correlation with THCA survival and staging

As Fig. 2A shows, the expression of ZCCHC12 was significantly higher in THCA than in normal controls. And then, we collected 501 THCA patients and categorized them into the high (376) and low (125) groups based on ZCCHC12 expression levels. The Kaplan-Meier curves revealed that those patients with high ZCCHC12 expression showed a higher survival prognosis (P = 1.0E-3; Fig. 2B). And the Sankey diagram (Fig. 2C) showed that the majority of THCA patients were younger than 60 years old, with more female patients than male patients. Patients at tumor stage I accounted for the majority of patients, whereas the smallest proportion of patients was identified in stage II. Notably, most patients with high ZCCHC12 expression are at stage I, while low ZCCHC12 expression patients are at stages II and III. Correspondingly, a higher proportion of patients with low ZCCHC12 expression were among the patients who died (Table S3). What needs to be clarified here is that the expression of ZCCHC12 is higher in THCA patients compared to normal individuals. However, after being diagnosed with THCA, the high expression of ZCCHC12 is actually beneficial for the survival of patients.

Fig. 2

ZCCHC12 expression in THCA and its correlation with survival and staging. AZCCHC12 expression in THCA. B K-M analysis of ZCCHC12 and survival in THCA. C Sankey diagram of age, sex, stage, ZCCHC12 expression, and survival status in THCA. Symbols “****” indicate the significant difference based on P < 0.0001

Correlation of ZCCHC12 with Tumor immunity, Tumor stemness, and TMB, and the PPI network with the ZCCHC gene family

In THCA, the ZCCHC12 expression was not correlated with B-cell expression (P = 0.112; Fig. 3A) but was negatively correlated with CD4 + T cells (P = 0.039; correlation coefficient = − 0.09; Fig. 3B) and positively correlated with CD8 + T cells (P = 1.07e-11; correlation coefficient = 0.30; Fig. 3C), neutrophils (P = 2.41e-04; correlation coefficient = 0.16; Fig. 3D), macrophages (P = 1.85e-13; correlation coefficient = 0.32; Fig. 3E), and dendritic cells (P = 0.02; correlation coefficient = 0.14; Fig. 3F).And the ZCCHC12 expression was also negatively correlated with both tumor stemness (P = 3.82e-15; correlation coefficient = − 0.34; Fig. 3G) and tumor mutational load (P = 1.59e-04; correlation coefficient = − 0.17; Fig. 3H). The PPI network analysis revealed a total of 10 proteins associated with ZCCHC12, including ASPRV1, LAMA3, LAMA4, PNMA1, PNMAL1, RGAG1, RGAG4, SMAD1, SMAD7, and SUMO1 (Fig. 3I, Table S4).

Fig. 3

Correlation analysis of ZCCHC12 with tumor immunity, tumor stemness, and TMB, and the PPI network for the ZCCHC gene family. A-F Correlation of ZCCHC12 with different types of immune cells in THCA. G Correlation of ZCCHC12 with tumor stemness score in THCA. H Correlation of ZCCHC12 with tumor mutational load in THCA. I PPI network of proteins interacting with ZCCHC12. Correlation analyses are performed using Spearman statistics, with the horizontal coordinates representing gene expression and the vertical coordinates representing the immune score distribution. The density curve (right) represents the trend in score distribution, the upper density curve is the trend in gene distribution, and the uppermost value represents the P-value of the correlation analysis

Chromosomal location and MOTIF structure of the ZCCHC gene family and its expression and correlation with THCA

Except for ZCCHC genes that were either not expressed in THCA or unstably expressed, a total of 12 genes were retained, including ZCCHC17 on chromosome 1, ZCCHC4 on chromosome 4, both ZCCHC9 and ZCCHC10 on chromosome 5, ZCCHC7 on chromosome 9, ZCCHC24 on chromosome 10, ZCCHC8 on chromosome 12, ZCCHC14 on chromosome 16, ZCCHC2 on chromosome 18, ZCCHC3 on chromosome 20, and both ZCCHC12 and ZCCHC18 on chromosome X (Fig. 4A). Motif analysis showed that these 12 genes consisted of 15 progenitors and belonged to the same gene family (Fig. 4B). The expressions of ZCCHC genes in THCA were all statistically different, with the most significant difference detected in ZCCHC12 (Fig. 4C). The top 25 genes associated with THCA were identified based on the GENECARDS database with the correlation between the ZCCHC gene family and these top 25 genes analyzed, showing that most of these genes were related to each other (Fig. 4D, Table S5).

Fig. 4

Chromosomal location and MOTIF structure of the ZCCHC gene family and their expressions and correlations with THCA. A Location of the ZCCHC gene family on the chromosomes. B Motif structure of the ZCCHC gene family. C Expression of the ZCCHC gene family in THCA. D Correlation of the ZCCHC gene family with THCA-related genes. *, P < 0.05; **, P < 0.01; ***, P < 0.005

Association of the ZCCHC gene family with THCA survival and staging

Analysis of the line graphs revealed that only ZCCHC3 and ZCCHC7 showed an effect on survival based on the univariate COX, but no effect on their association with the other 10 genes (Fig. 5A). Analysis of the expression of the ZCCHC gene family in different THCA stages revealed that all ZCCHC genes differed at various stages (Fig. 5B). ZCCHC3 shows a gradual decline in expression with the progression of THCA. However, it is similar to ZCCHC12 and exhibits higher expression levels in the THCA group compared to the control group, showed strong potential as a biomarker of disease progression (Table S6).

Fig. 5

Association of the ZCCHC gene family with THCA survival and staging. A Plot of one-way cox analysis of the ZCCHC gene family and survival in THCA. B Differential expression of the ZCCHC gene family at different stages in THCA. The horizontal coordinates represent the different groups of samples (G1 = stage I), the vertical coordinates represent the distribution of the gene expression, the different colors represent different groups, and the numbers in the upper left corner represent the significant p-values. *, P < 0.05; **, P < 0.01; ***, P < 0.005; ****, P < 0.001

Diagnostic efficiency of the ZCCHC gene family for THCA

To evaluate the diagnostic efficiency of the ZCCHC gene family for THCA, the receiver operating characteristic (ROC) curves were constructed (Table S7), showing that the AUC value of ZCCHC2 was 0.81 (Fig. 6A), 0.90 for ZCCHC3 (Fig. 6B), 0.64 for ZCCHC4 (Fig. 6C), 0.78 for ZCCHC7 (Fig. 6D), 0.80 for ZCCHC8 (Fig. 6E), 0.73 for ZCCHC9 (Fig. 6F), 0.64 for ZCCHC10 (Fig. 6G), 0.93 for ZCCHC12 (Fig. 6H), 0.63 for ZCCHC14 (Fig. 6I), 0.64 for ZCCHC17 (Fig. 6J), 0.89 for ZCCHC18 (Fig. 6K), and 0.90 for ZCCHC24 (Fig. 6L).

Fig. 6

ROC curves of the ZCCHC gene family members for THCA, including AZCCHC2, BZCCHC3, CZCCHC4, DZCCHC7, EZCCHC8, FZCCHC9, GZCCHC10, HZCCHC12, IZCCHC14, JZCCHC17, KZCCHC18, and LZCCHC24.

Correlation of the ZCCHC gene family with immunity and variations of Tumor stemness in THCA

In THCA, the immune scores were assessed based on the TIMER database (Table S8), and the correlations between the ZCCHC gene family and immune cells were investigated (Fig. 7A). The results revealed that ZCCHC9 showed positive correlations with B cells, macrophages, neutrophils, and CD8 + T cells; ZCCHC8 was positively correlated with B cells, macrophages, neutrophils, and CD8 + T cells, and negatively correlated with CD4 + T cells; ZCCHC7 was positively correlated with B cells, macrophages, neutrophils, and CD8 + T cells, and negatively correlated with CD4 + T cells; ZCCHC4 was positively correlated with B cells, macrophages, neutrophils, and CD8 + T cells, and negatively correlated with CD4 + T cells; ZCCHC3 was positively correlated with macrophages and CD8 + T cells, and negatively correlated with dendritic cells, neutrophils, and CD4 + T cells; ZCCHC24 showed positive correlation with B cells, macrophages, dendritic cells, neutrophils, and CD8 + T cells; ZCCHC2 showed positive correlation with B cells, macrophages, dendritic cells, neutrophils, CD8 + T cells, and CD4 + T cells; ZCCHC18 was positively correlated with B cells, macrophages, dendritic cells, neutrophils, and CD4 + T cells; ZCCHC17 was positively correlated with macrophages and CD8 + T cells and negatively correlated with dendritic cells, neutrophils, and CD4 + T cells; ZCCHC14 was positively correlated with B cells, macrophages, dendritic cells, neutrophils, and CD8 + T cells and negatively correlated with CD4 + T cells; ZCCHC12 was positively correlated with macrophages, dendritic cells, neutrophils, and CD8 + T cells and negatively correlated with CD4 + T cells; ZCCHC10 was positively correlated with B cells, macrophages, neutrophils, and CD8 + T cells and negatively correlated with CD4 + T cells. The tumor stemness was significantly higher in THCA compared to normal controls (P = 4.1e-154; Fig. 7B). ZCCHC12 was negatively correlated with THCA tumor stemness (P = 3.82e-15; correlation coefficient = − 0.34; Fig. 7C). A pattern plot of tumor stemness with ZCCHC12 expression, age, gender, and the stage was constructed by collating THCA clinical data and tumor stemness scores in TCGA. The results revealed lower ZCCHC12 expressions when the tumor stemness was elevated, more people aged > 60 years old, no significant difference in gender distribution, and a high number of patients at stages higher than I (Fig. 7D).

Fig. 7

Correlation of the ZCCHC gene family with immune cells in THCA and the correlation of tumor stemness differences in THCA with ZCCHC12. A Correlation of the ZCCHC gene family with immune cells. B Tumor stemness variations in THCA. C Tumor stemness variations in THCA and their correlations with ZCCHC12. D Tumor stemness distribution in THCA in relation to age, sex, stage, and ZCCHC12. *, P < 0.05; **, P < 0.01; ***, P < 0.005; **** P < 0.001

GO analysis and PPI network of the ZCCHC gene family

To analyze the function of the ZCCHC gene family (Table S9), the GO annotation analysis was performed based on the ZCCHC gene family, showing the differences in the zinc ion binding, transition metal ion binding, metal ion binding, cation binding, nucleic acid binding, RNA binding, Nucleolus, TRAMP complex, rRNA (adenine) methyltransferase activity, and rRNA methyltransferase activity pathways (Fig. 8A). The results showed that most of the ZCCHC gene family members except for ZCCHC18 were interlinked (Fig. 8B).

Fig. 8

GO analysis (A) and PPI network (B) of the ZCCHC gene family

Validation of ZCCHC gene family genes and proteins expression

We also verified the differential expression of ZCCHC gene family genes in commonly used cell lines for thyroid cancer, and all genes were differentially expressed except ZCCHC3 and ZCCHC18 (Fig. 9A). The results of immunohistochemistry confirmed that ZCCHC2, ZCCHC4, ZCCHC7, ZCCHC8, ZCCHC9, ZCCHC10, ZCCHC12, ZCCHC17, ZCCHC24 were higher in THCA patients than in normal controls; ZCCHC18 was lower in THCA patients than in normal controls (Table S10, Fig. 9B). Before experimental validation, we obtained the ZCCHCs gene expression matrix of THCA cell lines from the CCLE database (https://portals.broadinstitute.org/ccle/about). Through analysis, we found that differences in gene expression between cells and gene expression in patients are commonly observed. TPC-1 and Nthy-ori-3 cells are widely used cell lines for THCA research, and we believe it is reasonable to use them to validate the differential expression of the ZCCHC gene family.

Fig. 9

Validation of ZCCHC gene family genes and proteins expression. A ZCCHC gene family genes are differentially expressed in TPC1 and NTHYORI 3 − 1 cells. B Differential expression of immunohistochemistry of the ZCCHC gene family in The Human Protein Atlas database