Zinc Finger Proteins in Head and Neck Squamous Cell Carcinomas: ZNF540 May Serve as a Biomarker

Head and neck squamous cell carcinoma (HNSCC) is one of the ten most common cancers. Most cancer cases originate from alcohol and tobacco consumption. However, studies have demonstrated that human papillomavirus (HPV) infection, particularly HPV-16, may also significantly influence disease progression. The KRAB-ZNF family of genes is involved in epigenetic suppression, and its involvement in carcinogenesis is the subject of extensive studies. The available literature data demonstrate that they may play different roles, both as tumor suppressors and oncogenes. In this study, six ZNF genes, ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880, were tested using several in silico approaches based on the TCGA and GEO datasets. Our analyses indicate that the expression of the analyzed ZNFs was significantly downregulated in tumor tissues and depended on tumor localization. The expression levels of ZNFs differed between HPV-positive vs. HPV-negative patients depending on the clinical-pathological parameters. More specifically, the patients with higher levels of ZNF418 and ZNF540 showed better survival rates than those with a lower expression. In addition, the level of ZNF540 expression in HPV-positive (HPV(+)) patients was higher than in HPV-negative (HPV(−)) patients (p < 0.0001) and was associated with better overall survival (OS). In conclusion, we demonstrate that ZNF540 expression highly correlates with HPV infection, which renders ZNF540 a potential biomarker for HNSCC prognosis and treatment.


Introduction
Head and neck squamous cell carcinoma (HNSCC) contributed to the death of 450,000 people worldwide in 2018, which makes it the seventh most severe cancer. HNSCC is mainly associated with tobacco and alcohol abuse. However, human papillomavirus (HPV) infection, mostly with HPV-16, also appears as a crucial etiologic factor in HNSCC development [1,2]. The characteristics of HNSCC are a poor response to treatment and high mortality, where only about 50-60% of patients reach the 5-year survival rate. Thus, there is an urgent need to develop novel, more effective, personalized therapies and specific prognostic biomarkers which are based on genes with protein-coding and non-coding abilities [3][4][5][6]. However, knowledge of the exact molecular mechanisms driving HNSCC is still limited.
Zinc finger proteins (ZNFs) constitute the most numerous family of sequence-specific DNA-binding proteins encoded by 2% of human genes. They bind to their target DNA sequences through the zinc finger domain [7] and exert various functions, including transcriptional regulation, signal transduction, or protection against DNA double-strand breaks [8][9][10]. They may interact with DNA sequences, RNAs, proteins, and post-translational modifications [8][9][10]. ZNFs are divided into several subgroups based on their structural conformation, and C2H2 ZNFs are the most common [7,11]. Besides various zinc finger motifs, the C2H2 class contains additional domains involved in gene expression or cellular localization-such as the KRAB (Krüppel-associated box), SCAN, or BTB/POZ domain [12,13]. Due to the considerable complexity within the ZNF family, little is known about the exact molecular function of most of its members. Of note, many KRAB-ZNFs were shown to play an essential role in carcinogenesis, acting as oncogenes, suppressors, or both, depending on the cancer type [14]. The association with tumor biology was already described for several ZNFs in various cancers, including melanoma [15], colorectal [16,17], renal [18], gastric [19], and esophageal cancers [20], or lung adenocarcinomas [13]. Numerous KRAB-ZNFs show altered expression in various tumors, e.g., HNSCC, as was demonstrated in the transcriptomic profiling based on the TCGA datasets [21]. Nevertheless, the contribution to biological processes and the potential diagnostic utility of specific ZNFs in HNSCC remain undefined. Moreover, there isstill no data on ZNFs' involvement in head and neck cancers with HPV origin.
For this study, based on the preselection with the UALCAN database [22], we chose six ZNF genes: ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880. In our previous analysis, these factors were shown to be downregulated in multiple tumor types, including HNSCC [21]. Moreover, ZNF132 was reported to be epigenetically inactivated in laryngeal squamous cell carcinoma due to promoter hypermethylation [23,24]. In the same tumor type, ZNF418 promoter methylation was demonstrated as a potent diagnostic factor distinguishing between high-and low-risk groups of patients [25]. To the best of our knowledge, no other report has been published to date describing the involvement of ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 in HNSCC. Here, we hypothesize that these genes may be implicated in HNSCC biology and related HPV phenotypes. To test this hypothesis, we used the TCGA data and performed bioinformatics analyses of mRNA expression. We aimed to explore the correlation of ZNF expression with clinico-pathological parameters, their engagement in various cancer-associated processes, and their potential role as biomarkers in HNSCC.

Materials and Methods
In our study, we analyzed six ZNFs preselected based on the UALCAN database: ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880, using the RNA sequencing data downloaded from the TCGA [22]. The study's main steps included: (i) analysis of pathological and clinical features associated with ZNFs, (ii) functional enrichment analysis of genes correlated with selected ZNFs, (iii) analysis of infiltration of immune cells into tumor tissues, and (iv) validation of the selected results. The main steps of the methodology used by us are presented in Figure 1.

Pathological and Clinical Analysis
The differences between healthy and cancer tissues for analyzed ZNFs were obtained from the UALCAN database. To determine whether the expression level of each transcript allowed us to distinguish healthy from cancer samples, we applied a receiver operating characteristic (ROC) analysis with an area under the curve (AUC) estimation in a group of 43 adjacent-matched healthy and neoplastic tissues. We performed the Spearman correlation test to assess the correlation between expression levels of analyzed ZNFs. Next, expression levels of ZNFs were checked depending on localizations in the oral cavity (n = 316), pharynx (n = 90), and larynx (n = 116).
To determine the differences in overall survival (OS) and disease-free interval (DFI), the patients were divided into two groups depending on the expression level of the specified gene using the mean expression level as cut-off: ZNF418 for OS: n Low = 271, n High = 250 and for DFI: n Low = 68, n High = 62; ZNF540 for OS: n Low = 261, n High = 260 and for DFI: n Low = 55, n High = 75. The time of the observation was set as up to 5 years (1825 days).

Functional Enrichment Analysis of Genes Correlated with Selected ZNFs
Genes correlated with the analyzed transcripts were acquired from the cBioPortal for Cancer Genomics (www.cbioportal.org (accessed on 20 November 2020)). Positively (R > 0.3) and negatively (R < −0.3) correlated genes, according to the Spearman correlation, were used to study cellular involvement and interactions with a REACTOME database (https://reactome.org (accessed on 20 November 2020)). To determine statistically significant correlations, p < 0.05 was used for negative correlations with all analyzed transcripts and positive correlations with ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880.
Gene Set Enrichment Analysis (GSEA) software version 4.1 (http://www.gsea-msigdb. org/gsea/index.jsp (accessed on 12 January 2021)) was used to analyze the functional enrichment for HSNCC patients divided into two groups with low and high expression levels of specific ZNFs, using the mean expression levels as a cut-off (the same as for survival analysis). The input file contained expression data from TCGA for 20,530 genes and 520 patients. The groups were compared in terms of Hallmark gene sets (H) and Oncogenic Signatures (C) from the MSigDB collection using an analysis of 1000 gene permutations for testing the significance of the specified gene set enrichment. A nominal p-value of p ≤ 0.05 and a false discovery rate (FDR) ≤ 0.25 were considered significant, as described previously [26].

Infiltration of Immune Cells into Tumor Tissues
Analysis of the Immune and ESTIMATE scores (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) was downloaded from https: //bioinformatics.mdanderson.org/estimate/disease.html (accessed on 12 December 2020). These scores were used to define the infiltration of immune cells into tumor tissues and to infer tumor purity. The subpopulations of specific immune cells were estimated using supporting data presented by Thorsson et al. [27] and analyzed as described previously [26].

Statistical Analysis
All statistical analyses were performed using GraphPad Prism 6, 8, and 9 (GraphPad, San Diego, CA, USA). For two-group analysis, the t-test or Mann-Whitney U test were used for measuring ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 levels and gene expression depending on the Shapiro-Wilk normality test. Next, the expression levels of ZNFs were compared between tumor localizations using the one-way ANOVA and Tukey's multiple comparisons test or Dunn's multiple comparisons test. Receiver operating characteristic (ROC) curve analysis of ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 expression was used to compare adjacent normal and cancerous tissues obtained for 43 patients, and the AUC (Area Under Curve) with a 95% Confidence Interval (CI) was calculated.
For OS and DFI prognosis, the Log-Rank (Mantel-Cox) and Gehan-Breslow-Wilcoxon tests were used, and a 95% CI ratio was calculated. A heatmap was generated using MORPHEUS, an online visualization tool (https://software.broadinstitute.org/Morpheus (accessed on 20 January 2021)). For all of the analyses, p < 0.05 was indicated as statistically significant.

Validation of the Results
To validate the obtained results from the TCGA database, we used the Gene Expression Omnibus (GEO) data repository, with GSE65858 [28] set for HNSCC samples. The

Results
The ZNFs are downregulated in HNSCC and show a high potential to distinguish normal from cancer tissues.
Based on the UALCAN database, we observed a significant downregulation of ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 expression levels in primary tumor tissues compared to normal tissues ( Figure 2A). Moreover, the data indicate a low number of positive correlations between the expression of: ZNF426 and ZNF880 (R = 0.11, p = 0.0093), as well as ZNF426 and ZNF132 (R = 0.10, p = 0.025) genes. For the rest of the analyzed ZNFs, no significant (p > 0.05) correlations were indicated ( Figure 2B).
Next, we applied the ROC curve test to assess the potential of the analyzed ZNFs to discriminate between cancer and healthy tissues. To this end, we utilized paired normal and cancer tissues obtained from 43 HNSCC patients. The data indicate highly sensitive and specific discriminatory abilities for all six ZNFs, with the AUC ranging between 0.77 and 0.91 (p < 0.0001) ( Figure 2C). , and ZNF880 depending on the cancer localization in the oral cavity (n = 316), pharynx (n = 90), and larynx (n = 116) in HNSCC patients. One-way ANOVA was used to assess the statistical significance. Graphs with box and whiskers present 5-95 percentile; CI-confidence interval; ns-not significant; * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001; **** p ≤ 0.0001 are considered as significant.

Expression of the ZNFs Depends on the Tumor Localization and Clinical-Pathological Parameters
Next, the ZNFs expression was analyzed in the HNSCC tissues obtained from various tumor localizations, including the oral cavity, larynx, and pharynx. Significant differences were observed in the expression level of ZFP28, ZNF132, ZNF418, ZNF426, and ZNF540 between the oral cavity and pharynx localizations. The expressions of these factors were downregulated in the oral cavity, except for ZNF426, whose expression increased in the same localization. Moreover, we observed a significant difference in the expression of ZFP28, ZNF132, ZNF426, and ZNF540 between the pharynx vs. larynx localizations, and ZNF132, ZNF418, and ZNF540 in the oral cavity compared to the larynx. The expression levels of ZFP28 and ZNF426 were at the same level for the larynx localization and oral cavity. Moreover, the expression levels of ZNF418, ZNF540, and ZNF132 were more similar between tumors located in the larynx and pharynx than the oral cavity. No significant differences were observed in ZNF880 expression among the three cancer localizations ( Figure 2D).
The expression levels of ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 were also investigated in relation to the clinical-pathological parameters ( Table 1). The analysis demonstrated that the expression of ZFP28, ZNF132, ZNF418, and ZNF540 is significantly lower in patients over 60 years of age. Moreover, a lower transcription of ZNF132, ZNF426, and ZNF540 was noticed in women compared to men. The two important carcinogenic factors, namely, alcohol consumption and smoking, were positively associated with higher ZFP28 expression (p = 0.0480) and the elevated expression of ZFP28, ZNF132, ZNF418, and ZNF880, respectively.
Next, we correlated the expression of selected ZNFs with the TNM classification. In patients with higher T stages (3 and 4), the expression level of ZFP28, ZNF132, and ZNF540 weresignificantly lower than in patients with a less advanced T stage. In the patients with stages N2 and N3, we observed an increased expression of ZNF426 (p = 0.0012) and a decreased level of ZNF132 (p = 0.0064) compared to the patients with N0 and N1. Moreover, in the group of patients with lymph node neck dissection, we noticed lower expressions of ZFP28, ZNF132, and ZNF540. Furthermore, in patients with grade 1 and 2 compared to the group with grade 3 and 4, the expression levels of ZFP28, ZNF132, ZNF418, and ZNF540 were reduced, while that of ZNF426 was increased. We also observed that patients with perineural invasion had a significantly decreased level of ZNF540, whereas patients with angiolymphatic invasion had lower levels of ZNF426 and higher levels of ZNF418 and ZNF540.
We associated the HPV p16 status with the expression levels of ZNFs. We found that the lower level of ZNF426 and the increased levels of ZNF132 and ZNF540 were characteristic for HPV(+) HNSCC patients. All results are presented in Table 1.
Finally, based on hierarchical clustering for the expression levels of ZNFs depending on clinical-pathological parameters, we observed that the direction of the changes in the expression levels of ZNF540 and ZNF132 was very similar in comparison to other ZNFs. Moreover, the expression levels of ZNF426 were distinct compared to the rest of the analyzed genes for all clinical-pathological parameters (Supplementary Figure S1). Table 1. Expression levels of ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 depending on clinical-pathological parameters. t-test or Mann-Whitney U test was used to assess the statistical significance; n-number of cases. The differences with p < 0.05 were considered as significant and marked in bold in the specified cell.

Patients with Low ZNF418 and ZNF540 Expression Display Shorter Overall Survival
We further tested whether the expression of selected ZNFs may correlate with HNSCC patient outcome. To this end, we divided the patient cohort into two groups (high and low expression of each ZNF), with the mean expression used as a cut-off. We focused on the disease-free interval (DFI) and overall survival (OS). In the case of ZFP28, ZNF132, ZNF426, and ZNF880, no differences for OS or DFI were observed (p > 0.05). Moreover, there were no significant differences in DFI for ZNF418 and ZNF540. However, the OS of patients significantly differed. A low expression of ZNF418 and ZNF540 was associated with worse OS compared to the increased expression. The results are presented in Figure 3 and Supplementary Figure S2.

Expression of the ZNFs Is Connected with Critical/Essential Cellular Processes and Pathways
Next, the genes with negative and positive correlations identified via the cBioPortal domain (Spearman correlation: R < −0.3 and R > 0.3) were analyzed using the REAC-TOME online tool. We concentrated mainly on the genes associated with ZNF418 and ZNF540, as both factors showed a linkage to the patients' survival ( Figure 3). The correlated genes were classified into cellular processes and pathways. The studied ZNFs showed positive and negative correlations with various processes. The pathway that was negatively correlated with ZNF418 and ZNF540 included the formation of the cornified envelope, the keratinization process, gap junction trafficking, transport of connexons to the plasma membrane, microtubule-dependent trafficking of connexons from Golgi to the plasma membrane, and prefoldin-mediated transfer of substrate to CCT/TriC. The genes negatively associated only with ZNF418 were involved in gap junction assembly, as well as trafficking and regulation, recruitment of NuMA to mitotic centrosomes, and carboxyterminal post-translational modifications of tubulin. The genes negatively correlated exclusively with ZNF540 were connected with insulin-like growth factor-2 mRNA-binding proteins (IGF2BPs/IMPs/VICKZs) binding RNA, RHO GTPases activate IQGAPs, type I hemidesmosome assembly and signaling by MAPK mutants. For positively correlated genes, only those associated with voltage-gated potassium channels were common for these two ZNFs. Moreover, in this group, for ZNF418, the processes involved in elastic fiber formation, defective B3GALTL causes Peters-plus syndrome (PpS), O-glycosylation of TSR domain-containing proteins, molecules associated with elastic fibers, collagen chain trimerization, cGMP effects, defective CHST6 causes MCDC1, defective B4GALT1 causes B4GALT1-CDG (CDG-2d), and defective ST3GAL3 causes MCT12 and EIEE15 were indicated.
The most significant number of genes (n = 53) was correlated with the pathway involved with ZNF426, namely, the regulation of the expression of SLITs and ROBOs. Moreover, several pathways associated with the studied ZNFs were connected to cancerogenesis, such as defective base excision repair associated with OGG1, which positively correlated with ZNF540; RAS signaling downstream of NF1 loss-of-function variants positively correlated with ZNF426 and negatively with ZNF880 and more. All data are presented in Figure 4A and in Supplementary Table S2.
Next, the GSEA (Gene Set Enrichment Analysis) was carried out and patients with a high and low expression of ZNF418, as well as ZNF540, were compared to show potential differences in processes and pathways. In the group of patients with a low expression of ZNF418, most upregulated genes are MYC targets (V1 and V2), genes connected with KRAS, genes downregulated in primary keratinocytes with knockdown of RB1 and RBL1 genes. Surprisingly, it was observed that, in the case of the patients with a high expression of ZNF418, 121 different processes and pathways were indicated. The most downregulated genes were connected with KRAS changes and KRAS signaling, coagulation, IL2/STAT5 signaling, complement, genes associated with knockdown of PTEN, BRCA1, and RBBP8, and genes changed in cells after IL2 starvation and then stimulated by IL15 or IL21 ( Figure 4B and Supplementary Table S3). In the case of patients with lower levels of ZNF540, it was indicated that the upregulation of genes was observed in foreskin fibroblasts in early response to serum starvation (CSR_EARLY_UP.V1_UP) and genes defining the KRAS dependency signature. For patients displaying higher levels of ZNF540, 48 different processes and pathways were indicated as changed, including genes connected with KRAS, genes changed in cells after knockdown of SUZ12, RBBP8, CRX, as well as genes changed after starvation and later stimulation by IL2, or genes changed in cells after treatment with mTOR pathway inhibitor or dichloroacetate ( Figure 4B and Supplementary Table S3). For ZFP28, ZNF132, ZNF426, and ZNF880, no associations with DFI and OS were observed, but some processes indicated in the GSEA displayed similarities with ZNF418 and ZNF540. It was observed that the SING KRAS DEPENDENCY SIGNATURE process was changed for patients with lower levels of ZFP28 and ZNF880 similarly as for ZNF418 and ZNF540, and in the case of HALLMARK_MYC_TARGETS_V1 for ZFP28, ZNF418, and ZNF880, as well as RB_P107_DN.V1_DN for ZFP28 and ZNF418. Moreover, a similarity was observed for patients with a higher expression of ZFP28, ZNF880, ZNF418, and ZNF540 for over 30 different processes and pathways connected with changes in GSEA's gene sets. ZFP28, ZNF132, ZNF880, ZNF418, and ZNF540 are connected with changes in IL2_UP.V1_DN. All data are summarized in Supplementary Table S3.

The Expression Levels of ZNFs Are Associated with the Tumor Immunological Profile in HNSCC Patients
We further investigated the immunological profile of HNSCC tumors, depending on the low and high levels of ZNF418 and ZNF540 genes. First, the score of stromal cells, immune cells, and finally the ESTIMATE score were evaluated using the ESTIMATE analysis. An elevated fraction of stromal and immune cells were found in the samples with an upregulated expression of ZNF418 and ZNF540 ( Figure 5A). These factors clearly showed that the ESTIMATE scores of the HNSCC samples differ significantly (p < 0.05) depending on ZNF expression (except for ZNF426) ( Figure 5A and Supplementary Figure S3A).
It was indicated that patients with a higher expression of ZNF418 and ZNF540 genes displayed similar immunological profiles, which were manifested by significantly higher levels of lymphocytes and lower levels of mast cells and dendritic cells. Only for ZNF540 were differences in the fraction of macrophages observed, and patients with higher levels of ZNF540 displayed lower levels of infiltration of these cells in the tumor mass ( Figure 5B). Eosinophils and neutrophils did not show any significant changes (p > 0.05) depending on the expression of these two ZNFs ( Figure 5B).
Further analysis of specific subpopulations of immune cells indicated that patients with higher levels of ZNF418 had a higher fraction of CD4+ memory resisting and regulatory Treg cells, and a lower fraction of CD4 naïve T cells. Moreover, higher levels of CD8+, follicular helper cells, and regulatory Treg cells, as well as lower levels of CD4 naïve T cells for patients with higher levels of ZNF540, were observed ( Figure 5C). In the case of B cells, for patients with higher levels of ZNF418, only significantly higher levels of naïve B cells were observed. In contrast to that, higher levels of both naïve and memory B cell subpopulations were indicated in patients with higher levels of ZNF540( Figure 5C). The last analyzed subtypes were macrophages. In the case of higher expression levels of ZNF418, a significantly increased population of M2 macrophages was observed. In contrast to that, the reduced (p < 0.05) fraction of M0 and M2 populations was characteristic for the tumors with higher expressions of ZNF540 ( Figure 5C).
As mentioned above, for ZFP28, ZNF132, ZNF426, and ZNF880, no associations with DFI and OS were observed; however, immunological profiling depending on the expression levels of these ZNFs was carried out. Patients with higher levels of ZFP28, ZNF132, and ZNF880 displayed similar immune profiles with a higher fraction of lymphocytes and lower levels of macrophages and dendritic cells. Only for ZNF426 were significant differences (p < 0.05) indicated for mast cells, dendritic cells, as well as neutrophils (Supplementary Figure S3C). Moreover, the analysis of specific subpopulations of T cells showed that patients with higher levels of ZNFs displayed higher levels of resistant CD4+ memory cells. It should be noted that changes in the expression levels of ZNF426 were associated with differences in most subtypes of T cells. A higher fraction of naïve B cells was observed for ZFP28, ZNF132, and ZNF880, and only patients with higher levels of ZNF426 displayed a lower fraction of memory B cells. Surprisingly, macrophage profiles differed the most in the analyzed immunological cells depending on ZNF levels, and significant differences were observed only in the case of M1 and M0 subtypes for ZFP28 and ZNF132, respectively. All results are presented in Supplementary Figure S3C.

Validation of ZNF540 as a Potential Biomarker Using GEO Data
Next, we utilized the GSE65858 dataset to validate the possible association of ZNF540 expression with HPV status that was identified in our previous analysis based on the TCGA data.
The expression level of ZNF540 was also assessed in four different molecular cancer clusters. The highest expression of ZNF540 was observed in the "atypical (IR) 1" cluster compared to the "basal 4", "classical 2", and "mesenchymal 3" types of HNSCCs ( Figure 6A). Moreover, we tested the association between the levels of ZNF540 expression and various clinicopathological parameters and we found some significant differences only in the case of smoking. Interestingly, no differences in ZNF540 expression were observed between various tumor localizations (p > 0.05). All data are presented in Supplementary  Table S4.
First of all, our results from the GSE65858 dataset confirmed the upregulation of ZNF540 expression in the case of HPV(+) compared to HPV(−) HNSCC samples (6.51 ± 0.021 vs. 6.43 ± 0.009; p = 0.0002) ( Figure 6B). The ROC analysis also indicated the high ability of ZNF540 to distinguish HPV-positive and negative patients (AUC = 0.65; p = 0.0002) ( Figure 6C). Moreover, we observed that the patients infected with HPV type 16 had a higher expression of ZNF540 in comparison to other types of HPV (6.55 ± 0.023 vs. 6.36 ± 0.015; p = 0.0001) ( Figure 6D). Next, the expression levels of ZNF540 were examined depending on the virus activity. Its upregulation was demonstrated in the group of HPV(+) (DNA+/RNA+) vs. HPV(+) (DNA+/RNA−) (6.59 ± 0.027 vs. 6.48 ± 0.046; p = 0.008). The ROC analysis and estimation of the AUC revealed that the ZNF540 expression level may be utilized to discriminate between an active and inactive HPV infection status with high specificity and sensitivity (AUC = 0.72; 95% CI = 0.56 to 0.87; p = 0.0084) ( Figure 6E,F).
Patients from the GSE65858 dataset were divided based on the HPV status, and a between-group OS was calculated. The results suggest a slightly better outcome of the HPV(+) patients compared to the HPV(−) patients (p = 0.0552). For the verification of the possible usage of ZNF540 expression level as a prognostic marker, we divided HNSCC patients into the following subgroups: all patients (HPV(+) and HPV(−)), only HPV(−), and only HPV(+). Each of the subgroups was further divided into high and low expression groups based on the mean ZNF540 expression. For all patients and for HPV(−) patients, no differences in overall survival were observed between groups with a low and high expression of ZNF540 (p = 0.7060 and p = 0.6805, respectively). However, a significantly better OS was observed for the HPV(+) patients with a higher level of ZNF540 compared to the low expression group, with a median of survival of 1249 vs. 933 days, respectively, (95% CI = 0.3093 to 0.9909; p = 0.0351) ( Figure 6G). Mann-Whitney U test or one-way ANOVA test with post-test; nnumber of cases, CI-confidence interval; ns-not significant, ** p < 0.01, *** p < 0.001, **** p < 0.0001; p a -log-rank (Mantel-Cox) test, p b -Gehan-Breslow-Wilcoxon Test; p < 0.05 considered as significant.

Discussion
The zinc finger proteins (ZNFs) are one of the most abundant proteins encoded in the human genome. However, due to the vast complexity of this large family of transcriptional factors, the exact roles of ZNFs are still unexplored. In this study, we analyzed ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 in HNSCC, focusing on their biological role, association with various clinico-pathological parameters, and potential utility as biomarkers. The analysis was carried out using the TCGA data, followed by the validation with an alternative dataset from GEO.
First of all, our results based on the single-gene approach show that the expression of all analyzed ZNFs was lower in HNSCC samples compared to healthy controls, which is a favorable feature for diagnostic biomarkers. This approach confirmed the outcomes of our previous pan-cancer transcriptomic analysis utilizing the TCGA data [21]. Surprisingly, we found that the cross-correlation between the expression of ZNFs was only marginal. Of note, our ROC curve test showed a very good capacity for each ZNF to discriminate between tumor and normal samples, further pinpointing their potential applicability as biomarkers.
Although none of the existing studies explored the above transcripts in more detail in HNSCC, our data align with other published observations. For example, ZNF132 was shown to be epigenetically silenced via promoter hypermethylation in HNSCC [23,24], esophageal squamous cell carcinoma (ESCC) [20], and lung adenocarcinoma (LUAD) [29]. Phenotypically, the cells with reduced ZNF132 expression had decreased mobility in LUAD [29] and growth, migration, invasion, and tumorigenicity in ESCC [20]. These observations suggest the tumor suppressor function of ZNF132. High promoter methylation was also reported in the case of ZNF418 in laryngeal squamous cell carcinoma [25] and ZNF540 in clear cell renal cell carcinoma [18]. Moreover, the study by Hui et al. indicated that ZNF418 was significantly downregulated in gastric carcinoma patients [19]. In contrast, ZNF880 and ZFP28 were found upregulated in colorectal cancer [16] and melanoma [15], respectively. Interestingly, promoter hypermethylation frequently leads to the epigenetic inactivation of ZNF genes with the TSG features in various cancer types [14]. Thus, it is likely that at least some of the above-analyzed ZNFs may also become downregulated via the CpG methylation mechanism in HNSCC, and this possibility warrants further studies.
We further investigated whether the mRNA expression of selected ZNFs may differ depending on various clinico-pathological parameters. We observed particular similarities between ZFP28, ZNF540, and ZNF132 signatures. First, we found that the expression level of most ZNFs (apart from ZNF880) depended on tumor location. The pharynx was the site of the highest expression for ZFP28, ZNF540, and ZNF132, and the lowest for ZNF426. Of note, ZNF132 and ZNF540 expression differed in all three anatomical sites: oral cavity, larynx, and pharynx. Secondly, ZFP28, ZNF540, and ZNF132 were downregulated in tumors with a higher T stage, in older patients, and in the cohort that underwent lymph node dissection from the neck. In contrast, these factors demonstrated an increased expression in high-grade tumors (G3 + G4). For further explanation of whether these ZNFs were associated with more aggressive forms of HNSCC, we analyzed their expression profile in different molecular subtypes. We observed that all examined ZNFs were upregulated in the mesenchymal and less aggressive, atypical [30] tumors.
According to our knowledge, there are no comprehensive studies on the expression of ZNFs in the context of HNSCC and related risk factors. Here, we demonstrate that a higher ZFP28 expression was associated with alcohol consumption, whereas smoking was related to higher ZFP28, ZNF880, and ZNF418 levels and lower ZNF132 expression. Moreover, we observed that ZNF132 and ZNF540 were upregulated, and ZNF426 was downregulated in the HPV(+) group compared to the HPV(−) group. Although ZNF132 expression and promoter methylation were analyzed previously in HPV(+) HNSCC cases, no association with HPV was demonstrated [24]. Such a discrepancy between our study and [24] may reflect ethnic differences between populations analyzed or may be due to the lower number of patients included in [24].
Since their biological roles remain largely uncharacterized, we further sought to determine ZNFs' involvement in tumor-associated pathways. Using REACTOME and GSEA tools, we confirmed the correlation of ZNFs with various signaling pathways engaged in tumorigenesis and immune responses. In summary, those pathways include MAPK, NF-κB, TNF, JNK, and RAS signaling. For example, ZNF418 and ZNF540 expression was linked to KRAS signaling. In addition, ZNF540 expression positively correlated with defective base excision repair associated with OGG1 and with the expression of the TNF receptor superfamily (TNFSF) that mediates non-canonical NF-κB signaling, which is essential for immune response and cell growth regulation [31]. Our data also indicate that both ZNF540 and ZNF418 are associated with IL-2 signaling (and IL-15 and IL-21 in the case of ZNF418). Moreover, we found that both factors were related to altered immunological profiles in HN-SCC patients. Thus, it may be hypothesized that the decreased expression of ZNF540 and ZNF418 may affect tumor formation not only through various oncogene-related pathways but also via interfering with the immune response.
Furthermore, our study, for the first time, demonstrates that ZNF418 and ZNF540 expressions could be used as potential biomarkers in HNSCC. Notably, the patients with increased levels of ZNF418 showed significantly longer OS, while those with higher ZNF540 expression had prolonged disease-free intervals. Other reports may indirectly support our findings. For example, high ZNF418 expression correlated with improved OS in gastric carcinoma [18], whereas its promoter hypermethylation showed a good discriminatory potential between high-and low-risk patient cohorts [25]. Additionally, Arai et al. determined that ZNF540 is frequently methylated in clear renal cell carcinoma patients with worse survival [18]. Of note, our data reveal a link between ZNF540 expression and NF-κB, MAPK, and JNK pathways, which contribute to the epithelial-to-mesenchymal transition (EMT) [32,33]. Hypothetically, ZNF540 may interfere with the EMT events, thus lowering the chance of local and distant metastases and improving the prognosis of HNSCC patients. Nevertheless, these hypotheses require further clinical and wet-lab investigation.
Importantly, our detailed analysis of the potential clinical usage of ZNFs revealed that ZNF540 expression might serve as a prognostic marker in the context of HPV infection. Based on the TCGA and GEO data, we showed that the ZNF540 level is higher in HPV(+) patients than in HPV(−) patients. Moreover, we indicated that a higher expression of ZNF540 is observed in patients with an active HPV infection. So far, only one study revealed that ZNF540 is upregulated in HPV(+) active vs. HPV(+) inactive patients, as well as HPV(+) active and HPV(+) inactive in comparison to HPV(−) HNSCC patients [34]. However, no in vitro studies describe the biological role of ZNF540. As mentioned previously, the highest ZNF540 expression was observed in the HNSCC molecular subtype characterized as "atypical". It was indicated that the atypical subtype was a less aggressive type of HNSCC and was associated with a strong immune signature [35]. Based on our immunological profile results, we observed that patients with higher levels of ZNF540 displayed higher levels of CD8, follicular helper T cells, memory and naive B cells, and lower levels of M2 macrophages. Cillo et al. described the immune landscape of viral-and carcinogen-driven HNSCC and indicated that a higher level of follicular helper T cells was associated with longer progression-free survival [35]. We also observed in the TCGA data that patients with higher levels of ZNF540 displayed longer disease-free intervals. While the TCGA dataset comprised all patients, HPV(+) and HPV(−), further analysis based on the GEO dataset clearly showed that the HPV(+) patients with higher levels of ZNF540 had significantly longer overall survival.
In conclusion, we showed that ZFP28, ZNF132, ZNF418, ZNF426, ZNF540, and ZNF880 had reduced expression in HNSCC compared to healthy tissues. Moreover, their expression levels were associated with various clinical parameters, risk factors, and signaling pathways crucial for tumorigenesis and immune responses. We revealed that the high expression of ZFP540 and ZFP418 correlated with a favorable prognosis in HNSCC. Specifically, high ZFP540 levels were associated with improved survival of HPV(+) patients. Altogether, our findings emphasize the potential applicability of ZNF418 and ZNF540 as prognostic biomarkers in HNSCC. These promising data open new avenues for additional research to dissect the mechanisms responsible for ZNF downregulation. The limitation of our study is that it is based on the TCGA and GEO data, where we had no control over the quality of samples and their sequencing. However, in both data sets, different methodologies were implicated, and the results are similar, which confirms that ZNF540 is closely associated with HPV infection. More importantly, however, the molecular mechanisms contributing to the ZNF540 involvement in HNSCC biology are unknown and need to be clarified in the in vitro cell line models and in vivo based on large patient samples with known HPV status.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/curroncol29120779/s1, Figure S1. Heat map and clustering of mean expression levels of ZNFs depending on the specific clinical parameters; Figure S2. (A) Diseasefree interval (DFI) and (B) overall survival (OS) of HNSCC patients depending on the ZFP28, ZNF880, ZNF426, and ZNF132 expression levels, respectively. Results presented for 5 years of observation, low and high subgroups of patients divided on the mean of expression; p a -Long-rank (Mantel-Cox) test; p b -Gehan-Breslow-Wilcoxon test; p < 0.05 considered significant; Figure S3. The immunological profile of HNSCC patients depends on the low and high levels of ZFP28, ZNF132, ZNF426, and ZNF880 transcripts. (A) Stromal, immune, and ESTIMATE scores; (B) Infiltration of specific immune cells in tumor samples; (C) Differences in the fraction of T cells, B cells, and macrophages; ns-not significant; * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001; **** p ≤ 0.0001 considered as significant; Table S1. Number of patients' cases analyzed in the groups depending on the specific clinical parameters; n-number of cases; Table S2. Collected genes assigned to pathways positively and negatively correlated with ZFP28, ZNF880, ZNF540, ZNF418, ZNF426, and ZNF132; Table S3. Involvement of ZNFs transcripts in cellular processes based on GSEA analysis in HNSCC patients. Normalized enrichment scores for GSEA analysis of MSigDB gene sets for oncogenic and hallmark gene sets in the group of patients with low and high expression levels of specified ZNFs. Only results set with p ≤ 0.05 and FDR ≤ 0.25 were shown. NES (normalized enrichment score), p-Val (nominal p-value), and FDR q-val (false discovery rate); Table S4. Expression level of ZNF540 depending on clinical-pathological parameters in HNSCC patients from GSE65858; Mann-Whitney U test or one-way ANOVA, p < 0.05 considered as significant.

Institutional Review Board Statement:
The study is based on the analysis of freely available data sets and does not need any ethics committee's agreement and does not violate the rights of other persons or institutions. Informed Consent Statement: Not applicable. Study was based on fully available data from the TCGA and GEO databases. Patients' consents were not necessary for uses of the published data.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Raw data are available in the TCGA and GEO databases.