What this is
- Lung adenocarcinoma (LUAD) is a leading cause of cancer-related deaths, with a poor 5-year survival rate.
- This study identifies exosome-related (ER-) that correlate with LUAD prognosis and treatment response.
- A based on ER- was constructed, revealing significant differences in patient outcomes and immune characteristics.
Essence
- A based on exosome-related (ER-) predicts outcomes in lung adenocarcinoma patients. Higher ER-lncRNA-related risk scores correlate with poorer survival and greater sensitivity to certain chemotherapy drugs.
Key takeaways
- 134 differentially expressed ER- were identified, with 19 linked to LUAD prognosis. These ER- stratified patients into two subtypes, one associated with worse outcomes.
- A higher ER-lncRNA-related risk score (ERS) correlates with poor overall survival across multiple cohorts. The model demonstrated superior predictive power compared to traditional clinical factors.
- Differences in immune landscape were noted between high- and low-ERS groups, with the high-ERS group showing greater sensitivity to most chemotherapy drugs, except rapamycin and erlotinib.
Caveats
- The study relies on data from online databases, which may introduce biases. Further validation with additional data samples is necessary.
- The analysis is based on tissue samples, and the presence of ER- in or blood needs further exploration.
Definitions
- Exosomes: Tiny extracellular vesicles involved in intercellular communication, transferring nucleic acids and proteins.
- Long non-coding RNAs (lncRNAs): RNA molecules longer than 200 nucleotides that regulate gene expression without encoding proteins.
- Prognostic model: A statistical tool used to predict patient outcomes based on specific biomarkers or clinical features.
AI simplified
Introduction
Lung cancer is the most prevalent cancer and the major cancer-related cause of death [1], with lung adenocarcinoma (LUAD) being the most common histological type, accounting for 40% of lung cancers [2]. The incidence and mortality of LUAD are increasing year by year. In recent decades, great advances have been made in treatment strategies, such as surgery, chemotherapy, immunotherapy, radiotherapy or targeted therapy, but the 5-year survival rate for LUAD remains less than 20% [3–5]. In addition, due to the heterogeneity of tumors, criteria for identification of high-ERS patients and selection of precise treatment options are still lacking. Therefore, the exploration of robust predictor is crucial for better individualized precision medicine.
Exosomes, tiny extracellular vesicles of 30–100 nm, have been shown to play a pivotal function in intercellular communication through cell-to-cell transfer of nucleic acids and proteins, participating in a wide range of physiological processes [6, 7]. Recently, studies have revealed that exosomes also have a critical effect in tumorigenesis and progression [8, 9]. In particular, the role of exosomes in LUAD has been preliminarily elucidated, involving in LUAD progression [10, 11] and metastasis [12], drug resistance [13].
Long non-coding RNAs (lncRNAs), no less than 200 nucleotides in length, can regulate gene expression levels but are not involved in encoding proteins and play a regulatory role in physiology and disease development [14]. Evidence has suggested a potential role for lncRNAs in the onset and progression of LUAD [15]. Recent studies have revealed that lncRNAs can be encapsulated into exosomes and widely participate in intercellular signaling, which has important pathophysiological significance [16]. The exosomal lncRNAs can regulate antigen presentation, change the tumor microenvironment (TME) and participate in critical oncological behaviors, including tumor growth, metastasis, invasion [17–19]. In addition, studies have demonstrated that the dysregulation of exosomal lncRNAs are participating in multiple processes of cancer progression, making them diagnostic biomarkers and potential targets for cancer therapy [20].
However, there are few studies that focused on the role of exosomal lncRNAs in LUAD. In this study, we recognized the prognostic exosome-related lncRNAs (ER-lncRNAs) in LUAD and constructed an ER-lncRNAs signature to predict the prognosis and therapeutic response of LUAD patients. Our study may provide new directions for prognostic assessment and personalized therapy strategies for patients with LUAD.
Materials and methods
Data collection and preprocessing
Gene expression data (raw counts) and clinical data of the Cancer Genome Atlas (TCGA)-LUAD cohort were obtained from TCGA database (https://portal.gdc.cancer.gov/repository↗). Raw counts were transformed into transcripts per million (TPM) values for subsequent analysis. Gene expression data and clinical data of five independent Gene Expression Omnibus (GEO) datasets, namely, GSE30219↗ (n = 83), GSE3141↗ (n = 58), GSE50081↗ (n = 127), GSE31210↗ (n = 226) and GSE72094↗ (n = 398) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/↗). The ComBat algorithm in the “sva” R package was used for eliminating the batch effects between the TCGA-LUAD and GEO cohorts [21, 22]. The TCGA-LUAD, GSE30219↗ (n = 83), GSE3141↗ (n = 58), and GSE50081↗ (n = 127) cohorts were consolidated into a Meta cohort. Then, the gene expression data of Meta cohort were normalized using scale function of "limma" R package. GSE31210↗ (n = 226) and GSE72094↗ (n = 398) were used as the external validation cohorts. Besides, 121 exosome-related mRNAs (ER-mRNAs) [23] were collected from ExoBCD database (https://exobcd.liumwei.org/↗).
ER-lncRNAs
112 ER-mRNAs were matched in the Meta cohort, which were used to identify ER-lncRNAs. Differentially expressed (DE) analyses of ER-mRNAs and lncRNAs were performed using the “DESeq2” R package based on raw counts data. |logFC| >1 and FDR <0.05 were set as the screening threshold. The correlations between the expression levels of DE-ER-mRNAs and DE-lncRNAs by Spearman correlation analysis, and |cor| > 0.3 and p < 0.05 were screened as ER-lncRNAs.
Consensus clustering
Prognostic ER-lncRNAs were screened using univariate cox regression analysis by “survival” R package (p < 0.05). Subsequently, based on the expression of the prognostic ER-lncRNAs, unsupervised consensus clustering was performed to clustering patients by “ConsensusClusterPlus” R package [24]. Survival curves of the overall survival (OS) between different clusters were constructed by the Kaplan–Meier analysis. The principal component analysis (PCA) was applied to verify its validity.
Gene set enrichment analysis (GSEA)
GSEA was performed to assess biological progress between distinct groups using the GSEA software (3.0). Gene expression data in the Meta cohort were used for GSEA analysis. The reference gene set of c2.cp.kegg.v7.0.symbols.gmt or h.all.v7.4.symbols.gmt was obtained from Molecular Signatures Database [25] (http://www.gsea-msigdb.org/gsea/downloads.jsp↗). The statistically significant pathways were screened with p < 0.05.
Assessment of the immune landscape
The activity of immune cell, function, and pathway of each LUAD sample were calculated with a single sample gene set enrichment analysis (ssGSEA) algorithm using “GSVA” R package [26]. The “estimate” R package was used to assess the immune score, stromal score of each LUAD sample [27].
Construction and validation of an ER-lncRNA-related prognostic model
The differentially expressed genes (DEGs) between distinct clusters were recognized using the “limma” R package with |logFC|> 1 and a false discovery rate (FDR) < 0.05. The prognostic DEGs were identifed using univariable cox regression (p value < 0.05), which were extracted for constructed an ER-lncRNA-related prognostic model using least absolute shrinkage and selection operator (LASSO) Cox regression analysis by “glmnet” R package [28, 29]. Based on the constructed prognostic model, the risk score for each LUAD sample in Meta cohort, which we termed as ERS, was calculated as follows: ERS = n∑Coefi × Expi, with n, Coefi and Expi representing the number of model genes, gene i beta coefficient value and gene expression level, respectively. Patients were divided into high- and low- ERS groups by the median of the ERS. Kaplan–Meier curve and time‐dependent receiver operating characteristic (ROC) curve by the “survival” and “timeROC” R package were used to assess the OS difference and predictive accuracy of the ERS, respectively. The univariate and multivariate Cox regression analyses were applied to evaluate the independent prognostic value of the ERS. Two separate cohorts (GSE31210↗, GSE72094↗) were used to validate the ERS. Specifically, we used the same model genes to calculate the ERS in these two validation cohorts. Patients in these two validation cohorts were also classified into high and low ERS groups, respectively, based on the median ERS values of the Meta cohort. Similarly, we also used Kaplan–Meier curves and time-dependent ROC curves to assess OS differences and predictive accuracy of ERS, respectively. Univariate and multivariate Cox regression analyses were applied to assess the independent prognostic value of ERS in the two validation cohorts.
Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were employed to explore the enrichment in biological function and pathway of the DEGs using the “clusterProfiler” R package [30, 31]. p values < 0.05 were statistically significant.
Development and verification of nomogram
We established the nomogram in the Meta cohort using the constructed prognostic model and the above-collected clinicopathological parameter to anticipate the prognosis of LUAD patients by the ‘rms’ R package. Calibration curves were used to evaluate the fitness between actual and predicted survival statuses with the established nomogram. Decision curve analysis (DCA) and ROC curves were applied to identify whether this nomogram model was suitable for clinical utility using the “rms”, “rmda” and “timeROC” R packages [32–34].
Drug sensitivity
The “pRRophetic” R package [35] was used to calculate the half-maximum inhibitory concentration (IC50) for eight common chemotherapy drugs in LUAD, including cisplatin, docetaxel, methotrexate, paclitaxel, erlotinib, gemcitabine, rapamycin and erlotinib, to predict the chemotherapy sensitivity of each LUAD sample based on the Genomics of Drug Sensitivity in Cancer database [36] (GDSC, https://www.cancerrxgene.org/↗).
Cell lines, RNA interfering
NCI-H1975 and A549 Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) culture medium, supplemented with 10% fetal bovine serum (FBS) in a standard humidified incubator with 5% CO2 at 37 °C. The knockdown of thymidine kinase 1 (TK1) in LUAD cells was achieved via the transfection of the TK1 specific small Interfering RNA (siRNA) using Lipofectamine 3000 reagent (Invitrogen, Massachusetts, USA) according to the manufacturer’s protocol. The TK1 specific siRNA were synthesized from GenePharma (Shanghai, China) and the sequence of siRNA is as the following: 5′-AGAAACUCAGCAGUGAAAGCC-3′. The knockdown efficiency was evaluated by Real-Time quantitative PCR (RT-qPCR) after 48 h transfection.
Migration, invasion, and proliferation assay
Colony formation assays were employed to determine cell viability. Cells were seeded at a density of 6 × 102 cells/plate in 6-well culture plates. After transduction with TK1-specific siRNA or random control, cells were cultured in normal medium for 10 days. Surviving tumor cells were fixed with 4% paraformaldehyde and stained with crystal violet before counting the colonies.
For wound healing assays, cells (1 × 105 cells/well) were seeded into 6-well plates. After the formation of an adherent confluent cell monolayer, cancer cells were starved for 8 h. Wounds were created using a 10 μL pipette, photographed at 0 h, and then again after 24 h. ImageJ was used to calculate the percentage of areas covered by migrated cells.
For in vitro migration and invasion assays, 24-well transwell chambers (Transwell, Corning Costar) were utilized. In migration assays, 4 × 104 cells suspended in serum-free medium were added to the upper chamber, while the lower chamber contained medium with 20% FBS. In invasion assays, the upper membranes were pre-coated with 40 μL matrigel (Matrigel™ GFR Membrane Matrix, #356231, Corning, USA), and 8 × 104 cells were seeded to the upper chamber. After 24 h of culture, cells in the upper chamber were carefully cleaned with a cotton swab. Cells attached to the filter's lower surface were fixed with 4% paraformaldehyde and stained with crystal violet. The cells on the lower surface of the membrane filter were captured on camera under a microscope.
Additionally, western blot (WB) was used to detect cell cycle-related proteins. Radio-immunoprecipitation assay (RIPA) buffer was used to prepare whole cell lysates, and proteins (15 μg/well) were separated with sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) gel (10%) before being transferred to polyvinylidene fluoride (PVDF) membranes. The membranes were incubated overnight with corresponding primary antibodies. Primary antibodies against CDK4 (1:1000, #12790), cyclinD (1:1000, #2978), cyclinE (1:1000, #4129), c-Myc (1:1000, #18583) and beta-Actin (1:1000, #4970) were purchased from Cell Signaling Technology (Danvers, MA, USA). Subsequently, HRP-conjugated secondary antibodies (ab6721, ab205719, Abcam) were incubated, and chemiluminescent signals were detected using enhanced chemiluminescence (ECL) solution.
Statistical analysis
Categorical data were compared using the chi-square test, and continuous variables were compared using the Mann–Whitney test or Kruskal–Wallis test, as appropriate. All p values of statistical data were based on two-sided statistical tests. p < 0.05 was considered to be statistically significant. All the statistical analysis was conducted by the R software (version: 4.1.0).
Results
Identification of ER-lncRNAs
We identified 30 DE-ER-mRNAs (Fig. 1A) and 437 DE-lncRNAs in tumor and normal patients in the TCGA-LUAD cohort. Subsequently, spearman correlation analysis was conducted on DE-ER-mRNAs and DE-lncRNAs in tumor patients, and 134 ER-lncRNAs were identified. Then, univariate COX analysis defined 19 as prognostic ER-lncRNAs (Table S1). Figure 1B shows the 19 lncRNAs expression level in tumor and normal patients.
Selection of ER-lncRNAs.Heatmap of differentially expressed ER-mRNAs between tumor tissues and normal tissues.Heatmap of differentially expressed ER-lncRNAs between tumor tissues and normal tissues A B
Identification of molecular subtypes based on prognostic ER-lncRNAs
Based on 19 prognostic ER-lncRNAs, patients in Meta cohort were divided into two clusters with good internal consistency and stability (Figs. 2A–C). We termed them cluster A and B. PCA analysis showed significant difference in the two clusters (Fig. 2D). Figure 2E shows the expression of 19 prognostic ER-lncRNAs between cluster A and cluster B. Survival analysis demostrated that the prognosis of patients in cluster A was significantly worse than those in cluster B (Fig. 2F).
ER-lncRNA subtypes.Consensus clustering matrix when= 2.Consensus clustering CDF withvalued 2–9.Relative change in area under the CDF curve for= 2 through 9.PCA analysis showed that LUAD could be well differentiated into two subtypes based on the expression of ER-lncRNA.Heatmap of the 19 prognostic ER-lncRNAs between cluster A and cluster B.Kaplan–Meier curves of OS for two subtypes in LUAD A B C D E F k k k
Immune characterization between different ER-lncRNA-related subtypes and drug sensitivity analysis
There are significant differences in several immune cells, functions and pathways between cluster A and cluster B (Fig. 3A). While the two clusters had no significant differences in stromal and immune scores (Fig. 3A). In addition, we investigated the expression levels of six immune checkpoint genes (i.e., PDCD1, CD274, CTLA4, PDCD1LAG2, LAG3 and ERBB2) and discovered that most of them were expressed at higher levels in cluster A than in cluster B, while ERBB2 was expressed at higher levels in cluster B (Fig. 3B). Subsequently, we explored the sensitivity of classical chemotherapy drugs in Meta cohort patients. As the IC50 value indicated, the patients in the cluster A were more sensitive to cisplatin, docetaxel, methotrexate, paclitaxel, erlotinib, and gemcitabine, while the patients in the cluster B were more sensitive to rapamycin and erlotinib (Fig. 3C).
Immune characterization analysis and drug sensitivity analysis of ER-lncRNA-related subtypes.The expression of immune function between different subtypes.The expression levels of immune checkpoint genes between different subtypes.The boxplot of sensitivity of common chemotherapy drugs in different subtypes. ns, not significant, *< 0.05, **< 0.01, ***< 0.001 A B C p p p
Biological characteristics of ER-lncRNA-related subtypes
To explore the reasons for the different prognosis and immune characteristics, we investigated the hallmark pathways in cluster A and B. The GSEA results revealed that cluster A were enriched in DNA replication and the cell cycle related pathways, such as E2F target, G2M checkpoint, MYC target, MTORC1, and PI3K-Akt pathway (Fig. 4A). Cluster B has more metabolism processes, such as heme metabolism, bile acid metabolism, and fatty acid metabolism (Fig. 4B).
Biological characterization of ER-lncRNA-related subtypes. The GSEA pathway enrichment analysis in cluster Aand cluster B () A B
Establishment and verification of an ER-lncRNA‐related prognostic model in LUAD
Identification of the DEGs and functional annotation
To further explore the potential biological activity of ER-lncRNA modification subtypes, we recognized 286 ER-lncRNA-subtypes-related DEGs (Table S2), 261 of which were significantly relevant to the OS of LUAD (Table S3). GO and KEGG results indicated that the 261 prognostic genes were primarily concentrated in cell cycle, metabolism and cancer-related pathways (Fig. 5).
Enrichment analysis of DEGs between the two ER-lncRNA-related subtypes. Gene ontology (GO)analysis of the 261 prognostic DEGs. BP, biological process; CC, cellular component; MF, molecular function.Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the 261 prognostic DEGs.< 0.05 A B p
Establishment of ER-lncRNA‐related prognostic model
To accurately predict the ER-lncRNA-related subtype in each LUAD patients, we established a prognostic model based on the identified 261 prognostic ER-lncRNA‐related genes. Then, 13 out of the 261 prognostic ER-lncRNA‐related genes were selected as the most potent prognostic markers using LASSO Cox regression analysis. The formula for calculating the risk score is as follows: ERS = PRC1 × 0.005 + ANLN × 0.044 + SCN7A × (− 0.017) + MAMDC2 × (− 0.008) + TK1 × (0.003) + VIPR1 × (− 0.019) + CHRDL1 × (− 0.018) + GAPDH × 0.043 + LOXL2 × 0.115 + IRX3 × (− 0.015) + CPA3 × (− 0.019) + CXCL17 × (− 0.037) + FAM83A × 0.038. LUAD patients were divided into high- and low- ERS groups based on a median ERS of 0.821.
Evaluation and validation of ER-lncRNA‐related prognostic model
Figure 6A shows the ERS and survival time for each patient in the Meta cohort. The Kaplan–Meier curve indicated that patients in the high-ERS group had a significantly worse prognosis than those in the low-ERS group (Fig. 6B). To evaluate the validity of ERS, we established a ROC curve. The AUCs for 1, 3 and 5 years was 0.695, 0.707 and 0.696, respectively (Fig. 6C).
Furthermore, according to the same prognostic model, 227 samples from GSE31210↗ cohort (Fig. 6D–F) and 400 samples from GSE72094↗ cohort (Fig. 6G–I) were validated the superior prognostic prediction abilities of the ERS. Kaplan–Meier curves of the GSE31210↗ and GSE72094↗ validation cohorts confirmed that patients in the high-ERS group had shorter OS than those in the low-ERS group. The AUCs of the ERS in GSE31210↗ cohort for 1, 3 and 5 years were 0.686, 0.681 and 0.684. And the AUCs of the ERS in GSE72094↗ cohort for 1, 3 and 5 years were 0.661, 0.661 and 0.696. These results suggested that the prognostic model based on the 13 ER-lncRNA-related genes may be a reliable predictor of prognosis in LUAD patients.
Evaluation and validation of ER-lncRNA‐related prognostic model. Distribution of ERS and survival time, Kaplan–Meier curves of overall survival and ROC curves in TCGA-LUAD (A-C),(–) and(–) cohorts GSE31210 GSE72094 B F G I
ERS is an independent prognostic predictor for LUAD patients
Univariate cox analysis revealed that ERS, T-stage, N-stage, M-stage, and pTNM staging were significantly related to OS in patients with LUAD in the Meta cohort (Fig. 7A). Further multivariate cox analysis demonstrated that the ERS was an independent prognostic predictor with HR = 3.355 (95% CI 2.331–4.828), p < 0.001 (Fig. 7B). Similarly, the ERS proved to be an independent prognostic factor in GSE72094↗ cohort (Fig. 7C, D).
Predictive nomogram. Univariate () and multivariateanalysis of the clinicopathologic features and the ERS in Meta cohort. Univariate () and multivariateanalysis of the clinicopathologic features and the ERS incohort.Nomogram to predict the survival of the LUAD patients.Calibration curve for 1-, 3-, and 5-year survival. ROC () and DCAcurve of the clinicopathologic features and risk score A B C D E F G H GSE72094
Development and validation of the nomogram prediction model
Morever, we integrated the collected clinicopathological characteristics and the established ER-lncRNA-related prognostic model to develop an accurate prognostic nomogram model to predict the 1-, 3-, and 5- years OS of LUAD patients in the Meta cohort (Fig. 7E). The calibration curves of 1-, 3- and 5- year OS rates in the Meta cohort were close to the ideal reference line, which indicated that the nomogram model had good prediction performance (Fig. 7F). In addition, ROC and DCA results showed that the ERS and nomogram model had high reliability (Fig. 7G, H). The AUC of the 5-year OS predictions for the ERS and nomogram model were 0.722 and 0.750, respectively (Fig. 7G). The ERS and nomogram model had significantly higher benefits than extreme curves according to the DCA, and nomogram model had higher benefits than ERS (Fig. 7H).
Biological features of low- and high- ERS group
To explore the reasons for the different prognoses and immune landscapes, we investigated the hallmark pathways between the low- and high- ERS groups. The results showed that low-ERS group had more metabolism-related pathways (Fig. 8A). The high-ERS group had more cell cycle and cancer-related pathways (Fig. 8A).
Biological characterization, immunological characterization and drug sensitivity analysis of high- and low- ERS groups.The GSEA pathway enrichment analysis in low- (left) and high- (right) ERS groups.The expression of immune function between different groups.The expression levels of immune checkpoint genes between different groups.The boxplot of sensitivity of common chemotherapy drugs in different groups. ns, not significant, *< 0.05, **< 0.01, ***< 0.001 A B C D p p p
Assessment of immune characteristics and drug sensitivity in different ERS patients
There are significant differences in several immune cells, functions and pathways between low- and high- ERS groups (Fig. 8B). And the high-ERS group showed higher stromal scores than low-ERS group (Fig. 8B). In addition, except ERBB2, the immune checkpoint genes expression level (i.e., PDCD1, CD274, CTLA4, PDCD1LAG2, LAG3) were significantly higher in high-ERS group than in low-ERS group (Fig. 8C). Subsequently, the IC50 value indicated, the patients in the high-ERS group were more sensitive to classical chemotherapeutic drugs than those in low-ERS group, except for rapamycin and erlotinib (Fig. 8D).
Knockdown of TK1 inhibited cell proliferation and metastasis of LUAD
We further investigated the role of TK1 in LUAD cells, including NCI-H1975 and A549. The effectiveness of TK1 silencing in NCI-H1975 and A549 was confirmed through RT-qPCR validation (Fig. 9A). Results from colony formation assays indicated that TK1 knockdown led to a substantial suppression of cell proliferation (Fig. 9B). Subsequent wound healing assays demonstrated a significant inhibition of cell migration in NCI-H1975 and A549 upon TK1 inhibition (Fig. 9C). Moreover, the transwell migration assay illustrated that TK1 knockdown markedly diminished the migration capability of both NCI-H1975 and A549 (Fig. 9D), while the invasion ability of these cells was also compromised after TK1 knockdown, as evidenced in Fig. 9E. Additionally, there was a noteworthy reduction in the expression of cell cycle-related proteins in the TK1 knockout cell lines (Fig. 9F). Therefore, the silencing of TK1 had a substantial impact on suppressing the proliferation and metastasis of LUAD cells.
TK1 promote cell proliferation, migration, invasion and cell-cycle progression in LUAD.The efficiency of TK1 silencing was confirmed by RT-qPCR in NCI-H1975 and A549 cell lines.The colony assay revealed that the knockdown of TK1 significantly inhibited the proliferation ability of both NCI-H1975 and A549 cells.Wound healing assays demonstrated that TK1 knockdown exerted a pronounced inhibitory effect on the migration of LUAD tumor cell lines (NCI-H1975 and A549).Transwell migration assays indicated that TK1 knockdown compromised the migration ability of both NCI-H1975 and A549 cells.Transwell invasion assays illustrated that TK1 knockdown attenuated the invasion ability of NCI-H1975 and A549 cells.Analysis of cell cycle-related proteins in TK1 knockout cell lines revealed a significant decrease in their expression. The data presented herein are the result of at least three independent experiments A B C D E F
Discussion
Currently, numerous studies have demonstrated that many exosomal lncRNAs expression levels are significantly different from normal tissues in pathological states, especially in tumor tissues. This suggests that exosomes can selectively package, translocate and secrete specific lncRNAs to regulate the corresponding biological processes [37, 38]. Exosomal lncRNAs have multiple biological functions in tumor biology, affecting tumor cell proliferation, apoptosis, migration, promoting angiogenesis, and also play an essential role in regulating the TME and mediating drug resistance [20, 39, 40]. Exosomal lncRNAs are generally considered as novel biomarkers for early diagnosis, prognosis prediction and efficacy assessment for tumors [41, 42]. Therefore, systematic studies on the prognostic implications and drug response prediction of exosomal lncRNAs in LUAD are necessary.
To our knowledge, this is the first systematic elucidation of the role of ER-lncRNAs in LUAD prognosis and drug response. We identified 134 DE-ER-lncRNAs using spearman correlation analysis, 19 of which were associated with prognosis in LUAD. Then, based on the 19 prognostic ER-lncRNAs expression, we identified two subtypes of ER-lncRNAs related in LUAD—subtype A and subtype B. Significant prognostic differences were observed between the two ER-lncRNAs-related-subtypes. GSEA results revealed that subtype A was enriched in DNA replication and cell cycle-related pathways, while subtype B had more metabolic processes. This is consistent with previous studies showing that LUAD subtypes enriched in metabolism-related pathways have a good prognosis [43]. In addition, there is a significant difference in immune cell infiltration between the two ER-lncRNAs-related-subtypes. Meanwhile, mRNA expression levels of immune checkpoints (i.e., PDCD1, CD274, CTLA4, PDCD1LAG2, LAG3, and ERBB2) were higher in subtype A than in subtype B. Consistent with this, subtype A was more sensitive to common therapeutic drugs, except rapamycin and erlotinib. Therefore, we hypothesized that ER-lncRNAs are relevant to prognosis and treatment response in LUAD patients.
Prognostic scoring systems based on gene expression have been established as robust predictors of patient outcomes [44, 45]. Although several prognostic scoring systems based on lactate metabolism-related lncRNA [46] and pyroptosis-related lncRNA [47] have been developed in LUAD, there is a scarcity of prognostic scoring systems based on the features of ER-lncRNA. To accurately predict the ER-lncRNA modification pattern in individual patients, we developed a prognostic model based on DEGs of ER-lncRNA-related-subtypes identified in the Meta cohort and well validated in two external cohorts (GSE31210↗ and GSE72094↗). The prognostic model was used to calculate the ERS of each LUAD patient. The low-ERS patients showed better OS. Moreover, ERS was shown to be an independent prognostic factor for survival of LUAD patients and had better predictive indication value compared to other common clinical biomarkers such as T/N/M-stage and pTNM staging. To explore the reasons for the prognostic differences, we performed a GSEA analysis, which showed that the low-ERS group had metabolism-related pathways, while the high-ERS group had cell cycle and cancer-related pathways.
Since TME plays a crucial role in antitumor response [48, 49], we also investigated the immune features between high- and low-ERS groups. Tumor-infiltrating immune cells, as an indicator to assess TME and treatment response [50, 51], were significantly different between the high- and low-ERS groups. In addition, we observed higher stromal scores in the high-ERS group. Studies demonstrated that higher stromal groups are involved in tumor progression by reshaping antitumor immunity and influencing the response of immunotherapy [52, 53]. Immunotherapy has made great progress in recent decades, and immune checkpoint inhibitors (ICIs) have been approved as first-line agents for patients with advanced LUAD [54], which have shown a significant survival advantage in LUAD [55, 56]. In our study, immune checkpoint genes expression levels were significantly higher in the high-ERS group than in the low-ERS group, which means that patients in the high-ERS group may benefit more from ICI compared to those in the low-ERS group. Chemotherapy remains the first line of treatment for metastatic LUAD [57]. Due to the heterogeneity of tumor, LUAD patients exhibit different sensitivities to chemotherapy. There, we evaluated the predictive value of the ERS for chemotherapy sensitivity in LUAD patients. The above results suggested that ERS can predict the effects of the chemotherapeutic and immunotherapeutic treatments and may be useful in determining the most suitable drug regimen for individual LUAD patients based on the ERS.
TK1 participates in cell proliferation and plays a potential role as a diagnostic tool and prognostic factor in evaluating cancer treatment and progression [58]. Additionally, serum levels of TK1 have been shown to be associated with cancer staging, with higher levels of TK1 indicating a more severe prognosis [58]. A clinical study indicates that elevated levels of TK1 in the serum after chemotherapy for lung cancer signify treatment failure and a poorer prognosis [59]. Our experiments also preliminarily demonstrated that silencing TK1 had a significant impact on suppressing the proliferation, metastasis, and reducing the expression of cell cycle-related proteins in LUAD cells. Further research is needed to explore whether TK1 is included in exosomes of LUAD cells.
Our study has certain limitations. Firstly, our data were obtained from an online database. While the feasibility of our risk model was demonstrated by two independent external cohorts, additional data samples may be required to further validate the results and reduce potential biases. Secondly, the online data we utilized are all based on tissue samples. Further research is needed to explore whether TK1 and other genes in the ERS model are present in exosomes or blood, and to further investigate the potential biological mechanisms underlying the relationship among ERS model genes, exosome and LUAD prognosis. Further clinical studies are needed to explore the relationship between the expression of ERS model genes in exosomes and the prognosis and treatment response of LUAD patients.
Conclusions
We identified two ER-lncRNA-related subtypes based on the prognostic ER-lncRNAs we identified, with significant differences in clinical and immunological characteristics. In addition, based on the ER-lncRNA-related genes, we developed an ER-lncRNA-related prognostic model—ERS, and performed antitumor drug sensitivity analysis. Our findings can be used to design individualized treatment strategies for patients with different LUAD subtypes.
Supplementary Information
Below is the link to the electronic supplementary material.
Acknowledgements
Not applicable.
Author contributions
Wangping Li and Haitao Zhang contributed study design and revising the manuscript. Daifang Chu and Liulin Chen contributed data interpretation and statistical analysis. Daifang Chu contributed drafting manuscript. Daifang Chu and Liulin Chen contributed equally to the article. All authors contributed to data analysis, drafting or revising the article, gave final approval of the version to be published, and agree to be accountable for all aspects of the work.
Funding
This work was supported by grants from National Natural Science Foundation of China (No. 82103537).
Availability of data and materials
All the data used to support the results of this study can be found in the public databases: TCGA: https://portal.gdc.cancer.gov/repository↗. GEO: https://www.ncbi.nlm.nih.gov/↗. GDSC: https://www.cancerrxgene.org/↗. MSigDB: http://www.gsea-msigdb.org/gsea/downloads.jsp↗. ExoBCD: https://exobcd.liumwei.org/↗.
Declarations
Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Ethics approval and consent to participate
We performed this research though publicly available datasets. There was no ethical approval was required.
Consent for publication
The work has not been published previously, and it is not under consideration for publication elsewhere.
Footnotes
Contributor Information
Wangping Li, Email: qxd25@163.com.
Haitao Zhang, Email: zht_best@163.com.
References
Associated Data
Supplementary Materials
Data Availability Statement
All the data used to support the results of this study can be found in the public databases: TCGA: https://portal.gdc.cancer.gov/repository↗. GEO: https://www.ncbi.nlm.nih.gov/↗. GDSC: https://www.cancerrxgene.org/↗. MSigDB: http://www.gsea-msigdb.org/gsea/downloads.jsp↗. ExoBCD: https://exobcd.liumwei.org/↗.