Machine learning-based prognostic model of lactylation-related genes for predicting prognosis and immune infiltration in patients with lung adenocarcinoma

Abstract

Background

Histone lactylation is a novel epigenetic modification that is involved in a variety of critical biological regulations. However, the role of lactylation-related genes in lung adenocarcinoma has yet to be investigated.

Methods

RNA-seq data and clinical information of LUAD were downloaded from TCGA and GEO datasets. Unsupervised consistent cluster analysis was performed to identify differentially expressed genes (DEGs) between the two clusters, and risk prediction models were constructed by Cox regression analysis and LASSO analysis. Kaplan–Meier (KM) survival analysis, ROC curves and nomograms were used to validate the accuracy of the models. We also explored the differences in risk scores in terms of immune cell infiltration, immune cell function, TMB, TIDE, and anticancer drug sensitivity. In addition, single-cell clustering and trajectory analysis were performed to further understand the significance of lactylation-related genes. We further analyzed lactate content and glucose uptake in lung adenocarcinoma cells and tissues. Changes in LUAD cell function after knockdown of lactate dehydrogenase (LDHA) by CCK-8, colony formation and transwell assays. Finally, we analyzed the expression of KRT81 in LUAD tissues and cell lines using qRT-PCR, WB, and IHC. Changes in KRT81 function in LUAD cells were detected by CCK-8, colony formation, wound healing, transwell, and flow cytometry. A nude mouse xenograft model and a KrasLSL-G12D in situ lung adenocarcinoma mouse model were used to elucidate the role of KRT81 in LUAD.

Results

After identifying 26 lactylation-associated DEGs, we constructed 10 lactylation-associated lung adenocarcinoma prognostic models with prognostic value for LUAD patients. A high score indicates a poor prognosis. There were significant differences between the high-risk and low-risk groups in the phenotypes of immune cell infiltration rate, immune cell function, gene mutation frequency, and anticancer drug sensitivity. TMB and TIDE scores were higher in high-risk score patients than in low-risk score patients. MS4A1 was predominantly expressed in B-cell clusters and was identified to play a key role in B-cell differentiation. We further found that lactate content was abnormally elevated in lung adenocarcinoma cells and cancer tissues, and glucose uptake by lung adenocarcinoma cells was significantly increased. Down-regulation of LDHA inhibits tumor cell proliferation, migration and invasion. Finally, we verified that the model gene KRT81 is highly expressed in LUAD tissues and cell lines. Knockdown of KRT81 inhibited cell proliferation, migration, and invasion, leading to cell cycle arrest in the G0/G1 phase and increased apoptosis. KRT81 may play a tumorigenic role in LUAD through the EMT and PI3K/AKT pathways. In vivo, KRT81 knockdown inhibited tumor growth.

Conclusion

We successfully constructed a new prognostic model for lactylation-related genes. Lactate content and glucose uptake are significantly higher in lung adenocarcinoma cells and cancer tissues. In addition, KRT81 was validated at cellular and animal levels as a possible new target for the treatment of LUAD, and this study provides a new perspective for the individualized treatment of LUAD.

Background

Lung cancer is one of the most common cancers and the leading cause of cancer-related deaths worldwide [1]. Lung adenocarcinoma (LUAD) is the predominant histologic subtype, accounting for nearly 50% of all lung cancers [2]. Despite significant advances in the treatment of LUAD due to innovations in surgery, chemotherapy, targeted and immunotherapy, the prognosis of LUAD patients remains poor, with 5-year survival rates below 20% [3, 4]. The development of LUAD is a complex process that may be closely related to the aberrant expression of specific genes [5]. Therefore, searching for novel molecular markers and effective therapeutic targets for LUAD is essential to improve the poor prognosis and therapeutic outcome of LUAD patients.

Normal cells rely mainly on mitochondrial oxidative phosphorylation (OXPHOS) to generate energy for cellular processes. However, most cancer cells rely on aerobic glycolysis and secrete more lactate, a phenomenon known as the "Warburg effect" [6]. Lactate is involved in biological processes such as energy regulation, redox buffering, and regulation of fatty acid metabolism [7]. Lactate-derived lysine lactylation of histones is a novel epigenetic modification that directly stimulates gene transcription in chromatin [8]. Histone lactylation is involved in a variety of critical biological regulations. For example, Ricardo A. et al. found that the conversion of inflammatory macrophages to reparative macrophages is regulated by the B-cell junction of PI3K (BCAP) through histone lactylation [9]. Lactic acid-induced lactylation also recruits CAFs, macrophages, dendritic cells, and regulatory T cells in the tumor microenvironment, in order to remodel the tumor microenvironment and promote tumor development [10]. The results of xiong et al. highlight the importance of lactylation-driven METTL3-mediated RNA m6A modification in promoting the immunosuppressive capacity of tumor-infiltrating myeloid cells (TIMs) [11]. Metabolic reprogramming is a well-known feature of cancer [12]. Lactylation reflects the levels of lactate (an important metabolite), which in turn drives lactylation. This establishes an intrinsic link between lactylation and cellular metabolism. Histone lactylation is the link between reprogramming of cellular metabolism and transcriptome disruption in cancer cells [13].

In ocular melanoma, histone lactylation promotes tumorigenesis by facilitating YTHDF2 expression [14]. In hepatocellular carcinoma (HCC), lactylation affects various enzymes involved in metabolic pathways and promotes proliferation and metastasis of HCC cells [15]. Researchers have found that non-small cell lung cancer preferentially utilizes lactate over glucose to promote the TCA cycle and maintain tumor metabolism in vivo [16]. Upregulation of lactate dehydrogenase is associated with poor prognosis in lung adenocarcinoma, and lactate regulates cellular metabolism and promotes lung adenocarcinoma through histone lactylation-mediated gene expression [17]. Wang et al. found that BZW2 leads to the development of LUAD through glycolysis-mediated lactylation of the IDH3G through in vitro and in vivo experiments [18].

However, studies of lactylation-related genes in the LUAD tumor environment are still limited. In this paper, we analyzed the expression of lactylation-related genes in lung adenocarcinoma and developed a machine learning-based optimal prognostic model. LUAD patients were stratified using risk scores, and immune infiltration analysis, chemotherapeutic drug sensitivity, enrichment analysis, and clinical correlation analysis were performed for high- and low-risk populations using the R language. The results showed that LUAD patients with high lactylation scores had greater immune escape potential and lower immunotherapy response rates. Single-cell RNA sequencing analysis (scRNA-seq) verified the expression levels of model genes in LUAD. For trajectory analysis, MS4A1 decreased gradually as B cells differentiated into plasma cells. We further found that lactate content was abnormally increased in lung adenocarcinoma cells and cancer tissues, and glucose uptake capacity of lung adenocarcinoma cells was significantly enhanced. Knockdown of LDHA inhibited tumor cell proliferation, migration and invasion. In addition, we selected the model gene KRT81 for a series of in vivo and in vitro experiments and found that KRT81 was highly expressed in lung adenocarcinoma cell lines and tissues, and that knockdown of KRT81 significantly inhibited the proliferation, migration, and invasive ability of LUAD cells, leading to G0/G1 phase blockage and promoting their apoptosis. In this study, we elucidated the prognostic value of lactylation-related genes in LUAD patients through bioinformatics and experimental validation to provide new options for the individualized treatment of LUAD patients.

Methods

Data sources

We obtained mRNA matrix data (FPKM), mutation data, and RNA stemness scores (RNAss) for 541 LUAD samples and 59 normal tissue samples from TCGA (https://portal.gdc.cancer.gov/). From the TCGA database, we also obtained relevant clinical data of LUAD patients. The LUAD dataset GSE68465 was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), based on the GPL96 platform [HG-U133A] Affymetrix Human Genome U133A Array, which contains 442 lung adenocarcinoma samples gene matrix files and clinical information; The integrated ids were then converted to official gene symbols and the data were log2 processed. LUAD samples without survival time were excluded from this study. Expression profiles of LUAD samples in the TCGA database were converted from FPKM to TPM, and batch effects were removed by normalizing the TCGA and GEO datasets with the ComBat function of the SVA package in R. Based on previously published studies, we included 327 lactylation-related genes (Table S1) [8, 19, 20]. Ethics committee approval was not required as the clinical information of the patients participating in this study was obtained from the TCGA and GEO databases, and the publication guidelines were strictly adhered to.

The mutation frequency of lactylation differential genes (DE-LACAGs) was explored using the R package "maftools". The amplified and deleted copy number variations (CNVs) of DE-LACAGs were investigated based on the TCGA database. Protein–protein interaction (PPI) networks were constructed using the STRING tool (v11. 5, https://www.string-db.org/). The location of DE-LACAGs in chromosomes was explored using the R package "RCircos".

Consensus clustering based on lactylation-related genes

Unsupervised consensus clustering analysis was performed using "Consensus ClusterPlus" from the R package to categorize patients into multiple clusters. To ensure the stability of the supervised clustering, all processes were repeated 1000 times. Survival analyses were then compared across clusters.

Prognostic risk modeling

Differentially expressed genes (DEGs) were identified in 2 Clusters with |log2FC |≥ 1 and p-value < 0.05 using the "limma" package in R. Survival-associated DEGs were identified by univariate Cox regression analysis. To avoid overfitting, the least absolute shrinkage and selection operator (LASSO) was used to select genes with high prognostic values. Next, 1000 LASSO iterations were performed using the 'glmnet' package in R to construct the prognostic models and obtain their regression coefficients. LASSO regression is commonly used to fit selection variables in high-dimensional generalized linear models [21]. it is achieved by constructing a penalty function to obtain a finer model, effectively avoiding overfitting. Significant contributing variables are then entered into a stepwise descending logistic regression analysis. A stepwise approach was used to avoid covariance due to the omission of redundant variables. The final combined model was generated from the independent significant genetic variables and validated by a tenfold cross-validation scheme. When λ was smallest, it retained 21 genes, and finally, further screening using multivariate Cox regression identified the prognosis for 10 genes to be modeled. Ten genes were finally identified to participate in the prognostic risk modeling, including IGFBP1, CYP17A1, DKK1, KRT81, MS4A1, C11orf16, BCAN, FBN2, ANGPTL4, and SERPINB7. Risk scores were calculated as follows:

$${\text{Risk score}} = \Sigma {\text{expGenei}}*\beta {\text{i}}$$

where expGenei is the relative expression of the prognostic model gene, and βi is the regression coefficient.

Risk assessment models

Subjects were categorized into high-risk and low-risk groups based on the median risk score. Overall survival (OS) was analyzed using the Kaplan–Meier method. Patient survival status based on risk score was plotted using the pheatmap software package. Receiver Operating Characteristic (ROC) analysis and area under the curve (AUC) calculation were performed using the "timeROC" R software package. Multivariate Cox regression analysis was used to investigate whether the risk model was an independent risk factor excluding other clinical characteristics (age, sex, and stage) in the training and validation sets. Clinicopathologic factors were combined with our constructed risk scores to construct column-line plots to predict OS rates at 1, 3, and 5 years in patients with LUAD. We constructed line plots to predict the OS rates of LUAD patients at 1, 3, and 5 years, taking into account clinicopathologic factors. The R software packages “rms”, “regplot” and “survival” were used to construct the column plots and their corresponding calibration curves, and the closer the calibration curves were to the diagonal, the more accurate the prediction was. The closer the calibration curve is to the diagonal, the more accurate the prediction is.

Tumor Mutational Burden (TMB) analysis and immune-related functional analysis

We downloaded TMB-related data from TCGA samples and analyzed the number of mutations in both subgroups of LUAD patients using the R package (Maftools package). Use Expression Data (ESTIMATE) to estimate stromal and immune cells in tumor tissue to calculate tumor stroma score, immune score [22]. The scores of 22 immune cell subtypes in each tumor sample were identified by CIBERSORT (identification of cell types by estimation of relative isoforms of RNA transcripts) [23].

We used the “limma” R package to normalize the transcriptomic data of OC patients and uploaded the prepared data to the CIBERSORT algorithm (R script v1.03) to assess the distribution of different immune cells in each sample. The relative abundance of different immune cell types in the high- and low-risk groups was quantified and assessed to compare and predict immune cell infiltration between the two groups. Immune function correlation analysis was based on ssGSEA [24].

Predicting response to immunotherapy and chemotherapy

NSCLC immune dysfunction and rejection were obtained by TIDE prediction of response to immunotherapy and chemotherapy (http://tide.dfci.harvard.edu/), and scores were analyzed in the high- and low-risk groups using the "limma" and "ggpubr" R packages. TIDE scores accurately predicted the efficacy of immunotherapeutic agents received by patients [25], and higher TIDE scores predicted poorer response to immunotherapy. Correlations between risk scores and TMB were analyzed using the Spearman correlation. The Genomics of Drug Sensitivity in Cancer (GDSC) is a public dataset containing information on molecular markers of drug sensitivity and drug response in cancer cells [26]. The "oncoPredict" package is used to predict the sensitivity of high- and low-risk groups to various antitumor drugs [27].

Biological Process and Pathway Enrichment Analysis

KEGG analysis was performed using the R clusterProfiler software package [28, 29]. Gene set enrichment analysis (GSEA) was performed using GSEA software (version 4.1.0).

Single-cell sequencing analysis

The single-cell dataset was downloaded from https://codeocean.com/capsule/8321305/tree/v1. It contains 10 tumor samples and 10 paired normal samples. The scRNA-seq data were analyzed using the R software "Seurat" package [30]; Low-quality cells were excluded by calculating the percentage of mitochondrial or ribosomal genes. Data were normalized using the "NormalizeData" function. The "FindVariableFeatures" function was used to filter the top 2000 highly variable genes. The "RunPCA" function was used to perform principal component analysis (PCA) on 2000 genes. The "Harmony" function performs batch correction for different samples and uses uniform manifold approximation and projection (UMAP) [31] for downscaling and cluster identification. The "SingleR" package [32] is applied to annotate different clusters of cell subgroups. The R package "AddMouduleScore" was used to calculate the prognostic model score. The "monocle" package [33] was used for cell trajectory and pseudotemporal analysis, and the "DDRTree" function was used for dimensionality reduction. B-cell subpopulations were annotated according to Hao et al. [34].

Tissue specimen collection and lung adenocarcinoma cell culture

All tissue samples were obtained from the Department of Thoracic Surgery, Northern Jiangsu People’s Hospital and approved by the Medical Ethics Committee of the hospital. We obtained informed consent from each relevant patient before collection. Twelve pairs of samples were obtained from patients with lung adenocarcinoma who underwent tumor resection from January 2020 to December 2021, including tumor tissues (T) and paired normal tissues (N), and the pathological type of all LUAD cases was lung adenocarcinoma. All samples were stored at -80 °C. HBE, A549, H1975, H1299, and PC9 cell lines were obtained from the China Cell Resource Center (Shanghai, China). Cells were cultured in RPMI 1640 (Solarbio, 31800) medium supplemented with 10% fetal bovine serum (FBS) (Procell, 164210–50). Cells were incubated in a humidified incubator (Thermo Scientific, China) at 5% CO2, 37 °C.

RNA extraction and quantitative real-time polymerase chain reaction qRT-PCR

RNA was extracted from cells using TRIzol reagent (Vazyme). We measured RNA concentration using a spectrophotometer and stored the samples at -80 °C. cDNA was synthesized using HiScript Ill RT SuperMix for gPCR (+ gDNA wiper) (Vazyme, China). Hieff®qPCR SYBR Green Master Mix (High Rox Plus) (Yeasen Biotechnology, Shanghai, China) was used to perform real-time fluorescence quantitative PCR in the StepOne Plus real-time PCR System (Applied Biosystems). Primers Sequencing for qPCR.GAPDH:F:5'-AATGGGCAGCCGTTAGGAAA-3',

R:5'-GCGCCCAATACGACCAAATC-3'.LDHA: F:5'-GCCGTGATAATGACCAGCTT-3', R: 5'-TGGCAGCCTTTTCCTTAGAA-3'.

KRT81:F:5'-GCATTGGGGCTGTGAATGTCT-3',R:5'-ACCCAGGGAGCTGATACCAC-3'. The siRNA and siRNA negative control (si-NC) were purchased from GenePharma (Shanghai, China). siRNA sequences were as follows: LDHA: siNC:UUCUCCGAACGUGUCACGUTT-3';siRNA1:5'GCUGAUUUAUAAUCUUCUATT-3';siRNA2:5′-GAAUAAGAUUACAGUUGUUTT-3';siRNA3:5′-GCCCGAUUCCGUUACCUAATT-3'.KRT81:siNC:5'-GGCTCTAGAAAAGCCTATGC-3';siRNA1:5'-GCATCATTGCCGAGATTAA-3';siRNA2:5′-CACAGTATGACGACATTGT-3′;siRNA3:5'-GGCCAATTGAACACCACCT -3'. Cells were incubated in 6-well plates, and transfection was started when cell density reached 60%. Transfection was performed using gp-transfection-mate (GenePharma). Transfection efficiency was detected using qRT-PCR.

Western blot (WB) assay

Whole cell or tissue mixtures were separated using RIPA lysate (Solarbio, item number:R0020), PMSF, protease inhibitor, and protein phosphatase inhibitor mixtures, and equal amounts of proteins were separated on 10% SDS-PAGE and transferred to a PVDF membrane (Immobilon-P, IPVH00010). The membranes were closed with 5% skimmed milk and incubated with primary antibodies KRT81 (Proteintech,11342–1-AP), GAPDH (Proteintech, 10494–1-AP), MMP2 (Proteintech,10373–2-AP), MMP9 (Proteintech, 10375–2-AP), E-cadherin (Proteintech,20874–1-AP), N-cadherin (Proteintech, 22018–1-AP), Vimentin (Proteintech,10366–1-AP), p-PI3K (Cell Signaling, 17366), t-PI3K (Cell Signaling, 9655), p-AKT (Cell Signaling,13038S), t-AKT (Cell Signaling, 4685S), P53 (Proteintech, 10442–1-AP), Caspase3 (Cell Signaling, 9664 T) and Bcl2 (Cell Signaling, 4223 T) were combined overnight at 4 °C. They were then incubated in secondary IgG (ABclonal, AS014) for 1 h at room temperature. Protein bands were visualized using Super ECL detection reagent (Yeasen Biotechnology, Shanghai, China). Grayscale analysis of protein bands was performed using Image J software.

lactate assay in cells and tissues

Lactate accumulation in the cell culture medium was assessed using a commercially available lactate assay kit (Solarbio® BC2235) according to the manufacturer's instructions.

Glucose uptake measurement

Glucose uptake in the cell culture medium was assessed using a commercially available glucose uptake fluorometric test kit (Elabscience E-BC-F041) according to the manufacturer's instructions.

Immunohistochemical staining and scoring

Tissues cut into 4 μm thick paraffin-embedded tissues were dewaxed and hydrated and then antigenically repaired in pH 6.0 sodium citrate antigen repair solution in a microwave oven for 2 min on high and 15 min on bottom. The sections were incubated in endogenous peroxidase containment solution for 15 min and immunostaining containment solution for 15 min. Sections were incubated overnight at 4 °C with anti-KRT81 (Proteintech, 11342–1-AP 1:100), anti-Ki67 (Proteintech, 28074–1-AP 1:1000), and anti-PCNA (Cell Signaling, 13110S 1:500). After being left at room temperature, the nuclei were conjugated with a secondary IgG antibody (Servicebio, G1215-200 T) for 1 h, and the nuclei were stained with hematoxylin using DAB as a chromogen. HE staining was performed using Hematoxylin Eosin Staining Kit (Solarbio, G1120): Hematoxylin staining solution for 10 min, distilled water rinse. The differentiation solution was differentiated for 30 s and rinsed twice with distilled water. Eosin stain was placed for 1 min and then dehydrated, transparent, and sealed.

The positive rate of KRT81 in tissues was assessed as follows: staining intensity scores 0 (negative), 1 (weak), 2 (moderate), 3 (strong), and positive cell ratio scores 0 (0–5%), 1 (6–25%), 2 (26–50%), 3 (51–71%), and 4 (more than 75%), and the scoring reference plots are shown in the Supplementary Material (Fig. S1), and the final scores were designated as the product of staining intensity and positive rate cell score.

Cell proliferation assay

A cell proliferation assay was performed in 96-well plates, 1000 cells were added to each well after counting. 24, 48, 72, 96 h later, 10 μL CCK-8 solution (Yeasen) was added to each well, incubated for 1 h, and the absorbance at 450 nm (OD) was detected by an enzyme labeling instrument (Skanlt RE 7.0).

Colony formation assay

Cells were transfected with siNC, siRNA1 and siRNA2. 6-well culture plates were inoculated with 1000 cells per well and incubated for 14 days. Cells were fixed in 4% paraformaldehyde for 20 min, stained with 0.1% crystal violet solution for 8 min, air-dried, and photographed, and colonies were counted using Image J software.

Cell migration and invasion assays

We used 8 µm pore size transwell chambers (Corning, USA) in 24-well plates; 200ul of cell suspension containing 10,000 cells without FBS with matrix gel (BD Biocoat) or without matrix gel (Corning) was added to the upper chamber, and 500ul of cell suspension containing 10% FBS was added to the lower chamber. Cells were cultured in a cell incubator for 48 h. Cells suspended in the chamber were rinsed with phosphate buffered saline (PBS), cells were fixed with 4% paraformaldehyde for 20 min, stained with 0.1% crystal violet solution for 6 min, washed three times with PBS, and the cells on the upper surface of the bottom chamber were gently wiped with a cotton swab. Cells were located on the lower surface of the bottom of the drying chamber, and images were taken with an inverted microscope (OLYMPUS-CKX53). Cell counting using Image J.

Wound Healing Assay

Wound healing assay was used to assess the migration ability of the cells. Transfected cells (10 × 10⁴/well) were inoculated in 6-well plates with a monolayer of cells evenly distributed on the bottom of the plate. The cell layer was scraped with a 200 μl pipette tip, washed with PBS, and FBS-free medium was added to each well. Images were then taken under an inverted microscope at 0, 24, and 48 h (OLYMPUS-CKX53, China). The images were analyzed using Image J software。

Flow cytometry

Cell cycle: PBS washed 3 times, trypsin digested for 5 min, fixed in 70% ethanol at 4 °C for 30 min, then incubated in 500 μl propidium iodide staining solution (PI) (Beyotime) at 37 °C for 30 min. Cell apoptosis: PBS was washed 3 times, trypsin digested for 5 min, then annexin-FITC reagent (Beyotime) was added sequentially and incubated for 20 min away from light. Cell cycle and apoptosis were detected by flow cytometry (BD Biosciences, USA), and the results were analyzed using FlowJo software.

Xenograft model

A mouse xenograft model was established to explore the functional role of KRT81 in vivo. A549 cells were placed in six-well plates and transfected with siNC and siRNA2. 8 BALB/c nude mice were selected from Nanjing Gempharmatech Co., Ltd, and randomly divided into two groups, and the above-treated cells (9 × 10⁵) were injected into the axilla of nude mice. Tumor volume was monitored every 7 days and calculated as:V = (Length × Width²) × 0.5. Nude mice were executed after 7 weeks, tumor size was recorded after isoflurane anesthesia, and nude mice were dislocated from the cervical vertebrae until death.

KRAS lung cancer in situ mouse model

4–6-week-old B6-KrasLSL-G12D/ + mice were obtained from Nanjing Gempharmatech Co., Ltd and bred under specific pathogen-free (SPF) conditions. After anesthesia with isoflurane, mice were injected with 5 × 10⁷ PFU adenovirus (Ad-Cre, OBIO) by intratracheal drip and were executed 20 weeks after infection. Animal experiments were performed by animal care guidelines and approved by the Laboratory Animal Ethics Committee of Yangzhou University.

Statistical methods

Statistical analysis was performed using GraphPad Prism (8.4.3) and R software (4.2.2). The normal distribution variables were analyzed by t test. Non-normally distributed variables were analyzed using Wilcoxon rank sum test. P < 0.05 was considered statistically significant.

Results

Screening and characterization of genes associated with lactonization in lung adenocarcinoma

To identify genes associated with differential lactylation in lung adenocarcinoma (LUAD), we obtained 541 LUAD samples and 59 normal tissue samples from TCGA. We performed differential analysis using the limma package (threshold for differential genes |log2FC|≥ 1, p-value < 0.05). This analysis revealed 5145 genes (2593 up-regulated and 3027 down-regulated) differentially expressed in LUAD and normal tissues (Table.S2). We then intersected these results with 327 lactylation-associated genes, yielding 52 genes (Fig. 1A). Further univariate Cox analysis screened 26 lactylation prognosis-associated genes, and box plots showed the expression of lactylation-associated genes in lung adenocarcinoma and normal tissues (Fig. 1B). Our PPI analysis using STRING revealed four genes H2AX, CCNA2, GAPDH and ENO1 showed extensive association in the PPI network (Fig. 1E). In addition, somatic mutation characterization showed that 120 (30.77%) of the 390 samples in the TCGA database were significant for DE-LACAGs mutation frequency, with 12%, 7% and 4% for AHNAK, PRKDC and MKI67, respectively (Fig. 1C). We further analyzed the CNV of DE-LACAGs in LUAD, and the results showed that HDGF, NSUN2, CCT5, MNDA, and PABPC1 had high amplification rates; In contrast, CCNA2, LCP1, and GAPDH had high CNV deletions (Fig. 1D). The locations of DE-LACAGs on chromosomes were shown in Fig. 1F.

Fig. 1

Identification of lactylation prognosis-associated genes in lung adenocarcinoma. A 52 overlapping genes were identified as lactylation-associated DEGs. B Expression of 26 lactylation prognosis-associated genes in lung adenocarcinoma in Tumor and Normal. C 120 out of 390 patients with LUAD showed genetic alterations in prognosis-associated lactylation genes. D CNV mutations are widespread in 26 prognosis-associated lactylation genes. The column indicates the frequency of changes. Loss, green dots; GAIN, pink dots. E PPI network showing interactions between proteins encoded by lung adenocarcinoma lactylation prognosis-associated genes. F Locations of CNV alterations in lactylation prognosis-related genes on chromosome 23

Identification of subgroups of LACAGs in LUAD

We utilized unsupervised consensus clustering to construct a patient-specific classification of LUAD based on 26 LACAGs. The heatmap used the best classification with k = 2 to cluster the 507 LUAD samples into cluster C1 (n = 194, 38.3%) and cluster C2 (n = 313, 61.7%) (Fig. 2A). Kaplan–Meier survival analysis showed that LUAD patients in cluster C1 had a significantly worse prognosis compared to cluster C2 (Fig. 2B). The principal component analysis plot showed that LUAD samples in LACAGs clusters C1 and C2 could be clearly distinguished (Fig. 2C). In addition, clinicopathologic features and 26 LACAGs expression were shown in the heat map, and most DEGs were in gene cluster C1 (Fig. 2D).

Fig. 2

Identification of subgroups of LACAGs in LUAD. A Unsupervised consensus clustering analysis of LUAD. B Clinical survival analysis of C1 and C2 clusters. C Graph of PCA analysis of C1 and C2 clusters. D Characterization of clinicopathology and expression of LACAGs in LUAD subtypes

A total of 886 DEGs (threshold for differential genes |log2FC |≥ 1, p-value < 0.05) were identified between cluster C1 and cluster C2 (Table. S3). To reveal the biological pathways associated with LACAGs, we analyzed GO functional enrichment and KEGG pathway enrichment based on DEGs between cluster C1 and cluster C2. GO enrichment analysis showed that DEGs were associated with nuclear division, chromosome segregation, microtubule, chromosomal region, and tubulin binding (Fig. S2A). KEGG found that DEGs were enriched in the Cell cycle, ECM-receptor interaction, Motor proteins, and Oocyte meiosis (Fig. S2B). GSEA enrichment analysis found that Cell cycle, ECM-receptor interaction, Focal adhesion, regulation of actin cytoskeleton, and spliceosome were enriched in the high-risk group (Fig. S2C).

Risk modeling of DEGs through machine learning

To explore the prognostic value of DEGs associated with the LACAGs subtype for LUAD, we built a risk model to explore the impact on their prognosis. 359 of the DEGs were identified as prognostic genes by univariate Cox regression analysis. Then, these prognostic DEGs were examined using LASSO-Cox regression analysis (Fig. 3A), which screened out 21 genes (Fig. 3B). Finally, multivariate Cox regression was performed to model the prognosis of 10 genes. Score = (0.12822022848453 × IGFBP1)—(0.644037444131072 × CYP17A1) + (0.096333204249806 × DKK1) + (0.056508018166641 × KRT81)—(0.123234435555447 × MS4A1)—(0.12557360722698 × C11orf16) + (0.308812791755336 × BCAN) + ( 0.222019823099794 × FBN2) + (0.0936025202708982 × ANGPTL4) + ( 0.0944548850230676 × SERPINB7). We also constructed Sankey diagrams of LACAGs clusters, risk scores, and patient survival status (Fig. 3C). Among the LACAGs clusters, the risk scores of LUAD samples were higher in the C1 cluster, which has a poorer clinical prognosis (Fig. 3D). We then used TCGA data as the training set and GSE68465 as the validation set. Based on the risk score and survival status showed that mortality increased with increasing score (Fig. 3E). Based on the risk grouping, we then performed a Kaplan–Meier (KM) survival analysis (Fig. 3F), and both the training and validation sets illustrated that the high-risk group had a worse prognosis. In addition, we used ROC curves to assess the predictive accuracy of risk scores. The AUC values for the 1-, 3-, and 5-year risk scores were 0.744, 0.734, and 0.715 for the training set, and the AUC values for the 1-, 3-, and 5-year risk scores were 0.612, 0.598, and 0.566 for the test set, respectively (Fig. 3E, F).

Fig. 3

Construction of risk modeling. A, B LASSO regression analysis of selected prognostic genes. C Alluvial diagram showing the relationship between survival status, DEGs clusters, and risk scores. D Difference-in-difference analysis of cluster risk scores. E ROC curves of training set risk score plots, survival status plots, and 1-, 3-, and 5-year risk scores. F ROC curves for the validation set risk score plot, survival status plot, and 1-, 3-, and 5-year risk scores

Development and validation of column-line diagrams

The results of multivariate COX regression analysis of the training and validation sets showed that the risk score could be used as an independent prognostic indicator for the prognosis of LUAD patients (both P < 0.001) (Fig. 4A, B). A nomogram of LUAD patients was constructed based on gender, age, clinical stage, and risk score (Fig. 4C). The calibration curves showed that the predicted OS of the 1-, 3-, and 5-year. Nomogram was generally consistent with the corresponding observed OS of LUAD patients (Fig. 4D). ROC analysis was performed based on the clinical characteristics of LUAD in the TCGA database, and risk scores were found to have greater prognostic power than other clinical characteristics (Fig. 4E).

Fig. 4

Independent prediction of risk models and nomogram construction. A,B Multivariate COX regression analysis of training and validation sets. C Nomogram survival prediction of LUAD patients with risk scores. D Calibration curves assessing 1-, 3- and 5-year OS consistency. E ROC curves for common clinical parameters and risk scores

Risk model-based estimation of tumor immune microenvironment and immune-related genes

In order to determine the difference in infiltrating immune cells between the high-risk and low-risk groups, we found by ESTIMATE algorithm that the low-risk group had higher ESTIMATE scores and immune scores compared to the high-risk group (Fig. 5A). Correlation of 10 genes involved in constructing the risk model with immune cells (Fig. 5B). We further explored the relationship between risk scores and immune-related features. ssGSEA analysis showed that higher risk scores were significantly associated with reduced levels of most immune-related features (Fig. 5C), including immune cell infiltration (e.g., B cells, HLA, iDCs, Mast_cells, T_cell_co- stimulation, T_helper_cells, Tfh, TIL and Type_II_IFN_Reponse). Risk scores were correlated with Neutrophils (R = 0.18,p = 0.00014), Macrophages M0 (R = 0.24, p = 1.1e-07), NK cells (R = 0.15, p = 0.0013), and T cells CD4 memory activated (R = 0.15, p = 0.0018) were positively correlated, and negatively correlated with B cells memory (R = 0.26, p = 1.1e-0.8), Mast cells (R = -0.23, p = 6.9e-07) (Fig. 5D).

Fig. 5

Analysis of tumor immune microenvironment and immune infiltration. A StromalScore, ImmuneScore, and EstimatedScore between high and low-risk groups. B Correlation analysis of model genes with immune cells. C Differences in immune-related functions between high and low-risk groups. D Correlation between risk score and immune cells

Tumor mutation load analysis and drug treatment response prediction

In the expression of 26 LACAGs in high and low-risk groups, we found that except for PRAM, MNDA, LCP1, the rest were highly expressed in the high-risk group (Fig. 6A). We used the maftools algorithm to observe the mutations in the high-risk and low-risk groups, and the range of somatic mutations in the high-risk group was more extensive than that in the low-risk group (TP53: high-risk, 48%; low-risk, 39%, TTN: high-risk, 51%; low-risk, 36%, MUC16: high-risk, 41%; low-risk, 38%; and CSMD3: high-risk, 44%, low-risk, 32%) (Fig. 6B). The difference in TMB between the high-risk and low-risk groups was statistically significant (p < 0.05) (Fig. 6C). We also found that LACAGs clusters C1 and C2 increased TMB with increasing risk score (R = 0.18, p = 0.00096) (Fig. 6D).

Fig. 6

Tumor mutation load analysis. A Expression differences of 26 LACAGs between high and low-risk groups. B Waterfall plots showing the top 15 mutated genes of LUAD in the high-risk group (159 samples) and low-risk group (165 samples). C, D Differences in TMB between high- and low-risk groups and correlation analysis between risk scores and TMB

The tumor immune dysfunction and rejection (TIDE) algorithm calculated the risk of tumor immune escape. The results showed that the TIDE score was higher in the high-risk group than in the low-risk group, indicating a higher probability of immune escape (P < 0.001) (Fig. S3A). In addition, tumor stemness analysis showed that RNAss was increased in the high-risk group, and the risk score was positively correlated with RNAss (Fig. S3B). The RNAss values of LUAD samples are shown in (Table S4). We further explored the potential effective drugs for the treatment of LUAD patients using the "oncoPredict" R package. Among the commonly used NSCLC therapeutic agents, 5-Fluorouracil, Savolitinib sensitivity, Alisertib sensitivity, and Crizotinib sensitivity had lower IC50 values (50% inhibition of cell growth) in the high-risk group, suggesting that these agents may be more effective in high-risk patients may be more effective. However, low-risk patients may benefit more from Axitinib sensitivity, and Ribociclib sensitivity (Fig. S3C).

Single-cell RNA-seq analysis and trajectory analysis

Single-cell data were downscaled and clustered using Seurat, the harmony function was used for batch correction of data from different samples (Fig. S4A), UMAP downscaled and clustered tumor tissues into 19 clusters and normal tissues into 20 clusters (Fig. S4B), and cell types were annotated using SingleR (Fig. 7A). The expression of 10 lactylation-associated risk genes was calculated, and DKK1 was mainly expressed in Epithelial_cells. KRT81 was mainly expressed in T_cells. MS4A1 was mainly expressed in B_cells (Fig. 7B). AddMouduleScore in Seurat was applied to estimate the score of 10 relevant risk genes in each cell type. The results showed that the scores of B cells in tumor tissues were higher than those in normal tissues. This suggests our prognostic model may influence tumor progression by regulating B cells (Fig. 7C). B cells from tumor tissues were extracted and clustered into 12 subgroups (Fig. S4C). Pseudo-time series analysis using Monocle 2 identified 5 states of fate trajectories in their expression patterns (Fig. 7D). Interestingly, the relative expression of MS4A1 gradually decreased after B cells differentiated into plasma cells (Fig. 7E).

Fig. 7

Single-cell RNA-seq analysis in tumor and normal tissues. A Umap was used to downscale data from tumor tissues and normal tissues and annotate each cluster. B Expression of 10 lactylation-related genes in tumor tissues in different cells. C Scoring of lactylation-related genes in tumor tissues and normal tissues. D The 3 cell subtypes of the B cell differentiation process are shown. E Expression of MSA1 in the differentiation stages of the 3 B cell subtypes

Knockdown of LDHA inhibits lung adenocarcinoma cell proliferation, migration and invasion

To investigate the role of lactate in lung adenocarcinoma, we first examined the lactate content of lung adenocarcinoma cell lines and normal human bronchial epithelioid cell line (HBE) medium, and found that the lactate content was abnormally increased in lung adenocarcinoma cells (Fig. 8A), and further examined the lactate content of six pairs of cancerous tissues and paired normal tissues, which demonstrated that lung adenocarcinoma cancerous tissues had a significantly higher lactate content than that of the surrounding normal tissues (Fig. 8B). The glucose uptake capacity of lung adenocarcinoma cell lines A549 and H1975 was significantly higher than that of normal human bronchial epithelioid cell line (HBE) (Fig. 8C). Because of the higher lactate content of A549 and H1975 cells, siRNA for LDHA was transfected into A549 and H1975 cells and its knockdown efficiency was verified by qPCR (Fig. 8D). Selected siRNA1 was subjected to a series of functional experiments. CCK-8 and clone formation experiments showed that LDHA knockdown resulted in growth retardation of A549 and H1975 cells (Fig. 8E, F). transwell experiments showed that LDHA knockdown resulted in a significant reduction of cell migration and invasion ability (Fig. 8G).

Fig. 8

A Detection of lactate expression in HBE, A549, H1299, PC9 and H1975. B Determination of lactate concentration in human lung adenocarcinoma tissues and paired adjacent normal tissues (n = 6). C Glucose uptake capacity was measured in HBE, A549 and H1975. D The knockdown efficiency of LDHA in A549 and H1975 cells transfected with siRNA1 / 2 / 3 was detected using qRT-PCR. nc: no siRNA infection; siNC: negative control. E,F Proliferative capacity was detected by CCK8 and clone formation assay. G Transwell assay was used to detect the migration and invasion number of A549 cells and H1975 cells transfected with siNC and siRNA1

KRT81 expression was significantly elevated in LUAD cells and tissues

Based on GEPIA2 (http://gepia2.cancer-pku.cn), box plots showed that KRT81 expression was significantly up-regulated in 483 lung adenocarcinoma tissues compared with 347 normal tissues (Fig. 9A), and the KM survival curves showed that the prognosis of the KRT81 high-expression group was significantly worse than that of the KRT81 low-expression group (Fig. 9B). The differential expression of the model gene KRT81 in LUAD cell line and HBE cell line was verified by qRT-PCR and WB. It was found that the expression of LUAD cell line was higher than that of HBE (Fig. 9C, D). Further analysis of the protein levels of KRT81 in 12 pairs of lung adenocarcinoma tissues and paired normal tissues revealed that the expression of KRT81 was significantly higher in LUAD tissues (Fig. 9E). IHC results showed that the expression of KRT81 in lung adenocarcinoma tissues was significantly elevated and localized in the cytoplasm (Fig. 9F).

Fig. 9

Expression of KRT81 in LUAD cells and tissues. A We retrieved the expression profiles of KRT81 in 483 lung adenocarcinoma tissues versus 347 normal tissues on the GEPIA2 website. B KM curves of high and low KRT81 expression groups. C, D qRT-PCR and WB were performed to detect KRT81 expression in lung normal epithelial cell lines HBE and LUAD cell lines (A549, H1299, PC9, H1975). E WB detection of KRT81 expression in 12 pairs of LUAD tissues. F IHC detection of KRT81 expression in LUAD tissues. **p < 0.01, ***p < 0.001

Knockdown of KRT81 inhibits tumor cell proliferation, migration, invasion, and G1 phase inhibition and promotes apoptosis

A549 and H1975 cells highly expressed KRT81, so A549 and H1975 cells were transfected with siRNA, and the knockdown efficiency was verified by qPCR (Fig. 10A). The siRNA2 was selected for a series of functional experiments. CCK-8 and clone formation experiments showed that knockdown of the KRT81 gene resulted in growth retardation of A549 and H1975 cells (Fig. 10B, C). Wound healing and transwell assays showed that knockdown of KRT81 resulted in a significant reduction in cell migration and invasive ability (Fig. 10D,E). In addition, we found that knockdown of the KRT81 gene resulted in a significant decrease in the expression of the invasive proteins metalloproteinase 2 (MMP-2) and metalloproteinase 9 (MMP-9) proteins (Fig. 10E).

Fig. 10

KRT81 knockdown inhibits lung adenocarcinoma cell proliferation, migration, and invasion. A The knockdown efficiency of KRT81 in A549 cells and H1975 cells transfected with siRNA1 / 2 / 3 was determined by qRT-PCR. NC: no siRNA infection; siNC: negative control. B, C CCK8 and clone formation assay was used to detect the proliferation ability. D Transwell assay was used to detect the number of migrating and invading A549 cells and H1975 cells transfected with siNC and siRNA2. E Wound healing assay was used to detect the migration rate. WB assay was used to detect the expression levels of migration-associated proteins (MMP2, MMP9) in A549 cells transfected with siNC and siRNA2. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

We applied flow cytometry to explore whether KRT81 knockdown resulted in LUAD cell cycle arrest and increased apoptosis. In A549 cells, the proportion of cells in the G0/G1 phase was significantly increased, and the proportion of cells in the S phase was significantly decreased in the siRNA2 group compared with the siNC group (Fig. 11A). We further found that the apoptosis rate was significantly increased after knockdown of KRT81 (Fig. 11B). In addition, we analyzed the expression of apoptosis-related proteins. We found that Caspase3 and P53 were significantly up-regulated after KRT81 knockdown. However, Bcl-2 was significantly down-regulated (Fig. 11C). It was demonstrated that knockdown of KRT81 promoted apoptosis.

Fig. 11

Knockdown of KRT81 leads to G0/G1 phase block transition and apoptosis. A The effect of KRT81 knockdown on the cell cycle of A549 cells was detected by Cell Cycle and Apoptosis Analysis Kit. B The effect of KRT81 knockdown on apoptosis in A549 cells was detected by Annexin V Apoptosis Detection Kit. C WB detection of the expression levels of apoptosis-related proteins (Caspase3, P53, Bcl-2) in A549 cells. *p < 0.05, ***p < 0.001, ****p < 0.0001

KRT81 exerts its biological role through the EMT and PI3K/AKT pathway activation

Based on the effect of KRT81 on cell invasion, we further investigated the effect on Epithelial-Mesenchymal Transition (EMT) proteins after knocking down KRT81 and found that E-cadherin was significantly elevated and N-cadherin and Vimentin proteins were significantly decreased after knocking down KRT81 (Fig. 12A), which demonstrated that KRT81 regulates EMT to promote metastasis. Since changes in KRT81 expression had a significant impact on LUAD cell phenotype, we next explored the mechanism of KRT81's role in LUAD. To analyze the changes in gene expression and pathways that KRT81 may affect, differential gene analysis was first performed. We observed that the expression level of KRT81 could significantly alter hundreds of genes (Table S5). KEGG enrichment analysis was performed on the up-regulated genes. The KEGG results showed that the PI3K-AKT signaling pathway, Human T-cell leukemia virus 1 infection, and Regulation of actin cytoskeleton pathway were significantly enriched (Fig. 12B). In addition, the mRNA expression level of KRT81 in LUAD was divided into high and low expression groups at the median value for GSEA enrichment, and it was found that the upregulation of KRT81 was able to activate the PI3K/AKT signaling pathway (Fig. 12C). Protein blotting results showed that p-PI3K/P13K and p-AKT/AKT protein levels were reduced after knockdown of KRT81 (Fig. 12D). These results demonstrated that KRT81 may promote lung adenocarcinoma development through the PI3K/AKT pathway.

Fig. 12

KRT81 promotes EMT and PI3K/AKT pathway activation. A WB detection of EMT-associated proteins (E-cadherin, N-cadherin, and Vimentin) after KRT81 knockdown. B, C KEGG and GSEA enrichment analysis of up-regulated DEGs. D WB detection of p-PI3K, PI3K, p-AKT, and AKT proteins levels in A549 cells after siNC, siRNA2 transfections.*p < 0.05, **p < 0.01

KRT81 promotes tumorigenesis in vivo

A xenograft nude mouse model was established to elucidate the role of KRT81 in LUAD in vivo. Tumor changes were closely monitored after injection of A549 cells from transfected siNC, siRNA2. The results showed that KRT81 knockdown inhibited tumor growth (Fig. 13A). Tumor size and weight were significantly reduced (Fig. 13B). IHC results showed that the expression levels of KRT81, Ki-67 and PCNA proteins, which are closely related to tumor proliferation, were significantly lower in the siRNA2 group than in siNC (Fig. 13C). To further confirm this, we established an in situ model of lung cancer using B6-178 KrasLSL-G12D mice, and collected cancerous tissues and adjacent normal tissues. The results of IHC showed that the expression of KRT81 was significantly up-regulated in mouse lung adenocarcinoma tissues as compared with that in normal mouse lung tissues (Fig. 13D), indicating that KRT81 plays an important role in the progression of lung adenocarcinoma.

Fig. 13

KRT81 promotes in vivo tumorigenesis in nude mouse model and lung in situ model. A Xenograft tumors. B Tumor size and weight. C Expression levels of proliferative proteins KRT81, Ki67 and PCNA detected by IHC. D B6- 178 KrasLSL-G12D in situ lung cancer mouse model was constructed, and the expression of KRT81 was detected by IHC. *p < 0.05, **p < 0.01

Discussion

Lactylation, a newly identified post-translational modification induced by lactate, plays a key role in tumor development [35]. A recent study conducted a global analysis of lactylation in human lungs. Ultimately, 141 proteins modified by lactylation were identified [36]. Lactylation modification affects not only tumor cells in TME but also immune cells. Studies have shown that lactic acid has an immunosuppressive effect on most of the different types of immune cells in TME. Tumor-derived lactic acid promotes the development of myeloid-derived suppressor cells (MDSC) [37]. As an important participant in tumor immunity, the relationship between T cell lactylation and tumor progression has also received widespread attention. Studies have shown that lactate can promote tumorigenesis by regulating MOESIN lactylation and enhancing TGF-β signaling [38]. The oncogenicity of M2 macrophages was found to be positively correlated with the level of histone lactylation [8].

Lactylation modification plays an important role in cancer progression; however, to the best of our knowledge, lactylation-related prognostic gene signatures based on lactylation have not been reported in LUAD.

In this study, we first defined two clusters, C1 and C2, based on 26 LACAGs using an unsupervised consensus clustering method. We found that the two clusters differed significantly in terms of patient survival, with the C1 cluster having a significantly worse prognosis. This finding suggests that abnormal lactylation may contribute to the aggressiveness and progression of LUAD, thereby affecting patient prognosis. We further analyzed KEGG and GSEA enrichment of differential genes in the C1 and C2 cluster, and found that the Cell cycle and ECM-receptor interaction pathways were significantly enriched, and the above results preliminarily explain that the difference in prognosis between the C1 and C2 clusters may be attributed to these two pathways.

Ten genes were finally identified to construct prognostic models: IGFBP1, CYP17A1, DKK1, KRT81, MS4A1, C11orf16, BCAN, FBN2, ANGPTL4, and SERPINB7. Insulin-like growth factor binding protein 1 (IGFBP1), as a member of the family of secreted proteins, has been found in a recent study that IGFBP1 expression is upregulated in tumor cells, which promotes tumor metastasis by enhancing mitochondrial ROS accumulation [39]. The multifunctional enzyme cytochrome P450 (CYP17A1) plays a crucial role in androgen production and is closely related to the development of prostate and breast cancers [40,41,42]. DKK1, as an inhibitor of the Wnt/β-catenin signaling pathway, has emerged as a prognostic marker for a variety of cancers [43]. Yao et al. found that DKK1 significantly promoted NSCLC tumor cell migration, invasion, and angiogenesis [44]. BCAN-NTRK1 is a potent glioma driver gene and therapeutic target. Studies have modeled four relatively rare chromosomal rearrangements of unknown oncogenic potential in human gliomas and found that one of the chromosomal deletions leads to a fusion of BCAN and NTRK1, which promotes the formation of highly aggressive gliomas [45]. As a key component of human elastic fibers, it has been found that the methylation of FBN2 may be closely associated with NSCLC invasion and metastasis [46]. Angiopoietin-like protein 4 ANGPTL4 plays an important role in the tumor microenvironment. Xiao et al. found [47] that in NSCLC cells, ANGPTL4 could promote glutamine consumption and fatty acid oxidation, and knockdown of ANGPTL4 inhibited tumor growth in mice. Serine family B member 7 (SERPINB7) has been associated with the development of psoriasis and long island-type palmoplantar keratosis [48], and a recent study found that SERPINB7 was down-regulated in NSCLC cells after vitamin C (VitC)-treated human NSCLC cell lines. SERPINB7 may serve as a non-small cell lung intravenous vitamin C (IVC) therapy as a predictor of favorable response [49]. Then, we used the TCGA data as a training set and GSE68465 as a validation set, which was categorized into high-risk and low-risk groups based on the median risk score to validate the confidence level of our risk model. For the training set, the overall survival of LUAD patients in the high-risk group was significantly shorter than that in the low-risk group (P < 0.001). The results of the training set were similar to those of the test set (P = 0.012). Multivariate COX regression analysis showed that risk score and clinical stage were independent indicators affecting the prognosis of LUAD patients (P < 0.001). We constructed a bar chart combining risk score and clinical characteristics that could reliably predict the prognosis of LUAD patients. In addition, the AUC of our model ROC curve was more significant than the conventional clinical characteristics of patients.

We performed an immune cell infiltration analysis using CIBERSORT and examined the correlation between immune cell infiltration and risk scores. Resting macrophages (M0) are derived from bone marrow and are usually considered precursors of polarized macrophages [23]. In our study, we found that patients with higher high-risk scores had higher M0 macrophage infiltration scores, suggesting that tumors in high-risk patients may have higher M0 macrophage infiltration. Mast cells are capable of remodeling TME through direct cell-to-cell interactions within the tumor [50]. We found that patients with higher high-risk scores had lower mast cell infiltration scores, suggesting that tumors in high-risk patients may have lower mast cell infiltration.

TMB is defined as the number of somatic mutations per million bases and is often used as a predictive biomarker of immune checkpoint blockade in lung cancer [51]. We analyzed the TMB status of lung adenocarcinoma patients in the high-risk and low-risk groups. TMB was higher in the high-risk group than in the low-risk group. Some mutations were strongly associated with risk scores. For example, TP53, TTN, and CSMD3 mutations were the top three mutations in the high-risk group. TP53 is the most commonly mutated gene in lung adenocarcinoma patients, and TP53 mutations in lung adenocarcinoma can be used as biomarkers for immune checkpoint inhibitors [52]. TTN is associated with elevated TMB in a variety of solid tumors, and TTN mutant phenotypes may be a potential predictive marker for LUAD patients undergoing ICI [52, 53]. CSMD3 was identified as the second most common mutated gene (after TP53) in lung cancer [54]. TIDE is a computational method used to predict ICB response [25]. According to the TIDE prediction results, patients in the high-risk group had higher TIDE values. This finding suggests that high-risk patients have a higher likelihood of tumor immune escape.

This study further explored the potential role of lactylation-related genes by analyzing single-cell data. By single-cell clustering and trajectory analysis, MS4A1 plays a crucial role in B-cell differentiation. A significant correlation between MS4A1 and B cells was also confirmed in the graph of our model gene-immune cell correlation analysis. MS4A1 encodes the surface molecule CD20, which is widely recognized as a B cell lineage marker [55]. CD20 is expressed at different B cell developmental stages but is down-regulated after differentiation of the B cells into plasma cells [56]. This is consistent with our finding that the relative expression of MS4A1 is reduced after B cells differentiate into plasma cells. In lung adenocarcinoma, MS4A1 has been used in the construction of several prognostic models [57, 58].

Historically regarded as a waste product of glycolysis and an important carbon source for cellular metabolism, lactate has recently been found to be an important product affecting tumor proliferation and metastasis [59]. Lung cancer cells can take up lactic acid autonomously and use it as fuel in the body [16]. In our study we found found an abnormal increase in lactate content in lung adenocarcinoma cells and cancer tissues. Lung adenocarcinoma cell lines A549 and H1975 had significantly higher glucose uptake capacity than normal human bronchial epithelioid cell lines. Lactate dehydrogenase (LDHA) is the enzyme responsible for the interconversion of lactate and pyruvate. Increasing evidence also suggests a relationship between upregulation of lactate dehydrogenase A (LDHA) and tumor cell proliferation [60]. Our experimental results showed that LDHA knockdown resulted in a significant reduction in the proliferation, migration and invasion of lung adenocarcinoma cells. In future studies, we intend to use CRISPR technology to knock down or down-regulate specific lactylation sites in tumor cells and observe their specific effects on tumor biological behaviors.

Finally, we analyzed the role of the prognostic model gene KRT81 in LUAD. KRT81 is a type II hair keratin, one of the major hair proteins expressed in the hair cortex [61]. We confirmed that KRT81 expression was significantly elevated in LUAD tissues and cells. We also explored the proliferation, migration, and invasion of KRT81 by CCK-8, clone formation, wound healing, and transwell assays. We found that the knockdown of KRT81 inhibited proliferation, migration, and invasion of LUAD cells. Increases in matrix MMP2 and MMP9 led to degradation of the extracellular matrix (ECM) and basement membrane (BM), which allowed tumor cells to invade other tissues and metastasize tumor cells [62]. Our experimental results revealed that MMP2 and MMP9 were correspondingly reduced after knockdown of KRT81. By flow cytometry analysis, we also found that the downregulation of KRT81 led to G0/G1 phase cycle block and increased apoptosis in LUAD cells. Bcl-2 is a key marker of the endogenous apoptosis pathway, and the Caspase3 protein belongs to the key proteins of the exogenous apoptosis pathway [63]. P53 protein can promote apoptosis of tumor cells through transcriptionally nondependent promotion of apoptosis [64]. The WB results found that that the apoptotic proteins Caspase3 and p53 were significantly up-regulated, and Bcl-2 was significantly down-regulated after the knockdown of KRT81, demonstrating that the down-regulation of KRT81 led to increased apoptosis. EMT is the acquisition of mesenchymal characteristics by epithelial cells, which is closely related to tumor invasion and metastasis [65]. Xu et al. [66]. identified a new set of tumor cells associated with poor LUAD prognosis, namely KRT81 + tumor cells, which may be associated with activation of the EMT and hypoxia pathways. Therefore, we investigated the relationship between KRT81 and EMT proteins and found that knockdown of KRT81 resulted in a significant increase in E-cadherin and a significant decrease in N-cadherin and Vimentin proteins, demonstrating that KRT81 may promote metastasis by regulating EMT. The PI3K/AkT pathway is aberrantly activated in cancer and contributes to tumorigenesis and development [67]. Inhibitors of the PI3K/AkT/mTOR pathway have been widely used in clinical trials in non-small cell lung cancer [68]. Our findings were consistent with the fact that both KEGG and GSEA enrichment analyses showed that KRT81 up-regulation was able to activate the PI3K/AKT signaling pathway. The results of WB showed that knockdown of KRT81 resulted in a significant decrease in PI3K/AKT phosphorylation, and KRT81 may promote cell proliferation, migration, and invasion by activating the PI3K/AKT pathway. The PI3K/AKT pathway can be aberrantly activated by a variety of mechanisms, including different genomic alterations such as mutations in PIK3CA, phosphatase and tensin homolog (PTEN), AKT, TSC1 and mechanistic target of rapamycin (mTOR) [67]. KRT81 may be involved in the interaction of upstream and downstream proteins in the PI3K/AKT signaling pathway, and the specific mechanism needs to be further explored. In addition, whether KRT81 may affect the development of lung adenocarcinoma through other signaling pathways is also the direction of our further research. Further in vitro and in vivo experiments are needed to validate the mechanism of KRT81 in upstream and downstream protein regulation of EMT and PI3K/AKT pathways.

Finally we found that KRT81 knockdown significantly inhibited tumor growth by a xenograft nude mouse model. Ki67, PCNA proteins are markers of tumor cell proliferation and are commonly used to assess the growth fraction of tumor cells [69]. The IHC results of our xenograft nude mouse model showed that Ki67 and PCNA protein expression levels were significantly reduced after KRT81 knockdown. DuPage et al. produced lung adenocarcinomas in KrasLSL-G12D mice by dropping adenoviruses through tracheal intubation (70). We used the same method to establish an in situ model of lung cancer in B6-178 KrasLSL-G12D mice and collected cancerous tissues and adjacent normal tissues. The results of IHC showed that the expression of KRT81 was significantly up-regulated in tumor tissues, suggesting that KRT81 plays an important role in the progression of lung adenocarcinoma. In the future, we need to further investigate the mechanism by which KRT81 affects LUAD progression. Knockdown of KRT81 or application of specific downstream pathway-targeted drug inhibitors may be potential targets for the treatment of lung cancer. In conclusion, our study found that KRT81 can be used as a promising biomarker to predict the prognosis and tumor progression of LUAD, and may serve as a potential clinical therapeutic target to guide clinical decisions.

Although our findings were validated in an independent cohort, there are some limitations. First, our study was a retrospective study based on the publicly available TCGA database, and the prognostic model, although validated in GSE68465, requires more subsequent clinical prospective studies. We intend to collect more clinical samples for validation and extend our study to other tumors. Second, the detailed mechanisms of the interaction between the lactated prognostic model and TMB, the immune cell infiltration, are not fully understood, and more experiments are needed to explore and validate the relationship and molecular mechanisms. Third, further experiments are needed to explore the potential mechanisms of the other nine lactylation-related genes in the prognostic model. In our next study, we intend to focus on finding potential downstream target proteins of histone lactylation and how it affects key signaling pathways (e.g., PI3K/AKT, MAPK) that play a critical role in tumor growth and metabolic regulation through phosphorylated proteomics. Finally, although we have demonstrated in in vitro and in vivo experiments that KRT81 can inhibit the proliferation, migration, and invasion of lung adenocarcinoma cells, the specific signaling pathway mechanism and whether it affects tumor progression through lactylation modification need further investigation.

Conclusions

In summary, we constructed a prognostic model of 10 lactylation-associated genes in lung adenocarcinoma, which has high accuracy in predicting overall patient survival. The prognostic risk score was closely associated with common clinical features such as immune cell infiltration, immune-related function, TMB, and anticancer drug sensitivity. In addition, we preliminarily verified by in vitro in vivo experiments that KRT81 promotes lung adenocarcinoma proliferation, migration, and invasion and that knockdown of KRT81 blocks the G0/G1 phase cycle and increases apoptosis. KRT81 may play an oncogenic role in LUAD through EMT and PI3K/AKT pathways. These studies help to enhance clinical treatment strategies and guide future precision therapies.