Identification of metabolic driver genes and targeted drug screening for lung adenocarcinoma metastasis at the single-cell resolution

Introduction

Lung cancer is one of the most common types of cancer (1). Among all cancers globally, lung cancer stands out as particularly concerning due to its persistently rising incidence and mortality figures (2). In lung cancer, lung adenocarcinoma (LUAD) is the main subtype, accounting for about 40% of all lung cancer cases (3). LUAD carries a high risk of distant metastasis across all disease stage, with metastasis being a main cause of mortality and poor prognosis (4). Metabolic activation drivers LUAD metastasis, making it critical to decode the processes involved in metastasis and optimize prognosis prediction from a metabolic perspective.

Energy metabolism is the basis of tumor cell proliferation, invasion, and metastasis. Cancer cells need enough nutrients to adapt properly to the tumor microenvironment (TME) and create a favorable metabolic environment for themselves (5). Therefore, targeting tumor cell metabolism has become a strategy to explore new strategies for cancer therapy (6). However, progress in targeting cancer metabolism therapeutically in the past decade has been limited (7), due to metabolic compartmentalization and heterogeneity that exist in tumors. One of the difficulties in LUAD treatment is the heterogeneity in TME (8), which highlights the necessity of identifying different cellular and molecular factors of LUAD at single-cell level.

In this study, we employed single-cell RNA-sequencing (scRNA-seq) data analysis to reveal that metabolic alterations in tumor cells and their microenvironment are closely associated with the metastasis of LUAD. Simultaneously, we identified 43 pivotal metabolic driver genes (MDGs) and screened 307 drug molecules, which provides potential targets and candidate drugs for exploring therapeutic strategies against metastatic LUAD. We present this article in accordance with the MDAR reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-484/rc).

Methods

Data resource download

We obtained the single-cell data of LUAD from GSE131907 cohort study published by Kim et al., which encompasses the single-cell data of 44 patients (9). The quality control standards are detailed in the original literature (9). STAR-Counts data and clinical data for The Cancer Genome Atlas (TCGA)-LUAD were downloaded using the R package TCGAbiolinks (v2.28.3) (10). We retrieved lung cancer datasets with survival data from the Lung Cancer Explorer (LCE) database (11). Subsequently, the LUAD cohort datasets (GSE26939, GSE68465, GSE31210, and GSE50081) and a lung cancer cohort dataset GSE30219 were downloaded using the software package GEOquery (v2.64.2). The sample information for these datasets can be found in table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-1.xlsx. Additionally, we obtained differentially expressed gene (DEG) analysis results from RNA-sequencing (RNA-seq) data of LUAD 3D tumor spheroid models, which were derived from a study by De Vitis et al. (12). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Annotation, clustering, and dimensionality reduction of single cells

For the cell annotation of the single-cell dataset GSE131907, we utilized the annotation results initially provided by the dataset authors (9). Subsequently, we employed t-distributed stochastic neighbor embedding (t-SNE), implemented in the R package Seurat (v4.3.0), for data visualization (13).

Calculation of single-cell metabolic activity scores

The scMetabolism (v0.2.1) is an R tool for quantifying metabolism activity at the single-cell resolution (14). We used the sc.metabolism.seurat function to assess the activity scores of 76 metabolic pathways at the single-cell level. During the analysis, the method parameter was set to “VISION”. Subsequently, the data were visualized by DimPlot.metabolism and DotPlot.metabolism function.

Differential expression analysis at the single-cell level

To explore the expression differences of cell types in normal and disease tissues, we used the FindMarkers function from Seurat package (v4.3.0) to identify changes in expression of designated cells from different tissue sources at the gene level (13). Additionally, we performed an average expression analysis of specific types using the AverageExpression function. Genes with an absolute value of log₂fold change (log₂FC) greater than 1 and an adjusted P value less than 0.05 are considered to be differentially expressed.

Immune cell infiltration analysis

Single-sample gene set enrichment analysis (ssGSEA) (15) was utilized to assess the proportions of 28 tumor-infiltrating immune cells (TIICs) phenotypes within the TME.

Construction of the prognostic risk model and survival analysis

Using the TCGA-LUAD dataset, a univariate Cox regression analysis was performed to identify 43 MDGs with prognostic significance based on survival package (v3.3.1) in R. This analysis revealed 14 MDGs that could potentially serve as risk factors. To construct a potential independent prognostic model for LUAD based on MDGs, the 14 identified prognostic MDGs were further incorporated in a multivariate Cox regression calculation analysis. To mitigate multicollinearity, stepwise regression was employed to determine the optimal combination of prognostic factors. After establishing a formula to calculate the risk score, the risk score for each case was computed as follows:

Riskscore=∑coefi×Ei[1]

Here, coef and E represent the correlation coefficient and expression value of gene i, respectively. Using the median risk score of the TCGA-LUAD cohort as a cutoff value, patients in external datasets were stratified into high- and low-risk groups. Receiver operating characteristic (ROC) curves were generated using the timeROC R package (v0.4). The prognostic efficiency of the model was assessed by calculating the area under the curve (AUC) from the ROC analysis. Additionally, overall survival (OS) for specific groups was conducted the survival package (v3.3.1) in R.

Use of the non-negative matrix factorization (NMF) algorithm to identify molecular subtypes

A NMF clustering algorithm was used to cluster the patients based on the expression profiles of MDGs in the TCGA-LUAD cohort. The NMF (v0.26) package was employed to compute the average contour breadth of the communal member matrix, and a threshold of 10 members was set as the minimum for each of the subtypes. For the NMF methodology, we opted for the standard “brunet” configuration and executed 50 iterations. The range of clusters (k) was specified as 2 to 10. The dispersion, cophenetic, and silhouette indicators were used to determine the optimal clustering number, and the optimal clustering number was determined to be 3.

Differential gene expression and gene enrichment analysis

Differential expression analysis of genes was performed using the DESeq2 (v1.42.1) (16) package in R (v4.1.0). Genes with a |log₂FC| ≥1.5 and a false discovery rate (FDR) <0.01 were considered DEGs. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were conducted using the clusterProfiler package (v4.0.2) (17) in R (v4.1.0).

Drug screening

Drug and target data were downloaded from DrugBank database (https://go.drugbank.com/), a comprehensive resource that provides information on drugs and their targets. Subsequently, we utilized these data to screen drugs that target MDGs and compiled the relevant information. These operations were performed using the R software (v4.1.0).

Cell culture

The human LUAD cell line A549 was purchased from Wuhan Punosai Biotechnology Co., Ltd. (Wuhan, China). This cell was cultured in Dulbecco’s modified Eagle medium (DMEM) complete medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin-gentamicin solution. The cells were maintained in a cell culture incubator at 37 ℃ with 5% CO₂.

Cell wound healing assay

To assess the migration ability of LUAD cells treated with targeted drugs, a scratch assay was conducted. Cells were grown to confluence in six-well plates. Subsequently, cells were wounded with a pipette tip and treated with cholic acid (100, 150, 205 µM; Cat. No. HY-N0324; MCE, Monmouth Junction, NJ, USA) alone or dimethyl sulfoxide (DMSO). Cells were examined by light microscopy prior to the addition of experimental treatments (0 h) and at 24 and 48 h after treatment. The wound area was measured at each time point and the percent closure at 24 and 48 h was calculated. The area of each scratch wound was determined using ImageJ software to quantify the differences between 0, 24, and 48 h. The formula for calculating the wound healing rate (Shr) is as follows: Shrt(%)=Isw−FswtIsw×100%, where Shr_t is the wound healing rate at time point t, Isw is the area of the scratch at 0 h (initial area), Fsw_t is the mean distance between cells at time point t. Three wells per experimental treatment and three wounds per well were examined. Results reported are the mean ± standard error (SE).

Transwell migratory assay

The A549 cells (5×10⁴) were added into the upper chamber of the Corning Incorporated Transwell chamber (8 µm; Corning Inc., Corning, NY, USA) and DMEM containing 10% FBS was added to the lower chamber. The cells were then treated with either cholic acid (150 µM; Cat. No. HY-N0324; MCE) or DMSO as a control. After culturing for 24 h, the cells were fixed and stained with 0.5% crystal violet. The cells in randomly selected fields were photographed under an inverted microscope and then were roughly counted using ImageJ software.

Statistical analysis

Continuous variables that followed the normal distribution were compared using independent t-tests, while continuous variables with skewed distributions were compared with the Mann-Whitney U test. Survival curves for categorical variable prognostic analyses were generated using the Kaplan-Meier method, and the log-rank test was used to estimate statistical significance. The significance level was set at P<0.05, and all statistical tests were two-sided. All statistical data analyses were performed using R software (v4.1.0).

Results

Metabolic landscape in early, advanced, and metastatic LUAD

Processed scRNA-seq data and meta information were obtained from the original article (9). A total of 208,506 cells from 44 patients with treatment-naïve LUAD during endobronchial ultrasound/bronchoscopy biopsy or surgical resection were classified into nine distinct cell types, and these cells were visualized t using t-SNE plots (Figure 1A-1C). To understand the metabolic landscape of these cells, we first computed the metabolic score of 76 metabolic pathways in these cells. We found that most metabolic pathways exhibited high activation scores in the brain metastases (mBrain) tissues relative to other tissues (Figure 1D). In addition, we also noted that cells derived from normal lymph nodes (nLN) tissue exhibited low levels of activation scores in nearly all metabolic pathways (Figure 1D). Cells in metastatic lymph nodes (mLN) tissue had higher activation scores for many metabolic pathways compared with those in nLN (Figure 1D). Similarly, cells derived from advanced-stage lung cancers (tL/B) tissue also exhibited higher metabolic activation scores in most metabolic pathways relative to pleural fluids (PE) (Figure 1D). This phenomenon implied that the activation of metabolism may be associated with the progression and metastasis of LUAD. Additionally, we explored the activation of various metabolic pathways in different cell types. We found that epithelial cells, fibroblasts, myeloid cells, and oligodendrocytes are the major cell types involved in the activation of most metabolic pathways, whereas B/T lymphocytes and natural killer (NK) cells exhibit a very low fraction of metabolic activation relative to other cells (Figure 1E).

Figure 1 Metastatic LUAD single cell atlas. (A) Overview of tissue origins in the GSE131907 cohort study collection, including tLung, tL/B, PE, mLN, mBrain, nLung, and nLN. The red arrow indicates tumor metastasis from the primary locus to the brain parenchyma. (B,C) t-SNE plots illustrated the clustering results of the different origins (B) and the major cell lineages (C). (D,E) The average activation score of 76 metabolic related pathways from different samples (D) and cell subpopulations (E). LUAD, lung adenocarcinoma; mBrain, brain metastases; mLN, metastatic lymph nodes; nLN, normal lymph nodes; nLung, normal lung tissues; PE, pleural fluids; t-SNE, t-distributed stochastic neighbor embedding; tL/B, advanced-stage lung cancers; tLung, early-stage lung cancers.

Given the variations in cellular proportions derived from different tissues (Figure S1A-S1D), we focused our attention on elucidating the metabolic activation disparities within the same cell type but of diverse tissue origins, with the aim of better understanding the metabolic drivers involved in the metastatic process of LUAD. First, we conducted a comparative analysis of the metabolic activation disparities in immune cells between normal samples and metastatic samples of the same tissue type. We found that the metabolism of the four T lymphocytes subpopulations in normal lung tissues (nLung), early-stage lung cancers (tLung), and tL/B tissues displayed different forms (Figure S2A-S2D). Similarly, compared to nLung, the metabolic activation of B lymphocytes, NK cells, mast cells, myeloid cells, and their subpopulations in tLung and tL/B also exhibited significant differences (Figures S3-S6). Consistently, B lymphocytes, T lymphocytes, myeloid cells, and their subpopulations exhibited marked differences in most metabolic pathways between mLN and nLN (Figures S7-S9). Most metabolic pathways were differentially activated in B lymphocytes, T lymphocytes, myeloid cells, and NK cells, and their subpopulations from mLN compared to nLung (Figures S10-S13). Furthermore, T lymphocytes, B lymphocytes, NK cells, mast cells, myeloid cells, and their subpopulations from mBrain tissues showed significant differences in metabolic pathway activation scores versus nLung (Figures S14-S18). This collective evidence indicates that metabolic alterations in immune cells drive LUAD progression and metastasis. Next, we also assessed metabolic alterations in epithelial cells. Consistent with immune cells, epithelial cells metabolism differed significantly between nLung and tumor tissues (tL/B and tLung), and metastatic tissues (mBrain and mLN) (Figures S19-S21). Considering normal epithelial cells coexist in tumor tissues, we evaluated metabolic differences in epithelial cell subpopulations. Tumor cells [including transcriptional states (tS1, tS2, and tS3) and malignant cells] showed distinct metabolic profiles from normal epithelial cells [including alveolar types I (AT1) and II (AT2), and club cells and ciliated cells] (Figures S19-S21).

Taken together, this evidence suggests that the heterogeneity of cellular metabolism in the TME affects the progression and metastasis of LUAD.

Identification of MDGs in LUAD

To explore the potential of targeting metabolism-related genes for LUAD treatment, we first performed differential expression analysis to characterize gene expression changes in three original sites (nLung, tLung, and mBrain). For epithelial cells, compared with nLung, 232 and 53 DEGs were identified in mBrain and tLung, respectively (table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx); Compared to tLung, 95 DEGs were detected in mBrain tissue (table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx). In B lymphocytes, 96 and 43 DEGs were identified in mBrain and tLung versus nLung, respectively, while 42 genes were differentially expressed between mBrain and tLung (table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx). In T lymphocytes, 35 and 15 DEGs were found in mBrain and tLung compared to nLung, respectively, while 17 DEGs between in mBrain and tLung (table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx). Similarly, in myeloid cells, 88 and 21 DEGs were identified in mBrain and tLung versus nLung, respectively, while 15 genes were differentially expressed in mBrain and tLung (table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx). Subsequently, we focused on metabolism-related DEGs. As shown in Figure 2A, 43 DEGs were associated with metabolism, primarily originating from epithelial cells. Specifically, 28 metabolism-related DEGs identified in epithelial cells, with 26 genes from the mBrain versus nLung comparison (Figure 2A-2C, table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-2.xlsx). Consistently, expression changes of most genes aligned with TCGA-LUAD cohort (Figure 2D). Considering single-cell level heterogeneity, we analyzed these DEGs in epithelial cell subpopulations and found 14 genes differentially expressed between malignant and normal epithelial cells (Figure 2E). Additionally, we analyzed metabolism-related DEGs in B lymphocyte and myeloid cell subpopulations, revealing significant expression differences across samples but heterogeneity within subpopulations (Figure 2F,2G). A previous study using 3D tumor spheroid models of metastatic LUAD cultured under nutrient stress conditions demonstrated metabolic reprogramming in metastasis. Several metabolism-related genes (MSMO1, GAPDH, IDI1, DHCR24, CDA, PLA2G7, ALDH2, GPD1, GPX3) were differentially expressed in this model (Figure S22), supporting their role as metabolic drivers. In conclusion, through comprehensive differential expression analysis across cell types and tissue sites, we designated 43 metabolism-related DEGs as MDGs. These genes are primarily derived from epithelial cells and represent potential targets for investigating metabolic drivers in LUAD metastasis.

Figure 2 Identification of MDGs in LUAD. (A) Upset and Venn plots show the relationship between the DEGs of cell subtypes of the different origins and the metabolism related genes. (B) Heat map showing the expression of 28 MDGs of epithelial cells in nLung, tLung, and mBrain origins. (C) The upset diagram shows the source of 28 MDGs. (D) The expression levels of 28 MDGs in LUAD cohort based on the TCGA database. (E) Heat map showing the expression of 28 MDGs in epithelial cell subtypes. (F,G) Heat map showing the expression of MDGs in subtypes of B lymphocytes (F) and myeloid cells (G). ns, no significant difference (P>0.05); *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. AT, alveolar type; DCs, dendritic cells; DEGs, differentially expressed genes; LCs, Langerhans cells; LUAD, lung adenocarcinoma; Mac, macrophage; mBrain, brain metastases; MDGs, metabolic driver genes; mo-Mac, monocyte-derived macrophage; nLung, normal lung tissues; pDCs, plasmacytoid dendritic cells; TCGA, The Cancer Genome Atlas; tLung, early-stage lung cancers; TPM, transcripts per million.

Hazard scoring system based on MDGs

Based on the identified MDGs and observed metabolic alterations in diverse cell types across tissue, which collectively demonstrate their potential influence on LUAD progression and metastasis, we next developed a Cox-based hazard scoring system to investigate the clinical significance and predictive value of these gene expression changes. Following feature selection from the 43 identified MDGs, a Cox proportional hazards model comprising 14 MDGs was constructed, as outlined below:

Riskscore=0.1528×FBP1+0.3218×GAPDH+0.1677×ALDH3B1+0.2386×B4GALT1+0.1098×CA2+0.2639×ACOT7−0.4104×PNPLA6−0.162×GPX3+0.23×NME4+0.1427×PDE4D+0.1974×UPP1−0.1407×MGST1−0.1216×NNMT+0.181×HMOX1[2]

Subsequently, we calculated the patients’ risk score based on the model. The patients’ risk score distribution, survival status, and expression levels of the 14 MDGs were presented in Figure 3A-3C. Then we performed survival analysis between the high-risk group and low-risk group, showing that the OS of the low-risk group was better than that of the high-risk group (Figure 3D). The AUCs for 1-, 3-, and 5-year OS were 0.710, 0.723, and 0.683 in the TCGA-LUAD cohort (Figure 3E). To validate the predictive value of the model, risk scores for patients in GSE26939, GSE68465, GSE31210, and GSE50081 cohorts were calculated with the same formula. Consistently, the OS of the high-risk group was worse than that of the low-risk group (Figure 3F). Additionally, we found that our model was well-validated in the unstratified lung cancer cohort GSE30219 (Figure 3G). Collectively, our research results indicate that the MDG signature performed well for OS prediction.

Figure 3 Hazard scoring system based on MDGs. (A) Distribution of hazard scores of patients with LUAD in training cohort based on MDGs. Each data point represents an individual patient. (B) Training cohort scatter plots showing the association between the OS and the risk score in LUAD patients according to prognostic features of MDGs. Each data point represents an individual patient. (C) Training cohort heatmap showing fourteen unfavorable genes in high-risk patients, in contrast to the three favorable genes. Each column represents an individual patient. (D) Kaplan-Meier analysis of OS in the TCGA-LUAD cohort based on the risk scores of the 14-MDG signature. (E) ROC curves for 1-, 3-, 5-year OS prediction of the 14-MDG signature in the TCGA-LUAD cohort. (F) Kaplan-Meier analysis of OS in the LUAD validation cohorts of GSE26939, GSE68465, GSE31210, and GSE50081 based on the risk scores of the 14 MDGs signature. (G) Kaplan-Meier analysis of OS in the lung cancer cohorts GSE30219 based on the risk scores of the 14 MDGs signature. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LUAD, lung adenocarcinoma; MDGs, metabolic driver genes; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Molecular subtype based on MDGs

Building upon the identification of 43 MDGs and their potential as targets for investigating metabolic drivers in LUAD metastasis, we further explored the molecular subtypes based on these MDGs to gain deeper insights into LUAD heterogeneity. We employed the NMF algorithm to define a novel molecular subtype classification based on MDGs. Using cophenetic, dispersion, and silhouette indicators, an optimal clustering number of three was selected (Figure 4A, Figure S23). Principal component analysis confirmed clear separation among the three subtypes (Figure 4B). To assess clinical difference among subtypes, we analyzed their prognostic characteristics. The results showed that subtype 3 (C3) had the best prognosis, while C1 and C2 had the worst, with significant prognostic differences among the three subtypes (Figure 4C).

Figure 4 Molecular subtypes based on MDGs. (A) Consensus matrix based on NMF algorithm. (B) TCGA-LUAD populations identified after unsupervised clustering in (A). (C) Kaplan-Meier survival plot of C1, C2, and C3 groups. A log-rank test was conducted. LUAD, lung adenocarcinoma; MDGs, metabolic driver genes; NMF, non-negative matrix factorization; PC, principal component; TCGA, The Cancer Genome Atlas.

Immune infiltration, metabolism, and function enrichment of different molecular subtypes

Understanding the differences in immune infiltration, metabolic profiles, and functional enrichments among these subtypes can provide further insights into disease progression and prognosis. Thus, we first analyzed DEGs among the three molecular subtypes (Figure 5A-5C). We obtained a total of 153 DEGs among the subtypes (Figure 5D). GO biological process enrichment analysis showed that these DEGs were significantly enriched in metabolism-related pathways, such as olefinic compound metabolic process and retinol metabolic process (Figure 5E). Consistently, KEGG enrichment analysis also identified enrichment in metabolism-related pathways, including retinol metabolism (Figure 5F). These results indicate that tumor cell metabolic heterogeneity is associated with molecular subtypes. Subsequently, we evaluated metabolic enrichment scores across molecular subtypes and found that significant differences in most metabolic pathways (Figure S24). Finally, considering immune microenvironment heterogeneity, we assessed immune cell infiltration differences among subtypes. Most immune cells varied significantly across the three molecular subtypes, with C3 showing higher infiltration levels than C1 and C2 (Figure 5G). This observation may be associated with the better OS observed in the C3 subtype.

Figure 5 Immune infiltration, metabolism, and function enrichment of different molecular subtypes. (A-C) Volcano plot comparing DEGs between C1 and C2 (A), C3 versus C1 (B) and C3 versus C2 (C). (D) Venn plot constructed based on the DEGs presented in (A). (E,F) Enrichment analysis was performed to investigate GO-BP (E) and KEGG (F) categories of the DEGs. (G) The ssGSEA score for 28 TIICs across LUAD molecular subtypes. BP, biological process; DEGs, differentially expressed genes; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma; ssGSEA, single-sample gene set enrichment analysis; TIICs, tumor-infiltrating immune cells.

Drug screening targeting MDGs

The prognostic risk models and molecular subtypes established based on MDGs highlight the significant potential of these genes in LUAD metastasis therapy. We screened drugs targeting MDGs to explore therapeutic strategies for metabolism-related genes. To identify MDG-targeting drugs, we collected drug-target information from the DrugBank database and found that 34 MDGs were included in the DrugBank target pool (Figure 6A). Next, we summarized the corresponding targets-drugs information. Table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-3.xlsx lists the relationship between 307 drug molecules and these targets. We then counted the number of drug molecules corresponding to the targets. Among the down-regulated genes, the number of drug molecules targeting CA2 was the most, while among the up-regulated targets, the number of drug molecules targeting AKR1B1 was the most (Figure 6B,6C). Additionally, we also counted the number of MDG targets of these drugs. A total of 18 drug small molecules had at least two targets (Figure 6D, table available at https://cdn.amegroups.cn/static/public/tcr-2025-484-4.xlsx). Therefore, screened drugs targeting MDGs from these small molecule compounds may facilitate the development of therapeutic strategies for metastatic LUAD.

Figure 6 Drug screening targeting metabolism-related genes. (A) Venn plot shows MDGs in target pool of DrugBank database. (B,C) Statistical histogram shows the number of drugs in down regulated (B) and up-regulated MDGs (C). (D) Drugs targeting at least two MDGs. DEGs, differentially expressed genes; DEMGs, differentially expressed and methylated genes; MDGs, metabolic driver genes; NADH, nicotinamide adenine dinucleotide.

Effects of drugs targeting MDGs on cell metastasis

To evaluate the effect of screened compounds on LUAD metastasis, we treated A549 cells with different concentrations of cholic acid and performed wound healing assays and Transwell assays. The wound healing assay results clearly demonstrated that cholic acid significantly inhibited A549 cells migration. Specifically, compared with the untreated control, A549 cell migration showed significant concentration-dependent inhibition after 24 and 48 h of cholic acid treatment (Figure 7A,7B). To further validate the impact of cholic acid on metastatic ability, we used the Transwell assay to assess cells invasion. Results showed that as cholic acid concentration increased, the number of A549 cells migrating through the Transwell chamber membrane gradually decreased, confirming its inhibitory on cell metastasis (Figure 7C,7D). Collectively, these findings indicate that cholic acid, a screened drug molecule, inhibits LUAD metastasis. This also highlights that the drugs we screened can serve as candidate molecules for the development of therapeutic strategies against metastatic LUAD.

Figure 7 Effects of cholic acid on A549 cell metastasis. (A) Representative images depicting scratch wound healing at 24 and 48 h post-treatment with 50, 150, and 250 µM concentrations of cholic acid, with DMSO serving as the control. The cell-free area has been delineated. The scale bar corresponds to 50 µm. (B) The cell-free area was quantified from the microscopic images utilizing ImageJ software and is expressed as the percentage of wound closure. (C,D) Crystal violet staining (C) and cell count statistics (D) were conducted to assess the migratory capacity of A549 cells at 24 h after treatment with 50, 150, and 250 µM concentrations of cholic acid, with DMSO serving as the control. DMSO, dimethyl sulfoxide.

Discussion

Altered metabolism not only affects tumor cell proliferation but also affect their metastasis (18). Tumor cells often reprogram nutrient utilization to adapt growth and metastasis, gaining survival within their microenvironment (19). Metabolic interactions between tumor cells and microenvironmental cells (e.g., immune cells and fibroblasts) support tumor proliferation, progression, metastasis, and immune escape (20). Thus, targeting cell metabolism has emerged as a promising approach for cancer therapy development (6). However, tumor and microenvironmental heterogeneity pose significant challenges to this strategy. scRNA-seq has been widely used in exploring the heterogeneity of tumor cells (21,22). By measuring cell transcription at the single-cell resolution, scRNA-seq provides an unprecedented new view of the complexity of the TME (23). In this study, we characterized the metabolic landscape of LUAD across early- and advanced-stage at the single-cell level. Our findings may facilitate the discovery of new therapeutic targets or targeted drugs by leveraging LUAD metabolic heterogeneity.

By evaluating metabolic activation in LUAD across early to advanced stages, including primary and metastatic sites, we found that LUAD cell metabolism and that of microenvironmental cells changed during disease progression. This aligns with metabolic alterations being associated with tumor cell proliferation and metastasis (18). T/B lymphocyte, myeloid cell, NK cell, and their subpopulations exhibited distinct metabolic activation patterns across different disease sites. Understanding these differences is of great importance for developing LUAD therapies. In this study, we identified 43 MDGs in T lymphocytes, B lymphocytes, myeloid cells and epithelial cells from diverse disease sites. Leveraging these MDGs and observed metabolic alterations in various cell types across tissue, which influence LUAD progression and metastasis, we developed a hazard scoring system. Our model showed potential in evaluating LUAD patient prognosis. However, we acknowledged certain limitations in dataset size, which may require further model optimization for clinical predictive in future studies. Furthermore, by targeting the 43 MDGs in our analysis, we applied the NMF algorithm to define novel molecular subtypes, facilitating a comprehensive analysis of LUAD heterogeneity. Our stratified LUAD into three subtypes, whose further analysis revealed metabolic and immune infiltration heterogeneity. While survival analysis showed no statistically significant OS difference between C1 and C2—potentially confounded by uncharacterized factors—C3 exhibited superior OS compared to both C1 and C2. Additionally, DEGs among the three subtypes were significantly enriched in metabolic-related pathways, indicating that metabolism-based molecular subclassification has potential clinical significance.

Among these MDGs, we found that the CA2 gene product is the target of numerous drugs, followed by PDE4D. Previous studies have shown that CA2 inhibited liver cancer cells metastasis (24,25), but its role in LUAD remains unclear. Our analysis revealed that CA2 was downregulated in LUAD, which may be associated with LUAD metastasis. PDE4D has been reported to promote the growth of colorectal (26), lung (27), prostate (28), and pancreatic (29) cancer cells. However, our study showed that PDE4D was expressed at lower levels in LUAD, suggesting that PDE4D may have a more complex role in LUAD. Notably, studies have indicated that LUAD patients with low CD8⁺ T infiltration in anti-PD-1 therapy-insensitive group exhibited significant PDE4D upregulation (27,30). Collectively, these findings warrant further investigation into the roles of these genes in LUAD metastasis, which may provide insights for developing metabolism-related gene-targeted therapeutic strategies.

In this study, we sought to identify drug molecules targeting MDGs. From the DrugBank database, we screened 307 drug molecules, which included endogenous molecules and essential metal ions. Notably, nicotinamide adenine dinucleotide (NADH) targeted 10 MDGs. NADH supplementation has been shown to improve function in Alzheimer’s disease patients (31), reduce physical disability, and alleviate depression in Parkinson’s disease patients (32). Additionally, NADH regulates metabolic vulnerability in lung cancer (33), though the direct impact of nicotinamide adenine (NAD)-dependent pathway impairment in humans remains unclear due to spatiotemporal NAD/NADH heterogeneity in cells (34). We hypothesize this relates to tumor metabolic dysregulation. Thus, modulating these endogenous or exogenous molecules may improve LUAD prognosis. For example, excessive isoleucine supplementation inhibited lung cancer cell growth (35). Furthermore, some screened drugs have confirmed antitumor effects. Artemisinin, for instance, inhibits tumor proliferation and metastasis, induces apoptosis, and enhances immunity (36,37). In this study, cholic acid inhibited lung cancer cells migration, consistent with its reported suppression of H1299 cell proliferation (38). Oral cholic acid replacement is effective in children with primary bile acid synthesis defects (39). However, notably, cholic acid itself is a type of bile acid, different types of bile acids may promote the proliferation of lung cancer cells (38), suggesting cholic acid may benefit only specific LUAD patients. Considering the OS differences among our identified molecular subtypes, these drugs may be effective in specific subtype populations. Collectively, our analysis indicates that the investigating of these drugs can provide candidate molecules for the development of therapeutic strategies against metastatic LUAD.

Conclusions

Based on single-cell transcriptome data, we identified 43 metabolism-related genes driving LUAD metastasis, which have demonstrated robust performance in LUAD subtype classification and prognosis prediction. Additionally, we screened 307 drug molecules targeting these genes from the DrugBank database and experimentally validated that cholic acid inhibits migration of LUAD cells. Our study offers novel metabolism-related targets and candidate drug molecules for targeted therapy of metastatic LUAD.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-484/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-484/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-484/prf

Funding: This work was supported by grants from the High-Level Talent Initiation Project of Guizhou Medical University (No. J [2022] 051), the National Natural Science Foundation Cultivation Project of Guizhou Medical University (No. 22NSFCP34), the National Natural Science Foundation of China (Nos. 32160224, 62166009 and 32471369), the Guizhou Provincial Science and Technology Projects (Nos. ZK2024-167, ZK2022-397 and ZK2021-029, 2021-5637), the Excellent Young Talents Plan of Guizhou Medical University (No. 2021-101), and the Youth Guidance Project of Guizhou Provincial Science and Technology (No. 2024-245).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-484/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Diao X, Guo C, Jin Y, et al. Cancer situation in China: an analysis based on the global epidemiological data released in 2024. Cancer Commun (Lond) 2025;45:178-97. [Crossref] [PubMed]
2. Siegel RL, Kratzer TB, Giaquinto AN, et al. Cancer statistics, 2025. CA Cancer J Clin 2025;75:10-45. [Crossref] [PubMed]
3. Song L, Yu X, Wu Y, et al. Integrin β8 Facilitates Macrophage Infiltration and Polarization by Regulating CCL5 to Promote LUAD Progression. Adv Sci (Weinh) 2025;12:e2406865. [Crossref] [PubMed]
4. Liu Y, Lankadasari M, Rosiene J, et al. Modeling lung adenocarcinoma metastases using patient-derived organoids. Cell Rep Med 2024;5:101777. [Crossref] [PubMed]
5. Sung JY, Cheong JH. New Immunometabolic Strategy Based on Cell Type-Specific Metabolic Reprogramming in the Tumor Immune Microenvironment. Cells 2022;11:768. [Crossref] [PubMed]
6. Karta J, Bossicard Y, Kotzamanis K, et al. Mapping the Metabolic Networks of Tumor Cells and Cancer-Associated Fibroblasts. Cells 2021;10:304. [Crossref] [PubMed]
7. Stine ZE, Schug ZT, Salvino JM, et al. Targeting cancer metabolism in the era of precision oncology. Nat Rev Drug Discov 2022;21:141-62. [Crossref] [PubMed]
8. Huang K, Zhang Y, Gong H, et al. Inferring evolutionary trajectories from cross-sectional transcriptomic data to mirror lung adenocarcinoma progression. PLoS Comput Biol 2023;19:e1011122. [Crossref] [PubMed]
9. Kim N, Kim HK, Lee K, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun 2020;11:2285. [Crossref] [PubMed]
10. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
11. Cai L, Lin S, Girard L, et al. LCE: an open web portal to explore gene expression and clinical associations in lung cancer. Oncogene 2019;38:2551-64. [Crossref] [PubMed]
12. De Vitis C, Battaglia AM, Pallocca M, et al. ALDOC- and ENO2- driven glucose metabolism sustains 3D tumor spheroids growth regardless of nutrient environmental conditions: a multi-omics analysis. J Exp Clin Cancer Res 2023;42:69. [Crossref] [PubMed]
13. Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018;36:411-20. [Crossref] [PubMed]
14. Wu Y, Yang S, Ma J, et al. Spatiotemporal Immune Landscape of Colorectal Cancer Liver Metastasis at Single-Cell Level. Cancer Discov 2022;12:134-53. [Crossref] [PubMed]
15. Teng JS, Wang Y. An immune-related eleven-RNA signature-drived risk score model for prognosis of osteosarcoma metastasis. Sci Rep 2024;14:13401. [Crossref] [PubMed]
16. Huang A, Sun Z, Hong H, et al. Novel hypoxia- and lactate metabolism-related molecular subtyping and prognostic signature for colorectal cancer. J Transl Med 2024;22:587. [Crossref] [PubMed]
17. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
18. Zanotelli MR, Zhang J, Reinhart-King CA. Mechanoresponsive metabolism in cancer cell migration and metastasis. Cell Metab 2021;33:1307-21. [Crossref] [PubMed]
19. Espinosa-Rodriguez BA, Treviño-Almaguer D, Carranza-Rosales P, et al. Metformin May Alter the Metabolic Reprogramming in Cancer Cells by Disrupting the L-Arginine Metabolism: A Preliminary Computational Study. Int J Mol Sci 2023;24:5316. [Crossref] [PubMed]
20. Arner EN, Rathmell JC. Metabolic programming and immune suppression in the tumor microenvironment. Cancer Cell 2023;41:421-33. [Crossref] [PubMed]
21. Jovic D, Liang X, Zeng H, et al. Single-cell RNA sequencing technologies and applications: A brief overview. Clin Transl Med 2022;12:e694. [Crossref] [PubMed]
22. Aran D. Single-Cell RNA Sequencing for Studying Human Cancers. Annu Rev Biomed Data Sci 2023;6:1-22. [Crossref] [PubMed]
23. Zhao J, Shi Y, Cao G. The Application of Single-Cell RNA Sequencing in the Inflammatory Tumor Microenvironment. Biomolecules 2023;13:344. [Crossref] [PubMed]
24. Zhang C, Wang H, Chen Z, et al. Carbonic anhydrase 2 inhibits epithelial-mesenchymal transition and metastasis in hepatocellular carcinoma. Carcinogenesis 2018;39:562-70. [Crossref] [PubMed]
25. Li XJ, Xie HL, Lei SJ, et al. Reduction of CAII Expression in Gastric Cancer: Correlation with Invasion and Metastasis. Chin J Cancer Res 2012;24:196-200. [Crossref] [PubMed]
26. Pan R, Dai J, Liang W, et al. PDE4DIP contributes to colorectal cancer growth and chemoresistance through modulation of the NF1/RAS signaling axis. Cell Death Dis 2023;14:373. [Crossref] [PubMed]
27. Feng B, Pan B, Huang J, et al. PDE4D/cAMP/IL-23 axis determines the immunotherapy efficacy of lung adenocarcinoma via activating the IL-9 autocrine loop of cytotoxic T lymphocytes. Cancer Lett 2023;565:216224. [Crossref] [PubMed]
28. Powers GL, Hammer KD, Domenech M, et al. Phosphodiesterase 4D inhibitors limit prostate cancer growth potential. Mol Cancer Res 2015;13:149-60. [Crossref] [PubMed]
29. Jeong MH, Urquhart G, Lewis C, et al. Inhibition of phosphodiesterase 4D suppresses mTORC1 signaling and pancreatic cancer growth. JCI Insight 2023;8:e158098. [Crossref] [PubMed]
30. Andrews LP, Butler SC, Cui J, et al. LAG-3 and PD-1 synergize on CD8+ T cells to drive T cell exhaustion and hinder autocrine IFN-γ-dependent anti-tumor immunity. Cell 2024;187:4355-4372.e22. [Crossref] [PubMed]
31. Sorrentino V, Romani M, Mouchiroud L, et al. Enhancing mitochondrial proteostasis reduces amyloid-β proteotoxicity. Nature 2017;552:187-93. [Crossref] [PubMed]
32. Mantle D, Hargreaves IP. Mitochondrial Dysfunction and Neurodegenerative Disorders: Role of Nutritional Supplementation. Int J Mol Sci 2022;23:12603. [Crossref] [PubMed]
33. Weiss-Sadan T, Ge M, Hayashi M, et al. NRF2 activation induces NADH-reductive stress, providing a metabolic vulnerability in lung cancer. Cell Metab 2023;35:487-503.e7. [Crossref] [PubMed]
34. Amjad S, Nisar S, Bhat AA, et al. Role of NAD(+) in regulating cellular and metabolic signaling pathways. Mol Metab 2021;49:101195. [Crossref] [PubMed]
35. Wang H, Chen S, Kang W, et al. High dose isoleucine stabilizes nuclear PTEN to suppress the proliferation of lung cancer. Discov Oncol 2023;14:25. [Crossref] [PubMed]
36. Dai X, Zhang X, Chen W, et al. Dihydroartemisinin: A Potential Natural Anticancer Drug. Int J Biol Sci 2021;17:603-22. [Crossref] [PubMed]
37. Zhou Z, Lei J, Fang J, et al. Dihydroartemisinin remodels tumor micro-environment and improves cancer immunotherapy through inhibiting cyclin-dependent kinases. Int Immunopharmacol 2024;139:112637. [Crossref] [PubMed]
38. Ma L, Yang F, Wu X, et al. Structural basis and molecular mechanism of biased GPBAR signaling in regulating NSCLC cell growth via YAP activity. Proc Natl Acad Sci U S A 2022;119:e2117054119. [Crossref] [PubMed]
39. Gonzales E, Matarazzo L, Franchi-Abella S, et al. Cholic acid for primary bile acid synthesis defects: a life-saving therapy allowing a favorable outcome in adulthood. Orphanet J Rare Dis 2018;13:190. [Crossref] [PubMed]