Single-Cell RNA Sequencing Reveals Tumor Microenvironment Remodeling after Neoadjuvant Immunotherapy in Non-small Cell Lung Cancer

doi:10.21203/rs.3.rs-934132/v1

Download PDF

Resource

Single-Cell RNA Sequencing Reveals Tumor Microenvironment Remodeling after Neoadjuvant Immunotherapy in Non-small Cell Lung Cancer

https://doi.org/10.21203/rs.3.rs-934132/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Immunotherapy has revolutionized cancer treatment, but most patients are refractory to immunotherapy or acquire resistance. To explore immunotherapy resistance mechanisms, we characterized the transcriptomes of ~92000 single cells from 15 patients with non-small cell lung cancer (NSCLC) during neoadjuvant PD-1 blockade combined with chemotherapy. CD8+ T, natural killer, B, and dendritic cells were activated by therapy. Therapy also promoted differentiation of memory T cells into effector T cells. Macrophages were remodeled into an M0-like phenotype and neutrophils into an aged phenotype. Distinct therapy-induced cancer-cell transcriptomes were associated with clinical response. Major pathologic responders (MPRs) activated antigen presentation via major histocompatibility complex class II (MHC-II). Cancer cells of non-MPRs exhibited overexpression of estrogen metabolism enzymes and elevated serum estradiol. Elevated estradiol activated EGFR signaling and upregulated the expression of VEGFA, which promoted an immunosuppressive microenvironment. FCRL4+FCRL5+ memory B cells and CD16+CX3CR1+ monocytes were identified as biomarkers for positive immunotherapy response.

Oncology

Medical Genetics

Cancer Biology

neoadjuvant PD-1 blockade

effector T cells

macrophages

Lung cancer is the leading cause of cancer-related death worldwide ¹, with non-small cell lung cancer (NSCLC) representing approximately 85% of lung cancer cases ². Cancer immunotherapy through immune check point blockade (ICB) is the first-line treatment for advanced NSCLC without an identified driver-gene mutation ³. For resectable NSCLC, immunotherapy prior to surgery (neoadjuvant immunotherapy) is emerging as a promising therapeutic regimen ⁴. A common measure of neoadjuvant immunotherapy efficacy is the “major pathologic response” (MPR), which is defined as having no more than 10% residual viable tumor cells by routine hematoxylin and eosin (H&E) staining after therapy ⁵. Despite the benefits of immunotherapy, its efficacy by this measure is limited. The mean MPR rate of neoadjuvant anti-PD-1/PD-L1 immunotherapy is approximately 32% (range 18% to 63%) ⁶. Most patients are refractory to therapy or acquire resistance, and the underlying mechanisms remain to be explored.

The tumor microenvironment (TME) plays an important role in tumor development, progression, metastasis, and drug resistance ⁷. Immunotherapy remodels the TME, and the TME in turn influences the response to immunotherapy ⁸. Previous studies have characterized TME remodeling after ICB and associated the changes with clinical outcomes. ICB treatment overcomes T-cell dysfunction or exhaustion and promotes clonal expansion of T cells ⁹. Patients with clonotypic expansion of T cells respond best to ICB therapy ¹⁰. Studies also suggest that the abundance of CD8+TCF7+ T cells and EOMES+CD69+CD45RO+ T cells in the TME before treatment predicts ICB response and better survival ^11,12. Other studies have indicated that the formation of tertiary lymphoid structures (TLS), which are aggregations of B cells, T cells, dendritic cells (DC), and high endothelial venules, promote immunotherapy response ¹³. The myeloid component also has been reported to be associated with immunotherapy response, with a subset of CD73+ macrophage persistence marking therapy resistance ¹⁴. Most previous studies of immunotherapy and the TME focus on easily accessible cancer types such as melanoma or skin carcinoma. However, a recent study focusing on mutation-associated neoantigen (MANA) specific T cells in NSCLC reported that these cells express an incompletely activated cytolytic program in the TME during immunotherapy ¹⁵. MANA-specific clones from non-responders were largely confined to tissue-resident memory subsets, which exhibited upregulated checkpoints and inhibitory receptors. However, this study did not address the response of the entire TME in immunotherapy in NSCLC.

The dynamics of cancer cells and other immune cells, and a more comprehensive characterization of the TME of NSCLC during ICB treatment, is of interest for more accurately predicting patient response and providing novel therapeutic targets. Also of interest is characterization of immunotherapy in combination with other drugs, which is increasingly common in clinical use. To explore mechanisms of immunotherapy resistance and their relation to changes in the TME after PD-1 blockade combined with chemotherapy, we performed single-cell RNA sequencing (scRNA-seq), bulk RNA-seq and untargeted metabolomics from NSCLC samples.

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, Extended Data Fig. 1a and Supplementary Table 1). 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 ⁵. The dataset analyzed here also includes bulk RNA-seq from fresh biopsies from 24 independent TN patients (Supplementary Table 1).

Single-cell transcriptomes were generated using BD Rhapsody Single-Cell Analysis System (BD Biosciences). After quality control and removal of doublets, transcriptomes from 92330 cells with a median of 1256 genes per cell were used for further analyses. To mitigate possible batch effects from patients 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 Extended Data Fig. 1b), with no significant batch effects observed across different patients, PD-1 antibodies or response groups (Fig. 1c,d). Further, the average gene numbers and unique molecular identifiers (UMIs) were comparable between different clusters (Extended Data Fig. 1d). We then annotated the 26 clusters into T cells, NK cells, B cells, myeloid cells, neutrophils, plasma cells, plasmacytoid DC (pDC), mast cells, cycling immune cells, stromal cells (fibroblasts/endothelia) and epithelia, according to the expression of corresponding canonical marker genes (Fig. 1b and Extended Data Fig. 1c).

Combined Therapy Promotes Infiltration of Anti-tumor Cells and Decreases the Prevalence of Immunosuppressive Cells, Especially in MPR Patients

To characterize the TME remodeling in response to treatment, we calculated the fraction of different cell types in TN, MPR, and NMPR patients. We observed that the fraction of T cells, NK cells and B cells increased after therapy, especially for MPR patients (Fig. 1d,e and Extended Data Fig. 1e). This indicates that PD-1 blockade combined with chemotherapy promotes the infiltration of lymphoid cells into TME, consistent with previous studies of immunotherapy ^9,13 and chemotherapy ¹⁶.

We observed a significant increase in the share of myeloid cells in NMPR patients following therapy. The increase in the myeloid cell fraction in patients with poor response to therapy suggests an immunosuppressive function for myeloid cells. However, myeloid cells are known to have diverse and complex functions in the TME ¹⁷, which are further explored later in this study. Therapy resulted in reduction of neutrophils and mast cells in all post-treatment samples (Fig. 1d,e and Extended Data Fig. 1e). Neutrophils typically serve a pro-tumor role by interfering with the contacts between cytotoxic T/NK cells and tumor cells through extracellular traps (NETs) in immunotherapy ¹⁸. High mast cells infiltration either before or after treatment has been correlated with poor response in chemotherapy in breast cancer ¹⁹.

The fraction of epithelial cells (which includes cancer cells) and stromal cells decreased after therapy, because therapy induced the infiltration and expansion of lymphocytes and eliminated tumor cells. As described above, this infiltration was more pronounced in MPR patients, and accordingly MPR patients a had lower proportion of epithelial and stromal cells than TN or NMPR patients (Fig. 1d-e and Extended Data Fig. 1e).

Thus, PD1 blockade combined with chemotherapy significantly remodeled the TME by promoting infiltration of anti-tumor cells and decreasing the prevalence of immunosuppressive cells. This pattern was most pronounced in MPR patients.

Regeneration of Normal Lung Structures 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 (Fig. 2a and Extended Data Fig. 2a) and separated malignant and normal cells using the CopyKAT algorithm ²⁰ by copy number variations (CNVs) (Fig. 2b and Extended Data Fig. 2b,c). Clusters E0-DST, E3-PCNA, E4-TOP2A, E7-SERPINB9, and a subpopulation of E7-KRT17 had higher CNV scores than other clusters, and were inferred to be malignant cells (Fig. 2b and Extended Data Fig. 2c,d). The normal clusters were annotated as alveolar cells (E5-SFTPA2, type I: AGER, type II: SFTPA2), secretory club cells (E8-SCGB3A1), ciliated airway epithelial cells (E9-TPPP3) and basal airway epithelial cells (subpopulation of E1-KRT17) based on traditional markers (Fig. 2b and Extended Data Fig. 2a) ²¹.

The fractions of alveolar cells, club cells, and ciliated cells were significantly increased after therapy, especially for MPR patients (Extended Data Fig. 2e,f). We interpret this to reflect the regeneration of normal lung structure at the tumor site. We also observed that most alveolar cells were AT2 cells (Fig. 2b and Extended Data Fig. 2c). Previous studies reported that AT2 cells could act as stem-like progenitor cells, which become active and proliferate after diverse types of lung injury ²². Our results suggest that after therapy and the resulting anti-tumor inflammation, the rebuilding of normal lung structure is mediated by AT2 cells.

When comparing the cellular fraction of epithelial cells between different patients, we noted the enrichment of the malignant cluster E7-SERPINB9 in Patient6. This is unexpected because Patient6 was classified as having a pathologic complete response (pCR; Extended Data Fig. 2e,g), which is defined by a complete absence of viable tumor cells upon H&E staining ⁵. Although it is possible that this arose from sampling bias among the histopathologic slides, it is more likely that the sample from Patient6 contains malignant cells with genome alterations, but not morphological changes that can be detected by traditional histopathology. This observation suggests that scRNA-seq may be more sensitive in assessing residual tumor cells than conventional H&E staining. Our observation is consistent with previous reports that pCR patients may nonetheless experience tumor recurrence ^23,24. A consensus strategy for management of pCR patients after neoadjuvant immunotherapy has not been established due to lack of definitive clinical evidence, although one year of maintenance immunotherapy has been recommended ²⁵. Our observation of malignant cells by scRNA-seq that were not detectable by traditional histopathology provides support for maintenance immunotherapy 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. We found that the transcriptomes of malignant cells in MPR showed evidence of antigen presentation activation, immune cell chemotaxis, and anti-GrB activity. NMPR transcriptional signatures indicated aberrant estrogen metabolism. As further detailed below, these observations suggest novel combinatorial treatment strategies.

In response to therapy, malignant cells from MPR patients highly expressed SPARCL1, CX3CL1, CD74, SERPINB9 and major histocompatibility complex class II (MHC-II) genes (Fig. 2c, Extended Data Fig. 3b and Supplementary Table 2). Each key component of this MPR signature is addressed in turn below. In previous studies, SPARCL1 has been reported to suppress osteosarcoma metastasis and recruit M1 macrophages by activating WNT/β-catenin signaling ²⁶. Further supporting a connection between higher SPARCL1 expression and better outcomes, SPARCL1 was downregulated in lung adenocarcinoma (LUAD) tumors compared to normal lung tissues in TCGA cohorts (Extended Data Fig. 3a), and higher expression of SPARCL1 within tumors was associated with better prognosis (Extended Data Fig. 3d).

CX3CL1, another of the highly expressed genes in MPR samples, is the ligand of CX3CR1. Previous studies have reported that CX3CR1 is highly expressed in many immune cells including NK cells ²⁷ and monocytes ²⁸. Therefore, cancer cells expressing CX3CL1 in response to therapy may promote immune cell infiltration into the TME, thereby improving overall survival (Extended Data Fig. 3e).

CD74 and MHC-II genes, also components of the MPR signature, are required for tumor antigen presentation ²⁹. We observed that the gene signature of antigen presentation via MHC-II (Supplementary Table 2) was higher in cancer cells from MPR patients than TN or NMPR patients (Extended Data Fig. 3c). Previous studies have shown that MHC-II expression is associated with anti-PD-1 therapy response ³⁰, progression-free and overall survival in melanoma ³¹. Consistent with these observations, higher expression of CD74 and HLA-DRA was associated with a better prognosis in TCGA-LUAD cohorts (Extended Data Fig. 3f,g).

SERPINB9 was highly expressed in MPR patients and a cluster of E7-SERPINB9 cells was enriched in MPR patients (Extended Data Fig. 2f). SERPINB9 is an inhibitor of granzyme B (GrB) that protects cancer cells from GrB-induced apoptosis ³². Therefore, the enrichment of SERPINB9 among MPR patients was surprising to us. We hypothesize that the high SERPINB9 expression may arise from the residual cancer cells that survived in part by secreting SERPINB9 to escape the attack from GrB. This may represent a novel strategy to eliminate the residual cancer cells in MPR patients. For example, the combination of SERPINB9 inhibitors with ICB may further the survival of MPR patients after neoadjuvant immunotherapy, while the consensus recommends only ICB maintenance therapy ²⁵.

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 Supplementary Table 2). The AKR1B family has been previously reported to promote tumor metastasis and drug resistance ^33-35. The AKR1C family (hydroxysteroid dehydrogenases) has been previously reported to be involved in estrogen metabolism, catalyzing the reduction of estrone to estradiol ³⁶. Consistent with the function of the AKR1C family, Gene Set Variation Analyses (GSVA) revealed that following combined therapy, estrogen response pathways were up-regulated in malignant cells from NMPR patients (Fig. 2e and Extended Data Fig. 3h). There was no significant bias to female patients (Supplementary Table 2) and none of the 15 patients used any estrogen-related drugs during therapy. This indicates that the estrogen metabolism may be aberrantly high in NMPR patients following treatment.

To validate this, we used untargeted metabolomics to detect the abundance of steroids in serum from cells collected at baseline, on-treatment (after the first of 2-4 treatment cycles) and pre-surgery (after the last cycle) in 11 MPR and 14 NMPR patients (Supplementary Table 2). In confirmation of our previous result, levels of estradiol were significantly elevated in NMPR patients compared to baseline during therapy (Fig. 2f and Extended Data Fig. 3i). Estradiol has been reported to be an immunosuppressor in the TME ³⁷, through promoting the infiltration of M2 macrophages ³⁸, mobilization of myeloid-derived suppressor cells (MDSCs) ³⁹, and expansion of Tregs ⁴⁰. Estradiol has also been reported to drive the activation of the EGFR pathway and upregulate VEGFA to promote tumor progression in NSCLC ⁴¹. We consistently observed that the cancer cells in NMPR patients had significantly higher signature of EGFR signaling (Supplementary Table 2) and expression of VEGFA than TN and MPR patients (Fig. 2g,h). NSCLC harboring EGFR mutation was reported to have a poor response to ICB treatment ⁴². VEGFA has been reported to mobilize MDSCs and Tregs, and reduce the ability of antigen-presenting cells (APCs) and the anti-cancer responses of effector T cells ⁴³.

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 (Extended Data Fig. 3j), 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 ⁴⁴. Thus, elevated estrogen levels in serum could be a negative biomarker reflecting poor response to immunotherapy, and anti-estrogen drugs may improve the efficacy of immunotherapy among these patients.

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 cells 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 13 clusters (Fig. 3a). Each cluster was annotated based on marker genes (Fig. 3a,b, Extended Data Fig. 4a,b and Supplementary Table 4): including 2 subtypes of NK / natural killer T (NKT) cells (T6-FCGR3A and T9-KLRD1), 4 subtypes of CD8+ T cells (T0-GZMK, pre-effector T cells ; T1-GZMB, effector T cells; T7-HAVCR2, pre-exhausted T cells and T11-STMN1, cycling effector T cells), 3 subtypes of conventional CD4+ T cells (T4-CCR7, naïve T cells; T5-MAF cytotoxic CD4+ T cells; and T8-CXCL13, follicular helper T cells [Tfhs]), 2 subtypes of regulatory T (Treg) cells (T3-FOXP3 and T12-ILRA), 1 proliferating subtype (T10-MKI67) and 1 memory subtype (T2-IL7R) with both CD8A and CD4 signature.

We further calculated the cytotoxic and exhausted gene signatures for CD8+ T cells and NK/NKT 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 cellular factions of all CD8+ T clusters after therapy (Fig. 3e). Relative to TN patients, the fractions of pre-effector T cells (T0-GZMK) and effector T cells (T1-GZMB) were significantly increased in post-treatment samples, with the increase of effector T cells more pronounced in MPR patients (Fig. 3e). The faction of NK/NKT cells (T6 and T9) was increased only in MPR patients (Fig. 3e). This indicates that the cytotoxic T cell subsets and NK cells are possibly responsible for tumor cell killing and associated with a positive response to therapy.

We observed a higher fraction of pre-exhausted T cells (T7-HAVCR2) in both MPR and NMPR patients relative TN patients following combined therapy. To determine the source of pre-exhausted T cells, we performed differential expression analysis before and after therapy in pre-exhausted T cells (Extended Data Fig. 4c). We found that the transcription factors (TFs), NR4A1-2, that have been associated with T cell exhaustion ⁴⁵ were enriched in TN patients. This indicates that the T cells may have been exhausted before treatment, driven by NR4A1/2 during chronic T-cell dysfunction. Cytotoxic (GZMH, NKG7 and PRF1) and exhausted markers (LAG3 and TIGIT) were both highly expressed in post-treatment patients (Extended Data Fig. 4c). Pre-exhausted T 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 ⁴⁶.

Cluster T6-FCGR3A was most representative of cytotoxic cells and was distinguished from T9-KLRD1 cells by expression of FCGR3A, FGFBP2, and CX3CR1 (Fig. 3b and Extended Data Fig. 4a,b) ²⁷. Given the expression of CX3CL1 in cancer cells from MPR patients (Fig. 2c), it is possible that NK cells were recruited into the TME by CX3CL1. As expected, cell-cell interaction analysis using the CellPhoneDB algorithm ⁴⁷ showed the CX3CL1-CX3CR1 interaction between cancer cells and NK cells was significantly enriched in MPR patients (Fig. 3f).

We next focused on CD4+ T cells. Both naïve and cytotoxic CD4+ T cells decreased after therapy. CXCL13-expressing Tfhs in MPR patients showed an increasing trend (Fig. 3e). A previous report suggested that CXCL13- and CD40LG-expressing CD4+ T cells could attract B cells and form TLS ⁴⁸. TLS formation could then promote an immunotherapy response ¹³. Activated Tregs have been previously reported to have a stronger immunosuppressive function than naïve Tregs, and to be correlated with poor prognosis ⁴⁹. Consistent with this, activated Tregs (T12-IL2RA, expressing TNFRSF4 and TNFRSF9) decreased only in MPR patients. The proportion of naïve-like Tregs (T3-FOXP3+, also expressing SELL and LEF1) decreased in both MPR and NMPR patients relative to TN patients (Fig. 3e). NMPR patients consistently showed higher exhausted signature than TN and MPR patients (Fig. 3g).

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

Therapy Promote the Differentiation of Memory T cells to Effector T cells

After combined therapy, memory T cells (T2-IL7R) decreased while effector T cells increased. This suggested that treatment might directly induce the differentiation of T cells in the TME (Fig. 3e). To test this hypothesis, we performed trajectory analysis using RNA velocity algorithm ⁵⁰. One detected transition path went from pre-effector T cells (T0-GZMK) to effector T cells (T1-GZMB) to pre-exhausted T cells (T7-HAVCR2; Fig. 3h). This path confirmed the sequential activation and exhaustion of CD8+ T cells in the TME. The pre-effector T cells (T0-GZMK) can be resident in the TME, or may have infiltrated from peripheral blood ⁹ and lymph nodes ⁵¹. The analysis also showed that the effector cytotoxic cells (T1-GZMB) may differentiate directly from memory T cells (T2-IL7R). Two origins of effector cytotoxic T cells were confirmed using Monocle2 (Extended Data Fig. 5d), another trajectory analysis algorithm ⁵². Differential expression analysis between TN and post-treatment samples showed that the memory T cells in TN patients were enriched for NR4A1-3 (Fig. 3i). Relative to TN patients, cytotoxic-related genes such as GZMA, CCL5 and genes related to antigen presentation such as CD74 and HLA-DRA were up-regulated in memory T cells from post-treatment patients (Fig. 3i). Combined, these observations suggest that therapy could increase cytotoxic CD8+ T cells by promoting infiltration of pre-effector T cells and by differentiation of memory T cells.

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 ⁵³. To assess the B-cell diversity after therapy, we re-clustered B cells into 7 sub-clusters (Fig. 4a, Extended Data Fig. 5a,b and Supplementary Table 5), including 5 subgroups of memory B cells (CD20+CD27+GPR183+IGHD-, B0-MS4A1, B1-IGHM, B2-HSP1A1, B4-FCRL4 and B5-CD83), 1 naïve B cell (CD20+CD27-IGHD+, B3-IGHD) and 1 germinal center (GC) B cell (B6-RGS13) .

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 increased after treatment, while memory B cells were slightly reduced (Fig. 4b and Extended Data Fig. 5c). In MPR patients, GC B cells increased, and the proportion was significantly higher in MPR than NMPR patients (Extended Data Fig. 5c). We observed that GC B-cell marker RGS13 was associated with improved survival in TCGA-LUAD patients (Extended Data Fig. 5d). These analyses were consistent with previous studies showing that PD-1 blockade could increase germinal center formation ⁵⁴.

Although memory B cells in general decreased following combined treatment, the FCRL4+FCRL5+ B cells (B4-FCRL4), defined as “atypical memory B cells” ⁵⁵, were significantly increased in MPR patients but not NMPR patients (Fig. 4b,c). The FCRL4 and FCRL5 genes encode the Fc receptors for IgA and IgG, respectively, and are drivers of human memory B-cell activations ⁵⁵. B cells expressing FCRL4 have been previously reported to be associated with inflammation in rheumatoid arthritis ⁵⁶ and viral infections ⁵⁷. Among TCGA-LUAD patients, we consistently found that patients with high expression of FCRL4 and FCRL5 had a better prognosis (Extended Data Fig. 5e,f). 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) ⁵⁸. Consistently, immunofluorescence staining showed that FCRL4+FCRL5+ cells were much more enriched in MPR than NMPR samples (Fig. 4d). Interestingly, we noticed that CD20+ B cells aggregated in tumors and FCRL4+FCRL5+ B cells located in the center of the CD20+ B cells (Fig. 4d). Taken together, our analysis suggests that FCRL4+FCRL5+ B cells are associated with anti-tumor activity, and are associated with 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 (Supplementary Table 5) in our bulk RNA-seq data from pre-treatment patients. The signature scored significantly higher in MPR patients (Fig. 4e). We then performed similar analyses on published datasets from two independent melanoma cohorts with ICB treatment ^59,60. Again, 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 and Extended Data Fig. 5g). Higher B4-FCRL4 signature was associated with improved survival in previously published “melanoma dataset 2” ⁵⁹ (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 ⁶¹, which predicts ligands driving the transcriptomic changes of target cells. Several IFNα genes, IFNγ genes and tumor necrosis factor (TNF) were predicted as possible ligands driving the phenotype of B4-FCRL4 cells (Extended Data Fig. 5h). CellPhoneDB analysis revealed that FCRL4+ FCRL5+ B cells could interact with CD4+ T cells through ligand-receptor pairs: CXCL13-CXCR5, CD40-CD40LG and CD28-CD86 (Fig. 4H). CXCL13-CXCR5 interaction between B cells and CD4+ T cells is essential for the formation of TLS ⁴⁸. The data suggest that in the TME of immunotherapy responders, FCRL4+ FCRL5+ B cells are activated by IFNα, IFNγ, and TNF signals, and that these B cells interact with CD4+ T cells to enhance anti-tumor immunity through TLS formation.

Combined Therapy Expanded Tissue Regenerative Macrophages, Reprogramed TAM into a Neutral Phenotype, and Inhibited the Immunosuppressive Function in MPR Patients

The myeloid component in the TME exhibited remarkable heterogeneity and accordingly was categorized into 11 clusters (Fig. 5a,b, Extended Data Fig. 6a and Supplementary Table 6), including 2 subtypes of monocyte (C0-VCAN and C7-CX3CR1), 6 subtypes of macrophage (C1-SPP1, C2-CXCL9, C3-SELENOP, C4-FABP4, C6-MKI67 and C8-C1QA), and 3 subtypes of DC (C5-CD1C, C9-LAMP3 and C10-XCR1) according to canonical marker genes.

We first focused on the macrophages. C4-FABP4 cells are alveolar macrophages (AM) with the expression of the canonical AM markers (FABP4, MCEMP1 and MARCO) ⁶². Concordant with its tissue repair function ⁶³, C4-FABP4 AMs were significantly elevated in post-treatment patients, although to a greater extent in MPR patients (Fig. 5c). Given tissue repair function of AT2 cells and the increase of AT2 cells after therapy (Extended Data Fig. 2f), AMs may work together with AT2 cells to regenerate normal lung structure.

Other macrophage subtypes show similarity to known tumor-associated macrophages (TAM). C1-SPP1 macrophage has been previously reported to be associated with tumor angiogenesis ⁶⁴ and to facilitates immune escape by upregulating PD-L1 ⁶⁵. C3-SELENOP macrophage was previously reported to have anti-inflammation roles ⁶⁶. In contrast, C2-CXCL9 cells over-expressed pro-inflammatory factors (CXCL9 and CXCL10) that could attract T cells, NK cells and DCs ⁶⁷. The M2-like C1-SPP1 and C3-SELENOP macrophages decreased after therapy, while C1-SPP1 macrophages increased in NMPR patients (Fig. 5c).

Macrophages are traditionally classified into three subtypes: non-polarized or neutral (M0), pro-inflammatory or anti-tumor (M1), or anti-inflammatory or pro-tumor (M2). To better characterize the phenotype of macrophage subsets, we calculated M0, M1 and M2 signature scores based on the gene signatures from CIBERSORT ⁶⁸. As expected, the AM exhibited a M0-like phenotype, both C1-SPP1, C3-SELENOP and C8-C1QA had a stronger M2 signature. Only C2-CXCL9 cells showed a high M1 signature (Fig. 5d). Lineage tracing analysis suggested two distinct differentiation paths for monocytes. One path leads to M1-like C2-CXCL9 cells and another path differentiated into M2-like C3-SELENOP and C1-SPP1 cells (Fig. 5e). Comparing the macrophage signature to TN patients (Fig. 5f), the M0 signature increased significantly in both MPR and NMPR patients. However, the M2 signature decreased only in the macrophages of MPR patients. Surprisingly, M1 signature did not increase as expected, but even decreased significantly in the macrophages of MPR patients. We consistently observed that the fraction of M1-like CXCL9 cells did not increase in both MPR and NMPR patients (Fig. 5c). When comparing the three signatures in only TAMs, we found similar results (Extended Data Fig. 6b).

These analyses suggest that combined therapy induces expansion of tissue regenerative macrophages, and reprograms TAM into a neutral instead of an anti-tumor phenotype, but inhibits the immunosuppressive function only in MPR patients. Suppressing the M2-like function of TAMs may be more effective than promoting M1-like activity to enhance immunotherapy response.

Dendritic Cells Were Activated by Therapy and Expanded in MPR Patients

The DCs were classified into 3 subtypes, including conventional type I DCs (cDC1, C10-XCR1), conventional type II DCs (cDC2, C5-CD1C), and recently described LAMP3+ DCs (C9-LAMP3) (Fig. 5g and Extended Data Fig. 6a) ⁶⁹. cDC1 and cDC2 have been previously reported to activate CD8+ T cells and CD4+ T cells, respectively ⁶⁹. LAMP3+ DC was reported to be “mature DCs enriched in immunoregulatory molecules” (mregDC) ⁷⁰ due to expression of mature (LAMP3), migration (CCR7 and FSCN1) and immunoregulatory (CD274, PDCD1LG2 and CD200) markers, and downregulation of Toller-like receptors.

mregDCs have been reported to interact with tumor-infiltrating Treg cells or to inhibit CD8+ T cell-mediated tumor immunity by IL-4 stimulation ⁷⁰. Trajectory analysis indicated that mregDCs may be generated from both cDC1s and cDC2s (Fig. 5h), consistent with recent reports that mregDCs are derived from cDC1s and cDC2s upon uptake of tumor antigens ^70,71. After combined therapy, the fraction of cDC1s and cDC2s increased in MPR patients (Fig. 5c). In addition, antigen-presenting signature ⁷² of DCs was significantly increased after therapy, especially for MPR patients (Fig. 5i), suggesting that DCs were activated after therapy. mregDCs also showed an increasing trend in MPR patients, suggesting the conversion of cDC1 and cDC2 to mregDCs also occurred in MPR patients (Fig. 5c). However, the immunosuppressive signature did not show a significant difference between MPR and NMPR patients (Fig. 5i), possibly due to lack of IL4 receptor after treatment (Extended Data Fig. 6c,d). Thus, although mregDCs are increased in MPR patients, they do not contribute to immunosuppression in MPR patients due to lack of IL-4 signaling activation. The data indicates that the DCs are activated after therapy, and expanded in MPR patients, which could contribute to the activation of CD8+ and CD4+ T cells in the TME.

Patrolling Monocytes Predict Immunotherapy Response

Cluster C7-CX3CR1 highly expressed monocyte markers (FCN1, VCAN, S100A8 and S100A9), naïve marker (SELL) and lower MHC-II molecules (Fig. 5b), representing a “naïve-like” status. Cluster C0-VCAN had lower expression of monocyte markers, higher macrophage markers (MRC1, CD163 and MSR1) and MCH-II molecules than C7-CX3CR1 cells, suggesting a “pre-macrophage” status (Fig. 5b). Trajectory analysis validated that C0-VCAN was an intermediate phenotype between naïve monocytes and macrophages (Fig. 5e).

The C7-CX3CR1 monocytes closely resembled non-classical ‘‘patrolling monocytes (PMos)’’ (CD14+CD16+CX3CR1+, Fig. 6a) ²⁸. Previous studies show that PMos scavenge tumor material from the lung vasculature, and promote NK cell recruitment and activation ²⁸. C7-CX3CR1 cells were of increased abundance in MPR patients (Fig. 5c). C7-CX3CR1 cells also highly expressed CFP (complement factor properdin; Fig. 6b,c), which has been reported to suppress breast cancer growth by inducing apoptosis ⁷³. The expression of CFP was strongly downregulated in TCGA-LUAD patients compared to normal lung tissues (Extended Data Fig. 6e), and higher CFP expression was associated with better survival (Extended Data Fig. 6f).

We then investigate whether the gene signature from C7-CX3CR1 monocytes (Supplementary Table 6) could be used as a biomarker to predict ICB response using bulk RNA-seq data. Higher C7-CX3CR1 signature was observed in MPR patients in this study and responders from two independent melanoma cohorts (Fig. 6c,d and Extended Data Fig. 6g). Higher C7-CX3CR1 signature was associated with improved survival in Melanoma dataset 2 ⁵⁹ (Fig. 6f). Cell-cell interaction analysis showed that the CX3CL1-CX3CR1 pair was significantly enriched for C7-CX3CR1 cells (Fig. 6g), indicating that CX3CL1 expression in cancer cells attract C7-CX3CR1 monocytes in MPR patients. Our analysis identified CX3CR1+ monocytes as associated with anti-tumor activity, and a gene signature linked to immunotherapy response.

Therapy Induced Remodeling of Mature Neutrophils into Aged Phenotype

Neutrophils were abundant in our samples. They were divided into 6 sub-clusters (Fig. 7a, Extended Data Fig. 7a, and Supplementary Table 7), including 3 mature subsets (CD16+CXCR2^highCXCR4^low; including G0-PTGS2, G3-S100A12 and G4-FKBP5), and 3 aged subsets (CD16+CXCR2^lowCXCR4^high; including G1-CCL3, G2-IFIT3 and G5-CSTB) (Fig. 7b,c) ⁷⁴. Mature G0-PTGS2 neutrophils were characterized by high expression of PTGS2, a rate-limiting catalyzer for production of prostaglandin E2 that has been implicated in T cell suppression in cancer ⁷⁵. Mature G3-S100A12 neutrophils highly expressed S100A12, MMP9 and ARG1, which are associated with immune suppression ⁷⁴. Previous studies suggest that neutrophils could promote tumor growth and metastasis through NET ⁷⁶. NET pathway genes were up-regulated in mature G3-S100A12 and G4-FKBP5 cells (Fig. 7c). This suggests that mature neutrophils have strong tumor-promoting and immunosuppression functions in the TME.

Among the aged subsets, G1-CCL3 cells overexpressed multiple chemokines, including CCL3, CCL4 and CCL4L2. These cytokines exert both antitumor and pro-tumor effects ^77,78. G2-IFIT3 cells expressed interferon-stimulated genes (IFIT1-3, RSAD2 and MX1). CD74 and MHC-II genes were highly expressed in the three aged clusters (Fig. 7c), reflecting a previously described non-classical antigen-presenting phenotype ⁷⁹.

We calculated the previously described N1 (anti-tumor) and N2 (pro-tumor) signatures for all neutrophil subsets ⁷⁴. Nearly all clusters were dominated by the N2-phenotype. The aged neutrophils (G1, G2 and G5) showed slight enrichment in the N1-phenotype (Fig. 7d). We found that the expression of the ELANE (elastase) gene was negative in all neutrophils (Extended Data Fig. 7b). ELANE has been reported to have an anti-cancer function in human neutrophils ⁸⁰. Taken together, these results suggested that the neutrophils in the TME are dominated by immunosuppressive phenotypes, with aged neutrophils showed moderate anti-tumor and antigen-presenting functions.

To explore the reprogramming of neutrophils after combined therapy, we compared abundance of neutrophil subtype fractions before and after treatment (Fig. 7e). Most mature neutrophils decreased after treatment (G0 and G3), while aged neutrophils increased (G1, G2 and G5). This suggests that therapy induces neutrophils aging (Fig. 7e).

Finally, we explored the mechanisms of neutrophil subsets reprogramming during therapy. NicheNet analysis predicted SPP1, IFNγ and TNF ligands to drive the phenotype of aged G1-CCL3 cells. NicheNet analysis also predicted that IFNα ligand to drive the phenotype of aged G2-IFIT3 cells (Fig. 7f and Extended Data Fig. 7c). SPP1 was a key marker of C1-SPP1 macrophages. We evaluated the cell-cell interactions between aged G1-CCL3 neutrophils and C1-SPP1 macrophages. The SPP1-PTGER4 and SPP1-CD44 pairs were significantly enriched in the two cell-types (Fig. 7g). These receptor-ligand pairs are reported to mediate neutrophil phagocytosis and IL-8 production ^81,82. CCL3 and CCL4 secreted by aged G1-CCL3 cells were predicted to interact with CCR1 and CCR5 in C1-SPP1 macrophages, respectively. These indicates that C1-SPP1 macrophages drive the production of aged G1-CCL3 neutrophils, and that aged G1-CCL3 neutrophils in turn attract C1-SPP1 macrophages. This interaction may promote the increase of C1-SPP1 macrophages and G1-CCL3 neutrophils in NMPR patients (Fig. 5c and 7e).

Improving the response efficacy and identifying robust biomarkers are the major challenges for current immunotherapy. Although ICB therapy has been used in advanced NSCLC for years, many patients are refractory to the treatment. The transcriptional characteristics underlying ICB resistance have not been characterized due to the difficulty of sample acquisition. The advent of neoadjuvant immunotherapy for resectable NSCLC provides the opportunity to collect tumor tissues before and after treatment and pathologic assessment of resected tumor tissues enables more precise response information compared to traditional RECIST classifications. In this study, we examined single-cell transcriptomes from resectable NSCLC before and after combination treatment of PD-1 blockade and chemotherapy, and analyzed the TME across pathologic responses to investigate immune system and cancer responses to therapy.

We uncovered transcriptional signatures of cancer cells specific to different pathologic responses. MHC-II genes were highly expressed in malignant cells from MPR patients. The important role of antigen presentation via the MHC-II pathway in MPR patients was also demonstrated by CD20+ B cells and cDC2s presenting tumor antigens to CD4+ T cells by MHC-II. Further, CD4+ T cells were activated by CD40-CD40LG and CD28-CD86 interactions in MPR patients. This suggests that although MHC-II expression is usually restricted to APCs, however, MHC-II could also be expressed intrinsically in a subset of cancer cells ⁸³ or induced by IFNγ ²⁹. Given the low expression of MHC-II in TN patients, it is likely that the expression of MHC-II genes was induced by IFNγ secreted by effector T or NK cells as a result of therapy. Recent studies report that inhibition of histone deacetylases (HDAC) and mitogen-activated protein kinase kinase (MEK) enhanced MHC-II expression in NSCLC cell lines ⁸³. Thus, indicated that promoting antigen presentation via the MHC-II pathway may be a strategy to enhance response to immunotherapy. The upregulation of SERPINB9 after therapy suggests that anti-SERPINB9 drugs combined with ICB may a better adjuvant therapy regimen than ICB alone for MPR patients. Compound 3034, a SERPINB9 inhibitor, was reported to reduce melanoma and breast tumor growth in mice ³². However, the safety and efficacy of compound 3034 in combinatorial therapy remains to be established clinically.

In NMPR patients, we observed an overexpression of enzymes involved in estrogen metabolism in cancer cells, which resulted in the elevation of estradiol in serum. Estradiol led to an immunosuppressive environment, EGFR signaling activation, and upregulation of VEGFA in cancer cells. Previous studies have indirectly suggested a negative effect for estradiol in immunotherapy. A meta-analysis including 11,351 patients treated with ICB from 20 trials showed that ICB was significantly less effective in females males, and that females had no significant survival benefits in NSCLC ⁸⁴. Therefore, change of estradiol during therapy may be a biomarker for a poor response to therapy, and raises the possibility that a regimen that combined anti-estrogen with ICB may improve response. Indeed, anti-estrogen had been explored in EGFR mutation positive NSCLC, although not in combination with ICB. Phased II trials failed to show a survival benefit for anti-estrogen combined with targeted therapy ^85,86, and therefore a phase III trial was not pursued. Future studies are needed to explore the efficacy of anti-estrogen combined with ICB.

Patients with different pathological responses showed remarkable differences in the TME remodeling after therapy. We hypothesize that there may be a “positive feedback” immune response in good responders (MPR patients), and a “negative feedback” response in poor responders (NMPR patients). ICB activates CD8+ T cells and NK cells to release multiple cytokines, thereby inducing expression of MHC-II in cancer cells and activating B cells and DCs to present tumor antigens. The reprogrammed cancer cells then present the tumor antigens to CD4+ T cells via MHC-II and release CX3CL1 and SPARCL1 to recruit NK cells, PMos and M1 macrophages. In good responders, simultaneous with the recruitment and activation of cytotoxic cells into the TME, immunosuppressive cells (Tregs, neutrophils and M2 macrophages) are significantly reduced in the TME. In poor responders, therapy triggers the anti-tumor immunity and cytotoxic cells are activated at the beginning, but the immunosuppressive cells were not adequately remodeled and become more abundant in the TME. In poor responders, cancer cells limit the immune response by increasing estradiol levels. Adjunctive therapies that decrease the immunosuppressive cells or factors in the TME to allow a return to a positive feedback cycle may greatly improve response to immunotherapy.

In summary, scRNA-seq analysis of resectable NSCLC after neoadjuvant ICB combined with chemotherapy revealed that residual cancer cells and TME of good responders had properties that were distinct from poor responders. We identified two cell-types in the TME (FCRL4+ FCRL5+ memory B cells and CD16+CX3CR1+ monocytes) that could serve as biomarkers for treatment response. The data indicate that maintenance therapy is necessary for pCR patients, that a regimen of SERPINB9 inhibitor plus ICB may be effective for MPR patients after surgery, and that a new anti-estrogen plus ICB strategy improve response in NMPR patients.

Patient cohorts and sample collection

Patient data from three prospective clinical trials of the anti-PD-1 antibody plus chemotherapy in resectable non-small cell lung cancer (NSCLC) were used in this study: LungMate-001 (ChiCTR1900023758, http://www.chictr.org.cn/showproj.aspx?proj=39917), LungMate-002 (ChiCTR1900024014, http://www.chictr.org.cn/showproj.aspx?proj=40326) and LungMate-003 (NCT04379739, https://clinicaltrials.gov/ct2/show/NCT04379739). All the three investigator-initiated phase II trials were single-center and open-label, which aimed to explore the efficacy of neoadjuvant anti-PD-1 immunotherapy plus chemotherapy for NSCLC. Patients with stage IB-III NSCLC enrolled in the three trials were treatment-naïve and EGFR/ALK mutation negative. Institutional Review Board approval and informed consent were obtained prior to tissue acquisition and genomic sequencing for each patient. We obtained the primary tumor tissue by percutaneous pulmonary biopsy, bronchoscopy biopsy or endobronchial ultrasound (EBUS) biopsy before drug administration. The patients received 2 cycles of neoadjuvant therapy (PD-1 antibody + platinum-based chemotherapy), and then underwent surgery. The fresh tumor tissues were collected immediately after surgical resection.

Tissue dissociation and cell purification

Samples from biopsy or surgery were isolated and transported rapidly to the research facility. Tissues were transported in a sterile culture dish with 10 ml 1x Dulbecco's Phosphate-Buffered

Saline (DPBS; Thermo Fisher, Cat. no. 14190144) on ice (4 ℃) to remove the residual tissue storage solution, then minced into 1-3 mm³ pieces in another culture dish. We used 10 mg type Ⅰ collagenase (Sigma, Cat. no. C0130) dissolved in 10 ml RPMI 1640 medium (Thermo Fisher, Cat. no. 10-040-CM) with 10% Fetal Bovine Serum (FBS; Thermo Fisher, Cat. no. SV30087.02) to digest the tissues. Tissues were dissociated at 37 ℃ with a shaking speed of 50 r.p.m. Cell suspensions were filtered using a 70-um nylon cell strainer and red blood cells were removed by 1X Red Blood Cell Lysis Solution (Thermo Fisher, Cat. no. 00-4333-57). Dissociated cells were washed with 1x DPBS containing 2% FBS. Cells were stained with 0.4% Trypan blue (Thermo Fisher, Cat. no. 14190144) to check the viability on Countess® II Automated Cell Counter (Thermo Fisher).

Single cell RNA sequencing

To capture single-cell transcriptomic information of lung cancer samples, we used the BD Rhapsody Single-Cell Analysis System (BD Biosciences) according to the manufacturer’s protocol. Single-cell capture was achieved by the random distribution of a single cell suspension across ~200,000 microwells. Beads with unique molecular identifiers (UMIs) and cell barcodes were loaded close to saturation, so that each cell was paired with a bead in a microwell. After exposure to cell lysis buffer, polyadenylated RNA molecules hybridized to the beads. Beads were retrieved into a single tube for reverse transcription. On cDNA synthesis, each cDNA molecule was tagged on the 5’ end (that is, the 3’ end of a messenger RNA transcript) with UMI and cell label indicating its cell of origin. Whole transcriptome libraries were prepared using the BD Resolve single-cell whole transcriptome amplification workflow. Briefly, Rhapsody beads were then subject to second-strand cDNA synthesis, adaptor ligation, and universal amplification. Sequencing libraries were prepared using random priming PCR of the whole-transcriptome amplification products to enrich the 3’end of the transcripts linked with the cell label and UMI. Sequencing libraries were quantified using a High Sensitivity DNA Chip (Agilent) on a Bioanalyzer 2100 and the Qubit High Sensitivity DNA Assay (Thermo Fisher Scientific). The libraries were sequenced on NovaSeq6000 (Illumina) using 2x150 chemistry. The BD Rhapsody analysis pipeline was used to process raw sequencing data (FASTQ files).

Bulk RNA sequencing

RNA was isolated from fresh frozen tissues and perform RNA-seq. Total RNA was isolated with the RNeasy Mini Kit (Qiagen). The NEBNext Ultra RNA library (New England Biolabs) was used to construct the RNA-seq libraries according to the manufacturer’s protocol. Then, the quality-checked libraries were sequenced on the Illumina Novaseq 6000 platform.

Untargeted metabolomics

Metabolites Extraction

The samples (100 μL) were placed in the EP tubes and resuspended with prechilled 80% methanol and 0.1% formic acid by a well vortex. Then the samples were incubated on ice for 5 min and centrifuged at 15,000 g, 4°C for 20 min. Some of the supernatant was diluted to final concentration containing 53% methanol by LC-MS grade water. The samples were subsequently transferred to a fresh Eppendorf tube and then were centrifuged at 15000 g, 4°C for 20 min. Finally, the supernatant was injected into the LC-MS/MS system analysis.

UHPLC-MS/MS Analysis

UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (ThermoFisher, Germany) coupled with an Orbitrap Q Exactive^TM HF mass spectrometer (Thermo Fisher, Germany) in Novogene Co., Ltd. (Beijing, China). Samples were injected onto a Hypesil Goldcolumn (100×2.1 mm, 1.9μm) using a 17-min linear gradient at a flow rate of 0.2 mL/min. The eluents for the positive polarity mode were eluent A (0.1% FA in Water) and eluent B (Methanol). The eluents for the negative polarity mode were eluent A (5 mM ammonium acetate, pH 9.0) and eluent B (Methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2-100% B, 12.0 min; 100% B, 14.0 min; 100-2% B, 14.1 min; 2% B, 17 min. Q Exactive^TM HF mass spectrometer was operated in positive/negative polarity mode with a spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas flow rate of 40 arb and aux gas flow rate of 10 arb.

Multiplex immunohistochemistry

Multiplex immunofluorescence staining was performed using PANO 4-plex IHC kit (cat 10001100100, Panovue). We performed the fluorescent dyes by using the CD20 mouse anti-human antibody (Dako, M0755), CD86 rabbit anti-human antibody (CST, 91882S), and FCRL4 rabbit anti-human antibody (Abcam, ab239076). Different above primary antibodies were applied, followed by horseradish peroxidase-conjugated secondary antibody incubation and tyramide signal amplification. The slides were microwave heat-treated after each TSA operation. Nuclei were stained with DAPI (SIGMA-ALDRICH, D9542) after all the human antigens had been labelled. To obtain multispectral images, the stained slides were scanned using the Mantra System (PerkinElmer, Waltham, Massachusetts, US), which captures the fluorescent spectra at 20-nm wavelength intervals from 420 to 720 nm with identical exposure time; the scans were combined to build a single stack image.

Single cell RNA-seq data generation and quality control

scRNA-seq FASTQ files were processed using the BD Rhapsody Whole Transcriptome Analysis (WTA) Pipeline to get a unique molecular identifier (UMI) matrix for each sample. The matrix of read counts per gene per sample was further analyzed by the Seurat package (version 3.2.2) ⁸⁷in the R software (version 3.6.3). For each cell, we used four quality control (QC) measures. Cells meeting any of the following criteria were excluded: (1) < 500 expressed genes, (2) > 20% UMIs of mitochondria genes, (3) > 50% UMIs of ribosome genes and (4) housekeeping score (defined as the sum of the UMIs of three canonical housekeeping genes: ACTB, GAPDH and MALAT1) < 1. To exclude data from droplets containing more than one cell, doublet detection and removal were performed using Scrublet ⁸⁸. An expected doublet rate parameter of 0.025 was used, and doublet score thresholds were chosen manually to divide putative singlet and neotypic doublet modes in the score distribution. Predicted doublets were then removed from gene-barcode matrices.

Across-sample integration

The gene expression matrices were normalized by the NormalizeData function with default parameters. To adjust for biological sources of variation between samples, the standard anchor-based workflow for dataset integration in Seurat was used. As a previous study ⁴⁶, 3000 or 4000 (for B cells and neutrophils clustering) variable features for CCA (canonical correlation analysis) ⁸⁹ were chosen based on stabilized variance, and integration anchors were identified using the first 20 reduced dimensions. Integration-transformed expression values were used only for dimension reduction and clustering. The original log-normalized expression values were used for all differential expression and gene set level analyses.

Dimension reduction and unsupervised clustering

Principal components analysis (PCA) was performed on the integration-transformed expression matrix using the RunPCA function, and the first 15 principal components (PCs) were used in the FindNeighbors function. The resolution parameters of the FindClusters function were different for different cell types, with 0.6 for all cells, 0.8 for all immune cells, 0.4 for T and myeloid cells and 0.3 for epithelia, B cells and neutrophils. Uniform manifold approximation and projection (UMAP) was performed for visualization in two dimensions using the RunUMAP function with the same PCs and other default parameters. Major cell lineages were assigned to each cluster of cells using the abundance of canonical marker genes (Table S1), and marker genes for each cluster were found using the FindAllMarkers function with the parameter “min.pct = 0.25, thresh.use = 0.25”. Notably, when did cell clustering, we manually removed the clusters that expressed two or more major lineage markers (such as LYZ for myeloid cells and CD3E for T cells) on UMAP plot and probably were doublets that were nor recognized by Scrublet.

CNV estimation and identification of malignant cells

To identify malignant cells from epithelia, we used the CopyKAT algorithm ²⁰ to estimate the copy number variations (CNVs). The stromal cells (fibroblasts and endothelia) were used as normal reference, and the parameters were default. The sum of calculated CNV for each gene per cell was defined as the CNV score of the cell.

Differential expression analysis and gene set variation analysis (GSVA)

Differential expression analysis comparing cells from treatment exposure or response groups was performed using the FindAllMarkers function with the parameter “min.pct = 0.25, thresh.use = 0.25”.

To assign pathway activity estimates to individual cells, we applied GSVA using standard settings, as implemented in the GSVA R package (version 1.32.0), as described previously ²¹. The gene set of 50 hallmark pathways we investigated (h.all.v7.2.symbols.gmt) was downloaded from the Gene Set Enrichment Analysis (GSEA) website (https://www.gsea-msigdb.org/gsea/index.jsp). The differential activities of pathways between groups were calculated using limma R package (version 3.42.2). Significantly disturbed pathways were identified with Benjamini-Hochberg–corrected P value of ≤ 0.01.

Gene module enrichment analysis

To estimate the signature of MHC-II antigen presentation and EGFR signaling in cancer cells, we calculated the enrichment scores for each cell using the AddModuleScore function in Seurat with the gene list from the REACTOME_MHC_CLASS_II_ANTIGEN_PRESENTATION and REACTOME_SIGNALING_BY_EGFR_IN_CANCER pathway (c2.cp.reactome.v7.2.symbols.gmt, download from https://www. gsea-msigdb.org/gsea/index.jsp). To explore the cytotoxic and exhausted functions of T and NK cells, we calculated the cytotoxic score and exhausted score for each cell using the canonical cytotoxic (GZMA, GZMB, GZMK, GNLY, IFNG, PRF1 and NKG7) and exhausted (LAG3, TIGIT, PCCD1, HAVCR2, CTLA4, LAYN and ENTPD1) markers, respectively. To explore the N1 and N2 phenotype of neutrophils, we summarized the marker genes (Table S6) from two review articles ^74,90, and calculated the module score for each cell. With the same method, we used the gene list (Table S7) in “LM22.xls” from CIBERSORT ⁶⁸ to estimate the phenotype (M0, M1 or M2) of macrophage for each cell. We also calculated the antigen presentation score and immunosuppressive score for DCs with the previously reported markers (Table S7) ⁷².

Cellular fraction calculation

Due to few samples in TN (n=3) and MPR (n=4), we did not calculate the cellular fraction for each sample but the total cells in each group as previous studies ^67,91. Error bars indicate the 95% confidence interval for the calculated relative frequencies. For the fraction of B cell subsets, we removed the cells from Patient14 when calculating, because cells from Patient14 dominated cluster B2-HSPA1A, which had a strong impact on the fraction of B2-HSPA1A in MPR patients. The significance of differences among response groups for the fraction of each cell type was compared using the chi-square test.

Trajectory analysis

RNA velocity

The bam files from the WTA Pipeline were converted into sam files first using samtools (version 1.7). Then the cell barcode tag “MA” for each sequence was replaced by “UB” which could be recognized by the velocyto run function from RNA velocity (version 0.17.17) algorithm ⁵⁰. We removed the sequences without cell barcode, and the sam files were converted back into bam files and sorted using samtools. The sorted bam files were used to generate loom files using the velocyto run function with the genome annotation file “GRCh38.gtf”. The loom files for each sample were merged into one loom file. To calculate the velocity and visualize on plot, we used the scVelo ⁹² method using steady-state mode, following the Seurat to RNA-Velocity guides (https://github.com/basilkhuder/Seurat-to-RNA-Velocity).

Monocle2

We also used Monocle2 ⁵² to infer the cell lineage trajectory of T cells and myeloid cells with the top 400 signature genes calculated by differentialGeneTest function. The differentiation trajectory was inferred with the default parameters of Monocle after dimension reduction and cell ordering.

Cell-cell interaction analysis

We used CellPhoneDB ⁴⁷ to infer cell-cell interaction between different cell types. This method infers the potential interaction strength between two cell subsets based on gene expression level, and provides the significance through permutation test (1000 times). The enriched ligand-receptor interactions between two cell subsets were calculated based on the permutation test. We extracted significant ligand-receptor pairs with P-value < 0.01.

NicheNet analysis

NicheNet ⁶¹, a powerful tool that predicts ligands driving the transcriptomic changes of target cells, was used to identify potential ligands that drive the unique phenotype of B cell and neutrophil subsets. As described previously ⁷¹, we used all expressed genes of the B4-FCRL4 cells as the background of genes and the top 50 genes ordered by log2FC as gene sets of interest. Genes were considered as expressed when they have nonzero values in at least 10% of the cells in a cell type. We only used the expressed receptors in B4-FCRL4 cells to construct the expressed ligand-receptor interactions and calculate the ligand activity using the predict_ligand_activities function. For the two subsets of neutrophils, the same method was performed.

Bulk RNA-seq data processing and quantification

Raw FASTQ files were aligned on the hg38 genome reference using the STAR aligner with default parameters. Salmon was used to quantitate gene expression by transcripts per kilobase million (TPM). The TPM matrix was transformed by log2(TPM+1).

Cell signature estimation in bulk RNA-seq

The cell signatures of B4-FCRL4 and C7-CX3CR1 were estimated by ssGSEA method in the GSVA R package (version 1.32.0). The ssGSEA transforms specific gene expression patterns into quantities of cell populations in individual tumor samples at bulk level. The markers of the two cell types (Table S5-6) were used as the gene sets in the gsva function with the parameters “method=’ssgsea’, kcdf='Gaussian', abs.ranking=TRUE”. For melanoma dataset #2 ⁵⁹, we analyzed the prognostic value of the two cell signatures. The patients were divided into high and low signature groups by the median value, and the Kaplan–Meier survival curves with the cumulative number of events table and log-rank test were plotted by survminer (version 0.4.8) and survival (version 3.1-8) R package.

package.

TCGA analysis

The analyses of lung adenocarcinoma (LUAD) from the TCGA database were performed on the TCGA visualization web server, GEPIA2 (http://gepia2.cancer-pku.cn/), developed by Zeming Zhang ⁹³. The server provided several function modules, and we used the “Survival Analysis” module to explore the correlation between the expression of genes of interest and overall survival. The patients were divided into high and low signature groups and hazards ratio (HR) was calculated.

Metabolomic data processing and metabolite identification

The raw data files generated by UHPLC-MS/MS were processed using the Compound Discoverer 3.1 (CD3.1, ThermoFisher) to perform peak alignment, peak picking, and quantitation for each metabolite. The main parameters were set as follows: retention time tolerance, 0.2 minutes; actual mass tolerance, 5ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, et al. After that, peak intensities were normalized to the total spectral intensity. The normalized data was used to predict the molecular formula based on additive ions, molecular ion peaks and fragment ions. Then peaks were matched with the mzCloud (https://www.mzcloud.org/), mzVault and MassList database to obtain the accurate qualitative and relative quantitative results. When data were not normally distributed, normal transformations were attempted using of area normalization method.

Drug sensitivity analysis

The gene expression data (RPKM matrix) of NSCLC cell lines was downloaded from the (CCLE) database (https://portals.broadinstitute.org/ccle/), and the IC50 data of drugs was downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerRxgene.org). The gene expression of 65 NSCLC cell lines and the IC50 data of 16 drugs were used to analyze in this study. We first estimated the signature of the 5 AKR family genes (AKR1C1-3 and AKR1B1/10) for each cell line using ssGSEA method, then we calculated the pearson correlation coefficient between the signature and IC50 value for each drug.

Statistics

The cell signature or gene expression comparisons between two groups were performed using two-tailed Wilcoxon rank-sum tests. All statistical analyses and presentations were performed using R software (version 3.6.3). Statistical significance was set at P < 0.05.

Data and code availability

Raw single-cell RNA-seq and bulk RNA-seq data are deposited to Genome Sequence Archive of the BIG Data Center at the Beijing Institute of Genomics, Chinese Academy of Science, under accession number HRA001033 (accessible at http://bigd.big.ac.cn/gsa-human). The count matrix of scRNA-seq and TPM matrix of bulk RNA-seq can be obtained from Gene Expression Omnibus (GEO:) and Table S3, respectively. Two Melanoma treatment-related datasets ^59,60 are available from https://static-content.springer.com/esm/art%3A10.1038%2Fs41591-019-0654-5/MediaObjects/41591_2019_654_MOESM3_ESM.txt and GEO: GSE91061, respectively. Code used for all processing and analysis is available upon request.

Acknowledgements

We thank the patients who participated in this study. We thank X.Shirley Liu for giving helpful advice. This work was supported by the National Natural Science Foundation of China (Grant No. 81972172 and 81872290), Shanghai Science and Technology Committee (Grant No. 19XD1423200), Shanghai Hospital Development Center (Grant No. SHDC2020CR2020B), Shanghai Rising-Star Program (21QA1408200), The General Program Fund of Shanghai Natural Science Foundation (21ZR1467600), and Shanghai Pulmonary Hospital (Grant No. FKCX1904 and FKLY20004).

Author contributions

P.Z. and C.W. designed the study, X.Z., F.S., and H.X. cared for patients analyzed in the study and handled the patient samples, J.H., L.Z., and H.Y. processed scRNA-seq raw data and Y.Y. and L.S. processed bulk RNA-seq raw data. X.Z., F.S., H.X., J.Z., and D.B. collected and analyzed clinical metadata for all patient samples. J.H., L.Z., L.S., H.X., J.W. and Y.H., performed scRNA-seq analysis, and Y.Y., S.L., and D.L. performed bulk RNA-seq analysis. L.Z. and J.Z. analyzed H&E images. Y.C., P.F., and Q.M. provided help for statistic tests. G.J. provided administrative support and give critical advice. J.H., L.Z., and Y.Y., wrote the manuscript with input from all authors, and C.W. and P.Z. revised the manuscript. All authors read and approved the manuscript.

Competing interests

The authors declare no potential conflicts of interest.

1 Sung, H. et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin, doi:10.3322/caac.21660 (2021).

2 Goldstraw, P. et al. Non-small-cell lung cancer. Lancet (London, England) 378, 1727-1740, doi:10.1016/s0140-6736(10)62101-0 (2011).

3 Mok, T. S. K. et al. Pembrolizumab versus chemotherapy for previously untreated, PD-L1-expressing, locally advanced or metastatic non-small-cell lung cancer (KEYNOTE-042): a randomised, open-label, controlled, phase 3 trial. Lancet (London, England) 393, 1819-1830, doi:10.1016/s0140-6736(18)32409-7 (2019).

4 Forde, P. M. et al. Neoadjuvant PD-1 Blockade in Resectable Lung Cancer. The New England journal of medicine 378, 1976-1986, doi:10.1056/NEJMoa1716078 (2018).

5 Travis, W. D. et al. IASLC Multidisciplinary Recommendations for Pathologic Assessment of Lung Cancer Resection Specimens After Neoadjuvant Therapy. Journal of thoracic oncology : official publication of the International Association for the Study of Lung Cancer 15, 709-740, doi:10.1016/j.jtho.2020.01.005 (2020).

6 Jia, X. H. et al. Efficacy and safety of neoadjuvant immunotherapy in resectable nonsmall cell lung cancer: A meta-analysis. Lung Cancer 147, 143-153, doi:10.1016/j.lungcan.2020.07.001 (2020).

7 Altorki, N. K. et al. The lung microenvironment: an important regulator of tumour growth and metastasis. Nature reviews. Cancer 19, 9-31, doi:10.1038/s41568-018-0081-9 (2019).

8 Binnewies, M. et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nature medicine 24, 541-550, doi:10.1038/s41591-018-0014-x (2018).

9 Yost, K. E. et al. Clonal replacement of tumor-specific T cells following PD-1 blockade. Nature medicine 25, 1251-1259, doi:10.1038/s41591-019-0522-3 (2019).

10 Wu, T. D. et al. Peripheral T cell expansion predicts tumour infiltration and clinical response. Nature 579, 274-278, doi:10.1038/s41586-020-2056-8 (2020).

11 Sade-Feldman, M. et al. Defining T Cell States Associated with Response to Checkpoint Immunotherapy in Melanoma. Cell 175, 998-1013.e1020, doi:10.1016/j.cell.2018.10.038 (2018).

12 Gide, T. N. et al. Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer cell 35, 238-255.e236, doi:10.1016/j.ccell.2019.01.003 (2019).

13 Helmink, B. A. et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature 577, 549-555, doi:10.1038/s41586-019-1922-8 (2020).

14 Goswami, S. et al. Immune profiling of human tumors identifies CD73 as a combinatorial target in glioblastoma. Nature medicine 26, 39-46, doi:10.1038/s41591-019-0694-x (2020).

15 Caushi, J. X. et al. Transcriptional programs of neoantigen-specific TIL in anti-PD-1-treated lung cancers. Nature 596, 126-132, doi:10.1038/s41586-021-03752-4 (2021).

16 Gaudreau, P. O. et al. Neoadjuvant Chemotherapy Increases Cytotoxic T Cell, Tissue Resident Memory T Cell, and B Cell Infiltration in Resectable NSCLC. Journal of thoracic oncology : official publication of the International Association for the Study of Lung Cancer 16, 127-139, doi:10.1016/j.jtho.2020.09.027 (2021).

17 DeNardo, D. G. & Ruffell, B. Macrophages as regulators of tumour immunity and immunotherapy. Nature reviews. Immunology 19, 369-382, doi:10.1038/s41577-019-0127-6 (2019).

18 Teijeira, Á. et al. CXCR1 and CXCR2 Chemokine Receptor Agonists Produced by Tumors Induce Neutrophil Extracellular Traps that Interfere with Immune Cytotoxicity. Immunity 52, 856-871.e858, doi:10.1016/j.immuni.2020.03.001 (2020).

19 Reddy, S. M. et al. Poor Response to Neoadjuvant Chemotherapy Correlates with Mast Cell Infiltration in Inflammatory Breast Cancer. Cancer immunology research 7, 1025-1035, doi:10.1158/2326-6066.Cir-18-0619 (2019).

20 Gao, R. et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nature biotechnology, doi:10.1038/s41587-020-00795-2 (2021).

21 Xing, X. et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Science advances 7, doi:10.1126/sciadv.abd9738 (2021).

22 Nabhan, A. N., Brownfield, D. G., Harbury, P. B., Krasnow, M. A. & Desai, T. J. Single-cell Wnt signaling niches maintain stemness of alveolar type 2 cells. Science (New York, N.Y.) 359, 1118-1123, doi:10.1126/science.aam6603 (2018).

23 Shu, C. A. et al. Neoadjuvant atezolizumab and chemotherapy in patients with resectable non-small-cell lung cancer: an open-label, multicentre, single-arm, phase 2 trial. The Lancet. Oncology 21, 786-795, doi:10.1016/s1470-2045(20)30140-6 (2020).

24 Provencio, M. et al. Neoadjuvant chemotherapy and nivolumab in resectable non-small-cell lung cancer (NADIM): an open-label, multicentre, single-arm, phase 2 trial. The Lancet. Oncology 21, 1413-1422, doi:10.1016/s1470-2045(20)30453-8 (2020).

25 Liang, W. et al. Expert consensus on neoadjuvant immunotherapy for non-small cell lung cancer. Translational lung cancer research 9, 2696-2715, doi:10.21037/tlcr-2020-63 (2020).

26 Zhao, S. J. et al. SPARCL1 suppresses osteosarcoma metastasis and recruits macrophages by activation of canonical WNT/β-catenin signaling through stabilization of the WNT-receptor complex. Oncogene 37, 1049-1061, doi:10.1038/onc.2017.403 (2018).

27 Cózar, B. et al. Tumor-Infiltrating Natural Killer Cells. Cancer discovery, doi:10.1158/2159-8290.Cd-20-0655 (2020).

28 Hanna, R. N. et al. Patrolling monocytes control tumor metastasis to the lung. Science (New York, N.Y.) 350, 985-990, doi:10.1126/science.aac9407 (2015).

29 Axelrod, M. L., Cook, R. S., Johnson, D. B. & Balko, J. M. Biological Consequences of MHC-II Expression by Tumor Cells in Cancer. Clinical cancer research : an official journal of the American Association for Cancer Research 25, 2392-2402, doi:10.1158/1078-0432.Ccr-18-3200 (2019).

30 Alspach, E. et al. MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature 574, 696-701, doi:10.1038/s41586-019-1671-8 (2019).

31 Johnson, D. B. et al. Melanoma-specific MHC-II expression represents a tumour-autonomous phenotype and predicts response to anti-PD-1/PD-L1 therapy. Nature communications 7, 10582, doi:10.1038/ncomms10582 (2016).

32 Jiang, L. et al. Direct Tumor Killing and Immunotherapy through Anti-SerpinB9 Therapy. Cell 183, 1219-1233.e1218, doi:10.1016/j.cell.2020.10.045 (2020).

33 van Weverwijk, A. et al. Metabolic adaptability in metastatic breast cancer by AKR1B10-dependent balancing of glycolysis and fatty acid oxidation. Nature communications 10, 2698, doi:10.1038/s41467-019-10592-4 (2019).

34 Cheng, B. Y. et al. IRAK1 Augments Cancer Stemness and Drug Resistance via the AP-1/AKR1B10 Signaling Cascade in Hepatocellular Carcinoma. Cancer research 78, 2332-2342, doi:10.1158/0008-5472.Can-17-2445 (2018).

35 Wu, X. et al. AKR1B1 promotes basal-like breast cancer progression by a positive feedback loop that activates the EMT program. The Journal of experimental medicine 214, 1065-1079, doi:10.1084/jem.20160903 (2017).

36 Jin, Y. et al. Human cytosolic hydroxysteroid dehydrogenases of the aldo-ketoreductase superfamily catalyze reduction of conjugated steroids: implications for phase I and phase II steroid hormone metabolism. The Journal of biological chemistry 284, 10013-10022, doi:10.1074/jbc.M809465200 (2009).

37 Velez, M. A., Burns, T. F. & Stabile, L. P. The estrogen pathway as a modulator of response to immunotherapy. Immunotherapy 11, 1161-1176, doi:10.2217/imt-2019-0024 (2019).

38 Svensson, S. et al. CCL2 and CCL5 Are Novel Therapeutic Targets for Estrogen-Dependent Breast Cancer. Clinical cancer research : an official journal of the American Association for Cancer Research 21, 3794-3805, doi:10.1158/1078-0432.Ccr-15-0204 (2015).

39 Svoronos, N. et al. Tumor Cell-Independent Estrogen Signaling Drives Disease Progression through Mobilization of Myeloid-Derived Suppressor Cells. Cancer discovery 7, 72-85, doi:10.1158/2159-8290.Cd-16-0502 (2017).

40 Polanczyk, M. J. et al. Cutting edge: estrogen drives expansion of the CD4+CD25+ regulatory T cell compartment. Journal of immunology (Baltimore, Md. : 1950) 173, 2227-2230, doi:10.4049/jimmunol.173.4.2227 (2004).

41 Siegfried, J. M., Gubish, C. T., Rothstein, M. E., Henry, C. & Stabile, L. P. Combining the multitargeted tyrosine kinase inhibitor vandetanib with the antiestrogen fulvestrant enhances its antitumor effect in non-small cell lung cancer. Journal of thoracic oncology : official publication of the International Association for the Study of Lung Cancer 7, 485-495, doi:10.1097/JTO.0b013e31824177ea (2012).

42 Qiao, M. et al. Immune checkpoint inhibitors in EGFR-mutated non-small cell lung cancer: Dusk or Dawn? Journal of thoracic oncology : official publication of the International Association for the Study of Lung Cancer, doi:10.1016/j.jtho.2021.04.003 (2021).

43 Fukumura, D., Kloepper, J., Amoozgar, Z., Duda, D. G. & Jain, R. K. Enhancing cancer immunotherapy using antiangiogenics: opportunities and challenges. Nature reviews. Clinical oncology 15, 325-340, doi:10.1038/nrclinonc.2018.29 (2018).

44 Shi, Y. E. et al. Synuclein gamma stimulates membrane-initiated estrogen signaling by chaperoning estrogen receptor (ER)-alpha36, a variant of ER-alpha. The American journal of pathology 177, 964-973, doi:10.2353/ajpath.2010.100061 (2010).

45 Liu, X. et al. Genome-wide analysis identifies NR4A1 as a key mediator of T cell dysfunction. Nature 567, 525-529, doi:10.1038/s41586-019-0979-8 (2019).

46 Bi, K. et al. Tumor and immune reprogramming during immunotherapy in advanced renal cell carcinoma. Cancer cell, doi:10.1016/j.ccell.2021.02.015 (2021).

47 Efremova, M., Vento-Tormo, M., Teichmann, S. A. & Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nature protocols 15, 1484-1506, doi:10.1038/s41596-020-0292-x (2020).

48 Sautès-Fridman, C., Petitprez, F., Calderaro, J. & Fridman, W. H. Tertiary lymphoid structures in the era of cancer immunotherapy. Nature reviews. Cancer 19, 307-325, doi:10.1038/s41568-019-0144-6 (2019).

49 Guo, X. et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nature medicine 24, 978-985, doi:10.1038/s41591-018-0045-3 (2018).

50 La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498, doi:10.1038/s41586-018-0414-6 (2018).

51 Dammeijer, F. et al. The PD-1/PD-L1-Checkpoint Restrains T cell Immunity in Tumor-Draining Lymph Nodes. Cancer cell 38, 685-700.e688, doi:10.1016/j.ccell.2020.09.001 (2020).

52 Qiu, X. et al. Single-cell mRNA quantification and differential analysis with Census. Nature methods 14, 309-315, doi:10.1038/nmeth.4150 (2017).

53 Lu, Y. et al. Complement Signals Determine Opposite Effects of B Cells in Chemotherapy-Induced Immunity. Cell 180, 1081-1097.e1024, doi:10.1016/j.cell.2020.02.015 (2020).

54 Kim, S. S. et al. B Cells Improve Overall Survival in HPV-Associated Squamous Cell Carcinomas and Are Activated by Radiation and PD-1 Blockade. Clinical cancer research : an official journal of the American Association for Cancer Research 26, 3345-3359, doi:10.1158/1078-0432.Ccr-19-3211 (2020).

55 Li, H. et al. Fc receptor-like 4 and 5 define human atypical memory B cells. International immunology 32, 755-770, doi:10.1093/intimm/dxaa053 (2020).

56 Yeo, L. et al. Expression of FcRL4 defines a pro-inflammatory, RANKL-producing B cell subset in rheumatoid arthritis. Annals of the rheumatic diseases 74, 928-935, doi:10.1136/annrheumdis-2013-204116 (2015).

57 Siewe, B., Nipper, A. J., Sohn, H., Stapleton, J. T. & Landay, A. FcRL4 Expression Identifies a Pro-inflammatory B Cell Subset in Viremic HIV-Infected Subjects. Frontiers in immunology 8, 1339, doi:10.3389/fimmu.2017.01339 (2017).

58 Kreslavsky, T. et al. Essential role for the transcription factor Bhlhe41 in regulating the development, self-renewal and BCR repertoire of B-1a cells. Nature immunology 18, 442-455, doi:10.1038/ni.3694 (2017).

59 Liu, D. et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nature medicine 25, 1916-1927, doi:10.1038/s41591-019-0654-5 (2019).

60 Riaz, N. et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 171, 934-949.e916, doi:10.1016/j.cell.2017.09.028 (2017).

61 Browaeys, R., Saelens, W. & Saeys, Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nature methods 17, 159-162, doi:10.1038/s41592-019-0667-5 (2020).

62 Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nature immunology 20, 163-172, doi:10.1038/s41590-018-0276-y (2019).

63 Allard, B., Panariti, A. & Martin, J. G. Alveolar Macrophages in the Resolution of Inflammation, Tissue Repair, and Tolerance to Infection. Frontiers in immunology 9, 1777, doi:10.3389/fimmu.2018.01777 (2018).

64 Zhang, L. et al. Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer. Cell 181, 442-459.e429, doi:10.1016/j.cell.2020.03.048 (2020).

65 Zhang, Y., Du, W., Chen, Z. & Xiang, C. Upregulation of PD-L1 by SPP1 mediates macrophage polarization and facilitates immune escape in lung adenocarcinoma. Experimental cell research 359, 449-457, doi:10.1016/j.yexcr.2017.08.028 (2017).

66 Barrett, C. W. et al. Selenoprotein P influences colitis-induced tumorigenesis by mediating stemness and oxidative damage. The Journal of clinical investigation 125, 2646-2660, doi:10.1172/jci76099 (2015).

67 Zilionis, R. et al. Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals Conserved Myeloid Populations across Individuals and Species. Immunity 50, 1317-1334.e1310, doi:10.1016/j.immuni.2019.03.009 (2019).

68 Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nature methods 12, 453-457, doi:10.1038/nmeth.3337 (2015).

69 Wculek, S. K. et al. Dendritic cells in cancer immunology and immunotherapy. Nature reviews. Immunology 20, 7-24, doi:10.1038/s41577-019-0210-z (2020).

70 Maier, B. et al. A conserved dendritic-cell regulatory program limits antitumour immunity. Nature 580, 257-262, doi:10.1038/s41586-020-2134-y (2020).

71 Cheng, S. et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell 184, 792-809.e723, doi:10.1016/j.cell.2021.01.010 (2021).

72 Liu, Y. et al. Tumour heterogeneity and intercellular networks of nasopharyngeal carcinoma at single cell resolution. Nature communications 12, 741, doi:10.1038/s41467-021-21043-4 (2021).

73 Block, I. et al. CFP suppresses breast cancer cell growth by TES-mediated upregulation of the transcription factor DDIT3. Oncogene 38, 4560-4573, doi:10.1038/s41388-019-0739-0 (2019).

74 Jaillon, S. et al. Neutrophil diversity and plasticity in tumour progression and therapy. Nature reviews. Cancer 20, 485-503, doi:10.1038/s41568-020-0281-y (2020).

75 Sinha, P., Clements, V. K., Fulton, A. M. & Ostrand-Rosenberg, S. Prostaglandin E2 promotes tumor progression by inducing myeloid-derived suppressor cells. Cancer research 67, 4507-4513, doi:10.1158/0008-5472.Can-06-4174 (2007).

76 Papayannopoulos, V. Neutrophil extracellular traps in immunity and disease. Nature reviews. Immunology 18, 134-147, doi:10.1038/nri.2017.105 (2018).

77 Ntanasis-Stathopoulos, I., Fotiou, D. & Terpos, E. CCL3 Signaling in the Tumor Microenvironment. Advances in experimental medicine and biology 1231, 13-21, doi:10.1007/978-3-030-36667-4_2 (2020).

78 Mukaida, N., Sasaki, S. I. & Baba, T. CCL4 Signaling in the Tumor Microenvironment. Advances in experimental medicine and biology 1231, 23-32, doi:10.1007/978-3-030-36667-4_3 (2020).

79 Lok, L. S. C. et al. Phenotypically distinct neutrophils patrol uninfected human and mouse lymph nodes. Proceedings of the National Academy of Sciences of the United States of America 116, 19083-19089, doi:10.1073/pnas.1905054116 (2019).

80 Cui, C. et al. Neutrophil elastase selectively kills cancer cells and attenuates tumorigenesis. Cell, doi:10.1016/j.cell.2021.04.016 (2021).

81 Lu, C. H. et al. Intermediate Molecular Mass Hyaluronan and CD44 Receptor Interactions Enhance Neutrophil Phagocytosis and IL-8 Production via p38- and ERK1/2-MAPK Signalling Pathways. Inflammation 40, 1782-1793, doi:10.1007/s10753-017-0622-5 (2017).

82 Dey, I., Giembycz, M. A. & Chadee, K. Prostaglandin E(2) couples through EP(4) prostanoid receptors to induce IL-8 production in human colonic epithelial cell lines. British journal of pharmacology 156, 475-485, doi:10.1111/j.1476-5381.2008.00056.x (2009).

83 Neuwelt, A. J. et al. Cancer cell-intrinsic expression of MHC II in lung cancer cell lines is actively restricted by MEK/ERK signaling and epigenetic mechanisms. Journal for immunotherapy of cancer 8, doi:10.1136/jitc-2019-000441 (2020).

84 Conforti, F. et al. Cancer immunotherapy efficacy and patients' sex: a systematic review and meta-analysis. The Lancet. Oncology 19, 737-746, doi:10.1016/s1470-2045(18)30261-4 (2018).

85 Garon, E. B. et al. Randomized phase II study of fulvestrant and erlotinib compared with erlotinib alone in patients with advanced or metastatic non-small cell lung cancer. Lung Cancer 123, 91-98, doi:10.1016/j.lungcan.2018.06.013 (2018).

86 Mazieres, J. et al. Randomized Phase II Trial Evaluating Treatment with EGFR-TKI Associated with Antiestrogen in Women with Nonsquamous Advanced-Stage NSCLC: IFCT-1003 LADIE Trial. Clinical cancer research : an official journal of the American Association for Cancer Research 26, 3172-3181, doi:10.1158/1078-0432.Ccr-19-3056 (2020).

87 Stuart, T. et al. Comprehensive Integration of Single-Cell Data. Cell 177, 1888-1902.e1821, doi:10.1016/j.cell.2019.05.031 (2019).

88 Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell systems 8, 281-291.e289, doi:10.1016/j.cels.2018.11.005 (2019).

89 Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature biotechnology 36, 411-420, doi:10.1038/nbt.4096 (2018).

90 Keeley, T., Costanzo-Garvey, D. L. & Cook, L. M. Unmasking the Many Faces of Tumor-Associated Neutrophils and Macrophages: Considerations for Targeting Innate Immune Cells in Cancer. Trends in cancer 5, 789-798, doi:10.1016/j.trecan.2019.10.013 (2019).

91 Maynard, A. et al. Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing. Cell 182, 1232-1251.e1222, doi:10.1016/j.cell.2020.07.017 (2020).

92 Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nature biotechnology 38, 1408-1414, doi:10.1038/s41587-020-0591-3 (2020).

93 Tang, Z., Kang, B., Li, C., Chen, T. & Zhang, Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic acids research 47, W556-w560, doi:10.1093/nar/gkz430 (2019).

There is NO Competing Interest.

Download PDF

Version 1

posted

You are reading this latest preprint version

Single-Cell RNA Sequencing Reveals Tumor Microenvironment Remodeling after Neoadjuvant Immunotherapy in Non-small Cell Lung Cancer

Status:

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1