Identification of malignant cells from scRNA-seq dataFollowing rigorous quality control measures and data filtering, a total of 39,639 single cells were obtained from the scRNA-seq data. Utilizing canonical cell markers, various cell types were distinguished, encompassing immune cells (myeloids, NK cells, mast cells, B cells, and T cells), fibroblasts, endothelial cells, and epithelial cells (Fig. 1A,B). The UMAP plot illustrates that the distribution of most cell types is similar between EGFR-MU and EGFR-WT tumors (Fig. 1C). For cancer cells, 14 subclusters were further identified (Fig. 1D). EGFR-MU and EGFR-WT cancer cells were clustered according to their transcriptional patterns (Fig. 1E). The cancer cell transcriptional profiles seemed more distinguishable between individual EGFR-WT tumors than between individual EGFR-MU tumors (Fig. 1F), suggesting that EGFR-MU cancer cells are more transcriptionally homogenous. Then, to select malignant cells with high confidence for cell type-specific functional analysis, we obtained somatic large-scale chromosomal copy number variations (CNVs) and calculated CNV scores for each cancer cell cluster. Ciliated cells, endothelial cells, and fibroblasts were used as references (Fig. 1G,H). As expected, most cancer cells exhibited higher numbers of CNVs relative to the reference cells. To eliminate potential normal epithelial cells, we excluded three clusters—specifically clusters 6, 9, and 11—whose median CNV scores were below the average (Fig. 1H). The remaining 11 clusters were subsequently utilized for cancer cell-specific analysis. Alternatively, we considered all 14 clusters and excluded cells with CNV scores in the lowest 5%, as these were likely normal cells. Subsequent analyses using both methods to identify cancer cells yielded highly similar results. For conciseness, the results for the 11 clusters of cancer cells are presented in the main text, whereas the results derived from the alternative method are displayed in Supplemental Figs. 3–5.Figure 1Identification and clustering of single cells in LUAD samples. (A,C): UMAP plot of the major cell types combined or grouped by EGFR-mutation status. (B): Heatmap of the mean expression of canonical marker genes for ten major cell types. (D–F): UMAP plot of cancer cells grouped by cluster, EGFR-mutation status, and sample, with each group assigned a unique color. (G): Heatmap of copy number variation (CNV) for each cell. (H): CNV scores for each cancer cell cluster and referenced clusters.
EGFR mutation-related genes and modules as well as their functionsWe performed hdWGCNA analysis to explore gene modules that are functionally related to EGFR-mutation status. After optimization, the software threshold power was set at 26 to construct a coexpression network (Supplemental Fig. 1). Four modules were identified and each was assigned a unique color as an identifier. Then, the harmonized module eigengenes (hMEs) were calculated to determine the unique distribution area of each module (Fig. 2A). To compare the module expression patterns between EGFR-MU and EGFR-WT cancer cells, we performed differential module eigengene (DME) analysis and found that the yellow and brown modules were enriched in the EGFR-MU group, while the turquoise and blue modules were enriched in the EGFR-WT group (Fig. 2B). Next, the top 25 hub genes for each module were used to calculate signature scores, as depicted in the dimensionality reduction plot (Fig. 2C and Supplemental Fig. 2). No significant positive correlation was observed between the different modules, and gene network analysis suggested that these modules were linked to distinct biological functions (Fig. 2C,D). However, given the small sample size (n = 5 vs. 4; EGFR-MU and EGFR-WT) of heterogeneous lung cancer patients, it is possible that there are other confounding variables in this analysis other than EGFR mutation status.Figure 2Identification of modules associated with EGFR mutations by hdWGCNA. (A): Visualization of the different coexpression modules by network analysis. (B): UMAP plot of each coexpression module colored in each module’s uniquely assigned color. (C): Differential module eigengene (DME) analysis for different EGFR-mutation groups. (D): Visualization of the correlation between each module.Subsequently, enrichment analysis (KEGG and HALLMARK) was performed to study the potential biological functions of each module (Fig. 3A). The expression levels of the top 20 hub genes for each module in EGFR-MU and EGFR-WT cancer cells are indicated in Fig. 3B and Supplement Fig. 2B. As mentioned before, the yellow and brown modules were more engaged in EGFR-MU cancer cells. The results reveal that the yellow module regulates MAPK or KRAS signaling and the epithelial–mesenchymal transition (EMT) process, both of which play critical roles in cancer cell proliferation, migration, and metastasis. Meanwhile, the brown module was found to be associated with the regulation of various immune-related pathways such as inflammatory response, cytokine-cytokine signaling, and chemokine signaling (Fig. 3A).Figure 3Hub genes and enrichment analysis of each module. (A): HALLMARK and KEGG enrichment analysis of modules. (B): Comparative analyses of the normalized expression levels of the top 20 hub genes from each module between EGFR-MU and EGFR-WT cancer cells. (C): Gene module scores in EGFR-MU or EGFR-WT tumors in the TCGA-LUAD datasets. A comparative analysis of the normalized expression levels of the top 20 hub genes.The turquoise module was found to be engaged in immunoregulatory processes, such as antigen-processing and presentation, as well as response to IFN-γ. Notably, HLA-B, a hub gene in the turquoise module, seemed to be downregulated in EGFR-MU tumors (Fig. 3B). Although the turquoise module was least differential with respect to EGFR mutation status (Supplemental Fig. 3), we further analyzed the expression of other MHC-I molecules, including B2M and HLA-A. Our findings indicate that their expression is generally reduced in EGFR-MU tumors relative to EGFR-WT tumors (Fig. 5A). Therefore, we hypothesize that EGFR-MU cancer cells may evade CD8+ cytotoxic T-cell-mediated attack by inhibiting MHC-I-dependent antigen presentation. The turquoise module is also involved in the regulation of various metabolic pathways, such as glutathione metabolism, phenylalanine metabolism, and hypoxia (Fig. 3A).The TNF-α/NF-κB pathway was enriched in the turquoise and blue modules (Fig. 3A). Given that TNF-α is mainly produced by intratumoral immune cells, these data suggest that EGFR-WT tumors may possess a more inflamed or immunologically “hot” tumor milieu. However, immune cells were not excluded from EGFR-MU tumors, as shown in Fig. 1C. Furthermore, distinct inflammation-related pathways were enriched in the brown module (Fig. 3A), implying that immune landscapes differ significantly between EGFR-MU and EGFR-WT tumors.Finally, although the hdWGCNA analysis revealed gene modules potentially responsible for the suppressive TME, the sample sizes were small. Therefore, we used the TCGA-LUAD bulk RNA-seq dataset to calculate module scores for individual samples based on the expression of hub genes specific to each module. The differential enrichment of these modules based on EGFR-mutation status aligns with previous findings and strengthens the validity of the scRNA-seq analysis results, as depicted in Fig. 3C.Suppressed T-cell responses and altered MHC-I/II expression in EGFR-MU tumorsTo study how EGFR mutations affect the T-cell immune response, we initially characterized major T-cell subpopulations based on canonical marker genes (Fig. 4A,B). While not achieving statistical significance, EGFR-MU tumors numerically exhibited a reduced proportion of CD8+ T cells—including that of effector CD8+ T cells (clusters 1 and 9) and exhausted CD8+ T cells (cluster 5)—relative to EGFR-WT tumors (Fig. 4C). Conversely, two Treg subtypes, Treg1 and Treg2 (clusters 6 and 8), also were numerically more prevalent in EGFR-MU tumors (Fig. 4C) while the difference was statistically not significant. These results are consistent with prior research11,14,15. We then analyzed cell–cell communication using the CellChat software. In brief, cancer cells and T-cell subpopulations were grouped first by EGFR-mutation status. Then, pairwise cell communication was analyzed, and ligand/receptor interactions as well as their strength were compared between EGFR-MU and EGFR-WT groups. The results indicated that although EGFR-WT tumors, in general, had more cell–cell communication, the intensity of this communication was lower than that in EGFR-MU tumors (Supplemental Fig. 6A). Moreover, there is a notable disparity in the strength of particular cell–cell communication pathways, such as MHC-I/II and CXCL, between EGFR-MU and EGFR-WT groups, which will be elaborated upon in the following discussion. Then, specific molecular interactions between cancer cells and T cells were examined. Interactions that were enhanced or suppressed in EGFR-MU tumors relative to EGFR-WT tumors are shown in Fig. 4D–E and Supplemental Fig. 6B. Analysis of the interaction between cancer cells and CD8+ effector T cells (CD8eff_1 and CD8eff_2) indicated that the MHC-I-CD8 interaction was decreased in the EGFR-MU group (Fig. 4D), consistent with the reduced expression of MHC-I genes in EGFR-MU cancer cells (Fig. 5A). These results further support that decreased MHC-I gene expression is essential for the ability of EGFR-MU cancer cells to evade immune attacks from CD8+ cytotoxic T cells.Figure 4Identification of T-cell subtypes and cell–cell communication analysis. (A): UMAP plot of all T cells colored by subtype and EGFR-mutation status. (B): Heatmap of the mean expression of canonical marker genes for different T-cell subtypes. (C): Proportions of each subcluster in the EGFR-MU and EGFR-WT groups. (D,E): Significantly altered ligand-receptor interactions between cancer cells and effector CD8+ T cells. (D) or Tregs (E) in EGFR-MU or EGFR-WT tumors.Figure 5Examining HLA-DR positive cancer cells by mIHC. (A,B): Violin plots of HLA-I and HLA-II expression in EGFR-MU and EGFR-WT tumors. (C): Representative mIHC images of EGFR-MU (EGFR 19del and EGFR L858R) and EGFR-WT tumors. (D): Clinical features of the LUAD cohort for mIHC analysis. (E): Percentages of HLA-DR+ cancer cells (CD45-PanCK+ HLA-DR+) among the total cancer cells (CD45-PanCK+) in EGFR-WT and EGFR-MU tumors.Next, we analyzed the interaction between cancer cells and Tregs, and found that the interaction between MHC class II molecules (MHC-II) derived from cancer cells and CD4 derived from Tregs was significantly enhanced in EGFR-MU tumors (Fig. 4E). Consistent with this finding, the expression of MHC-II genes, including HLA-DR, HLA-DP, and HLA-DQ, was significantly upregulated in EGFR-MU cancer cells (Fig. 5B and Supplemental Fig. 7A). To verify this finding, we conducted multiplex immunohistochemistry (mIHC) analysis on 48 surgically resected tumor samples from patients with LUAD, including 22 EGFR-WT and 26 EGFR-MU samples (Supplemental Table 1). The results demonstrated that HLA-DR expression was significantly elevated in EGFR-MU cancer cells compared with that in EGFR-WT cancer cells in the LUAD cohort (Fig. 5C–E). Notably, the expression of costimulatory molecules, including CD80, CD86, and ICAM1, was significantly lower in cancer cells than in myeloid cells in EGFR-MU tumors (Supplemental Fig. 7B). Prior research indicates that insufficient CD28-mediated costimulatory signaling can lead to the differentiation of naïve T cells into Tregs16,17,18. Therefore, we speculate that EGFR-MU cancer cells might function as nonclassical MHC-II antigen-presenting cells and thus skew the differentiation of tumor-infiltrating CD4+ T cells towards Tregs.Previous studies have reported that airway type II (AT-II) epithelial cells express MHC-II and regulate local immune responses in lung tissues19,20,21. To assess whether EGFR-MU cancer cells express MHC-II because of their AT-II cell origin, we analyzed the expression of several surfactant protein genes that indicate the AT-II cell lineage (https://panglaodb.se), including SP-A (SFTPA1, SFTPA2), SP-B (SFTPB), and SP-C (SFTPC)22. In comparison with EGFR-WT cancer cells, EGFR-MU cancer cells expressed high levels of these genes, which positively correlated with HLA-DRB1 expression (Supplemental Fig. 7C). The positive correlation between the expression of surfactant protein and HLA-DRB1 was validated in the TCGA-LUAD dataset (Supplemental Fig. 7D).Roles of chemokines/chemokine receptors in recruiting Tregs into tumorsCell–cell communication analysis indicated that several chemokine/chemokine receptor interactions, such as CXCL1-CXCR2, CXCL3-CXCR2, and CXCL16-CXCR6 (Fig. 4E), between cancer cells and Tregs were elevated in EGFR-MU tumors, consistent with the increased expression of these chemokines in EGFR-MU cancer cells (Supplemental Fig. 8A). To further explore the roles of chemokine/chemokine receptors in Treg recruitment as well as their cellular sources, we analyzed the single-cell expression profile of several chemokine receptors that have been reported to be expressed by Tregs and to mediate Treg recruitment, including CCR423,24, CCR725,26, CCR827,28, CCR1029,30, and CXCR431,32 as well as CXCR2 and CXCR6 identified in our analysis. As shown in Fig. 6A and Supplemental Fig. 8B, CCR4 was expressed at higher levels by Treg2 cells in EGFR-MU tumors, although CCR4 expression was not restricted to Tregs. In contrast, CCR8 expression was more restricted to Tregs, and its expression by Treg2 cells in EGFR-WT tumors was more abundant. CCR10 and CXCR2 were generally expressed at low levels by Tregs, while CXCR4 was highly and universally expressed by Tregs, with higher expression in Treg1 cells than in Treg2 cells. CXCR6 was highly expressed by Treg2 cells, and its levels were comparable between EGFR-MU and EGFR-WT tumors. These results reveal a complex pattern of chemokine receptor expression by Tregs in vivo, and suggest that multiple chemokine receptors might be involved in the increased infiltration of Tregs in EGFR-MU tumors, particularly of Treg2 cells that express CCR4.Figure 6Role of chemokines in recruiting Tregs in EGFR-MU LUAD. (A): Violin plots representing chemokine receptor expression in the Treg subtypes. (B): Volcano plot of differentially expressed genes between Treg1 and Treg2 cells; Delta Percent refers to the percentage of (Treg1 vs. Treg2) expression of a specific gene. (C): Developmental trajectory of Treg1, Treg2, and naïve CD4+ T cells. (D): Violin plots representing chemokine ligand expression in all cell types in EGFR-MU and EGFR-WT tumors. (E, F): Association between CCL17/CCL22 and FOXP3 in EGFR-MU LUAD in the TCGA database. (G): Boxplot showing the differential expression of FOXP3 in the TCGA-LUAD samples grouped by median expression of CCL17 and CCL22. (H–I): CCL17/CCL22 expression in EGFR-MU cells (PC9, NCIH3255, NCIH1975, NCIH1650, HCC827, and HCC4006) compared with that in EGFR-WT cells (NCIH69, NCIH522, NCIH520, NCIH460, NCIH292, NCIH1299, A549, and A427) in the CCLE database.To examine whether the Treg1 and Treg2 subtypes were functionally different in terms of immune suppression, we analyzed differentially expressed genes (DEGs) between Treg1 and Treg2 cells. Some genes known to be involved in Treg-mediated immune suppression, including FOXP3, CTLA4, TNFRSF9, and TIGIT, were more highly expressed by Treg2 cells than by Treg1 cells, suggesting that Treg2 cells might be functionally more immunosuppressive (Fig. 6B and Supplemental Fig. 8C). Trajectory analysis of CD4+ naïve T, Treg1, and Treg2 cells suggested that Treg2 cells were terminally differentiated and derived from naïve CD4+ T cells, whereas Treg1 cells were in an intermediate state (Fig. 6C). These findings imply that Treg1 cells might be differentially unstable or exhibit a transitional differentiation status.Analysis of TF activity revealed that certain TFs involved in cytokine-mediated signal transduction, including STAT1, STAT2, and STAT5, exhibited increased activity in Treg2 cells (Supplemental Fig. 9A). In addition, NFATC2 and SMAD3, which reportedly enhance FOXP3 expression, were also activated in Treg2 cells33. Pathway enrichment analysis revealed that the cytokine/cytokine receptor signaling pathway was significantly enriched in Treg2 cells relative to Treg1 cells (Supplemental Fig. 9B–C). Interestingly, Treg2 cells in EGFR-MU tumors seemed to express less POXP3 and CTLA4 than their counterparts in EGFR-WT tumors (Supplemental Fig. 8C), suggesting that Treg2 cells in EGFR-MU tumors might be relatively less active; this may be attributable to attenuated feedback-triggered activation of Tregs resulting from weaker CD8+ T-cell reactions.To investigate the chemokine ligands involved in Treg recruitment and their sources, we comprehensively analyzed the expression of major ligands for the critical chemokine receptors (Fig. 6A). Interestingly, CCL1734, which binds to CCR4, was predominantly expressed by myeloid cells and cancer cells in EGFR-MU tumors. CCL22, another CCR4 ligand34, was expressed in a similar pattern as CCL17 (Fig. 6D). Analysis of the TCGA datasets verified that the CCL17 and CCL22 levels positively correlated with FOXP3 expression in EGFR-MU LUAD tumors (Fig. 6E–F). TCGA-LUAD tumors with high levels of both CCL17 and CCL22 (CCL17–CCL22high) expressed significantly more FOXP3 than CCL17–CCL22low LUAD tumors (Fig. 6G). In addition, analysis of the CCLE database (https://sites.broadinstitute.org/ccle/) verified that CCL17 and CCL22 expression levels were higher in EGFR-MU lung cancer cells than in EGFR-WT cells (Fig. 6H–I). These results suggest that CCL17 and CCL22 produced by EGFR-MU cancer cells, at least in part, play a role in Treg recruitment.Ligands for CCR7, CCR8, and CCR10 were weakly expressed by various cell types, except for CCL18, which was mainly expressed by myeloid cells with similar expression levels observed in EGFR-MU and EGFR-WT tumors (Fig. 6D). Ligands for CXCR2 were mainly expressed by myeloid cells, neutrophils, and cancer cells (Fig. 6D). Considering the weak expression of CXCR2 by Tregs (Fig. 6A), the roles of CXCR2 ligands in Treg recruitment are unclear. CXCL16, the ligand for CXCR6, was mainly expressed by myeloid cells, neutrophils, and mast cells (Fig. 6D). CXCR6 levels were higher in Treg2 than in Treg1 cells (Fig. 6A); however, CXCL16 was not correlated with FOXP3 expression in EGFR-MU LUAD in the TCGA datasets (Supplemental Fig. 8D). Collectively, our results reveal that multiple chemokine/chemokine receptors might mediate Treg recruitment—myeloid and cancer cell-derived CCR4 ligands, CCL17 and CCL22 might be critical in Treg recruitment in EGFR-MU tumors, whereas CCR8 and its ligand CCL18 expressed by myeloid cells might play a more important role in Treg recruitment in EGFR-WT tumors.CAF and tumor endothelial cell (TEC) subtypes are associated with EGFR-mutation statusPrevious studies have revealed that a subpopulation of cancer-associated fibroblasts (CAFs) with antigen-presenting capability (apCAFs) can induce the expansion of Tregs in pancreatic cancer35. We analyzed major CAF subtypes—including myofibroblasts (mFibs), myoCAFs (mCAFs), inflammatory CAFs (iCAFs), apCAFs, and proliferating CAFs (pfCAFs)—based on the expression of classical marker genes (Fig. 7A–C). The findings indicated that the proportion of apCAFs (MHC-II+) was increased while that of iCAFs was reduced in EGFR-MU tumors, relative to EGFR-WT tumors (Fig. 7D). We also analyzed DEGs among CAFs (Fig. 7E) and performed an enrichment analysis (Fig. 7F). The results demonstrated that different CAF subtypes carried out distinct functions—specifically, apCAFs actively participated in antigen processing and presentation, whereas in iCAFs, cytokine- and chemokine-mediated effects, such as TNF/NF-κB and IL-17 signaling, were activated, reflecting that their functions were closely associated with the establishment of an inflammatory environment.Figure 7Subtypes, functionality, and developmental trajectory of CAFs. (A, B): UMAP plot of all CAF subtypes and EGFR-mutation status. (C): Heatmap of the mean expression of canonical marker genes for different CAF subtypes. (D): Proportion of each CAF subtype in the EGFR-MU and EGFR-WT groups. (E): Differentially expressed genes in each CAF subtype relative to the other CAF subtypes. (F): KEGG and GO enrichment analysis for each CAF subtype. (G): Semisupervised trajectory of CAFs, colored by subtype and pseudotime. (H): Heatmap showing inferred transcription factor activity in different CAF subtypes. (I,J): Boxplot representing the fractions of apCAFs or iCAFs in EGFR-MU and EGFR-WT tumors in the TCGA-LUAD datasets.Trajectory analysis revealed that apCAFs and iCAFs were derived from mFibs, following distinct differentiation routes (Fig. 7G) with the preferential generation of apCAFs or iCAFs in EGFR-MU or EGFR-WT tumors, respectively. Next, we found that compared with other CAF subtypes, iCAFs exhibited a very unique TF activity profile featuring enhanced activity of TFs engaged in NF-κB and cytokine signaling pathways, such as RELA/RELB, NFKB1, and STAT3, suggesting the critical roles of inflammatory cytokines in iCAF differentiation and function (Fig. 7H). In apCAFs, several TFs related to cellular stemness, such as FLI, SIX2, and SOX10, as well as RFX5, which has been reported to promote MHC-II expression, were highly active (Fig. 7H). To overcome the limitation of the small scRNA-seq sample size, we examined whether EGFR mutations affected CAF differentiation by quantifying apCAFs and iCAFs in the TCGA-LUAD dataset with cell type-specific matrixes derived from the scRNA-seq data and CIBERSORT algorithm. Consistently, the results indicated that apCAFs were enriched in the EGFR-MU TCGA-LUAD samples, whereas iCAFs were more abundant in the EGFR-WT samples (Fig. 7I–J). To validate these findings, we examined CAF subsets in the LUAD cohort with mIHC analysis. The results confirmed enrichment of apCAFs and iCAFs in EGFR-MU and EGFR-WT tumors, respectively (Fig. 8A–C). Furthermore, we utilized normalized deviation to evaluate the composition of apCAFs or iCAFs in all CAFs, and found that the proportion of apCAFs was higher than that of iCAFs in EGFR-MU tumors compared with EGFR-WT tumors (Fig. 8D).Figure 8Analysis of apCAFs and iCAFs by mIHC. (A): Representative mIHC images of EGFR-MU (EGFR 19del and EGFR L858R) and EGFR-WT tumors. (B): Fractions of apCAFs (Pan-CK-CD45-αSMA+ HLA-DR+ IL-6−) among total CAFs (Pan-CK-CD45-αSMA+) in EGFR-MU and EGFR-WT tumors. (C): Fractions of iCAFs (Pan-CK-CD45-αSMA+ HLA-DR-IL-6+) among total CAFs in EGFR-MU and EGFR-WT tumors. D: Normalized deviation between apCAFs and iCAFs formulated by [(% apCAFs—% iCAFs)/% CAFs] in EGFR-MU and EGFR-WT tumors.Tumor endothelial cells (TECs) are critical tumor stromal cells and can serve as APCs in response to inflammatory stimuli. We identified four vascular TEC clusters (Vascular_end_1–4) and one lymphatic TEC cluster (Fig. 9A–C) based on canonical marker genes. Cluster 3 TECs (vascular_end_3) were more abundant in EGFR-MU tumors than in EGFR-WT tumors (Fig. 9D). Furthermore, both DEG and pathway enrichment analyses demonstrated that their function was related to MHC-II antigen presentation (Fig. 9E–F). Notably, the expression of CD86 was markedly lower in TECs and CAFs than in myeloid cells (Supplemental Fig. 7B), suggesting that, similar to apCAFs, vascular_end_3 TECs lack sufficient costimulatory molecules and might be involved in immune suppression instead of immune activation.Figure 9Subtypes, functionality, and developmental trajectory of ECs. (A–B): UMAP plot of all EC subtypes and EGFR-mutation status. (C): Heatmap of the mean expression of canonical marker genes for different EC subtypes. (D): Proportion of EC subtypes in the EGFR-MU and EGFR-WT groups. (E): Differentially expressed genes in each EC subtype relative to the other EC subtypes. (F): KEGG and GO enrichment analysis for each EC subtype.We also evaluated the fractions of various immune cell types in each TCGA-LUAD sample and found that CD8+ T cells and NK cells were decreased in EGFR-MU samples relative to EGFR-WT samples, whereas Tregs and myeloid cells showed an inverse tendency (Fig. 10). These findings were again consistent with the scRNA-seq analysis results.Figure 10Analyzing gene modules and immune infiltration with TCGA bulk RNA-seq data. Heatmap showing infiltration of multiple immune cell types in EGFR-MU and EGFR-WT tumors in the TCGA-LUAD datasets.
Integrative transcriptome analysis reveals the molecular events underlying impaired T-cell responses in EGFR-mutant lung cancer
