Cuproptosis-associated CDKN2A is targeted by plicamycin to regulate the microenvironment in patients with head and neck squamous cell carcinoma

Head and neck squamous cell carcinoma (HNSCC), the most common malignancy of the head and neck, has an overall 5-year survival rate of <50%. Genes associated with cuproptosis, a newly identified copper-dependent form of cell death, are aberrantly expressed in various tumours. However, their role in HNSCC remains unknown. In this study, bioinformatic analysis revealed that the cuproptosis-related gene CDKN2A was correlated with the malignant behaviour of HNSCC. Kaplan-Meier (KM) curves showed that patients with high CDKN2A expression had a better prognosis. Multiomic analysis revealed that CDKN2A may be associated with cell cycle and immune cell infiltration in the tumour microenvironment and is important for maintaining systemic homeostasis in the body. Furthermore, molecular docking and molecular dynamics simulations suggested strong binding between plicamycin and CDKN2A. And plicamycin inhibits the progression of HNSCC in cellular assays. In conclusion, this study elucidated a potential mechanism of action of the cuproptosis-associated gene CDKN2A in HNSCC and revealed that plicamycin targets CDKN2A to improve the prognosis of patients.

Background

Head and neck squamous cell carcinoma (HNSCC) is the most common malignant tumour of the head and neck that develops in the mucosal epithelium of the mouth, pharynx, larynx, nasal cavity, and sinus cavities (Taberna et al., 2017). It is the sixth most common malignancy worldwide, with approximately 890,000 new cases and 450,000 deaths owing to HNSCC reported worldwide in 2018 (Bray et al., 2018; Ferlay et al., 2019). Although significant progress has been made in the combined use of surgery, radiotherapy, chemotherapy and targeted therapy for the treatment of HNSCC, 40%–50% of patients have post-treatment relapse, and the overall 5-year survival rate is <50% (Canning et al., 2019). Therefore, biomarkers that can improve HNSCC prognosis need to be identified.

Tsvetkov et al., (2022) recently described a new copper-dependent type of cell death, cuproptosis, closely related to mitochondrial respiration. Copper binds to fatty acylated proteins in the tricarboxylic acid cycle during mitochondrial respiration, resulting in the aggregation of fatty acylated proteins and low expression of iron–sulfur cluster proteins, thereby inducing proteotoxic stress and eventually leading to cell death (Tsvetkov et al., 2022). There have been four genes identified for cuproptosis that are positively regulated (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB and ATP7B) and eight genes that are negatively regulated (MTF1, GLS, cyclin-dependent kinase inhibitor 2A [CDKN2A] and SLC31A1) (Tsvetkov et al., 2022). Biological processes such as tumour cell proliferation, vascular growth, and metastasis have been shown to be significantly affected by copper, an essential element for mitochondrial respiration and iron uptake (Ruiz et al., 2021; Oliveri, 2022). There are a number of factors that regulate and maintain the body’s intake, excretion, and metabolism of copper. Copper-induced cell death or abnormal copper metabolism can occur when copper homeostasis is disrupted in the body. Previous research reported that the cuproptosis-related gene LIPT1 is significantly correlated with prognosis and immune infiltration in melanoma (Lv et al., 2022). Similar findings have been reported by Yun Y et al. indicating that SLC31A1, DX1 and TP7B0 are associated with lung cancer (Yun et al., 2022). Several studies have shown that cuproptosis-related genes are associated with poor prognosis, reduced drug sensitivity, and tumor microenvironment in renal clear cell carcinoma (Ji et al., 2022). However, further studies are needed to determine the effects of cuproptosis on HSNCC.

Bioinformatic tools have become increasingly popular with the rapid development of high-throughput technology for assessing prognostic markers and studying mechanisms (Wei et al., 2019; Shi et al., 2021; Luo et al., 2022a; Luo et al., 2022b; Chen et al., 2022; Lin et al., 2022; Mei et al., 2022; Xuan et al., 2022; Zhao and Jiang, 2022). In melanoma, lung cancer, and renal clear cell carcinoma, cuproptosis-related genes have been studied. In recent studies, researchers examined ten cuproptosis-related lncRNAs associated with immune function and prognosis in HNSCC (Li et al., 2022). However, the functions and mechanisms of action of cuproptosis-related genes in HNSCC warrant further investigation. Molecular dynamics simulation serves as one of the important tools for assessing the stability of drug-targeting ligands for the development of oncological therapeutics (Chikan and Vipperla, 2015; Thai et al., 2015). The aim of this study was to identify cuproptosis-related genes associated with poor prognosis of HNSCC and screen for potential candidate drugs for its treatment using molecular docking and molecular dynamics simulations.

Materials and methods

Gene identification and data acquisition

TCGA-HNSCC data were analyzed to identify differentially expressed genes (DEGs) between tumours and healthy tissues as previous researches (Tomczak et al., 2015; Kang et al., 2021). Based on the intersection of tumour-associated DEGs and cuproptosis-associated genes, CDKN2A was identified as a key gene. We investigated the correlation between CDKN2A expression and clinical staging of HNSCC. CDKN2A prognostic significance was evaluated using KM curves. In order to investigate the performance of CDKN2A for predicting 1, 3 and 5 year overall survival (OS), ROC curves were plotted using the time ROC package. Furthermore, ROC curves were used to assess the relationship between CDKN2A and HNSCC clinical characteristics.

CDKN2A function assessment

For further analysis, RNA-seq data were extracted from TCGA, and log2-transformed gene expression data were obtained (Tomczak et al., 2015). Based on the median CDKN2A expression, tumour samples were divided into high and low-expression groups for survival analysis. In the high and low CNKN2A expression groups, DEGs were screened using the limma R package. Adjusted p-values of <.05 and |logFC| values of >1 were used as the screening criteria for significant DEGs (Ritchie et al., 2015). The ggplot2 package was used to plot volcano and heat maps to visualise the expression of 17 significant DEGs. Data on CDKN2A mutations was downloaded from TCGAbiolinks and visualized using track Viewer (Colaprico et al., 2016; Ou and Zhu, 2019).

Gene set variation analysis (GSVA)

The “GSVA” package was used to analyze all CDKN2A-associated DEGs, followed by the “limma” package to identify high and low CDKN2A expression levels (Hänzelmann et al., 2013; Lin et al., 2021).

GSEA and KEGG enrichment analysis

GSEA was performed on DEGs associated with CDKN2A, KEGG enrichment analysis was performed with gseKEGG, and pathways of interest were visualized using ggplot2 (Subramanian et al., 2005; Ito and Murphy, 2013).

Immuno-infiltration analysis

In order to calculate 28 immune cells, ssGSEA was used on gene expression profile data (Barbie et al., 2009). Based on CDKN2A expression levels, 28 immune cells were compared between groups with high and low expression levels of the protein. The Pearson correlation coefficient was used to examine the correlation between CDKN2A and immune cell infiltration further. Pearson correlation coefficient was used to analyze the correlation between each type of immune cell.

Acquisition and optimization of FDA structures

Food and Drug Administration (FDA) approved 2,568 small molecules (as of 2022-01-04). We downloaded the small molecule structures from Drug Bank (http://www.drugbank.com/) in SDF format (Luo et al., 2022b). In RDKit, the Experimental-Torsion Basic Knowledge Distance Geometry (ETKDG) algorithm was used to generate 3D conformations based on the modified distance geometry algorithm, while the MMFF94 stand was used to optimize small molecule structure and energy using the MMFFOptimize Molecule module (Li et al., 2022).

The structures of the proteins were obtained from the Uniport website (https://www.uniprot.org/). PDB structure 1A5E obtained under CDKN2A entry (P42771) using gene name query. Docking was carried out using Smina Chikan and Vipperla, (2015). Protein-Ligand Interaction Profiler (Plip, https://plip-tool.biotec.tu-dresden.de/plip-web/plip) was used to analyze docking results.

Analysis of PPI networks and gene networks

By using the STRING website and cytoscape software, we explored CDKN2A’s protein-protein interaction network (Shannon et al., 2003; Szklarczyk et al., 2021). Genemania was used to analyze the CDKN2A network (Warde-Farley et al., 2010).

Construction of ceRNA network

Firstly, the multiMiR package was used to find miRNAs related to CDKN2A (Huang et al., 2019a). The LncRNADisease database was then used to identify HNSCC-related LncRNAs, and the miRTarBase database was used to identify shared miRNAs with HNSCC-related LncRNAs (Bao et al., 2019; Huang et al., 2019b). Cytoscape is used for the final visualization (Shannon et al., 2003).

Drug analysis

To identify drugs that may act on DEGs (99 highly expressed genes and 56 lowly expressed genes) between high and low expression groups of CDKN2A, the cmap website was accessed (http://clue.io/).

Molecular dynamics simulation (MDS)

The lowest energy conformation was selected as the kinetic initial conformation. The quantitative software Orca was used to perform quantum chemical optimization for small molecules under B3LYP/6-31G* basis set conditions, involving corrections for bond lengths, bond angles, dihedral angles, and calculations of RESP2 at 0 fixed charges. Gromacs 2019.6 was chosen as the kinetic simulation software, amber14sb was chosen as the protein force field, and Gaff2 force field was chosen for small molecules, and the TIP3P water model was used to add TIP3P water model to the complex system to build a water box and add sodium ions to equilibrate the system. Under elastic simulation by Verlet and cg algorithms, Particle-mesh Ewald (PME) deals with electrostatic interactions using the steepest descent method for energy minimization for the maximum number of steps (50,000 steps). The Coulomb force cutoff distance and van der Waals radius cutoff distance were both 1.4 nm. Finally the system was equilibrated with the regular system (NVT) and isothermal isobaric system (NPT), and the MD simulations were performed at room temperature and pressure for 100 ns. During the MD simulation, the hydrogen bonds involved were constrained using the LINCS algorithm with an integration step of 2 fs. The PME method was calculated with a cutoff value set to 1.2 nm and a non-bond interaction cutoff value set to 10 Å. The V-rescale temperature coupling method was used to control the simulation temperature at 300 K and the Berendsen method to control the pressure at 1 bar. Additionally, 30 ps of NVT and NPT equilibrium simulations were performed at 300 K. Finally, 50 ns of finished MD simulations were performed for the protein–ligand complex system. Root mean square fluctuations (RMSF) were used to observe the local loci variation structure of the system during the simulation (the fluctuation cutoff was set to 0.2). The radius of gyration (Rg) was used to evaluate the tightness of the system structure. The RMSF can observe the local loci variation of the system during the simulation.

Calculation of free energy of binding of proteins and small molecules

The MD trajectory operation is calculated by the following equation:

${\Delta G}_{bind}={\Delta G}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}–\left({\Delta G}_{\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}\mathrm{o}\mathrm{r}}+{\Delta G}_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{a}\mathrm{n}\mathrm{d}}\right)$

${=\Delta \mathrm{E}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}+{\Delta \mathrm{E}}_{\mathrm{V}\mathrm{D}\mathrm{W}}+{\Delta \mathrm{E}}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}+{\Delta \mathrm{G}}_{\mathrm{G}\mathrm{B}}+{\Delta \mathrm{G}}_{\mathrm{S}\mathrm{A}}$

In the above equation, ${\Delta \mathrm{E}}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}$ represents internal energy, ${\Delta \mathrm{E}}_{\mathrm{V}\mathrm{D}\mathrm{W}}$ represents van der Waals interaction, and ${\Delta \mathrm{E}}_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}}$ represents electrostatic interaction. Internal energy includes bond energy (Ebond), angular energy (Eangle), and torsion energy (Etorsion). ${\Delta \mathrm{G}}_{\mathrm{G}\mathrm{B}}$ and ${\Delta \mathrm{G}}_{\mathrm{S}\mathrm{A}}$ are collectively referred to as the solvation free energy. ${\mathrm{G}}_{\mathrm{G}\mathrm{B}}$ is the polar solvation free energy and ${\mathrm{G}}_{\mathrm{S}\mathrm{A}}$ is the non-polar solvation free energy. We used the GB model for calculating ${\Delta \mathrm{G}}_{\mathrm{G}\mathrm{B}}$ (Weiser et al., 1999). The ${\Delta \mathrm{G}}_{\mathrm{S}\mathrm{A}}$ was calculated based on the product of surface tension (γ) and solvent accessible surface area (SA): ${\Delta \mathrm{G}}_{\mathrm{S}\mathrm{A}}=0.0072\times \Delta \mathrm{S}\mathrm{A}$. The entropy variation was neglected in this study due to the high consumption of computational resources with low accuracy.

Cell culture and wound-healing assay

Human laryngeal cancer cell lines (TU212) were obtained from the Cell Bank of the Chinese Academy of Sciences. For cell culture, Dulbecco’s modified Eagle’s medium (DMEM) containing 10% fetal bovine serum (FBS) was used. In an incubator with 5% CO2 and 37°C, 10%FBS and 1% penicillin all cell lines were cultured. In TU212 cells, plicamycin was applied for 14 days at 10 nmol/L. Incubation at 37°C was carried out for both treated and untreated cells plated on 10-cm culture dishes. With a plastic pipette tip, a lane was scratched through the confluent monolayers, followed by addition of DMEM, 1% FBS, and 10 nmol/L plicamycin. Several wounded areas were observed and then photographed through a microscope 24 h after the scratch.

Statistical analysis

Standard error of the mean is indicated by bars on figures, and was calculated using Microsoft Office Excel 2016. All experiments were performed with at a minimum of triplicate samples, and all p-values were calculated with two-tailed t-tests.

Results

Data acquisition and screening of DEGs

After intersecting tumour-associated DEGs with cuproptosis-related genes, CDKN2A was identified as a key gene (Figure 1A), which may play a critical role in tumour development through cuproptosis. CDKN2A expression was significantly higher among patients with T2-stage HNSCC than among patients with T3-and T4-stage HNSCC (Figure 1B), indicating that CDKN2A expression is high in early malignant tumours and low in advanced malignant tumours. KM curves demonstrated that high CDKN2A expression was associated with a better clinical prognosis (Figure 1C). The area under the ROC curves indicated a more significant protective effect of CDKN2A on 1-, 3- and 5-year OS (Figure 1D). In addition, age; sex; T, N, and M stages and CDKN2A were found to have predictive significance for the prognosis of HNSCC (Figure 1E).

FIGURE 1

FIGURE 1. Identification of key genes. (A) Venn diagram demonstrating the results of intersection analysis of differentially expressed genes and cuproptosis-related genes; (B) Bar graph demonstrating the difference in CDKN2A expression among patients with different T stages; (C) KM curve demonstrating a better prognosis for patients with head and neck squamous cell carcinoma with high CDKN2A expression; (D) ROC curve demonstrating the predictive performance of CDKN2A for 1-, 3- and 5-year OS; (E) ROC curves demonstrating the relationship between clinical characteristics and the predictive performance of CDKN2A for clinical prognosis.

Characterisation of CDKN2A expression and enrichment analysis

HNSCC samples in TCGA cohort were divided into the high- and low-CDKN2A-expression groups based on the median CDKN2A expression. The scatter plot demonstrated differences in the survival status of patients between the two groups (Figures 2A, B). A total of 155 DEGs were identified in the two groups; of which, 99 were upregulated and 56 were downregulated (Figure2D). The top 17 genes with the most significant differences in expression were visualised on a heat map (Figure 2C). GSEA suggested that these genes were significantly enriched in pathways associated with cell cycle and DNA replication in the high-CDKN2A-expression group (Figure 2E). In addition, GSVA suggested that the genes were enriched in pathways associated with DNA repair, E2F targets and the G2M checkpoint in the high-CDKN2A-expression group and pathways associated with coagulation, apical junction and inflammatory response in the low-CDKN2A-expression group (Figure 2F). These results suggest that high CDKN2A expression is associated with activation of the cell cycle in HNSCC.

FIGURE 2

FIGURE 2. Relevant functions of CDKN2A. (A) Tumour samples in TCGA database were divided into the high- and low-expression groups according to median CDKN2A expression; (B) Scatter plot demonstrating the difference in survival status between the two groups; (C) Top 17 differentially expressed genes; (D) Volcano plot demonstrating the distribution of differentially expressed genes; (E) GSEA; (F) Bar graph demonstrating the differentially enriched pathways between the two groups.

Immune infiltration analysis

ssGSEA revealed that the infiltration of memory B cell was higher in the low-CDKN2A-expression group; however, no significant difference was observed in the infiltration of other immune cells between the two groups (Figure 3A and Supplementary Figure S1A). Correlation analysis revealed a co-expression relationship among 28 immune cell types (Supplementary Figure S1B), and the correlation between CDKN2A and immune cells indicated that CDKN2A promoted immune cell infiltration (Figure 3B). CDKN2A was significantly positively correlated with activated B-cell (cor = 0.182, p < 0.001) and activated CD4 T-cell (cor = 0.160, p < 0.001) (Figures 3C, D) and significantly negatively correlated with neutrophils (cor = −0.262, p < 0.001) and gamma-delta T-cell (cor = −0.166, p < 0.001) (Figures 3E, F). This suggests that CDKN2A is involved in tumor immune microenvironment homeostasis.

FIGURE 3

FIGURE 3. Correlation between CDKN2A and immune cells. (A) Differences in immune cell infiltration between the high- and low-CDKN2A-expression groups were examined via ssGSEA; (B) Correlation between CDKN2A and the infiltration of immune cells; (C–F) Correlation plots demonstrating the relationship between the top four immune cell types and CDKN2A.

Gene network and mutation analysis of CDKN2A

A protein–protein interaction (PPI) network is shown in Figure 4A. The gene interaction network of CDKN2A was visualised using geneMANIA (Figure 4B). In the ceRNA network of CDKN2A, HOTAIR was found to be closely associated with has-miR-34a and has-miR-125a (Figures 4C, D). HOTAIR may regulate has-miR-34a and has-miR-125a in HNSCC, thus exerting a regulatory effect on CDKN2A expression (Figure 4D). In addition, CDKN2A had a higher mutation rate at the p. R80 locus than at the p. R58 and p. W110 loci (Supplementary Figure S2).

FIGURE 4

FIGURE 4. Single gene network diagram. (A) PPI network analysis of CDKN2A; (B) Gene interaction analysis of CDKN2A using geneMANIA; (C) miRNAs targeted by CDKN2A were predicted based on data extracted from seven databases; (D) Visualization of the ceRNA network maps using Cytoscape.

Molecular docking results

The distribution of docking fractions for virtual screening is shown in Figure 5. The binding affinity of proteins and small molecules can be obtained after molecular docking. The docking fraction for each ordinal number corresponds to −5.42 kcal/mol, with a mean value of −5.42 kcal/mol. Based on a potential screening threshold of −7 kcal/mol, the screening rate is 12.13% (308/2,539), i.e., 12.13% of small molecules have the potential for precursor optimization as Potential dead drug set (PLDS). The mid-docking effects of PLDS were examined to avoid analytical bias caused by transient false-positive docking results. There were 104 residue sites contacted (Supplementary Figure S3). A cut-off of 20 was selected to exclude transient positive contacts. The results are shown in Figure 6. The orange colour represents hydrogen bonds; the number share is the main interaction force [47.41% (1427/3010)] and the amino acid residue contact sites include 46-ARG, 47-ARG, 54-MET, 84-ASP, 87-ARG, 88-GLU, 105-ASP, 111-GLY, 116-ASP, 131-ARG, 138-ARG, 139-GLY, 142-HIS, 144-ARG and 147-ALA. The blue colour represents hydrophobic interactions, which are the main interaction forces [38.27% (1152/3010)], and the amino acid residue contact sites include 21-ALA, 44-TYR, 51-VAL, 77-THR, 79-THR, 107-ARG, 110-TRP, 112-ARG, 113-LEU, 117-LEU, 121-LEU, 137-THR, 148-ALA, 149-GLU, 151-PRO and 154-ILE. Molecular docking analysis suggests that plicamycin binds CDKN2A most strongly.

FIGURE 5

FIGURE 5. Different amino acid binding modes. (A) The X-axis represents the serial number of FDA-approved drugs, and the Y-axis represents the docking fraction (unit: kcal/mol); (B) Distribution of docking fractions. The horizontal coordinate represents the docking fraction (unit: kcal/mol), and the vertical coordinate represents the distribution probability.

FIGURE 6

FIGURE 6. Different amino acid binding modes. Blue represents hydrophobic interactions (HI), orange represents hydrogen bonds (HdB), green represents π-π stacking (pSp), red represents π-cation interactions (pCI), purple represents halogen bonds (HIB) and brown represents salt bridge (SB).

Molecular dynamics and cellular effects of plicamycin

Molecular dynamics simulation (MDS) is an important method for assessing the stability of complexes in an aqueous solution (Hollingsworth and Dror, 2018; Hildebrand et al., 2019; Huang et al., 2022; Kang et al., 2022; Zhang et al., 2022). Stability of a system can be measured by the atomic root-mean-square deviation (RMSD). In Figure 7A, RMSD values fluctuated between 20 and 30 ns due to transient instability and then remained stable at 30–50 ns. Root-mean-square fluctuations (RMSF) reflect changes in the local sites of the system during MDS. According to Figure 7B, amino acids at positions 1–15, 33–44, 58, 90, 107, and 124–156 fluctuated more than other amino acids. Figure 7C shows that the radius of rotation (Rg) is an important measure of architecture tightness. In Figure 7D, the solvent-accessible surface area (SASA) of the protein decreases steadily over 0–100 ns, indicating favourable binding and progressive protein tightening. Figure 7E shows the potential-energy curve of hydrogen bonded complexes in their steady state. Figures 8A, B show the binding site, docking pose, and overall protein structure.

FIGURE 7

FIGURE 7. Results of molecular dynamics simulation. (A) Root-mean-square deviation (RMSD) of complex MDS; (B) Root-mean-square fluctuation (RMSF) of complex MDS; (C) Variations in the radius of gyration (Rg) of complex MDS; (D) Variations in SASA of the protein in 0–100 ns of complex MDS; (E) Variations in hydrogen bonding in the stable complex.

FIGURE 8

FIGURE 8. Binding site of plicamycin on the CDKN2A protein. (A) Docking pose; (B) Amino acids and docking sites.

We compared the binding free energies of the two solvated molecules in their bound and free states, as well as the binding free energies of various solvated conformations of plicamycin (Figure 9A). Analysis of the variation of the binding free energy with MDS revealed that the total free energy (G_total) was <0, indicating a likely interaction between the CDKN2A protein and plicamycin. The binding energy of both van der Waals and electrostatic interactions was <0, indicating that hydrophobic interactions and electrostatic energy contribute to the binding between CDKN2A and plicamycin. Non-polar interactions favored binding, whereas polar solvation didn’t. Positive values of polar solvation energy (EGB) indicate that non-polar interactions favored binding. MET-52, MET-53, MET-54, and ASP-84 contact residues of the complex have free energies of zero, indicating they are the major binding sites for CDKN2A and plicamycin (Figure 9B). There was a positive free energy difference between ARG-46 and ARG-87, indicating that plicamycin binds poorly to CDKN2A at these sites. In the stable complex, the contact residue MET-52 contributes to the major binding force (Figure 9B). Furthermore, we found that plicamycin treatment reduced the mobility of TU212 cells in scratch experiments (Figure 10A, B). Accordingly, plicamycin inhibited HNSCC progression in cellular assays. Based on molecular dynamics simulations, this study shows that plicamycin targets CDKN2A to improve patient outcomes.

FIGURE 9

FIGURE 9. Energy analysis of molecular dynamics simulations. (A) Binding free energy of protein–ligand complexes. VDWAALS, van der Waals energy; EEL, electrostatic energy; EGB, polar solvation energy; ESURF, non-polar solvation energy; GGAS, total gas phase free energy; GSOLV, total solvation free energy; G_total free energy = GSOLV + GGAS. (B) The relationship between each contact residue and the binding energy.

FIGURE 10

FIGURE 10. Pplicamycin treatment reduced the mobility of TU212 cells. (A) Wound-healing assay for migration activity of untreated (upper panel) and treated (lower panel) TU212 cells. After 24 h, we captured representative images of the migrating cells. (B) Quantitative analysis of untreated TU212 invasion effects compared to treated TU212 invasion effects. Invasive cells were counted on average.

Discussion

Cuproptosis, a newly discovered form of cell death, is involved in the accumulation of intracellular copper and possesses the common features of multiple regulated cell death (RCD) mechanisms. Genes associated with cuproptosis play an important role in the prognosis of various tumours such as melanoma and lung cancer. (Shanbhag et al., 2021; Yun et al., 2022). In this study, CDKN2A, a cuproptosis-associated gene, was identified as a novel biological marker for predicting the prognosis of patients with T2-stage HNSCC compared with T3-and T4-stage HNSCC. GSEA and GSVA revealed that CDK2NA overexpression was associated with the cell cycle. Positive correlations were found between CDKN2A expression and activated B and CD4 T cells, but negative correlations were found between CDKN2A expression and neutrophils and gamma-delta T cells. Molecular docking and MDS analysis revealed that plicamycin inhibits HNSCC progression by acting on CDKN2A.

CDKN2A contributes to cell proliferation and angiogenesis by being a member of the INK4 family of tumour suppressor genes (Zhang et al., 2016; Zhao et al., 2016). CDKN2A contains four exons, namely, exons 1α, 1β, 2 and 3. It possesses independent promoters for their respective proteins: exons 1, 2 and 3 encode P16INK4a and exons 1, 2 and 3 encode P14ARF, thus they are both independent proteins (Li et al., 2011). The cyclin-dependent kinase (CDK) inhibitor p16INK4a inhibits the phosphorylation of retinoblastoma protein (pRb) by binding to CDK4 and CDK6. Hypophosphorylated pRb inhibits the transcription factor family E2F, impairing cell cycle progression from G1 to S phase (Rubin, 2013; van den Heuvel and Dyson, 2008). By forming a trimeric complex with MDM2 and p53, P14ARF inhibits the degradation of P53 by MDM2, resulting in G1 and G2 arrest (Pollice et al., 2008). CDKN2A expression is correlated with the development and prognosis of various tumours, including hepatocellular carcinoma, pancreatic cancer and melanoma (Zeng et al., 2018; Kimura et al., 2021; Luo et al., 2021). Copy number deletion in CDKN2A and low expression of p16INK4a indicate a poor prognosis and can be used as independent prognostic predictors in HNSCC (Chen et al., 2018; Cury et al., 2021). A meta-analysis on HNSCC reported that methylation of CDKN2A is significantly correlated with tumorigenesis, progression and metastasis and can be used as a potential diagnostic marker (Zhou et al., 2018). Similar results were obtained in this study, indicating that CDKN2A plays an important role in regulating the tumour cell cycle and may serve as a prognostic marker in HNSCC.

Tumor-infiltrating immune cells are critical components of the tumour microenvironment and are highly predictive of prognosis and treatment outcomes (Nguyen et al., 2016; Guo et al., 2020). In this study, CDKN2A was positively correlated with the infiltration of activated B and CD4 T cells in the tumour microenvironment, leading to a better prognosis of HNSCC. The IGJ gene encodes CD19 and the J-chain gene, typical markers of B cells. In HNSCC, CD19 and J-chain expression were studied, and B cell infiltration in the tumor microenvironment indicated a better prognosis (Kim et al., 2020). Moreover, anti-blockade of PD-1/PD-L1 immune checkpoints in an AT-84-E7 murine model of HNSCC led to tumor enlargement by depleting B cells in the tumor microenvironment by regulating B cell activation and germinal center formation (Kim et al., 2020). This suggests that B cell infiltration and activity play an important role in the treatment of tumors. Hladíková et al., (2019) observed a high abundance of tumor-infiltrating B cells (TIL-Bs) in the HNSCC tumor microenvironment, which may interact with CD8⁺ T cells to promote tumor growth (). They produce tumour-associated antibodies and cytokines, which exert cytotoxic effects on tumours, as well as presenting tumour-associated antigens (TAA) (Candolfi et al., 2011; Mahmoud et al., 2012; Shi et al., 2013; Germain et al., 2014). Infection with the human papillomavirus (HPV) is causing an increase in HNSCC. HNSCC with HPV-positive microenvironment contain more CD4⁺ T cells, which result in a better prognosis and response to radiotherapy than those with HPV-negative microenvironment. However, whether this better prognostic performance is related to CD4⁺ T cell infiltration remains unclear (van Kempen et al., 2016). Multiple studies have shown that neutrophils can promote tumour development and angiogenesis through the secretion of cytokines such as vascular endothelial growth factor, hepatocyte growth factor, IL-6 and IL-8 (McCourt et al., 1999; McCourt et al., 2001; Schaider et al., 2003). Low CDKN2A expression was associated with neutrophil infiltration and poorer prognosis in the present study. Yang et al. reported in a previous study that high neutrophil infiltration can result in the generation of reactive oxygen species, arginase and nitric oxide, causing lymphocytes to be suppressed, T cells to be activated, and, ultimately, tumour growth (Gooden et al., 2011; Yang et al., 2019). Watermann et al., (2021) also found a link between neutrophil infiltration and tumour malignancy in the microenvironment of recurrent HNSCC.

There is a positive correlation between CDKN2A expression and activation of B cells and CD4⁺ T cells. There is evidence that CD8⁺ and CD4⁺ T cells infiltrate the microenvironment to ensure the beneficial effects of gamma-delta T cells on prognosis (Nielsen et al., 2017; Lu et al., 2020). Gamma-delta T cells are important immune cells in the mucosal region responsible for removing pathogens and maintaining the integrity of the epithelium. Lu et al. found a high abundance of gamma-delta T cells infiltration in the tumour microenvironment of patients with HNSCC with a better prognosis. However, Bas et al. examined gamma-delta T cell in the peripheral blood of patients with HNSCC and found that their high abundance was associated with the recurrence of HNSCC (Bas et al., 2006). Therefore, CDKN2A, an immune cell related gene, can aid in improving the prognosis and microenvironment of tumors by inhibiting immune cells that activate B and CD4 T cell, neutrophils and gamma-delta T cells.

Plicamycin, also known as mithramycin A, is a natural polycyclic aromatic polyketide compound that inhibits SP1 transcription factor binding to DNA, which interferes with biological processes like tumour cell proliferation, apoptosis, angiogenesis, invasion, and metastasis (Beishline and Azizkhan-Clifford, 2015; Schweer et al., 2021). Plicamycin is a natural polycyclic aromatic polyketide that inhibits SP1 transcription factor binding to DNA, impeding apoptosis, angiogenesis, invasion, and metastatic processes in cancer cells (Saha et al., 2015). In mouse models of lung cancer, low SP1 expression can effectively inhibit tumour growth and nicotine-induced lung cancer cell growth (Brown et al., 2013). The histone methyltransferase gene SETDB1 is an important gene in the development and metastasis of melanoma in vivo. Several studies have shown that plicamycin can effectively target the activity of SP-1 protein on the SETDB1 promoter to inhibit SETDB1 expression, thus offering a beneficial therapeutic strategy for melanoma (Federico et al., 2020). Plicamycin, an inhibitor of SP1, can inhibit tumour progression in several cancer types and has been used for the treatment of lung, breast and gastrointestinal tract cancers in phase II clinical trials, with good efficacy. However, the role and mechanisms of action of plicamycin in HNSCC remain to be elucidated. In this study, strong binding was observed between plicamycin and CDKN2A through molecular docking and MDS, suggesting that plicamycin improves the prognosis of HNSCC by targeting CDKN2A at the molecular level. However, the relationship between CDKN2A and SP1 warrants further investigation. The results of this study are mainly based on bioinformatics analysis and MDS, and they aren’t validated in vivo or in vitro. Additionally, plicamycin’s effects on HNSCC should be studied in large clinical trials.

Conclusion

CDKN2A, which is closely related to the maintenance of copper metabolic homeostasis in the body, is a biomarker of HNSCC and may improve its prognosis by regulating the cell cycle and immune cell infiltration. Further, plicamycin targets and binds CDKN2A, offering a novel strategy for the treatment of HNSCC in the future.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: Head and neck squamous cell carcinoma data were obtained from TCGA database (https://portal.gdc.cancer.gov).

Author contributions

KF and YL conceived and designed the analysis; KF, TL, and YL collected the data; TL, YL, and YD contributed data or analysis tools; KF, TL, YL, and YD performed the analysis; KF, TL, YL, and YD wrote the paper.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.1036408/full#supplementary-material

Supplementary Figure S1 | Correlation among immune cell types. (A) Heat map demonstrating immune cell infiltration in the high- and low-CDKN2A-expression groups; (B) Correlation among different immune cell types.

Supplementary Figure S2 | CDKN2A mutation sites.

Supplementary Figure S3 | PLDS shows a medium docking effect with a total of 104 residue sites.

References

Bao, Z., Yang, Z., Huang, Z., Zhou, Y., Cui, Q., and Dong, D. (2019). LncRNADisease 2.0: An updated database of long non-coding RNA-associated diseases. Nucleic Acids Res. 47, D1034–D1037. doi:10.1093/nar/gky905

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Bas, M., Bier, H., Schirlau, K., Friebe-Hoffmann, U., Scheckenbach, K., Balz, V., et al. (2006). Gamma–delta T-cells in patients with squamous cell carcinoma of the head and neck. Oral Oncol. 42, 691–697. doi:10.1016/j.oraloncology.2005.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Beishline, K., and Azizkhan-Clifford, J. (2015). Sp1 and the ‘hallmarks of cancer. FEBS J. 282, 224–258. doi:10.1111/febs.13148

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA A Cancer J. Clin. 68, 394–424. doi:10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, K. C., Perry, H. E., Lau, J. K., Jones, D. V., Pulliam, J. F., Thornhill, B. A., et al. (2013). Nicotine induces the up-regulation of the α7-nicotinic receptor (α7-nAChR) in human squamous cell lung cancer cells via the Sp1/GATA protein pathway. J. Biol. Chem. 288, 33049–33059. doi:10.1074/jbc.M113.501601

PubMed Abstract | CrossRef Full Text | Google Scholar

Candolfi, M., Curtin, J. F., Yagiz, K., Assi, H., Wibowo, M. K., Alzadeh, G. E., et al. (2011). B cells are critical to T-cell—mediated antitumor immunity induced by a combined immune-stimulatory/conditionally cytotoxic therapy for glioblastoma. Neoplasia 13, 947–960. doi:10.1593/neo.11024

PubMed Abstract | CrossRef Full Text | Google Scholar

Canning, M., Guo, G., Yu, M., Myint, C., Groves, M. W., Byrd, J. K., et al. (2019). Heterogeneity of the head and neck squamous cell carcinoma immune landscape and its impact on immunotherapy. Front. Cell Dev. Biol. 7, 52. doi:10.3389/fcell.2019.00052

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W. S., Bindra, R. S., Hayman, T., Husain, Z., Contessa, J. N., Mo, A., et al. (2018). CDKN2A copy number loss is an independent prognostic factor in HPV-negative head and neck squamous cell carcinoma. Front. Oncol. 8, 95. doi:10.3389/fonc.2018.00095

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Sun, Y., Luo, Z., Lin, J., Qi, B., Kang, X., et al. (2022). Potential mechanism underlying exercise upregulated circulating blood exosome miR-215-5p to prevent necroptosis of neuronal cells and a model for early diagnosis of alzheimer’s disease. Front. Aging Neurosci. 14, 860364. doi:10.3389/fnagi.2022.860364

PubMed Abstract | CrossRef Full Text | Google Scholar

Chikan, N. A., and Vipperla, B. (2015). KAISO inhibition: An atomic insight. J. Biomol. Struct. Dyn. 33, 1794–1804. doi:10.1080/07391102.2014.974072

PubMed Abstract | CrossRef Full Text | Google Scholar

Colaprico, A., Silva, T. C., Olsen, C., Garofano, L., Cava, C., Garolini, D., et al. (2016). TCGAbiolinks: An R/bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 44, e71. doi:10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

Cury, S. S., Miranda, P. M. d., Marchi, F. A., Canto, L. M. d., Chulam, T. C., Petersen, A. H., et al. (2021). Germline variants in DNA repair genes are associated with young-onset head and neck cancer. Oral Oncol. 122, 105545. doi:10.1016/j.oraloncology.2021.105545

PubMed Abstract | CrossRef Full Text | Google Scholar

Federico, A., Steinfass, T., Larribere, L., Novak, D., Moris, F., Nunez, L. E., et al. (2020). Mithramycin A and mithralog EC-8042 inhibit SETDB1 expression and its oncogenic activity in malignant melanoma. Mol. Ther. - Oncolytics 18, 83–99. doi:10.1016/j.omto.2020.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., ColoMbetM., , , Soerjomataram, I., Mathers, C., Parkin, D. M., PinerosM., , et al. (2019). Estimating the global cancer incidence and mortality in 2018: GLOBOCAN sources and methods. Int. J. Cancer 144, 1941–1953. doi:10.1002/ijc.31937

PubMed Abstract | CrossRef Full Text | Google Scholar

Germain, C., Gnjatic, S., Tamzalit, F., Knockaert, S., Remark, R., Goc, J., et al. (2014). Presence of B Cells in tertiary lymphoid structures is associated with a protective immunity in patients with lung cancer. Am. J. Respir. Crit. Care Med. 189, 832–844. doi:10.1164/rccm.201309-1611OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Gooden, M. J. M., de Bock, G. H., Leffers, N., Daemen, T., and Nijman, H. W. (2011). The prognostic influence of tumour-infiltrating lymphocytes in cancer: A systematic review with meta-analysis. Br. J. Cancer 105, 93–103. doi:10.1038/bjc.2011.189

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, L., Li, X., Liu, R., Chen, Y., Ren, C., and Du, S. (2020). TOX correlates with prognosis, immune infiltration, and T cells exhaustion in lung adenocarcinoma. Cancer Med. 9, 6694–6709. doi:10.1002/cam4.3324

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hildebrand, P. W., Rose, A. S., and Tiemann, J. K. S. (2019). Bringing molecular dynamics simulation data into view. Trends Biochem. Sci. 44, 902–913. doi:10.1016/j.tibs.2019.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hladíková, K., Koucky, V., Boucek, J., Laco, J., Grega, M., Hodek, M., et al. (2019). Tumor-infiltrating B cells affect the progression of oropharyngeal squamous cell carcinoma via cell-to-cell interactions with CD8+ T cells. J. Immunother. cancer 7, 261. doi:10.1186/s40425-019-0726-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollingsworth, S. A., and Dror, R. O. (2018). Molecular dynamics simulation for all. Neuron 99, 1129–1143. doi:10.1016/j.neuron.2018.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H.-Y., Lin, Y. C. D., Li, J., Huang, K. Y., Shrestha, S., Hong, H. C., et al. (2019). miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res. 48, D148. doi:10.1093/nar/gkz896

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Lin, W., Sun, Y., Wang, Q., He, S., Han, Z., et al. (2022). Quercetin targets VCAM1 to prevent diabetic cerebrovascular endothelial cell injury. Front. Aging Neurosci. 14, 944195. doi:10.3389/fnagi.2022.944195

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Z., Shi, J., Gao, Y., Cui, C., Zhang, S., Li, J., et al. (2019). HMDD v3.0: A database for experimentally supported human microRNA–disease associations. Nucleic Acids Res. 47, D1013–D1017. doi:10.1093/nar/gky1010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, K., and Murphy, D. (2013). Application of ggplot2 to pharmacometric graphics. CPT Pharmacometrics Syst. Pharmacol. 2, 79. doi:10.1038/psp.2013.56

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, Z.-H., Ren, W.-Z., Wang, H.-Q., Gao, W., and Yuan, B. (2022). Molecular subtyping based on cuproptosis-related genes and characterization of tumor microenvironment infiltration in kidney renal clear cell carcinoma. Front. Oncol. 12, 919083. doi:10.3389/fonc.2022.919083

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, X., Chen, Y., Yi, B., Yan, X., Jiang, C., Chen, B., et al. (2021). An integrative microenvironment approach for laryngeal carcinoma: The role of immune/methylation/autophagy signatures on disease clinical prognosis and single-cell genotypes. J. Cancer 12, 4148–4171. doi:10.7150/jca.58076

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, X., Sun, Y., Yi, B., Jiang, C., Yan, X., Chen, B., et al. (2022). Based on network pharmacology and molecular dynamics simulations, baicalein, an active ingredient of yiqi qingre ziyin method, potentially protects patients with atrophic rhinitis from cognitive impairment. Front. Aging Neurosci. 14, 880794. doi:10.3389/fnagi.2022.880794

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. S., Shen, S., Miyauchi, S., Sanders, P. D., Franiak-Pietryga, I., Mell, L., et al. (2020). B cells improve overall survival in HPV-associated squamous cell carcinomas and are activated by radiation and PD-1 blockade. Clin. Cancer Res. 26, 3345–3359. doi:10.1158/1078-0432.CCR-19-3211

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimura, H., Klein, A. P., Hruban, R. H., and Roberts, N. J. (2021). The role of inherited pathogenic CDKN2A variants in susceptibility to pancreatic cancer. Pancreas 50, 1123–1130. doi:10.1097/MPA.0000000000001888

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Poi, M. J., and Tsai, M.-D. (2011). Regulatory mechanisms of tumor suppressor P16 ^INK4A and their relevance to cancer. Biochemistry 50, 5566–5582. doi:10.1021/bi200642e

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Jun, L., Yan, H., Zhang, Q., and Wei, S. L. (2022). The prognostic value and immune landscape of a cuproptosis-related lncRNA signature in head and neck squamous cell carcinoma. Front. Genet. 13, 942785. doi:10.3389/fgene.2022.942785

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, W., Liu, J., Oualha, N., and Balesdent, M. (2021). Study on the effect of composite nanoparticles on corneal epithelial cell immune mechanism based on dectin-1 signaling pathway. Oxid. Med. Cell Longev. 2021, 1148–1153. doi:10.1166/jnn.2021.18703

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, W., Wang, Q., Chen, Y., Wang, N., Ni, Q., Qi, C., et al. (2022). Identification of a 6-RBP gene signature for a comprehensive analysis of glioma and ischemic stroke: Cognitive impairment and aging-related hypoxic stress. Front. Aging Neurosci. 14, 951197. doi:10.3389/fnagi.2022.951197

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, H., Dai, W., Guo, J., Wang, D., Wen, S., Yang, L., et al. (2020). High abundance of intratumoral γδ T cells favors a better prognosis in head and neck squamous cell carcinoma: A bioinformatic analysis. Front. Immunol. 11, 573920. doi:10.3389/fimmu.2020.573920

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, J., Wang, J., and Huang, J. (2021). CDKN2A is a prognostic biomarker and correlated with immune infiltrates in hepatocellular carcinoma. Biosci. Rep. 41, BSR20211103. doi:10.1042/BSR20211103

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Q., Ma, H., Guo, E., Yu, L., Jia, L., Zhang, B., et al. (2022). MicroRNAs promote the progression of sepsis-induced cardiomyopathy and neurovascular dysfunction through upregulation of NF-kappaB signaling pathway-associated HDAC7/ACTN4. Front. Neurol. 13, 909828. doi:10.3389/fneur.2022.909828

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Z., He, Z., Qin, H., Chen, Y., Qi, B., Lin, J., et al. (2022). Exercise-induced IL-15 acted as a positive prognostic implication and tumor-suppressed role in pan-cancer. Front. Pharmacol. 13, 1053137. doi:10.3389/fphar.2022.1053137

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, H., Liu, X., Zeng, X., Liu, Y., Zhang, C., Zhang, Q., et al. (2022). Comprehensive analysis of cuproptosis-related genes in immune infiltration and prognosis in melanoma. Front. Pharmacol. 13, 930041. doi:10.3389/fphar.2022.930041

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahmoud, S. M. A., Lee, A. H. S., Paish, E. C., Macmillan, R. D., Ellis, I. O., and Green, A. R. (2012). The prognostic significance of B lymphocytes in invasive carcinoma of the breast. Breast Cancer Res. Treat. 132, 545–553. doi:10.1007/s10549-011-1620-1

PubMed Abstract | CrossRef Full Text | Google Scholar

McCourt, M., Wang, J. H., Sookhai, S., and Redmond, H. P. (2001). Activated human neutrophils release hepatocyte growth factor/scatter factor. Eur. J. Surg. Oncol. (EJSO) 27, 396–403. doi:10.1053/ejso.2001.1133

PubMed Abstract | CrossRef Full Text | Google Scholar

McCourt, M., Wang, J. H., Sookhai, S., and Redmond, H. P. (1999). Proinflammatory mediators stimulate neutrophil-directed angiogenesis. Arch. Surg. 134, 1325–1331. doi:10.1001/archsurg.134.12.1325

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, S., Li, Y., and Kang, X. (2022). Prognostic and functional analysis of NPY6R in uveal melanoma using bioinformatics. Dis. Markers 2022, 4143447–4143514. doi:10.1155/2022/4143447

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, N., Bellile, E., Thomas, D., McHugh, J., Rozek, L., Virani, S., et al. (2016). Tumor infiltrating lymphocytes and survival in patients with head and neck squamous cell carcinoma: Tumor infiltrating lymphocytes. Head. Neck 38, 1074–1084. doi:10.1002/hed.24406

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, M. M., Witherden, D. A., and Havran, W. L. (2017). γδ T cells in homeostasis and host defence of epithelial barrier tissues. Nat. Rev. Immunol. 17, 733–745. doi:10.1038/nri.2017.101

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveri, V. (2022). Selective targeting of cancer cells by copper ionophores: overview. Front. Mol. Biosci. 9, 841814. doi:10.3389/fmolb.2022.841814

PubMed Abstract | CrossRef Full Text | Google Scholar

Ou, J., and Zhu, L. J. (2019). trackViewer: a Bioconductor package for interactive and integrative visualization of multi-omics data. Nat. Methods 16, 453–454. doi:10.1038/s41592-019-0430-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollice, A., Vivo, M., and Mantia, G. L. (2008). The promiscuity of ARF interactions with the proteasome. FEBS Lett. 582, 3257–3262. doi:10.1016/j.febslet.2008.09.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubin, S. M. (2013). Deciphering the retinoblastoma protein phosphorylation code. Trends Biochem. Sci. 38, 12–19. doi:10.1016/j.tibs.2012.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz, L. M., Libedinsky, A., and Elorza, A. A. (2021). Role of copper on mitochondrial function and metabolism. Front. Mol. Biosci. 8, 711227. doi:10.3389/fmolb.2021.711227

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, S., Mukherjee, S., Mazumdar, M., Manna, A., Khan, P., Adhikary, A., et al. (2015). Mithramycin A sensitizes therapy-resistant breast cancer stem cells toward genotoxic drug doxorubicin. Transl. Res. 165, 558–577. doi:10.1016/j.trsl.2014.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaider, H., Oka, M., Bogenrieder, T., Nesbit, M., Satyamoorthy, K., Berking, C., et al. (2003). Differential response of primary and metastatic melanomas to neutrophils attracted by IL-8. Int. J. Cancer 103, 335–343. doi:10.1002/ijc.10775

PubMed Abstract | CrossRef Full Text | Google Scholar

Schweer, D., McCorkle, J. R., Rohr, J., Tsodikov, O. V., Ueland, F., and Kolesar, J. (2021). Mithramycin and analogs for overcoming cisplatin resistance in ovarian cancer. Biomedicines 9, 70. doi:10.3390/biomedicines9010070

PubMed Abstract | CrossRef Full Text | Google Scholar

Shanbhag, V. C., Gudekar, N., Jasmer, K., Papageorgiou, C., Singh, K., and Petris, M. J. (2021). Copper metabolism as a unique vulnerability in cancer. Biochimica Biophysica Acta (BBA) - Mol. Cell Res. 1868, 118893. doi:10.1016/j.bbamcr.2020.118893

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, J.-Y., Gao, Q., Wang, Z. C., Zhou, J., Wang, X. Y., Min, Z. H., et al. (2013). Margin-infiltrating CD20+ B cells display an atypical memory phenotype and correlate with favorable prognosis in hepatocellular carcinoma. Clin. Cancer Res. 19, 5994–6005. doi:10.1158/1078-0432.CCR-12-3497

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Q., Yan, X., Wang, J., and Zhang, X. (2021). Collagen family genes associated with risk of recurrence after radiation therapy for vestibular schwannoma and pan-cancer analysis. Dis. Markers 2021, 7897994–7898015. doi:10.1155/2021/7897994

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: Customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49, D605–D612. doi:10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

Taberna, M., MenaM., , Pavon, M. A., Alemany, L., Gillison, M. L., and Mesia, R. (2017). Human papillomavirus-related oropharyngeal cancer. Ann. Oncol. 28, 2386–2398. doi:10.1093/annonc/mdx304

PubMed Abstract | CrossRef Full Text | Google Scholar

Thai, K.-M., Le, D. P., Tran, N. V. K., Nguyen, T. T. H., and Tran, T. D. (2015). Computational assay of Zanamivir binding affinity with original and mutant influenza neuraminidase 9 using molecular docking. J. Theor. Biol. 385, 31–39. doi:10.1016/j.jtbi.2015.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomczak, K., Czerwińska, P., and Wiznerowicz, M. (2015). Review the cancer genome atlas (TCGA): An immeasurable source of knowledge. Contemp. Oncol. Pozn. 19, 68–77. doi:10.5114/wo.2014.47136

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 375, 1254–1261. doi:10.1126/science.abf0529

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Heuvel, S., and Dyson, N. J. (2008). Conserved functions of the pRB and E2F families. Nat. Rev. Mol. Cell Biol. 9, 713–724. doi:10.1038/nrm2469

PubMed Abstract | CrossRef Full Text | Google Scholar

van Kempen, P. M. W., Noorlag, R., Swartz, J. E., Bovenschen, N., Braunius, W. W., Vermeulen, J. F., et al. (2016). Oropharyngeal squamous cell carcinomas differentially express granzyme inhibitors. Cancer Immunol. Immunother. 65, 575–585. doi:10.1007/s00262-016-1819-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Warde-Farley, D., Donaldson, S. L., Comes, O., Zuberi, K., Badrawi, R., Chao, P., et al. (2010). The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 38, W214–W220. doi:10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

Watermann, C., Pasternack, H., Idel, C., Ribbat-Idel, J., Bragelmann, J., Kuppler, P., et al. (2021). Recurrent HNSCC harbor an immunosuppressive tumor immune microenvironment suggesting successful tumor immune evasion. Clin. Cancer Res. 27, 632–644. doi:10.1158/1078-0432.CCR-20-0197

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, C.-Y., Zhu, M. X., Lu, N. H., Peng, R., Yang, X., Zhang, P. F., et al. (2019). Bioinformatics-based analysis reveals elevated MFSD12 as a key promoter of cell proliferation and a potential therapeutic target in melanoma. Oncogene 38, 1876–1891. doi:10.1038/s41388-018-0531-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiser, J., Shenkin, P. S., and Still, W. C. (1999). Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 20, 217–230. doi:10.1002/(sici)1096-987x(19990130)20:2<217::aid-jcc4>3.0.co;2-a

CrossRef Full Text | Google Scholar

Xuan, Z., Ma, T., Qin, Y., and Guo, Y. (2022). Role of ultrasound imaging in the prediction of TRIM67 in brain metastases from breast cancer. Front. Neurol. 13, 889106. doi:10.3389/fneur.2022.889106

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Huang, Y., Zhou, L., Dai, Y., and Hu, G. (2019). High pretreatment neutrophil-to-lymphocyte ratio as a predictor of poor survival prognosis in head and neck squamous cell carcinoma: Systematic review and meta-analysis. Head Neck 41, 1525–1535. doi:10.1002/hed.25583

PubMed Abstract | CrossRef Full Text | Google Scholar

Yun, Y., Wang, Y., Yang, E., and Jing, X. (2022). Cuproptosis-related gene – SLC31A1, FDX1 and ATP7B – polymorphisms are associated with risk of lung cancer. PGPM 15, 733–742. doi:10.2147/PGPM.S372824

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, H., Jorapur, A., Shain, A. H., Lang, U. E., Torres, R., Zhang, Y., et al. (2018). Bi-Allelic loss of CDKN2A initiates melanoma invasion via BRN2 activation. Cancer Cell 34, 56e9–68. doi:10.1016/j.ccell.2018.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Zeng, M., and Fu, B. M. (2016). Inhibition of endothelial nitric oxide synthase decreases breast cancer cell MDA-MB-231 adhesion to intact microvessels under physiological flows. Am. J. Physiology-Heart Circulatory Physiology 310, H1735–H1747. doi:10.1152/ajpheart.00109.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zhang, J., Sun, C., and Wu, F. (2022). Identification of the occurrence and potential mechanisms of heterotopic ossification associated with 17-beta-estradiol targeting MKX by bioinformatics analysis and cellular experiments. PeerJ 9, e12696. doi:10.7717/peerj.12696

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B., and Jiang, X. (2022). hsa-miR-518-5p/hsa-miR-3135b regulates the REL/SOD2 pathway in ischemic cerebral infarction. Front. Neurol. 13, 852013. doi:10.3389/fneur.2022.852013

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, R., Choi, B. Y., Lee, M.-H., Bode, A. M., and Dong, Z. (2016). Implications of genetic and epigenetic alterations of CDKN2A (p16 INK4a ) in cancer. EBioMedicine 8, 30–39. doi:10.1016/j.ebiom.2016.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Shen, Z., Ye, D., Li, Q., Deng, H., Liu, H., et al. (2018). The association and clinical significance of CDKN2A promoter methylation in head and neck squamous cell carcinoma: A meta-analysis. Cell Physiol. Biochem. 50, 868–882. doi:10.1159/000494473

PubMed Abstract | CrossRef Full Text | Google Scholar