Introduction

Lung cancer is the most common cancer among all human tumor types, with more than 1.7 x 10⁶ new cases worldwide each year. According to the Global Cancer Report data, lung adenocarcinoma (LUAD) accounts for most lung cancers (1). The application of adjuvant or neoadjuvant chemotherapy (NCT) has significantly improved the long-term survival of LUAD patients. At present, for most LUADs that need chemotherapy after being assessed, chemotherapy will be used before and after surgery (2). However, chemotherapy drugs are highly toxic and can often become ineffective (3). In addition, continued ineffective chemotherapy will lead to the generation of drug-resistant tumor cell clones (4, 5) and a delay in tumor removal. Almost all cancer patients show inherent or acquired drug resistance, leading to treatment failure and unsatisfactory overall survival. Therefore, to accurately develop therapies that can overcome drug resistance, it is essential to understand the alterations in the tumor microenvironment driven by chemotherapy.

Many studies have increasingly proved the tumor microenvironment (TME) to be an essential source of intratumoral heterogeneity (6). The heterogeneity within the tumor microenvironment (TME) encompasses not only the variations between different tumor cells but also among various stromal and immune cell types. Investigating the dynamic changes in multiple cell populations within the TME of LUAD following chemotherapy may provide crucial insights into overcoming chemotherapy resistance in LUAD. In this study, we demonstrated the changes in the microenvironment of lung adenocarcinoma with chemotherapy. In particular, we focused on the effect of chemotherapy on the metabolic reprogramming of tumor cells, stromal cells, and immune cells.

Formerly, it was generally believed that consuming glucose in TME by cancer cells may promote nutritional competition, a metabolic mechanism of immunosuppression (7). However, recent studies have shown that tumor-infiltrating immune cells rely on glucose for their energy needs and functionality, with immune cells, particularly macrophages, consuming more glucose than malignant cells. The impaired immune cell metabolism in the tumor microenvironment (TME) helps tumor cells escape immunity (8). The internal metabolic changes in the cells drive immune cells and cancer cells to preferentially obtain glucose and glutamine. It is believed that the selective cellular allocation of these nutrients can be used to develop therapeutic and imaging strategies to enhance or monitor the metabolic processes and activities of specific cell populations in TME (9). Metabolic reprogramming in various cell types in the tumor microenvironment after undergoing chemotherapy may be an essential feature that affects chemotherapy. Our research fully demonstrated the metabolic reprogramming landscape of tumor cells, stromal cells, and immune cells before and after chemotherapy.

Materials and methods

Patients

All patients included in this study understood and signed written informed consent, (Approval number: B2019-436). The clinical samples of scRNA-seq came from patients diagnosed with LUAD, of which 4 cases received no treatment before surgery, and 5 cases received chemotherapy (Pemetrexed + Cisplatin). These samples were donated by inpatients in the Department of Thoracic Surgery, Zhongshan Hospital of Fudan University. After the lung adenocarcinoma tissue sample was taken, a small part was cut for paraffin sections, and the remaining tissue was dissociated into a single-cell suspension. 1×10⁶ cells were drawn from the single-cell suspension for single-cell RNA sequencing.

Preparation of single-cell suspensions

For each patient, as described above, we dissociated the lung adenocarcinoma tumor sample into a single-cell suspension and then took 1×10⁶ cells for single-cell RNA sequencing. We used the Tumor Dissociation Kit (Miltenyi Biotec, Gladbach, Germany) to digest tumor tissues with enzymes according to the manufacturer’s instructions. In short, we first cut the lung adenocarcinoma tissue sample into small tissue pieces about 1cm³ with a surgical scalpel. We then transferred these small tissue pieces to the MACS C tube containing 4.7 mL DMEM serum-free medium, 200 µL Enzyme H, 100 µL Enzyme R, and 25 µL Enzyme A. After the tissue was incubated and digested in a constant temperature incubator 37 [for 1 hour, the tissue was mechanically separated by the MACS™ instrument. This procedure was repeated twice. After the tissue sample was dissociated, the sample was filtered with a 40 μm filter to remove the remaining large particles from the single-cell suspension. Centrifuge the suspension at 300 × g for 7 minutes, then discard the supernatant.

Next, we used red blood cell lysate (10×) (Sigma-Aldrich, St. Louis, MO, USA) to remove red blood cells from the single-cell suspension. In short, add 1x Lysis Buffer to the centrifuge tube containing the single-cell pellet described above. The cell suspension was then incubated at room temperature for 15 minutes. To improve the quality of our samples, we also used a Dead Cell Removal Kit (Miltenyi Biotec) to ensure that the cell survival rate of our sequencing samples was >90%.

The 10x scRNA-seq data analysis

The R version used in our scRNA-seq data analysis study is 3.6.1. The cell quality control criteria are as follows: 1) The number of expressed genes is less than 300 or greater than 5000; 2) 10% or more of UMI is localized to mitochondrial or ribosomal genes. If they meet one of the criteria, the cells are excluded. After quality standardization, we applied the Seurat R package (10) to analyze the scRNA-seq data. First, we convert the scRNA-seq data into Seurat identifiable objects, and then we use the “FindVariableFeatures” function to find the first 2000 highly variable genes. After that, we applied principal component analysis (PCA) to reduce the dimensionality of scRNA-seq data. The “RunTSNE” function is used to perform t-distributed random neighborhood embedding (TSNE) to visualize various types of cells. The “FindClusters” and “FindAllMarkers” functions are used for cluster analysis of cell subclusters and detection of marker genes of cell subclusters.

Finally, according to the SingleR package (11), the CellMarker (http://bio-bigdata.hrbmu.edu.cn/CellMarker/) data set, and a previous report (12), we annotated different cell types. Simultaneously, some new potential marker genes were verified through experiments.

Analysis of Sub-Clusters of Cells in LUAD

After preliminary classification and annotation of all cells, epithelial cells, stromal cells, and immune cells are extracted through the “SubsetData” function. Then, we apply the “FindClusters” and “FindAllMarkers” functions to find the marker genes of each cell and perform dimensionality reduction clustering on each extracted cell through TSNE. The sub-clusters are annotated by dominantly expressed cell markers published by previous researchers. To select the marker genes that meet the requirements, we set the following cut-off thresholds to reveal the marker genes of each cluster: adjusted P-value <0.01 and multiple Log2FC >0.5.

Estimation of the copy number variations

To estimate the initial copy number variation (CNV) of each region, the R package “scCancer” (13) was applied. The expression level of each cell was used as the original input file for calculating CNV. Immune cells served as a background reference for calculating the CNVs scores of other cells. In addition, the R package “inferCNV” was used to quantify CNV in tumor cells as described previously (14).

Definition of cell scores and signature

To evaluate the M1/M2 polarization state and pro-/anti-inflammatory potential of macrophages, we performed a GSVA (Gene Set Variation Analysis) analysis. We retrieved gene sets related to the above functions from previous studies (15) and used them as references in this analysis.

We used the average expression of a published list of characteristic genes for T cell toxicity and exhaustion to define T cells’ cytotoxicity, exhaustion, and costimulation scores.

Identification of gene markers of malignant cells

We used the identified malignant cell marker genes in tumor cells to identify gene expression characteristics in malignant cells. Then, we performed unsupervised NMF (Non-negative Matrix Factorization) to reveal the malignant characteristics of tumor cells through the NMF R package (16).

Trajectory analysis

We used the monocle2 R package to analyze the trajectory of all cells to explore the trajectory progression of various types of cells in a single cell (17). First, apply the function “newCellDataSet” to construct a data object that the monocle 2 R package can recognize. Afterward, the differentially expressed genes identified by the Seurat R package were selected for cell trajectory analysis. The “reduceDimension” function was used to reduce the dimensionality. We used the “orderCells” function to project cells on a pseudo-time trajectory to show the trend of cell evolution. A state consisting of cells mainly derived from nonmalignant tissues in a cluster identified as epithelial cells was defined as “root cells.”

Analyses of metabolic pathways

To evaluate the activity of various metabolic pathways of each cell type, we applied the algorithm developed by Xiao et al. (18). In short, the analysis of metabolic programs is based on the average expression level of metabolic genes across cell types to indirectly reflect the metabolic activity of cells.

A variety of environmental factors may potentially affect the metabolic reprogramming of tumors, such as chemotherapy, nutrient supply, and the environment where the cells are located. Therefore, exploring these factors and the cross-conversion between glycolysis and mitochondrial activity in various cells in the tumor microenvironment is essential for understanding the metabolic reprogramming of tumors.

We calculated the average gene expression levels in glycolysis and OXPHOS as indicators of glucose supply and mitochondrial activity, respectively. The data of genes that were responsive to the two groups of genes (known to be responsive to glycolysis and OXPHOS) used in the calculations were retrieved from the MsigDB database. At the same time, the cells were sorted by flow cytometry, and the contents of various metabolites were tested, in turn, to verify whether they were consistent with gene expression levels.

Cell Interaction Network analysis

To study the cell-to-cell interactions between tumors and nonmalignant cells, immune cells, and stromal cells, we applied the R package “CellChat” (19) and “CellPhoneDB” Python package for analysis (20). The crosstalk analysis between cells through the “CellChat” package was as follows: (1) First, use the “createCellChat” function to create a data set object that can be identified by “CellChat”; (2) Then use “aggregateNet”, “computeCommunProbPathway”, and “computeCommunProb” function to automatically infer the possible cellular communication network between cells; (3) Finally, the “netVisual_aggregate”, “netVisual_bubble” and “netVisual_signalingRole” functions were used to visualize the interaction between these cells. Then use the built-in parameters to apply the “CellPhoneDB” R package.

Immunohistochemistry and immunofluorescence

The paraffin-embedded lung cancer tissue sections were deparaffinized with xylene and rehydrated. Discard the blocking solution, add the primary antibody, and incubate overnight at 4 degrees. After removing the primary antibody and washing thoroughly, add the secondary antibody to incubate for 1 hour, and then add DAB chromogenic reagent (Gene Tech, China) for color development. Finally, hematoxylin is used for nuclear dyeing.

As mentioned in the above immunohistochemistry experiment, the steps before incubating the primary antibody are the same. Incubate with the corresponding primary and secondary antibodies with green and red fluorescent dyes, respectively, and then use DAPI to stain the nuclei.

Flow cytometry assay

Cells and APC-conjugated mouse anti-human CD45, FITC-conjugated mouse anti-human CD11b, BV421-conjugated mouse anti-human ARG1, as well as pe-cy-conjugated mouse anti-human CD86 (5 μL/10⁶ cells; BD Biosciences) were incubated on ice for 30 minutes. Then, FACSAria III (BD Biosciences) was used to quantify the required cells, and FlowJo software (TreeStar, Woodburn, OR, USA) was used to analyze the results.

Animal experiments

All animals involved in this study were treated humanely and received standard care. The animal experimental procedures were approved by the Institutional Review Board of Zhongshan Hospital of Fudan University (Shanghai, China). In this experiment, we housed male athymic nude mice (BALB/cASlac-nu) in a specific pathogen-free environment. We mixed treated A549 cells and TAMs to make a 1:1 cell mixture at a cell concentration of 5×10⁶ cells/ml. Take 0.05ml of the mixed suspension of cells and Matrigel, and implant them into the lung thoracic cavity of nude mice for in situ tumor formation experiments.

Animals were sacrificed when one of the following signs of disease was observed: tumor ulceration (greater than 0.5 cm); inability to move or eat; or serious injury. Changes in tumor size were detected using an optical imaging system for in vivo small animals (IVIS Spectrum, PerkinElmer, USA).

Statistical Analysis

The statistical tools, methods, and thresholds of each analysis are clearly described in the results or detailed in the legend or materials and methods.

Results

Single-cell transcriptomic profiling of LUAD

A total of 9 patients with non-metastatic LUAD underwent lobectomy with curative intent in the Department of Thoracic Surgery, Zhongshan Hospital of Fudan University. Among them, five received three cycles of preoperative neoadjuvant combination chemotherapy with cisplatin plus pemetrexed (defined as NCT group), while others only received surgery (defined as the Control group). Following resection, a malignant lung tumor sample was obtained from each patient, rapidly digested to a single-cell suspension, and analyzed using 10X scRNA-seq (Figure 1a). After quality control, a total of 83,622 cells that met the inclusion criteria were subjected to subsequent analyses, with 33,567 and 50,055 cells derived from the control and NCT groups, respectively (Figure 1a-c, S1a). Next, we classified cell types through dimensional reduction and unsupervised clustering using the Seurat package and relative maker genes.

Figure 1.

Using the SingleR package, the CellMarker dataset, and our previous studies (21, 22), we identified cell clusters that could be assigned to known cell lineages: epithelial cells (marked by SFTA2 and KRT8), T cells (marked by CD3D and TRBC2), B cells (marked by CD79A and CD19), endothelial cells (marked by EMCN and CXorf36), mast cells (marked by TPSB2 and TPSAB1), macrophages (marked by CD68 and APOE), monocytes (marked by FGL2 and LGALS2), fibroblasts (marked by LUM and DCN), neutrophils (marked by FCGR3B and CMTM2). Meanwhile, the consensus clustering of these cells also exhibited the consistency and homogeneity of the expression profile within each identified cell type (Figure 1b). For instance, clusters 1, 3, 6, 7, and 15, all designated as epithelial cells, were adjacent to each other in the consensus heatmap. This result confirms the robustness and reliability of our data pre-processing. Detailed distributions of these marker genes in each cluster are depicted in Figure S1.

By comparing the composition of different types of cells in each group, we noticed tumor microenvironment heterogeneity: the proportion of cells other than tumor cells, especially immune cells (mainly T and B), was significantly higher in the NCT (Figure 1e). Therefore, to identify subclusters within each of these nine major cell types, we observed a complex cellular ecosystem containing eight different epithelial subclusters and 43 non-epithelial clusters. Interestingly, the epithelial subclusters, mainly composed of cancer cells, were highly patient-specific, while the immune cell subclusters mostly consisted of cells derived from four or more patients (Figure 1f). This observation demonstrated the substantial variation and heterogeneity of tumor microenvironment among groups and individuals. Therefore, we further explored these alterations associated with the therapeutic regimen in greater detail for the primary cell types in subsequent analyses.

Metabolic reprogramming in lung adenocarcinoma driven by neoadjuvant chemotherapy

Metabolic reprogramming is a hallmark of malignant tumors. Recent studies have also shown that tumors’ metabolic characteristics and preferences change during cancer progression (23). In each type of cell derived from the Control and NCT groups, more significantly up-regulated metabolic pathways were enriched in cancer cells, nonmalignant epithelial cells, fibroblasts, and macrophages (Figure 2a). The enrichment of oxidative phosphorylation, glycolysis, pyruvate metabolism, and the tricarboxylic acid cycle indicates active glucose metabolism in these four cell types. By analyzing the activity of metabolic pathways in cells from different sources, we found that the activity scores of the metabolic pathways of tumor cells and macrophages were significantly higher than those of other types of cells. Notably, the metabolic pathway activity of macrophages and malignant cells increased after chemotherapy (Figure 2b).

Figure 2.

Changes in metabolism and gene expression of tumor cells after neoadjuvant chemotherapy

To accurately analyze the effect of chemotherapy on the cancer cells, we first re-clustered the epithelial cells, and 12 clusters were identified (Figure 3a). Copy number variations (CNVs) (Figure S2a) and marker genes were used to accurately separate malignant and nonmalignant epithelial cells in Control and NCT samples. They were finally defined as MalignantSA cells (Marker genes: FOXL2/MET/CD74), MalignantNCT cells (Marker genes: RAC1/MAF/CXCL1), and Nonmalignant cells (Marker genes: ABCA3/SFTPB/LPCAT) (Figure 3d). These marker genes were further confirmed by immunofluorescence experiments (Figure S2d). We found that the proportion of malignant cells was significantly reduced after chemotherapy (Figure 3b, c, Figure S2b). Although malignant cells were significantly reduced after chemotherapy, genetic aberrations by CNVs analysis revealed that MalignantNCT cells exhibited significantly higher malignant scores compared to MalignantSA cells (Figure S2a).

Figure 3.

We performed trajectory analysis to track the reprogramming of epithelial cells across the three groups. Nonmalignant cells evolved in two directions and developed into two clusters of cells (Figure 3e). In this evolutionary process, glycolysis-related genes (ENO1, LDHB, GAPDH), oxidative phosphorylation-related genes (NDUFA4), mitochondrial repair-related genes (TOMM7), glucose and lipid metabolism regulation genes (S100A16), ATPase activity-related genes (CCT6A), tumor immune regulation-related genes (CCL20, CXCL1, PAEP, PPP1R14B), hypoxia response regulation genes (CHCHD2), apoptosis regulation genes (MEG3, CEACAM5), mRNA alternative splicing-related genes (LSM5), and Ras-related protein (RAC1, RALA) gradually increased over time in the pseudotime analysis. These findings indicate that these genes play an essential role in the transformation of epithelial cells into tumor cells (Figure 3f, Figure S2e). Correspondingly, during the process of epithelial cells transforming into malignant tumor cells, the activity of the glycolysis pathway, oxidative phosphorylation pathway, angiogenesis pathway, DNA repair pathway, mTORC1 signaling pathway gradually increased over time. However, P53 pathway, apoptosis signaling pathway activity then steadily decreased (Figure S2f).

Similarly, we performed GSVA analysis on malignant and nonmalignant cells from the three groups. We found that the glycolysis pathway, oxidative phosphorylation pathway, MYC-targets, E2F-targets, DNA repair pathway, and mTORC1 signaling pathway were significantly enriched in MalignantNCT cells derived from the NCT group (Figure 3g). The metabolic reprogramming enables cancer cells to resist anti-cancer drugs, thereby developing chemoresistance (24). To find the hub genes that cause the malignant transformation of epithelial cells, through Single-Cell Regulatory Network Inference and Clustering (SCENIC) analysis, we found that E2F1, BRCA1, PURA, NKX2-1, NFIC, ETV7, STAT1, EGR1, and CEBPD were highly expressed in malignant cells (cluster 2, 6, 11) from the Control group. In contrast, the malignant cells from the NCT group (clusters 1, 7, 8, 9) have high expression of transcription factors (TFs) such as PATZ1, SIX5, BATF, IRF1, FOXA1, and CEBPG. After neoadjuvant chemotherapy, the increased expression of these TFs promoted the occurrence of lung adenocarcinoma complex phenotypic remodeling (Figure 3h).

Tumor cells have significant heterogeneity. We re-clustered the malignant cells and obtained 13 sub-clusters (Figure 3i). Cluster 1, 2, 5, 6, 8, and 10 were derived from the NCT group (Figure 3i). Through the analysis of the metabolism of these cell subclusters, we found that clusters 5 (marker genes: PCP4/NPW/VSIG1), 6 (marker genes: ERG1/HSPA6), and 10 (marker genes: C9orf172/SLC39A10) from the NCT group showed high levels of glycolysis, oxidative phosphorylation and pyruvate metabolism (Figure 3k). GSVA analysis also showed that the glycolysis and oxidative phosphorylation signaling pathway-related genes were significantly enriched in clusters 5, 6, and 10 (Figure 3l). Similarly, we re-clustered nonmalignant cells to obtain 16 sub-clusters, of which clusters 1, 3, 8, 11, 12, 13, 15 were from the NCT group, and the rest were from the Control group (Figure S2g, h). Clusters 4 and 7 from the Control group showed high levels of glycolysis and oxidative phosphorylation (Figure S2i, j), which contrasts with the glucose metabolism observed in malignant cells from the Control group.

Changes in stroma cells resulted from neoadjuvant chemotherapy

To investigate stromal cell dynamics in the tumor microenvironment (TME), we obtained 8944 presumed stromal cells, as shown in Figure 1c. We re-clustered them into five sub-populations, including COL14A1-positive fibroblasts, endothelial-1, endothelial-2, myofibroblasts, pericytes, and smooth muscle cells (SMC) (Figure 4a-c) (25–27). Detailed expression of the marker genes in each cell type is outlined in Figure 4c. Herein, we noticed a significant difference between the distribution of each of these five clusters in patients receiving varied types of treatment. The COL14A1-positive fibroblasts comprised the main fibroblast types in NCT groups, in which both endothelial 1 & 2 were mainly found. Pericyte and SMC were presented in all three groups. In contrast, myofibroblasts exclusively originated from the control group. According to previous research, myofibroblasts have been described as cancer-associated fibroblasts that participate in extensive tissue remodeling, angiogenesis, and tumor progression (25, 26). Therefore, this finding revealed that NCT and immunotherapy significantly altered the stromal cell composition in the tumor microenvironment.

Figure 4.

To explore the activity of known biological pathways in these stromal cells, we performed functional enrichment analysis. In particular, GSVA analysis exhibited that endothelial 1 & 2 shared several up-regulated pathways related to cell proliferation and fate regulation, including IL6-JAK-STAT3, TGFβ, and WNT-β catenin signaling. Besides, pathways associated with energy metabolisms such as glycolysis and hypoxia were up-regulated in myofibroblast, whereas pericyte was characterized by enriched oxidative phosphorylation and adipogenesis (Figure 4d). Meanwhile, when comparing the GSVA scores of these biological processes between patients from control or NCT groups, we noted that the stromal cells exhibited enhanced metabolic levels after NCT, as represented by up-regulated glycolysis, oxidative phosphorylation, and fatty acid metabolism pathways (Figure 4e).

Considering the essential role of fibroblasts and their complicated function in shaping the tumor microenvironment, we further re-clustered them into ten subgroups (Figure 4f-h). As shown in Figure 4i-j, the GSVA score of the metabolic pathways, including glycolysis and oxidative phosphorylation, and pyruvate metabolism and citrate cycle (TCA cycle), were up-regulated in clusters 5, 6, and 9. Intriguingly, the upregulation of these pathways was mainly observed in NCT groups (Figure S3a). The three clusters were represented by distinct gene expression profiles, such as overexpressed MYH11 in cluster 5, RGS5 in cluster 6, and TOP2A in cluster 9. Since the potential involvement of these genes in the manipulation of fibroblast metabolism has never been proposed yet, they might serve as new specific markers of the fibroblast subtype with such a high metabolic rate in the tumor microenvironment. Besides, the SCENIC analysis demonstrated that MEF2C, NFIA, and RAD21 might drive the formation of these clusters, respectively (Figure S3c). Further in vitro studies are required to elucidate these notable fibroblasts’ potential function and driver genes in LUAD’s development and response to NCT. Conclusively, cellular dynamics in stromal cells support a consistent phenotypic shift of fibroblasts towards an increased metabolic level after preoperative chemotherapy.

Chemotherapy drove tumor-associated macrophages to turn more into phenotypes that promote tumor progression

In the process of cancer formation, tumor-associated macrophages (Tumor-Associated Macrophages, TAM) have an essential influence on the inflammatory response in the tumor microenvironment (28). To study the effects of chemotherapy on TAMs, we first extracted all macrophages (10526 cells) and re-clustered them into ten cell clusters (Figure 5a). From Figure 1e, we can see that the proportion of macrophages after chemotherapy was reduced.

Figure 5.

The cell clusters derived from the Control group were 1/2/3/5/7 clusters, those from the NCT group were mainly 0/4/8 clusters, and the number of cells in the 6/9 clusters from the Control group and the NCT group was similar (Figure 5a). The proportion of cells in cluster 0 (marker genes: CXCL8/ CCL20/CHIT1), 4 (marker genes: CCL3/ CCL4/ SEPP1), 8 (marker genes: ARG2/ S100A2) decreased after chemotherapy, while the remaining cell clusters increased (Figure 5b, c). Through the GSVA analysis, we found that glycolysis, angiogenesis, PI3K-AKT-mTOR-signaling, IL6-JAK-STAT3-signaling, hypoxia, TGF-beta-signaling, and other signaling pathways were significantly enriched in cluster 0/1/8. Promoting inflammation-related signaling pathways such as TNF-signaling-via-NFKB, inflammatory-response, Notch-signaling, fatty-acid-metabolism, and oxidative-phosphorylation were increased dramatically in clusters 2/4/7/9 (Figure 5d).

Similarly, we found that glycolysis/gluconeogenesis, amino sugar and nucleotide sugar metabolism, alanine, aspartate, and glutamate metabolism were more active in the 0/1/8 cluster. In contrast, oxidative phosphorylation, citrate cycle, pyruvate metabolism, fatty acid elongation, fatty acid biosynthesis, etc., were more active in clusters 2/4/7/9 (Figure 5e). According to the GSVA analysis, the glycolysis/gluconeogenesis signaling pathway was significantly enriched in macrophages from the NCT group. In contrast, macrophages from the Control group showed a high activity in oxidative phosphorylation, fatty acid elongation, fatty acid degradation, fatty acid biosynthesis, and citrate cycle (TCA cycle) (Figure 5f). These results indicate that significant metabolic reprogramming occurred in tumor-associated macrophages after chemotherapy, and different TAMs cell clusters also showed huge metabolic differences. In general, our results revealed that chemotherapy could promote glycolysis of TAMs and inhibit fatty acid metabolism.

To explore the key genes that regulate the differences in the metabolism of each subcluster of macrophages, we performed a SCENIC analysis. We found that HES, PPARG, SPI1, CEBPB, and IRF7 were highly expressed in cluster 0/1, which may be the key genes that regulate the conversion of macrophages into M2-like TAMs, while clusters 2/4/7/9 highly expressed STAT1, STAT2, NFKB1, JUN, and FOS that regulate the conversion of macrophages to M1-like TAMs (Figure 5g).

According to the gene expression of macrophages, we divided these 10 clusters of cells into three subtypes of macrophages through cluster analysis (Figure 5h). We scored the expression levels of pro-inflammatory and anti-inflammatory genes in all macrophages. We displayed each color-coded macrophage subtype’s M1 and M2 scores (left) and pro-inflammatory and anti-inflammatory scores (right) through a scatter plot.

Similarly, we found that 0/1/8 cluster cells exhibited M2-like polarization and anti-inflammatory properties, while 2/4/7/9 exhibited M1-like polarization and pro-inflammatory properties (Figure 5i). Based on these analyses, we divided these 10 clusters of macrophage subtypes into three categories: M1-like polarized phenotype was defined as Anti-mac; M2-like polarized phenotype was defined as Pro-mac; those without obvious polarized phenotype were defined as Mix (Figure 5j). We found that the proportion of Pro-mac in the tumor microenvironment increased after chemotherapy, especially in the case of NCT-1 (Figure 5k). Interestingly, via trajectory analysis we found that two subtypes, Anti-mac and Mix, can be converted to Pro-mac. In this evolution process, the high expression of LYZ, FBP1, ALOX5AP, MARCO, S100A9, FN1, CXCL8, APOC1CTSL, and other genes may have played an essential role in promoting the conversion of Anti-mac to Pro-mac (Figure 5l,m). This suggests that we can change the phenotype of TAMs in the tumor microenvironment by altering the expression of these genes.

Chemo-driven Pro-mac and Anti-mac metabolic reprogramming exerted diametrically opposite effects on tumor cells

To further verify the remodeling effect of chemotherapy on the functional phenotype of TAMs in the tumor microenvironment, we first used the FindAllMarkers function in the Seurat package to find the marker genes of Pro-mac, Anti-mac, and Mix cells. Pro-mac was mainly characterized by high expression of CXCL8, ARG1, CREM, CD206, STAT6, CCL22, MMP7, and CCL3L3, while Anti-mac was mainly characterized by high expression of CD86, HLA-DR, PLAC8, CXCL10, COX2, IL15R, and SCGB3A1 (Figure 6a). Based on these marker genes, we sorted out Anti-mac cells (CD45+CD11b+CD86+) and Pro-mac cells (CD45+CD11b+ARG+) by flow cytometry (Figure 6b). To verify whether the cells we sorted were the cell population we wanted, we re-verified the positive rates of Pro-mac and Anti-mac cells by flow cytometry (Figure 6c). Our results showed that the proportion of Pro-mac cells in lung adenocarcinoma tissues after neoadjuvant chemotherapy increased significantly (Figure 6d). In fact, by performing immunofluorescence staining on lung adenocarcinoma tissue samples derived from surgery alone and neoadjuvant chemotherapy, we also found that the proportion of cells marked by the marker gene CD206 of M2-like TAMs increased significantly after chemotherapy (Figure 6e). Macrophages can promote tumor progression by secreting many cytokines. By analyzing the differentially expressed genes of Pro-mac and Anti-mac cells, we found IL10, PDCD1LG2, PDGF, VEGF, MMP9, CXCL9, CXCR4, IL22, KLF4, and TGF-β were highly expressed in Pro-mac cells that promote tumor growth, angiogenesis and suppress tumor immunity (Figure 6f). We obtained the Pro-mac and Anti-mac cells from 12 cases (6 cases of surgery alone, 6 cases of surgical samples after neoadjuvant chemotherapy) by flow cytometry. We named them Control Anti-mac, Control Pro-mac, NCT Anti-mac, NCT Pro-mac. After placing them in a cell culture flask for 24 hours, the content of some key cytokines in the supernatant of the culture medium was detected by enzyme-linked immunosorbent assay (ELISA). The levels of MMP9, EGF, and VEGF secreted by Pro-mac after neoadjuvant chemotherapy were significantly higher than those of Pro-mac from the surgery alone group. MMP9, EGF, VEGF, and IL10 secreted by Pro-mac were significantly higher than Anti-mac (Figure 6g). Similarly, when Control Anti-mac, Control Pro-mac, NCT Anti-mac, and NCT Pro-mac were inoculated subcutaneously with A549 cells at a ratio of 1:1 (Reinjection of macrophages two weeks later), we also found that NCT Pro-mac can significantly promote tumor growth. Interestingly, NCT Anti-mac in the tumor microenvironment after chemotherapy can significantly inhibit the growth of tumor cells, and this inhibitory ability was stronger than Control Anti-mac (Figure 6h, i).

Figure 6.

Our previous analysis found that Pro-mac glycolysis-related signaling pathways were significantly enriched, while in Anti-mac, oxidative phosphorylation and fatty acid metabolism signaling pathways were greatly enhanced (Figure 6j). In vitro experiments show that NCT Pro-mac’s ability to take up glucose and produce lactate was considerably more potent than other cells (Figure 6k). It was worth noting that the glycolysis level of NCT Anti-mac was markedly higher than that of Control Anti-mac (Figure 6l). When we placed the Pro-mac and Anti-mac in a 24-well plate and co-cultured with A549 cells in the Transwell chamber, we found that NCT Pro-mac can significantly enhance the invasion ability of A549 cells. At the same time, NCT Anti-mac showed a stronger ability to inhibit tumor invasion than Control Anti-mac (Figure 6m). However, when we used 2-DG (800uM, the concentration determined in pre-experiment) to inhibit the glycolysis of TAMs, the ability of Pro-mac to promote tumor progression was significantly weakened, and the power of NCT Anti-mac to suppress tumors was also considerably reduced (Figure 6m). By mixing these cells with macrophages for 3D culture, we found that the ability of NCT Anti-mac to inhibit tumor proliferation was significantly weakened when inhibiting its glycolytic activity. This showed that glycolysis could enhance the ability of Pro-mac to promote tumor progression and increase the capacity of Anti-mac to inhibit tumors (Figure 6n). Finally, through in vivo experiments, we inoculated a mixture of TAMs and A549 to nude mice and obtained the same experimental results as in Figure 6m/n (Figure 6o).

Chemotherapy treatment-induced remodeling of T and B cells

Considering the essential role of the tumor microenvironment, especially the immune infiltration level, in tumor development and response to therapy, we next investigated the characteristics of T and B cells. In our study, 22530 T cells were detected, which accounted for 26.9% of the total. We noticed that the re-clustered T cells could not be visibly distinguished among patients receiving different therapeutic regimens (Supplementary Figure 6a-b). According to the expression of a series of canonical markers of T cell subtypes, the T cells were divided into CD4+ T (marked by LTB, CD45RO, etc.), CD8+ T (marked by NKG7, GZMA, GZMB, CD8A, etc.), and Tregs (marked by FOXP3, CTLA4, etc.) (12, 21, 29) (Supplementary Figure 6c-d). The detailed expression profile of these marker genes is exhibited in Supplementary Figure 6e. Meanwhile, aside from these previously published T-cell markers, we also noted the specific upregulation of several genes in a particular cluster. At the same time, their expression specificity has not been elucidated yet.

As the major executor of tumor immunology, CD8+ T cells are thought to differentiate into cytotoxic T cells (CTLs) and specifically recognize endogenous antigenic peptides presented by the major histocompatibility complex I, thereby eliminating tumor cells (30). By comparing the composition of T cell subtypes in LUAD cells derived from different groups, we found that the proportion of CD8+T cells in the NCT group was significantly higher than those in patients receiving only surgical treatment (Supplementary Figure 6d). Therefore, we focused on CD8+ T cells for subsequent analyses and re-clustered them into five new subgroups, in which clusters 1-4 were mainly derived from the NCT group. In contrast, cluster 5 was predominantly enriched in the control group (Supplementary Figure 6f-k).

We next explored the expression profile of genes associated with T cell’s function in each CD8+ T sub-cluster. As depicted in Supplementary Figure 6i, clusters 1 and 2 were characterized by up-regulated naïve T cell markers, such as TCF7, LEF1, and CCR7, whereas genes associated with immune inhibition, like TIGIT, CTLA4, PDCD1, and HAVCR2, were enriched explicitly in cluster 3. Cytotoxic function-related genes, including GZMA GNLY, PRF1, GZMP, and GZMK, IFNG, IL2, were respectively overexpressed in clusters 4 and 5. Based on this evidence, we defined clusters 1 and 2 as naive T, three as regulatory/exhausted T, and 4 & 5 as cytotoxic T cells. Intriguingly, regarding both the sample origins and expression profiles of CD8+ T cells in clusters 4 and 5, we can reasonably hypothesize that NCT treatment potentially induces the reprogramming of CD8+ cytotoxic cells. To further verify this statement, we performed pseudotime-ordered trajectory analysis to monitor the dynamic view of CD8+ T cells’ reprogramming process via Monocle. As shown in Supplementary Figure 6l-p, three phases were detected in these clusters. Cluster 1, which exhibited the lowest cytotoxicity, was designated as the “root” state according to pseudotime.

In contrast, the immune inhibition-related genes like LAG3, TIGIT, and PDCD1, and cytotoxicity-related genes such as GZMB and IFNG were respectively activated in phases 2 and 3. This phenomenon is consistent with our T cell phenotype classification mentioned above. Then, our results showed differentiation paths from naive T to Treg/exhausted cells and cytotoxic cells. Considering the transcriptional changes associated with T cell reprogramming, naive T cells (phase 1) expressing high CCR6 and TCF7 differentiate into two distinct fates, clusters 4 and 5, in phase 3. Notably, the cells positioned at the cluster 4 branch were characterized by higher cytotoxicity than in cluster 5 (Supplementary Figure 6l, m, o). Regarding the sample origins of the two clusters, these findings demonstrated that NCT treatment ignites a relatively more robust immune cytotoxic response towards tumor cells, which could be partly explained by the excessive production of neoantigen caused by NCT-induced DNA damage.

SCENIC analyses suggested that distinct transcriptional mechanisms drove the differentiation of naive T cells to either cluster 4 or 5. As revealed in Supplementary Figure 6q, the cytotoxic cells derived from NCT-treated LUAD patients (cluster 4) were characterized by increased activation of FOSL2-extended, REL, YBX1, and NF-KB pathways. In contrast, those from the control group (cluster 5) had up-regulated JUN, FOSB, and ELF3 extended pathways. Together, our results revealed that preoperative chemotherapy prompts the naïve T cells to differentiate towards a more cytotoxic phenotype.

As for B cells, only 3902 (4.6%) cells were detected. 475 cells were derived from the control group, while 3427 were from the NCT group (Supplementary Figure 4a). Herein, we re-clustered the B cells into two sub-clusters. Based on canonical cell markers, class-switched memory B-cells (marked by CD19, CD37, and HLA-DRA) and plasma cells (marked by IGHA2, IGHG4, and CD38) were defined (Supplementary Figure 4a-c). The former compromised the majority of the total B cells (80.7%). Notably, the sample origins of the B cells demonstrated that a higher proportion of plasma cells characterized the control groups. In contrast, the class-switched memory B cells were significantly enriched in preoperatively treated patients.

Meanwhile, we performed GSVA analysis to explore several key biological pathways in the B cells derived from different groups. As depicted in Supplementary Figure 4d, B cells from the control group exhibited significant activating ways associated with metabolism and energy supply, including glycolysis and oxidative phosphorylation. However, the B cells derived from the NCT group exerted essential roles in most of the pathways, including glycolysis, fatty acid metabolism, apoptosis, and hypoxia. Overall, our observations demonstrated that NCT not only induced T cell reprogramming but also extensively impacted the composition and function of B cells in the tumor microenvironment.

Crosstalk among tumor and immune cells

The tumor microenvironment consists of numerous cell types, and the importance of crosstalk between cancer and immune cells has been implicated in various biological processes associated with tumor development (21, 29, 31). As depicted in Supplementary Figure 7a-b and Supplementary Figure 5a, the interactions between malignant cells and macrophages exhibited the strongest activity in both control and NCT groups, highlighting the important role of the macrophage in tumor immunology. Notably, we noted that the cell-to-cell communications among different cell types, especially between tumoral and immune cells such as cytotoxic CD8+ T, Treg, and memory B, were significantly strengthened in the NCT group. Specifically, we further investigated the ligand-receptor atlas within and between tumor cells and immune cells, which seemed to be quite reshaped by NCT (Supplementary Figure 7c-d, Supplementary Figure 5b). For example, MIF-CXCR4, whose activation usually promotes leukocyte recruitment (32), was increasingly activated in the NCT group between malignant and memory B, CD4+ T, and cytotoxic CD8+ T, whereas inhibited in macrophages. Meanwhile, MDK-NCL exhibited a similar activating phenotype with MIF-CXCR4, but its function in shaping the tumor microenvironment has never been reported. So, it might serve as a potential target of immune checkpoint inhibitor treatment in the future.

Given the above-mentioned NCT-induced immune activation, which was characterized by CD8+ T with higher cytotoxicity and an increased proportion of class-switched memory B cells, these findings further clarified that NCT could ignite a strong intrinsic immune response towards tumor cells. However, the inhibitory interaction pairs LGALS9-CD44 and LGALS9-HAVCR2 was abnormally activated in the NCT group between malignant and several T cells or macrophages (33). Its exact role in such conditions still requires further exploration.

In summary, our study revealed that the LUAD tissues that have experienced NCT had a distinct landscape of intracellular interactions, which might provide new ideas for future research focusing on implementing immunotherapy in the comprehensive anti-tumor therapeutic regimen.

Discussion

Although important advances in chemotherapy have reduced the mortality of cancer patients, the 5-year survival rate is still low, mainly due to the inherent or acquired mechanism of anti-tumor drug resistance (34). Chemoresistance results from complex reprogramming processes, such as drug export/import, drug detoxification, DNA damage repair, and cell apoptosis. Recently, the correlation between metabolic regulation and chemoresistance has received great attention. More efforts are devoted to targeting cell metabolism to overcome chemoresistance (35). The classic mechanism is to target the transport of anti-cancer drugs by increasing the activity of the efflux pump, such as the adenosine triphosphate (ATP) binding cassette (ABC) transporter. Cancer cells exhibit a special metabolic phenotype-aerobic glycolysis, quickly transporting and consuming glucose to produce ATP and promote drug efflux. PI3K/AKT pathway is activated by producing 3’-phosphorylated phosphoinositol, which is an important signaling pathway for lung cancer MDR (36). Glycolysis is beneficial to cancer cells by producing ATP faster, providing many intermediates for violent biosynthesis, maintaining redox balance, and creating a microenvironment with low immunity (24). The combination therapy of shikonin+2-DG could inhibit glycolytic phenotype, migration, and invasion by regulating the Akt/HIF1α/HK-2 signal axis (37).

Normal and healthy cells mainly produce energy through OXPHOS. However, due to rapid cell growth and frequent division, cancer cells face impressive metabolic challenges, which force them to adjust their energy metabolism to meet these needs (38) It is generally believed that cancer cells mainly obtain energy through glycolysis, which is named the Warburg effect. After chemotherapy, cancer cells change their metabolism from glycolysis to OXPHOS. This process is regulated by the SIRT1-PGC1α signaling pathway, thus increasing the resistance of cells to chemotherapy (38) Drug-resistant cancer cells can often be re-sensitized to anti-cancer treatments by targeting the metabolic pathways of import, catabolism, and synthesis of basic cell components (40). Recent studies have determined the cancer-promoting function of mitochondrial oxidative phosphorylation (OXPHOS) by regulating cell growth and redox homeostasis (41). Our study also found that after chemotherapy, the glycolysis and oxidative phosphorylation of tumor cells was enhanced. This metabolic reprogramming may enable cancer cells to have higher proliferation, invasion, and metastasis capabilities.

Tumor endothelial cells (ECs) have high glycolytic metabolism, shunting intermediates to nucleotide synthesis. Blocking of the glycolysis activator PFKFB3 in EC cells does not affect tumor growth. Still, it reduces cancer cell invasion, intravascular, and metastasis by normalizing tumor blood vessels, thereby improving blood vessel maturation and perfusion. PFKFB3 inhibition tightens the vascular barrier by reducing VE-cadherin endocytosis in endothelial cells and reduces glycolysis to make cells more quiescent and adherent (by up-regulating N-cadherin); it also reduces NF-κB signaling to reduce the expression of cancer cell adhesion molecules in ECs. PFKFB3 blockade therapy also improves chemotherapy for primary and metastatic tumors (42).

Due to rapid cell growth and frequent division in tumor cells, cancer cells face impressive metabolic challenges, which force them to adjust energy metabolism to meet these needs, namely metabolic reprogramming (43). However, studies have shown that metabolic plasticity in tumors is contributed by the glycolytic phenotype (as explained by Warburg) and that mitochondrial energy reprogramming has recently been identified as a feature of tumors (44). Chemotherapy can increase SIRT1/PGC1α-dependent oxidative phosphorylation (OXPHOS) in tumor cells, thereby promoting the survival of colorectal tumors during treatment. This phenomenon was also observed in chemotherapy-exposed liver metastases, which strongly suggests that chemotherapy causes long-term changes in tumor metabolism, which may interfere with drug efficacy (39). In addition, elevated glycolysis and OXPHOS promote epithelial-mesenchymal transition and cancer stem cell (CSC) phenotype in tumor cells (45). Therefore, recent research emphasizes the mixed glycolysis/OXPHOS phenotype rather than the phenotype that relies excessively on glycolysis to meet cellular energy requirements, thereby significantly promoting aggressiveness and treatment resistance (44). Chemotherapy has a significant effect on the metabolic reprogramming of tumor cells and profoundly affects stromal and immune cells’ metabolism in the tumor microenvironment.

Theoretically, chemical drugs can inhibit tumorigenesis by blocking the proliferation of tumor cells or depositing in tumor cell apoptosis, but this unintentionally causes “tissue damage.” The body will mistake this tumor-specific damage for normal tissue damage and then inevitably activate the tissue damage repair mechanism dominated by TAM (46). The result of this effect is that tumors will grow rapidly, and patients will develop resistance to anti-tumor chemotherapy.

Meanwhile, as suggested by Parra et al., neoadjuvant chemotherapy exerted PD-L1 upregulation in NSCLC patients. It increased the density of CD68+ macrophages, which were associated with better outcomes in both univariate and multivariate analyses (47). However, opposite results were also reported by Talebian et al. that NSCLC patients treated with radiotherapy, rather than a platinum-based standard-of-care chemotherapy, displayed a decrease in lymphoid cells and a relative increase in macrophages (48). Therefore, the role of TAMs in LUAD cells’ response to chemotherapy still requires further investigation.

Our research reveals the remodeling effect of chemotherapy on tumor microenvironment. However, this study still has many limitations. Firstly, only 9 samples were included in our study, and the number of samples is a defect of our study. Secondly, our study did not detect other causes of tumor heterogeneity, such as EGFR-mutant or ALK-translocated. We explored the possible impact of these factors and found that there was no significant difference in the expression of EGFR and other genes between the NCT group and the Control group (Supplementary Table 1). Thirdly, our data can only reflect the change in gene expression of various types of cells in the tumor microenvironment after chemotherapy. We can not draw a direct conclusion on whether chemotherapy will benefit patients or not. These need further study in the future.

Availability of data and material

The single-cell sequencing data used in this study and meta data (including) can be obtained from the figshare (10.6084/m9.figshare.24797265).

Acknowledgements

We thank OE Biotechnology Co., Ltd. (Shanghai, China) for helping us perform single-cell sequencing analysis of lung adenocarcinoma surgical samples. We thank International Science Editing Co. for providing language editing services for our papers.

Additional information

Author contributions

S.Y.Y, Y.W.H., G.S.B, C.Z., and G.J.W. designed the research; Y.W.H., L.C., J.Q.L., T.L., M.N.Z., G.S.B., M.L., J.Q.L., Z.Y.H., Y.S.Z., J.J.X, C.Z., Z.W.L., W.J., Q.W., and L.J.T. performed the research; Y.W.H., G.J.W., C.Z., L.X., and G.S.B. analyzed the data; Y.W.H., G.J.W., G.S.B., J.Q.L., L.X., and C.Z. wrote the paper. Final approval of the manuscript: All authors.

Declarations

A total of 9 patients were included in this study. The Ethics Committee of Zhongshan Hospital of Fudan University approved our research (B2021-230R), and the content of the study complies with the Helsinki Declaration.

A total of 60 animals were involved in this study. They were treated humanely and received standard care. The animal experimental procedures were approved by the Institutional Review Board of Zhongshan Hospital of Fudan University (Shanghai, China).

Competing interests

The authors have no conflicts of interest to declare.

Consent for publication

All of the authors have consented to publication.

Funding

This research was funded by the Sailing Program of Shanghai (No. 22YF107300), the National Natural Science Foundation of China (No. 82203645), the China Postdoctoral Science Foundation (No. 2023M730655), the Natural Science Foundation of Shanghai, China (No. 22ZR1411900), and the Special Foundation for Supporting Biomedical Technology of Shanghai, China (No. 22S11900300).

Supplementary figures

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Abbreviations

• LUAD: lung adenocarcinoma

• NCT: neoadjuvant chemotherapy

• TME: tumor microenvironment

• TSNE: t-distributed random neighborhood embedding

• CNV: copy number variation

• SCENIC: Single-Cell Regulatory Network Inference and Clustering

• SMC: smooth muscle cells

• TCA cycle: pyruvate metabolism and citrate cycle

• TAM: Tumor-Associated Macrophages

• CTLs: cytotoxic T cells

• OXPHOS: oxidative phosphorylation

• CSC: cancer stem cell

• GSVA: Gene Set Variation Analysis

• NMF: Non-negative Matrix Factorization