Abstract

Background. Anoikis, a form of programed cell death, plays a pivotal role in the invasion and metastasis of various tumors, including lung squamous cell carcinoma (LUSC). This study aims to construct a prognostic model for LUSC, leveraging anoikis-related genes (ARGs). Methods. A total of 357 ARGs were extracted from the GeneCards database and Harmonizome portals. Subsequently, ARGs influencing LUSC patients’ prognosis were identified using univariate Cox regression analysis. Unsupervised clustering analysis was carried out utilizing the “consensusplus” R package, and LASSO regression was deployed to craft a risk regression model. The ‘IBOR’ R package quantified the immune cell infiltration abundance. Moreover, the “maftools” R package, paired with the GISTIC online tool, facilitated the assessment of gene copy number variations. Experimental validation was conducted through RT-PCR, evaluating the differential expression of eight pivotal genes, and cellular functional assays discerned the influences of the CHEK2 and SDCBP genes on LUSC cells’ migratory and invasive capabilities. Results. Fifteen survival-associated ARGs delineated three molecular subtypes within the TCGA-LUSC cohort. An eight ARG-based risk prognostic model was constructed, delineating significant survival disparities between high and low-risk groups. Notably, the low-risk group manifested a diminished propensity for immune therapy evasion and gene mutations. A comprehensive nomogram, incorporating risk scores and clinical attributes, was fashioned, exemplifying remarkable predictive acumen. Cellular functional assays substantiated that the modulation of CHEK2 and SDCBP expressions conspicuously influenced the migratory and invasive propensities of LUSC cells. Conclusions. This rigorous study unveils novel anoikis-related biomarkers integral to LUSC prognostication. The meticulously constructed risk prognostic model, underscored by these biomarkers, augurs a potent predictive tool for enhancing the LUSC patient prognosis and therapeutic strategies.

1. Introduction

Lung cancer stands as the most prevalent malignancy globally, with staggering statistics revealing 230,000 lung cancer patients and 15,300 associated fatalities in 2020 [1]. It’s crucial to note that LUSC constitutes 20%–30% of non-small cell lung cancers (NSCLCs), positioning itself as the second most common NSCLC subsequent to lung adenocarcinoma (LUAD). Contrary to LUAD, diagnosing LUSC in its early stages proves to be a formidable challenge, coupled with a scarcity of comprehensive treatment options. There has been noteworthy progress in the realm of LUSC diagnosis and treatment in the recent years. Despite these advancements, the prognosis remains disheartening, with a mere 5-year survival rate of 17.1% [2]. In light of these realities, an exigent call is made for the unveiling of innovative biomarkers pertinent to LUSC. Such a discovery could be monumental in facilitating early diagnosis, honing the precision of prognosis determinations, and tailoring suitable treatment strategies, thereby enhancing the survival outlook for LUSC patients [3].

Anoikis, a term first coined in 1994, is used to describe a unique form of apoptosis triggered by the detachment of normal epithelial cells from the extracellular matrix (ECM) [4]. This physiological process plays a pivotal role in organism development, tissue homeostasis, disease onset, and tumor metastasis [5]. Remarkably, the vital steps of cancer progression and metastasis, such as anchorage-independent growth and epithelial–mesenchymal transition (EMT), have been identified as hallmarks of anoikis-resistance [6]. Anoikis operates as a formidable barrier against tumor metastasis by obstructing the migration of tumor cells from their native ECM to the other organs. It has been substantiated through research that a tumor cell’s resistance to anoikis is instrumental in its survival postdetachment from the ECM, a fundamental requisite for malignant tumors to achieve invasive and metastatic colonization [7, 8]. Consequently, anoikis emerges as a prospective therapeutic target for mitigating tumor aggressiveness [5]. Tumor cells deploy a myriad of strategies to the counteract anoikis. These include leveraging autocrine growth factors or paracrine signaling by adjacent cells to trigger the activation of prosurvival pathways, modulating integrin expression patterns, and employing reactive oxygen species- (ROS-) mediated activation of growth factor receptors to eschew apoptosis, and instigating EMT [9, 10]. An array of pathways orchestrating the regulation of anoikis resistance in lung cancer has been unveiled in various studies. For instance, PLAG1 is known to activate the CamKK2–AMPK signaling pathway, which in turn fosters anoikis resistance in LKB1-deficient lung cancer through the regulation of the GDH1 expression [11]. Additionally, the intricate interplay of the TGFβ1-SH2B3 axis with the JAK2/STAT3 and SHP2/grb2s signaling pathways has been implicated in modulating lung cancer’s anoikis resistance [12].

Recent revelations underscore the significance of fibronectin (FN) as a crucial component underpinning resistance to anoikis following cell detachment in lung cancer, marking it as a potential therapeutic target [13, 14]. In alignment with this, a spectrum of targeted pharmacological agents, such as the third-generation EGFR inhibitor WZ4002 and the TMPRSS4 serine protease inhibitor KRT1853, have been heralded for their capacity to induce anoikis and curb lung cancer metastasis by modulating various signaling pathways [15]. Despite these strides, a discernible gap persists in the literature, particularly concerning the intricate regulatory mechanisms of anoikis-resistance in LUSC and the exploration of anoikisdiagnostic and evaluative potential in the context of LUSC. This study pioneers a comprehensive evaluation, shedding light on the impactful role of anoikis in influencing the prognosis and immunotherapeutic approaches in LUSC patient care.

In this study, we conducted a comprehensive bioinformatics analysis on LUSC samples from the TCGA and GEO databases, based on anoikis-related genes (ARGs). We identified eight potential biomarkers and established a risk prognosis model and related nomogram with superior predictive performance, comparisons were made regarding the survival prognosis, immune infiltration abundance, and somatic mutation frequency of high and low-risk patients. Lastly, we validated the differential expression of eight hub genes through qPCR experiments, and verified the impact of two biomarkers, CHEK2 and SDCBP, on the migration and invasion abilities of LUSC cells through cellular functional experiments. This study offers new insights for the subsequent research related to the occurrence and development of LUSC concerning ARGs.

2. Materials and Methods

2.1. Data Source

RNA-seq transcriptome data, along with associated clinical data, were extracted for analysis in this study. The data comprised 470 LUSC samples (each with a survival time exceeding 30 days) and 49 normal samples, all obtained from the TCGA database (The Cancer Genome Atlas Program (TCGA) - NCI). Additionally, two microarray datasets (GSE73403 and GSE30219), complete with corresponding clinical data, were retrieved from the Gene Expression Omnibus database (Home - GEO - NCBI (nih.gov)). These datasets encompass 69 and 61 LUSC samples, respectively. For the purpose of creating a structured methodology, the TCGA-LUSC cohort was designated as the training cohort. Meanwhile, the GSE73403 and GSE30219 datasets were amalgamated, serving as a validation cohort subsequent to the elimination of batch effects. All extracted data underwent a transformation process; the transcript sequencing data from TCGA, along with the series matrix data from GSE73403 and GSE30219, were converted using log2 (data + 1). A collection of 357 ARGs was meticulously curated from two databases: the GeneCards database (GeneCards—Human Genes | Gene Database | Gene Search) and the Harmonizome database (Harmonizome (maayanlab.cloud)), as documented in Table S1.

2.2. Screening Prognositic ARGs

We utilized the “limma” R package, applying stringent criteria (|log2FC > 1| and a false discovery rate (FDR) < 0.05), to discern differentially expressed ARGs between tumor and normal tissues in the TCGA-LUSC cohort. Pursuing precision, a Cox regression analysis honed in on ARGs with a significant association with the survival prognosis of LUSC patients. Enhancing our analytical rigor, the “maftool” R package was employed, facilitating a nuanced exploration and visualization of the mutational landscapes characterizing prognostic ARGs. Culminating our analysis, we meticulously examined the expression correlations among the prognostic ARGs, achieving a comprehensive visualization utilizing the “Circlize” R package.

2.3. Unsupervised Clustering Analysis

Building upon the identification of prognostic ARGs, an unsupervised clustering analysis was meticulously executed using the “ConsensusClusterPlus” R package. Utilizing the “hc” arithmetic and Pearson distances, this analysis involved 1,000 iterations of repeated sampling, each encompassing specimens (≥80%) derived from the TCGA-LUSC cohort. The strategy involved evaluating up to a maximum of eight clusters (k), with the optimal cluster number (k) being ascertained based on the relative alteration in the area under the cumulative distribution function (CDF) curve. To corroborate the clustering analysis results, principal component analysis (PCA) was employed as a verification measure. Subsequently, a Kaplan–Meier (K–M) analysis was conducted to scrutinize the disparities in survival curves across diverse clusters. The culmination of this rigorous analytical process was marked by the visualization of a heat map depicting gene expression, paired with a comparative overview of the clinical information pertinent to various molecular subtypes, facilitated through the adept use of the R package “ComplexHeatmap”.

2.4. Construction and Validation of an Anoikis-Related Risk Model

In the initial phase of our analytical procedure, the application of LASSO analysis was indispensable, utilized via the R package “glmnet” to meticulously exclude overfitting ARGs. Through an exhaustive process involving 1,000 iterations of tenfold cross-validation, the penalty regularization parameter λ was decisively determined. Following this, we embarked on the construction of a risk model, denoted as the ARG score, employing prognostic ARGs, and corresponding coefficients ascertained from the optimal λ value.

The calculation of the ARG score for each patient was meticulously performed using the equation as follows:where βi and Expi symbolize the risk coefficient and the expression of each ARG, respectively. In a subsequent stage of analysis, samples within the TCGA-LUSC cohort were stratified into high- and low-risk categories based on median risk scores. Employing the –M survival analysis in conjunction with the log-rank test, facilitated through the R package “survminer”, allowed for the precise determination of predictive accuracy concerning overall survival (OS). Additionally, the tool time ROC played a crucial role in computing the area under the receiver operating characteristic (AUC) curve, thereby enabling the prediction of survival rates at 1-, 3-, and 5-year intervals. Through rigorous univariate and multivariate Cox regression analyses, the ARG score emerged as an independent prognostic risk factor. In the concluding phase of our study, a validation process was conducted by harmonizing the GSE73403 and GSE30219 datasets (integrated and labeled as the GSE73403&30219 dataset), subsequent to the elimination of batch effects. This step was crucial in reaffirming the prognostic potency of our meticulously crafted risk model.

2.5. Immune Microenvironment Landscape and Drug Response

We utilized seven algorithms from R packages “IBOR” and “CIBERSORT” R script to evaluate the abundance of immune cell infiltration in the tumor microenvironment (TME) of high- and low-risk patients. And visualized the immune cells with significant differences in the infiltration grades. Subsequently, we further conducted a correlation analysis between the infiltration abundance of 22 types of immune cells in the “CIBERSORT” algorithm, risk scores, and eight hub genes. Lastly, based on the ESTIMATE algorithm, we calculated and compared the ESTIMATE, IMMUNE, Stromal scores, and tumor purity of high and low-risk patients, and evaluated the possibility of immune escape in high and low-risk patients using the “TIDE” online tool (http://tide.dfci.harvard.edu/login/). The drug responsivity of the LUSC samples in the high- and low-risk groups was predicted using the R package OncoPredict [16].

2.6. Functional Enrichment Analysis

We employed the R package “limma” to discern differentially expressed genes (DEGs) among various molecular subtypes and risk groups, adopting the criteria of |logFC > 1| and a p-value < 0.05. Transitioning into gene set enrichment analysis (GSEA), gene set annotation files, name, “c2.cp.kegg.v7.4.symbols.gmt” and “c5.go.bp.v7.4.symbols.gmt,” were retrieved from the MSigDB database. This facilitated a comprehensive investigation into the intrinsic functions of the DEGs. Utilizing the R package “clusterProfiler,” our analysis was meticulously based on the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, ensuring a nuanced exploration of the underlying biological processes and pathways associated with the identified DEGs.

2.7. Construction of a Predictive Nomogram

We commenced our analysis with univariate and multivariate Cox regression analyses, emphasizing the TCGA-LUSC cohort. This initial phase involved a meticulous examination of clinical variables—pathological staging, age, gender, and risk scores—to gauge their impact on the prognostic outcomes. The identification of independent prognostic risk factors allowed for the construction of pertinent nomograms, thereby enhancing the visualization of prognostic data. Moving forward, to assess the predictive accuracy of the formulated nomograms, we utilized calibration and ROC curves. These tools facilitated an accurate evaluation of the 1-, 3-, and 5-year predictive outcomes of the nomograms, bolstering the credibility of our predictive models in forecasting clinical paths. Concluding our analysis, we applied decision curve analysis (DCA), a crucial step to appraise the clinical utility of our developed nomograms. This final stage emphasized the practical significance of our models, showcasing their potential to influence and guide clinical decision-making processes [17].

2.8. Quantitative Real-Time PCR

The human normal lung epithelial cell line BEAS-2B, and the LUSC cell lines SK-MES-1 and NCL-H226, were obtained from Fuheng Biotechnology (Shanghai, China). Each cell line was cultivated under the optimized conditions: BEAS-2B and NCL-H226 in 90% RPMI-1640 medium with 10% FBS, and SK-MES-1 in 90% DMEM, enriched with 10% FBS and 1% P/S (PB180120). For molecular analysis, total RNA was extracted using a Trizol Kit (Gibco, USA), ensuring RNA integrity and purity. The RNA was then carefully reverse transcribed into complementary DNA (cDNA) sequences using a specific cDNA reverse transcription kit (GeneCopoeia, USA). Subsequently, the expression levels of eight essential hub genes were analyzed in the mentioned cell lines. The Bio-Rad IQ5 Real-Time PCR System (Life Medicine, California, USA), used with the All-in-One™ qPCR mix (GeneCopoeia, USA) and specific primers, allowed for precise gene expression detection. The exact forward and reverse primer sequences for each gene are provided in Table 1 for reference.

2.9. Transient Transfection

The SK-MES-1 and NCI-H226 cells were first cultivated in 6-well plates, allowing them to reach a cell density ranging from 40% to 60%. Aligning with the Lipofectamine 2000 transfection reagent guidelines, transfections were meticulously executed. After an interval of 48–72 hr succeeding the transfection, the cells were carefully harvested, paving the way for ensuing experiments.

2.10. Transwell Experiment

Matrigel was first dissolved overnight at 4°C and subsequently diluted in a 1 : 8 ratio with prechilled, serum-free culture medium. A 100-μL aliquot of this mixture was carefully added to a prechilled transwell chamber, then incubated at 37°C for an hour to allow the Matrigel to solidify. Excess fluid was diligently removed, and the chambers were prepared by adding 100 μL of serum-free medium to the upper chamber and 600 μL of 5% serum culture medium to the lower chamber, followed by an overnight incubation at 37°C for equilibration. 1 × 106 cells were meticulously counted, resuspended in 100 μL of serum-free DMEM, and positioned in the upper chamber, while 600 μL of complete medium was introduced to the lower chamber. Following a 24-hr incubation period at 37°C in a 5% CO2 environment, the cells from the upper chamber were expertly removed. The remaining cells were subjected to a process of fixation, washing, and staining. Finally, the migrated cells were thoroughly examined, counted, and photographed for the detailed analysis.

2.11. Cell Scratch Experiment

Uniform horizontal lines were meticulously drawn on the back of a 12-well plate using a marker pen and ruler, preparing it for the experiment. Cells, during their logarithmic growth phase, were carefully adjusted to a concentration of 2 × 106/ mL, inoculated into the plate, and then incubated at 37°C in a 5% CO2 environment. For the scratching process, a sterile 200-µL micropipette tip was precisely used, ensuring alignment with the predrawn lines, promoting consistency. Postscratching, the plate was gently washed with PBS two–three times, facilitating the removal of dislodged cells, before adding the complete culture medium for further incubation. Throughout the incubation, regular observations were made at specified intervals. Photographs were captured to diligently monitor the scratch healing process under an inverted microscope. Quantitative analysis was performed using ImageJ software, where scratch widths were meticulously measured and analyzed to ascertain the migratory capabilities of each cell group. A specialized formula, percentage = (0 hr−other time points)/0 hr, was utilized for precise calculations, followed by a thorough statistical analysis to validate and consolidate the experimental observations.

2.12. Statistical Analysis

Statistical analyses were diligently conducted using R software v4.2.0. Various tests were employed to ensure a thorough and nuanced analysis of the data. Pearson’s test was utilized to effectively assess the correlations between two continuous variables, ensuring a comprehensive evaluation of their relationship. For the comparison of continuous variables between two independent groups, the Student’s t-test was meticulously applied, facilitating a robust assessment of differences. To discern disparities in survival times between two independent groups, the log-rank test was judiciously used. A p-value of less than 0.05 was established as the threshold for statistical significance, ensuring the reliability and validity of the results.

3. Results

3.1. Screening and Investigating Prognostic ARGs Related to LUSC

In our study involving 357 ARGs, 138 were identified as differentially expressed in LUSC tissues—with 60 showing high expression and 78 exhibiting low expression (Figure 1(a)). Following a univariate Cox regression analysis of these ARGs, 15 were singled out as significantly impacting the survival prognosis of LUSC patients (Figure 1(b)). Five of these 15 pivotal ARGs—FADD, CHEK2, NTRK2, TUBB3, and SPINK1—manifested high-expression levels, while the remaining 10 were characterized by low expression (Figure 1(c)). An analysis focusing on copy number variations revealed that SPINK1 possessed the most frequent deletion mutations, in contrast, ITGA3 had the highest amplification mutation frequency (Figure 1(d)).Regarding the overall somatic mutation frequency, HGF emerged as the most mutated among the prognostic-related ARGs (Figure 1(e)). The specific chromosomal mutation locations of these 15 ARGs are depicted in Figure 1(f). A closer look at the expression correlation among the selected ARGs unveiled a notable negative correlation between CHEK2 and NTRK2 and other genes, enhancing our understanding of their relational dynamics (Figure 1(g)).

3.2. Identification of Molecular Subtypes

In the TCGA-LUSC cohort, samples were clustered utilizing 15 prognostic ARGs as a basis. We observed that when k equals 3, the cumulative distribution function (CDF) value diminished more gradually (Figure 2(a)), and the relative alteration of the area beneath the CDF curve was minimized (Figure 2(b)), which facilitated maximal consensus within clusters. Consequently, the samples were categorized into three distinct molecular subtypes: A, B, and C. PPCA revealed significant variations among patients across the three molecular subtypes, illustrating a clear distinction (Figure 2(c)). The K–M analysis further corroborated these findings, demonstrating notable differences in OS across the subtypes, with a p-value of 0.017 (Figure 2(d)). Notably, Subtype A exhibited the most favorable survival prognosis, while Subtype B had the least favorable outcome. Utilizing a heat map, we illustrated the expression of prognostic ARGs alongside the clinicopathological features of samples across the three subtypes (Figure 2(e)). This visual representation suggested that NTRK2 might act as an indicator of positive prognosis. Further exploration through GSEA enabled us to scrutinize the differential enrichment across various pathways, including Hallmark, KEGG, and GO biological processes (BP), particularly between Subtypes A and B (Figure 2(f)2(h)). The GSEA unveiled that gene sets pertinent to pathways such as allograft rejection, complement, interferon response, and immune cells were predominantly expressed in Subtype B. Conversely, pathways linked to cytochrome metabolism and keratinization were predominantly prevalent in Subtype A, indicating a discernible variation in gene expression patterns between the subtypes.

3.3. Immune Landscape in Different Subtypes

The ssGSEA algorithm was employed to meticulously analyze the infiltration patterns of 28 immune cell types in Subtypes A and B samples, as depicted in Figure 3(a). Notably, Subtype B manifested a higher infiltration of most immune cells, with natural killer CD56 cells being the sole exception, when compared to the other two subtypes. Subsequent to this, a comprehensive evaluation involving 45 immune checkpoint genes was undertaken (Figure 3(b)). The results elucidated that, out of these genes, 24 exhibited significant differences in expression across the three subtypes, with nine being predominantly expressed in subtype A. Leveraging the capabilities of the R package ESTIMATE, a calculated analysis was conducted to ascertain the Stromal, Immune, ESTIMATE Scores, and tumor purity across the various sample subtypes. A careful examination of the findings revealed that subtype B conspicuously possessed the highest scores in Stromal, Immune, and ESTIMATE categories. Contrastingly, subtype A emerged with the highest score in tumor purity, as illustrated in Figure 3(c)3(f).

3.4. Establishment and Validation of the ARG Score

A Lasso-penalized Cox regression analysis was performed in the TCGA-LUSC cohort, selecting eight ARGs based on the optimum λ value (Figures 4(a) and 4(b)) to construct a risk model termed the ARG score. The ARG score formula was delineated as follows: Riskscore = 0.006738897 ExpressionPDK4 + 0.117685875ExpressionSDCBP + 0.034158724ExpressionEDA2R + 0.007896589ExpressionSERPINA1-0.195125554 ExpressionCHEK2 + 0.095678261ExpressionFADD + 0.054168418ExpressionITGA3 + 0.009860555ExpressionSPINK1 (Figure 4(c)). Based on the median ARG score, samples in the TCGA cohort were categorized into high- and low-risk groups. The K–M analysis revealed a significant difference in OS between these groups (Figure 4(d)). A scatter-dot plot further underscored the close association between the ARG scores and the patients’ survival statuses (Figure 4(e)). For validation, the GSE73403 and GSE30219 datasets were merged post removal of the batch effects. The predictive accuracy of the ARG score for patient prognosis was assessed within these validation cohorts, revealing consistent findings with the TCGA cohort in both the K-M analysis and scatter-dot plots (Figures 4(g) and 4(h)). The ROC curve affirmed the risk model’s commendable predictive efficacy across both the TCGA and validation cohorts (Figures 4(f) and 4(i)).Visual representations, such as alluvial diagrams and heat maps, were utilized to elucidate differences in molecular subtypes, risk groups, and survival statuses (Figures 5(a) and 5(b)), where notably, CHEK2 exhibited a heightened expression in the low-risk group. In summation, the ARG score demonstrated robust predictive accuracy for determining the prognosis of LUSC patients.

3.5. Functional Enrichment Analysis

DEGs between high- and low-risk groups were filtered using stringent criteria: |log2FC > 1| and FDR < 0.05. Subsequent enrichment analyses, involving GO and KEGG pathways, were conducted to uncover the underlying biological functions of the DEGs. The analysis revealed a significant enrichment in specific pathways. For instance, in the GO pathways, the DEGs prominently featured in ECM organization, and collagen-containing ECM (Figure 5(c) and 5(d)). KEGG pathway analysis further indicated that the DEGs were notably associated with ECM organization and neuroactive ligand–receptor interaction. Moreover, GSEA illustrated that in the high-risk group, DEGs were predominantly enriched in pathways related to EMT, allograft rejection, and the inflammatory response (Figure 5(e)). Conversely, in the low-risk group, there was a significant enrichment in pathways such as E2F targets, MYC targets, and G2M checkpoint (Figure 5(f)). These findings align with the established research [18, 19], corroborating the crucial roles of EMT and ECM in fostering resistance to anoikis in cellular structures, thus underscoring their significance in the underlying molecular mechanisms.

3.6. Immune Infiltration in Different Risk Groups

Prior research has illuminated the complex interplay between tumors and the immune TME, a pivotal arena influencing tumor invasion and the efficacy of immunotherapies [20]. In our study, utilizing the “IOBR” package for immune cell infiltration abundance calculations, a conspicuous increase was observed in high-risk patients compared to their low-risk counterparts (Figure 6(a)). To elaborate, the relative proportions of 22 immune cells infiltrating the high- and low-risk groups were meticulously evaluated using the “CIBERSORT” R script (Figure 6(b)). Among these, M0 macrophages prominently surfaced as the most abundant (Figure 6(c)). A notable divergence in the infiltration of CD8+ T cells, follicular helper T cells, resting memory CD4+ T cells, regulatory T cells (Tregs), gamma-delta T cells, and neutrophils was observed between the two risk groups, underscoring the potential pivotal role of T cells in dictating the prognosis of LUSC patients [21]. Advancing our exploration, we scrutinized the correlation landscapes, revealing intriguing associations between risk scores and immune cell infiltrations (Figure 6(d)). A discernible positive correlation emerged with the infiltration abundance of resting memory CD4+ T cells and risk scores, while an inverse correlation was apparent with follicular helper T cells (Figure 6(e)). In a comparative analysis of ESTIMATE scores between risk factions, high-risk patients manifested elevated Stromal, Immune, and ESTIMATE scores, diverging significantly from the low-risk counterparts (Figure 6(f)). This was inversely mirrored in tumor purity, which was markedly higher in the low-risk group (Figure 6(g)). Employing the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm, known for its precision in evaluating tumor immune escape trajectories [22], we discerned lower TIDE scores in the low-risk group, exhibiting a robust positive correlation with ARG scores (Figure 6(h)). This ensemble of findings collectively paints a picture where low-risk patients exhibit a heightened propensity to derive benefits from the immunotherapy interventions.

3.7. Mutation and Variation Landscape in Different Risk Groups

We explored the gene mutation landscapes in the two risk groups. As shown in Figures 7(a) and 7(c), a higher mutation rate of the 14 genes whose mutations were closely related to the survival of patients with LUSC [23] was observed in the low-risk group than in the high-risk group (96.49 vs. 89.24%), and the mutation rate of TP53 differed significantly between the two groups. However, there were no significant differences in TMB between the high- and low-risk groups (Figure 7(b)), and there was also no significant correlation between the risk score and the TMB (Figure 7(d)). Subsequent to exploring the variations in copy number loci and frequency differences between high and low-risk patient groups, we discovered significant amplification mutations at loci 3q26.33 and 3q26.32 in both groups, whereas at locus 9p21.3, there were notable deletion mutations, with no significant differences in mutation frequency and the number of mutated genes observed (Figures 7(e) and 7(g)). However, in the high-risk group, there were loci, 8p11.23 and 11q13.3, with a higher frequency of copy number variations (Figures 7(a) and 7(h)), possibly indirectly affecting the survival prognosis of the high-risk patients.

3.8. Establishment of a Prognostic Nomogram

In assessing the impact of clinicopathological features and risk score on LUSC prognosis, we integrated the risk score with clinical data, conducting univariate and multivariate Cox regression analyses (Figures 8(a) and 8(b)) using TCGA cohort data. This revealed N stage and risk score as independent prognostic risk factors. Subsequently, a nomogram was constructed utilizing the ARG score model and N stage, aiming to forecast the 1-, 3-, and 5-year survival probabilities (Figure 8(c)). The calibration curve affirmed the prediction accuracy for OS (Figure 8(d)). The time-dependent ROC curve exhibited the nomogram’s robust predictive efficacy (Figure 8(e)), and the DCA curve underscored its substantial clinical utility for 1-, 3-, and 5-year OS projections (Figure 8(f)). These findings suggest that the nomogram, built upon the ARG score and N stage, holds promise for precise prognostic prediction in LUSC, demonstrating potential clinical applicability.

3.9. Quantitative Real-Time PCR

Through quantitative real-time PCR assays, we evaluated the expression of eight hub ARGs in the human normal lung epithelial cell line BEAS-2B, and the LUSC cell lines SK-MES-1 and NCL-H226. We observed that, except for CHEK2, FADD, and SPINK1, these ARGs were more prominently expressed in the BEAS-2B cell line compared to the LUSC cell lines SK-MES-1 and NCL-H226 (Figure 9).

3.10. Functional Roles of SDCBP and CHEK2 in LUSC Cell Lines

We selected two ARGs with the largest impact coefficients, SDCBP and CHEK2, to further study their effects on the migration and invasion abilities of LUSC cell lines. Since SDCBP is lowly expressed in LUSC tissues, and the CHEK2 gene is highly expressed, First, we transfected CHEK2 siRNA into LUSC cell lines NCI-H226 and SK-MES-1, and transfected SDCBP plasmids into the same cell lines, western blot results showed that, compared with the blank control group, the expression of CHEK2 was significantly downregulated (Figures 10(a) and 10(b)), while the expression of SDCBP was significantly upregulated (Figures 10(h) and 10(i)). Next, we conducted cell scratch assays and transwell experiments. The results showed that in both LUSC cell lines, compared with the control group, the scratches continued to heal over time. Knocking down CHEK2 promoted the cells’ ability to migrate into the scratch area (Figures 10(c), 10(e) and 10(f)), and the overexpression of SDCBP also enhanced this ability (Figures 10(j), 10(l), and 10(m)). At the same time, transwell analysis revealed that knocking down CHEK2 or overexpressing SDCBP significantly increased the number of NCI-H226 and SK-MES-1 cells penetrating to the lower chamber (Figures 10(d), 10(g), 10(k), and 10(n)).

3.11. Drug Sensitivity Analysis

Chemotherapy is a critical strategy for the treatment of LUSC. Recent research efforts have focused on the development of more sensitive chemotherapy drugs [24]. OncoPredict is an R package used predict drug response in cancer patients [16]. We calculated drug sensitivity to 198 potential antineoplastic drugs in the low- and high-risk groups. Patients in the high-risk group had a higher sensitivity to AZD5991_1720, Acetalax_1804, Carmustine_1807, Ibrutinib_1799, and MIRA.1_1931 than patients in the low-risk group, whereas the sensitivity to Selumetinib_1736, PRIMA.1MET_1131, SB216763_1025, and Nutlin.3a_1047 was higher in the low-risk group (Figure 11). These results indicated potential sensitivity to drugs for the treatment of patients in different risk groups, and may provide guidance on the choice of chemotherapy drugs for the treatment of patients with LUSC in the high- and low-risk groups.

4. Discussion

LUSC is an extraordinarily aggressive malignancy for which targeted therapies are currently nonexistent [25]. Despite recent advancements in the treatment strategies of LUSC, the prognosis predominantly remains unfavorable [26, 27]. A considerable challenge is the typical late-stage diagnosis, which leads to a median OS of merely 17.1 months for patients undergoing immunotherapy and chemotherapy [28]. Thus, there is a pressing need to unveil innovative diagnostic and therapeutic strategies that promise an enhanced prognosis for LUSC patients [29]. Contemporary research suggests that early intervention in LUSC could be facilitated through predictive models built upon prognostic genes [30]. Various models, incorporating elements like aging-related, immune-related, and pyroptosis-related gene signatures, have been presented. However, the current biomarkers and prognostic models lack comprehensive sufficiency and robustness [3133].

Anoikis significantly influences oncological outcomes by maintaining tissue integrity, curbing abnormal cellular proliferation, and restricting undesirable adherence within an anomalous ECM [34]. Numerous recent studies have meticulously delineated the resistance mechanisms to anoikis in LUSC among other tumors [12, 35]. A specialized gene signature pertinent to anoikis has been elucidated, finding relevance in various tumors including head and neck squamous cell carcinoma, clear cell renal cell carcinoma (ccRCC), and glioblastoma multiforme (GBM) [3638]. These investigative efforts underscore the pivotal prognostic implications of ARGs in oncology, illuminating pathways for the advent of precision therapeutics in cancer management.

This study utilized the RNA-seq transcriptome profile from the TCGA database to identify prognostic ARGs, resulting in an eight-ARG signature and a predictive risk model for forecasting LUSC patients’ prognosis. Additionally, a comprehensive analysis was conducted to explore the correlation between the immune landscape, mutation landscape, immunotherapy response, and chemotherapy response in LUSC patients alongside risk stratification. The findings revealed that prognostic ARGs are instrumental in classifying LUSC patients into distinct molecular subtypes and risk stratifications, showing significant variations in OS, immune infiltration levels, mutation loads, and therapeutic sensitivity.

The potential biological mechanisms of the eight ARGs identified in our study were diverse. Previous studies have indicated that inhibition of glycolysis by CEMIP-mediated downregulation of PDK4 impaired the migration and invasion of prostate cancer cells resistant to anoikis by attenuating the expression of MMP2 and VEGF [39]. Furthermore, PDK4 has been reported to potentiate cisplatin-induced cell death [40]. Our study was the first to provide evidence supporting PDK4 as a potential biomarker for LUSC.

Syntenin/MDA-9 (SDCBP) is prevalently overexpressed in various tumors, marking it a significant therapeutic target to thwart tumor metastasis [41]. Additionally, SDCBP orchestrates protective autophagy and modulates anoikis resistance via BCL-2 phosphorylation in GBM [42]. Kim et al. [43] revealed that SDCBP fosters the migration and endothelial formation of lung cancer cells by managing the secretion of small extracellular vesicles (sEVs) [44]. Knocking down SDCBP expression unveiled that sEVs from syntenin-1-knockdown cells negated the migratory promotion of cancer and endothelial cells, and inhibited endothelial tube formation.

This study corroborated SDCBP’s role in enhancing the migration and invasion of LUSC cells through the construction of overexpressed LUSC cell lines, using diverse approaches. EDA2R has been reported to be a putative tumor suppressor and a direct target of TP53 that prevents tumor progression through regulation of apoptosis and anoikis. EDA2R is often downregulated in the colorectal cancer tissues due to epigenetic alterations or through TP53 gene mutations [45]. We showed that the expression of EDA2R was correlated with a worse prognosis in LUSC patients and was higher in patients without TP53 mutations.

This finding agreed with the results of previous studies and indicated that EDA2R may be associated with the TP53 mutation and resistance to anoikis in LUSC. SERPINA1 has been shown to play an active role in the pathogenesis of NSCLC [46]. CHEK2, a tumor suppressor protein, has been reported to be closely related to the progression of several types of tumors [47] and can induce cell death to regulate anoikis by activating PRAS40 in papillary thyroid cancer (PTC) [48]. Furthermore, CHEK2 has been shown to be an important component in an anoikis-related risk model of ccRCC [37].

Although CHEK2 is highly expressed in LUSC cells, it still plays a role as a tumor suppressor gene in LUSC. Our research has also confirmed this through cellular functional experiments. Additionally, our study showed that CHEK2 was significantly associated with the ESTIMATE score and tumor purity, indicating that future research focusing on the effects of CHEK2 on the immune microenvironment in LUSC may be valuable. FADD plays an important role in cancer progression [49] and has been reported to be a regulator of cell life and death [50]. FADD can block anoikis through the formation of a Fas–FADD complex [51]. ITGA3 can mediate miR-124-associated regulation of tumor cell anoikis. IGTA3 is highly expressed IS NSCLC and has been shown to be associated with invasion and metastasis in LUAD and LUSC [52]. Serine peptidase inhibitor kazal type I (SPINK1) has been shown to promote anoikis resistance in ovarian cancer through a distinct mechanism involving protease inhibition [53].

The immune TME, a crucial determinant in the trajectory of tumor progression and metastasis, substantially influences the efficacy of immunotherapy. In our meticulous analysis, we explored the correlations between the eight pivotal ARGs and the infiltration abundance of 22 immune cell types, yielding insightful findings. Specifically, it was discerned that follicular helper T cells exhibited an inverse relationship with the eight central ARGs, with the exception of CHEK2, FADD, and SPINK1. Conversely, memory CD4+ T cells demonstrated a positive correlation with these ARGs, barring CHEK2, FADD, and SPINK1. Such discoveries underscore the potential significance of concentrating on the relationship between follicular helper T cells and memory resting CD4+ T cells, suggesting that this focus could unveil valuable insights in the realm of immunotherapy efficacy and tumor progression. This nuanced understanding could potentially be instrumental in enhancing therapeutic strategies and interventions in the future.

In the present study, we constructed and validated an eight ARGs signature and a nomogram-based model for patient outcomes. The model had excellent predictive performance in patients with LUSC and may help in clinical decision-making and the development of personalized LUSC treatments. However, studies at the single cell level may improve the predictive power of the model due to significant heterogeneity between the tumor cells. Furthermore, a larger sample size is needed to calibrate the risk model. The addition of in vitro experiments may further strengthen the findings in our study. In summary, our study provided a new perspective and a powerful tool for the early diagnosis and development of precision drugs for the treatment of LUSC.

Data Availability

The datasets TCGA-LUSC, GSE30219, and GSE73403 for this study can be found in the TCGA database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) and GEO database(https://www.ncbi.nlm.nih.gov/geo/), respectively.

Disclosure

Preprint of the manuscript [54]: https://assets.researchsquare.com/files/rs-3121381/v1/3f1c1959-9af5-4744-967c-b1bc9a38eb61.pdf?c=1688916571.

Conflicts 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.

Authors’ Contributions

JL performed the data analysis, qPCR experiments, and drafted the manuscript. HZ completed the cellular function experiments. LW designed this study. All authors contributed to the article and approved the submitted version.

Acknowledgments

I would like to express my gratitude to Hui Zheng for her support and encouragement in my research work. This research was supported by Key Science and Technology Projects in Henan Province (No. 202102310457, No. 212102310669) and Henan Province medical science and technology research plan joint construction project (No. LHGJ20200034), and “23456 Talent Project” of Henan Provincial People’s Hospital.

Supplementary Materials

Table S1: the 357 anoikis-related genes. (Supplementary Materials)