Article Text
Abstract
Background Lung adenocarcinoma (LUAD) is a highly heterogeneous disease, posing significant challenges to accurate prognosis prediction. Mitochondria play a central role in the energy metabolism of eukaryotic cells and can influence programmed cell death (PCD) mechanisms, which are critical in tumorigenesis and cancer progression. However, the prognostic significance of the interplay between mitochondrial function and PCD in LUAD requires further investigation.
Methods We analyzed data from 1231 LUAD patients across seven global cohorts to develop a mitochondrial-related PCD signature (MPCDS) using machine learning. Validation was done using six immunotherapy cohorts (LUAD, melanoma, clear cell renal cell carcinoma; n=935) and a pan-cancer cohort of 21 tumor types. An in-house LUAD tissue cohort (n=100) confirmed the prognostic significance of nucleoside diphosphate kinase 4 (NME4). In vivo and in vitro experiments explored NME4’s role in immune exclusion.
Results The MPCDS demonstrated strong predictive performance for prognosis in LUAD patients, surpassing 114 previously published LUAD signatures. Additionally, MPCDS effectively predicted outcomes in immunotherapy patients (including those with LUAD, melanoma, and clear cell renal cell carcinoma). Biologically, MPCDS was significantly associated with immune features, with the high MPCDS group exhibiting reduced immune activity and a tendency towards cold tumors. NME4, a key gene within the MPCDS (correlation=0.55, p<0.05), was associated with poorer prognosis in LUAD patients with high expression, particularly in CD8 desert phenotypes, as validated by our in-house cohort. Multiplex immunofluorescence confirmed the spatial colocalization and exclusion relationship between NME4 and immune cells such as CD3+ T cells and CD20+ B cells. Further experiments revealed that NME4 regulated the proliferation and invasion of LUAD cells both in vitro and in vivo. Importantly, inhibiting NME4 increased the abundance and activity of CD8+ T cells and enhanced the antitumor immunity of anti-programmed cell death protein-1 therapy in vivo.
Conclusion The MPCDS provides personalized risk assessment and immunotherapy interventions for individual LUAD patients. NME4, a key gene within the MPCDS, has been identified as a novel oncogene associated with immune exclusion and may serve as a new target for LUAD intervention and immunotherapy.
- Tumor Microenvironment
- Immunotherapy
- Biomarker
- Lung Cancer
Data availability statement
All datasets pertinent to this study are accessible through the TCGA database (http://cancergenome.nih.gov/), GEO database (https://www.ncbi.nlm.nih.gov/geo/), Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa-human/), UCSC Xena database (https://xenabrowser.net/datapages/), the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/), TIP (http://biocc.hrbmu.edu.cn/TIP/index.jsp), or the data availability sections of the relevant publications. All data relevant to this investigation, whether generated or analyzed, are comprehensively detailed in this manuscript and its supplementary materials. For further inquiries or data requests, interested parties are advised to reach out to the corresponding authors.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
Mitochondrial function and programmed cell death (PCD) play crucial roles in tumor progression, but their interplay and prognostic significance in lung adenocarcinoma (LUAD) remain understudied, particularly in predicting immunotherapy response.
WHAT THIS STUDY ADDS
This study developed a mitochondria-related PCD signature (MPCDS) demonstrating superior prognostic performance across multiple cohorts, and identified nucleoside diphosphate kinase 4 (NME4) as a key gene, revealing its mechanism in promoting tumor immune evasion by significantly inhibiting CD8+ T cell infiltration and activation, providing new insights into the regulation of the immune microenvironment in LUAD.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
This research provides a novel risk assessment tool for LUAD patients and may promote the development of personalized immunotherapy strategies, particularly by enhancing anti-programmed cell death protein-1 therapy efficacy through targeting NME4, while offering an important theoretical foundation and potential targets for future clinical trial designs and new therapeutic strategy development.
Introduction
Malignant lung tumors are the most common malignant tumors of the respiratory system and the leading cause of cancer-related deaths in both men and women.1 The WHO classifies malignant lung tumors based on histological type into non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC),2 with the former accounting for approximately 85% and the latter 15%. NSCLC primarily includes adenocarcinoma and squamous cell carcinoma. Currently, lung adenocarcinoma (LUAD) is the most common type of NSCLC, comprising about 40% of lung cancer cases.3 Although surgery is the most effective treatment for stage I to II and some stage IIIA LUAD, the invasive nature of these tumors often makes complete surgical resection challenging. The advent of immunotherapy has ushered in a new era for LUAD patients, with immune checkpoint inhibitors significantly improving clinical outcomes and offering a neoadjuvant immunotherapy option for early resectable disease. Unfortunately, research indicates that while immunotherapy significantly improves the prognosis for some patients, the 5-year survival rate remains only 4%–10%.4 Therefore, there is an urgent need to elucidate the underlying molecular mechanisms and develop reliable molecular classification models to effectively assess prognosis and guide personalized treatment strategies for LUAD patients.
Programmed cell death (PCD), also known as regulated cell death, is a crucial physiological process that plays a pivotal role in maintaining tissue homeostasis and eliminating damaged or unnecessary cells. PCD can occur through various mechanisms, including apoptosis, anoikis, autophagy, alkaliptosis, cuproptosis, entosis, entotic cell death, immunogenic cell death, ferroptosis, lysosome-dependent cell death, methuosis, necroptosis, netotic cell death, NETosis, oxeiptosis, pyroptosis, parthanatos, and paraptosis. PCD involves a regulated program of intrinsic clearance mechanisms to remove damaged or unnecessary cells, thereby maintaining tissue health.5 Among these PCD mechanisms, several are associated with mitochondrial dysfunction. Mitochondria play a key role in providing energy for cellular functions, regulating cell signaling pathways, and controlling PCD. Recent studies suggest that cancer is primarily a mitochondrial metabolic disease, with mitochondrial dysfunction characterized by changes in structure, function, and dynamics being implicated in tumorigenesis.6 Apoptosis is a well-recognized PCD mechanism that plays a significant role in maintaining tissue homeostasis and eliminating damaged or unnecessary cells, characterized by a series of biochemical and morphological changes.4 Inhibition or resistance to apoptosis often contributes to cancer development. Pyroptosis is a form of PCD that occurs following inflammasome activation and caspase-1 cleavage, characterized by cell swelling, membrane rupture, and the release of pro-inflammatory cytokines. Ferroptosis is a recently discovered PCD characterized by iron-dependent cell death and lipid peroxidation. Autophagy is a cellular mechanism that maintains cellular homeostasis by degrading damaged proteins and organelles. Autophagy can promote cell survival or induce cell death, depending on the specific context in which it occurs. Previous studies have shown that USP15 may induce autophagy to negatively regulate lung cancer progression. Necroptosis is a form of PCD characterized by necrosis-like cell death, triggered by the activation of RIPK1 and RIPK3. Cuproptosis is a PCD induced by copper overload, characterized by lipid peroxidation and mitochondrial dysfunction. Entotic cell death occurs in viable cells and their neighboring regions, distinct from traditional apoptotic pathways, as it does not require the activation of apoptotic execution pathways. Netotic cell death is another form of PCD resulting from the release of neutrophil extracellular traps, commonly observed during infections or injuries. Parthanatos is a tightly regulated form of cell death induced by the overactivation of the nuclear enzyme PARP-1. Lysosome-mediated cell death involves the action of hydrolytic enzymes that permeate the cell membrane and enter the cytoplasm. Currently, ferroptosis, cuproptosis, and pyroptosis have garnered widespread discussion in the context of LUAD.7–9 Additionally, alkaliptosis is an emerging form of PDC regulated by intracellular alkalinization processes. Oxeiptosis, leveraging the reactive oxygen species sensing ability of KEAP1, is a recently identified cellular pathway that likely interacts with other cell death pathways.
In the context of LUAD, tumor cells can evade the activation of PCD mechanisms through various strategies. They may activate specific signaling pathways or alter their morphology to avoid being eradicated by the host. As our understanding of PCD mechanisms deepens, numerous drugs targeting these pathways have been developed. For instance, a novel HDAC1/2 inhibitor can effectively inhibit colorectal cancer by regulating the cell cycle and inducing apoptosis.10 GSDME-induced pyroptosis is a unique form of PCD that has shown potential as an antitumor immunotherapy. Additionally, research has indicated that obstructing ferroptosis can lead to resistance against programmed cell death protein-1 (PD-1)/programmed death-ligand 1 (PD-L1) therapy. These findings underscore the significance of PCD research in enhancing our understanding of LUAD and developing new anti-LUAD therapies.
Disruption of mitochondrial morphology, such as changes in shape, size, or cristae organization, can impair normal mitochondrial function and trigger PCD. Structural abnormalities may affect the release of pro-apoptotic factors from mitochondria, leading to caspase activation and subsequent apoptosis. Mitochondrial function is also closely associated with PCD mechanisms. Dysfunctional mitochondria with impaired oxidative phosphorylation and ATP production can induce cellular stress and initiate PCD pathways. To better link mitochondrial function with the aforementioned PCD mechanisms, this study screened a series of relevant biomarkers. Currently, several molecular markers have been identified as clinically significant for the diagnosis and prognosis of LUAD.11 Markers such as PD-L1 expression levels, tumor mutational burden (TMB), and hematological biomarkers play a crucial role in predicting treatment efficacy and are closely related to the prognosis of LUAD. Besides genetic factors, various environmental and lifestyle factors contribute to the development of LUAD. Carcinogens in tobacco smoke can induce differentiation of bronchial epithelial cells, thereby promoting tumor formation.12
In the process of screening and validating myriad biomarkers, various confounding factors—including sample selection bias, intratumoral and intertumoral heterogeneity, analytical variability, and publication bias—can significantly impact the accuracy and generalizability of results. Consequently, the identification of survival-associated genes through transcriptomic profiling is paramount for prognostic stratification and the selection of targeted therapeutic modalities. Over recent decades, a substantial body of evidence has elucidated the pivotal role of mitochondrial dysfunction and PCD mechanisms in the pathogenesis and metastatic progression of malignant neoplasms. Neoplastic cells must circumvent diverse modalities of cellular demise and mitigate mitochondrial dysfunction to sustain proliferation and invasion. Notwithstanding these advances, a comprehensive elucidation of the intricate interplay between mitochondrial perturbations and PCD pathways in LUAD remains elusive. Moreover, there is a paucity of detailed functional studies delineating these processes specifically within the context of LUAD tumorigenesis and progression. To bridge these knowledge gaps, we introduced a mitochondrial-related PCD signature (MPCDS) to predict the efficacy of therapeutic interventions and prognosis in LUAD. Through our research, we have uncovered the high heterogeneity13 among LUAD patients and evaluated their clinical outlooks.
Method
Dataset source
The transcriptomic, copy number variation (CNV), mutation, H&E slides, and clinical data of LUAD were successfully retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) and used as the training set for model construction. Six transcriptome datasets were obtained from the GEO database for model validation, including GSE13213 dataset14 (n=119), GSE26939 dataset15 (n=115), GSE29016 dataset16 (n=39), GSE30219 dataset17 (n=86), GSE31210 dataset18 (n=227), and GSE42127 dataset19 (n=134). Three LUAD immunotherapy cohorts (OAK dataset,20 POPLAR dataset,21 NG dataset22) and three pan-cancer immunotherapy cohorts (GSE91061 dataset,23 phs000452 dataset,24 and Braun_2020 dataset25) were used to validate the model’s ability to predict the prognosis of immunotherapy patients. A pan-cancer cohort covering 21 tumors was used to evaluate the model’s broad applicability (data sourced from UCSC, https://xena.ucsc.edu/). Two single-cell datasets were sourced from GEO (GSE189357 dataset26) and the Genome Sequence Archive at the BIG Data Center (HRA001130 dataset27). All datasets used in this study are sourced from those listed in online supplemental table 1. Mitochondrial and PCD genes were derived from previous studies28; see online supplemental table 2 for details.
Supplemental material
To ensure data uniformity and comparability, gene expression data were converted to transcripts per million format. The “combat” function in the “sva” package29 was applied to mitigate potential batch effects. Additionally, all datasets obtained from TCGA and GEO underwent log transformation to establish a standardized data format from the onset of the analysis. Principal component analysis (PCA) was used to evaluate batch effects between datasets.
Identification of prognostic MPCDS
The limma package was used to analyze differentially expressed mitochondrial PCD-related genes (MPRGs) between normal and tumor samples, using selection criteria of false discovery rate (FDR) <0.05 and log2 fold change (FC)>1. Univariate Cox regression analysis was conducted to evaluate genes with prognostic value for LUAD patients. Subsequently, using 10-fold cross-validation, we applied 101 combinations of 10 machine learning algorithms, including stepwise Cox, Lasso, Ridge, partial least squares regression for Cox (plsRcox), CoxBoost, random survival forest (RSF), generalized boosted regression modeling (GBM), elastic net (Enet), supervised principal components (SuperPC), and survival support vector machine (survival-SVM). The goal was to identify the most valuable MPCDS characterized by the highest concordance index (C-index).
Our machine learning approach is similar to the combination of 10 machine learning algorithms used in Liu et al’s article.30 We used the same algorithms: stepwise Cox, Lasso, Ridge, plsRcox, CoxBoost, RSF, GBM, Enet, SuperPC, and survival-SVM. Both studies employed 101 algorithm combinations and used the C-index as the metric for evaluating model performance. The main difference lies in the cross-validation method: we used 10-fold cross-validation, while Liu et al used leave-one-out cross-validation.
The accuracy of MPCDS was evaluated using receiver operating characteristic (ROC) curves and PCA. Subsequently, we comprehensively retrieved 114 prognostic signatures from published papers, including lncRNA and mRNA, and assessed the performance of MPCDS using the C-index as the evaluation metric.
Investigating potential mechanisms and pathways
To investigate the biological functions and pathway processes associated with MPRGs, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the “clusterProfiler” package. Specifically, we used MPRGs as inputs, converted them to Entrez IDs, and conducted GO and KEGG enrichment analyses with an adjusted p value <0.05 as the criterion. Additionally, the potential mechanisms were investigated using Gene Set Variation Analysis (GSVA). The cancer immunity cycle and pathways predicted to respond to immunotherapy were examined according to previously established methods.31 32 A comprehensive collection of pathway gene sets was sourced from the Molecular Signatures Database.
Comprehensive analysis of genomic alterations and immune cell infiltration
Genomic alterations, including recurrently amplified and deleted regions, were identified using GISTIC 2.0 analysis. The TMB was calculated using the R package “maftools.” Previous studies have evaluated immune cell infiltration abundance in TCGA data, with IFNG scores sourced from the TIDE website (http://tide.dfci.harvard.edu/). The Cancer Immunome Atlas (https://tcia.at/home) assesses the immunophenoscore (IPS) of LUAD patients and identifies populations suitable for immunotherapy. Additionally, the ssGSEA algorithm evaluates the abundance of immune cell infiltration and the activity of immune-related pathways in tumor samples. The TIMER2.0 database provides a summary of immune cell infiltration abundance in TCGA, incorporating results from multiple algorithms (see online supplemental table 3 for details). The IOBR package33 includes a variety of immune infiltration algorithms that are used to infer the abundance of immune cells within tumor samples.
Preprocessing and integration of single-cell gene expression data
The initial gene expression matrix underwent preprocessing using the Seurat R package (V.4.2.0).34 Genes expressed in at least 10 cells within each sample were included in the analysis. Subsequently, low-quality cells were excluded based on specific criteria: cells with more than 5000 or fewer than 200 expressed genes, or cells with more than 10% of unique molecular identifiers originating from the mitochondrial genome. Sample integration was performed using the harmony R package. Next, a set of highly variable genes was selected for PCA, and the top 30 significant principal components were chosen for Uniform Manifold Approximation and Projection dimension reduction to visualize gene expression patterns. Differentially expressed genes within each cell subpopulation were identified using the “FindAllMarker” function and annotated based on previously published reliable cell type markers.27
Cell–cell communication
We used CellChat35 to integrate gene expression data and analyze differences in cell–cell communication modules. The default CellChatDB was used as the ligand-receptor database, following the standard CellChat protocol. By detecting overexpressed ligands or receptors within cell groups, we inferred cell type-specific interactions and identified enhanced ligand-receptor interactions when ligands or receptors were highly expressed.
Analysis of drug sensitivity in cancer cell lines and LUAD patients
Expression profile data of human cancer cell lines were sourced from the Broad Institute Cancer Cell Line Encyclopedia database (https://sites.broadinstitute.org/ccle). CERES scores for 18333 genes in 739 cell lines were obtained from the DepMap portal (https://depmap.org/portal/). CERES scores measure the dependency of genes in specific cell lines, with lower scores indicating a higher likelihood that the gene is essential for cell growth and survival. Drug sensitivity data for cancer cell lines were obtained from the Cancer Therapeutics Response Portal (CTRP,36 https://portals.broadinstitute.org/ctrp.v2.1/) and the PRISM database37 (https://www.theprismlab.org/). These datasets provide the area under the dose-response curve (AUC) as a measure of drug sensitivity, with lower AUC values indicating greater sensitivity. Missing AUC values were imputed using the K-nearest neighbor method, and compounds with more than 20% missing data, as well as hematopoietic and lymphoid tissue cell lines, were excluded.38 The sensitivity of each LUAD patient to chemotherapeutic agents was assessed based on the IC50 values quantified using the pRRophetic package from the Genomics of Drug Sensitivity in Cancer database.39 A ridge regression model was applied to generate drug sensitivity estimates for each patient.
Cell culture and small interfering RNA transfection
The A549 and H1299 LUAD cell lines were obtained from the Institute of Biochemistry and Cell Biology at the Chinese Academy of Sciences in Shanghai, China. The culture medium, RPMI 1640, was supplemented with 10% FBS and 1% antibiotics (100 U/mL penicillin and 100 mg/mL streptomycin) and was used to maintain the A549 and H1299 cells, respectively. Small interfering RNA (siRNA) transfection was performed using the Lipo2000 reagent (Invitrogen, Shanghai, China) strictly following the manufacturer’s protocol. Typically, A549 and H1299 cells were seeded onto coverslips in six-well plates, and siRNA transfection was carried out the following day. qRT-PCR was used to assess knockdown efficiency. The primer information for nucleoside diphosphate kinase 4 (NME4) and GAPDH is provided in online supplemental table 4.
Colony formation
5000 cells were seeded into each well of a six-well plate for the colony formation assay, with conventional growth medium added and replaced after 1 week. After 2 weeks, once the colonies had matured, they were fixed with methanol for 15 min and then stained with 0.1% crystal violet (Sigma) for 30 min. Following staining, the colonies were counted to assess their formation capability.
Invasion and migration assays
Invasion and migration assays were conducted using Corning’s Transwell system (24 wells, 8 µm pore size, New York, USA). For migration assays, 5×104 transfected cells were placed in the upper chambers with 350 µL of serum-free medium, while 700 µL of medium containing 10% FBS was added to the lower chambers. The invasion assays used Transwell membranes precoated with Matrigel (Sigma-Aldrich). After a 16-hour incubation, cells on the upper surface were removed, and those that migrated to the lower surface were stained with methanol and 0.1% crystal violet. Images were captured using an Olympus inverted microscope (Tokyo, Japan).
Clinical and pathological characteristics of LUAD patients
All paraffin-embedded tissue sections were obtained from the Department of Pathology at Tianjin Medical University Cancer Institute and Hospital. All case slides were histopathologically confirmed as LUAD, and none of the patients received radiotherapy, chemotherapy, or targeted therapy prior to surgery. The clinical and pathological characteristics of the patients are shown in online supplemental table 5.
Multiplex immunohistochemistry protocol for NME4 and immune cell marker detection
According to the standard protocol,40 multiplex immunohistochemistry (mIHC) was used to examine the expression relationship between NME4 and immune cell markers, with simultaneous detection of DAPI. The primary antibodies used for mIHC included: NME4 (1:200 dilution, Cat# ab228005, Abcam, USA), CD3 (1:150 dilution, Cat# ab21703, Abcam, USA), and CD20 (1:100 dilution, Cat# ab64088, Abcam, USA). The immunohistochemistry protocol involved deparaffinizing the LUAD tissue sections using xylene, followed by rehydration in a graded series of ethanol solutions. To inhibit peroxidase activity, the sections were treated with 3% hydrogen peroxide and then blocked with 5% normal goat serum. Subsequently, the sections were incubated overnight at 4°C with primary antibodies, including NME4 (1:200 dilution, Cat# ab228005, Abcam, USA) and CD8A (1:2000 dilution, Cat# ab217344, Abcam, USA). Following this, HRP-conjugated secondary antibodies were applied and incubated for 30 min at the same temperature. The sections were then stained with DAB (3,3′-diaminobenzidine) and counterstained with hematoxylin for visualization. Two independent pathologists independently evaluated all stained sections. For the semi-quantitative evaluation of NME4 staining, the H-score criterion was used. Additionally, tumors were classified into three phenotypes based on the spatial distribution of CD8+ T cells: inflamed, excluded, and desert subtypes.
Subcutaneous xenograft experiments in C57BL/6 mice
C57BL/6 mice were purchased from Jiangsu Gempharmatech (Jiangsu, China). For subcutaneous xenograft experiments, 2×106 LLC cells (sh-NME4) were resuspended in 100 µL PBS and injected subcutaneously into the mice. Once the tumors reached approximately 100 mm3, the mice were treated with either mouse PD-1 monoclonal antibody (mAb) (BioXcell, BE0146; New Haven, CT, USA) or IgG isotype control (BioXcell, BE0089) via intraperitoneal injection (100 mg per mouse in 100 µL D-PBS buffer). Tumor size was measured every 2–3 days. Tumor weight was recorded, and volumes were estimated using the formula 1/2×(length×width2).
Flow cytometry analysis of tumor-infiltrating immune cells in C57BL/6 mice
After euthanizing the C57BL/6 mice, fresh tumor tissue was excised and digested with collagenase IV (Gibco 17104019) and DNase I (Roche 10104159001) for 2 hours. The digested tissue was then filtered through 70 µm cell strainers to obtain a single-cell suspension. The cells were washed, resuspended in 40% Percoll, and centrifuged against 70% Percoll to isolate lymphocytes. For intracellular cytokine staining, the cells were stimulated for 4 hours with a leukocyte activation cocktail (BD Pharmingen, 550583; San Diego, CA, USA). Non-specific antibody binding was blocked with CD16/CD32 antibody (553141), followed by surface staining with anti-CD45 (553080), anti-CD3e (553064), and anti-CD8α (551162) antibodies (BD Pharmingen). After fixation and permeabilization with a Fixation/Permeabilization kit (BD Biosciences, 554714), intracellular GZMB was stained with an anti-GZMB antibody (Biolegend, 372204; San Diego, CA, USA). Stained cells were analyzed using a BD FACS Canto II flow cytometer, and the data were processed with FlowJo software.
Statistical analysis
Data manipulation, statistical analyses, and visualization were performed using R software V.4.2.0. Kaplan-Meier analysis and log-rank tests were used to estimate and compare subtype-specific overall survival (OS). Differences in continuous variables between two groups were assessed using the Wilcoxon test or t-test, while categorical variables were evaluated with χ2 tests or Fisher’s exact tests. P values were adjusted using the FDR method. Pearson correlation analysis was employed to explore relationships among variables. All p values were calculated using a two-tailed test, with p<0.05 considered statistically significant.
Result
Elucidating the role of mitochondrial and PCD-related biomarkers
Tumor cells can evade PCD mechanisms through various strategies, affecting the tumor microenvironment and thereby limiting immune response, which has become a significant obstacle in current cancer treatments. Numerous studies have demonstrated the close connection between mitochondrial function and PCD. To better elucidate this relationship, our study screened a series of mitochondrial and PCD-related biomarkers specifically for LUAD. Figure 1 presents the overall framework of this study. To ensure the consistency of data across different databases, we removed batch effects from normal samples in both the GTEx and TCGA databases (online supplemental figure S1A). We integrated normal lung tissue data from GTEx and LUAD samples from TCGA to perform differential expression analysis and identify differentially expressed MPRGs (figure 2A). Next, we conducted univariate Cox analysis to identify MPRGs associated with prognosis (figure 2B), ultimately identifying 92 prognosis-related MPRGs (p<0.05, online supplemental table 6). Further GO and KEGG enrichment analyses revealed that MPRGs are mainly enriched in mitochondrial components, regulation of autophagy, apoptosis, PI3K-Akt, other signaling pathways, and so on (figure 2C,D). The chromosomal locations and expressions of MPRGs are shown in figure 2E. To ensure comparability between datasets, we removed batch effects from seven LUAD transcriptome cohorts (figure 2F). CNV status analysis indicated frequent changes in MPRGs. Notably, STYXL1 and STAT1 exhibited the most widespread CNV amplification (figure 2G).
Supplemental material
Development and validation of an MPCDS for prognosis and immunotherapy suitability in patients with cancer
Using the expression profiles of 92 prognosis-related MPRGs, we developed an MPCDS through a machine learning combinatorial algorithm. The TCGA dataset served as the training cohort, while six GEO datasets were used for validation. The average C-index across the six validation cohorts was the criterion for model selection. Ultimately, the Lasso+Stepcox (forward) algorithm emerged as the optimal model (figure 3A). The MPCDS score distinguished patient prognosis across all seven cohorts (figure 3B–H). Patients in the high MPCDS group exhibited worse outcomes compared with those in the low MPCDS group. Interestingly, we extrapolated the MPCDS scores for the immunotherapy cohorts using the model’s formula and found consistent results in three LUAD immunotherapy cohorts (figure 3I–L) and three pan-cancer cohorts (including two melanoma cohorts and one clear cell renal cell carcinoma cohort, figure 3M–P). The low MPCDS group appeared to be more suitable for immunotherapy, as they demonstrated better prognoses following treatment (figure 3I–P). Furthermore, IPSs calculated using the TCIA website also suggested that the low MPCDS group is more appropriate for treatment with PD-1, CTLA-4, or their combination therapy (figure 3Q–T).
Superior predictive performance of MPCDS in prognosis
To evaluate the predictive efficacy of MPCDS, we integrated clinical features from seven datasets. MPCDS demonstrated higher C-index values than any other clinical feature (such as age, gender, stage, EGFR status, etc) (figure 4A). Subsequently, we conducted a PCA based on the expression levels of the model genes of MPCDS across all datasets. The results showed that the expression levels of the model genes could effectively distinguish LUAD patients, with high and low MPCDS groups forming distinct clusters (figure 4B). ROC curves further confirmed that MPCDS scores could reliably predict the prognosis of LUAD patients, which was validated across all datasets (with 1-, 3-, and 5-year AUC values generally above 0.65, figure 4C). Next, we compared MPCDS against 114 previously published LUAD signatures and found that MPCDS consistently demonstrated the best predictive performance across all datasets, achieving the highest C-index values (figure 4D). Further exploration of the correlation between MPCDS and clinical indicators revealed that patients with high MPCDS scores had a higher number of deaths, and tended to have more advanced T stage, N stage, and overall stage (online supplemental figure S1B). Additionally, the model genes LDHD, CTSG, and ADRB2 were highly expressed in the low MPCDS group, while the remaining model genes were highly expressed in the high MPCDS group. These observations indicate that MPCDS is an effective prognostic predictor and is significantly associated with clinical characteristics.
Complex relationship between MPCDS and key cancer immunity and biological pathways
To explore the potential pathway mechanisms underlying the differences between different MPCDS groups, GSVA enrichment analysis was conducted. The primary pathways enriched in the high MPCDS group include DNA replication checkpoint signaling, regulation of mitotic cytokinesis, positive regulation of DNA-templated DNA replication, DNA unwinding involved in DNA replication, phosphorylation of EMI1, and MYC gene upregulation-related target genes (online supplemental figure S2A–G). The abnormal activation or suppression of these pathways leads to uncontrolled cell cycle progression, increased DNA replication stress, accumulation of genetic variations, and significant chromosomal structural changes. These collectively promote the unlimited proliferation of tumor cells, genomic instability, and enhanced invasive and metastatic potential. Next, we explored the associations between MPCDS scores and immune therapy-related pathways as well as cancer immunity cycle steps. First, MPCDS showed significant positive correlations with various cancer immunity cycle steps (online supplemental figure S2H), such as cancer antigen release (step 1) and the recruitment of immunosuppressive cells (eg, MDSCs and neutrophils). This suggests that high MPCDS score groups may exhibit heightened activity in these immune steps, leading to more effective antigen release and recruitment of immunosuppressive cells. Second, MPCDS demonstrated strong correlations with several key biological pathways, such as DNA replication, cell cycle regulation, and mismatch repair. These pathways may be abnormally activated or suppressed in high MPCDS score groups, thereby affecting tumor progression and immune response. For example, abnormally active DNA replication and cell cycle regulation can result in rapid tumor cell proliferation and genomic instability, while alterations in mismatch repair mechanisms may increase mutation rates. Surprisingly, when we extended the analysis of MPCDS to a pan-cancer level, GSVA confirmed our previous results. MPCDS scores were significantly positively correlated with cell cycle-related pathways such as E2F targets, MYC targets V1/V2, and the G2M checkpoint across 21 types of tumors (online supplemental figure S5A,B). Overall, patients with high MPCDS scores exhibit more pronounced activity in these key pathways and immunity cycle steps, which are closely related to tumor invasiveness and prognosis. This indicates a complex relationship between MPCDS and multiple cancer-related biological pathways and immunity cycle steps, where high MPCDS scores may predict poorer prognosis and higher tumor aggressiveness.
Exploring the relationship between tumor mutations and patient prognosis
Given the significant correlation between tumor mutations, immune response, and patient prognosis, this analysis explores these relationships in detail. Online supplemental figure S3A,B vividly shows distinct chromosomal alterations between high- and low-MPCDS groups. The heatmap highlights a markedly elevated TMB in the high-MPCDS group, with TP53, TNN, CSMD3, RYR2, and LRP1B being the most frequently mutated genes (online supplemental figure S3C). Analyses reveal more frequent chromosomal deletions or amplifications within the high-MPCDS group (online supplemental figure S3D). Furthermore, online supplemental figure S3E indicates the worst prognosis in the L-TMB+ high-MPCDS subgroup.
Immune infiltration and regulatory gene expression in high and low MPCDS groups
Genomic-level comparisons revealed that the proportions of leukocytes, lymphocytes, and NK cells were higher in the low MPCDS group than in the high MPCDS group (p<0.05) (figure 5A–D). Additionally, we used the tumor-infiltrating lymphocyte (TIL) fraction data provided by Saltz et al41, who employed deep learning methods to estimate TILs on H&E-stained slides. Strikingly, the H&E-estimated TIL fraction was higher in the low MPCDS group (p<0.05) (figure 5E). The IFNG levels estimated by the TIDE website were also higher in the low MPCDS group (figure 5F). Immune regulatory genes play a significant role in tumor immune response, and we further compared the relative expression levels of immune regulatory genes between the two groups. The low MPCDS group appeared to favor the activation of MHC class II molecules and co-stimulatory molecules, whereas the high MPCDS group showed activation of MHC class I molecules and certain co-inhibitory molecules(figure 5G). Furthermore, we used the ssGSEA algorithm to calculate the abundance of immune infiltrating cells and immune-related pathway activity between the high and low MPCDS groups. The results were consistent with previous findings: the low MPCDS group had higher immune cell abundance, such as B cells and helper T cells, and the activation of type II interferon response and antigen presentation-related immune pathways (figure 5H). The resultsN calculated using the TIMER2.0 database further corroborated these findings (figure 5I). To confirm the above results, we examined H&E slides from the TCGA database for the high and low MPCDS groups and found more lymphocyte infiltration in the low MPCDS group (figure 5J).
Detailed exploration of MPCDS distribution and cellular interactions in different cell types
We explored the detailed distribution of MPCDS across different cell types using single-cell RNA transcriptome data (GSE189357 and HRA001130). We annotated the major cell types in these datasets (figure 6A,D) and examined the distribution and intensity of MPCDS scores across different cell populations (figure 6B,E). Based on the median MPCDS score, we divided the cell populations into high and low MPCDS groups to observe cell type proportions (figure 6C,F). We found that the proportion of cycling cells was higher in the high MPCDS group, while immune cells such as B cells, plasma cells, and NK/T cells were lower. This finding is consistent with the bulk transcriptome analysis results. Violin plots also demonstrated that cycling cell types had the highest MPCDS scores (figure 6G,H), indicating that higher MPCDS scores may be associated with greater proliferative and stem cell potential. Furthermore, we examined cell interactions between high and low MPCDS groups (figure 6I–K) and found stronger cell–cell communication and more robust outgoing and incoming signals in the high MPCDS group (figure 6L,M). The high MPCDS group is likely to have stronger pro-tumor signaling pathways, including VEGF, EGF, and PDGF pathways (online supplemental figure S4).
Identification of therapeutic targets and candidate drugs for high MPCDS score patients
Proteins that show a high positive correlation with MPCDS may have potential therapeutic implications for patients with high MPCDS scores. However, most human proteins remain undruggable because they lack clear active sites for small molecule compounds to bind, or they are located within cells, making them inaccessible to biological agents. Therefore, to identify potentially druggable therapeutic targets for patients with poor prognosis in the high MPCDS group, we calculated the correlation coefficients between the expression levels of druggable proteins and MPCDS, identifying protein targets with correlation coefficients greater than 0.30 (p<0.05). Subsequently, using the CERES scores and MPCDS scores from LUAD cell lines, we identified genes with Spearman’s r<−0.3 and p<0.05. Through these two analyses, we identified seven genes: PARP1, HK2, P4HTM, PLOD2, LDHA, GPX8, and RFK (online supplemental figure S6A–O). Notably, the CERES scores of PLOD2 were greater than zero in most LUAD cell lines, suggesting that PLOD2 may not be essential in LUAD. The remaining six genes are considered potential therapeutic targets, indicating that inhibiting these genes in patients with high MPCDS scores may achieve favorable treatment outcomes.
To identify candidate drugs with higher sensitivity in patients with high MPCDS scores, we employed two different methods, using drug response data derived from CTRP and PRISM. First, we conducted differential drug response analysis between the high MPCDS group (top decile) and the low MPCDS group (bottom decile) to identify compounds with lower estimated AUC values in the high MPCDS group (log2FC>0.10). Next, we performed Spearman correlation analysis between AUC values and MPCDS scores to select compounds with negative correlation coefficients (Spearman r<−0.50 for CTRP or Spearman r<−0.3 for PRISM, online supplemental figure S6P–S). These analyses identified six CTRP-derived compounds (Ispinesib, Paclitaxel, Volasertib, Docetaxel, NVP-AUY922, Litronesib) and three PRISM-derived compounds (SB-743921, BI-2536, Paclitaxel). All these compounds showed lower estimated AUC values in the high MPCDS group and negative correlations with MPCDS, indicating their potential as therapeutic drugs for patients with high MPCDS scores.
NME4 as a potential oncogenic biomarker
The previous analyses validated that MPCDS stratification effectively distinguishes differences in patient prognosis and the tumor immune microenvironment. To identify potential biomarkers associated with MPCDS, we focused on a key model gene, NME4. Currently, there is a lack of sufficient literature elucidating its role in tumor progression. Prior analyses have demonstrated that NME4 is a risk gene (HR>1, p<0.05). Further analysis revealed a significant positive correlation between NME4 and MPCDS (correlation=0.55, p<0.05, figure 7A). Consistent results from the kmplotter website showed that LUAD patients with high NME4 expression had poorer prognosis (figure 7B). The HPA database confirmed that NME4 is highly expressed in LUAD tumor samples compared with normal lung tissues (figure 7C,D). We further investigated whether NME4, similar to MPCDS, could differentiate the immune microenvironment. Heatmap results indicated that NME4 was negatively correlated with the majority of chemokines and receptors, MHC molecules, and immune co-stimulatory and co-inhibitory molecules (figure 7E), suggesting a potential role for NME4 in regulating tumor immunity. To explore the functional role of NME4 in LUAD tumorigenesis, we knocked down NME4 expression in A549 and H1299 cells using specific siRNA (online supplemental figure S7A). CCK8 assays demonstrated that the OD values of LUAD cells decreased with NME4 knockdown, indicating reduced proliferation capacity (figure 7F,G). Colony formation assays showed a significant reduction in the number of colonies in NME4 knockdown LUAD cells (figure 7H). Transwell assays confirmed that NME4 knockdown inhibited tumor cell invasion and migration (figure 7I,J). These findings indicate that NME4 functions as an oncogene, regulating the proliferative capacity of LUAD cells.
NME4 expression and immune cell infiltration
NME4, as a marker of MPCDS, was further explored for its immune-related properties. The IOBR package predicted immune cell abundance in LUAD tumor samples, and figure 8A shows that as NME4 expression increased, immune cells gradually decreased (eg, CD8 T cells, B cells, etc). This observation was repeatedly validated by multiple algorithms (figure 8A). We used in-house paraffin sections for staining NME4 and CD8 and found that patients with high NME4 expression had reduced CD8 cell infiltration (figure 8B). Based on the spatial distribution of CD8+ T cells, we classified the sections into inflamed, excluded, and desert phenotypes, finding that patients with the inflamed phenotype had the best prognosis (figure 8C), which also corroborates previous findings.42 The inflamed phenotype showed the lowest NME4 expression, the desert type the highest, and the excluded type in between (figure 8D). Further, by categorizing patients into high and low NME4 expression groups, we observed that patients with high NME4 expression had a poorer prognosis(figure 8E). To investigate the distribution of NME4 with other immune cells, we performed multiplex immunofluorescence staining, revealing that patients with low NME4 expression had more infiltration of CD3+ T cells and CD20+ B cells(figure 8F,G), further proving the relationship between NME4 and immune cell exclusion.
NME4 knockdown enhances CD8+ T cell infiltration and boosts anti-PD-1 therapy efficacy
To evaluate the immune regulatory effect of NME4 in vivo, we used the LLC mouse LUAD cell line to construct subcutaneous xenografts in immunocompetent C57BL/6 mice, which were then treated with PD-1 mAb or IgG isotype control (IgG2a) (figure 9A, online supplemental figure S7B). The results showed that, compared with the control group, NME4 knockdown significantly inhibited tumor growth, consistent with the results of the PD-1 mAb treatment. Additionally, the combination of NME4 inhibition and PD-1 blockade resulted in the most significant reduction in tumor burden, with no observed weight loss or other common toxic effects (figure 9B–E). At the end of the treatment, tumor samples were collected for further analysis. Flow cytometry analysis indicated that both PD-1 mAb and NME4 knockdown enhanced CD8+ T cell infiltration and activation. Notably, the combination of PD-1 mAb and NME4 knockdown led to the most significant enrichment of activated CD8+ T cells in the tumor regions (figure 9F–H). These findings demonstrate that NME4 knockdown enhances CD8+ T cell infiltration and function, thereby significantly improving the efficacy of anti-PD-1 therapy. Overall, these results indicate that NME4 promotes tumor immune escape by inhibiting CD8+ T cell infiltration, and that NME4 knockdown can reverse this process, thereby enhancing the effectiveness of anti-PD-1 therapy.
Discussion
LUAD is a complex and highly heterogeneous disease, making prognostic prediction for LUAD patients challenging.43 Prognostic factors play a crucial role in predicting disease progression, selecting appropriate treatment strategies, and estimating OS.44 Among the significant prognostic factors, histological subtype and grading are paramount. Higher grading is associated with more aggressive tumor behavior, shorter progression-free survival, and lower OS. Molecular markers have become powerful predictors of prognosis in LUAD. Age at diagnosis is another key factor influencing prognosis, with younger patients generally exhibiting longer OS compared with older patients. The extent and location of the tumor are also prognostically significant, with smaller tumor size and more accessible locations associated with better outcomes. Residual disease after surgical resection is another critical factor affecting prognosis. As our understanding of LUAD biology expands, emerging prognostic factors and novel molecular markers are being investigated to refine prognostic predictions and personalized treatment approaches.
Machine learning techniques are becoming increasingly prevalent in predicting patient survival. However, effectively implementing these methods in clinical practice while maintaining accuracy remains a challenge. Two important considerations are why a specific machine learning algorithm should be used and which solution is the best. Researchers’ choice of algorithms may largely depend on their preferences and biases.45 In this study, we collected expression profiles from 1231 LUAD patients across seven global multicenter cohorts and developed an MPCDS using a novel computational framework. The MPCDS is based on the expression of 10 genes that are most effective in predicting patient survival. Enrichment analysis revealed that high MPCDS group are primarily enriched in various cellular processes, environmental information processing, and organismal systems. Notably, our analysis uncovered a significant correlation between the high MPCDS group and biological processes related to DNA replication, the cell cycle, and other proliferation-related pathways. These findings partly explain the poorer prognosis observed in this group. Our results suggest that MPCDS can serve as a valuable tool for guiding treatment decisions and improving patient outcomes. Importantly, MPCDS can effectively differentiate the prognosis of patients in immunotherapy cohorts, providing clinicians with more reliable patient stratification and helping to identify those most suitable for immunotherapy.
NME4, also known as non-metastatic clone 23 human subtype 4 (nM23-H4), is a member of the nM23 family. NME4 is a mitochondria-associated gene whose protein contains a mitochondrial targeting sequence. Increasing evidence suggests that NME4 is involved in mitochondrial autophagy, inflammatory responses, invasion potential, and apoptosis. Research on NME4 has primarily focused on tumors such as oral cancer, breast cancer, gastric cancer, renal tumors, colorectal cancer, and testicular germ cell tumors. However, its role in LUAD remains underexplored. Studies indicate that reducing the n-6/n-3 fatty acid ratio can downregulate the expression of cell adhesion/invasion-related molecules (such as NME4), thereby inhibiting the invasion of human lung cancer cells. Wang et al46 suggested that NME4 may promote NSCLC progression by overcoming cell cycle arrest and promoting cell proliferation. In this study, we identified NME4 as a potential oncogenic biomarker for LUAD and highlighted its significant role in tumor progression and immune regulation. Our analysis validated that the MPCDS effectively distinguishes patient prognosis and the tumor immune microenvironment. Specifically, we found that NME4 is a key gene within this signature, showing a significant positive correlation with MPCDS and poor patient prognosis. High NME4 expression was associated with increased proliferation, colony formation, and invasion of LUAD cells, indicating its oncogenic role. Additionally, we explored the immune-related properties of NME4 and found that high NME4 expression was associated with reduced immune cell infiltration, particularly CD8+ T cells. Our in vivo experiments further demonstrated that NME4 knockdown enhances CD8+ T cell infiltration and activation, significantly improving the efficacy of anti-PD-1 therapy. These findings suggest that NME4 promotes tumor immune escape by inhibiting CD8+ T cell infiltration, and targeting NME4 can enhance the effectiveness of immunotherapy. In summary, NME4 is an important regulatory factor in the progression and immune evasion of LUAD, providing a new target for therapeutic intervention.
A notable limitation of this study is the lack of multicenter cohort validation to assess the accuracy of our model. While we have conducted extensive bioinformatics analyses and preliminary experimental validations, further in-depth investigation of the mechanisms involving NME4 is needed. Future research should aim to include a more diverse patient population and validate our findings across multiple independent cohorts to ensure robustness and applicability. The potential of targeting NME4 for the treatment of LUAD requires comprehensive functional studies and clinical trials.
In summary, the MPCDS can effectively predict LUAD patient prognosis and response to immunotherapy. As a key gene within the MPCDS, NME4 is a crucial oncogene in LUAD, significantly influencing tumor progression and immune regulation. High NME4 expression is associated with poor prognosis and immune cell exclusion, particularly CD8+ T cells. Inhibiting NME4 enhances CD8+ T cell infiltration and improves the efficacy of anti-PD-1 therapy, suggesting its potential as a therapeutic target.
Data availability statement
All datasets pertinent to this study are accessible through the TCGA database (http://cancergenome.nih.gov/), GEO database (https://www.ncbi.nlm.nih.gov/geo/), Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa-human/), UCSC Xena database (https://xenabrowser.net/datapages/), the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/), TIP (http://biocc.hrbmu.edu.cn/TIP/index.jsp), or the data availability sections of the relevant publications. All data relevant to this investigation, whether generated or analyzed, are comprehensively detailed in this manuscript and its supplementary materials. For further inquiries or data requests, interested parties are advised to reach out to the corresponding authors.
Ethics statements
Patient consent for publication
Ethics approval
All human experiments conducted in this study were approved by the Medical Ethics Committee of the Tianjin Medical University Cancer Institute and Hospital (approval number: bc2022232), and informed consent was obtained from all patients. For animal experiments, all procedures were reviewed and approved by the Animal Ethical and Welfare Committee of Tianjin Medical University Cancer Institute and Hospital (as per the official stamp of the committee), with the approval number 2023078. Both human and animal studies were conducted in strict adherence to the respective guidelines and regulations established by these committees, ensuring ethical standards were maintained throughout the research process.
Acknowledgments
This work was supported by the Tianjin Natural Science Foundation under Grant/Award Number 21JCYBJC01020, and funded by Tianjin Key Medical Discipline (Specialty) Construction Project (TJYXZDXK-011A).
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
LZ, YC, GZ and PZ contributed equally.
Contributors The study was conceived and designed by PZ. Data collection was conducted by YC. GZ performed the statistical analysis. The first draft of the manuscript was written by PZ. The experiment was performed by YC and PZ. The final approval of the submitted version was given by LZ and ZZ. All authors contributed to the manuscript and approved the submitted version. ZZ is responsible for the overall content as the guarantor.
Competing interests It is hereby declared by the authors that the research was carried out without the presence of any potential conflict of interest arising from commercial or financial relationships.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.