Identification of mitochondrial-related genes to evaluate the immune infiltration and prognosis of lung adenocarcinoma

Data collection

RNA-seq data and the clinical information of 550 LUAD samples (V36. 0, December 12, 2022) were downloaded from The Cancer Genome Atlas (https://portal.gdc.cancer.gov/), 14 samples were excluded due to incomplete clinical information or survival less than 30 days. In total, 536 samples, comprising 478 tumor samples and 58 healthy samples, were analyzed in the present study. Samples from two datasets, GSE31210 and GSE30219, were applied to validate the hub genes and risk score (https://www.ncbi.nlm.nih.gov/geo/). GSE31210 is annotated by GPL570 and includes 226 tumor samples and 20 healthy samples. GSE30219 is annotated by GPL570, comprising 293 tumor samples and 14 healthy samples. The list of mitochondrial related genes was collected from MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta). The research flowchart was shown in Fig. 1.

IHC images were retrieved from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/), a publicly available, peer-reviewed resource for comprehensive human protein expression data (Uhlén et al., 2015).

Differential expression genes analysis

To identify differentially expressed genes (DEGs), we utilized the “DESeq2” (version 1.47.3) package of R (version 4.3.3). Specifically, DEGs were screened under two comparison settings: (1) between normal and tumor samples, and (2) between high-risk and low-risk subgroups derived from the training set. DEGs were defined using the thresholds of |log₂ fold change| > 1.5 and adjusted P-value < 0.01. For visualization of the identified DEGs, a volcano plot was generated with the assistance of R’s “ggpubr” (version 0.6.0) and “ggthemes” (version 5.1.0) packages. Additionally, to visualize the overlapping DEGs in both DEGs groups and mitochondrial-related genes, a Venn diagram was constructed using R’s “Ggplot2” (version 3.4.4) and “VennDiagram” (version 1.7.3) packages.

Development and validation of a prognostic mitochondrial-related set

The selection of Mitochondrial-Related (MTR) genes was accomplished through bioinformatics approaches in R, including univariate Cox regression, LASSO regression (Kang et al., 2021; Yi et al., 2022), and multivariate Cox regression. Specifically, the “survival” package (version 3.5.7) was employed to perform univariate and multivariate Cox proportional hazards regression analyses, the “glmnet” package (version 4.1.8) was used to implement Lasso regression for further variable screening to refine the MTR gene set, and the “forestplot” package (version 1.12.1) was utilized to visualize the results of the Cox regression for intuitive presentation of the associations between the selected MTR genes and survival outcomes. The risk score for each individual sample was ascertained. In the selected cohorts, participant samples were stratified into low-risk and high-risk categories utilizing the median risk score as the threshold. The determination of the optimal cutoff point and the generation of Kaplan-Meier (K-M) survival curves were facilitated by the R packages “survival” (version 3.5.7) and “survminer” (version 0.4.97). The predictive accuracy was quantified through the receiver operating characteristic (ROC) curve, risk stratification plot, and the concordance index (C-index). Comprehensive data regarding the prognostic genes were sourced from the Human Gene Database (https://www.genecards.org/), the Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/).

Development and validation of a predictive nomogram

A comprehensive analysis was conducted to determine the risk score of clinical parameters such as age, gender, T stage, N stage, M stage, tumor stage. These variables were subjected to univariate Cox regression and multivariate Cox regression to refine the candidate predictors that exhibited a significant relationship with survival (P < 0.05). Utilizing these predictors, nomograms were developed, and corresponding scores were allocated for each variable within the nomogram framework. The aggregate score for each patient was derived by summing the individual scores of the included predictors. The performance of the nomogram model was evaluated using the calibration plots to assess its discriminative power and predictive accuracy. “Survival” and “rms” (version 6.7.0) packages in R were employed.

Analyses of Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, Gene Set Enrichment Analyses

For gene set functional enrichment analysis, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Gene Ontology (GO) annotations of genes from the R package “org.Hs.eg.db” (version 3.18.0) were used as the background. We mapped the genes to the background set and performed enrichment analysis using the R package “clusterProfiler” (version 4.10.1) to obtain the gene set enrichment results. For Gene Set Enrichment Analysis (GSEA), the GSEA 4.2.1 software was obtained from the GSEA website (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005). Samples were divided into two groups according to risk levels. Following subsets were downloaded from Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/index.jsp) (Liberzon et al., 2011), including c5.go.cc.v7.4.symbols.gmt, c2.cp.kegg.v7.4.symbols.gmt, c7.immunesigdb.v7.4.symbols.gmt, c7.all.v7.4.symbols.gmt, and c5.go.mf.v7.4.symbols.gmt. These subsets were used to assess related pathways and molecular mechanisms based on gene expression profiles and phenotype groupings. The minimum gene set size was set to 5 and the maximum to 5,000. One thousand re-samplings were performed. P < 0.05 and false discovery rates (FDR) < 0.1 were considered statistically significant.

Analyses of immune landscape in LUAD

TME was characterized by computing stromal scores, immune scores, ESTIMATE scores and tumor purity through R package “estimate” (version 1.0.13). Lists of TME-related biomarkers was extracted from Gene set enrichment analyses (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp). Tumor immune dysfunction and exclusion (TIDE) acted as a vital biomarker for immunotherapy response. Additionally, we assessed whether there was any link between TIDE scores and risk scores.

Immune infiltration analysis

To better characterize the immune cell features in high-risk and low-risk populations, we employed gene expression profiles from microarray processed as previously described in Zhang et al. (2023). Based on these profiles, we applied the single-sample gene-set enrichment analysis (ssGSEA) algorithm to estimate the differential infiltration abundance of 28 distinct immune cell types between the two groups. Spearman correlation analysis was used to elucidate the relationship between risk stratification and immune cell distribution patterns.

The dataset of immune cell signature was retrieved from the Cibersort data (https://cibersortx.stanford.edu/). The relative frequencies of 22 distinct tumor-infiltrating immune cell (TIIC) subtypes were ascertained via the CIBERSORT algorithm. A Spearman test was applied to assess the differences in TIICs’ relative frequencies between high-risk and low-risk groups.

Tumor mutation burden and immune checkpoint gene differences analysis

The somatic mutation data were downloaded from TCGA GDC data portal 2.0 database (https://portal.gdc.cancer.gov/). The R package “maftools” (version 2.18.0) was then used to draw a waterfall plot to illustrate the mutation landscape in LUAD patients with the high- and low-risk group and calculate the TMB score for each sample. In addition, identifying the relationship between high and low risk scores and the expression levels of common immune checkpoints.

Cell culture and transfection

Human LUAD cell lines A549, PC9, and H1299, together with the normal human pulmonary epithelial cell line BEAS-2B, were sourced from the China Center for Type Culture Collection (CCTCC, Shanghai, China). All cell lines were authenticated using short tandem repeat (STR) profiling and confirmed to be free of mycoplasma contamination before experimentation. Cell maintenance was performed in RPMI 1,640 medium supplemented with 10% fetal bovine serum (FBS; Vazyme, Nanjing, China). Cultures were kept in a humidified incubator at 37 °C with a 5% CO₂ atmosphere.

siRNA targeting specific genes were synthesized by Beijing Tsingke Biotechnology Co., Ltd. (see Table S1 for sequences). Transfection was performed using Lipo8000™ Transfection Reagent (Beyotime, Jiangsu, China) according to the manufacturer’s protocol. Lentivirus infection was screened with puromycin (1 μg/ml).

RNA extraction and quantitative real-time PCR

Total RNA was isolated from LUAD cells and BEAS-2B cells using an RNA extraction kit (Beyotime, Jiangsu, China). With 1 μg of total RNA as template, reverse transcription into complementary DNA (cDNA) was carried out using HiScript III RT SuperMix for qPCR (Vazyme, Nanjing, China). The synthesized cDNA was then amplified using the same qPCR SuperMix from Vazyme. For normalization of gene expression, β-actin served as the internal control. Relative RNA levels of the target genes were determined via the 2^−ΔΔCT method. All experiments were run in triplicate, and primer sequences are provided in Table S2.

Flow cytometry

Cells were resuspended in pre-cooled Staining Buffer (Proteintech, Wuhan, China) and adjusted to 10⁶ cells/mL. A 100 μL aliquot of the cell suspension was transferred to polypropylene tubes, followed by adding an appropriate volume of fluorochrome-labeled antibody to each tube; the cell-antibody mixture was then incubated for 30 min at room temperature in the dark. After incubation, cells were washed twice with Staining Buffer via centrifugation (1,200 rpm, 5 min at 4 °C). Fluorescence signals were detected using a Beckman CytoFLEX S flow cytometer, and subsequent data were analyzed and visualized using FlowJo software.

Mitotracker, TMRM, ROS and Hoechst staining

Cells were stained with Mitotracker (Invitrogen, Waltham, MA, USA), Tetramethylrhodamine (TMRM) (Thermo Fisher Scientific, Waltham, MA, USA), ROS (Beyotime, Jiangsu, China) and Hoechst 33342 (Yeasen, Shanghai, China) yielding a final concentration of 200 nM for 30 min, 100 nM for 30 min, 10 µM for 20 min, 0.5 µg/mL for 10 min before live cell imaging or flow cytometry respectively. The cells were then collected and washed with PBS three times, followed by FACS or fluorescence microscope analysis. Fluorescence in live cells was captured by confocal microscopy or plotted using the FlowJo software. Perform at least three biological replicates for each condition.

mtDNA localization

Cells were stained with Mitotracker (Invitrogen, Waltham, MA, USA), PicoGreen (Yeasen, Shanghai, China), and Hoechst 33342 (Yeasen, Shanghai, China) at final concentrations of 200 nM for 30 min, 3 µL/mL for 1 h, and 0.5 µg/mL for 10 min, respectively. Wash with PBS three times and then analyze with confocal microscopy. Perform at least three biological replicates in each case.

Quantitative PCR and mtDNA copy number analysis

Extract DNA components from cytoplasm using a DNA extraction kit, and detect the expression levels of mtDNA markers ND1 using qPCR; Use TERT as internal reference. The primer sequences are listed in Table S2.

Statistical analysis

All statistical analyses were performed utilizing R statistical software (version 4.3.3; R Foundation for Statistical Computing) and GraphPad Prism 8 (GraphPad Software Inc., San Diego, CA, USA). Continuous variables with normal distribution were compared between groups using Student’s independent t-test, while categorical variables were analyzed through Pearson’s chi-square test. Nonparametric Spearman’s rank correlation coefficient was employed to assess bivariate associations between continuous variables. Survival outcomes were evaluated through Kaplan-Meier methodology with log-rank testing for univariate analysis, followed by multivariable Cox proportional hazards regression modeling to adjust for potential confounding factors. All statistical tests were conducted using two-tailed comparisons, with a probability threshold of P < 0.05 defining statistical significance.

Analytical workflow for mitochondrial-related prognostic model construction in LUAD

The general workflow employed in this study was illustrated in Fig. 1. Our workflow began with processing RNA-seq data from the TCGA-LUAD cohort (n = 550), and identifying 3,358 differentially expressed genes (DEGs) with statistical significance. By intersecting these DEGs with mitochondrial genes in MitoCarta 3.0 (n = 1,136), 71 mitochondrial-associated DEGs were obtained. Next, we developed a prognostic signature using univariate Cox, LASSO regression, and multivariate Cox analyses. Subsequently, we constructed a clinical nomogram integrating the risk score and clinicopathological factors. Model validation was performed via survival analysis, ROC curves, and external GEO datasets (GSE31210 and GSE30219). Functional enrichment and TME analyses were conducted in risk-stratified subgroups. Finally, mitochondrial functions of identified biomarkers were experimentally validated.

Development of a mitochondrial-related risk model in LUAD

Through rigorous quality control of TCGA-derived LUAD cohorts (n = 550), we established an analytical cohort comprising 478 tumor samples and 58 adjacent normal tissues with complete clinical annotations and minimum 30-day survival data. Systematic differential expression analysis (|log2FC| > 1.5, adj. P < 0.01) identified 3,358 protein-coding DEGs, including 2,214 upregulated and 1,134 downregulated transcripts in LUAD vs normal controls (Fig. 2A). Bioinformatic integration with 1,136 mitochondrial genes from MitoCarta3.0 database revealed 71 putative mitochondrial-associated DEGs in LUAD (Fig. 2B). Univariate Cox proportional hazards regression (P < 0.01 threshold) identified 14 DEGs with prognostic significance for overall survival (OS) (Fig. 2C). Subsequent LASSO regression and multivariate Cox analysis refined the prognostic signature to three core mitochondrial-related (MTR) genes: SFXN1, CPS1, and MTFR2 (Figs. 2D–2F, Table 1).

Transcriptomic analysis confirmed elevated mRNA levels of SFXN1, CPS1, and MTFR2 in LUAD (Figs. 2G–2I). Immunohistochemical analysis based on HPA datasets revealed concordant protein overexpression patterns in LUAD specimens compared to adjacent normal controls (Figs. 2J–2L). Kaplan-Meier survival curves with log-rank testing demonstrated significantly reduced OS in LUAD patients with higher MTFR2 (log-rank P < 0.001), or with higher SFXN1 (log-rank P = 0.005) (Figs. S1A, S1B). Furthermore, although the traditional median cutoff failed to demonstrate the prognostic value of CPS1 (log-rank P = 0.294), patients with CPS1 expression in the top 20% had significantly shorter overall survival than those in the bottom 20% (log-rank P = 0.018). Additionally, we used X-tile (3.6.1) software to objectively determine the optimal cutoff value for CPS1 expression (CPS1 = 7.33) to stratify patients into high and low expression groups, and the log-rank test showed a P value of less than 0.001 (Figs. S1C–S1E).

Patients were stratified into high-risk and low-risk subgroups based on the median risk score, which in both training and validation cohorts was computed based on the following formula:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=0.3067\ast \mathrm{S}\mathrm{F}\mathrm{X}\mathrm{N}1+0.0798\ast \mathrm{C}\mathrm{P}\mathrm{S}1+0.17\ast \mathrm{M}\mathrm{T}\mathrm{F}\mathrm{R}2.$$

Figs. 3A–3D presents an integrated visualization of clinical-molecular correlations, featuring: (a) risk score distribution, (b) survival status/time alignment, (c) risk stratification patterns, and (d) heatmap of the three-gene expression profile. Kaplan-Meier analysis confirmed significantly shorter OS in high-risk patients (log-rank P < 0.001; Fig. 3E). Time-dependent ROC analysis demonstrated robust accuracy for 1-year (AUC = 0.712), 3-year (AUC = 0.647), and 5-year (AUC = 0.615) survival prediction (Fig. 3F). External validation in GSE31210 and GSE30219 cohorts consistently supported the prognostic stratification (P < 0.001, P = 0.001; Figs. S3A–S3E). The model maintained predictive performance in validation sets, with GSE31210 showing AUCs of 0.652 (1-year), 0.600 (3-year), and 0.725 (5-year) (Fig. S2F), and GSE30219 achieving AUCs of 0.827 (1-year), 0.755 (3-year), and 0.784 (5-year) (Fig. S3F). These multi-cohort validations confirm the clinical generalizability of our MTR risk stratification framework.

Clinical correlation and prognostic nomogram development

Clinical correlation analysis revealed significant relationships between risk scores and key demographic-clinical parameters, which were visualized using circular visualization and bar charts (Fig. 4A, Figs. S4A–S4F, Table 2). Male patients exhibited higher risk scores than females (Mann-Whitney P < 0.001), and advanced-stage (III–IV) tumors were markedly enriched in the high-risk subgroup ( ${\chi }^2$ P < 0.001). However, no significant age-dependent risk stratification was observed (P = 0.734). Univariable and multivariable Cox regression adjusted for TNM stage, age, and gender confirmed the risk score as an independent prognostic indicator (HR = 1.476, 95% CI [1.088–2.002]; P = 0.0123), outperforming conventional clinical parameters (Fig. S4J, Fig. 4B, Tables 3, 4). We developed an innovative prognostic nomogram integrating the MTR risk signature with clinicopathological variables, demonstrating strong concordance between predicted and observed survival probabilities at 1-, 3-, and 5-year intervals (Figs. 4C, 4D). The model achieved a C-index of 0.680 (95% CI [0.656–0.703]), signifying robust clinical applicability.

Functional enrichment analysis of DEGs in high risk and low risk groups

Comparative transcriptomic profiling between high-risk and low-risk subgroups identified 600 DEGs meeting stringent thresholds (|log2FC| > 1.5, Benjamini-Hochberg adj. P < 0.01), comprising 387 upregulated and 213 downregulated transcripts in high-risk patients (Figs. 5A, 5B). Functional enrichment analysis demonstrated striking pathway polarization. GO analysis revealed high-risk-specific enrichment in fibrinogen complex and extracellular region at the cellular component level. Molecular function terms predominantly mapped to receptor regulator activity and signal release, while biological process centered on regulation of hormone levels and receptor ligand activity. KEGG pathway analysis identified eight significantly enriched pathways, including Neuroactive ligand-receptor interaction, complement and coagulation cascades, GABAergic synapse, vascular smooth muscle contraction, steroid hormone biosynthesis, biosynthesis of amino acids, NOD-like receptor signaling pathway and Toll-like receptor signaling pathway (Figs. 5C, 5D). GSEA results showed that the high-risk status was significantly associated with extracellular membrane bounded organelle, fibrinogen complex, mitochondrial protein containing complex, outer mitochondrial membrane protein complex and immune related pathways (Figs. S5A–S5D).

Association between MTR risk scores and immune landscape in LUAD

We performed bioinformatics analysis to determine the effects of mitochondrial function on TME. The high-risk group exhibited higher tumor purity but lower stromal score, immune score, and Estimate score than the low-risk group. Moreover, the MTR risk score was significantly associated with Tumor Immune Dysfunction and Exclusion (TIDE) scores (Fig. 6A), suggesting impaired immune surveillance in high-risk tumors. To further confirm the correlation of MTR signatures expression with the immune microenvironment, the proportion of tumor infiltrating immune subsets was analyzed using CIBERSORT algorithm, and the immune cell profiles in LUAD samples were constructed.

We found through rainbow plots, CIBERSORT-based bulk immune cell infiltration analysis, and single immune cell correlation analysis respectively that levels of activated CD4+ T cells, memory B cells, effector memory CD4+ T cells, and Type 2 T helper cells were elevated in the high-risk group (Figs. 6B, 6C, Fig. S6). By contrast, the risk score was negatively associated with the abundance of innate immune cells, including macrophages, activated B cells, activated dendritic cells, central memory CD4+ T cells, Natural killer T cells, and Type 17 helper cells (Figs. 6B, 6C, Fig. S7). A heatmap was generated to display the correlation between 21 kinds of tumor-infiltrating immune cells. The numeric value in each tiny box represented the P value of correlation between two types of cells (Fig. S8).

Mutation status and immune checkpoint markers of LUAD patients in high risk and low risk groups

As shown in Figs. 7A, 7B, the top 20 most frequently mutated genes in high-risk and low-risk group were presented. TP53, TTN, MUC16, CSMD3, RYR2, ZFHX4, LRP1B, USH2A, KRAS, SPTA1, XIRP2, FLG and NAV3 were the common high-frequency mutation genes in both groups (Figs. 7A, 7B).

Importantly, the high-risk group showed remarkedly elevated expression of immune checkpoint markers such as CD274, CD276, PDCD1, and RELT, indicating a close relationship between MTR index and immune therapy tolerance (Fig. 7C).

Experimental verification of MTR genes expression and mitochondrial functions in LUAD

Previously, through bioinformatics analysis, we initially screened and identified three key MTR genes—SFXN1, CPS1, and MTFR2—that are potentially associated with the occurrence and development of LUAD, immune regulation, and mitochondrial function. However, the accuracy of bioinformatics predictions still needs further confirmation through experimental verification. Moreover, there is a lack of direct experimental evidence regarding how these genes specifically regulate mitochondrial morphology, function, and genomic stability in LUAD cells. To substantiate our bioinformatic findings, we systematically evaluated MTR gene expression and mitochondrial functional alterations in LUAD models. Quantitative RT-PCR analysis revealed consistent upregulation of SFXN1, CPS1, and MTFR2 across malignant cell lines (A549, H23, H1975, H1299, PC9) compared to normal bronchial epithelium (BEAS-2B) (Figs. 8A–8C). To elucidate the functional consequences of MTR gene dysregulation, we performed siRNA-mediated knockdown followed by comprehensive mitochondrial phenotyping (70–85% efficiency at 72 h post-transfection, Figs. 8D–8I, Fig. S9). CPS1 knockdown induced fission-dominant morphology with reduced ΔΨm↓, but elevated ROS and enhanced mtDNA extrusion. SFXN1 knockdown resulted in mitochondrial fusion, ΔΨm↑, and ROS suppression. MTFR2 knockdown increased mitochondrial fusion, ΔΨm, and ROS accumulation. All perturbations significantly reduced the copy number of mtDNA, with CPS1 knockdown uniquely promoting cytosolic mtDNA release (Figs. 9, 10). These results demonstrate gene-specific modulation of mitochondrial dynamics and genome maintenance.