Comprehensive Analysis of a Zinc Finger Protein Gene–Based Signature with Regard to Prognosis and Tumor Immune Microenvironment in Osteosarcoma

Osteosarcoma is the most common malignant bone tumor that seriously threatens the lives of teenagers and children. Zinc finger (ZNF) protein genes encode the largest transcription factor family in the human genome. Aberrant expressions of ZNF protein genes widely occur in osteosarcoma, and these genes are therefore attractive biomarker candidates for prognosis prediction. In this study, we conducted a comprehensive analysis of ZNF protein genes in osteosarcoma and identified prognosis-related ZNF protein genes. Then, we constructed a prognostic signature based on seven prognosis-related ZNF protein genes and stratified patients into high- and low-risk groups. The seven genes included MKRN3, ZNF71, ZNF438, ZNF597, ATMIN, ZNF692, and ZNF525. After validation of the prognostic signature in internal and external cohorts, we constructed a nomogram including clinical features such as sex and age and the relative risk score based on the risk signature. Functional enrichment analysis of the risk-related differentially expressed genes revealed that the prognostic signature was closely associated with immune-related biological processes and signaling pathways. Moreover, we found significant differences between the high- and low-risk groups for the scores of diverse immune cell subpopulations, including CD8⁺ T cells, neutrophils, Th1 cells, and TILs. Regarding immune function, APC co-inhibition, HLA, inflammation promotion, para-inflammation, T-cell co-inhibition, and the type I IFN response were significantly different between the high- and low-risk groups. Of the seven ZNF protein genes, lower expressions of ATMIN, MKRN3, ZNF71, ZNF438, and ZNF597 were correlated with a high risk, while higher expressions of ZNF525 and ZNF692 were associated with a high risk. The Kaplan–Meier survival analysis suggested that lower expressions of ATMIN, ZNF438, and ZNF597 and the higher expression of ZNF692 were associated with worse overall survival in osteosarcoma. In conclusion, our ZNF protein gene–based signature was a novel and clinically useful prognostic biomarker for osteosarcoma patients.

Introduction

Osteosarcoma is the most common malignant bone tumor that seriously threatens the lives of teenagers and children, with a second incidence peak in the elderly over the age of 60 years (Ritter and Bielack, 2010). In the last few decades, huge progresses in the diagnosis and management of osteosarcoma have been made (Kager et al., 2017; Chen et al., 2020). Surgical resection combined with neo-adjuvant chemotherapy is currently the standard treatment for patients with osteosarcoma, and it increases the 5-year overall survival from 20 to 70% for patients with localized disease (Kansara et al., 2014; Harrison et al., 2018). However, the 5-year overall survival is still less than 20% for patients with metastatic or recurrent osteosarcoma and has remained almost unchanged over the past 30 years (Kansara et al., 2014). Additionally, resistance to chemotherapy regimens often occurs due to the high heterogeneity in the genome of osteosarcoma, and it seriously impacts the prognosis of patients with osteosarcoma (Gianferante et al., 2017; Sayles et al., 2019). Despite the existence of some prognostic biomarkers in osteosarcoma, it is still unachievable to accurately predict the prognosis of patients with osteosarcoma. Therefore, the development of novel risk stratification methods, as well as identification of novel prognostic biomarkers, are clinically useful.

The zinc finger (ZNF) protein genes encode the largest family of regulatory proteins with zinc ion–binding finger-like domains in the human genome (Laity et al., 2001). The functions of these ZNF proteins are extremely diverse, and they can interact with RNA, DNA, or proteins to participate in various cellular processes including transcriptional regulation, mRNA stability, and protein degradation (Leon and Roth, 2000; Fu and Blackshear, 2017; Wang and Zheng, 2021). Over the last decades, increasing number of researches have revealed the vital role of ZNF proteins in oncogenesis and its progression (Jen and Wang, 2016; Cassandri et al., 2017; Ye et al., 2019). Aberrant expressions of ZNF proteins or genes have been found in malignant tumors and might serve as diagnostic or prognostic biomarkers. For example, Poirier et al. (2021) reported that ZNF768 was elevated in lung adenocarcinoma when compared to normal lung tissue and was associated with clinicopathological features. In vitro experiments suggested that depletion of ZNF768 severely impairs proliferation in several lung cancer cell lines. A further example includes ZNF471, which acts as a tumor suppressor in cervical cancer. Hypermethylation of the ZNF471 gene promoter and low ZNF471 expression was correlated with poor overall and recurrence-free survival in cervical cancer (Bhat et al., 2021).

In this study, we conducted a comprehensive analysis of ZNF protein genes in osteosarcoma. Prognosis-related ZNF protein genes were identified and further used to construct a prognostic signature in osteosarcoma. We found that our novel ZNF protein gene–based signature could act as a risk stratification method and a prognostic biomarker. Meanwhile, we also developed a prognostic nomogram comprising clinical features such as gender and age and the relative risk score based on the risk signature. Additionally, we found that our ZNF protein gene–based signature was closely associated with immune cell infiltration and immune function. Taken together, our ZNF protein gene–based signature was a novel and clinically useful prognostic biomarker for osteosarcoma patients, and our findings shed light on the potential role of ZNF proteins in the pathogenesis of osteosarcoma.

Materials and Methods

Data Collection and Zinc Finger Protein Genes

The RNA-seq data and corresponding clinical information of patients with osteosarcoma were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The GSE21257 data set, containing the gene expression matrix data from a total of 55 osteosarcoma samples and clinical data, was collected from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) for external validation. Patients without follow-up time or incomplete clinical data were excluded in the subsequent analysis. The list of ZNF protein genes was obtained from the UniProt database (https://www.uniprot.org/), and the ZNF protein genes were adopted to a unified terminology in the gene expression matrixes from the TCGA and the GEO databases.

Identification of Prognosis-Related Zinc Finger Protein Genes in Osteosarcoma

Prognosis-related ZNF protein genes were identified by performing univariate Cox regression analysis. Genes with a p-value <0.05 in the univariate Cox regression analysis were regarded as prognosis-related ZNF protein genes.

Construction and Validation of a Zinc Finger Protein Gene–Based Signature

The TCGA osteosarcoma cohort, hereafter referred to as the entire cohort, was randomly separated into a training cohort and a testing cohort at a ratio of approximately 1:1. Then, in the training cohort, the prognosis-related ZNF protein genes were subjected to the least absolute shrinkage and selection operator (LASSO) regression analysis to eliminate collinearity among these genes. Subsequently, the multivariate Cox regression analysis was conducted using the survival package in R to calculate the regression coefficients, and we ultimately constructed a ZNF protein gene–based signature in training cohort. The risk scores of the patients in the testing cohort, entire cohort, and GSE21257 cohort were calculated based on the signature and the formula as follows:

$\text{risk score }={\displaystyle \sum }_{i=1}^n\left(Coefi\ast Expi\right)$

The Coefi is the regression coefficients of the selected genes, and the Expi is the expression levels of the selected genes. Using the median risk score value in the training cohort, patients in the testing cohort, entire cohort, and GSE21257 cohort were classified into high- and low-risk groups. The clinicopathological characteristics of the osteosarcoma patients in the TCGA entire cohort and the GSE21257 cohort are shown in Table 1. The Kaplan–Meier survival analysis was performed using the survival and survminer packages in R to compare the overall survival of patients in the high- and low-risk groups. Time-dependent receiver operating characteristic (ROC) curve analysis was conducted using the survivalROC package in R to verify the accuracy of the ZNF protein gene–based signature in predicting the prognosis of patients with osteosarcoma.

TABLE 1

TABLE 1. Clinicopathological characteristics of osteosarcoma patients in high- and low-risk groups.

Construction and Validation of a Prognostic Nomogram

A prognostic nomogram, comprising clinical features such as gender and age, and the relative risk score based on the risk signature were constructed using the data set from the entire cohort. Calibration curves were plotted in both the entire cohort and GSE21257 cohort to evaluate the accuracy of the nomogram in predicting 1-, 2-, and 3-year overall survival.

Functional Enrichment Analysis

The differentially expressed genes (DEGs) between the high- and low-risk groups in the gene expression matrix of the TCGA entire cohort were identified using the limma package in R. The genes with a threshold of | log2 (fold change) | > 0.5 and false discovery rate (FDR) <0.05 were regarded as DEGs and were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the clusterProfiler and org.Hs.eg.db packages in R.

Gene Set Enrichment Analysis

The gene set enrichment analysis (GSEA) was performed between the high- and low-risk groups in the gene expression matrixes from the TCGA osteosarcoma cohort and the GSE21257 cohort using the GSEA 4.0.2 software. The enriched pathways were filtered with a threshold of nominal p-value <0.05 and normalized enrichment score (NES) >1.0.

Single-Sample Gene Set Enrichment Analysis

The single-sample gene set enrichment analysis (ssGSEA) was performed using the GSVA package in R to quantify the level of immune cell infiltration and immune function. The patients in the entire cohort and GSE21257 cohort were imported to ssGSEA, and the scores of the immune cell and immune function in each sample were calculated and were further compared between the high- and low-risk groups.

Statistical Analysis

All the statistical analyses were performed using R software version 4.1.0. The differences were considered to be significant at a p-value <0.05 if not specified.

Results

Identification of Prognosis-Related Zinc Finger Protein Genes in Osteosarcoma

The flowchart of the present study is exhibited in Figure 1. First, we performed univariate Cox regression analysis to identify prognosis-related ZNF protein genes in osteosarcoma. As shown in Figure 2A, a total of 41 ZNF protein genes were significantly associated with the prognosis of osteosarcoma patients in the univariate Cox regression analysis, and these genes were identified as prognosis-related ZNF protein genes in osteosarcoma. Of these 41 ZNF protein genes, 16 were risk factors (hazard ratio (HR) > 1), while 25 were protective factors (HR < 1) in osteosarcoma. The expression profile of these 41 prognosis-related ZNF protein genes in osteosarcoma and the correlation of them are exhibited in Figures 2B,C.

FIGURE 1

FIGURE 1. Workflow diagraph of data collection and analysis.

FIGURE 2

FIGURE 2. Identification of prognosis-related ZNF protein genes in osteosarcoma. (A) Prognosis-related ZNF protein genes were identified by performing univariate Cox regression analysis. (B) The expression profiles of the prognosis-related ZNF protein genes. (C) The correlation of the prognosis-related ZNF protein genes based on gene expression.

Construction of a Zinc Finger Protein Gene–Based Signature in Training Cohort

In the training cohort, we performed LASSO regression analysis to eliminate overfitting and then conducted multivariate Cox regression analysis to calculate the regression coefficients of the selected genes (Figures 3A,B). Finally, we developed a novel prognostic signature based on seven ZNF protein genes in osteosarcoma. The seven genes included MKRN3, ZNF71, ZNF438, ZNF597, ATMIN, ZNF692, and ZNF525, and the coefficients of these are shown in Figure 3C. The risk scores of the patients in the training cohort were calculated using the following formula: risk score = MKRN × (−4.202) + ZNF71 × (−2.800) + ZNF438 × (−2.102) + ZNF597 × (−1.675) + ATMIN × (−0.343) + ZNF692 × 0.420 + ZNF525 × (0.652), and it allowed patients to be classified into high- and low-risk groups based on the median value of the risk score. The risk score distribution of patients in the high- and low-risk groups is exhibited in Figure 3D. The survival time of the patients in the high-risk group tended to be worse than that in the low-risk group (Figure 3E). Moreover, the Kaplan–Meier survival analysis revealed that the survival probability of the patients in the high-risk group was dramatically lower than that in the low-risk group (Figure 3F). The time-dependent ROC curves were plotted, and the areas under the curve (AUC) for 1-, 2-, and 3-year survival were 0.890, 0.992, and 0.991, respectively (Figure 3G).

FIGURE 3

FIGURE 3. Construction of a ZNF protein gene–based signature in osteosarcoma training cohort. (A, B) Lasso regression and multivariable Cox regression analyses were performed to identify genes for the prognostic signature. (C) The coefficients of the seven genes. (D) Risk score distribution. (E) Survival time and survival status of patients in high- and low-risk groups. (F) Kaplan–Meier survival analysis in high- and low-risk groups. (G) Time-dependent ROC curves for 1-, 2-, and 3-year survival.

Validation of the Zinc Finger Protein Gene–Based Signature in Internal and External Cohorts

To verify the credibility of our ZNF protein gene–based signature in predicting the prognosis of patients with osteosarcoma, the internal cohorts, including the testing cohort and entire cohort, were employed. Using the same aforementioned formula, the risk scores were calculated, and it allowed patients to be stratified into the high- and low-risk groups with the median risk score in the training cohort as the cutoff value. The risk score distribution of the high- and low-risk patients in the testing and entire cohorts is exhibited in Figures 4A,B. The survival time and status of the patients in the testing and entire cohort are shown in Figures 4C,D, and it seems that the mortality rate of the patients in the high-risk group was higher than in the low-risk group. The expression profiles of the seven ZNF protein genes are displayed in Figures 4E,F. The Kaplan–Meier survival curves showed that the prognosis of the patients in the high-risk group was worse than in the low-risk group in both the testing and entire cohorts (Figures 4G,H). The AUC values for 1-, 2-, and 3-year survival in the testing cohort were 0.751, 0.748, and 0.839, respectively (Figure 4I), and were 0.788, 0.888, and 0.914, respectively, in the entire cohort (Figure 4J). To further evaluate the predictive power of the ZNF protein gene–based prognostic signature, the Kaplan–Meier survival analyses were performed in various subgroups stratified by gender and age, and the results suggested that in all these subgroups, the overall survival of the patients in the high-risk group was worse than in the low-risk group (Figures 5A–D).

FIGURE 4

FIGURE 4. Validation of the ZNF protein gene-based signature in osteosarcoma testing and entire cohorts. (A, B) Risk score distribution of patients in testing and entire cohorts. (C, D) Survival time and survival status of patients in testing and entire cohorts. (E, F) The expression profile of the seven ZNF protein genes in testing and entire cohorts. (G, H) Kaplan–Meier survival analysis of high- and low-risk patients in testing and entire cohorts. (I, J) Time-dependent ROC curves for 1-, 2-, and 3-year survival in testing and entire cohorts.

FIGURE 5

FIGURE 5. Kaplan–Meier survival analysis in subgroups stratified by gender (A, B) or age (C, D) in the TCGA cohort.

To assess the general applicability of the ZNF protein gene–based prognostic signature, the external GSE21257 cohort was enrolled in the subsequent analysis. The risk scores of the patients in the GSE21257 cohort were calculated, and the patients were then divided into the high- and low-risk groups according to the median risk score in the training cohort (Figure 6A). The survival time and survival status of the patients in the GSE21257 cohort are exhibited in Figure 6B, and the result indicated that the patients of the high-risk group harbored a higher mortality rate than those in the low-risk group. The expression profiles of the seven ZNF protein genes in the GSE21257 gene expression matrixes are shown in Figure 6C. The Kaplan–Meier survival analysis suggested that the overall survival of the patients in the high-risk group continued to be worse than in the low-risk group (Figure 6D). The AUC values for 1-, 2-, and 3-year survival were 0.811, 0.831, and 0.737, respectively (Figure 6E). Taken together, these internal and external cohorts confirmed the satisfactory accuracy and generalizability of the ZNF protein gene–based prognostic signature.

FIGURE 6

FIGURE 6. Validation of the ZNF protein gene-based signature in GSE21257 cohort. (A) Risk score distribution of patient in the GSE21257 cohort. (B) Survival time and survival status of patients in the GSE21257 cohort. (C) The expression profile of the seven ZNF protein genes in GSE21257 cohort. (D) Kaplan–Meier survival analysis of high- and low-risk patients in the GSE21257 cohort. (E) Time-dependent ROC curves for 1-, 2-, and 3-year survival in the GSE21257 cohort.

The Prognostic Independence of the Zinc Finger Protein Gene–Based Signature and Construction of a Nomogram in Osteosarcoma

The univariate and multivariate Cox regression analyses were utilized to evaluate the independent predictive power of the ZNF protein gene–based signature. As shown in Table 2, the results indicated that the risk score based on the seven-gene signature was the only independent prognostic biomarker in osteosarcoma. Moreover, we constructed a prognostic nomogram using the gene expression matrixes of the entire cohort and corresponding clinical information. The total point score was calculated based on the relative risk score and clinical features including gender and age. The probability of the 1-, 2-, and 3-year overall survival was reflected by the total points and is shown in Figure 7A. To evaluate the validity of the nomogram, the calibration plots were visualized in both TCGA entire cohort and the GSE21257 cohort. The results showed that the predicted overall survival at 1, 2, and 3 years was consistent with the actual overall survival at 1, 2, and 3 years (Figures 7B,C), suggesting that the nomogram had great accuracy in predicting the overall survival in osteosarcoma.

TABLE 2

TABLE 2. Univariable and multivariable analysis of the ZNF protein gene–based signature and clinical factors in the TCGA osteosarcoma cohort.

FIGURE 7

FIGURE 7. Construction and validation of a nomogram in TCGA and GSE21257 cohorts. (A) Nomogram of clinical feature including age and gender, and relative risk based on the signature. (B, C) Calibration plots of the nomogram for predicting the probability of overall survival at 1, 2, and 3 years in TCGA and GSE21257 cohorts.

Identification of Risk-Related Differentially Expressed Genes and Functional Enrichment Analysis

To reveal the underlying mechanism related to the ZNF protein gene–based signature, we identified genes that were differentially expressed between the high- and low-risk groups using the limma package in R. A total of 1,819 DEGs were screened out with the criteria of | log2 (fold change) | > 0.5 and an FDR < 0.05 (Figure 8A). Of the 1,819 DEGs, 1,178 genes were upregulated and 641 genes were downregulated in the high-risk group. The expression profiles of these 1,819 genes in the high- and low-risk groups are exhibited in Figure 8B. Then, we performed GO and KEGG enrichment analyses on the 1,819 DEGs. In biological processes, DEGs mainly participate in the MHC class II protein complex, immune response–activating signal transduction, humoral immune response, lymphocyte-mediated immunity, and receptor-mediated endocytosis. In terms of molecular functions, DEGs are primarily enriched in signaling receptor activator activity, receptor ligand activity, antigen binding, glycosaminoglycan binding, and structural constituent of the ribosome. In cellular components, enriched terms were mainly related to the external side of the plasma membrane, collagen-containing extracellular matrix, and immunoglobulin complex (Figure 8C). The KEGG enrichment analysis revealed that the DEGs were primarily enriched in cytokine–cytokine receptor interaction, the ribosomes, the phagosomes, transcriptional mis-regulation in cancer, and cell adhesion molecules (Figure 8D).

FIGURE 8

FIGURE 8. Identification of risk-related differentially expressed genes and functional enrichment analysis. (A) Vocal plot exhibited risk-related differentially expressed genes. (B) Expression profiles of genes differentially expressed between high- and low-risk groups in gene matrixes of the TCGA cohort. (C, D) GO and KEGG enrichment analysis of risk-related differentially expressed genes.

Gene Set Enrichment Analysis and Immune Score Analysis

The gene set enrichment analysis (GSEA) was performed between the high- and low-risk groups in TCGA and GSE21257 gene expression matrixes with a threshold of NES >1 and FDR <0.05. As shown in Figures 9A,B, the immune-related pathways including antigen processing and presentation, chemokine signaling pathway, complement and coagulation cascades, cytokine–cytokine receptor interaction, natural killer cell–mediated cytotoxicity, and ribosomes were significantly enriched in the two gene expression matrixes. These analyses connected the ZNF protein gene–based signature with the immune-related pathways. We further compared the enrichment scores of various immune cell subpopulations and related immune functions in the high- and low-risk groups by performing ssGSEA. As shown in Figures 10A,B, the scores of diverse immune cell subpopulations including the CD8⁺ T cells, neutrophils, Th1 cells, and TILs were significantly lower in the high-risk group than in the low-risk group in the two gene expression matrixes. As for the related immune functions, there were significant differences between the high- and low-risk groups for APC co-inhibition, HLA, inflammation promoting, para-inflammation, T-cell co-inhibition, and type I IFN response (Figures 10C,D).

FIGURE 9

FIGURE 9. Gene set enrichment analysis between high- and low-risk groups in gene expression matrixes of the TCGA cohort (A) and GSE21257 cohort (B).

FIGURE 10

FIGURE 10. Comparison of ssGSEA scores between the high- and low-risk groups in TCGA and GSE21257 cohorts. (A, B) Scores of 16 immune cell types in TCGA and GSE21257 cohorts. (C, D) Scores of related immune functions in TCGA and GSE21257 cohorts.

Expression and Kaplan–Meier Survival Analysis of the Seven Zinc Finger Protein Genes

Finally, we performed expression and Kaplan–Meier survival analyses of the seven ZNF protein genes using the data set from the TCGA entire cohort. As shown in Figures 11A–E, the expressions of ATMIN, MKRN3, ZNF71, ZNF438, and ZNF597 were lower in the high-risk group than in the low-risk group. Conversely, the expressions of ZNF525 and ZNF692 were higher in the high-risk group than in the low risk-group (Figures 11F,G). Moreover, the Kaplan–Meier survival analysis suggested that the lower expressions of ATMIN, ZNF438, and ZNF597 and higher expression of ZNF692 were associated with worse overall survival in osteosarcoma, while the expressions of MKRN3, ZNF71, and ZNF525 were not significantly related to the prognosis of patients (Figures 12A–G).

FIGURE 11

FIGURE 11. Comparison of the expression of ATMIN (A), MKRN3 (B), ZNF71 (C), ZNF438 (D), ZNF597 (E), ZNF525 (F), and ZNF692 (G) in high- and low-risk groups.

FIGURE 12

FIGURE 12. Kaplan–Meier survival analysis of ATMIN (A), MKRN3 (B), ZNF71 (C), ZNF438 (D), ZNF597 (E), ZNF525 (F), and ZNF692 (G) in the TCGA cohort.

Discussion

Osteosarcoma is a highly heterogeneous malignant tumor that mainly occurs in the metaphysis of the long bone (Saraf et al., 2018). The clinical outcome of patients with osteosarcoma varies even when patients are treated under the standard management (Rickel et al., 2017). Predicting the prognosis of patients in osteosarcoma accurately is of important clinical value, which may be helpful for risk stratification and clinical decision-making and for guiding individual treatment. Up to now, the most widely used prognostic markers for osteosarcoma include tumor features such as the AJCC-TNM classification and clinical factors, including gender, age, and site (Hattinger et al., 2018). However, the sensitivity and accuracy of these factors in predicting the prognosis of patients with osteosarcoma are limited. Recently, reanalyzing data sets from public databases using bioinformatic methods provides an avenue for the identification of prognostic biomarkers in malignant tumors. For example, Zhang et al. (2021) constructed an autophagy-related long noncoding RNA signature in osteosarcoma using data sets from the TCGA database and confirmed the accuracy of the signature in predicting the clinical outcomes of patients with osteosarcoma.

Given the critical role of ZNF proteins in oncogenesis and progression, here we conducted a comprehensive analysis of the ZNF protein genes in osteosarcoma and identified prognosis-related ZNF protein genes. Furthermore, these genes were utilized to constructed a seven-gene prognostic signature that allowed patients to be stratified into high- and low-risk groups. The Kaplan–Meier survival analyses in the internal and external cohorts confirmed that the patients of the high-risk group had worse overall survival than these of the low-risk group, suggesting the satisfactory accuracy and generalizability of the ZNF protein gene–based signature in predicting the prognosis of the osteosarcoma patients. Moreover, time-dependent ROC curve analysis confirmed the sensitivity and specificity of the ZNF protein gene–based signature. Compared with the previously reported prognostic signatures in osteosarcoma (Zhang et al., 2021), we found that the specificity and accuracy of our ZNF protein gene–based signature were superior to these prognostic signatures, since the AUC values were higher in our research. Moreover, a prognostic nomogram comprising clinical features such as gender and age and the relative risk score based on the risk signature was constructed, and it might be helpful for clinical decision-making and designing personalized management. To identify the underlying mechanism associated with the ZNF protein gene–based signature, genes differentially expressed between the high- and low-risk groups were screened out and subjected to functional enrichment analyses. The results suggested that the DEGs participated in immune-related biological processes such as MHC class II protein complex, immune response–activating signal transduction, humoral immune response, and lymphocyte-mediated immunity. Interestingly, the gene set enrichment analysis also connected immunity with the ZNF protein gene–based signature. We also found that the risk-related DEGs were enriched in COVID-19–related pathways. Subsequently, we further compared enrichment scores of various immune cell subpopulations and related immune functions between the high- and low-risk groups. We found that the scores of diverse immune cell subpopulations, including CD8⁺ T cells, neutrophils, Th1 cells, and TILs, were significantly lower in the high-risk group than in the low-risk group. As for the related immune functions, there were significant differences between the high- and low-risk groups for APC co-inhibition, HLA, inflammation promoting, para-inflammation, T-cell co-inhibition, and type I IFN response. These analyses indicated the immune suppressive states of the high-risk group, and we might speculate that patients in the high-risk group may not respond to immune checkpoint inhibitors (Zhang et al., 2020). Therefore, our ZNF protein gene–based signature could be utilized to predict the outcome of patients treated with immune therapy.

Of the seven genes comprising our prognostic signature, the expressions of ATMIN, ZNF438, and ZNF597 were lower in the high-risk group, and lower expressions of ATMIN, ZNF438, and ZNF597 predicted worse overall survival in patients with osteosarcoma, suggesting that ATMIN, ZNF438, and ZNF597 might act as protective factors in osteosarcoma. Up to now, few researches reported the role of ZNF438 and ZNF597 in tumors. ATMIN, also named as ZNF822, functions as a cofactor of ATM and is required for the activation of the ATM signaling pathway (Kanu et al., 2010; Liu et al., 2017). ATM participates in the maintenance of genomic stability and DNA damage repair (Prochazkova et al., 2015; Li et al., 2019) and has been implicated as a tumor suppressor in various cancer types including lung adenocarcinoma (Foster et al., 2019), B-cell lymphoma (Loizou et al., 2011), and colorectal cancer (Li et al., 2021). However, little is known about the role of ATMIN in osteosarcoma, and based on our analysis, we can speculate that ATMIN functions as a tumor suppressor in osteosarcoma. The effect of ATMIN on osteosarcoma cell proliferation, migration, and invasion and the underlying mechanism need to be further investigated. MKRN3, also known as ZNF127, acts as an E3 ligase and is responsible for the degradation of target proteins through ubiquitination processes (Abreu et al., 2015). MKRN3 is frequently mutated in non–small-cell lung cancers (NSCLCs), and genomic aberrations of MKRN3 are found to be enriched in NSCLC samples harboring KRAS mutations (Li et al., 2021). Lower expression of MKRN3 is associated with worse overall survival in patients with NSCLCs. Reconstitution of MKRN3 in NSCLC cells directly abrogates tumor growth in vitro and in vivo, suggesting MKRN3 functions as a tumor suppressor in NSCLCs. Here, we found MKRN3 expression was lower in the high-risk group, and the lower expression of MKRN3 had a tendency to be associated with a worse overall survival in osteosarcoma. Therefore, MKRN3 might be a protective factor in osteosarcoma. ZNF692, also known as AREBP and Zfp692, was reported to have an important role in cell proliferation, invasion, and metastasis of lung and colon adenocarcinomas (Zhang et al., 2017; Xing et al., 2019). ZNF692 was upregulated in cervical cancer tissues, and its overexpression was associated with poor clinicopathological characteristics in patients with cervical cancer (Zhu et al., 2019). These reports indicate that ZNF692 might be an oncogene in tumors, and consistently, we found that ZNF692 was higher in the high-risk group than in the low risk-group, and the higher expression of ZNF692 was associated with worse overall survival in osteosarcoma. Therefore, it is reasonable to speculate that ZNF692 might have an effect on the malignant behaviors of osteosarcoma cells.

There are some limitations in our study. First, the ZNF protein gene–based signature was only validated using cohorts from the public database, the accuracy and generalizability of the prognostic signature should also be verified in using prospective data, and we are collecting osteosarcoma patients to get our own cohort. Second, the role of the seven ZNF protein genes in osteosarcoma cells should be further explored. We will design specific siRNAs and construct overexpression plasmid for these ZNF protein genes and perform functional experiments to explore the effect of these genes on cell proliferation, migration, invasion, cell cycle, cell apoptosis, and stemness.

In all, we here comprehensively analyzed ZNF protein genes in osteosarcoma and constructed a prognostic signature based on seven ZNF protein genes. The signature is a novel and clinically useful prognostic biomarker for osteosarcoma patients.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, and further inquiries can be directed to the corresponding author.

Author Contributions

WG designed the study. DZ conducted bioinformatic analysis. XS wrote the manuscript and was responsible for language revisions. All authors reviewed the manuscript.

Funding

This study was funded by grants from the Science and Technology Major Project of Hubei Province (2020BCB058).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We sincerely thank the TCGA project for using their data.

References

Abreu, A. P., Macedo, D. B., Brito, V. N., Kaiser, U. B., and Latronico, A. C. (2015). A New Pathway in the Control of the Initiation of Puberty: the MKRN3 Gene. J. Mol. Endocrinol. 54, R131–R139. doi:10.1530/jme-14-0315

CrossRef Full Text | Google Scholar

Bhat, S., Kabekkodu, S. P., Adiga, D., Fernandes, R., Shukla, V., Bhandari, P., et al. (2021). ZNF471 Modulates EMT and Functions as Methylation Regulated Tumor Suppressor with Diagnostic and Prognostic Significance in Cervical Cancer. Cell Biol Toxicol 37, 731–749. doi:10.1007/s10565-021-09582-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cassandri, M., Smirnov, A., Novelli, F., Pitolli, C., Agostini, M., Malewicz, M., et al. (2017). Zinc-finger Proteins in Health and Disease. Cell Death Discov. 3, 17071. doi:10.1038/cddiscovery.2017.71

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Cao, J., Zhang, N., Yang, B., He, Q., Shao, X., et al. (2020). Advances in Differentiation Therapy for Osteosarcoma. Drug Discov. Today 25, 497–504. doi:10.1016/j.drudis.2019.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, H., Ruiz, E. J., Moore, C., Stamp, G. W. H., Nye, E. L., Li, N., et al. (2019). ATMIN Is a Tumor Suppressor Gene in Lung Adenocarcinoma. Cancer Res. 79, 5159–5166. doi:10.1158/0008-5472.Can-19-0647

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, M., and Blackshear, P. J. (2017). RNA-binding Proteins in Immune Regulation: a Focus on CCCH Zinc finger Proteins. Nat. Rev. Immunol. 17, 130–143. doi:10.1038/nri.2016.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Gianferante, D. M., Mirabello, L., and Savage, S. A. (2017). Germline and Somatic Genetics of Osteosarcoma - Connecting Aetiology, Biology and Therapy. Nat. Rev. Endocrinol. 13, 480–491. doi:10.1038/nrendo.2017.16

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, D. J., Geller, D. S., Gill, J. D., Lewis, V. O., and Gorlick, R. (2018). Current and Future Therapeutic Approaches for Osteosarcoma. Expert Rev. Anticancer Ther. 18, 39–50. doi:10.1080/14737140.2018.1413939

PubMed Abstract | CrossRef Full Text | Google Scholar

Hattinger, C. M., Patrizio, M. P., Tavanti, E., Luppi, S., Magagnoli, F., Picci, P., et al. (2018). Genetic Testing for High-Grade Osteosarcoma: a Guide for Future Tailored Treatments? Expert Rev. Mol. Diagn. 18, 947–961. doi:10.1080/14737159.2018.1535903

PubMed Abstract | CrossRef Full Text | Google Scholar

Jen, J., and Wang, Y.-C. (2016). Zinc finger Proteins in Cancer Progression. J. Biomed. Sci. 23, 53. doi:10.1186/s12929-016-0269-9

CrossRef Full Text | Google Scholar

Kager, L., Tamamyan, G., and Bielack, S. (2017). Novel Insights and Therapeutic Interventions for Pediatric Osteosarcoma. Future Oncol. 13, 357–368. doi:10.2217/fon-2016-0261

PubMed Abstract | CrossRef Full Text | Google Scholar

Kansara, M., Teng, M. W., Smyth, M. J., and Thomas, D. M. (2014). Translational Biology of Osteosarcoma. Nat. Rev. Cancer 14, 722–735. doi:10.1038/nrc3838

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanu, N., Penicud, K., Hristova, M., Wong, B., Irvine, E., Plattner, F., et al. (2010). The ATM Cofactor ATMIN Protects against Oxidative Stress and Accumulation of DNA Damage in the Aging Brain. J. Biol. Chem. 285, 38534–38542. doi:10.1074/jbc.M110.145896

CrossRef Full Text | Google Scholar

Laity, J. H., Lee, B. M., and Wright, P. E. (2001). Zinc finger Proteins: New Insights into Structural and Functional Diversity. Curr. Opin. Struct. Biol. 11, 39–46. doi:10.1016/s0959-440x(00)00167-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Leon, O., and Roth, M. (2000). Zinc Fingers: DNA Binding and Protein-Protein Interactions. Biol. Res. 33, 21–30. doi:10.4067/s0716-97602000000100009

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C. G., Mahon, C., Sweeney, N. M., Verschueren, E., Kantamani, V., Li, D., et al. (2019). PPARγ Interaction with UBR5/ATMIN Promotes DNA Repair to Maintain Endothelial Homeostasis. Cell Rep. 26, 1333–1343. e1337. doi:10.1016/j.celrep.2019.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K., Zheng, X., Tang, H., Zang, Y.-S., Zeng, C., Liu, X., et al. (2021a). E3 Ligase MKRN3 Is a Tumor Suppressor Regulating PABPC1 Ubiquitination in Non-small Cell Lung Cancer. J. Exp. Med. 218, 1. doi:10.1084/jem.20210151

CrossRef Full Text | Google Scholar

Li, Y.-J., Yang, C.-N., Kuo, M. Y.-P., Lai, W.-T., Wu, T.-S., and Lin, B.-R. (2021b). ATMIN Suppresses Metastasis by Altering the WNT-Signaling Pathway via PARP1 in MSI-High Colorectal Cancer. Ann. Surg. Oncol. 28, 8544–8554. doi:10.1245/s10434-021-10322-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., King, A., Hoch, N. C., Chang, C., Kelly, G. L., Deans, A. J., et al. (2017). ASCIZ/ATMIN Is Dispensable for ATM Signaling in Response to Replication Stress. DNA Repair 57, 29–34. doi:10.1016/j.dnarep.2017.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Loizou, J. I., Sancho, R., Kanu, N., Bolland, D. J., Yang, F., Rada, C., et al. (2011). ATMIN Is Required for Maintenance of Genomic Stability and Suppression of B Cell Lymphoma. Cancer Cell 19, 587–600. doi:10.1016/j.ccr.2011.03.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Poirier, A., Gagné, A., Laflamme, P., Marcoux, M., Orain, M., Plante, S., et al. (2021). ZNF768 Expression Associates with High Proliferative Clinicopathological Features in Lung Adenocarcinoma. Cancers 13, 4136. doi:10.3390/cancers13164136

PubMed Abstract | CrossRef Full Text | Google Scholar

Prochazkova, J., Sakaguchi, S., Owusu, M., Mazouzi, A., Wiedner, M., Velimezi, G., et al. (2015). DNA Repair Cofactors ATMIN and NBS1 Are Required to Suppress T Cell Activation. Plos Genet. 11, e1005645. doi:10.1371/journal.pgen.1005645

PubMed Abstract | CrossRef Full Text | Google Scholar

Rickel, K., Fang, F., and Tao, J. (2017). Molecular Genetics of Osteosarcoma. Bone 102, 69–79. doi:10.1016/j.bone.2016.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritter, J., and Bielack, S. S. (2010). Osteosarcoma. Ann. Oncol. 21 (Suppl. 7), vii320–vii325. doi:10.1093/annonc/mdq276

PubMed Abstract | CrossRef Full Text | Google Scholar

Saraf, A. J., Fenger, J. M., and Roberts, R. D. (2018). Osteosarcoma: Accelerating Progress Makes for a Hopeful Future. Front. Oncol. 8, 4. doi:10.3389/fonc.2018.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sayles, L. C., Breese, M. R., Koehne, A. L., Leung, S. G., Lee, A. G., Liu, H.-Y., et al. (2019). Genome-Informed Targeted Therapy for Osteosarcoma. Cancer Discov. 9, CD–17. doi:10.1158/2159-8290.Cd-17-1152

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G., and Zheng, C. (2021). Zinc finger Proteins in the Host-Virus Interplay: Multifaceted Functions Based on Their Nucleic Acid-Binding Property. FEMS Microbiol. Rev. 45, 1. doi:10.1093/femsre/fuaa059

CrossRef Full Text | Google Scholar

Xing, Y., Ren, S., Ai, L., Sun, W., Zhao, Z., Jiang, F., et al. (2019). ZNF692 Promotes colon Adenocarcinoma Cell Growth and Metastasis by Activating the PI3K/AKT Pathway. Int. J. Oncol. 54, 1691–1703. doi:10.3892/ijo.2019.4733

CrossRef Full Text | Google Scholar

Ye, Q., Liu, J., and Xie, K. (2019). Zinc finger Proteins and Regulation of the Hallmarks of Cancer. Histol. Histopathol 34, 1097–1109. doi:10.14670/hh-18-121

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Zhang, G., Sun, N., Zhang, Z., Zhang, Z., Luo, Y., et al. (2020). Comprehensive Molecular Analyses of a TNF Family-Based Signature with Regard to Prognosis, Immune Features, and Biomarkers for Immunotherapy in Lung Adenocarcinoma. EBioMedicine 59, 102959. doi:10.1016/j.ebiom.2020.102959

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Ding, R., Wu, T., Jia, J., and Cheng, X. (2021). Autophagy-Related Genes and Long Noncoding RNAs Signatures as Predictive Biomarkers for Osteosarcoma Survival. Front. Cell Dev. Biol. 9, 705291. doi:10.3389/fcell.2021.705291

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Q., Zheng, X., Sun, Q., Shi, R., Wang, J., Zhu, B., et al. (2017). ZNF692 Promotes Proliferation and Cell Mobility in Lung Adenocarcinoma. Biochem. Biophysical Res. Commun. 490, 1189–1196. doi:10.1016/j.bbrc.2017.06.180

CrossRef Full Text | Google Scholar

Zhu, B., Pan, Y., Zheng, X., Zhang, Q., Wu, Y., Luo, J., et al. (2019). A Clinical, Biologic and Mechanistic Analysis of the Role of ZNF692 in Cervical Cancer. Gynecol. Oncol. 152, 396–407. doi:10.1016/j.ygyno.2018.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar