The RNA-sequencing (RNA-seq) expression profiles of LUAD were downloaded from TCGA database (https://portal.gdc.cancer.gov) (25), comprising 530 LUAD tissue samples (T) and adjacent tissues samples from 59 patients with LUAD (N). Utilizing the rjson R package (https://github.com/alexcb/rjson) through R studio (2023.09.1 494) (26), the gene expression matrices were integrated, duplicate data were removed and gene names were converted to common identifiers, thereby acquiring LUAD-associated genes. CRGs were sourced from prior literature studies (9,14,15). SangerBox 3.0 (http://vip.sangerbox.com/) is a comprehensive, user-friendly bioinformatics analysis platform; this software was used to perform image visualization (27).

To analyze TCGA data (T=530; N=59), differential expression analysis was performed using R packages limma-voom (28), DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) and edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html) to identify differentially expressed genes (DEGs) between the selected groups and their respective controls. Genes with expression values of zero in >50% of samples were excluded. To perform differential expression analysis using DESeq2 on the obtained TCGA expression dataset, the DESeqDataSetFromMatrix and DESeq functions were utilized to input the matrix and normalize the data, followed by the results function to determine the significance of differential expression for each gene. For differential expression analysis using edgeR, the DGEList and calcNormFactors functions were adopted to input the matrix and normalize the data, with the exactTest function used to assess the significance of differential expression for each gene. In the limma approach for differential expression analysis, the data were transformed using the voom function, and multivariate linear regression analysis was performed using the lmFit function. The eBayes function was then applied to calculate moderated statistical values, yielding the significance of differential expression for each gene. Hypothesis testing approaches were utilized to set P-values, false discovery rate (FDR) values and fold changes for screening DEGs. In the present study, P<0.05, FDR <0.05 and fold change >1.5 were applied to obtain statistically significant differences for each gene. Common DEGs from the three methods were intersected using Venn diagrams to ensure accuracy. These DEGs were then intersected with CRGs to identify 10 candidate LUAD CRGs.

The STRING database (version 11.5; http://string-db.org/) (29) was employed to construct a PPI network among the differentially co-expressed genes. When constructing the gene interaction network using STRING to illustrate the PPI associations among the 10 selected genes, a confidence score of ≥0.150 was applied as the threshold. The resulting network was visualized using Cytoscape (v3.9.1) (30). Among the network analysis methods, node degree and maximal clique centrality (MCC) algorithms were applied via the cytoHubba plug-in in Cytoscape to identify potential hub nodes (31). Based on these analyses, CDKN2A, a LUAD CRG, was selected as the most credible hub gene.

A total of five databases were used to predict the potential TFs for CDKN2A and the nine other LUAD CRGs identified (AOC3, ULK2, SLC31A2, PDK1, CP, GCSH, COA6, LOXL2 and H3C1): ENCODE (https://www.encodeproject.org/), hTFtarget (http://bioinfo.life.hust.edu.cn/hTFtarget), Cistrome (http://cistrome.org/db/), TCGA-LUAD (https://portal.gdc.cancer.gov) and Genotype-Tissue Expression (GTEx) Lung (https://gtexportal.org/home). Based on Pearson correlation analysis, the correlation was calculated between the expression of target genes and TFs to ultimately obtain the required TFs. In the present study, the prediction of target gene TFs across multiple databases was performed using the website https://jingege.shinyapps.io/TF_predict/. Subsequently, the overlapping TFs predicted by ≥2 databases were identified through Venn diagrams. Among them, Spi-1 proto-oncogene (SPI1) was selected for prognostic analysis based on its prognostic significance.

The Gene Expression Profiling Interactive Analysis (GEPIA) 2.0 database (http://gepia.cancer-pku.cn/) (32) encapsulates a vast collection of data from 9,736 tumor samples and 8,587 normal samples sourced from TCGA and GTEx projects. In the present study, the GEPIA 2.0 database facilitated the analysis of differential gene expression levels between LUAD tissues and normal lung tissues, as well as correlations among genes. The Human Protein Atlas 24.0 (HPA 24.0; images available from v24.proteinatlas.org) (33) was utilized to further corroborate the protein expression levels of specific genes in LUAD vs. normal lung tissues by visually inspecting the protein expression of relevant genes in pathological sections of normal lung and LUAD tissues, thereby verifying the differences in protein expression.

The UALCAN database (http://ualcan.path.uab.edu/) (34), an online analysis and mining tool based on TCGA data, was used to verify the reliability of hub genes by analyzing the differences in mRNA expression levels between LUAD tissues and normal tissues, as well as their prognostic significance. Additionally, the Kaplan-Meier (KM) plotter database (http://kmplot.com/analysis/) (35) was employed to examine the effect of gene expression on overall survival (OS), first progression (FP) and post-progression survival (PPS) in patients with LUAD; the patients were divided according to median values, and the dataset names were 211156-at, 207039-at and 209644-at. The test method used log-rank test to compare between groups and to calculate P-values, and then Cox multivariate analysis was performed to calculate the hazard ratio to determine the independent prognostic role of genes; the follow-up threshold was 120 months.

LinkedOmics 1.2 (http://www.linkedomics.org/login.php) (36), a visualization platform, was used to explore gene expression profiles. Using LinkedOmics, co-expressed genes of CDKN2A were identified through Spearman's correlation coefficient analysis. The results were visualized using heatmaps and volcano plots. Subsequently, Gene Set Enrichment Analysis was performed to investigate the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with CDKN2A and its co-expressed genes. The data were downloaded from LinkedOmics and used HelixLife 2.0.28 (https://www.helixlife.cn/) to visualize the GO graph. The KEGG database (https://www.genome.jp/kegg/) (37), which interconnects disease genes, pathways, drugs and diagnostic markers, served as the primary resource for signal pathway enrichment analysis.

To assess the credibility of core genes, the normal dataset (accession number: 000319563900002) from GTEx was incorporated into the present analysis. A comprehensive analysis of RNA-seq data in transcripts per million format was performed from TCGA and GTEx, which were uniformly processed using the Toil pipeline, in the University of California Santa Cruz Xena database (https://xenabrowser.net/datapages/) (38). The final analysis included 515 LUAD samples from TCGA, as well as 59 normal tissue samples from TCGA and an additional 288 normal tissue samples from GTEx. The Wilcoxon rank-sum test was selected for statistical analysis, and the ggplot2 3.4.4 package (https://ggplot2.tidyverse.org) was utilized for data visualization to validate the significance of the core genes.

The normal human lung epithelial BEAS-2B cell line and the LUAD H1975 cell line were obtained from The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences. Cells in the logarithmic growth phase were seeded into 6-well plates containing RPMI-1640 (cat. no. PM150110; Wuhan Pricella Biotechnology Co., Ltd.; Wuhan Elabscience Biotechnology Co., Ltd.) medium supplemented with 10% fetal bovine serum (Wuhan Pricella Biotechnology Co., Ltd.; Wuhan Elabscience Biotechnology Co., Ltd.) and 1% penicillin-streptomycin (Beyotime Institute of Biotechnology). When the cells reached 60–70% confluence, plasmids (1.5 µg) were transfected into the cells using Lipofectamine^® 2000 (Thermo Fisher Scientific, Inc.). The transfection ratio was 1:1.5, Lipofectamine 2000:plasmid. The plasmid backbone names and sequences are presented in Table I. The short hairpin (sh)RNA-CDKN2A (sh-CDKN2A), sh-negative control (NC; empty plasmid), overexpression (oe)-SPI1 and oe-NC (empty plasmid) plasmids (1.5 µg) were transfected into the H1975 cell line and cultured at 37°C in a 5% CO₂ incubator for 48 h. The transfection efficiency was assessed by reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

Total RNA was extracted from BEAS-2B and H1975 cells using RNAiso Plus (cat. no. 9109; Takara Bio, Inc.) and RNA concentrations were quantified with an ultra-micro spectrophotometer. According to the manufacturer's protocol, RT of RNA was performed using HiScript II Q RT SuperMix for qPCR (+gDNA wiper; cat. no. R223; Vazyme Biotech Co., Ltd.). SYBR qPCR Master Mix (cat. no. Q311; Vazyme Biotech Co., Ltd.) was used for qPCR analysis. According to the manufacturer's protocol the following thermocycling conditions were performed: Pre-denaturation at 95°C for 30 sec; followed by 40 cycles of denaturation at 95°C for 10 sec, annealing at 60°C for 30 sec and extension at 72°C for 30 sec. Relative gene expression levels were analyzed using the 2^−ΔΔCt method (39). GAPDH served as the internal reference gene with the following primer sequences: Forward, 5′-TGCACCACCAACTGCTTAGC-3′ and reverse, 5′-GGCATGGACTGTGGTCATGAG-3′. The primer sequences for CDKN2A were: Forward, 5′-CTTCCTCGGGTGCCGATAC-3′ and reverse, 5′-ACCCCTTCATTGCTACTCGAT-3′. For SPI1, the primer sequences were: Forward, 5′-GTGCCCTATGACACGGATCTA-3′ and reverse, 5′-AGTCCCAGTAATGGTCGCTAT-3′.

Transfected cells were lysed on ice using RIPA lysis buffer (Beyotime Institute of Biotechnology) to extract proteins, and protein concentration was determined using the BCA kit (cat. no. P0012; Beyotime Institute of Biotechnology). Proteins (20 µg) were subjected to sodium dodecyl sulfate-polyacrylamide gel electrophoresis on a 12% gel and transferred to a PVDF membrane using conventional wet transfer methods. The membrane was blocked with 5% skimmed milk powder in TBS-0.1% Tween (TBST) at room temperature for 2 h. An overnight incubation with primary antibodies at 4°C followed. The next day, the membrane was washed three times with TBST and incubated with secondary antibodies for 2 h at 4°C. After another three washes with TBST, enhanced chemiluminescence solution (Biosharp Life Sciences) was added for exposure. ImageJ software (1.54f; National Institutes of Health) was used to measure the gray values of each band, with GAPDH serving as the internal control. Relative protein expression levels were semi-quantified as the ratio of gray values of the target protein to gray values of GAPDH. The antibodies used were in the present study were as follows: Rabbit anti-CDKN2A antibody (1:1,000; cat. no. BS40808; Bioworld Technology, Inc.), rabbit anti-p53 antibody (1:2,000; cat. no. 21891-1-AP; Proteintech Group, Inc.), rabbit anti-GAPDH antibody (1:3,000; cat. no. AP0063; Bioworld Technology, Inc.) and a Goat Anti-Rabbit IgG secondary antibody (1:6,000; cat. no. BS13278; Bioworld Technology, Inc.).

Cell viability was assessed using the CCK-8 assay kit (Wuhan Elabscience Biotechnology Co., Ltd.). Transfected H1975 cells (sh-NC, sh-CDKN2A, oe-NC and oe-SPI1) were seeded into 96-well plates at a density of 2×10³ cells/well and incubated at 37°C. Each group contained four replicate wells. The cells were incubated for 0, 24, 48 and 72 h, with 10 µl CCK-8 solution added to each well at each time point, followed by 2 h incubation. Absorbance at 450 nm was measured using a microplate reader.

Cells transfected with sh-NC, sh-CDKN2A, oe-NC and oe-SPI1 were seeded into 6-well plates at a density of 5×10⁵ cells/well and cultured to form a monolayer. A scratch was created using a 10-µl pipette tip, and the cells were maintained in RPMI-1640 medium containing 1% FBS. The scratch area was observed and images were captured using an optical microscope at 0, 24 and 48 h. The scratch areas were marked using Photoshop software (version 23.0.0; Adobe Systems, Inc.).

The LUAD cell line H1975 was selected and maintained in RPMI-1640 medium under standard humidified incubator conditions at 37°C with 5% CO₂. Three distinct groups were established within the experiment: A control group, an ES-CuCl₂ group and a tetrathiomolybdate (TTM)-ES-CuCl₂ group. ES, an efficacious Cu²⁺ ionophore that enhances cellular apoptosis, was sourced from MedChemExpress. TTM, a copper chelator that functions as a copper antagonist, was procured from Shanghai Yuanye Biotechnology Co., Ltd. ES was dissolved in 100% dimethyl sulfoxide (Beijing Solarbio Science & Technology Co., Ltd.) and the control group was treated with an equal amount of dimethyl sulfoxide. CuCl₂ (Shanghai Macklin Biochemical Co., Ltd.) was dissolved in sterile water. ES-CuCl₂ was prepared by mixing ES and CuCl₂ in a 1:1 ratio. In the ES-CuCl₂ group, H1975 cells were treated with 30 nM ES-CuCl₂ for 2 h. By comparison, cells in the TTM-ES-CuCl₂ group were pretreated overnight with 20 µM TTM before exposure to 30 nM ES-CuCl₂ for 2 h. After the aforementioned treatment, the cells of each group were cultured in fresh medium at 37°C and 5% CO₂ for 48 h, and RNA was extracted. Gene expression was quantitatively detected by RT-qPCR.

For gene pathway analysis (n=515) and core gene validation and correlation analysis (T=483; N=347), Pearson correlation analysis was used. When analyzing the gene expression differences between the normal group and the cancer group (T=483, N=347) using the GEPIA database, an independent samples t-test was employed. When comparing differences between normal and cancer groups in both TCGA and GTEx datasets (T=515; N=347), non-parametric Wilcoxon rank-sum test were performed, using disease status as the variable for calculating differential expression. For patient survival analysis, the log-rank test and Cox univariate regression were utilized. In experimental validation, differences between groups were assessed using one-way ANOVA with Tukey's post hoc test. A unpaired Student's t-test was used to compare differences between two groups. Results are expressed as the mean ± standard deviation. P<0.05 was considered to indicate a statistically significant difference. All calculated results were statistically analyzed using GraphPad Prism 9.1 software (Dotmatics). All experiments were repeated three times.

Differential gene expression analysis was performed on 59,427 genes associated with LUAD, retrieved from TCGA database using the R packages limma-voom, DESeq2 and edgeR. The datasets employed in all three methods excluded genes with >50% zero expression values from TCGA dataset. All analyses were executed within the SangerBox software platform. Genes with P<0.05, FDR <0.05 and fold change >1.5 were selected and sorted based on their P-values. According to the limma-voom method, 13,651 significantly DEGs were identified between LUAD and normal lung tissues (Fig. 1A). A total of 16,559 DEGs were detected using DESeq2 (Fig. 1B), whereas edgeR analysis yielded 12,151 DEGs (Fig. 1C). The intersection of these three gene sets revealed 10,252 common DEGs (Fig. 1D). These common DEGs were subsequently intersected with a set of 45 CRGs derived from the literature (Fig. 1E), resulting in 10 candidate genes (Fig. 1F). Among these candidates, three genes (AOC3, ULK2 and SLC31A2) were downregulated, whereas seven genes (PDK1, CP, GCSH, COA6, LOXL2, CDKN2A and H3C1) were upregulated (Fig. 1G). The statistical values associated with genes in the heatmap all adhered to P<0.05, FDR <0.05 and fold change >1.5.

The STRING database was utilized to investigate potential interactions among candidate genes, and the results were visualized using Cytoscape. A PPI network was constructed by importing 10 candidate genes into the STRING database, which expanded to 15 nodes and 34 edges (Fig. 2A). Hub genes were identified using the MCC algorithm from the cytoHubba plug-in within Cytoscape. Genes with high scores are represented by red nodes, whereas those with low scores are denoted by yellow nodes (Fig. 2B). Based on image analysis, the LUAD CRG CDKN2A, with a score of 252, was selected as the hub gene for further investigation (Fig. 2C). The UALCAN database was employed to predict the effect of CDKN2A expression on OS of patients with LUAD, samples were categorized into two groups: High-expression (with TPM values above upper quartile) and low/medium expression (with TPM values below upper quartile), revealing that patients in the CDKN2A high-expression group exhibited worse OS compared with the low/medium expression-group (P=0.0023; Fig. 2D). The KM plotter database was used to assess the influence of CDKN2A (also known as multiple tumor suppressor 1) expression on OS, PPS and FP in LUAD, thereby validating the findings (Fig. 2E-M). The results indicated that of the nine datasets analyzed, only two showed non-significant PPS (P=0.14, Fig. 2H; P=0.1, Fig. 2J), whereas all other datasets demonstrated statistical significance (P<0.05). Due to the relatively small sample size of ~300 cases for PPS compared with ~1,000 cases for OS and FP, there may be potential biases or errors in the PPS results. Nevertheless, high CDKN2A expression was associated with unfavorable patient prognosis, which supports the potential of CDKN2A levels as a predictive marker for long-term survival in patients with LUAD.

The GEPIA database predicted the clinical specimen tissue expression levels of CDKN2A between LUAD tissues and normal tissues. The results indicated that CDKN2A was significantly upregulated in LUAD tissues compared with in normal tissues (Fig. 3A). RT-qPCR validation revealed that the mRNA expression levels of CDKN2A in H1975 cells were significantly higher compared with those in BEAS-2B cells (Fig. 3B). The HPA database was utilized to assess and verify the presence of protein differences between normal lung tissues and LUAD tissues for CDKN2A, further confirming its role in LUAD tissues; the results revealed that CDKN2A protein expression was increased in LUAD tissues than in normal tissues (Fig. 3C). Following the knockdown of CDKN2A in H1975 cells using sh-CDKN2A, RT-qPCR validation of knockdown efficiency demonstrated that the expression levels of CDKN2A in the sh-CDKN2A group were significantly lower compared with those in the sh-NC group (Fig. 3D). The CCK-8 assay showed that CDKN2A knockdown inhibited cell proliferation (Fig. 3E). The wound healing assay assessing cell migration demonstrated reduced migratory ability of H1975 cells after CDKN2A knockdown (Fig. 3F). These findings indicated that CDKN2A knockdown may inhibit the proliferative and migratory capacity of H1975 cells.

To investigate the potential functions and mechanisms of CDKN2A in LUAD, analysis of CDKN2A co-expressed genes was performed using LinkedOmics on RNA-seq data from patients with LUAD in TCGA database. The resulting gene correlation heatmap analysis revealed that CDKN2A expression was positively correlated with the top 50 genes, including MTAP, CDKN2B, CCNE1 and RFC4 (Fig. 4A), and was negatively correlated with the top 50 genes, such as EPAS1, ATF7, MTUS1 and SIAE (Fig. 4B). Subsequently, the KEGG pathways and GO terms enriched in the co-expressed genes of CDKN2A were analyzed. The KEGG pathway analysis indicated that the co-expressed genes of CDKN2A primarily participated in the ‘cell cycle’, ‘DNA replication’, ‘p53 signaling pathway’ and ‘cAMP signaling pathway’ (Fig. 4C), All mentioned KEGG pathways were P<0.05. The GO functional analysis revealed that CDKN2A co-expressed genes participated in biological processes such as ‘DNA replication’ and ‘DNA-templated DNA replication’; cellular components such as ‘nuclear chromosome’, ‘chromosomal region’ and ‘spindle’; and molecular functions including ‘catalytic activity, acting on DNA’, ‘protein kinase regulator activity’ and ‘single-stranded DNA binding’ (Fig. 4D). All mentioned GO terms were P<0.05. The p53 pathway related to cuproptosis was selected for further study with FDR<0.05. Further exploration of the KEGG database demonstrated that CDKN2A exerts its effects by regulating different genes through the p53 pathway (Fig. 4E).

Utilizing the hTFtarget database, 17 potential co-TFs were predicted for nine CRGs, including AOC3, ULK2, SLC31A2, PDK1, CP, GCSH, COA6, LOXL2 and H3C1 (Fig. 5A). A total of 25 common TFs were identified for CDKN2A by integrating data from the ENCODE, hTFtarget, Cistrome, TCGA-LUAD and GTEx Lung databases (Fig. 5B). The intersection of these two predicted TF sets suggested that SPI1 and SRF might serve roles in the transcriptional regulation of CDKN2A (Fig. 5C). Prognostic survival analysis of SPI1 and SRF using the UALCAN database revealed that SRF did not show prognostic significance (P=0.74; Fig. 5D). By contrast, the prognostic survival analysis of SPI1 revealed a statistically significant difference between the high and low expression groups; notably, the prognosis of the low expression group was poor (P=0.016; Fig. 5E).

The tissue expression level predictions from clinical specimens revealed significant differences in SPI1 expression between normal and LUAD tissues, with LUAD tissues exhibiting lower SPI1 expression levels compared with normal tissues (Fig. 6A). RT-qPCR validation demonstrated that SPI1 expression in H1975 cells was significantly lower compared with that in BEAS-2B cells (Fig. 6B). Furthermore, the HPA database was utilized to assess and validate the differential SPI1 protein levels between normal lung tissues and LUAD tissues, highlighting its role in LUAD; the results indicated that SPI1 protein expression was lower in LUAD tissues than in normal tissues (Fig. 6C). Following the overexpression of SPI1 in H1975 cells using oe-SPI1, RT-qPCR analysis revealed a significant increase in SPI1 expression levels in the oe-SPI1 group compared with those in the oe-NC group (Fig. 6D). CCK-8 assays evaluating cell viability indicated that SPI1 overexpression inhibited cell proliferation (Fig. 6E). Wound healing assays assessing cell migration demonstrated reduced migration of H1975 cells in the oe-SPI1 group compared with that in the oe-NC group (Fig. 6F). These findings indicated that high SPI1 expression may inhibit cell proliferation and migration.

A cuproptosis induction model was established (Fig. 7A) to investigate whether the expression of SPI1 and CDKN2A is influenced by cuproptosis. The results indicated that ES-CuCl₂ induced cuproptosis leading to significant upregulation of SPI1 expression and downregulation of CDKN2A expression. However, when H1975 cells were pretreated with TTM, the upregulation of SPI1 was inhibited and the downregulation of CDKN2A was also suppressed (Fig. 7B). Furthermore, the sensitivity of H1975 cells to cuproptosis activators was assessed after knocking down CDKN2A or overexpressing SPI1 to explore the roles of SPI1 and CDKN2A in LUAD cuproptosis. CCK-8 assay results demonstrated that SPI1 overexpression (Fig. 7C) and CDKN2A knockdown (Fig. 7D) promoted ES-CuCl₂-induced cuproptosis in H1975 cells.

Through integrated analysis of lung tissue data from TCGA-GTEx, boxplots showed statistically significant differences in SPI1 (Fig. 8A), CDKN2A (Fig. 8B) and TP53 (Fig. 8C) expression (P<0.001) between normal lung tissues and LUAD tissues. The expression levels of SPI1 were significantly lower in LUAD tissues than those in normal tissues, whereas the expression levels of CDKN2A and p53 were significantly higher in LUAD tissues than those in normal tissues. The GEPIA database was used to predict the correlations between genes to demonstrate the regulatory association within the SPI1/CDKN2A/p53 signaling axis. The results showed that SPI1, as a TF, was negatively correlated with CDKN2A (P<0.01; Fig. 8D) and TP53 (P<0.01; Fig. 8E), whereas CDKN2A was positively correlated with TP53 (P<0.01; Fig. 8F). Notably, the correlations between SPI1 and CDKN2A, and between CDKN2A and p53 were considered weak (r-values <0.3/-0.3). Therefore, western blot analysis was performed for verification. The overexpression of SPI1 in H1975 cells resulted in decreased protein levels of CDKN2A and p53 (Fig. 8G), and the knockdown of CDKN2A decreased p53 protein levels (Fig. 8H). These findings suggested that SPI1 may promote the cuproptosis of LUAD cells. Whereas CDKN2A and p53 as anti-cuproptosis genes can inhibit the cuproptosis of LUAD, SPI1 can decrease the content of CDKN2A and p53 thereby affecting the cuproptosis of LUAD cells. Based on the discovered gene pathway association, SPI1 may promote cuproptosis in LUAD cells by negatively regulating CDKN2A and p53.