Four cohorts were included for evaluating the clinicopathologic, genomic, transcriptomic, and immune profiles of SPA: ① LUAD-TCGA: 199 untreated LUAD patients with genomic data (197 with transcriptomic data, 148 with proteomic data) and confirmed histologic subtypes from The Cancer Genome Atlas (TCGA) repository (12); ② LUAD-MSKCC: 604 untreated LUAD patients with next-generation sequencing data and definite histologic subtypes from the MSKCC cohort (6); ③ LUAD-Singapore: 173 untreated LUAD patients with detailed histologic subtypes (172 with genomic data, 141 with transcriptomic data) from Singapore (13); ④ LUAD-NCC: 437 untreated LUAD patients (103 with genomic and proteomic data, 51 with transcriptomic data) from our center (14). LUAD histologic subgroups were classified into 2 categories: SPA versus Non-solid predominant adenocarcinoma (Non-SPA).

In addition, a cohort of 48 LUAD patients who received neoadjuvant immunotherapy in our center was used to validate the better suitability of immunotherapy for SPA. A flow diagram was drawn to illustrate the study design ( Figure 1 ). The baseline information of the main data sets included is summarized in Table S1 . This study was approved by the Ethics Committee of National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College (Approval No. 2016YJC-01). Written informed consent was obtained from all participants.

Whole exome sequencing and RNA sequencing were performed in LUAD-TCGA, LUAD-Singapore and LUAD-NCC (12–14). MSK-IMPACT next-generation sequencing was performed in LUAD-MSKCC (6). The protein quantitation was based on reverse phase protein array (RPPA) in LUAD-TCGA and mass spectrometry in LUAD-NCC (12, 14). The detailed methods and experimental procedures regarding DNA, RNA, and protein extraction from tumors, library preparation, sequencing, quality control, and subsequent data processing have been previously reported (6, 12–14). The updated clinical, genomic, transcriptomic and proteomic data of LUAD-TCGA was retrieved from the TCGA data portal (https://gdc.cancer.gov/) using the R package “TCGAbiolinks” (15) and the corresponding histologic subtypes of patients were obtained from the supplementary material of the study conducted by TCGA Research Network (12). The clinical and genomic data of LUAD-MSKCC was obtained from the cBioPortal database (https://www.cbioportal.org/study/summary?id=luad_mskcc_2020). The clinical, genomic, and transcriptomic data of LUAD-Singapore was derived from the Singapore Oncology Data Portal (https://src.gisapps.org/OncoSG/) and the histologic subtype of cases were obtained from the supplementary material of the previously published study (13). The genomic, transcriptomic and proteomic data of LUAD-NCC was collected from the supplementary material of our previously published research (14). The genomic and transcriptomic data of LUAD-Singapore and LUAD-NCC were integrated into a combined cohort (LUAD-Asia) due to similar ethnic backgrounds. For rare mutational events, LUAD-TCGA, LUAD-MSKCC, and LUAD-Asia were combined into a larger cohort (LUAD-Combined).

The mutation status of each gene was inferred from the mutation MAF files by the maftools package (16). OncoPrints constructed by the ComplexHeatmap R package (17) were used to depict the mutation landscape. TMB was defined as the total nonsynonymous somatic mutation counts in coding regions (megabases). Neoantigen load, T cell receptor diversity, intratumor heterogeneity scores and homologous recombination deficiency scores of patients in LUAD-TCGA were retrieved from the supplementary material of a previous research (18). Deleterious NOTCH mutations were determined by PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/), which could classify a given mutation as damaging or benign (19, 20). The deconstructSigs R package (21) with default parameters was used to derive COSMIC mutational signatures (22), which account for the observed mutational profile in each patient. Therapeutically targetable alterations were annotated using OncoKB (https://www.oncokb.org/) (23), which collects information from the published research and defines drug actionability according to the clinical evidence.

The gene expression levels were measured as fragments per kilobase per million mapped reads and were then log2 transformed. For combination and comparison of expression data from different datasets, the expression level of each gene was further normalized by a z-score with mean value = 0 and standard deviation = 1. Heatmap was used to depict the gene expression levels by the pheatmap R package. Protein quantitation was evaluated by a label-free quantification algorithm (iBAQ) in LUAD-NCC and values of iBAQ were log2 transformed if necessary. Protein quantitation in LUAD-TCGA was measured by RPPA. The TCGA molecular subtypes were determined using the published LUAD 506-gene nearest centroid classifier (24). The TCGA immune subtypes and the Karolinska proteomic subtypes were obtained from previously published studies (18, 25).

Unbiased gene set enrichment analysis (GSEA) was performed to identify pathways enriched in SPA and Non-SPA using the javaGSEA Desktop Application (26). HALLMARK gene sets from the MSigDB database (26) were selected for GSEA. The threshold was set at false discovery rate (FDR) <0.1. The enrichment score in GSEA was calculated by first ranking the genes from the most to least significant with respect to the two phenotypes (SPA versus Non-SPA), the entire ranked list was then used to assess how the genes of each gene set were distributed across the ranked list. Functional gene expression signatures manually curated (18, 27–32) were used to investigate the correlation between pathological subtypes and other relevant biological processes ( Table S2 ). For each signature, we performed gene set variation analysis and assigned a signature score to each patient using the GSVA R package (33). The 18-gene T cell-inflamed gene expression profile (GEP) score was calculated as weighted sum of the normalized expression values of the genes (27) and the weightings for each gene in the signature ( Table S3 ) were obtained from the published patent (34).

The xCell tool (https://xcell.ucsf.edu) was used to analyze the relative levels of infiltration of 64 immune and non-immune cell types in the tumor immune microenvironment from the RNA sequencing data (35). This method combined the advantages of gene enrichment analysis and deconvolution approaches to evaluate the enrichment of immune cells. The spatial fractions of tumor regions with infiltrating immune cells estimated by analysis of the mapped TCGA digitized hematoxylin-eosin (H&E)-stained slides were obtained from a published study (36).

A total of 103 LUAD patients from LUAD-NCC were assessed by immunohistochemistry (IHC) staining with PD-L1 (clone SP263; Roche Ventana), CD4 (ZA-0519; Zsbio Tech), CD8 (ZA-0508; Zsbio Tech), CD19 (ZM-0038; Zsbio Tech), CD68 (ZM-0060; Zsbio Tech) and CD163 (ZM-0428; Zsbio Tech) antibodies. Formalin-fixed, paraffin-embedded LUAD specimens were cut into 4 mm slides for IHC staining. Slides were stained using an automated Leica Bond staining system according to the manufacturer’s protocol. PD-L1 expression scores were defined as the percentage of tumor cells with membranous staining. Staining status was stratified into three subgroups: negative (PD-L1 <1%), intermediate (PD-L1 1-49%), and high (PD-L1 ≧50%). For CD4 and CD8, the proportion of positive cells was assessed as low density (≦25%), intermediate density (25-50%), and high density (> 50%). For CD19, CD68 and CD163, the percentage of stained cells was classified as follows: low density (≦10%), intermediate density (10-20%), and high density (> 20%). All slides were evaluated by two pathologists and any disagreement was resolved by consensus using a multi-head microscope.

Wilcoxon rank-sum test (or Kruskal-Wallis test) and χ2 test (or fisher’s exact test) were used to assess differences between continuous and categorical variables, respectively. Survival analysis was conducted using the Kaplan–Meier method and log-rank test. Hazard ratios (HRs) were determined by Cox regression analyses. Forest plots showing HRs and confidence intervals were drawn. All statistical analyses were performed with the STATA 16.0 and R 3.6.3. All reported P values were two-sided, and P <0.05 was considered statistically significant unless otherwise specified. For multiple testing, P-values were adjusted using Benjamini–Hochberg FDR correction and FDR q value < 0.10 was considered significant.

The differences in clinicopathologic characteristics between SPA and Non-SPA are detailed in Table S4 . SPA had distinct clinicopathologic characteristics as compared with Non-SPA. SPA was more common in patients with a history of smoking and tended to have a higher rate of lymph node metastasis, a higher pathological stage and a lower degree of differentiation. Next, the survival of SPA and Non-SPA was compared. SPA showed significantly shorter recurrence-free survival (P<0.001; LUAD-MSKCC) and overall survival (P<0.001; LUAD-NCC) than did Non-SPA ( Figures 2A, B ). Multivariate cox regression analyses confirmed that SPA was an unfavorable prognostic factor for both recurrence-free survival (HR: 1.74, P=0.029; Figure 2D ) and overall survival (HR: 1.86, P=0.029; Figure 2E ). In LUAD-Singapore, we also observed a worse overall survival in SPA (P=0.026; Figure 2C ). To uncover the underlying mechanisms, we performed a multi-omics analysis.

We first depicted comprehensive genomic landscapes of SPA and Non-SPA in LUAD-Asia, LUAD-MSKCC and LUAD-TCGA ( Figures 3A , S1A, C ). Despite the genetic diversity of cancer genome between races, SPA and Non-SPA showed distinct genomic profiles across the three cohorts. Consistently higher levels of TMB was observed in SPA ( Figure 3B ) although the median TMB differed markedly between races ( Figure 3D ). Similarly, SPA had higher neoantigen load than Non-SPA ( Figure 3E ), indicating potential suitability for immunotherapy.

Next, we compared the mutation rates of the curated genes (19, 37–42) which were associated with response to immunotherapy between SPA and Non-SPA. Three cohorts were merged into a cohort (LUAD-Combined) due to the low mutation rates of these genes. The results showed that SPA had significantly higher frequency of SMRCA4, ARID1A, ARID2, POLE, PTPRD, and EPHA7 mutations ( Figure 3G ). In terms of co-occurring mutations (12, 38, 43–45), NF1/TP53 co-mutation, KRAS/TP53 co-mutation, and co-mutations in DNA damage response (DDR) pathways were enriched in SPA. A recent study demonstrated that TRAF2 loss and CCND1 amplification were associated with response and resistance to immunotherapy, respectively (46). Here, we found that TRAF2 loss was significantly enriched in SPA while no difference was found in CCND1 amplification between SPA and Non-SPA ( Figure 3F ).

We investigated the alteration frequencies of known driver mutations and those amenable to targeted therapies. Even though SPA had higher number of driver mutations as compared with Non-SPA (mean: 2.27 versus 1.91, P=0.002; Figure 4A ), lower frequencies of therapeutically targetable driver alterations were observed in SPA across the three cohorts ( Figure 4B ). Similar to TMB, higher number of driver mutations was seen in LUAD-TCGA and LUAD-MSKCC as compared with LUAD-Asia ( Figure 4C ).

Next, we compared the mutation rates of five of the most common driver mutations in LUAD, including EGFR, TP53, KRAS, STK11, and KEAP1, which were reported to be correlated with tumor antigenicity and immunogenicity (43, 45, 47, 48). Two genes (TP53, EGFR) were found to be statistically significantly altered in SPA, as compared with Non-SPA across the three cohorts, with significantly higher frequency of TP53 and lower frequency of EGFR mutation in SPA ( Figure 4D ). Our group previously revealed that LUAD patients with EGFR/TP53 co-mutation had poorer overall survival than those with EGFR mutation alone in a small population (14). Here, with larger cohorts, we confirmed that the co-occurrence of EGFR and TP53 mutation was associated with earlier recurrence (P=0.005; Figure 4E ) and poorer survival (P=0.017; Figure 4F ) in EGFR-mutated patients. Importantly, in EGFR-mutated patients, we found that EGFR/TP53 co-mutation was enriched in SPA (P<0.001; Figure 4G ).

Finally, we evaluated the alteration frequencies of 10 hallmark oncogenic pathways (49) and 30 somatic mutational signatures (22) in SPA and Non-SPA. Six pathways were statistically significantly altered in SPA as compared with Non-SPA: Cell cycle, Hippo, Myc, Notch, p53, and PI3K ( Figures S1F , 3C , S1B, D ). SPA had significantly higher number of pathways altered (NPA) than did Non-SPA (P<0.0001; Figure S1E ). Tumors with three or more NPA were enriched in SPA across the three cohorts ( Figure 4H ). Mutational signature analysis revealed that the frequencies of APOBEC signature (signature 2 and 13) and smoking signature (signature 4) were significantly higher in SPA than in Non-SPA ( Figure 4I ). Meanwhile, a notable increase in the rate of transversion/transition was seen in SPA ( Figure S1G ). In addition, SPA showed higher level of homologous recombination deficiency signature (signature 3) with a marginal statistical difference (FDR q value = 0.10). Consistently, higher homologous recombination deficiency score was observed in SPA ( Figure S1H ).

TTF-1 and Napsin-A are two canonical histologic markers used to demonstrate LUAD differentiation (1). In this study, consistent with the clinicopathologic characteristics of SPA (lower degree of differentiation and poorer prognosis), both the mRNA expression and protein abundance of these two markers showed significantly lower levels in SPA ( Figures S2A, B ), indicating that SPA might be atypical LUAD.

Unbiased GSEA of SPA versus Non-SPA showed that SPA was enriched in hallmark pathways associated with proliferation (E2F targets, G2M checkpoint, mTORC1 signaling, Mitotic spindle, MYC targets V1 and V2), immune response (Inflammatory response, IL-6 JAK STAT3 signaling, IFN-γ response, complement, allograft rejection, and IL-2 STAT5 signaling) and hypoxia ( Figure S2C , Table S5 ), consistent with the recent finding that proliferation and immune axes segregated LUAD transcriptomic subgroups (13). However, Non-SPA was not enriched in any biological processes.

To validate the results of GSEA, several functional gene expression signatures were curated ( Table S2 ) and individual signature scores were calculated for each signature. Using the proliferation-related gene signature (18), we confirmed that SPA was more aggressive with higher proliferation scores ( Figure S2D ), which was further supported by the increased Ki67 expression (both mRNA and protein level) in SPA ( Figure S2E ). In terms of hypoxia, SPA displayed higher hypoxia scores and increased HIF-1α expression ( Figures S2F, G ), indicating a potential hypoxic microenvironment in this subtype.

Since SPA showed a close relationship with immune-related pathways in GSEA, we explored the expression of several clusters of metagenes previously reported to be associated with response to immunotherapy (27, 31, 32), including an effector CD8 T cell signature (CD8 T effector), an immune checkpoint signature, a six-gene IFNγ signature (IFNγ-6), a chemokine signature and antigen presenting machinery (APM) signature. Heatmaps depicting the expression levels of these gene signatures revealed clear differences between SPA and Non-SPA. Specifically, SPA demonstrated higher levels of mRNA expression in these immunity markers ( Figure 5A ). To reinforce the robustness of our findings, scores of each signature were compared between these two subtypes. The results showed significantly higher scores across all these signatures in SPA, indicating preexisting immunity within this subtype ( Figure 5B ). Consistently, we observed significantly higher GEP scores (a pan-cancer predictor of response to immunotherapy (27)) in SPA ( Figure 5C ). Individual immune gene markers were also investigated. Both the LUAD-TCGA and LUAD-Asia cohorts showed significantly higher PD-L1 mRNA and protein expression in SPA than in Non-SPA ( Figures 5D, E ). IHC analysis confirmed that SPA had higher PD-L1 expression in the LUAD-NCC cohort (P<0.001, Figures 5F, G ). In addition, SPA showed greater T cell receptor diversity ( Figure 5H ) and higher CXCL9/CXCL13 expression ( Figure S2H ).

Next, we evaluated the infiltrating immune cells in SPA and Non-SPA. SPA had a distinct tumor immune microenvironment as compared with Non-SPA, with higher infiltration of M1 macrophages, Th1 and Th2 cells, CD8+ T cells and macrophages as well as higher immune scores in SPA ( Figures 6A , S3A ). No difference in infiltration of M2 macrophages, neutrophils and CD4+ T cells was observed between SPA and Non-SPA. The analysis of mapped TCGA digitized H&E-stained images also demonstrated higher proportion of tumor-infiltrating lymphocytes (TILs) and macrophages in SPA ( Figure S3B ). IHC analysis of the 103 LUADs in LUAD-NCC confirmed that SPA had higher density of CD8+ T cells (P=0.007) and macrophages (P=0.023) infiltration ( Figures 6B , S3E ) while no difference in M2 macrophages, CD4+ T cells, or B cells infiltration was found between SPA and Non-SPA. Th17 cell is generally associated with improved overall survival (50). In this study, using the Th17 cell signature ( Table S2 ) (18), we found that SPA had lower Th17 cells infiltration as compared with Non-SPA ( Figure S3C ).

Further, we explored the correlation between pathological subtypes and the four types of tumor immune microenvironment classified based on CD8A and PD-L1 expression (51). Positive CD8A and PD-L1 were defined as expression above the median level. As expected, we observed that SPA had a higher proportion of type I tumor immune microenvironment (PD-L1+/CD8A+) than Non-SPA ( Figure 6C ), indicating adaptive immune resistance in SPA. A recent study revealed that tumors with type I tumor immune microenvironment may have more complex genomic intratumor heterogeneity (52), which was associated with postoperative recurrence and poor survival in lung cancer (53, 54). Therefore, genomic intratumor heterogeneity was compared between SPA and Non-SPA. As expected, SPA was associated with increased intratumor heterogeneity ( Figure 6D ).

β-Catenin protein level was reported to inversely correlate with T cell infiltration (55). Therefore, we compared the abundance of β-Catenin protein between SPA and Non-SPA using the proteomic data in LUAD-TCGA and LUAD-NCC. We noted significantly lower levels of β-Catenin in SPA ( Figure 6E ). Finally, to explore why adjuvant chemotherapy was unsatisfactory in SPA at transcriptome level, we curated two gene expression signatures named chemoresistence signature and chemotherapy response signature ( Table S2 ) (29, 30), which were associated with chemotherapeutic outcome, and calculated signature score for each patient. As expected, SPA showed significantly higher chemoresistence signature scores and lower chemotherapy response signature scores in both LUAD-TCGA and LUAD-Asia ( Figure S3D ).

To validate the aforementioned findings, we assessed the correlation between pathological subtypes (SPA versus Non-SPA) and different transcriptomic and proteomic subtypes previously published. First, we evaluated the TCGA molecular subtypes. We found that SPA was enriched for the “Proximal-inflammatory” subtype while Non-SPA was enriched for the “Terminal respiratory unit” subtype ( Figure 7A ). Next, the TCGA immune subtypes were assessed. The analysis revealed that SPA consisted primarily of the “IFN-γ dominant” subtype ( Figure 7B ), consistent with the result of unbiased GSEA which showed that SPA was enriched for the IFN-γ response pathway ( Figure S2C and Table S5 ). Then, SPA and Non-SPA were classified into the Karolinska proteomic subtypes. Consistent with higher immune infiltration in SPA, SPA showed enrichment of the “Immune-hot” subtype in both LUAD-TCGA and LUAD-NCC ( Figure 7C ). Finally, analysis of the proteomic subtypes recently proposed by our group (NCC proteomic subtypes) showed that SPA was enriched for the “Proliferation and proteasome” subtype ( Figure 7D ).

The aforementioned multi-omics analyses revealed that SPA was more suitable for immunotherapy as compared with Non-SPA. To validate this finding, we investigated the non-small cell lung cancer patients who received neoadjuvant immunotherapy (alone or in combination with chemotherapy or angiogenesis inhibitor) from October 2018 to May 2022 in our center, of which 48 LUAD patients with detailed pathological results were included for further analysis ( Table S6 ). We first compared the pathological regression rates between these two subtypes. As expected, SPA had significantly higher regression rates than Non-SPA ( Figures 7E, F ). Next, we stratified the pathological regression rates into three groups: marked regression (regression rate: 80-100%), moderate regression (regression rate: 40-80%), and mild regression (regression rate: 0-40%). Patients with marked regression were enriched in SPA while patients with mild regression were enriched in Non-SPA ( Figure 7G ). Further, we compared the major pathological response (MPR) between these two subtypes. When the cut-off value was set as ≦10% viable malignant cells, SPA had a higher proportion of MPR (14.3% versus 5.9%) than Non-SPA although the difference was not significant ( Figure 7H ). When using the newly proposed cut-off value of MPR for LUAD (≦65% viable malignant cells) (56), patients with MPR were significantly enriched in SPA ( Figure 7I ).