Construction of a novel radioresistance-related signature for prediction of prognosis, immune microenvironment and anti-tumour drug sensitivity in non-small cell lung cancer.

Introduction
Lung cancer is the second most common cancer worldwide after breast cancer [1]. It is the malignant neoplasm with the highest mortality rate, the majority subtype of which is non-small cell lung cancer (NSCLC), accounting for 80–85% [2]. NSCLC is composed of multiple histological subtypes, including lung adenocarcinoma, lung squamous cell carcinoma and large-cell lung cancer [3]. Management strategies for NSCLC include local treatment (surgery, radiotherapy, etc.) and systemic treatment (chemotherapy, targeted therapy, immunotherapy, etc.) [4]. Despite significant advancements in treatment strategies, the 5-year survival rate for NSCLC patients is still less than 30%, and an important factor contributing to this situation is treatment resistance, including radioresistance [5].
Radiotherapy represents a crucial local treatment for NSCLC. Previous studies have indicated that upwards of 60% of lung cancer patients present with indications for radiotherapy at a specific stage of the disease, targeting chest diseases or extrathoracic metastases [6]. Over the past century, there have been notable advancements in radiotherapy technology. In recent years, advanced conformal radiotherapy techniques, such as IMRT and proton therapy, have been increasingly applied in clinical practice [7, 8], which may allow for the delivery of efficacious doses in a shortened schedule, reduce radiation damage to adjacent organs, and enhance the therapeutic ratio [9, 10]. However, the radioresistance of tumour cells remains the main obstacle to the application of radiotherapy.
Preclinical studies have demonstrated that mutations in some genes, such as EGFR and ERCC1, lead to radioresistance [11–13]. Additional pathological characteristics and molecular biomarkers, such as DNA mutations, DNA methylation, tumour immune microenvironment and metabolites, also play a role in the effectiveness of radiotherapy [14–16]. Nonetheless, the heterogeneity of the tumour and corresponding tumour microenvironment (TME) limited the sensitivity and specificity of individual genes or biomarkers, which finally induced the failure to predict the prognosis of the tumour [17, 18]. Correspondingly, radiologists have observed that there are individual differences in the efficacy and side effects of radiotherapy in cancer patients in clinical practice. The discrepancy in treatment response suggests that the mining of biomarkers related to radioresistance is particularly important for the treatment of NSCLC [19, 20].
In summary, there are gaps in NSCLC radioresistance that need to be elucidated, and it is still unclear which type of NSCLC patients may benefit more from radiotherapy. The objective of this research was to detect radioresistance-related genes (RRRGs) with prognostic value in NSCLC. Consequently, we conducted a comprehensive investigation of the key gene modules related to the radiation response in the expression profiles of NSCLC patients undergoing radiotherapy, as well as the differentially expressed genes (DEGs) detected in radiosensitive and radioresistant lung cancer cell lines. Additionally, the relationships between gene signatures and immune cell infiltration, as well as the potential underlying mechanisms, were examined to identify specific biomarkers that can predict the sensitivity to radiotherapy. We believe our findings provide perspectives on the mechanisms underlying radioresistance and propose prospective therapeutic targets.

Materials and methods

Study design and project extraction
The flowchart of the present study is depicted in Figure 1. To obtain the necessary data for this study, we used ‘lung cancer’, ‘radiation’ and ‘resistance’ as keywords and searched for the required items in the public functional genomics database Gene Expression Omnibus. We ultimately selected the GSE197236 series as our analytical dataset. In this project, the gene expression profiles of radioresistant A549 lung cancer cells and control cells were described.
The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) was used to extract gene expression profiles and clinical data from NSCLC patients. Specifically, we obtained expression profiles from TCGA for 1153 patients and clinical data for 998 patients, of which 994 matched the expression profiles and clinical data. We excluded 38 samples with missing survival information or overall survival of less than one month, and the final remaining 956 samples were used for subsequent analyses in this study. Seventy-seven patients in this TCGA-NSCLC cohort received radiotherapy and were evaluated for efficacy. Patients with complete or partial response were classified as the radiosensitive (RS) group (n = 33), while patients with progressive or stable disease were classified as the radioresistant (RR) group (n = 44). All TCGA data, GEO data, mRNA expression data, and clinical details were downloaded using R software (version 4.2.2) [21, 22].

Weighted gene coexpression network analysis (WGCNA)
WGCNA serves as a methodology used to cluster genes based on their similar expression patterns, enabling the exploration of the correlation between various clustering modules and clinical phenotypes. Therefore, it is widely utilized to explore the association of clinical phenotypes with gene expression. In this study, we selected patients who received radiotherapy in the TCGA-NSCLC cohort and constructed a gene expression coexpression network targeting the RS and RR groups utilizing the ‘WGCNA’ package in R [23]. The optimal soft threshold was determined to aim at maximizing the adherence of gene interactions to a scale-free distribution. A clustering tree diagram was constructed based on the calculation of gene adjacency and similarity. The dynamic tree cutting algorithm was used to further segment modules and merge similar modules. Each colour represents a module, each containing genes with similar expression patterns. The association of each gene module with sample traits was evaluated by Pearson correlation, and several modules with higher absolute values were selected for subsequent analysis.

Identification and functional enrichment analysis of DEGs
DEGs between radiosensitive and radioresistant cell lines in the GSE197236 dataset were identified using the ‘limma’ package [24]. The screening thresholds were set at |log 2FC| ≥ 1 and p value ≤ 0.05. To visualize the DEGs, the ‘ggplot2’ package was utilized to generate the Volcano plot, while the software package ‘pheatmap’ was utilized for the generation of the heat map [25,26]. The R software package ‘clusterProfiler’ was utilized to carry out enrichment analyses on the DEGs, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) [27].

Construction and validation of gene risk score model
After the radioresistance-related genes (RRRGs) were screened, we merged the expression data with survival data and used the createDataPartition function in the ‘caret’ package to randomly categorize 956 patients from TCGA into training and validation groups at a 7:3 ratio [28]. Subsequently, the ‘glmnet’ package and the ‘randomForestSRC’ package were used for least absolute shrinkage selection operator (LASSO) regression algorithm and random survival forest (RSF) algorithm, respectively, to identify hub RRRGs associated with the survival of patients with NSCLC [29, 30]. The predictive performance of the LASSO and RSF algorithms is evaluated by means of a consistency index (C-index).
In LASSO regression, the RRRGs associated with prognosis with the minimum lambda value were selected. In the case of the RSF algorithm, the 30 most important RRRGs were selected for subsequent analysis. Thereafter, the LASSO regression results intersected with those of the RSF to obtain hub RRRGs. A risk score model was established on the basis of these RRRGs through multivariate Cox regression analysis, utilizing the TCGA training set (n = 670). Using the following formula, the risk score was computed:

In the formula, β represents the corresponding Cox regression coefficient, while X denotes the expression of a specific gene. The median risk score was used to classify patients into low-risk (LS) and high-risk (HS) groups. To compare the overall survival (OS) of patients in the LS and HS groups, Kaplan-Meyer (KM) curves were plotted. The effectiveness of the model was estimated utilizing the receiver operating characteristic (ROC) curves constructed by the software package ‘SurvivalROC’ [31]. The risk score model was validated in the TCGA validation set and the entire TCGA set. Subsequently, to investigate whether risk score was an independent prognostic factor for NSCLC, univariate regression analyses were conducted. A nomogram was utilized to make predictions for the survival rates of patients with NSCLC at 1, 2 and 5 years, and its effectiveness was assessed through calibration curves. When performing multivariate COX regressions, we tested whether the proportional hazards (PH) assumption was met (p > 0.05) to ensure that the factors investigated have a constant effect on the risk over time.

Analysis of immune cell infiltration
Based on RNA-seq data from NSCLC patients, the relative proportions of 22 infiltrating immune cells in different groups were inferred to utilize the CIBERSORT algorithm [32]. Furthermore, the association among immune cells was investigated through Pearson correlation analysis. The immune checkpoint expression levels of different risk score subgroups were compared based on TCGA-NSCLC patient expression profiles. Furthermore, to ascertain the potential for immune escape in tumours between the HS and LS groups, the Tumour Immune Dysfunction and Exclusion (TIDE) algorithm was conducted [33].

Assessment of drug sensitivity
The Cancer Immune Atlas (TCIA) database (https://tcia.at/) was used to obtain the immune phenotype score (IPS). It is a reliable predictor of CTLA-4 (cytotoxic T lymphocyte-associated antigen-4) and PD-1 (programmed death-1) reactivity and is utilized to predict the response of LS and HS groups during immune checkpoint inhibitor treatment [34].
A comprehensive pretrained prediction model was generated using the package ‘oncoPredict’ based on Genomics of Drug Sensitivity in Cancer (GDSC) [35]. The gene expression profiles from the TCGA-NSCLC cohort were utilized to calculate the 50% maximum inhibitory concentration (IC50) for each sample to accurately predict drug reactions [36].

Cell culture
A549, a human NSCLC cell line obtained from the Shanghai Institute for Biological Sciences (Shanghai, China), was cultured in a humidified incubator at 37 °C and 5% carbon dioxide. The culture medium used was Roswell Park Memorial Institute-1640 (RPMI-1640) supplemented with 10% fetal bovine serum (FBS) and 100 U/ml penicillin–streptomycin.

Plasmids and short hairpin RNA (shRNA) transfection
We purchased the overexpression plasmids from GenePhama (Suzhou, China) to specifically overexpress LBH. The Lipo8000TM (Beyotime, Shanghai, China) transfection reagent was utilized to transfect the plasmids into A549 cells. The lentiviruses containing short hairpin RNA targeting LBH (shLBH) and negative control (shNC) were obtained from GenePhama (Suzhou, China).

Quantitative real-time polymerase chain reaction (qRT-PCR)
Total RNA from lung cancer cells was extracted with TRIzol reagent (Invitrogen, USA), which was then reverse transcribed into cDNA utilizing a reverse transcription kit (YEASEN, Shanghai, China). The Hieff® qPCR SYBR Green Master Mix kit (YEASEN, Shanghai, China) was used to conduct the qRT-PCR. The primers’ sequences utilized were as follows: LBH, forward 5′- GCCCCGACTATCTGAGATCG−3′, reverse 5′-GCGGTCAAAATCTGACGGGT−3′. The 2-ΔΔCt method was utilized to evaluate the expression levels of specific genes, with the β-actin gene used as an endogenous control.

Cell irradiation and colony formation assay
A medical linear accelerator (Varian Medical Systems, Inc.) was used to process the cells at a dose rate of 2 Gy/min. Appropriate cell densities were inoculated into six-well culture dishes and incubated until they adhered to the wall. Subsequently, the cells were subjected to different doses of X-rays (0, 2, 4, or 6 Gy). Phosphate buffer solution (PBS), 4% paraformaldehyde, and crystal violet were used for cell washing, fixation, and staining, respectively, after about 12-14 days of the initial inoculation. The colonies (at least 50 cells) were then numbered. The following linear-quadratic equation was used to plot the survival curve:

The experiment was conducted in triplicate.

Cell Counting kit-8 (CCK-8) assay
Cells were planted into 96-well plates (3000 cells/well) and irradiated with 4 Gy X-rays after the cells adhered. At 0, 24, 48, and 72 h after irradiation, 20 μl CCK-8 reagent (Sunview, Shenzhen, China) was added to every well. After incubation for 2 h, a microplate reader (Tecan Infinite 200 M; Männedorf, Switzerland) was used to measure absorbance at a wavelength of 450 nm.

Flow cytometry for apoptosis and cell cycle
To measure apoptosis, we used the Annexin V-FITC/PI cell apoptosis analysis kit (YEASEN, China) on the basis of the manufacturer’s instruction and IDEAS Application v6.0. Software was used to analyse the samples. To detect the cell cycle, we harvested cells at an anticipative time and fixed cells with 75% ethanol overnight at −20 °C. After rehydration with PBS and centrifugation, the fixed cells were stained with 500 µL of DNA stain (Liankebio, Hangzhou, China) for 30 min in the dark. The samples were analysed by Modfit 5.0 software. The number of cells collected for apoptosis and cell cycle analyses in each sample was 10,000 and 15,000, respectively.

Wound healing assay
The wound healing assay is a simple and efficient way to measure the migration ability of cells. Briefly, cells were planted into 6-well plates and scratched with a sterile 200 μl pipette tip, and then the serum-free medium was added to each plate. The scratch area was acquired using a light microscope (Olympus, Tokyo, Japan) at 0 h and 24 h, and analysed by ImageJ-2.1.0 Software.

Statistical analysis
This study utilized R software (version 4.2.2) and its associated software packages to process, statistically analyse and visualize the data. The Student’s t-test and one-way ANOVA were used for two-group comparisons and multiple-group comparisons for normally distributed variables, respectively. The Wilcoxon test/Mann–Whitney U-test and the Kruskal–Wallis test were used for two-group comparisons and multiple-group comparisons for nonnormally distributed variables, respectively. The Pearson correlation analysis and Spearman’s correlation analysis were developed for analysing correlations for normally and nonnormally distributed variables, respectively. *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001. Statistical significance was set to p ≤ 0.05, while ns represented not significance.

Results

Recognition of radioresistance-related gene modules by WGCNA
To investigate the key gene modules significantly associated with radiation response in NSCLC patients, we selected expression profiles from 77 patients receiving radiotherapy from the TCGA-NSCLC cohort for WGCNA. The 77 patients were clustered using R software (Figure 2(A)), after which a soft threshold of β = 6 was determined to be appropriate for establishing a scale-free network (scale-free R2 = 0.90) (Figure 2(B)). Subsequently, analogous gene modules were consolidated to generate a dynamic cut tree (Figure 2(C)). The correlation heatmap between the gene modules and clinical features on the basis of the Pearson correlation coefficient was shown in Figure 2(D). Among the 12 gene modules identified, four key modules exhibiting a high correlation (|cor| ≥ 0.2) with radioresistance were selected for subsequent analysis (Figure 2(C)). The turquoise module comprised 969 coexpressed genes, the pink module 175, the yellow module 444, and the yellow-green module 75 (Supplementary Table 1).

Screening and functional enrichment analysis of DEGs
To detect the genes that play pivotal roles in radioresistance, we selected the GSE197236 dataset as the research object. By analysing the expression profiles of the radiosensitive and radioresistant groups, a total of 1067 DEGs were detected according to the screening thresholds of |log 2FC| ≥ 1 and p value ≤ 0.05, of which 717 genes were up-regulated, and 350 genes were down-regulated (Figure 3(A)). Figure 3(B) illustrates the expression patterns of the DEGs through a heatmap presentation. Supplementary Table 2 contains a complete list of up- and down-regulated genes. Functional enrichment analyses, including GO and KEGG, were carried out to provide further insight into the functions and pathways of the DEGs. GO analysis indicated that the DEGs were significantly enriched in molecular function (MF), including cytokine activity, CXCR chemokine receptor binding, and receptor ligand activity (Figure 3(C)). Cytokine–cytokine receptor interaction, coronavirus disease, and Kaposi sarcoma were the top three enriched pathways in the KEGG analysis (Figure 3(D)), indicating that immune inflammation might play an important role in radioresistance. After intersecting 1067 DEGs and 1663 key module genes, 99 RRRGs were selected for the identification of hub RRRGs associated with the survival of NSCLC patients (Figure 3E, Supplementary Table 3).

Identification of signature RRRGs and construction of an eight-gene prognostic risk score model
The expression data of the 99 genes mentioned above and the survival data of NSCLC patients were merged to identify RRRGs associated with survival using machine learning algorithms. To ascertain the most optimal lambda value in the LASSO regression analysis, a 10-fold cross-validation was performed. To improve the performance of the RSF algorithm, the optimal parameters (ntree = 1000, mtry = 10, nodesize = 50, nsplit = 2) were determined using manual tuning and the ‘tune.rfsrc’ function in the randomForestSRC package [30]. The 16 genes identified by the LASSO algorithm (Figure 4(A,B) and Supplementary Table 4) and the top 30 genes in terms of importance in the RSF algorithm (Figure 4(C,D) and Supplementary Table 5) were retained for subsequent studies. The C-index of the LASSO and RSF algorithms was 0.64 and 0.66, respectively. Finally, the intersection of the two sets yielded 11 genes: APOBEC3B, DOCK4, FAM167A, IER5L, KITLG, LBH, LY6K, RERG, RMDN2, SPP1, and TSPAN2 (Figure 4(E) and Supplementary Table 6).
Having determined that they all met the PH assumption (Supplementary Table 7), the 11 genes were subsequently subjected to multivariate Cox regression and a risk score model consisting of 8 key RRRGs (APOBEC3B, DOCK4, IER5L, LBH, LY6K, RERG, RMDN2, and TSPAN2) was established (Figure 5(A)). The risk score calculated according to the method described above was used to classify the TCGA-NSCLC cohort into two risk groups (the LS group (n = 335) and the HS group (n = 335)), the distribution of risk scores in the NSCLC samples was showed in Figure 5(B). Figure 5(C) illustrates the distinct expression patterns of eight key genes in the two groups. Additionally, patients in the LS group were found to have a prolonged OS compared to those in the HS group (p < 0.001) (Figure 5(D)). To evaluate the effectiveness of the prediction model, we constructed a time-varying ROC curve. The risk score model displayed moderate accuracy, as evidenced by the area under the curve (AUC) of the TCGA training set was 0.676, 0.692, and 0.614 at 1, 3 and 5 years, respectively (Figure 5(E)).
We used the TCGA validation set (the LS group (n = 143), the HS group (n = 143)) (Figure 6(A)) and the TCGA all set (the LS group (n = 478), the HS group (n = 478)) (Figure 7(A)) to further validate the performance of prognostic signature. Figures 6(B) and 7(B) illustrate the distinct expression patterns of eight key genes in different risk groups for the TCGA validation set and the TCGA all set, respectively. Consistently, K-M analysis showed that patients in the LS group had a significantly better OS than those in the HS group (For the TCGA validation set: p = 0.0018, Figure 6C; For the TCGA all set: p < 0.001, Figure 7(C)). Furthermore, AUC for 1-, 3- and 5-year OS were 0.722, 0.673, and 0.623 in the TCGA validation set, respectively (Figure 6(D)), and 0.683, 0.687, and 0.626 in the TCGA all set, respectively (Figure 7(D)).

The risk score can independently predict prognosis in patients with NSCLC
The univariate and multivariate Cox regression analyses were performed to investigate the effectiveness of the risk score and other clinical features in predicting survival in NSCLC patients. Univariate Cox regression analysis demonstrated that both risk score and TNM stage were significantly related to the prognosis of NSCLC patients (Figure 8(A)). Furthermore, the outcomes of the multivariate Cox regression analysis implied that the risk score was an independent prognostic factor for NSCLC patients (HR: 1.958, 95%CI: 1.586–2.416, p < 0.001) (Figure 8(B)). We performed a PH test prior to multivariate COX regression, the results are shown in Supplementary Table 8.
The nomogram based on the combination of the risk score and TNM stage, which can calculate the total score of each patient, was utilized to estimate the survival of NSCLC patients. The total score was inversely correlated to patient survival (Figure 8(C)). As shown in Figure 8(D), the 1-, 2- and 3-year OS rates predicted by the calibration curve were not significantly different from the observed OS rate.

Comparative analysis of the immune microenvironment in the LS and HS groups
Due to the significant impact of the immune microenvironment on the advancement of NSCLC, we sought to elucidate the disparities in the immune landscape between the LS and HS groups [37]. The infiltration of 22 immune cell types was calculated in both groups utilizing the CIBERSORT algorithm. The results demonstrated that in the LS group, there was greater infiltration of memory B cells, immature B cells, resting dendritic cells, resting mast cells, monocytes, resting memory CD4+ T cells, CD8+ T cells and regulatory T cells (Tregs), while in the HS group, there was greater infiltration of M0 macrophages, activated mast cells, neutrophils, and resting NK cells (Figure 9(A)). Figure 9(B) illustrates the correlation among immune cells in NSCLC. To investigate the expression levels of 44 immune checkpoints, a comparative analysis was conducted between the LS and HS groups. The results revealed significant expression disparities between the two groups for 34 immune checkpoint molecules. With the exception of CD44, CD276, and TNFRSF18, the remaining 31 immune checkpoint genes, including CD244 and CTLA4, exhibited increased expression in the LS group (Figure 9(C)).
Furthermore, TIDE was utilized to evaluate the degree of immune exclusion between the LS and HS groups. Higher TIDE scores indicate a greater likelihood of immune escape and less effective ICI treatment. According to the results, the LS group exhibited lower TIDE and T-cell exclusion scores in comparison to the HS group, while the LS group had a higher T-cell dysfunction score than the HS group. When comparing the microsatellite instability (MSI) status of the LS and HS groups, no significant differences were observed (Figure 9(D)).

Assessment of treatment response on the basis of the risk score model
Evidence indicated that the treatment effect of immune checkpoint inhibitor (ICI) therapy can be precisely obtained by IPS, with the higher the IPS score, the better the response to immunotherapy [34]. Therefore, the IPS was calculated for patients in the TCGA-NSCLC cohort, and the results showed that compared to the HS group, the LS group owned a higher likelihood of therapeutic benefits in ICI treatment (Figure 10(A–C)).
The ‘Oncopredict’ tool was employed to predict the response of NSCLC patients to 198 distinct anti-tumour drugs, including both chemotherapeutic drugs (e.g. cisplatin and paclitaxel) and targeted drugs (e.g. gefitinib and axitinib) currently in clinical use. The lower the sensitivity score is, the more sensitive the clinical response, moreover, patients with different risk scores had different responses to 156 drugs (Figure 10(D) and Supplementary Figure 1). In most cases, the HS group exhibited greater sensitivity to these drugs.

Verification of the relationship between RRRGs and radioresistance in NSCLC cell line
To validate the role of key RRRGs in NSCLC radioresistance, we selected the LBH gene, which has the greatest impact on clinical outcomes, based on the hazard ratio values of the multivariate Cox and verified its biological function in vitro experiments. We overexpressed LBH in A549 cells and verified the results by qRT-PCR (Figure 11(A)). After irradiation, the colony formation assay and CCK-8 assay showed that cell proliferation was lower in the LBH-overexpressing group compared to the control group (Figure 11(B–D)), the survival fractions were 31.28% ± 2.94% and 16.22% ± 0.79% for the parental A549 cells and LBH-overexpression A549 cells after 2 Gy X-rays irradiation, respectively (p < 0.01). Furthermore, LBH overexpression exerted a suppressive impact on the malignant phenotype of A549 cells, as evidenced by increased apoptosis, increased G2/M phase arrest, and decreased cell migration, as observed in the apoptosis, cell cycle, and wound healing assays (Figure 11(E–J)).
We also used shRNA to knockdown LBH and verified the silencing efficiency by qRT-PCR (Figure 12(A)). CCK-8 assay showed that the knockdown of LBH significantly enhanced the proliferation of A549 cells, which can also damage the cytotoxic effects caused by radiotherapy (p < 0.0001) (Figure 12(B)). Similarly, the colony formation ability of A549 cells was enhanced with the knockdown of LBH (Figure 12(C,D)). As shown in Figure 12(E–F), flow cytometry demonstrated LBH markedly enhanced the apoptosis induced by radiotherapy (apoptosis rate: X-rays group: 15.57% ± 0.61%, shLBH + X-rays group: 12.66% ± 0.98%, p < 0.05). the G2/M arrest was decreased in the shLBH combined with radiotherapy group compared with the radiotherapy group (Figure 12(G–H)). Finally, we found that the knockdown of LBH significantly promoted the migration ability of A549 cells (percent wound closure: shNC group: 72.37% ± 1.85%, shLBH group: 79.32% ± 0.88%, p < 0.01) (Figure 12(I–J)).

Discussion
Radiotherapy serves as a crucial local treatment for NSCLC patients. Despite the significant progress in radiotherapy technology, radioresistance is inevitable, regardless of the type of radiation (X-ray, proton, or carbon ions) or the dose division method (IMRT, SBRT) used [20, 38–41]. Therefore, it is imperative to fully understand the heterogeneity of NSCLC, identify radiosensitive and radioresistant groups, understand the mechanism of radioresistance, and develop new strategies to overcome radioresistance. These endeavours will contribute to the enhancement of clinical benefits for NSCLC patients.
In the present study, we analysed the expression profiles of radiosensitive and radioresistant cell lines in the GEO database, along with the RNA-seq data of patients receiving radiotherapy in the TCGA-NSCLC cohort. Subsequently, we screened eight key radioresistance-related genes and used them to construct a risk model for predicting the survival of NSCLC patients. In addition, our outcomes revealed that the immune landscape and response to treatment differed between patients in different risk groups, providing a reference for clinical decision-making.
We conducted a risk score model consisting of 8 key RRRGs. Briefly, APOBEC3B is a component of the innate immune system and belongs to the apolipoprotein B mRNA editing enzyme family, functioning as a cytidine deaminase. In numerous tumour types, there is a marked upregulation of APOBEC3B expression, which is strongly implicated in the poor prognosis of cancer patients [42–44]. TSPAN2 (tetraspanin 2) is a member of the tetraspanin family that contains four transmembrane domains of membrane proteins and is present in almost all mammalian cells. Previous research has indicated that high expression of TSPAN2 predicts a poor prognosis for cancer patients [45, 46]. Dedicator of cytokinesis 4 (DOCK4) is a member of the DOCK family. As an agonist of Rho family GTPases, it has biological functions such as nerve cell development and angiogenesis [47, 48]. DOCK4 contributes to tumour progression and metastasis in gastric cancer and breast cancer, but high levels of DOCK4 predict better survival in glioblastoma [49–51]. Many studies have shown that immediate early response 5-like (IER5L), similar to other types of immediate early genes, such as Jun and c-myc, has been extensively documented to exert biological functions in regulating the proliferation and migration of cancer cells and has been associated with poor prognosis in a variety of cancers [52–54]. Human LY6K (lymphocyte antigen 6 complex locus K) is a member of the Ly-6/urokinase type plasminogen activator receptor (uPAR) family. It can increase the expression of matrix metalloproteinases and activate the Raf-1/MEK/ERK signaling pathway, thus promoting the invasion of cancer cells [55, 56]. Ras-related and estrogen-regulated growth inhibitor (RERG) was identified for the first time in breast cancer [57]. RERG is expressed at low levels in a variety of cancers. It can inhibit the malignant phenotype of tumour cells by blocking the activation of the Ras/mitogen-activated protein kinase (MAPK)/ERK signaling pathway and attenuating the proliferation and migration of tumour cells [58–60]. The regulator of microtubule dynamics 2 (RMDN2) gene encodes a microtubule-related protein that is involved in regulating microtubule growth [61]. Andrade et al. have reported that high RMDN2 expression was linked to better clinical outcomes in breast cancer patients [62]. There are few studies on RMDN2 in cancer, and its function and mechanism in cancer need to be further explored.
Located on vertebrate chromosome 2p23.1, the limb-bud and heart development (LBH) is a highly conserved tissue-specific transcription factor that plays an important role in human embryogenesis and disease [63–65]. Deng et al. demonstrated that LBH inhibited cancer cell growth and invasion in lung adenocarcinoma and that low expression of LBH was associated with poor prognosis in lung adenocarcinoma patients [66]. Conversely, other studies have demonstrated that LBH exhibits high expression levels in malignant tumours, including breast cancer, hepatocellular carcinoma, and gastric cancer. Its involvement in boosting cancer cell proliferation and angiogenesis contributes to the initiation and advancement of cancer [67–69]. In conclusion, the role of LBH in cancer inhibition or promotion remains inconclusive. Our findings demonstrated that LBH can suppress the malignant phenotype and enhance the sensitivity to radiation of NSCLC cells in vitro.
Radiation not only kills tumours locally but also exerts antitumour effects by stimulating local or systemic immune responses. Radiation can directly induce an immune response by triggering DNA damage and immunogenic cell death in tumour cells, releasing tumour-associated antigens and damage-associated molecular patterns [70–72]. These molecules facilitate the maturation and activation of antigen-presenting cells and the activation of tumour-specific T cells. Second, radiotherapy can also lead to the occurrence of the abscopal effect, which is a phenomenon of tumour regression in areas not directly irradiated [73, 74]. Radiation-induced tumour cell death leads to the release of new antigens, which boost the migration of effector T cells to nonirradiated areas and exert antitumour effect [75–77]. Furthermore, radiation can increase the expression of MHC-I in tumour cells, increase the infiltration and recognition ability of immune cells, and improve the degree of inflammation in the tumour microenvironment by increasing the levels of the chemokines CXCL10 and CXCL16 in the tumour microenvironment to promote the recruitment of effector T cells to the tumour site [78–80]. Nevertheless, radiation may also exert a certain degree of immunosuppressive effects. In our study, GO and KEGG analysis showed that radioresistance-related genes might influence the immune status of TME. Studies implied that irradiation could result in the death of T lymphocytes in the circulating bloodstream [81], as well as create a hypoxic environment in the irradiated area [82], which can facilitate the aggregation of immunosuppressive cells, such as myeloid-derived suppressor cells (MDSCs) and M2 macrophages, and the production of inhibitory cytokines, such as TGF-β. Concurrently, radiotherapy can induce fibrosis, vascular damage, chronic inflammation, and immunosuppression in the TME, which may lead to radioresistance and thus contribute to tumour survival [82, 83]. Radiation-induced immunosuppression may be one of the reasons why some patients do not respond to radiotherapy.
Our research observed that the risk score established by using key RRRGs was significantly related to immune infiltration and the response to immunotherapy. Specifically, the LS group exhibited heightened immune activity and a more favourable response to immunotherapy. This finding suggests a potential mechanism of mutual sensitization to radioimmunotherapy and has important reference value for the development of personalized radioimmunotherapy strategies for NSCLC patients. Notably, this study demonstrated that patients in the HS group exhibited greater sensitivity to chemotherapy, while patients in the LS group demonstrated greater sensitivity to immunotherapy. These findings have the potential to inform clinical decision-making.
It should be noted that this study has certain limitations which require further improvement. In this study, the ability of the risk score model constructed from key genes to predict and validate the clinical outcome, immune invasion status, and anti-tumour drug response of NSCLC patients depended on TCGA clinical data and in vitro experiments. However, the samples sizes of these clinical datasets are small, which may lead to inherent bias to the interpretation of the outcomes, and the in vitro experiments cannot fully represent the heterogeneity and complexity of tumours in the human body. In vivo experiments, large sample clinical trials, and real-world studies should be conducted in the future, which could further validate our findings. Secondly, we predicted drug sensitivity in NSCLC patients using the ‘OncoPredict’ and ‘TCIA’ tools, which needs to be further validated by designing prospective studies in the future. Finally, the potential mechanisms underlying the key genes have yet to be fully elucidated, which may limit the comprehensibly understanding of our conclusions. Further research is needed to help us better understand the mechanism of key genes in radioresistance through genomic sequencing, tumour biopsy, and other technologies.

Conclusion
In summary, we screened radioresistance-related genes and constructed an 8-gene risk score model that could predict clinical outcomes in NSCLC patients. Furthermore, we constructed a prognostic nomogram with clinicopathological features and the risk score model and established calibration curves. In addition, we evaluated the associations between the risk score and immune infiltration and the response to antitumor therapy in NSCLC patients, providing a reference for clinical decision-making. These key molecules could be potential targets for overcoming radioresistance in NSCLC.

Supplementary Material

Supplemental Material