A total of 4 LUAD transcriptome-related datasets, including the scRNA-seq dataset [GSE131907 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131907) (15)] and three bulk transcriptomic cohorts [GSE30219 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30219) (16), platform GPL570; GSE26939 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE26939) (17), platform GPL9053; and GSE72094 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72094) (18), platform GPL15048], were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. Detailed sample information of these datasets is presented in Table SI. In addition, the gene expression profiles and corresponding clinical information of the TCGA-LUAD cohort (https://portal.gdc.cancer.gov/projects/tcga-luad) (19) were downloaded from the TCGA (https://tcga-data.nci.nih.gov/tcga/) database. TCGA-LUAD contains an expression matrix of 585 tumor samples and the overall survival (OS) data of 572 tumor samples. Specifically, 513 patients with both gene expression data and OS information were used to construct the prognostic model, and the GSE30219, GSE26939 and GSE72094 datasets were used for external validation.

For bulk transcriptomic analyses, TCGA-LUAD and the GEO cohorts were included if they comprised histologically confirmed LUAD and provided genome-wide gene expression profiles together with OS data and basic clinicopathological information. Cohorts without survival data, with incomplete clinical annotations or with data quality issues were not considered. For scRNA-seq analysis, GSE131907 was selected as it profiled treatment-naïve patients with LUAD with well-annotated tumor, normal lung and lymph node samples, and provided sufficient T cell coverage for robust identification of Treg-related transcriptional programs.

In detail, TCGA-LUAD contains data for 276 female and 237 male patients with a median age of 66 years (range, 33–88 years). GSE30219 contains data for 43 female and 250 male patients with a median age of 63 years (range, 15–84 years). GSE26939 contains data for 62 female and 53 male patients with a median age of 65 years (range, 41–90 years). GSE72094 contains data for 222 female and 176 male patients with a median age of 70 years (range, 38–89 years). GSE131907 contains data for 17 female and 27 male patients with a median age of 63 years (range, 39–81 years).

The ‘Seurat’ (v.4.1.0) R package (R Foundation for Statistical Computing, Vienna, Austria; http://www.r-project.org/) was employed for the data processing of GSE131907. Low-quality cells and genes were filtered out using 200<nFeature_RNA <6,000, nCount_RNA <60,000 and percent.mt <20% as the criteria (15), otherwise the default was used for the other parameters. Subsequently, distinct cell clusters were generated using uniform manifold approximation and orojection clustering (resolution=1). These cell clusters were then annotated according to corresponding marker genes (15).

Treg markers were selected in GSE131907 according to the following criteria: avg_|log2FoldChange| ≥1.5 (Tregs vs. the other cells) and pct.m >10% (the markers expressed in ≥10% of Tregs). From these differentially expressed genes, the feature genes associated with prognosis were selected in TCGA-LUAD using univariate Cox analysis with P<0.05 as the significance cut-off. Subsequently, LASSO regression was performed to select the prognostic genes using the ‘glmnet’ (4.1–4) R package. The risk score formula was follows: Risk score=(Expression_1 × Coefficient_1) + (Expression_2 × Coefficient_2) + … + (Expression_n × Coefficient_n), where n is the number of prognostic genes.

A total of 513 patients with LUAD were then divided into the high- and low-risk groups based on the median risk score. A Kaplan-Meier curve was generated to compare the OS difference between the high- and low-risk groups using the ‘survminer’ (v. 0.4.9) R package (P<0.05). Meanwhile, the ‘survivalROC’ (v. 1.0.3.1) R package was used to construct the 1-, 2- and 3-year receiver operating characteristic (ROC) curves. The prognostic model was then verified using the GSE30219, GSE26939 and GSE72094 datasets using Kaplan-Meier and ROC curves. In consideration that the external validation cohorts were profiled on different platforms, which introduce distributional and scale shifts in gene expression and in the resulting risk scores, using a single absolute cut-off across heterogenous cohorts can be poorly calibrated. Therefore, a cohort-specific optimal cut-off within each validation set was used, which is standard for platform-heterogenous validations. All bulk RNA-seq/microarray cohorts were analyzed independently without cross-cohort merging. Therefore, batch-effect correction across different datasets or platforms was not applicable.

To assess the effect of clinical indicators on the prognosis of patients with LUAD, 5 clinical indicators [age, sex, tumor (T) stage, lymph node (N) stage and metastasis (M) stage] and risk score were evaluated. Subsequently, independent prognostic indicators were identified using univariate and multivariate Cox analyses (P<0.05). Finally, the differences in risk score in different clinical subgroups were analyzed (Wilcox rank-sum test; P<0.05).

To evaluate the TME of patients with LUAD, the CIBERSORT algorithm was utilized. The infiltration ratio differences of 22 immune cells were compared in TCGA-LUAD [adjusted (adj.) P<0.05]. Subsequently, 7 immune-system-related metagene clusters [hemopoietic cell kinase, IgG, interferon, lymphocyte-specific kinase (LCK), major histocompatibility complex (MHC)-I, MCH-2 and STAT1] were obtained from the published literature (20) and the immune scores of each sample were calculated. In addition, the T-cell receptor (TCR) richness and TCR Shannon diversity were determined in the low- and high-risk groups. TCR richness refers to the number of distinct productive TCR clonotypes detected in a sample; it reflects how many different T cell clones are present (breadth of the repertoire). TCR Shannon diversity summarizes both the number of clonotypes and how evenly they are distributed. Higher values indicate that numerous different clones are present and none of them overwhelmingly dominates (greater breadth and evenness) (21).

In TCGA-LUAD, genes that were moderately and highly correlated with risk scores were selected based on the Spearman's correlation (|cor| >0.5; adj.P<0.05). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to evaluate the function of these genes using the ‘clusterProfiler’ (v. 4.2.2) R package (adj.P<0.05) (22).

The tumor and para-carcinoma tissues of 5 patients with LUAD (2 male and 3 female patients; median age, 60 years; age range, 52–65 years) were collected from the Shanghai General Hospital Jiuquan Branch (Jiuquan People's Hospital; Jiuquan, China) from March 2025 to April 2025. The inclusion criteria for patient selection were as follows: i) Histopathologically confirmed primary LUAD; ii) patients who underwent curative-intent surgical resection and could provide paired LUAD and adjacent non-tumorous lung tissues; iii) no neoadjuvant chemotherapy, radiotherapy, targeted therapy or immunotherapy before surgery; and iv) basic clinicopathological information were available and the patient provided written informed consent. The exclusion criteria were as follows: i) History of other malignant tumors; ii) recurrent or metastatic LUAD at the time of surgery; iii) severe concomitant infectious, autoimmune or systemic inflammatory diseases that could markedly affect immune status; and iv) insufficient tissue quality or quantity for RNA and protein extraction.

Total RNA from all tissues was extracted using the E.Z.N.A.^® Total RNA Kit I (cat. no. R6834-01; Omega Bio-tek, Inc.). The synthesis of cDNA was then performed using the 5X Evo M-MLV RT Master Mix (cat. no. AG11728; Hunan Accurate Bio-Medical Technology Co., Ltd.) according to the manufacturer's instructions. Subsequently, qPCR was performed using the 2X SYBR Green Premix Pro Taq HS qPCR Tracking Kit (cat. no. AG1170; Hunan Accurate Bio-Medical Technology Co., Ltd.) on the StepOne Plus™ Real-Time PCR System (Thermo Fisher Scientific, Inc.) under the following thermocycling conditions: Initial denaturation at 95°C for 30 sec, followed by 40 cycles of denaturation at 95°C for 5 sec and annealing at 60°C for 30 sec. The primer sequences for all prognostic genes are listed in Table SII. The data was then exported and analyzed using the 2^−ΔΔCq method (23), with β-actin as the endogenous control.

Fresh-frozen LUAD tumor and adjacent normal tissues were homogenized on ice in RIPA buffer (MilliporeSigma). Lysates were centrifuged at 8,000 × g at 4°C for 15 min and supernatants were collected. Protein concentrations were determined using the BCA assay following the manufacturer's instructions (Thermo Fisher Scientific, Inc.) and quantified on a microplate reader at 562 nm. Aliquots containing equal amount of protein (20 µg per lane) were mixed with 5X loading buffer (Beyotime Biotechnology), heated at 98°C for 5 min, rapidly cooled on ice and loaded onto 4–12% SDS-PAGE gels (Thermo Fisher Scientific, Inc.). Proteins were then transferred to PDVF membranes, followed by blocking with 5% BSA (MilliporeSigma) at room temperature for 1 h. Primary antibodies against centromere protein M (CENPM) (cat. no. 19840-1-AP), zinc finger protein 101 (ZNF101) (cat. no. 25599-1-AP), melanoma-associated antigen H1 (MAGEH1) (cat. no. 12424-1-AP), Ikaros family zinc finger protein 4 (IKZF4) (cat. no. 31094-1-AP) and β-actin (cat. no. 20536-1-AP) purchased from Proteintech Group, Inc., were diluted to 1:1,000 and incubated with the membranes overnight at 4°C. β-actin was used as the loading control. After washing with PBST, membranes were incubated with HRP-conjugated goat anti-rabbit recombinant secondary antibody (H+L) (cat. no. RGAR001; Proteintech) diluted to 1:5,000 at room temperature for 1 h. After washing with PBST, signals were developed using ECL chemiluminescence (Beyotime Biotechnology) and imaged on a chemiluminescence system. Band intensities were quantified in ImageJ (version 1.53a; National Institutes of Health) and target proteins were normalized to β-actin to obtain relative expression.

All statistical analyses were performed using R software and GraphPad Prism 9.0 (Dotmatics). Continuous variables between two groups in the TCGA and GEO cohorts were compared using the Wilcoxon rank-sum test. For RT-qPCR and western blot experiments, expression levels of key genes are presented as mean ± SD and differences between tumor and adjacent tissues were assessed using the paired Student's t-test (n=5 for each group). Survival curves were estimated by the Kaplan-Meier method and compared using the log-rank test. P<0.05 was considered to indicate a statistically significant difference. For multiple testing in the differential expression and enrichment analyses, the Benjamini-Hochberg method was used to control the false discovery rate and adj. P<0.05 was considered to indicate a statistically significant difference.

First, quality control of all samples in GSE131907 was performed (Table SIII). Subsequently, the top 2,000 variable genes and the top 50 principal components were used for dimensionality reduction and all high-quality cells were grouped into 29 clusters (Fig. 1A). After annotation using canonical cell markers, 13 cell types, comprising epithelial cells, fibroblasts, T cells, Tregs, natural killer (NK)T cells, NK cells, B cells, plasma cells, monocytes, macrophages, MAST cells, classical dendritic cells (DCs) and tumor specific DCs, were identified (Fig. 1B). The expression levels of the canonical marker genes in each cell type are displayed in a bubble plot shown in Fig. 1C. The average proportion of Tregs in the nLNs of patients with LUAD was 2.4346%, which was notably higher than that of the other five groups (Fig. 1D and Table SIV). In addition, the tumor lung (average proportion, 2.0831%) of patients with LUAD had markedly more Tregs than normal lung (average proportion, 0.6636%). By comparing gene expressions between Tregs and the other 12 cell types, a total of 326 Treg markers were obtained for the further analysis (Table SV).

A total of 17 Treg markers [C-C motif chemokine receptor 6 (CCR6), zinc finger CCCH-type containing 12D (ZC3H12D), ankyrin repeat and SOCS box protein 2 (ASB2), IKZF4, ZNF101, T cell receptor β constant 1, signaling lymphocytic activation molecule family member 1, CD28, cytokine inducible SH2 containing protein (CISH), CD5, MAGEH1, glucocorticoid induced 1 (GLCCI1), lymphocyte transmembrane adaptor 1, baculoviral IAP repeat containing 3 (BIRC3), interleukin 1 receptor type 2 (IL1R2), pituitary tumor-transforming gene 1 protein (PTTG1) and CENPM] associated with prognosis were identified based on the univariate Cox analysis (Fig. 2A). To obtain more robust Treg markers, LASSO analysis was employed and 13 prognostic Treg features with non-zero coefficients, including CENPM, PTTG1, IL1R2, GLCCI1, BIRC3, MAGEH1, CD5, CISH, ZNF101, IKZF4, ASB2, ZC3H12D and CCR6, were selected (λ.min=0.005779793) (Fig. 2B and Table SVI). The risk coefficients of CENPM, PTTG1, IL1R2 and BIRC3 were positive and were thus risk factors in the model. The risk coefficients of GLCCI1, MAGEH1, CD5, CISH, ZNF101, IKZF4, ASB2, ZC3H12D and CCR6 were negative and were thus protective factors in the model. In addition, the expression levels of these 13 genes in the Treg cells of metastatic (m)LN & tumor lung of advanced-stage LUAD (tL/B), pleural effusion, nLN, nLung, tLung and mBrain groups were detected. It was found that, except CENPM, the other 12 genes were differentially expressed across different anatomical groups. For example, the expression levels of MAGEH1 and CCR6 were the highest in the tLung group, but the abundances of IKZF4, BIRC3 and PTTG1 were the highest in the mLN & tL/B group (Fig. S1). Subsequently, based on the median risk score, 513 patients with LUAD were grouped into the high- (n=256) and low- (n=257) risk subgroups. As shown in Fig. 2C, the high-risk group had more deaths than the low-risk group. The survival probability of patients in the low-risk groups was significantly higher than that in the high-risk group (Fig. 2D). Meanwhile, the area under ROC curves (AUC) for 1-, 2- and 3-years of the risk score model were 0.719, 0.715 and 0.711, respectively, demonstrating the predictive power of the prognostic model (Fig. 2E).

To assess the universality of the prognostic model, three independent datasets were examined, including GSE30219, GSE26939 and GSE72094. In GSE30219, 193 patients with LUAD were classified into the high- (n=150) and low- (n=143) risk subgroups based on the optimal threshold of −5.823151. The OS of the low-risk subgroup was significantly higher than that of the high-risk subgroup (P=0.01; Fig. 3A), and the 1-, 2- and 3-year AUC values of the risk score model were 0.611, 0.626 and 0.615, respectively (Fig. 3B). In GSE26939, 115 patients with LUAD were grouped into the high- (n=58) and low- (n=57) risk subgroups based on the optimal threshold of −0.932411. The patients with LUAD in the low-risk subgroup had a higher OS rate than those in the high-risk group (Fig. 3C) and the AUC values of the risk score model for 1-, 2- and 3-year were 0.684, 0.647 and 0.660, respectively (Fig. 3D). In GSE72094, 398 patients with LUAD were divided into the high- (n=200) and low- (n=198) risk subgroups based on the optimal threshold of −2.738428. The OS of the low-risk subgroup was significantly higher than that of the high-risk subgroup (P=0.00036; Fig. 4A) and the 1-, 2- and 3-year AUC values of the risk score model were 0.616, 0.657 and 0.663, respectively (Fig. 4B).

In the TCGA-LUAD cohort, the risk score, T stage, N stage and M stage were determined to be associated with the prognosis of patients with LUAD using univariate Cox analysis (Fig. 5A). Moreover, the risk score, T stage and N stage were selected as independent prognostic indicators (Fig. 5B). In addition, an association was demonstrated between the risk score and clinical characteristics of patients with LUAD. Specifically, the risk score of patients with LUAD with stage T3-T4 and N1-N3 were markedly higher than those patients with stage T1-T2 and N1-N3, respectively (Fig. 5C). This indicates that the risk score may be relevant to LUAD development and progression. Moreover, male patients with LUAD had notably higher risk scores than female patients with LUAD (Fig. 5C); however, there was no significant difference in the risk score between the M0 and M1 stages (Fig. 5C).

To assess differences in the TME of the low- and high-risk groups, the infiltration of 22 immune cells was first evaluated (Fig. 6A). The results revealed that the infiltration of M0 macrophages and activated mast cells in the high-risk group was notably higher, whilst the infiltration of plasma cells and Tregs was markedly lower, compared with the low-risk group (Fig. 6B). Moreover, the high-risk group had lower immune, stromal and ESTIMATE scores than the low-risk group (Fig. 6C), further supporting that the immune environment of the two groups was different.

To further assess the potential effect of the prognostic Treg markers on immunotherapeutic treatment, the correlations between risk score and inflammatory metagenes were then determined. The findings demonstrated that risk score had significant weak correlations with 4 metagene clusters, including IgG (cor=−0.22), interferon (cor=0.22), LCK (cor=−0.31) and MHC-II (cor=−0.17) (Fig. 7A). In addition, the TCR repertoire analysis revealed that the high-risk group had lower TCR Shannon diversity (P<0.0001; Fig. 7B) and lower TCR richness (P<0.0001; Fig. 7C) than the low-risk group, suggesting a reduced capacity to recognize neoantigens in the high-risk group. Furthermore, considering the association between tumor mutation and immune response, the tumor mutation burden (TMB) score was assessed and the results revealed that the TMB of the high-risk group was markedly higher than that of the low-risk group (Fig. 7D). These findings indicate that the prognostic Treg-related markers may affect immune response and therapeutics in patients with LUAD.

To further evaluate the role of prognostic Treg-related genes, genes that were strongly correlated with the Treg-related risk score were identified. After correlation analysis, 33 genes were determined to be negatively correlated with the risk score and 4 genes were positively correlated (|cor| >0.5; adj.P<0.05; Fig. 8A and Table SVII). These genes were mainly associated with metabolic biological processes, such as ‘ADP catabolic process’, ‘pyridine nucleotide catabolic process’ and ‘purine nucleoside diphosphate catabolic process’ (Fig. 8B and Table SVIII) and KEGG pathways such as ‘Glycolysis/Gluconeogenesis’, ‘HIF-1 signaling pathway’, ‘Carbon metabolism’ and ‘Pyruvate metabolism’ (Fig. 8C and Table SIX). These results suggest the potential role of the prognostic Treg markers in metabolism.

The expression levels of these genes in the TCGA-LUAD cohort were examined and it was found that the CENPM, PTTG1, ZC3H12D and ZNF101 expression levels were significantly upregulated, while CISH, GLCCI1, IKZF4 and MAGEH1 were significantly downregulated in the tumor tissues (Fig. 9A). Clinical samples were then collected to assess the expression levels of these 8 genes and it was found that the levels were altered consistently with the TCGA-LUAD results (Fig. 9B). The protein expression levels of CENPM, ZNF101, MAGEH1 and IKZF4 in the clinical samples were also examined. The results showed that CENPM and ZNF101 were significantly increased, while MAGEH1 and IKZF4 were reduced in the tumor tissues (Fig. 9C and D). These findings indicate that the expression levels of these genes were reproducible in the independent clinical LUAD samples.