Four datasets—GSE32863 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32863) (10), GSE43458 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE43458) (11), GSE75037 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE75037) (12), and GSE72094 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72094)—were acquired from the Gene Expression Omnibus (GEO) repository. The GSE32863 dataset comprised microarray expression profiling of 58 pairs of LUAD and adjacent non-tumor lung tissues on the Illumina platform, which were log2-transformed and subjected to robust spline normalization with the lumi package (13). Gene expression profiles in the GSE43458 (containing 80 LUAD and 30 normal lung tissues) and GSE72094 (containing 442 LUAD specimens) datasets on the Affymetrix platform were normalized, summarized, and adjusted for background with robust multichip average, together with log2 transformation utilizing integrated BRB-Array Tools. The GSE75037 dataset included expression profiling of 83 pairs of LUAD and matched adjacent non-malignant lung tissues on the Illumina platform. Bead-level data were background-adjusted and summarized with model-based background correction for the BeadArrays algorithm (14), followed by quantile normalization and log2 transformation. Moreover, RNA sequencing data of 510 primary LUAD and 58 normal tissues were obtained from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/), and count data were normalized with Trimmed Mean of M-values utilizing the edgeR package (15). Gene methylation data of 473 LUAD and 34 normal tissues were also downloaded from TCGA on the basis of the Infinium HumanMethylation450 BeadChip array.

The expression levels of genes between LUAD and normal lung tissues in each dataset (GSE32863, GSE43458, GSE75037, and TCGA-LUAD) were compared with Student’s t-test utilizing the limma package. The p-value was adjusted through the false discovery rate. Differentially expressed genes (DEGs) were identified in accordance with adjusted p < 0.05 (16).

The co-expression network of DEGs shared by four datasets was constructed via the weighted gene correlation network analysis (WGCNA) package (17). The soft-thresholding power β of 4 (scale-free R² of 0.9) was utilized for achieving an adjacency matrix with scale-free topology. Then, the topological overlap matrix was established for determining the connectivity and dissimilarity of the co-expression network. With the use of the dynamic tree cut approach, co-expression modules were clustered. Pearson’s test was utilized for computing the relationships between module eigengene (ME), which is defined as the first principal component of a given module and clinical traits.

Based on the markers of immune infiltrates, the enrichment score computed by single-sample gene set enrichment analysis (ssGSEA) was utilized for denoting the relative infiltration of each immune cell type (18). Univariate Cox regression analysis was implemented for inferring the relationships between immune infiltrates with overall survival (OS).

In accordance with the infiltration levels of risky immune cell types, a hierarchical clustering dendrogram was adopted for classifying LUAD cases as the optimal number of clusters in TCGA-LUAD and GSE72094 datasets. The Elbow method was used to validate the number of clusters. OS between clusters was evaluated via Kaplan–Meier approach together with the log-rank test.

Functional enrichment analysis of module genes in certain Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was evaluated with the clusterProfiler package (19). GO terms comprised biological processes (BPs) and cellular components (CCs) together with molecular functions (MFs). Gene set enrichment analysis (GSEA) software was adopted for determining gene sets with significant differences between groups (20). The pathway activity in each sample was estimated utilizing gene set variation analysis (GSVA) in accordance with the enrichment of gene sets (18). Terms with adjusted p < 0.05 were considered significant enrichment.

The IC50 value of therapeutic agents was estimated with the pRRophetic approach (21), which can denote the effect of an agent in mitigating specific biological or biochemical functions.

Differentially methylated sites were determined with the ChAMP package (22). First, quality control was implemented with the champ.QC function. Second, the β value was normalized via champ.norm function. Differentially methylated genes were screened according to the CpG site with p ≤ 0.05.

TMB denotes the number of somatic mutations in each megabase of the query genome sequence (23), which was computed from the number of variants out of the total length of the human exons (38 million). A waterfall diagram was conducted with the maftools package for assessing the number of somatic mutations across TCGA-LUAD specimens (24).

Module genes with methylation were imported into the online tool STRING (https://string-db.org) (25). The protein–protein interaction (PPI) network with a combined score of ≥0.7 (high confidence) was constructed and visualized with Cytoscape software (26). Genes showing more interactions with others were defined as key module genes.

Key module genes were utilized to establish a multivariate Cox regression model with a survival package, and the risk score was computed by combining the coefficient and mRNA level of each key module gene, following the formula: risk score = β₁x₁ + β₂x₂ + ··· + β_ix_i (where β_i represents the coefficient of gene i and x_i represents the mRNA level of gene i). With the median value, LUAD cases were classified as high- and low-risk groups. Survival state was compared between groups. Afterward, a key module gene-based nomogram was developed via the rms package for estimation of 5- and 8-year OS. The prediction accuracy of the nomogram was assessed via calibration curves.

Human lung normal epithelial cell line (BEAS-2B) and LUAD cell lines (NCI-H1975, NCI-H1395, A549, and PC9) were cultivated in Roswell Park Memorial Institute (RPMI) 1640 medium (Thermo Scientific, MA, USA) plus 10% fetal bovine serum (FBS) at 37°C in a humidified incubator with an atmosphere of 5% CO₂. The coding sequence of full-length human EXO1 was subcloned into pLenti-CMV-GFP-Puro (#17448; Addgene, MA, USA) following the manufacturer’s instructions. A549 cells were treated with pLenti-CMV-GFP-Puro-EXO1 or matched empty vector as negative control (NC) under 8 µg/ml of polybrene for 12 h. NCI-H1975 cells were transfected with 10 nmol/L short interfering RNAs (siRNAs) against EXO1 (siRNAs 1-5, Hanyin, Shanghai, China) or a scrambled negative control siRNA (siRNA NC; Hanyin) via Lipofectamine 2000 (Invitrogen, MA, USA) in line with the manufacturer’s instructions. The sequences included the following: EXO1 siRNA1, 5′-GAAGAGAAGTTTCGTTACA-3′; EXO1 siRNA2, 5′-GTTGGCCTATCTTAACAAA-3′; EXO1 siRNA3, 5′-GAAGTAGAGAGATCTAGAA-3′; EXO1 siRNA4, 5′-TGACTACAATCCAGACACT-3′; EXO1 siRNA5, 5′-CCATTTCACCACCCACTTT-3′. Quantitative real-time PCR was adopted to examine EXO1 expression after transfection.

Total RNA was isolated with TRIzol reagent (Thermo Scientific) following the manufacturer’s instructions and reverse-transcribed into cDNA with a reverse transcription kit (Takara, Dalian, China). Quantitative PCR was conducted with SYBR Green PCR Mix (Invitrogen, MA, USA) following the primers: EXO1 forward, 5′-TGAGGAAGTATAAAGGGCAGGT-3′; EXO1 reverse, 5′-AGTTTTTCAGCACAAGCAATAGC-3′; GAPDH forward, 5′-GGAGCCAAAAGGGTCATCATCTC-3′; GAPDH reverse, 5′-GAGGGGCCATCCACAGTCTTCT-3′; PD-L1 forward, 5′-TGGCATTTGCTGAACGCATTT-3′; PD-L1 reverse, 5′-TGCAGCCAGGTCTAATTGTTTT-3′. EXO1 mRNA level was computed with the 2^−ΔΔCt approach and normalized to GAPDH.

Transfected cells were seeded into 96-well plates (1 × 10³ cells/well) and cultivated for 5 days. Cellular proliferation was measured at 0, 1, 2, 3, and 4 days with Cell Counting Kit-8 (CCK-8; Beyotime, Shanghai, China). Cultures were exposed to CCK-8 reagent for 2 h. Afterward, the optical density (OD) value at 450 nm was measured.

Migration and invasion were examined utilizing transwell chambers (Corning, MA, USA). For migration assay, transfected cells were cultivated in serum-free RPMI 1640 medium for 15 h and then loaded into the upper chamber. The lower chamber was filled with RPMI 1640 medium with 20% FBS. For invasion assay, the upper chamber was coated with 25 µg of Matrigel (BD, CA, USA) and then seeded with 5 × 10⁴ cells. RPMI 1640 medium was added to the lower chamber. Following 24-h incubation, migratory and invasive cells were fixed with glutaraldehyde and dyed with crystal violet.

The cells were lysed in radioimmunoprecipitation assay (RIPA) buffer (Beyotime, China) plus a protease inhibitor cocktail, followed by protein quantification utilizing the bicinchoninic acid (BCA) method. The lysates were subjected to sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and transferred onto the polyvinylidene difluoride (PVDF) membrane. The membrane was immersed in 5% skim milk and probed with PD-L1 (1/1000; ab213480; Abcam, Waltham, MA, USA) and GAPDH (1/2500; ab9485; Abcam) overnight at 4°C. Afterward, a secondary antibody (1/5000; ab7097; Abcam) was applied to the membrane. Protein bands were visualized utilizing an enhanced chemiluminescence (ECL) substrate (Solarbio, Beijing, China).

In accordance with the manufacturer’s specification, soluble PD-L1 (sPD-L1) level was detected utilizing a human PD-L1 ELISA kit (ED-16408; Xiamen Lunchangshuo Biotechnology Co., Ltd., China).

Data were expressed as mean ± standard deviation. Statistical analysis was conducted with GraphPad Prism 8.0.1 (GraphPad Software, CA, USA) or R language 4.0.5. The difference was evaluated for significance with Student’s t-test, Wilcoxon rank-sum test, and one-way or two-way analysis of variance. Pearson’s or Spearman’s test was utilized for correlation analysis. Receiver operating characteristic (ROC) curve of each key module gene was drawn, and the area under the curve (AUC) was computed for estimating the prediction accuracy utilizing the survivalROC package. p < 0.05 denoted statistical significance.

The current study retrospectively collected four LUAD datasets (GSE32863, GSE43458, GSE75037, and TCGA-LUAD). A total of 16,667 DEGs were selected by comparing LUAD with controls in four datasets ( Figure 1A ). Among them, 3,150 upregulated genes and 2,545 downregulated genes were shared by four datasets ( Figures 1B, C ). After detection of outliers ( Figure 1D ), the WGCNA approach was adopted to determine different modules for shared DEGs in accordance with scale-free fit index β = 4 (scale-free R² = 0.9; Figure 1E ). Average linkage clustering generated nine co-expression modules containing 4,849 shared DEGs ( Figures 1F, G ). Module genes were notably in relation to immunity (neutrophil-mediated immunity, regulation of T-cell activation, IL-17 signaling pathway, etc.) and tumorigenic (autophagy, Wnt/cAMP/ErbB/mTOR/p53/PI3K-Akt/Ras/TGF-beta signaling pathways, etc.) processes and pathways ( Figures 1H, I ). The detailed information is listed in Supplementary Table 1 .

Infiltration levels of immune cell types notably differed between LUAD and control lung tissues. B cells, follicular helper T (Tfh) cells, regulatory T cells (Tregs), and T helper 2 (Th2) cells were abundant in LUAD, with other immune cell types (neutrophils, mast cells, etc.) abundant in normal lungs ( Figure 2A ). Eight immune cell types prominently correlated to LUAD cases’ OS, of which B cells, eosinophils, immature dendritic cells (iDCs), mast cells, T cells, effective memory T (Tem), and Tfh cells acted as risk factors of OS, with Th2 cells as a protective factor ( Figure 2B ). On the basis of infiltration of the above seven risky immune cell types, we classified LUAD patients in TCGA-LUAD cohort as two molecular subtypes, namely, C1 and C2 ( Figure 2C ). The Elbow method confirmed the accuracy of this molecular classification ( Figure 2D ). In Figure 2E , OS outcomes observably differed between subtypes, with the survival advantage in C2. In addition, we externally verified this immune cell infiltration-based classification in the GSE72094 dataset. As expected, similar results were observed, as shown in Figures 2F–H . Further analysis found that the MEblue module exhibited the strongest association with LUAD and C1 subtype, with the strongest correlation between the MEturquoise module and normal lung tissues and the C2 subtype ( Figure 2I ). For DEGs in the MEblue module, module membership dramatically correlated to gene significance for the C1 subtype ( Figure 2J ). Moreover, for DEGs in the MEturquoise module, module membership was markedly associated with gene significance for the C2 subtype ( Figure 2K ). Thus, MEturquoise and MEblue were regarded as the key modules.

GSEA demonstrated that base excision repair, cell cycle, DNA replication, Fanconi anemia pathway, and one carbon pool of folate exhibited a prominent activation in TCGA-LUAD than in normal lung tissues ( Figure 3A ). In both TCGA-LUAD and GSE72094 datasets, in comparison to the C2 subtype, higher activity of basal transcription factors, Fanconi anemia pathway, mismatch repair, nucleotide excision repair, and spliceosome was observed in the C1 subtype ( Figure 3B and Supplementary Figure 1 ). MHC molecules displayed different expressions between two subtypes across TCGA-LUAD samples ( Figure 3C ). Moreover, we compared the difference in the abundance of immune cells between the two subtypes. Overall, immune cell infiltrations were lower in C1, but with higher infiltration levels of Tgd, NK CD56dim cells, and Th2 cells in C2 ( Figure 3D ; and Supplementary Table 2 ). Figure 3E depicts the notable interactions between distinct immune cell types and their prognostic implications across TCGA-LUAD. Afterward, we separately analyzed the correlations between immune cell types in each subtype. As a result, prominent relationships between immune cell types were observed in both C1 and C2 subtypes ( Figure 3F ). TMB may infer the survival benefit of LUAD patients with ICB treatment (27). A higher TMB score was observed in the C1 subtype, which was predictive of clinical benefit with ICB ( Figure 3G ). Immune checkpoints ADORA2A, BTLA, CD4, and TIM3 exhibited higher activity in C2, with lower activity of LAG3, PD-L1, and IDO1 in C1 ( Figure 3H ), indicating the prominent differences in immune checkpoints between subtypes. Moreover, patients in the C1 subtype were more likely to respond to JNK inhibitor VIII, BI-D1870, RO-3306, A-443654, VX-680, SL-0101-1, S-trityl-l-cysteine, BI-2536, itomycin C, CCT018159, KU-55933, and docetaxel ( Figure 3I ).

In TCGA-LUAD dataset, 148,187 differentially methylated sites were determined between LUAD and normal lung tissues ( Figure 4A ). Figure 4B visualizes the distribution of differentially methylated sites on the chromosomes. Among differentially methylated genes, 3,340 were module genes whose expression was altered by methylation ( Figure 4C ). Additionally, the genetic mutation was evaluated. Among 567 LUAD samples, 537 (94.71%) occurred genetic mutation, indicating the widespread mutation ( Figure 4D ). According to the mutation frequency, we depicted the top 30 module genes across LUAD samples, including MUC16, NAV3, COL11A1, ANK2, PCDH15, DNAH9, PTPRD, CDH10, HMCN1, RELN, MXRA5, CACNA1E, VCAN, CSMD2, PKHD1L1, ADGRL3, ASPM, SYNE1, BRINP3, RIMS2, COL6A3, LAMA2, NEB, TSHZ3, AHNAK, MYH2, ADGRB3, NLRP3, DNAH3, and PCDH10.

Among methylated genes in two key modules MEturquoise and MEblue, 945 were significantly correlated to LUAD cases’ OS. GSVA revealed the top five upregulated pathways (mismatch repair, base excision repair, nucleotide excision repair, nucleotide excision repair, and one carbon pool by folate) and top five downregulated pathways (adherens junction, cAMP signaling pathway, complement, and coagulation cascades, TGF-beta signaling pathway and signaling pathway regulating pluripotency of stem cells) ( Figure 5A ), which were linked to methylated genes ( Figure 5B ). The PPI network showed their interactions ( Figure 5C ). The expression of the above genes between LUAD and controls is depicted in Figure 5D . The five genes with the highest connectivity were defined as key module genes, of which GAPDH and EXO1 were upregulated in LUAD, while FYN, CFTR, and KLF4 were downregulated ( Figure 5E ).

To facilitate the clinical application, the current study established a multivariate Cox regression model of key module genes, and the risk score of each LUAD case was computed. All cases were classified as high- and low-risk groups with the median value ( Figure 6A ). More death cases were observed in the high-risk group. All key module genes had >0.7 AUCs, demonstrating their excellent efficacy in inferring LUAD prognosis ( Figure 6B ). Moreover, a key module gene-based nomogram was established for inferring LUAD cases’ 5- and 8-year OS outcomes ( Figure 6C ). Calibration curves demonstrated that the nomogram performed well in contrast with the ideal model ( Figure 6D ). We also analyzed the differences in clinicopathological features and risk score between C1 and C2 subtypes. In both TCGA-LUAD and GSE72094 datasets, the proportion of stage I–II in C2 was higher than that in C1, indicating better prognostic outcomes for patients in the C2 subtype ( Figures 6E, F ). However, no notable difference in risk score was observed between C1 and C2 subtypes.

Among key module genes, we further focused on EXO1. As illustrated in Figure 7A , EXO1 was notably linked to multiple immune cell types. In addition, it prominently correlated to most MHC molecules ( Figure 7B ) and immune checkpoints such as PD-L1 ( Figure 7C ). The above evidence indicated that EXO1 correlated to antitumor immunity in LUAD. To characterize the function of EXO1 in LUAD cells, we first measured EXO1 mRNA level in human lung normal epithelial cell line (BEAS-2B) and LUAD cell lines (NCI-H1975, NCI-H1395, A549, and PC9). Among LUAD cell lines, EXO1 exhibited the highest level in NCI-H1975 and the lowest in A549 ( Figure 7D ). Thus, A549 cells were utilized to overexpress EXO1, with NCI-H1975 cells for its knockdown ( Figures 7E, F ). Overexpressed EXO1 dramatically reinforced the growth of A549 cells, with the opposite findings in NCI-H1975 cells with knockout EXO1 ( Figures 7G, H ), indicating that EXO1 might be essential for LUAD growth.

The current study further observed whether EXO1 impacted the migration and invasion of LUAD cells. Overexpressed EXO1 notably reinforced the number of migratory A549 cells, while EXO1 knockdown produced the opposite findings ( Figures 8A–D ). In addition, the number of invasive A549 cells was observably heightened by EXO1 upregulation, with opposite findings in NCI-H1975 cells with knockout EXO1 ( Figures 8E–H ). Altogether, EXO1 may motivate migratory and invasive features of LUAD cells.

Next, we observed whether EXO1 affected PD-L1 expression in LUAD. The results showed that overexpressed EXO1 enhanced PD-L1 expression in A549 and NCI-H1975 cells according to RT-qPCR and Western blotting ( Figures 9A–E ). Oppositely, EXO1 knockdown lowered PD-L1 expression in two LUAD cells ( Figures 9F–J ). In addition, sPD-L1 expression was examined. As a result, EXO1 overexpression elevated sPD-L1 levels in LUAD cells ( Figures 9K, L ), with opposite results when EXO1 was knocked out ( Figures 9M, N ).