Tumor microenvironment remodeling after neoadjuvant immunotherapy in non-small cell lung cancer revealed by single-cell RNA sequencing

scRNA-seq analysis of NSCLC during PD1 blockade combined with chemotherapy

We prospectively collected fresh tumor samples from a total of 15 patients with clinical stage IIIA NSCLC for analysis by scRNA-seq (Fig. 1A, Additional file 2: Fig. S1A and Additional file 1: Table S1). For three patients, samples were collected by biopsy before treatment and classified as treatment naïve (TN; n = 3). For the remaining 12 patients, samples were taken from surgical resections after PD-1 antibody combined with chemotherapy treatment. The 12 post-treatment samples were categorized into two groups: MPR (n = 4) and NMPR (non-major pathologic response; n = 8) based on pathologic assessment [5]. The dataset analyzed here also includes bulk RNA-seq from fresh biopsies from 21 independent TN patients (Additional file 1: Table S1).

The fresh tissues were rapidly digested to a single-cell suspension, and all single-cell transcriptomes were generated using commercial BD Rhapsody platform. After quality control and removal of doublets, transcriptomes from 92,330 cells with a median of 1256 genes per cell were used for further analyses. To mitigate batch effects from patients (Additional file 2: Fig. S1B) and allow for joint analysis of malignant and non-malignant cells, we performed canonical correlation analysis (CCA) and aggregated cells from different patient samples. Unsupervised clustering of all cells identified 26 clusters (Fig. 1B and Additional file 2: Fig. S1C), with no significant batch effects observed across different patients, PD-1 antibodies, or response groups (Additional file 2: Fig. S1D-E). Further, the average gene numbers and unique molecular identifiers (UMIs) were comparable between different clusters (Additional file 2: Fig. S1F). We then annotated the 26 clusters into T cells, NK cells, B cells, myeloid cells, neutrophils, plasma cells, plasmacytoid DC (pDC), mast cells, stromal cells (fibroblasts/endothelia), and epithelial cells, according to the expression of corresponding canonical marker genes (Fig. 1B and Additional file 2: Fig. S1G).

To characterize the TME remodeling in response to treatment, we calculated the fraction of different cell types in TN, MPR, and NMPR patients (Fig. 1C, D). We observed that the fraction of T cells, NK cells, and B cells were increased in MPR patients, although we did not get positive P values due to limited sample sizes (Fig. 1D). To further validate this, we performed immunohistochemical (IHC) staining in 10 post-treatment surgical tumor tissues (3 MPR and 7 NMPR, corresponding to scRNA-seq samples) and another 5 treatment-naïve surgical tumor tissues (Additional file 1: Table S1). IHC staining verified the increased abundance of T (CD3+) and B (CD20+) cells in MPR samples, except NK (CD56+) cells (Fig. 1E, F, and Additional file 2: Fig. S2). This was consistent with previous reports that T and B cell expansion was associated with better response to ICB [10, 13]. Myeloid cells were enriched in the TME, but showed no obvious differences among TN, MPR, and NMPR patients (Fig. 1D). However, myeloid cells are known to have diverse and complex functions in the TME [36], which are further explored later in this study. We also identified 5–20% cells as neutrophils in the TME of NSCLC [37], which are usually absent in previous scRNA-seq studies using 10X Genomics, reflecting the advantage of BD Rhapsody in capturing neutrophils.

Increase of normal lung epithelial cells and detection of residual cancer cells in pCR patients after combined therapy

We next investigated populations of epithelial cells. We first re-clustered the epithelial cells into 10 populations and separated malignant and normal cells using the CopyKAT algorithm [20] based on copy number variations (CNVs) (Fig. 2A, B and Additional file 2: Fig. S3A-B). Clusters E0_DST, E3_PCNA, E4_TOP2A, E7_SERPINB9, and a subpopulation of E1_KRT17 had higher CNV scores than other clusters and were inferred to be malignant cells (Fig. 2B and Additional file 2: Fig. S3C-D). The normal clusters were annotated as alveolar cells (E5_SFTPA2, type I: AGER, type II: SFTPA2), secretory club cells (E8_SCGB1A1), ciliated cells (E9_TPPP3), and basal epithelial cells (subpopulation of E1_KRT17) based on traditional markers (Fig. 2A, B) [21]. The fractions of alveolar cells (E5_SFTPA2), club cells (E8_SCGB1A1), and ciliated cells (E9_TPPP3) were increased after therapy, especially for MPR patients (Additional file 2: Fig. S3E-F). This indicated that combined therapy promoted expansion of normal epithelial cells after eliminating malignant cells. The normal epithelial cells may contribute to reconstruct normal lung structure in the previous tumor bed. In addition, it also has been reported that SCGB1A1+ club cells could increase the efficacy of ICB in lung cancer by promoting infiltration of cytotoxic cells [38].

Fig. 2figure 2

Epithelial cell reprograming after therapy. A UMAP plot of epithelia colored by clusters. The cells within the black dash line were malignant cells based on copy number variations (CNVs) inferred by the CopyKAT algorithm. B UMAP plots of epithelia colored by CopyKAT and normal lung epithelial markers. In the left top panel, the cells in red were predicted to be malignant cells and blue were normal cells. The cells in cluster E1_KRT17 contained both malignant and normal cells. C Boxplots of the average expression of CX3CL1, CD74, and HLA-DRA in malignant cells in TN (n = 3), MPR (n = 4), and NMPR (n = 7, one sample with less than 10 malignant cells was removed) patients. Center line indicates the median, lower, and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5× interquartile range. One-sided t-test was used, and the P values were adjusted by the FDR method. D Boxplots of the average expression of ARK1C1-3 in malignant cells in TN (n = 3), MPR (n = 4), and NMPR (n = 7) patients. One-sided t-test was used, and the P values were adjusted by the FDR method. E Boxplot of the average signature score in malignant cells in TN (n = 3), MPR (n = 4), and NMPR (n = 7) patients. One-sided t-test was used. F Longitudinal serums were collected from 24 patients (10 were assessed as MPR, and 14 as NMPR after surgery) at baseline, on-treatment and post-treatment timepoint. Non-targeted Metabolomic was conducted to detect the abundance of β-estradiol. G Boxplots of the β-estradiol abundance relative to baseline in 24 patients (10 patients were assessed as MPR and 14 as NMPR after surgery) at on-treatment and post-treatment timepoint. Two-sided unpaired Wilcoxon test was used. H Correlation analysis between the signature of AKR family genes and IC50 values in NSCLC cell lines under the condition of different drugs. P values were determined by two-sided Pearson correlation test

When comparing the cellular fraction of epithelial cells between different patients, we noted the enrichment of the malignant cluster E7_SERPINB9 in P06. This is unexpected, because P06 was classified as having a pathologic complete response (pCR; Additional file 2: Fig. S3E,G), which is defined by a complete absence of viable tumor cells upon H&E staining [5]. Although it is possible that this arose from sampling bias among the histopathologic slides, it is more likely that the sample from P06 contains malignant cells with genome alterations, but not morphological changes that can be detected by traditional histopathology. Our observation is consistent with previous reports that pCR patients may nonetheless experience tumor recurrence [39]. This reflects the presence of molecular residual disease (MRD) in pCR patients, a rising biomarker for NSCLC immunotherapy [40]. MRD is generally detected by circulating tumor DNA in serum. This case suggests that scRNA-seq has the ability in assessing MRD, which may be necessary to detection even for pCR patients.

Distinct molecular characteristics of malignant cells distinguish MPR and NMPR

To better characterize the malignant cell transcription programs activated in response to therapy, we performed differential expression analysis among TN, MPR, and NMPR patients. We were particularly interested in expression patterns that may drive interactions with the immune system and perform as signatures of therapy response. In response to therapy, malignant cells from MPR patients highly expressed CX3CL1, CD74, and major histocompatibility complex class II (MHC-II) genes (Fig. 2C and Additional file 2: Fig. S4A). Each key component of this MPR signature is addressed in turn below.

CX3CL1 is the ligand of CX3CR1. Previous studies have reported that CX3CR1 is highly expressed in many immune cells including NK cells [41] and monocytes [42]. CX3CL1 was downregulated in lung adenocarcinoma (LUAD) tumors compared to normal lung tissues in TCGA cohorts (Additional file 2: Fig. S4B), indicating tumor immune evasion. Therefore, cancer cells expressing CX3CL1 in response to therapy may promote immune cell infiltration into the TME, thereby improving overall survival (Additional file 2: Fig. S4D).

CD74 and MHC-II genes, also components of the MPR signature, are required for tumor antigen presentation [43]. We observed that the gene signature of antigen presentation via MHC-II was higher in cancer cells from MPR patients than TN or NMPR patients (Additional file 2: Fig. S4C). Previous studies have shown that MHC-II expression is associated with anti-PD-1 therapy response [44], progression-free, and overall survival in melanoma [45]. Consistent with these observations, higher expression of CD74 and HLA-DRA was associated with a better prognosis in TCGA-LUAD cohorts (Additional file 2: Fig. S4D).

Compared to TN and MPR patients, we observed that enzymes in the Aldo-Keto Reductase family (AKR1B1/10 and AKR1C1-3) were highly expressed in cancer cells from NMPR patients (Fig. 2D and Additional file 2: Fig. S4E-F). The AKR1B family has been previously reported to promote tumor metastasis and drug resistance [46,47,48], and the AKR1C family (hydroxysteroid dehydrogenases) was involved in estrogen metabolism, catalyzing the reduction of estrone to β-estradiol [49]. Consistent with this, Gene Set Variation Analyses (GSVA) revealed that following combined therapy, estrogen response pathways were upregulated in malignant cells from NMPR patients (Fig. 2E and Additional file 2: Fig. S4G). Only one of the eight (12.5%) NMPR patients was female, while two of the four (50%) patients in MPR group were female (Additional file 1: Table S1). None of the 15 patients used any estrogen-related drugs during therapy. Thus, the data suggest that estrogen metabolism may be aberrantly high in NMPR patients following treatment.

To validate this, we used non-targeted metabolomics to detect the abundance of steroids in serum from cells collected at baseline (before neoadjuvant therapy), on-treatment (at the first or second cycle, 3 weeks per cycle, total 2-4 cycle) and post-treatment (4 weeks after the last drug administration, blood samples were collected before surgery) in 10 MPR (30% female) and 14 NMPR (7% female) patients (Fig. 2F and Additional file 1: Table S1). In confirmation of the previous result, levels of β-estradiol were significantly elevated in NMPR patients compared to baseline during therapy (Fig. 2G, Additional file 2: Fig. S4H and Additional file 4: Table S3). When removing the patients from scRNA-seq cohort, the results were similar (Additional file 2: Fig. S4I). Thus, elevated estrogen levels in serum could reflect poor response to immunotherapy. Estradiol has been reported to be an immunosuppressor in the TME [50], through promoting the infiltration of M2 macrophages [51], mobilization of myeloid-derived suppressor cells (MDSCs) [52], and expansion of Tregs [53]. These suggested that estradiol may generate an immunosuppressive TME in the NMPR patients.

To identify potential drugs that may be effective on cancer cells in NMPR patients, we explored data in NSCLC cell lines from the Genomics of Drug Sensitivity in Cancer (GDSC) database. We found that the NMPR signature was negatively correlated with the IC50 (half the maximal inhibitory concentration) of 17-AAG (Fig. 2H), an inhibitor of HSP90, suggesting that cancer cells in NMPR patients may be sensitive to 17-AAG. Notably, 17-AAG is reported to inhibit estrogen signaling by disrupting HSP90 [54].

The degree of cytotoxic T/NK cell expansion and reduction of suppressive Tregs after combined therapy is positively associated with pathologic response

Next, we explored the dynamics of immune cell lineages in the TME in response to therapy. Since T cells are the most abundant tumor-infiltrating lymphocytes in the TME, we re-clustered T/NK cells and identified 14 clusters (Fig. 3A, B and Additional file 2: Fig. S5A). These includes 2 subtypes of NK cells (NK_FCGR3A and NK_KLRD1), 5 subtypes of CD8+ T cells (CD8_IL7R, memory T [Tm]; CD8_GZMK, effector memory T [Tem]; CD8_GZMB, Trm; CD8_HAVCR2, exhausted T [Tex] and CD8_STMN1, cycling effector T), 4 subtypes of conventional CD4+ T cells (CD4_CCR7, naïve T; CD4_IL7R, memory T; CD4_MAF, mature follicular helper T [Tfhs] [55]; and CD4_CXCL13, naïve Tfhs), 2 subtypes of regulatory T (Treg) cells (Treg_SELL, naïve-like Treg; Treg_CTLA4, activated Treg), and 1 proliferating subtype (T_MKI67).

Fig. 3figure 3

T/NK cell remodeling after therapy. A UMAP plot of T/NK cells colored by clusters. B Heatmap of normalized expression of canonical T/NK cell marker genes among clusters. TRM, tissue-resident memory. C Boxplots of the average cytotoxic and exhausted signature scores for CD8+ T cells in TN (n = 3), MPR (n = 4), and NMPR (n = 8) patients. Center line indicates the median, lower, and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5× interquartile range. One-sided t-test was used. D Boxplots of the average cytotoxic and exhausted signature scores for NK cells in TN (n = 3), MPR (n = 4), and NMPR (n = 8) patients. One-sided t-test was used. E Boxplot showing cellular fractions of each T/NK cluster in TN (n = 3), MPR (n = 4), and NMPR (n = 8) patients. All differences with adjusted P < 0.10 are indicated. One-sided unpaired Wilcoxon test was used and the P values were adjusted by the FDR method. F Summary of selected ligand-receptor interactions from CellPhoneDB between cancer cells and CD16+ NK cells in MPR patients. G Boxplots of the average exhausted signature scores for Tregs in TN (n = 3), MPR (n = 4), and NMPR (n = 8) patients. One-sided t-test was used. H The developmental trajectory of CD8+ T cells inferred by Monocle2. The memory CD8+ T cells (CD8_IL7R) and effector memory (CD8_GZMK) T cells were the roots of differentiation, and the exhausted CD8+ T cells (CD8_HAVCR2) were in the end-point state. I Heatmap of the top differential genes in memory (CD8_IL7R) cells along the pseudo-time (lower panel). The distribution of CD8_IL7R cells during the transition (divided into 2 phases: resting and activated) in TN, MPR, and NMPR patients, along with the pseudo-time (upper panel)

We further calculated the cytotoxic and exhausted gene signatures for CD8+ T cells and NK cells for TN, MPR, and NMPR patients. The cytotoxic and exhausted signatures were both significantly increased after therapy (Fig. 3C, D). This corresponds to an increase in the factions of all CD8+ T clusters after therapy (Fig. 3E). Relative to TN patients, the fractions of Tem (CD8_GZMK) and Trm (CD8_GZMB) and cycling effector T cells (CD8_STMN1) were increased in post-treatment samples, with the increase more pronounced in MPR patients (Fig. 3E). The CD8_GZMK cells can be resident in the TME and then locally expand after ICB, or newly infiltrate from peripheral blood [56]. The increase of Trm (CD8_GZMB) was consistent with expansion of neoantigen-specific T cells in NSCLC after immunotherapy [15].

We observed a higher fraction of Tex (CD8_HAVCR2) in both MPR and NMPR patients relative to TN patients following combined therapy. Recent studies have reported that exhausted T cells are specifically derived from tumor-specific T cells [57], and an increase in exhausted-like T cells is associated with ICB response [58]. To determine the source of exhausted T cells, we performed differential expression analysis before and after therapy in exhausted T cells (Additional file 2: Fig.S5B). We found that the transcription factors (TFs), NR4A2-3, that are associated with T cell exhaustion [59] were enriched in TN patients. This indicates that the T cells may have been exhausted before treatment, driven by NR4A2/3 during chronic T cell dysfunction. Cytotoxic (GZMH, NKG7, and PRF1) and exhausted markers (LAG3 and TIGIT) were both highly expressed in post-treatment patients. Tex cells that remain after treatment may arise from the coupled activation, expansion, and exhaustion process for cytotoxic T cells, which has been reported to be more evident in responders [18].

Cluster NK_FCGR3A was most representative of cytotoxic cells and was distinguished from NK_KLRD1 cells by expression of FCGR3A (CD16a), FGFBP2, and CX3CR1 (Fig. 3B) [41]. Given the expression of CX3CL1 in cancer cells from MPR patients (Fig. 2C), it was possible that NK_FCGR3A cells were recruited into the TME by CX3CL1. As expected, cell-cell interaction analysis using the CellPhoneDB algorithm [28] showed the CX3CL1-CX3CR1 interaction between cancer cells and NK_FCGR3A cells was significantly enriched in MPR patients (Fig. 3F).

We next focused on Tregs. Activated Tregs have been previously reported to have a stronger immunosuppressive function than naïve Tregs, and to be correlated with poor prognosis [60]. Consistently, activated Tregs (Treg_CTLA4, expressing TNFRSF4 and TNFRSF9) decreased only in MPR patients. The proportion of naïve Tregs (Treg_SELL, expressing SELL and LEF1) decreased in both MPR and NMPR patients relative to TN patients (Fig. 3E). MPR patients consistently showed lower Treg exhausted signature than NMPR patients (Fig. 3G).

Our analysis revealed the expansion and activation of cytotoxic T cells and CD16+ NK cells, and reduction of immunosuppressive Tregs after treatment. The strength of these trends was associated with positive response to combined therapy.

Therapy promotes the differentiation of memory T cells into an effector phenotype

After combined therapy, memory CD8+ T cells (CD8_IL7R) decreased while effector T cells increased (Fig. 3E). This suggested that treatment might directly induce the activation of memory T cells into a cytotoxic phenotype. To test this hypothesis, we performed trajectory analysis using Monocle2 [27]. One detected transition path went from Tem (CD8_GZMK) to Trm (CD8_GZMB) to Tex cells (CD8_HAVCR2; Fig. 3H). This path confirmed the sequential activation and exhaustion of CD8+ T cells in the TME. The analysis also showed that the cytotoxic cells may differentiate directly from Tm cells (CD8_IL7R). Two origins of cytotoxic T cells were confirmed using RNA velocity algorithm (Additional file 2: Fig. S5C), another trajectory analysis algorithm [25]. When we delineated the distribution of the CD8_IL7R cells in pseuso-time, we noticed that the CD8_IL7R cells could be categorized into 2 phases (Fig. 3I). The CD8_IL7R cells in resting phase highly expressed NR4A1-3, while cytotoxic-related genes (GZMA, NKG7, and CCL5) and MHC-II genes (CD74 and HLA-DRA) were upregulated in the activated phase (Fig. 3I). The proportion of activated cells was increased after therapy, and more CD8_IL7R cells were activated in MPR than NMPR samples (Fig. 3I and Additional file 2: Fig. S5D). Combined, these observations suggest that therapy could activate memory CD8+T cells into an effector phenotype, and the activation was most pronounced in MPR patients.

FCRL4+FCRL5+ memory B cells predict response to ICB and boost immunotherapy through activating CD4+ T cells

Studies indicate that B cells are actively involved in anti-tumor immunity after neoadjuvant chemotherapy [61]. To assess the B cell diversity after therapy, we re-clustered B cells into 7 subclusters (Fig. 4A and Additional file 2: Fig. S6A-B), including 5 subgroups of memory B cells (CD27+GPR183+IGHD-, B0_MS4A1, B1_IGHM, B2_HSP1A1, B4_FCRL4 and B5_CD83), 1 naïve B cell (CD27-IGHD+, B3_IGHD), and 1 germinal center (GC) B cell (B6_RGS13).

Fig. 4figure 4

B cell remodeling after therapy. A UMAP plot of B cells colored by clusters. B Boxplot showing cellular fractions of each B cluster in TN (n = 3), MPR (n = 4), and NMPR (n = 8) patients. Center line indicates the median, lower, and upper hinges represent the 25th and 75th percentiles, respectively, and whiskers denote 1.5× interquartile range. All differences with adjusted P < 0.10 are indicated. One-sided unpaired Wilcoxon test was used, and the P values were adjusted by the FDR method. C Violin plots of marker genes of B4_FCRL4 cells across clusters. D In situ multiplex immunofluorescence images of B4_FCRL4 cells in MPR and NMPR tumor tissues. E Violin and box plots of B4_FCRL4 signature in our validation cohort (9 patients were assessed as MPR and 12 as NMPR after surgery) before ICB + chemotherapy. One-sided unpaired Wilcoxon test was used. F Violin and box plots of B4_FCRL4 signature in responders (R) and non-responders (NR, removing SD patients) in advanced melanoma cohorts. Two-sided unpaired Wilcoxon test was used. G Kaplan–Meier survival curves of the signature of B4_FCRL4 in advanced melanoma cohort 2. Survival curves were compared by the log-rank test. H Summary of selected ligand-receptor interactions from CellPhoneDB among B4_FCRL4 cells, Tfhs and CD8+ T cells

To characterize the function of different B cell subsets, we compared their cell-type fractions in TN, MPR, and NMPR patients. Naïve B cells were increased after treatment, while memory B cells were slightly reduced (Additional file 2: Fig. S6C). Although memory B cells in general decreased following combined treatment, the FCRL4+FCRL5+ B cells (B4_FCRL4), defined as “atypical memory B cells” [62], exhibited a slightly increasing trend in MPR patients (Fig. 4B). The FCRL4 and FCRL5 genes encode the Fc receptors for IgA and IgG, respectively, and are drivers of human memory B cell activations [62]. B cells expressing FCRL4 have been previously reported to be associated with inflammation in rheumatoid arthritis [63] and viral infections [64]. Among TCGA-LUAD patients, we consistently found that patients with high expression of FCRL4 and FCRL5 had a better prognosis (Additional file 2: Fig. S6D-E). Also highly expressed in B4_FCRL4 cells were interferon-stimulated genes (CCR1, STAT1, and GBP4), co-stimulatory molecule (CD86), and activated follicular B cell marker (BHLHE40) [65]. Consistently, immunofluorescence staining showed that FCRL4+FCRL5+ cells were much more enriched in MPR than NMPR samples (Fig. 4D and Additional file 2: Fig. S6F). Interestingly, we noticed that CD20+ B cells aggregated in TLS and FCRL4+FCRL5+ B cells located in the center of the TLS (Fig. 4D and Additional file 2: Fig. S6G). Taken together, our analysis suggests that FCRL4+FCRL5+ B cells are associated with anti-tumor activity and a positive response to combined therapy.

We investigated whether the signature from FCRL4+FCRL5+ B cells could serve as a positive biomarker for immunotherapy. We first evaluated the B4_FCRL4 gene signature (Additional file 3: Table S2) in our validation cohort. The signature scored significantly higher in MPR patients before ICB combined with chemotherapy (Fig. 4E). We then performed similar analyses on published datasets from two independent melanoma cohorts with ICB treatment [34, 66]. Although the melanoma cohorts were not in the neoadjuvant setting, the B4_FCRL4 signature also perform well in predicting immunotherapy response. The B4_FCRL4 signature was higher in responders (complete response or partial response) than non-responders (stable disease or progressive disease) before and after therapy in both cohorts (Fig. 4F). Higher B4_FCRL4 signature was associated with improved survival in previously published “melanoma cohort 2” [34] (Fig. 4G). These results indicate that the signature of FCRL4+FCRL5+ B cells can be used as biomarker for predicting response to ICB.

To explore the underlying mechanisms for the activation and function of FCRL4+FCRL5+ B cells, we performed NicheNet analysis [29], which predicts ligands driving the transcriptomic changes of target cells. Several IFNα genes, tumor necrosis factor (TNF), and IL27 were predicted as possible ligands driving the phenotype of B4_FCRL4 cells (Additional file 2: Fig. S6H). CellPhoneDB analysis revealed that FCRL4+ FCRL5+ B cells could interact with Tfhs through ligand-receptor pairs: CXCL13-CXCR5, CD40-CD40LG, and CD28-CD86 (Fig. 4H). CXCL13-CXCR5 interaction between B cells and Tfhs is essential for the formation of TLS [

留言 (0)

沒有登入
gif