Exploring ribosome biogenesis in lung adenocarcinoma to advance prognostic methods and immunotherapy strategies

Abstract

Background

Lung adenocarcinoma (LUAD) presents a considerable danger to human health and has evolved into a major public health concern. Ribosome biogenesis (RiboSis) is a critical process for synthesizing ribosomes, closely associated with cancer initiation, progression, and treatment resistance, potentially serving as a target for future cancer therapies.

Methods

Utilizing single-cell RNA sequencing (scRNA-seq) technology, a single-cell atlas of LUAD was delineated, focusing on the analysis of T cell subpopulations. Cells were scored based on the expression patterns of 331 genes associated with RiboSis across different cell types, and monocle2 was employed to analyze the developmental trajectory of CD4⁺ T cells. Employing various machine learning algorithms, a ribosome biogenesis-related signature (RBS) was constructed and compared to 140 published LUAD prognostic models. The relationship between RBS risk scores and various factors in LUAD patients, including prognosis, the tumor immune microenvironment (TIME), responsiveness to immunotherapy, and sensitivity to pharmacological treatments was specifically analyzed. Immunohistochemistry was utilized to validate the expression levels of immune markers in the high- and low- RBS groups, and in vitro experiments were performed to validate the functional role of the pivotal gene KIF23 in the progression of LUAD.

Results

Using single-cell analysis, two distinct T cell subtypes were identified: CD8⁺ interferon (IFN) response T cells and CD4⁺ stress response T cells. It was observed that CD4⁺ naive-like T cells exhibit high expression of RiboSis-related genes, with a gradual decrease in RiboSis activity as CD4⁺ T cells develop. Compared to other prognostic models, RBS demonstrated superior performance in prognosis prediction. The low-RBS group exhibited a tumor microenvironment (TME) more favorable for efficient immune monitoring and reaction, higher responsiveness to immunotherapy, and a better prognosis. Immunohistochemistry confirmed higher expression levels of immune markers in the low-RBS group, while in vitro experiments validated the promoting role of KIF23 in LUAD cell proliferation, migration and invasion.

Conclusion

This study delves into the relationship between RiboSis and LUAD cell subpopulations, identifying a potent prognostic biomarker for LUAD. This biomarker aids in assessing immunotherapy efficacy in LUAD patients, ultimately enhancing their prognosis and guiding clinical decision-making.

Introduction

Lung cancer (LC), characterized by its high mortality rate, is a prevalent malignant tumor [1]. It is broadly classified into two main types: small cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC) [2]. Within the NSCLC subtypes, lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) are the most frequent [2]. For most early-stage LUAD patients, early lesion resection through surgery is the most effective treatment [3, 4]. The cure rate for advanced-stage LUAD remains low, with many patients at this stage having lost the opportunity for surgery and can only extend survival through methods such as chemotherapy and targeted therapy [4, 5]. In recent years, immunotherapy, specifically immune checkpoint blockade (ICB) specifically targeting programmed cell death protein 1 (PD-1) and cytotoxic T-lymphocyte associated antigen 4 (CTLA-4) pathways has emerged prominently, which has revolutionized cancer treatment, bringing hope to advanced-stage LUAD patients [6, 7]. However, not all patients respond favorably to immunotherapy [6]. Therefore, identifying specific subsets of LUAD suitable for immunotherapy is a crucial focus in contemporary oncology [8].

Ribosome biogenesis (RiboSis) is a complex and crucial cellular process that initiates in the nucleolus and concludes in the cytoplasm, generating ribosomes essential for cell growth and proliferation [9, 10]. Ribosomes are the mRNA translation machinery responsible for protein synthesis [11]. RiboSis involves multiple steps, including rRNA transcription, rRNA splicing, rRNA modifications, ribosome assembly, and output of ribosomal precursors [9, 10, 12]. Recent studies indicate that ribosomes in tumors exhibit heterogeneity, giving rise to the concept of “cancer ribosomes“ [9]. Increasing evidence suggests that RiboSis is closely associated with cancer initiation, metastasis, and therapy resistance, thereby positioning it as a potential target for the development of innovative cancer treatment approaches [9, 12,13,14]. Existing research has highlighted the crucial role of RiboSis in the progression of NSCLC [15, 16]. Wang et al. discovered that SOD1 in the cell nucleus promotes the proliferation of KRAS-driven NSCLC cells by regulating ribosome biogenesis and proliferation [15]. Lian et al. found that T cell differentiation protein 2 enhances mTOR-mediated ribosome biogenesis to stimulate NSCLC cell proliferation [16]. However, the impact of RiboSis on LUAD metastasis, treatment resistance, and clinical outcomes has yet to be fully elucidated.

In recent years, the continuous advancement of single-cell RNA sequencing (scRNA-seq) technology and bioinformatics has pioneered a new field in the study of LUAD [17]. ScRNA-seq aids in revealing cellular heterogeneity, enabling precise identification and differentiation of diverse cell subpopulations, and discovering specific subgroups that have a vital impact on the emergence and evolution of LUAD [18]. Identifying biomarkers for LUAD helps in evaluating patients’ disease status, predicting disease progression, and assessing treatment responses [19, 20]. Biomarkers play a significant role in personalized medicine, facilitating treatment planning, optimizing therapeutic outcomes, and ultimately augmenting both survival rates and the quality of life for patients with LUAD [19].

In this study, we utilized a scRNA-seq dataset to identify cell subpopulations in LUAD samples and extensively investigated the roles of RiboSis in various cell types. Subsequently, we developed a ribosome biogenesis-related signature (RBS) aimed at improving the prognosis of LUAD patients, assessing immunotherapy efficacy, and providing new insights for the purpose of clinical diagnosis and management of LUAD.

Materials and methods

Data collection and initial processing

The data derived from scRNA-seq for nine primary LUAD samples, were obtained from the Gene Expression Omnibus (GEO) database with the accession number GSE189357, accessible at https://www.ncbi.nlm.nih.gov/geo/ [21]. These data were retrieved specifically for subsequent single-cell analysis. Additionally, bulk RNA sequencing (bulk RNA-seq) data, along with its corresponding clinical information of 1,498 LUAD patients were obtained from various public databases for the purpose of constructing and validating RBS. The sources of these data include The Cancer Genome Atlas (TCGA) database, accessible via its website at https://portal.gdc.cancer.gov/, provides a wealth of data, including bulk RNA-seq data, copy number variation (CNV) data, mutation data, and clinical data of 596 LUAD patients. Furthermore, GEO databases provided additional data, with accession numbers GSE30219 (n = 83) [22]; GSE31210 (n = 246) [23]; GSE42127 (n = 131) [24]; and GSE68465 (n = 442) [25], encompassing RNA-seq data and clinical data of 902 LUAD patients. To facilitate comparative analysis, bulk RNA-seq data of 288 normal lung tissue samples were downloaded from the Genotype-Tissue Expression (GTEx) database (website: https://gtexportal.org/home/). Two immunotherapy cohorts specific to LUAD, namely the OAK dataset (n = 241) [26] and the POPLAR dataset (n = 55) [27], were utilized in this study to assess the predictive capability of the model in forecasting the prognosis of patients undergoing immunotherapy. Additionally, three pan-cancer immunotherapy cohorts, which included the Melanoma-GSE91061 dataset (n = 109) [28], Melanoma-phs000452 dataset (n = 153) [29], and Melanoma-PRJEB23709 dataset (n = 91) [30], were also employed for this purpose. A cohort encompassing 21 distinct types of tumors, constituting a pan-cancer dataset, was utilized to evaluate the broad applicability of the model (accessible at https://xena.ucsc.edu/). In the subsequent analyses conducted, samples that lacked clinical information and with a total survival period (OS) of 0 were excluded. To ensure data comparability, all bulk RNA-seq data were normalized to transcripts per million (TPM) except for differential analysis, which used counts data format. To address potential batch effects, the “combat” function from the “sva” R package was utilized for potential batch effect removal [31]. Prior to analysis, a log2 transformation was applied to all bulk RNA-seq data to ensure that the data were on a more appropriate scale for statistical analysis.

To elucidate the role of RiboSis in the biology and prognosis prediction of LUAD, we constructed a gene set related to RiboSis consisting of 331 genes based on the Gene Ontology (GO) term of the Molecular Signatures Database (MSigDB) and the characterization by Nerurkar et al. (Table S1) [32, 33].

Preliminary analysis of scRNA-seq data

In the comprehensive analysis of the GSE189357 dataset, the scRNA-seq data from nine LUAD samples were processed using the Seurat package, version 5.0.1, within the R software environment, version 4.3.0 [34]. Initially, a rigorous cell quality control procedure was implemented, which involved the systematic retaining of genes and cells of high quality based on a variety of criteria. These criteria encompassed genes that were expressed in a minimum of three cells, the number of genes expressed per cell ranging between 300 and 7000, the expression levels of mitochondrial genes kept below 10%, and the expression levels of red blood cell-related genes maintained under 3%, and cells with less than 100,000 unique molecular identifiers (UMIs). After following the quality control procedures,114,341 high-quality cells were chosen for a more in-depth examination. In the next phase, a series of data preprocessing steps were executed, including normalization, scaling, and the identification of highly variable genes. These tasks were efficiently handled using the specific functions named “NormalizeData”, “ScaleData”, and “FindVariableFeatures”. To further enhance the quality of the data, potential batch effects were addressed by utilizing the “RunHarmony” function, the “RunPCA” function was applied to identify anchors, subsequently employing the “RunUMAP” function for dimensionality reduction [35]. Clustering analysis was executed with the “FindNeighbors” and “FindClusters” functions (resolution = 0.5). Marker genes were pinpointed by contrasting cells within specific clusters against all others clusters using the “FindAllMarkers” function. Annotation of data was done according to representative marker genes for the primary cell types found in LUAD. RiboSis activity scores for individual cells were computed using the “AUCell” and “AddModuleScore” R packages [36, 37]. Cells were then divided into two groups: those with high RiboSis-AUC scores and those with low RiboSis-AUC scores, using the median AUC value as the cutoff. The “cellchat” R package was employed to assess the variations in cell communication between the high and low RiboSis-AUC groups [38].

Construction of T cell atlas and trajectory analysis

From the annotated data, we selected 39,149 T cells and 6,071 NK cells for further investigation. The data processing methods were similar to the previous steps, involving normalization, identification of highly variable genes, batch effect removal, dimensionality reduction, and clustering. When performing clustering using the “FindClusters” function here, the resolution was set to 1.0. Subsequently, cell clusters were annotated conducted on the basis of differentially expressed genes (DEGs) and typical immune markers. RiboSis activity scores were computed for each cell using the “AddModuleScore” R package. Cell communication within T cell subpopulations was assessed using the “cellchat” R package. Furthermore, Monocle2 (version 2.30.0) was employed to infer the cellular lineage trajectory of CD4⁺ T cells [39, 40]. After dimensionality reduction and cell ordering, differentiation trajectories were deduced using Monocle’s predefined settings. Signature genes with q-values < 0.001, calculated using the “differentialGeneTest” function, were used for RiboSis activity scoring of the inferred temporal results. Finally, enrichment analysis of the inferred pseudotime results was conducted using “clusterProfiler” and the “org.Hs.eg.db” R package.

Screening key RiboSis-related genes

Utilizing the “FindMarkers” function, we were able to pinpoint DEGs between the high and low RiboSis-AUC groups, employing the criteria of |log2 fold change (FC)| > 0.7 and p < 0.05. Following this, a Spearman correlation analysis was performed, ultimately selecting the top 150 genes that demonstrated the strongest correlation with RiboSis. By merging these highly correlated genes with DEGs, we obtained a comprehensive list of 830 RiboSis-related genes for further analysis. The bulk RNA-seq data originated from 288 normal lung tissues, sourced from the GTEx database and merged with the TCGA dataset. Normalization was performed applieding a log2(count + 1) transformation and mitigated batch effects. With the normalized and cleaned dataset, we focused our attention on the 830 RiboSis-related genes, resulting in 369 DEGs between LUAD and normal tissues (P < 0.05, |log2FC| > 0.8). Subsequently, univariate Cox regression analysis was carried out using the TCGA-LUAD dataset, leading to the identification of 83 RiboSis-related genes that demonstrated prognostic significance (P < 0.05).

Construction of RBS through combinations of machine learning algorithms

To establish a comprehensive and reliable prognostic signature known as RBS, we integrated 10 machine learning (ML) algorithms: Lasso, elastic net (Enet), stepwise Cox, partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), random survival forest (RSF), CoxBoost, generalized boosted regression modeling (GBM), Ridge, and survival support vector machine (survival-SVM) [41,42,43]. The TCGA-LUAD dataset served as the training set, whereas GSE30219, GSE31210, GSE42127, and GSE68465 with batch effects mitigated, were chosen as validation sets. The outline of RBS construction is as follows: (a) Leveraging a machine learning framework, 81 prognostic models were constructed using 101 algorithm combinations based on key genes. (b) The performance of these models was assessed using four validation datasets. (c) The concordance index (C-index) was computed for each model across all four validation datasets, and the one with the highest average C-index was deemed the optimal model. (d) Further analysis was conducted using the selected optimal model to pinpoint the involved genes, followed by risk score computation for each patient. Subsequently, within each dataset, patients were categorized into high-RBS group and low-RBS group based on the median risk score. This rigorous process guarantees the creation of a validated, reliable, and robust signature across multiple datasets.

Assessment of RBS and construction of a nomogram

Kaplan-Meier survival curves were produced using the “survminer” R package to investigate the disparities in overall survival (OS) and progression-free survival (PFS) between the high and low RBS groups. For evaluating the precision of RBS in predicting LUAD patients’ prognosis, receiver operating characteristic (ROC) analysis was executed with the “timeROC” R package, calculating the area under the curve (AUC). Dimensionality reduction of data was executed via principal component analysis (PCA), visualizing sample distribution differences between groups. A clinical correlation circle plot was generated to explore the relationship between RBS and clinical characteristics. The “CompareC” R package facilitated calculating and comparing the C-index of RBS and clinical features. Following a PubMed literature search, we obtained 140 signatures designed to forecast the prognosis of patients with LUAD. Subsequently, we compared the C-index of the RBS with the C-index of these signatures across multiple datasets. Finally, by combining RBS with TNM staging, a nomogram was constructed to forecast the survival rates for patients with LUAD. To comprehensively evaluate the predictive accuracy and overall performance, we utilized calibration curves, decision curve analysis (DCA), and ROC curves as graphical representation tools.

Mutation landscape

Using GISTIC 2.0 analysis, genomic alterations including recurrent amplifications and deletions were identified. Additionally, the “maftools” R package was utilized to calculate the tumor mutational burden (TMB) [44]. LUAD patients were divided into two groups based on their TMB, using the median TMB score as the cutoff point. Finally, integrating clinical data from TCGA-LUAD dataset, a survival analysis was performed to assess and compare the prognostic outcomes between these two patient groups.

Pathway enrichment

Utilizing the MsigDB database and R software packages such as “GSVA”, “ClusterProfiler”, “org.Hs.eg.db” and “GseaVis”, gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) were conducted [45]. Enrichment analysis of DEGs between the high and low RBS groups was carried out utilizing the comprehensive tools available on the Metascape website. (https://metascape.org/) [46]. Enrichment scores for each sample were calculated based on enriched pathways and were then visualized through the application of t-distributed stochastic neighbor embedding (tSNE) dimensionality reduction.

Immune infiltration

The TIMER2.0 database (accessible at http://timer.comp-genomics.org/timer/) consolidates findings from multiple algorithms and provides a comprehensive summary of immune cell infiltration levels in TCGA datasets [47]. We leveraged this resource to examine variations in immune cell infiltration between the high and low RBS groups. Utilizing the “ESTIMATE” R package, we computed the immune score and stromal score, and ESTIMATE score for each LUAD patient in the TCGA-LUAD dataset, and conducted correlation analysis to assess their association with RBS risk scores. To evaluate potential differences in the Tumor Immune Microenvironment (TIME) between high- and low- RBS groups, a total of 8 unique algorithms were employed, including TIMER, xCell, CIBERSORT, CIBERSORT-ABS, quanTIseq, MCP-counter, EPIC, and ssGSEA [47,48,49], the correlation between the RBS risk score and the cancer-immunity cycle was subsequently investigated. Ultimately, a comparative assessment of the expression of immune checkpoint (IC) molecules and immune modulators was conducted.

Immunotherapeutic response prediction

The relationship between the tumor immune evasion score and the RBS risk score was evaluated utilizing the Tumor Immune Dysfunction and Exclusion (TIDE) website, accessible at http://tide.dfci.harvard.edu/. Additionally, the Immune Phenotype Score (IPS) for patients with LUAD was retrieved from The Cancer Immunome Atlas (TCIA) database (http://tcia.at/home) to evaluate their immunogenicity and predict their response to ICB therapy drugs. We employed the model’s formula to calculate the RBS scores for five immunotherapy cohorts (two LUAD immunotherapy cohorts (OAK, POPLAR) and three pan-cancer cohorts (Melanoma-GSE91061, Melanoma-phs000452, and Melanoma-PRJEB23709)) and conducted survival analysis.

Predicting drug sensitivity and screening small molecule drugs

The data for drug sensitivity were acquired from various reputable sources: the Cancer Therapeutics Response Portal database (CTRP 2.0, accessible via the URL: https://portals.broadinstitute.org/ctrp.v2.1/), the PRISM database (accessible via the URL: https://www.theprismlab.org/), the Cancer Drug Sensitivity Genomics of Cancer Cell Lines (GDSC) database (accessible via the URL: https://www.ancerrxgene.org/). In parallel, the cell line expression profile data were extracted from the Cancer Cell Line Encyclopedia (CCLE) database (accessible via the URL: https://sites.broadinstitute.org/ccle/). A ridge regression model was constructed leveraging the transcriptome data provided by TCGA. By performing a comprehensive analysis, the Spearman correlation between the RBS risk score and the AUC score was calculated (CTRP 2.0: Spearman’s r < -0.4; PRISM: Spearman’s r < -0.55), a total of 24 small molecular drugs were identified. To further understand the sensitivity of high and low RBS groups to common anticancer drugs, Utilizing the GDSC database, drug sensitivity data was obtained and the “pRRophetic” R package was employed to compute the half-maximal inhibitory concentration (IC50) values.

Constructing molecular docking models

Obtain the two-dimensional (2D) structures of SB-743,921 and Ispinesib from the PubChem database (accessible via the URL: http://pubchem.ncbi.nlm.nih.gov/). After successfully acquiring these 2D structures, proceed to input them into the ChemOffice 20.0 software program in order to generate their corresponding three-dimensional (3D) structures. Following this step, the next task involves retrieving the protein target sequences of KIF23 (Uniprotkb: Q02241) and TPX2 (Uniprotkb: Q9ULW0) from the Uniprot database (accessible via the URL: https://www.uniprot.org/), and construct protein crystal models by utilizing the AlphaFold tool, which is available at the web address https://alphafold.ebi.ac.uk. After constructing the protein crystal models, perform energy minimization on the compounds using the Molecular Operating Environment 2019 (MOE 2019) software. Subsequently, preprocess the target proteins to identify active pockets. Conclude this step by performing molecular docking using MOE 2019, with a specified number of iterations set at 50. Evaluating the binding activities of the compounds based on the binding energies and visualize the results using PyMOL 2.6.0 and Discovery Studio 2019 software.

Clinical specimen collection and RNA sequencing

The process of collecting tissue samples has been granted ethical approval from the Medical Ethics Committee of the First Affiliated Hospital of Nanjing Medical University. On the surgery day, 10 LUAD samples, which have been categorized into three distinct groups: AIS, MIA, or IAC by pathology experts, are procured and dispatched to Oncocare Inc. located in Suzhou, China, for the purpose of undergoing RNA sequencing. We collected transcriptome sequencing data from 10 LUAD samples and calculated their respective RBS risk scores. Based on these scores, we classified a sample with the highest RBS scores into the high-RBS group, and a sample with the lowest scores into the low-RBS group.

Immunohistochemistry (IHC)

Paraffin-embedded tissue sections at 37℃ underwent an incubation process for a duration of 120 min. The primary antibodies utilized were anti-CD4 (Catalog number: UMAB64; ZSGB-BIO, Beijing, China), anti-CD8 (Catalog number: SP16; ZSGB-BIO, Beijing, China), and PD-L1 (Catalog number: MXR003; MXB, Fuzhou, China). Subsequent to this primary step, horseradish peroxidase-conjugated secondary antibodies were introduced and allowed to incubate for an additional period of 30 min. The temperature conditions were maintained throughout. Following this, the sections underwent staining with DAB (3,3’-diaminobenzidine) and were subsequently counterstained with haematoxylin to enhance visualization.

Cell lines culture

We acquired the 16HBE human bronchial epithelial cell line, the A549 and the H1299 human LUAD cell line from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences. It was necessary to culture all the cells in DMEM medium (Gibco; Thermo Fisher Scientific, Inc.) containing 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin (Gibco; Thermo Fisher Scientific, Inc.). The culture conditions were maintained to ensure optimal cell growth and viability. This involved setting the incubator temperature at 37 °C, maintaining a CO₂ concentration of 5% and ensuring a humidity level of 95%.

RNA extraction and qRT-PCR

Total RNA was extracted from the cells utilizing the TRIzol^® reagent, which was sourced from Invitrogen; a brand belonging to Thermo Fisher Scientific, Inc. according to the manufacturer’s protocol. Reverse transcription was subsequently carried out with the PrimeScript Reverse Transcription Kit (Takara Bio, Inc.). Quantitative real-time polymerase chain reaction assays (qRT-PCR) was performed using the Archimed-X4 Real-Time PCR System (Applied ARCHIMED X; RocGene Technology, Inc.). The specific PCR conditions employed during this process comprised an initial step involving denaturation, which was conducted at a temperature heated to 95 °C for 30 s during the initial denaturation phase, succeeded by 40 cycles of denaturation at 95 °C for 10 s and annealing at 60 °C for 30 s. Relative mRNA expression levels were quantified employing the 2-ΔΔCt method, utilizing GAPDH as the reference gene for normalization. Triplicates were conducted for each reaction. Below are the primers utilized for the genes investigated in this study:

GAPDH.

Sequence (5’ -> 3’).

Forward Primer TACCCATTTGAATCGTGAGTCA.

Reverse Primer CTCTGGTCCGGTTAGTTCTTTC.

KIF23.

Sequence (5’ -> 3’).

Forward Primer AGTCAGCGAGAGCTAAGACAC.

Reverse Primer GGTTGAGTCTGTAGCCCTCAG.

Transfection of plasmid DNA and small interfering RNA

The KIF23 cDNA was inserted into the pcDNA expression vector. Transfection of plasmids was achieved employing the X-tremeGEN™ HP reagent for DNA transfection sourced from Nanjing, China, whereas small interfering RNA (siRNA) transfection was facilitated by Lipo2000 reagent from Invitrogen, USA, both conducted in accordance with the manufacturers’ protocols. Coverslips positioned in six-well plates were inoculated with A549 cells, and transfection with either plasmid or siRNA was carried out 48 h post-seeding. The specific target sequences utilized are detailed as follows:

si-KIF23-1# Sequence (5’ -> 3’).

Sense GACUAUAUCUAGAUCAUGUCU.

Antisense ACAUGAUCUAGAUAUAGUCUU.

si-KIF23-2# Sequence (5’ -> 3’).

Sense GAAGUGAUCAAUAAUACAACU.

Antisense UUGUAUUAUUGAUCACUUCUA.

Cell counting Kit-8 (CCK‐8) assay

We performed the Cell Counting Kit- 8 (CCK-8) assay utilizing CCK-8 solution sourced from Dojindo Laboratories, Inc., adhering strictly to the manufacturer’s guidelines. Cells were plated in 96-well plates at a density of 2 × 10³ cells per 100 µl. Subsequently, CCK-8 solution was introduced into the wells at various time points: 0, 24, 48, 72, 96, and 120 h. Cell viability was evaluated after an additional 2-hour incubation period, with absorbance utilized a microplate reader to measure at a wavelength of 450 nm.

Colony formation assay

In the colony formation assay, LUAD cells were plated at a density of 1 × 10³ cells per well within six-well plates. These plates were then incubated at a temperature of 37˚C with an atmosphere containing 5% carbon dioxide, for a period of 14 days. During this incubation period, the growth medium was replaced every four days to ensure optimal conditions for colony formation. Once the incubation period was completed, the colonies were fixed using a 4% solution of paraformaldehyde. Subsequently, they were stained with a 0.1% solution of crystal violet, which was sourced from the reputable Beyotime Institute of Biotechnology. Colonies containing more than 50 cells were meticulously counted in each well to ensure accuracy. The experiments were conducted in triplicate.

Wound-healing assay

Once the LUAD cells reached a 90% confluence, they were plated into 6-well dishes. A linear scratch wound was then created across the monolayer using the tip of a 20 µl pipette. After culturing the cells in a medium devoid of FBS, the monolayer was gently rinsed with PBS to eliminate debris and any floating cells. Images of each wound were documented using an inverted microscope immediately after wounding and again after 24 h.

Transwell assay

For the LUAD cell invasion and migration assay, a Transwell chamber (8 μm pore size; Corning, Inc.) featuring 24 pores with a diameter of 8 mm was employed. The matrix gel invasion analysis necessitates the use of Transwell membranes (Sigma Aldrich) that are pre-coated with matrix gel. In each upper Transwell chamber, a total of 4 × 104 cells were carefully cultured and maintained in precisely 300 µl of serum-free medium, ensuring optimal conditions for cellular growth and viability. Simultaneously, the lower chamber was filled with 700 µl of Dulbecco’s Modified Eagle Medium (DMEM), which was further supplemented with 10% fetal bovine serum (FBS). Following a 16-hour incubation period, cells on the upper surface were removed and stained with methanol and 0.1% crystal violet after migrating through the membrane to reach the lower surface. Photographic records were captured utilizing an advanced inverted microscope, specifically manufactured by the reputable company Olympus, in Tokyo, Japan.

Statistical analysis

The initial data processing, coupled with statistical analysis and visualization of results, was carried out utilizing R software (version 4.3.0) and Perl software (version 5.30.0). Both Pearson and Spearman correlation analyses were utilized to evaluate the relationship between two continuous variables. For categorical variables, comparisons were made using the Chi-square test. When dealing with continuous variables, either the Wilcoxon rank-sum test or the t-test was utilized for comparisons. All statistical tests were conducted with a two-faceted approach, with statistical significance ascertained by a P-value < 0.05.

Results

The characterization of RiboSis-related genes in LUAD scRNA-seq data

Figure 1 depicts the workflow of this study. Initially, scRNA-seq data from the GSE189357 dataset underwent cellular quality control, resulting in the selection of 114,341 high-quality cells for subsequent analysis. After mitigating batch effects across samples, the cellular distributions of the 9 samples became relatively consistent (Fig. S1A). Following that, dimensionality reduction was carried out using the “RunUMAP” function, clustering all cells into 24 cell clusters (Fig. 2A). Annotation based on typical marker genes of major LUAD cell types was carried out, resulting in the identification of 9 cell subtypes comprising epithelial cells, endothelial cells, T lymphocytes, NK cells, B lymphocytes, plasma cells, myeloid cells, mast cells, and fibroblasts (Fig. 2B). A bubble plot displayed the typical marker genes of major LUAD cell types (Fig. 2C), a UMAP plot depicted the expression distribution of these marker genes in the data (Fig. 2D), and a proportional plot illustrated the differences in cellular type distributions among the 9 samples (Fig. S1B).

Fig. 1

The flowchart of this study. This figure illustrates the workflow of the study exploring the role of RiboSis in LUAD. The study utilized machine learning to construct RBS for predicting prognosis and immunotherapy response in LUAD patients. Immunohistochemistry experiments confirmed elevated expression levels of immune markers in the low RBS group, and in vitro experiments validated the promoting role of KIF23 in LUAD cell proliferation, migration, and invasion. RiboSis, ribosome biogenesis; LUAD, lung adenocarcinoma; RBS, ribosome biogenesis-related signature; ROC, receiver operating characteristic; TIME, tumor immune microenvironment; qRT-PCR, quantitative real-time polymerase chain reaction assays; CCK-8, cell counting kit-8 experiment

Fig. 2

Characterization of RiboSis-related genes in LUAD scRNA-seq data. (A) Following dimensionality reduction and clustering, cells from 9 LUAD samples were grouped into 24 cell clusters. (B) Cell annotation yielded 9 cell subtypes. (C) Bubble plot displaying typical marker genes for the 9 cell subtypes. (D) Expression distribution of 9 typical marker genes across samples. (E) Violin plots comparing RiboSis activity scores of the 9 cell subtypes. B cells and T cells exhibit higher RiboSis activity scores, while endothelial cells, plasma cells, and epithelial cells show lower RiboSis activity scores. (F) UMAP plot illustrating the cell distribution of high and low RiboSis-AUC groups. (G) Network plot comparing cell communication quantity and intensity between high and low RiboSis-AUC groups. (H) Contrasting incoming and outgoing interaction strengths among different cell subtypes between high and low RiboSis-AUC groups. (I, J) Comparison of signaling pathway strengths between high and low RiboSis-AUC groups. (K) Upregulated and downregulated ligand-receptor pairs in different cell subtypes of high and low RiboSis-AUC groups

Subsequently, based on the expression levels of 331 RiboSis-related genes in the dataset, the RiboSis activity scores for each cell were calculated. The results indicated that B lymphocytes and T lymphocytes exhibited higher RiboSis activity, while endothelial cells, plasma cells, and epithelial cells showed lower RiboSis activity (Fig. 2E, Fig. S1C). Cells were classified into high RiboSis-AUC and low RiboSis-AUC groups based on their average AUC scores (Fig. 2F). Cell communication variances between the two groups were then assessed, the results revealed that the high RiboSis-AUC group exhibited increased quantity and intensity of cell communication (Fig. 2G, Fig. S1D). Additionally, we observed that the high RiboSis-AUC group demonstrated higher incoming and outgoing interaction strengths compared to the low RiboSis-AUC group (Fig. 2H). Furthermore, most signaling pathways showed elevated expression levels in the high RiboSis-AUC group (Fig. 2I, J). Bubble plots illustrated that the high RiboSis-AUC group displayed stronger intercellular communication through a greater number of ligand-receptor pairs (Fig. 2K). Specifically, interactions between myeloid cells and epithelial cells, as well as endothelial cells, were upregulated for ligand-receptor pairs such as RETN-CAP1 and LGALS9-CD44. In summary, these findings suggest a positive correlation between RiboSis activity in cells and the intensity of intercellular communication.

Construction of T cell atlas and evaluation of ribosis activity

Tumor-infiltrating T cells (TILs) are a critical component of TIME and represent a major subset responsive to immunotherapy in LUAD, exhibiting higher RiboSis activity (Fig. 2E). Therefore, we conducted in-depth analysis of T cells and constructed a T cell atlas for LUAD. Following the removal of batch effects between samples, the distribution of T and NK cells across the 9 samples was relatively consistent (Fig. S1E). Through dimensionality reduction and clustering, 20 cell clusters were identified (Fig. 3A). Subsequently, based on DEGs and typical immune markers, we defined 11 T cell subtypes and 1 NK cell subtype (Fig. 3B) [50].

Fig. 3

Construction of T cell atlas and evaluation of RiboSis activity. (A) Dimensionality reduction clustering was performed on T cells and NK cells, identifying 20 cell clusters. (B) Cell annotation defined 11 T cell subtypes and 1 NK cell subtype. (C) Bubble plot illustrating the association between typical marker genes of various T cell subtypes and the 20 cell clusters. (D) Violin plots comparing RiboSis activity scores of 12 T/NK cell subtypes. Higher RiboSis activity was observed in CD4_TN, CD4_TFH, CD8_TISG, and CD4_Treg. (E) Cell communication quantity and intensity among the 12 T/NK cell subtypes. (F-H) Developmental trajectories of CD4 + T cells inferred by Monocle2. (I) RiboSis activity scoring based on pseudo-time results. (J) Expression changes of two hallmark RiboSis-related genes, RPL11 and RPS8, as differentiation progresses. (K) Heatmap displaying expression patterns of genes with the highest differential expression at various stages of CD4 + T cell differentiation

Firstly, we distinguished NK cells (C4, C11, and C19) characterized by high expression of genes such as NKG7, GNLY, and KLRD1. Next, we identified γδT cells (Tgd; C13) showing high expression of TRDC and cytotoxicity-related genes. Subsequently, we differentiated CD8⁺ T cells (high expression of CD8A, CD8B, etc.) and CD4⁺ T cells (high expression of CD4, CD40LG, etc.). For CD8⁺ T cells, based on the characteristic high expression of cytotoxicity-related genes (GZMA, GZMB, GZMH, GZMK, etc.), we defined CD8⁺ effector T cells (CD8_TEFF; CD8_TEFF_GZMK: C1 and CD8_TEFF_GZMH: C8). Furthermore, according to the feature of high expression of exhaustion-related genes (FASLG, PDCD1, LAG3, etc.), we defined CD8⁺ exhausted T cells (CD8_TEX; C16). For CD4⁺ T cells, based on the characteristic high expression of naive genes (CCR7, SELL, TCF7, etc.), we defined CD4⁺ naive-like T cells (CD4_TN; C0, C2, C6, C7, and C17). We also defined CD4⁺ helper T cells (CD4_TFH; C12), characterized by high expression of CXCL13, TOX, ICOS, etc. CD4⁺ regulatory T cells (CD4_Treg; C9 and C10), demonstrating high expression of FOXP3, IL2RA, TNFRSF4, and others. For memory T cells, we defined CD8⁺ resident memory T cells (CD8_TRM; C3 and C15; high expression of ITGAE, ITGA1, ZNF683, etc.) and CD4⁺ central memory T cells (CD4_TCM; C14). Additionally, we identified two distinct T cell subtypes. CD8⁺ interferon (IFN) response T cells (CD8_TISG; C18) exhibited high expression of interferon-stimulated genes (IFIT1, MX1, etc.). CD4⁺ stress response T cells (CD4_TSTR; C5) were characterized by the abundant expression of heat shock proteins and stress response-related genes (such as HSPA1B, JUN, FOS, etc.). Typical immune markers are shown in Fig. 3C.

Subsequently, we calculated the RiboSis activity scores for each T cell subtype. The results indicated higher RiboSis activity in CD4_TN, CD4_TFH, CD8_TISG, and CD4_Treg (Fig. 3D). Within the marker genes of the four cell clusters annotated as CD4_TN (C0, C2, C6, C7, and C17), we identified numerous RiboSis-related genes such as RPS8, RPS6, RPL11, among others. Furthermore, we observed strong intercellular communication between CD4_TN and CD8_TEFF_GZMK, as well as between CD4_TN and CD4_Treg (Fig. 3E). To explore the developmental trajectory of CD4⁺ T cells, we conducted pseudotime analysis. The results revealed two differentiation paths, with CD4_TN as a common starting point. The endpoint of trajectory 1 was CD4_TFH and CD4_Treg, while the endpoint of trajectory 2 was CD4_TSTR (Fig. 3F-H, Fig. S1F), unveiling the sequential activation and exhaustion processes of CD4⁺ T cells in tumor microenvironment (TME). The RiboSis activity scores from the pseudotime analysis indicated a gradual decrease in RiboSis activity as differentiation progressed, with CD4_TN exhibiting higher RiboSis activity scores compared to CD4_TFH, CD4_Treg, and CD4_TSTR (Fig. 3I). The expression levels of two hallmark RiboSis-related genes, RPL11 and RPS8, also decreased as differentiation proceeded (Fig. 3J). A heatmap illustrated the expression patterns of genes showing the highest differential expression at various stages of CD4⁺ T cell differentiation, emphasizing the high expression of RiboSis-related genes (RPS8, RPS11, RPS16, etc.) in the early stages of CD4⁺ T cell development (Fig. 3K). Enrichment analysis results highlighted a significant enrichment of protein synthesis-related pathways (response to unfolded protein, protein refolding, etc.) in the early stages of CD4⁺ T cell development (Fig. S1G).

Developing RBS utilizing integrative machine learning techniques

By integrating DEGs between high and low RiboSis-AUC groups and the top 150 genes most strongly correlated with RiboSis, a total of 830 genes were identified (Fig. 4A). Subsequent differential analysis revealed 369 DEGs between LUAD and normal lung tissues (P < 0.05, |log2FC| > 0.8) (Fig. 4B). Univariate Cox regression analysis identified 83 RiboSis-related genes with prognostic significance (P < 0.05, Table S2). In the model construction process, the TCGA-LUAD dataset was chosen as the training dataset, while GSE30219, GSE31210, GSE42127, and GSE68465 were selected as validation sets, with batch effects removed (Fig. 4C).

Fig. 4

Construction and Evaluation of RBS. (A) Top 150 genes most strongly correlated with RiboSis. (B) Differentially expressed genes (DEGs) between high and low RiboSis-AUC groups (P < 0.05, |log2FC| > 0.7). (C) Sample distribution of five datasets used for model construction before and after batch correction. (D) A total of 81 prediction models were established using a machine learning framework, with the C-index values calculated for each model across all validation datasets. Models constructed using CoxBoost and SuperPC algorithms demonstrated the highest average C-index value (0.69). (E-I) OS-based K-M curves for high and low RBS groups in TCGA-LUAD (P < 0.0001) (E), GSE30219 (P < 0.0001) (F), GSE31210 (P < 0.0001) (G), GSE42127 (P < 0.0001) (H), and GSE68465 (P < 0.0001) (I). (J) PFS-based K-M curves for high and low RBS groups in TCGA-LUAD (P < 0.05). (K-O) Time-dependent ROC curves for RBS in TCGA-LUAD (K), GSE30219 (L), GSE31210 (M), GSE42127 (N), and GSE68465 (O). (P) Patient distribution of high and low RBS groups visualized in the PCA plot of the TCGA-LUAD dataset

To develop RBS, a total of 81 models were constructed using 10 machine learning algorithms. The average C-index value for each model across all validation datasets was computed. The model built using the CoxBoost and SuperPC algorithms demonstrated the highest average C-index value (0.69) (Fig. 4D). This model encompassed 20 key genes, including AKAP12, ANGPTL4, CA4, CAV1, CCDC68, CD69, CD79A, COL4A3, DAAM2, EFNB2, KIF23, LAMB3, METTL7A, MKI67, PTPRC, RPS19, RRAS, SCNN1B, SLC7A11, and TPX2. Using this algorithm, a risk score was calculated for each patient, subsequently categorizing patients into high-RBS and low-RBS groups in each dataset based on the median risk score.

RBS evaluations

In the aforementioned five datasets, survival analysis was conducted for the high- and low-RBS groups, generating Kaplan-Meier survival curves. The results revealed that the OS of the high-RBS group was significantly lower than that of the low-RBS group (P < 0.0001, Fig. 4E-I). In the training dataset (TCGA-LUAD), the high-RBS group exhibited lower PFS rates (P < 0.05, Fig. 4J). Subsequently, time-dependent ROC curves were plotted for RBS in each of the five datasets, and the AUC values were calculated for 1-year, 3-year, and 5-year intervals as follows: 0.7, 0.68, 0.65 (TCGA-LUAD); 0.79, 0.78, 0.76 (GSE30219); 0.68, 0.68, 0.72 (GSE31210); 0.86, 0.74, 0.72 (GSE42127); 0.71, 0.67, 0.62 (GSE68465) (Fig. 4K-O). These findings underscore the robust predictive capabilities of RBS in forecasting the prognosis of LUAD patients. Utilizing PCA analysis for dimensionality reduction and visualization, noticeable disparities in sample distribution exist between the two groups were observed (Fig. 4P, Fig. S2A). Subsequently, we compared the efficacy of RBS and clinical features in predicting the prognosis of LUAD. The risk score of RBS demonstrated a strong correlation with T stage, N stage, and clinical stage in the TCGA-LUAD dataset (P < 0.05, Fig. 5A). We plotted the ROC curves related to clinical features, and the results showed that the AUC value for RBS was significantly higher than that for clinical features (Fig. S2B). RBS exhibited the second-highest C-index value across the five datasets, following only clinical stage (Fig. 5B).

Fig. 5

Evaluation of RBS Prognostic Predictive Ability. (A) Correlation of RBS risk score with age, stage, T stage, N stage, and M stage in the TCGA-LUAD dataset. (B) Comparison of C-index between RBS risk score and clinical features (age, gender, and stage) in the TCGA, GSE30219, GSE31210, GSE42127, and GSE68465 datasets. (C) Calculation and comparison of C-index between RBS and previously published 140 LUAD prognostic models in the TCGA, GSE30219, GSE31210, GSE42127, and GSE68465 datasets

Comparison of RBS and previous signatures

In the past few decades, with the rapid advancement of bioinformatics technologies, an increasing number of prognostic signatures have been established to predict the prognosis of LUAD patients [51, 52]. To assess the superiority of RBS, we conducted an extensive PubMed search and gathered a total of 140 published prognostic signatures for LUAD. Subsequently, by calculating their C-index values in aforementioned five datasets, we compared the predictive performance of RBS with other LUAD prognostic signatures. The results revealed that RBS ranked highest in terms of C-index values among four datasets, while placing sixth in the GSE31210 dataset. This implies that RBS surpasses other signatures in predicting the prognosis of LUAD patients (Fig. 5C).

The development and validation of a nomogram

Given the widespread adoption of the TNM staging system in the diagnosis and treatment of LUAD, we integrated RBS with the TNM staging system to construct a comprehensive and detailed nomogram. This integration aims to significantly enhance the predictive performance [53]. Through the utilization of this newly constructed nomogram, we were able to calculate specific scores for each individual patient. These scores allowed for a more precise prognosis assessment, providing a clearer picture of the patient’s potential outcome (Fig. S3A). The Cox analysis, utilizing this nomogram, pinpointed age, clinical stage, and the RBS risk score as significant prognostic factors (P < 0.05, Fig. S3B). Subsequently, TCGA-LUAD was chosen as the training set and GSE30219 as the test set to assess the precision of this nomogram’s predictions. Calibration curves demonstrated concordance between predicted and observed outcomes (Fig. S3C). The AUC values for 1-year, 3-year, and 5-year intervals in the time-dependent ROC curves were as follows: 0.749, 0.749, 0.709 (TCGA-LUAD); 0.538, 0.739, 0.738 (GSE30219) (Fig. S3C). Finally, DCA indicated the potential of this nomogram to accurately predict survival probabilities of LUAD patients at various time intervals (Fig. S3D). These results suggest that this nomogram exhibits strong prognostic predictive performance.

Genomic alterations

Assessing genomic variations and TMB aids in determining molecular subtypes, predicting patient prognosis, and identifying oncogenic driver genes [54, 55]. Fig. S4A illustrates significant chromosomal alterations between high- and low-RBS groups. In the high-RBS group, genes such as TP53, TTN, CSMD3, MUC16, RYR2, and LRP1B exhibit high mutation rates, with a noticeable increase in chromosomal instability (CIN) and more frequent chromosomal deletions or amplifications (Fig. S4B, C). The high-RBS group demonstrates a higher TMB score, with a positive correlation between RBS risk score and TMB score (Fig. S4D, E). To explore the impact of TMB on prognosis, all LUAD patients were divided into high TMB group (H-TMB) and low TMB group (L-TMB) based on TMB scores. Survival analysis indicates that H-TMB group shows a better prognosis (Fig. S4F, P < 0.05). Subsequently, survival analysis was conducted by combining TMB scores and RBS risk scores. The results show that the H-TMB + low-RBS group exhibits the best prognosis, while the L-TMB + high-RBS group has the poorest prognosis (Fig. S4G, P < 0.001).

Underlying biological mechanisms associated with RBS

The results of GSVA reveal that pathways enriched in the high-RBS group primarily include epithelial-mesenchymal transition (EMT), DNA replication, cell cycle regulation, tricarboxylic acid cycle (TCA cycle), glycolysis, PI3K-AKT-mTOR signaling, and oxidative phosphorylation, among various pro-oncogenic pathways (Fig. 6A, B). Aberrant activation or inhibition of these pathways leads to dysregulated cell cycle control, abnormal DNA replication, and metabolic disturbances, collectively promoting uncontrolled proliferation, genetic instability, enhanced invasion, and metastatic potential of tumor cells, thereby driving rapid progression and deterioration of tumors [56]. This may be a contributing factor to the higher disease risk and poorer prognosis observed in the high-RBS group. Conversely, the low-RBS group demonstrates a favorable prognosis, possibly due to less aggressive tumor biology. This group exhibits lower enrichment levels of pathways associated with unlimited cell proliferation and genetic instability. Metascape analysis indicates that DEGs between the high- and low-RBS groups are enriched in pathways related to cell cycle regulation and DNA replication (Fig. 6C, D). Furthermore, GSEA reveals that major pathways enriched in the high-RBS group include ribosome and ribosome biogenesis in eukaryotes, while the low-RBS group is enriched in cell adhesion molecules pathways (Fig. 6E).

Fig. 6

Underlying Biological Mechanisms Associated with RBS. (A) GSVA analysis describing the biological characteristics of high and low RBS groups. (B) T-SNE plot revealing differences in pathway activity between high and low RBS groups. (C, D) Network plot (C) and bar graph (D) showing enrichment analysis of differentially expressed genes (DEGs) between high and low RBS groups using Metascape. (E) GSEA analysis of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms specifically associated with RBS

The correlation between RBS and immune infiltration

The TIME is closely associated with the occurrence and development of tumors [57, 58]. Tumor immunity plays a critical role in the early stages and metastasis of tumor development, including immune surveillance, immune regulation, and tumor immune escape [59]. In order to acquire a more profound understanding of immune-related characteristics of the high-RBS group and low-RBS group, we investigated their potential biological mechanisms. The results obtained demonstrated that the risk score associated with RBS demonstrated an inverse relationship with the immune score, stromal score, and also the ESTIMATE score. Concurrently, it displayed a positive correlation with tumor purity (Fig. S5A). An in-depth analysis utilizing seven distinct algorithms further revealed a notably higher level of immune cell infiltration within the low-RBS group, which emcompasses infiltration of CD4⁺ and CD8⁺ T cells (Fig. S5B). ssGSEA results demonstrated that the group with low RBS risk scores showed elevated ssGSEA scores across numerous immune cell types (including activated B cells, activated dendritic cells, eosinophils, macrophages, mast cells, neutrophils, natural killer cells, helper T cells, memory T cells, and memory B cells, among others) (Fig. S5C). Furthermore, the low-RBS group also showed higher ssGSEA scores in immune processes such as antigen-presenting cell (APC) co-stimulation, T cell co-stimulation, type II interferon response, and others (Fig. S5C). Additionally, the expression levels of ICs were increased in the low-RBS group (Fig. S5D), suggesting that patients in the low-RBS group are more likely to benefit from ICB therapy. The risk score of RBS was negatively correlated with key steps in the cancer-immunity cycle, including cancer antigen presentation, recruitment of B cells, CD4⁺ and CD8⁺ T cell recruitment, macrophage recruitment, NK cell recruitment, helper T cell recruitment, Treg cell recruitment, T cell infiltration into the tumor stroma, and killing of cancer cells (Fig. S5E). In the low-RBS group, the levels of immune regulatory factors are elevated, such as MHC molecules HLA-DPA1, HLA-DPB1, and costimulatory molecules such as BTNL8, CD27, SLAMF1 (Fig. S5F). This suggests that a lower RBS risk score may help immune cells enter the TME, enhance their anti-tumor activity, which could be beneficial in immunotherapy. In summary, the low-RBS group exhibits higher immune infiltration and immune regulatory factor levels, creating a TME that is more conducive to effective immune surveillance and response, potentially better controlling tumor growth and spread. This enhanced immune activity may be one of the reasons for the better prognosis of the low-RBS group, as a strong presence of immunity within tumors is frequently linked to better treatment outcomes and extended patient lifespan [59].

Response to immunotherapy and drug sensitivity

Immunotherapy represents a significant breakthrough in cancer treatment in recent years, with notable examples including ICB therapy targeting PD-1 and CTLA-4, as well as chimeric antigen receptor T-cell (CAR-T) therapy [6, 60, 61]. For certain advanced LUAD patients, immunotherapy can help improve their prognosis and extend their survival. Therefore, we further explored the potential role of RBS in predicting immunotherapy response. The results showed that the high-RBS group had higher TIDE scores, indicating a stronger tendency towards immune evasion, which may lead to poorer response to immunotherapy (Fig. 7A, p < 0.001). The low-RBS group exhibited higher IPS, suggesting that this group may have a better response to ICB drugs (Fig. 7B, p < 0.001). We extrapolated the RBS risk scores for the immunotherapy cohorts by employing the formula and subsequently identified consistent findings across two distinct LUAD immunotherapy cohorts (OAK and POPLAR, Fig. 7C) and three pan-cancer cohorts (three melanoma cohorts, Fig. 7D). The group with lower RBS risk scores exhibited greater suitability for immunotherapy, showing improved prognoses post-treatment (Fig. 7C, D).

Fig. 7

Response to Immunotherapy and Drug Sensitivity. (A) Difference in TIDE scores between high and low RBS groups. (B) Prediction of IPS in TCGA-LUAD patients using cancer immunogenomic landscape, consistently indicating higher IPS and greater immunotherapy sensitivity in the low RBS group. (C) Calculation of RBS risk scores in two LUAD immunotherapy cohorts (OAK, POPLAR) based on the model formula, followed by survival analysis showing consistent prognostic stratification. (D) Evaluation of RBS in three pan-cancer immunotherapy cohorts (Melanoma-GSE91061, Melanoma-phs000452, and Melanoma-PRJEB23709), followed by survival analysis showing consistent prognostic stratification. (E) Construction of a drug resistance prediction model based on the CTRP2.0 and PRISM datasets, identifying 11 CTRP derivatives and 13 PRISM derivatives. (F) Comparison of sensitivity to four commonly used clinical drugs (Paclitaxel, Osimertinib, Gemcitabine, and Sorafenib) between high and low RBS groups using the PRISM database. (G) Molecular docking models illustrating the interaction between the KIF23 protein and SB-743,921, as well as the interaction between the TPX2 protein and ispinesib

Considering the poor response of the high-RBS group to immunotherapy and the fact that in clinical practice, combination therapy involving chemotherapy and targeted therapy persists the primary approach for treating advanced NSCLC, we targeted the high-RBS group to predict potentially effective small molecule drugs using CTRP2.0 and PRISM databases [4, 5, 62]. The aim was to provide new treatment options and improve their prognosis. Ultimately, 24 small molecule drugs were identified (CTRP2.0: 11 drugs, including SB-743921 and paclitaxel; PRISM: 13 drugs, including ispinesib and docetaxel), with lower AUC values in the high-RBS group (Fig. 7E). Paclitaxel, Osimertinib, Gemcitabine, and Sorafenib are four commonly used anticancer drugs. In the PRISM database, the high-RBS group exhibited higher sensitivity to these four drugs (Fig. 7F, p < 0.001). The molecular docking of two small molecules, SB-743,921 and ispinesib, with the KIF23 and TPX2 proteins (proteins encoded by two RBS modeling genes) was observed (Fig. 7G, Fig. S6A, B). The binding energies from molecular docking were as follows: -7.3895 kcal/mol (KIF23 + SB-743921), -7.1194 kcal/mol (KIF23 + ispinesib); -5.7081 kcal/mol (TPX2 + SB-743921), -5.9550 kcal/mol (TPX2 + ispinesib). Low binding energies signify a high affinity between the compounds and their protein targets. All binding energies were below − 5 kcal/mol, indicating good binding activity between the small molecule drugs and the proteins encoded by the modeling genes. Analysis using the GDSC database revealed significantly higher IC50 values for certain chemotherapy drugs (such as Oxaliplatin, 5-Fluorouracil, Cytarabine, and Cyclophosphamide) and targeted therapy drugs (such as Gefitinib, Crizotinib, and Sorafenib) in the high-RBS group. This suggests that the high-RBS group may respond better to these drugs, potentially compensating for their poor response to immunotherapy (Fig. S6C).

Immunohistochemistry and in vitro experiments

Previous studies demonstrated CD8 and CD4’s correlation with CD8⁺ T cells and CD4⁺ T cells respectively, both capable of anti-tumor activity and predicting favorable immunotherapeutic outcomes [63, 64]. Additionally, heightened expression of the immune checkpoint protein programmed death-ligand 1 (PD-L1) may signify improved immunotherapy efficacy [65]. We collected transcriptome sequencing data from 10 LUAD samples and calculated their respective RBS risk scores. Based on these scores, we classified a sample with the highest RBS scores into the high-RBS group, and a sample with the lowest scores into the low-RBS group. Immunohistochemical analysis revealed that, in contrast to the high RBS group, the low RBS group displayed elevated expressions of CD8, CD4, and PD-L1 (Fig. 8A).

Fig. 8

The Effect of KIF23 on Proliferation, Migration, and Invasion of Lung Adenocarcinoma (LUAD) Cells. (A) Validation of the expression differences of CD8, CD4, and PD-L1 between high and low RBS groups through immunohistochemical staining experiments. (B) Differential expression of 20 model genes at the pan-cancer level in tumor and normal tissues. (C) Correlation analysis demonstrating the relationship between KIF23 and RBS risk score. (D) Survival analysis showing a significant negative correlation between the expression level of KIF23 and the prognosis of LUAD patients. (E) Evaluation of KIF23 expression using qRT-PCR in LUAD cell lines (including A549 and H1299) and human bronchial epithelial cell line 16HBE. (F) Regulation of KIF23 expression in A549 cells using specific siRNA and overexpression plasmids. (G-J) Determination of the effects of KIF23 knockdown and overexpression on proliferation of A549 cells using CCK-8 assay (G) and colony formation assay (H). Exploration of the effects of KIF23 knockdown and overexpression on migration and invasion of A549 cells through scratch healing assay (I) and Transwell assay (J)

The RBS consists of 20 model genes. Utilizing transcriptome data from TCGA and GTEx databases, we compared the differential expression of these 20 model genes between LUAD and normal lung tissues (Fig. S7A). Using pan-cancer datasets, we analyzed the differential expression of the 20 model genes in tumor and normal tissues at a pan-cancer level (Fig. 8B). IHC images obtained from the Human Protein Atlas (HPA) database further validated the differential expression of these genes at the protein level (Fig. S7B). We evaluated the correlation between the 20 model genes and the RBS using ROC curves and correlation analysis, where KIF23 exhibited a significant positive correlation with RBS (correlation = 0.9, AUC = 0.958, Fig. S7C, Fig. 8C). KIF23, a member of the kinesin family, involved in cell division regulation, is highly expressed in lung cancer, associated with lung cancer prognosis, and may serve as a novel target for lung cancer therapy [66]. Survival analysis indicated a significant negative correlation between KIF23 expression and prognosis of LUAD patients (Fig. 8D). The expression of KIF23 was assessed using qRT-PCR in LUAD cell lines (including A549 and H1299) and human bronchial epithelial cell line 16HBE. Results indicated upregulated expression of KIF23 in LUAD cell lines (Fig. 8E). To investigate the functional role of KIF23 in the development of LUAD, we altered the expression levels of KIF23 in A549 cells by employing specific siRNA for knockdown and overexpression plasmids (Fig. 8F). Results from CCK-8 assays and colony formation experiments suggested that KIF23 may promote the proliferation of LUAD cells (Fig. 8G, H). Subsequently, we investigated the potential impact of KIF23 on the migration and invasion of LUAD cells. Wound healing experiments showed that knocking down KIF23 significantly hindered wound closure, while overexpression of KIF23 notably promoted wound healing (Fig. 8I). The results of Transwell assays revealed that the silencing of KIF23 impeded the migration and invasion of LUAD cells, whereas its overexpression markedly boosted their migratory and invasive capabilities (Fig. 8J). These observations validate KIF23’s role in stimulating the proliferation, invasion, and migration of A549 cells, suggesting its potential in driving LUAD progression.

Discussion

Lung cancer poses a significant threat to human health, with a persistently high mortality rate, making it a major public health concern [4, 67]. Identifying reliable and effective biomarkers for lung cancer has become an urgent task. Although various lung cancer biomarkers have been discovered using bioinformatics methods, such as the LUAD prognostic model based on B cell marker genes established by Song et al., the reliability and clinical translational potential of these prognostic models remain to be thoroughly investigated due to the lack of model comparisons and clinical validations [68]. In recent years, scRNA-seq technology has been widely employed in cancer research, with researchers aiming to explore the heterogeneity of tumor cells, map pan-cancer single-cell atlases, identify specific cell subpopulations, and investigate their functions within TME [17]. For instance, Tang et al. constructed a pan-cancer single-cell panorama of human natural killer cells and elucidated specific NK cell subtypes [69]. However, these studies tend to focus more on elucidating biological mechanisms and less on exploring the relationship between target cell subpopulations and cancer prognosis and treatment outcomes. RiboSis is a core step in the process of tumor initiation and progression, with considerable evidence suggesting a close association between RiboSis and cancer occurrence, metastasis, and treatment resistance, holding potential clinical translational value as a target for future cancer therapies [9, 12, 14]. Some studies have investigated the role of RiboSis in proliferation and metastasis in LUAD using laboratory techniques [15, 16]. However, there is currently no literature elucidating whether RiboSis-related genes can serve as biomarkers for LUAD and play a role in prognosis prediction and treatment efficacy assessment.

In this study, we first utilized LUAD scRNA-seq data to construct a single-cell atlas, focusing on describing the cell subpopulations of T cells. We identified two relatively novel T cell subtypes: CD8⁺ interferon (IFN) response T cells and CD4⁺ stress response T cells. Both are associated with stress responses and play distinct roles within the TIME. Subsequently, we delved into the association between cell subpopulations and RiboSis, revealing that CD4⁺ naive-like T cells exhibit high expression of RiboSis-related genes, with RiboSis activity gradually decreasing as development progresses. This discovery paves the way for investigating naive-like T cells and RiboSis in LUAD. Building upon LUAD bulk RNA-seq data and RiboSis-related genes, we constructed a prognostic model (RBS). Comparing it with 140 previously published LUAD prognostic models, we confirmed the superior predictive performance of RBS in prognosis forecasting.

RBS comprises 20 key genes, including AKAP12, ANGPTL4, CA4, CAV1, CCDC68, CD69, CD79A, COL4A3, DAAM2, EFNB2, KIF23, LAMB3, METTL7A, MKI67, PTPRC, RPS19, RRAS, SCNN1B, SLC7A11, and TPX2. AKAP12, a type A kinase anchor protein scaffold, loss of which increases cancer susceptibility, has been reported in studies to be improved in LUAD by miR-338-3p inhibition of AKAP12 [70]. CD69 serves as a leukocyte and NK cell activation marker, playing a crucial role in the pathogenesis of tissue damage under activation of different leukocyte subpopulations and various inflammatory states, with previous research confirming its close association with lung cancer [71]. The gene encoding ribosomal protein S19 (RPS19) is one of the pathogenic genes for Diamond-Blackfan anemia (DBA), with no reported function in lung cancer to date [72, 73]. Targeting protein for Xklp2 (TPX2) is associated with the spindle pole localization of the dynamic protein-dependent Xklp2. TPX2 is highly expressed in various cancers, including lung cancer, with studies confirming its regulation of macrophage polarization in LUAD through the NF-κB pathway [74].

After validation, RBS has demonstrated strong predictive abilities for the prognosis of LUAD patients. The high-RBS group is associated with multiple pro-oncogenic pathways, showing fewer immune cell infiltrates and a lower response to immunotherapy, potentially contributing to higher disease risk and poorer prognosis. Conversely, the low-RBS group exhibits a TME more conducive to effective immune surveillance and response, displaying a higher response to immunotherapy and better prognosis. At the cellular and molecular biology level, we selected one of the modeling genes, KIF23, for in vitro experiments, confirming its role in promoting proliferation and invasion of LUAD cells. This provides foundational biological evidence for the clinical translation of RBS. KIF23, a member of the kinesin motor protein family involved in cell division regulation, is highly expressed in lung cancer and associated with lung cancer prognosis, potentially serving as a novel target for lung cancer treatment [66]. Previous studies have shown that RNA interference-mediated depletion of KIF23 inhibits in vivo lung tumor formation and induces apoptosis in lung cancer cell lines [75].

The RBS developed in this study is highly associated with the prognosis of LUAD patients, demonstrating superior performance in prognostic prediction and exhibiting significant immune relevance, highlighting its potential as a biomarker in clinical practice. Immunotherapy has been a major breakthrough in cancer treatment in recent years, improving the prognosis and extending the survival of certain advanced LUAD patients [6, 7]. RBS can predict the differential therapeutic responses to immunotherapy among LUAD patients, thereby facilitating the identification of LUAD patient subgroups that are sensitive to immunotherapy. Integrating RBS into future LUAD clinical practice may facilitate personalized risk assessment and immunotherapeutic intervention, ultimately improving patient outcomes and guiding clinical decision-making. However, this study has several limitations. First, it is a retrospective study based on public databases. The limited sample size in the validation datasets may restrict the generalizability of the findings, necessitating further independent validation through robust prospective multicenter cohort studies. Additionally, more LUAD immunotherapy datasets should be incorporated to comprehensively validate the predictive power of RBS in assessing immunotherapy response. Second, the single-cell analysis was not comprehensive, as it primarily focused on T cell subpopulations. Finally, we only conducted in vitro validation experiments for KIF23, leaving the exploration of the biological functions of RiboSis-related genes in LUAD incomplete. Future studies should aim to include a more diverse patient population and validate our findings across multiple independent cohorts to ensure their robustness and applicability. Furthermore, additional experimental investigations, including in vivo studies, are required to elucidate the precise mechanisms by which KIF23 influences LUAD cell proliferation and migration, as well as its downstream signaling pathways.

Conclusion

This study employed scRNA-seq technology to construct a single-cell atlas of LUAD and identified that CD4⁺ naive-like T cells exhibit high expression of RiboSis-related genes. Notably, RiboSis activity gradually decreased as CD4⁺ T cells underwent differentiation and maturation. RiboSis, a fundamental process for ribosome synthesis, is closely associated with tumor initiation, progression, and therapeutic resistance, and holds promise as a potential target for future cancer treatments. Accordingly, we developed a RiboSis-based biomarker—referred to as RBS—by integrating multiple machine learning algorithms, including Lasso, SuperPC, and CoxBoost. Through a comprehensive multi-omics analytical framework, we systematically investigated the relationship between the RBS risk score and patient prognosis, TIME, response to immunotherapy, and drug sensitivity in LUAD. Compared to 140 existing prognostic biomarkers, RBS demonstrated superior performance in predicting clinical outcomes in LUAD patients. Moreover, RBS exhibited significant immune relevance. The low-RBS subgroup was characterized by a TME more favorable for effective immune surveillance and response, enhanced responsiveness to immunotherapy, and improved overall prognosis. The ability of RBS to stratify LUAD patients according to differential responses to immunotherapy facilitates the identification of immunotherapy-sensitive subpopulations. Incorporating RBS into future clinical practice may enable personalized risk assessment and immune-based therapeutic strategies, ultimately improving patient outcomes and guiding clinical decision-making. As a core gene within the RBS model, KIF23 was experimentally validated in vitro to promote LUAD cell proliferation, migration, and invasion, supporting its role as a key oncogene in LUAD and highlighting its potential as a therapeutic target. We believe the findings of this study will support future research into the biological mechanisms and clinical translational value of RiboSis in LUAD, offering novel insights and guidance for clinical decision-making.