Skip to main content

Comparative profiling of single-cell transcriptome reveals heterogeneity of tumor microenvironment between solid and acinar lung adenocarcinoma

Abstract

Background

The diversity of histologic composition reflects the inter- and intra-tumor heterogeneity of lung adenocarcinomas (LUADs) macroscopically. Insights into the oncological characteristics and tumor microenvironment (TME) of different histologic subtypes of LUAD at the single-cell level can help identify potential therapeutic vulnerabilities and combinational approaches to improve the survival of LUAD patients.

Methods

Through comparative profiling of cell communities defined by scRNA-seq data, we characterized the TME of LUAD samples of distinct histologic subtypes, with relevant results further confirmed in multiple bulk transcriptomic, proteomic datasets and an independent immunohistochemical validation cohort.

Results

We find that the hypoxic and acidic situation is the worst in the TME of solid LUADs compared to other histologic subtypes. Besides, the tumor metabolic preferences vary across histologic subtypes and may correspondingly impinge on the metabolism and function of immune cells. Remarkably, tumor cells from solid LUADs upregulate energy and substance metabolic activities, particularly the folate-mediated one-carbon metabolism and the key gene MTHFD2, which could serve as a potential therapeutic target. Additionally, ubiquitination modifications may also be involved in the progression of histologic patterns. Immunologically, solid LUADs are characterized by a predominance of exhausted T cells and immunosuppressive myeloid cells, where the hypoxic, acidified and nutrient-deprived TME has a non-negligible impact. Discrepancies in stromal cell function, evidenced by varying degrees of stromal remodeling and fibrosis, may also contribute to the specific immune phenotype of solid LUADs.

Conclusions

Overall, our research proposes several potential entry points to improve the immunosuppressive TME of solid LUADs, thereby synergistically potentiating their immunotherapeutic efficacy, and may provide precise therapeutic strategies for LUAD patients of distinct histologic subtype constitution.

Background

Invasive lung adenocarcinomas (LUADs) account for almost 70–90% of all surgically resected lung cancers [1]. The morphologic manifestations of invasive LUADs have been well characterized microscopically and are mainly differentiated into lepidic, papillary, acinar, micropapillary, and solid growth patterns [1, 2]. The diversity of histologic composition macroscopically reflects the inter-tumor and intra-tumor heterogeneity of LUADs, with most LUADs manifesting as a successive tissue transition between two or more histologic patterns [2]. Solid and acinar are two histologic subtypes of LUAD with high frequency, the solid type was identified as a histologic pattern stronger in aggressiveness, higher in grade, and worse in prognosis than the acinar type [2]. Identification of patients who may benefit from additional treatment after curative surgery for early LUAD has been a focus in the field of adjuvant therapy. An early study suggested that the solid type benefited from adjuvant chemotherapy in terms of disease-free survival (DFS) and specific DFS, while the acinar type did not [3]. Moreover, the latest grading system, introduced with 20% or more of high-grade patterns (including the solid pattern) as the cut-off for histologically high-grade LUADs, consistently demonstrated that those patients with high-grade LUADs could benefit from adjuvant chemotherapy [4, 5].

To our knowledge, no study as yet has conclusively reported any association between the histologic subtype and the response to targeted or immune therapies, both in the adjuvant and neoadjuvant scenarios (comprehensive histological assessment relies on radical surgical excision). Although the application of targeted therapies depends on driver gene detection, it has been shown that the solid type has a lower frequency of targetable alterations when compared to the acinar type, as evidenced by rare EGFR mutations and a higher frequency of KRAS mutations, which implicates that tyrosine kinase inhibitors (TKIs) may not be satisfactorily effective in solid-type tumors [6, 7]. More importantly, several studies have proposed that the solid type may benefit from additional immunotherapy, especially the perioperative immune checkpoint blocker (ICB) treatment, which is supported by its higher TMB and PD-L1 expression [6,7,8]. Nevertheless, the efficacy of immunotherapy can be hampered by various factors, including hypoxic, acidic conditions and tumor metabolic reprogramming, which is precisely what makes the prospects for combinational treatment so attractive [9]. Whether the histologic constitution of LUAD serves to guide the determination of future combinational therapeutic regimens is largely unknown.

Since previous research addressing the molecular underpinnings of histologic subtypes was mostly based on bulk tumor sequencing [6, 7, 10, 11]. In-depth comprehension of the oncological characteristics as well as the tumor microenvironment of different histologic subtypes of LUAD at the single-cell level can help identify potential therapeutic vulnerabilities and combinational approaches to improve the survival of LUAD patients. Here we performed scRNA-seq analysis on solid-, and acinar-prevalent LUAD tumors, and adjacent normal lung samples. Through the comparative profiling of scRNA-seq data-defined cell communities in the tumor microenvironment (TME), we brought some new insights into the molding of tumor cells to their surrounding circumstances and the metabolic preferences of tumor cells from distinct histologic subtypes. Relevant results were further validated in multiple bulk transcriptomic and proteomic datasets as well as in an independent immunohistochemical validation cohort. Moreover, the constitutional and functional differences of immune and stromal cell subpopulations were also addressed. Our study promises to provide some attractive clues and supporting evidences for histologic subtype-oriented therapy in patients with LUADs.

Methods

Data resources

The scRNA-seq data of untreated, primary LUAD samples were collected from two of our previously published researches [12, 13]. The histologic constituents of each tumor were re-accessed and recorded semi-quantitatively by two experienced pathologists (L.H. and S.L.) as suggested by the World Health Organization [2]. The clinicopathological information of the enrolled samples was summarized in Additional file 1: Table S1. For the complementation of scRNA-seq data, three additional LUAD cohorts comprising mRNA or protein expression datasets were also collected. The RNA-seq data for The Cancer Genome Atlas (TCGA) cohort was downloaded from UCSC Xena (http://xena.ucsc.edu/) and the corresponding histopathologic information was obtained from the work presented by The Cancer Genome Atlas Research Network [14]. The RNA-seq data and histopathological information for the East Asian ancestry (EAS) cohort were derived from Chen et al.'s work [15] and obtained from OncoSG (https://src.gisapps.org/OncoSG/)(16, 17). The RNA-seq and proteomic data, along with histopathological information for the Clinical Proteomic Tumor Analysis Consortium (CPTAC) cohort were obtained from Gillette et al.’s research [18]. The histologic categories of these LUAD samples were determined by their predominant histologic pattern. In addition, for exploring the correlations between tumor cell metabolic pathway activities and the ratio of exhausted CD8 + T cells, an independent LUAD scRNA-seq dataset was used [19], and the processed data (including the normalized log2TPM matrix together with the cell annotation table) was downloaded from the Gene Expression Omnibus database (Accession Code: GSE131907).

We also retrospectively collected an additional cohort from our center which included 16 patients with untreated, surgically resected LUADs, and the formalin-fixed paraffin-embedded (FFPE) samples from each patient were obtained. The basic clinical and histologic information for this validation cohort was provided in Additional file 1: Table S1. The study was conducted following the principles of the Declaration of Helsinki, and the study protocol was approved by the ethics committee of Shanghai Pulmonary Hospital. Because of the retrospective nature of the study, patient consent for inclusion was waived.

scRNA-seq data processing and analysis

For all enrolled samples with scRNA-seq data, the specific procedures for single cell suspension preparation, library construction and sequencing were detailed in the original article [12, 13]. Raw sequence files in fastq format were obtained and processed with Cell Ranger (version 4.0.0) coupled with the GRCh38 human reference genome to generate a gene count matrix for each sample. Then Seurat R package took over the downstream analytic procedures [20, 21]. First of all, cells met any following criteria were removed: (1) Cells with extreme feature counts (< 500 or > 6000); (2) Cells with extreme RNA counts (> 40,000); (3) > 10% reads aligned to mitochondria; (4) < 10% reads aligned to ribosome. Subsequently, we performed data normalization, identification of highly variable genes and principal component analysis (PCA) with Seurat’s classic workflow. Then, harmony algorithm was used to correct the potential batches among sample resources. Uniform Manifold Approximation and Projection (UMAP) was facilitated for dimension reduction with the parameter “reduction” set as “harmony”. The Seurat “FindNeighbors” and “FindClusters” functions were applied to detect communities and find cell clusters. And then the “FindAllMarkers” function with default settings was used to find markers for each identified cluster.

Cell annotation and doublets removal

All cells were preliminarily labeled as epithelium, immune and stromal cells (roughly, PTPRC for marking immune cells, EPCAM for epitheliums, VWF and COL1A2 for stromal cells). Annotation of primary immune cell types was implemented according to the expression of canonical cell markers and by inspecting the top marker genes of each cluster. Notably, cell clusters expressing two or more major cell lineage markers were manually removed due to their potential doublet identity (such as the co-expression of LYZ for myeloid cells and CD3E for T cells in one single cluster).

CNV estimation

Large-scale chromosomal copy number variations within the cancer cells at single-cell level were identified by inferCNV (https://github.com/broadinstitute/inferCNV) with stromal cells (fibroblasts and endothelium) as normal reference. The Hidden Markov Model (HMM) model was used to infer the variation degrees by setting “HMM = TRUE, denoise = TRUE”.

Differential expression analysis and gene set variation analysis (GSVA)

The Seurat “FindMarker” function was used to identify differentially expressed genes between cells from solid and acinar samples under the default Wilcoxon rank-sum test. To estimate pathway activity at single-cell level, we applied GSVA with default parameters, as implemented in the GSVA R package (version1.32.0), as previously described [22]. The gene sets of hallmarks (h.all.v7.2.symbols.gmt) and KEGG pathways (c2.cp.kegg.v7.2.symbols.gmt) used in this study were acquired from the GSEA website (https://www.gsea-msigdb.org/gsea/index.jsp) [23, 24]. The limma R package (version 3.42.2) was used to calculate the differential activities of pathways between groups. A Benjamini-Hochberg-corrected P value of ≤ 0.05 was used to identify significantly altered pathways.

Scoring of gene signatures in scRNA-seq

For exploring the functional state of immune cells from tumors of different histologic types, we curated gene expression signatures including the cytotoxic, exhausted and hypoxic signatures for T/NK cells [25]; the M1 and M2 phenotype signatures for macrophages [26]; the antigen presentation and immunosuppressive signatures for DCs [27]; the fibrillar collagens and epithelial mesenchymal transition signatures for fibroblasts [28]. The detailed list of genes for each signature and their reference sources were provided in the corresponding supplemental tables. The normalized weighted mean expression of those signatures was calculated by Seurat’s “AddModuleScore” function, and signatures expression among subgroups was compared by two-sided Wilcoxon rank-sum test.

Cellular fraction calculation

To compare the variation in the percentage composition of immune cell subpopulations between different histologic subgroups, we calculated the fraction of cell numbers of the resulting cell subpopulations to the total number of clustered cells at the individual sample level. The significance of differences among histologic subgroups for the fractions was compared using two-sided unpaired Wilcoxon rank-sum test.

Analyses on the GEPIA2 web server

The prognostic relevance of molecular expression in the TCGA LUAD cohort was analyzed on the GEPIA2 (http://gepia2.cancer-pku.cn/) web server [29]. Specifically, the “Survival Analysis” module was selected. Patient groups were divided by the median of molecular expression, the Kaplan–Meier curves were plotted to visualize survival differences, the cox proportional hazard ratio (HR) was calculated and Log-rank test was used for hypothesis test. Additionally, correlation analysis between the expression of MTHFD2 and UBE2S in the TCGA LUAD cohort was performed on the “Correlation Analysis” module by Pearson correlation test.

Molecular expression in bulk transcriptome and proteomic datasets

Differential expression analyses of molecules of interest among histologic types were performed on three RNA-seq datasets (log2-transformed FPKM expression values for TCGA and EAS cohorts, and log2-transformed RPKM for CPTAC cohort), and the CPTAC proteomic dataset (Two-component-normalized log2-transformed protein expression values). Pairwise comparisons of the solid type and other histological types were performed by two-sided Wilcoxon rank-sum tests.

Gene signatures enrichment in bulk RNA-seq

The enrichment analyses of gene signatures in the TCGA bulk RNA-seq dataset were performed by gsva function in the GSVA R package (version 1.42.0) with parameters “method = ssgsea, kcdf = Gaussian, abs.ranking = TRUE”. The detailed list of genes for each signature and their reference sources are provided in Additional file 2: Table S2. Comparisons of signature enrichment scores across histologic subtypes were performed using the Kruskal–Wallis test. Violin plots were drawn by the ggpubr R package (version 0.4.0) for visualization.

Metabolic pathway analysis

For comprehensive dissecting the metabolic differences of epitheliums across histologic subgroups in our scRNA-seq dataset, we quantified the metabolic activities of every single epithelium by applying the scMetabolism R tool pipeline with parameters “method = VISION, metabolism.type = KEGG” [30, 31]. The median absolute deviation value was then used to measure the degree of inter-subgroup variation in metabolic pathway activity, only the top 50% variable metabolic pathways were selected for heatmap plotting.

For the independent LUAD scRNA-seq dataset, only samples from tumor lung (tLung) were included [19]. We also removed two samples with less than 50 tumor cells, and nine samples remained for further analysis. CD8 + T cells in the dataset were annotated as CD8 low T, Cytotoxic CD8 + T, Exhausted CD8 + T and Naive CD8 + T, as described in the original manuscript [19]. The exhausted CD8 + T cell ratio of each sample was calculated. We next quantified the metabolic activities of every single tumor cell in the same manner as described above. For each sample, the score on a specific metabolic pathway was calculated as the average score of tumor cells in that sample on that metabolic pathway. Finally, Pearson’s correlation tests were applied to examine the correlations between exhausted CD8 + T cell ratio and tumor metabolic pathway activities, with correction for multiple testing by Benjamini–Hochberg method, only metabolic pathways with P values less than 0.1 were selected for heatmap mapping.

Immunohistochemistry

Tissues were fixed in 4% paraformaldehyde, embedded in paraffin, cut into sections, and placed on adhesion microscope slides. Sections were subjected to immunohistochemical (IHC) staining according to standard procedures. We performed the IHC by using the MTHFD2 mouse anti-human antibody (Abcam, ab56772). The primary antibody was incubated at 4 °C overnight followed by using the BOND™ Polymer Refine Detection Kit (Leica, DS9800) according to the manufacturer’s instructions. Whole slide scanning was performed using panoramic MIDI under a 40 × objective lens. For each slide, the histologic patterns were firstly identified according to the cellular structure of the tumor, then three to five non-overlapping fields of view for each histologic region were randomly captured at 100 × magnification, and the staining intensity of MTHFD2 was finally semi-quantified using the Image J software (1.53q) by transforming it into mean optical density [32]. The statistical difference in staining intensity of MTHFD2 between solid and lepidic/acinar was determined by the Wilcoxon rank-sum test.

Statistics

The statistical analyses involved in this study were described in the corresponding method section. All statistical analyses and data presentations were performed by the R program (versions 3.6.3 and 4.0.2). All reported P values were two-tailed, and P < 0.05 was considered statistically significant.

Results

Analysis of scRNA-seq data from histologically annotated LUAD samples

The present study was a repurposing of scRNA-seq data from two of our previously published researches [12, 13]. All surgically excised samples came from patients with untreated, primary non-metastatic LUADs. The histologic constituents of each tumor sample were assessed and recorded semi-quantitatively [1, 2]. We attempted to single out tumor samples with high histologic purity for the purpose of dissecting subtype-specific oncological and immunological characteristics. Collectively, four solid-type, four acinar-type LUAD samples, and five adjacent normal lung samples were enrolled in this study. The representative hematoxylin–eosin (HE) stained images clearly visualized the microscopic structure of the acinar and solid patterns (Additional file 7: Fig. S1A–B). There was a "near-pure" tumor with the solid pattern covering more than 70% of the whole tumor in each solid LUAD sample [33]. With regards to the acinar type, the proportion of acinar pattern in each sample was greater than 50%, with the content of solid/micropapillary patterns limited to less than 10%, allowing to minimize the impacts of high-grade histologic components. The clinicopathological information for all enrolled samples was summarized in Additional file 1: Table S1.

The single-cell transcriptomic profiles generated by each sample were then combined for integrated analysis. Following strict quality control procedures, a sparse matrix with 97,875 cells and 25,233 genes was obtained (Methods). Before performing unsupervised graph-based clustering analysis, potential batch effects between samples were assessed and eliminated. Subsequently, all cells were labeled preliminary based on the expression of canonical cell markers (roughly, PTPRC for immune cells, EPCAM for epitheliums, VWF and COL1A2 for stromal cells; Additional file 7: Fig. S1C-D). Among these cells, 30,208 (30.86%) originated from solid samples, 25,250 (25.80%) originated from acinar samples, and 42,417 (43.34%) originated from adjacent lung tissues.

Tumor cells from solid LUADs create a more anoxic and acidic TME

We then committed to comparing the transcriptional characteristics of tumor cells derived from solid or acinar samples. By inferring large-scale copy number variations from transcriptome information, extensive chromosomal aberrations were observed in tumor-derived epitheliums relative to stromal cells (Additional file 7: Fig S1E). Comparing solid and acinar samples using gene set variation analyses (GSVA) [34] revealed that hallmarks associated with aggressiveness and metabolic activity, such as G2M checkpoint, angiogenesis, epithelial-mesenchymal transition (EMT), MYC targets V1 and PI3K/AKT/mTOR signaling, were up-regulated in tumor cells from solid samples (Fig. 1A, Additional file 2: Table S2), which was consistent with a more aggressive histopathological phenotype of solid LUADs. Notably, immune response-related hallmarks (such as TNFα signaling via NF-κB, IL2-STAT5 signaling, and IL6-JAK-STAT3 signaling) were also significantly enriched in solid samples. These findings emphasized the invasiveness of tumor cells from solid LUADs as well as their adept immune evasion capabilities.

Fig. 1
figure 1

Tumor cells from solid LUADs create a more anoxic and acidic tumor microenvironment. A. Differentially enriched hallmarks (top) and KEGG pathways (bottom) between tumor cells from solid and acinar samples revealed by GSVA. B. Violin plots showing enrichment scores of tumor proliferating rate and hypoxia signatures by histologic subtypes in the TCGA LUAD cohort. Global differences were measured by the Kruskal–Wallis test. C. Violin plots of upregulated genes in solid LUAD tumor cells. D. Boxplots showing mRNA expression of HIF1A, LDHA, UBE2S and UBE2C by histologic subtypes in the TCGA LUAD cohort. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Comparisons were performed using two-sided Wilcoxon rank-sum test (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, n s not significant). E–H. Violin plots showing enrichment scores of glycolysis (E), lactate transmembrane transporter activity (F), checkpoint molecules (G) and ubiquitin mediated proteolysis (H) signatures by histologic subtypes in the TCGA LUAD cohort. Global differences were measured by the Kruskal–Wallis test

It has been well established that hypoxia and acidification characterized the tumor microenvironment [35]. When comparing tumor cells from solid and acinar samples, hypoxia and glycolysis hallmarks were found to be more prominent in the former (Additional file 2: Table S2). Consistently, by applying single sample enrichment analysis (ssGSEA) in bulk RNA-seq data from the TCGA LUAD cohort, we found the enrichment scores of the tumor proliferating rate [36] and hypoxia [25] signatures were increased stepwise with histologic progression (Fig. 1B; Additional file 2: Table S2). Moreover, the expression levels of hypoxia-inducible factor-1 alpha (HIF1A) and lactate dehydrogenase A (LDHA), were observably upregulated in tumor cells from solid LUADs (Fig. 1C, Additional file 2: Table S2). As the key mediator of hypoxic response, HIF1A was intimately linked to multiple aspects of antitumor immunobiological processes [37]. And LDHA was identified as a transcription target of the oncogene MYC and was required for enhanced glycolysis and malignant potential of tumor cells [38]. LDHA was also an essential enzyme in lactate production, which was reported to create an acidic tumor microenvironment and mediate immunotolerance [35]. The upregulation of HIF1A and LDHA in solid LUADs was further corroborated in three publicly available bulk RNA-seq cohorts for which the histologic annotation data were available (Fig. 1D, Additional file 8: Fig. S2A, B; Methods). Moreover, enhanced glycolysis and lactate transmembrane transporter activity were also observed in solid LUADs in the TCGA cohort (Fig. 1E, F; Additional file 2: Table S2). These findings supported the hypothesis that tumor cells of solid LUADs constructed a more anoxic and acidic peritumor niche, which not only facilitated malignant progression but also impaired antitumor effector cell function. Indeed, immune suppression by checkpoint molecules was found to be most intense in solid LUADs (Fig. 1G; Additional file 2: Table S2).

The role of epigenetic mechanisms in the differentiation of LUAD morphological subtypes

Previous research established the critical role of epigenetic mechanisms in determining the histologic identity of LUAD cells [10]. Intriguingly, we noticed the activity of the ubiquitin-mediated proteolysis pathway, as well as the expression of multiple ubiquitin-conjugating enzymes including UBE2S, UBE2E3 and UBE2C, were upregulated in tumor cells from solid samples (Fig. 1C; Additional file 2: Table S2). A consistent trend in the activity of the ubiquitin-mediated proteolysis across histologic subtypes was also observed in the TCGA cohort (Fig. 1H; Additional file 2: Table S2). Indeed, ubiquitin-conjugating enzymes have been shown to be overexpressed in various cancer types and involved in the regulation of a variety of cancer-associated biological processes [39]. Overexpression of UBE2S, for instance, has been shown to reduce the stability and activity of the p53 protein by increasing its ubiquitination [40]. Furthermore, UBE2S was demonstrated to be negatively correlated with the von Hippel-Lindau tumor-suppressor (pVHL), which mediated the ubiquitin-dependent proteolysis of HIF1A, in multiple tumor cell lines [41, 42]. Actually, ubiquitination engaged in metabolic reprogramming of cancer cells by modifying metabolic signaling pathways, transcription factors and metabolic enzymes [43]. The solid LUAD-enriched PI3K/AKT/mTOR and c-Myc signaling, could be enhanced by ubiquitination modifications in tumor cells [43]. Importantly, the differential expression of UBE2S and UBE2C across histologic subtypes were further corroborated in bulk transcriptome from the TCGA and the East Asian ancestry (EAS) cohorts, as well as proteome from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) cohort (Fig. 1D, Additional file 8: Fig. S2C–D; Methods). Briefly, our results hinted that post-transcriptional modification represented by ubiquitination might play a crucial part in the acquisition of morphological phenotypes and the development of adaptive metabolic reprogramming in LUADs.

Metabolic preferences vary between histologic subtypes

We noticed a striking difference in the enrichment of multiple metabolic pathways between tumor cells derived from solid and acinar samples (Fig. 1A, Additional file 2: Table S2). Tumor cells in solid samples, for instance, had enhanced glycolysis and oxidative phosphorylation, while tumor cells from acinar samples had increased nitrogen metabolism. We next quantified a comprehensive collection of metabolic pathways [30] to dissect the metabolic landscape of epitheliums across groups (Fig. 2A). Remarkably, activities of some metabolic pathways, such as glycolysis/gluconeogenesis, purine and pyrimidine metabolism, and the citrate cycle (TCA cycle), showed progressive increases in lung epitheliums/tumor cells from normal to acinar and then solid samples. Nevertheless, some metabolic pathways exhibited histologic pattern preference. Specifically, tumor cells from solid samples were distinguished by increased activity in sulfur, glycerolipid, glutathione, selenocompound metabolism, and folate biosynthesis, whereas tumor cells from acinar samples were recognized by enhanced activity in D-glutamine and D-glutamate, nitrogen, vitamin B6 metabolism, and arginine biosynthesis. From a global perspective, tumor cells from solid LUADs upregulated pathways in both energy and substance metabolism (including nucleotide, amino acid, lipid metabolism and TCA cycle), which corresponds to their requirement for robust growth, proliferation and invasion. Among all these pathways, folate-mediated one-carbon metabolism supports a range of anabolic processes essential for cancer cell survival and growth [44]. Importantly, we discovered that the methylenetetrahydrofolate dehydrogenase 2 (MTHFD2), a key enzyme in the transformation of folate metabolites, was markedly upregulated in tumor cells from solid LUADs, along with the enhancement of folate biosynthesis and one-carbon pool by folate (Fig. 1C, Fig. 2A; Additional file 2: Table S2). Encouragingly, the expression differences of MTHFD2 across histologic subtypes were further confirmed in the TCGA and EAS bulk transcriptome cohorts, as well as the CPTAC proteomic cohort (Fig. 2B–D). More importantly, immunohistochemical staining further verified at the clinical sample scale that the density of MTHFD2 was increased in tumor regions of the solid pattern compared to the lepidic/acinar pattern (Fig. 2E, Additional file 9: Fig. S3A). To our knowledge, overexpression of MTHFD2 promotes tumorigenesis and malignant progression of various malignancies by virtue of its dual functions in metabolism and epigenetic modification [44,45,46]. The overexpression of MTHFD2 was associated with a worse prognosis in the TCGA LUAD cohort (Additional file 9: Fig. S3B). Strikingly, MTHFD2 raises the basal and IFN-γ-stimulated PD-L1 expression in both basal and IFN-γ-stimulating conditions, which could be the epitome of the linking of onco-metabolic genes and tumor immunoresistance [45]. Intriguingly, a strong positive correlation between the expression of MTHFD2 and UBE2S was observed in the TCGA cohort (Additional file 9: Fig. S3C), suggesting a possible link between ubiquitination modifications and metabolic reprogramming in LUADs.

Fig. 2
figure 2

Metabolic preferences vary between histologic subtypes. A. Heatmap showing relative metabolic pathways activities in the three histologic subgroups. B-D. Boxplots showing MTHFD2 and C15orf48 mRNA expression in the TCGA (B) and the EAS cohorts (C), and protein expression in the CPTAC cohort (D). Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Comparisons were performed using two-sided Wilcoxon rank-sum test (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, n s not significant). E. Representative immunohistochemical images of MTHFD2 in lepidic, acinar and solid tumor regions. Scale bar, 50 μm. F. Violin plots showing the reactive oxygen species pathway score across the three histologic subgroups. Comparisons were performed by two-sided Wilcoxon rank-sum test. G. Kaplan–Meier survival curves showing the prognostic difference between the low and high C15orf48 expression groups in the TCGA LUAD cohort. H. Correlations between tumor metabolic pathway scores and the ratio of exhausted CD8 + T cells. P values were determined by Pearson's correlation test (*P < 0.1, **P < 0.05, ***P < 0.01)

What’s more, we noticed the enhancement of glutathione metabolism and the enrichment of reactive oxygen species (ROS) pathway in tumor cells from solid samples (Fig. 2A, F; Additional file 2: Table S2). This suggested that tumor cells in the solid type might confront a more intense oxidative stress reaction. Consistently, the apoptosis hallmark was significantly enriched in tumor cells from solid LUADs (Additional file 2: Table S2). ROS are a complex class of molecules with both pro-tumor and anti-tumor effects [47]. Excessive ROS can induce cell death via various approaches, such as apoptosis, necrosis and autophagy, thus limiting the progression of cancer [48]. Interestingly, we discovered that the gene C15orf48 (also known as modulator of cytochrome C oxidase during Inflammation; MOCCI), which may be involved in the production of reactive oxygen species, was significantly upregulated in tumor cells from solid LUADs (Fig. 1C). C15orf48 is located in the cytochrome C oxidase subunit (Complex IV) in mitochondria and has been reported to reduce mitochondrial membrane potential and ROS production during inflammation, resulting in cellular protection and attenuated immunoreaction [49]. The expression differences of C15orf48 among LUAD subtypes were verified in the TCGA and EAS bulk transcriptome cohorts, as well as the CPTAC proteomic cohort (Fig. 2B–D). Survival analysis further demonstrated that high expression of C15orf48 was associated with a poor prognosis in LUAD patients in the TCGA cohort (Fig. 2G).

We wondered whether subtype metabolic preferences affect antitumor immune function. Using a publicly available scRNA-seq dataset of LUAD patients [19], we investigated the correlations between cancer cell metabolic pathway activity and the ratio of exhausted CD8 + T cells (Fig. 2H, Methods). The results showed unequivocally that the activation of metabolic pathways such as galactose metabolism, glycolysis/gluconeogenesis, and pentose phosphate pathway, was positively correlated with the ratio of exhausted CD8 + T cells. Most of these pathways were upregulated in tumor cells from solid samples except for vitamin B6 metabolism, which was relatively more active in acinar samples. Remarkably, through the Leloir pathway, galactose can be metabolized to produce G6P as an alternative fuel for glycolysis [50]. For the complementary of scRNA-seq, we quantified the enrichment scores of one-carbon pool by folate, pyrimidine metabolism and galactose metabolism in the TCGA LUAD cohort with bulk RNA-seq data (Additional file 9: Fig. S3D–F; Additional file 2: Table S2). The results showed that these metabolic pathways exhibited cascade enhancement with histological progression and reached their highest in the solid subtype, which is consistent with our previous findings.

In summary, the metabolic preferences of different histologic subtypes were consistent with their corresponding biological behavioral traits, especially the upregulation of folate-mediated one-carbon metabolism and the key gene MTHFD2 by tumor cells of solid LUADs to satisfy their cellular replication and transcriptional requirements. More importantly, we also gained some observations on the influence of tumor metabolic preferences on anti-tumor immune responses, which might provide new insights for understanding the heterogeneity of LUAD subtypes and exploring combinatorial therapeutic strategies.

Compromised anti-tumor effects in T cells from solid samples

We next investigated the landscape of the tumor immune microenvironment across distinct histologic subtypes. Immune cells (n = 81,136) were re-clustered and identified as T cells, NK cells, myeloid cells, mast cells, plasmacytoid dendritic cells (pDCs), B cells, and plasma cells based on the expression of canonical markers (Additional file 10: Fig. S4A–D; Additional file 3: Table S3). One of our primary interests was the variation in the composition and functional status of T/NK subsets associated with distinct histologic subtypes. Of the 33,946 T/NK cells obtained after removal of doublets (Methods), 4204 (12.4%) cells were from acinar samples and 18,974 (55.9%) from solid samples. Re-clustering all of these cells revealed ten distinct subsets: four CD4 + T cell subsets (CD4-C1 to C3; CD3D + CD4 +) including one subset of regulatory T cell (Treg; FOXP3 +), three CD8 + T cell subsets (CD8-C1 to C3; CD3D + CD8A +), one cycling T cell subset (Cycling T; MKI67 + TOP2A +), and two subsets of NK cells (NK-C1 and C2, CD3D − TYROBP +) (Fig. 3A–C, and Additional file 10: Fig. S4E–G).

Fig. 3
figure 3

Compromised anti-tumor effects in solid LUAD-derived T cells. A. Cluster-colored UMAP plot of T/NK cells. B. UMAP plots showing the expression of chosen canonical T/NK cell marker genes. C. Heatmap depicting the normalized expression of canonical T/NK cell marker genes among clusters. D. Boxplot showing cellular fractions of each T/NK cluster in normal (n = 5), acinar (n = 4) and solid (n = 4) samples. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Line segments of different lengths indicate which two groups were compared. Two-sided Wilcoxon rank-sum test. E–G. Violin and box plots showing the cytotoxic (E), exhausted (F) and hypoxic (G) signature scores for CD8 + T cells in the three subgroups. Comparisons were performed by two-sided Wilcoxon rank-sum test. H. Violin plots showing immune checkpoints expression of T cells from solid and acinar samples. I. Differentially enriched KEGG pathways between CD8 + T cells from solid and acinar samples revealed by GSVA

The proportion of the CD8-C1 subset increased incrementally from normal to acinar and solid samples (Fig. 3D). CD8-C1 highly expressed ZNF683 and ITGAE, represented a subpopulation of tissue-resident memory T cells (TRMs), which was reported as a considerable origin of tumor neoantigen specific T cells (Fig. 3B, C) [51, 52]. Besides, the elevated expression of exhaustion-related markers (e.g., LAG3, PDCD1, HAVCR2 and ENTPD1) in CD8-C1 indicated that a considerable fraction of these cells was induced to an exhausted phenotype in the TME (Fig. 3C, Additional file 10: Fig. S4G). Indeed, shared TCRs were frequently observed between ZNF683 + and exhausted T cells in lung tumors [53]. CD8-C2 infiltrated relatively less in tumors, and exhibited an effector memory signature (high GZMK, IFNG and CCL4 but lack of GZMB) and the previously defined “precursor exhausted” T cell signature (high EOMES, NR4A2, and low TCF7) (Fig. 3C, Additional file 10: Fig. S4G) [52, 54]. CD8-C3 was designated as cytotoxic T cells based on their expression of GZMA, GZMB, GNLY and PRF1, and their relative abundance was reduced in tumors. Notably, we observed increased HIF1A expression in the exhausted CD8-C1 subpopulation (Fig. 3C). We further quantified the cytotoxicity, exhaustion, and hypoxia signatures to compare the functional status of CD8 + T cells from different origins (Fig. 3E–G, Additional file 4: Table S4) [25]. The results showed that the exhausted score and hypoxic score increased progressively from normal to acinar and solid samples, while correspondingly, the cytotoxic score of CD8 + T cells from solid samples was lower than that of normal samples but still higher than that of acinar samples. This might suggest that, despite the considerable cytotoxic activity of solid LUAD-derived CD8 + T cells, the intensification of exhaustion induction led to a compromised oncocidal effect, thus contributing to immune escape and poor prognosis of solid tumors. Moreover, multiple immune checkpoint molecules including LAG3, TIGIT and ENTPD1, were significantly upregulated in solid LUAD-derived CD8 + T cells (Fig. 3H, Additional file 4: Table S4). Meanwhile, consistent with our previous finding of a more hypoxic and acidic TME in solid LUADs, the significantly elevated hypoxic score of solid LUAD-derived CD8 + T cells might reflect the involvement of hypoxic mechanism in T cell dysfunction.

GSVA was further conducted to dissect the differences in pathway activities between CD8 + T cells of solid and acinar sample origins (Fig. 3I, Additional file 4: Table S4). The results revealed that cytokine-cytokine receptor interaction was significantly activated in solid LUAD-derived CD8 + T cells. Correspondingly, various cytokines and chemokines, such as CCL4L2, CCL3L1, CCL3, CXCL13 and IFNG, were significantly upregulated in solid LUAD-derived CD8 + T cells (Additional file 10: Fig. S4H, Additional file 4: Table S4H). In particular, multiple metabolic pathways were upregulated in solid LUAD-derived CD8 + T cells (Fig. 3I, Additional file 4: Table S4), some of which were also found to be upregulated in tumor cells from solid samples, such as glycolysis/gluconeogenesis, cysteine and methionine, taurine, hypotaurine, and pyruvate metabolism. This suggested that the metabolic patterns favored by histologic subtypes might have a substantial impact on the metabolism and function of tumor-infiltrating CD8 + T cells. We were gaining the recognition that cellular metabolism constituted an important modulator of T cell replication and function [55]. This prompted us to speculate that solid LUAD-derived T cells were forced to compete for nutrients and oxygen with malignant cells, which could result in metabolic deficiencies and diminished antitumor effects [56]. In comparison, CD8 + T cells derived from acinar samples were significantly enriched in leukocyte transendothelial migration, focal adhesion and T cell receptor signaling pathways (Fig. 3I).

It is widely acknowledged that CD4 + T cells play a prominent role in cancer immunosurveillance and immunotherapy [57]. For the three conventional CD4 + T cell populations, CD4-C1 was identified as effector memory CD4 + T cells based on their expression of ANXA1, CD40LG and GZMA, which infiltrated relatively less in solid samples compared to acinar samples (Fig. 3B–D, Additional file 10: Fig. S4G). Of note, CD4-C2 was designated as a subtype of follicular helper T cells by expressing CXCL13 and TOX2, as well as exhaustion markers including TIGIT, PDCD1, CTLA4 and MAF. This exhausted population, as well, upregulated HIF1A expression and accounted for a relatively higher proportion of solid samples (Fig. 3C, D). CD4-C3 was marked by naïve T cell features such as CCR7, SELL and LEF and was relatively more abundant in solid samples. The infiltration percentage of Treg was elevated in tumors, but comparable in acinar and solid samples (Fig. 3D). The GSVA further revealed that pathways including glycolysis/gluconeogenesis, porphyrin, and chlorophyll metabolism, were significantly upregulated in solid LUAD-derived CD4 + T cells (Additional file 10: Fig. S4I). Dissimilarly, acinar LUAD-derived CD4 + T cells were enriched in natural killer cell mediated cytotoxicity, antigen processing and presentation, and allograft rejection.

To summarize, when compared to acinar LUADs, more T cells in solid LUADs were reprogramed to an exhausted phenotype, and importantly, the anti-tumor effects of T cells in solid LUADs were at least partially circumscribed by the highly nutrient- and oxygen-deprived TME. On the other hand, this also suggested that the therapeutic regimen combining immunotherapy and vascular normalization therapy to improve the TME might be more beneficial for solid LUADs.

Immunosuppressive phenotype dominates myeloid cells in solid LUADs

Myeloid cells were further subdivided into monocytes/macrophages and dendritic cells based on the expression of canonical markers (Additional file 11: Fig. S5A–B). Sub-clustering of 36,443 monocytes/macrophages revealed nine subsets, including four subsets of alveolar macrophages (AMs; Macro-C1 and C6-C8: MACRO, FABP4), three subsets of monocyte-derived macrophages (Mo-Macs; Macro-C2 and C3: MAFB, CSF1R; Macro-C4: MAF, SPP1), one subset of actively cycling macrophage (Macro-C5: MKI67, TOP2A) and one subset of monocyte (VCAN, FCN1) (Fig. 4A, B, Additional file 11: Fig. S5C, D). AMs were dominated by the Macro-C1 subpopulation, which was abundant in adjacent normal lungs and had the lowest proportion in solid samples (Fig. 4C). Macro-C6-C8 constituted only a minor fraction of AMs, with Macro-C8 being marked by metallothionein genes, such as MT2A and MT1E, which might be vital in heavy metals binding and processing (Additional file 11: Fig. S5D) [58]. As for the Mo-Macs, the proportions of Macro-C2-C4 were increased in tumors relative to normal samples (Fig. 4C). Macro-C4 was highlighted by high SPP1 expression and absent MHC II expression, which was regarded as a marker of tumor-associated macrophages (TAMs) and was involved in angiogenesis and metastasis promotion (Fig. 4B, C; Additional file 11: Fig. S5D) [59]. Macro-C2 exhibited a conspicuous gene signature of M2 polarization (including MAF, APOE, CCL18 and CCL13), while Macro-C3 was dominated by the M1 signature and highly expressed MHC II molecules (Fig. 4B, Additional file 11: Fig. S5D). Gene ontology annotation of marker genes further revealed that Macro-C2 was associated with neutrophil activation and myeloid leukocyte migration (Additional file 11: Fig. S5E). While Macro-C3 was enriched in antigen processing and presentation, response to IFN-γ and regulation of lymphocyte activation (Additional file 11: Fig. S5F). To compare the general functional characteristics of macrophages between acinar and solid samples, we quantified the phenotypic polarization scores of all macrophages [26]. Intriguingly, we discovered that solid LUAD-derived macrophages had both significantly higher M1 and M2 scores (Fig. 4D, E; Additional file 5: Table S5). Nevertheless, the “M2-M1” score was also significantly higher in solid LUAD-derived macrophages, implying that the immunosuppressive M2 phenotype was predominant in macrophages from solid samples (Fig. 4F). Indeed, we discovered that solid LUAD-derived monocytes/macrophages increased expression of suppressive immune checkpoints including HAVCR2, ENTPD1 and VSIR (Fig. 4G, Additional file 5: Table S5). Beyond that, comparison of gene expression profiling revealed that genes implicated in immunosuppression (FOLR2, TMEM176B, CLEC10A) and chemotaxis (CXCL9, CXCL10, CCL13, CCL18) were upregulated in solid LUAD-derived macrophages (Fig. 4H, Additional file 5: Table S5). FOLR2 encoded the macrophage-specific folate receptor β, which was responsible for folate transportation and was regarded as an attractive therapeutic target for TAMs [60]. The tumor-promoting cytokine CCL18 was also found to be intimately involved in the induction of the M2 phenotype in macrophages [61, 62]. Complementarily, we dissected the functional characteristics of macrophages in bulk data from the TCGA LUAD cohort and obtained consistent results. Those were, signals regulating macrophage and dendritic cell traffic into tumor, as well as enrichment scores for M1 phenotype and myeloid-derived immunosuppression signature increased incrementally with histologic progression (Additional file 11: Fig. S5G–I; Additional file 2: Table S2) [36].

Fig. 4
figure 4

M2 phenotype dominates solid LUAD-derived macrophages. A. Cluster-colored UMAP plot of monocytes and macrophages. B. Heatmap of normalized expression of monocyte/macrophage marker genes among clusters. C. Boxplot showing cellular fractions of each monocyte/macrophage cluster in normal (n = 5), acinar (n = 4) and solid (n = 4) samples. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Line segments of different lengths indicate which two groups were compared. Two-sided Wilcoxon rank-sum test. D-F. Violin and box plots of M1 (D), M2 (E), and M1-M2 (F) signature scores for macrophages in the three subgroups. G. Violin plots showing immune checkpoints expression of monocytes/macrophages from solid and acinar samples. H. Volcano plot showing the upregulated genes in solid LUAD-derived macrophages

Dendritic cells are critical in initiating and maintaining anti-tumor T cell immunity. All 2129 dendritic cells were re-clustered and annotated as conventional dendritic cells, type I (cDC1s; CLEC9A); conventional dendritic cells, type II (cDC2s; CD1C); mature regulatory dendritic cells (mregDCs; LAMP3), and plasma dendritic cell (pDCs; LILR4) (Fig. 5A, B). Both cDC1 and pDC were slightly overrepresented in solid samples when compared to acinar samples, and the opposite was the case in cDC2 (Fig. 5C). mregDCs are DC cells in a specific differentiated state associated with tumor antigens, and their immunomodulatory program may be exploited by cancer cells to facilitate immune escape [63]. Although the proportion of mregDCs was comparable between the two subtypes, mregDCs in solid samples were induced to differentiate to a more immunosuppressive phenotype, which corresponded with the overexpression of immunoregulatory genes such as PD-L1, CD200, CMTM6, IDO1, SOCS1, SOCS2, EBI3 and IL4I1 (Fig. 5D, Additional file 5: Table S5). We further quantified the antigen-presentation and the immunosuppressive signatures to ascertain the functional difference between all infiltrating dendritic cells in the three subgroups (Fig. 5E, Additional file 5: Table S5) [27]. Remarkably, solid LUAD-derived DCs have both the highest antigen-presentation score and immunosuppressive score. As previously reported, solid LUADs exhibited a higher rate of somatic mutations, which implies a greater abundance of potential tumor neoantigens [6]. This might partially explain why solid LUAD-derived DCs exhibited increased antigen presentation activity. The tumor microenvironment of solid LUADs, on the other hand, might counteract the anti-tumor effects of DCs by inducing an immunosuppressive phenotype, yet the underlying mechanism of which remained unknown.

Fig. 5
figure 5

The immunosuppressive phenotype of solid LUAD-derived DCs. A UMAP plot colored by DC subtypes. B Heatmap of normalized expression of DC marker genes among subtypes. C. Boxplot showing cellular fractions of each DC subtype in normal (n = 5), acinar (n = 4) and solid (n = 4) samples. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Line segments of different lengths indicate which two groups were compared. Two-sided Wilcoxon rank-sum test. D. Violin plots of immunosuppressive molecules expression in mregDC from solid and acinar samples. E. Violin and box plots of antigen presentation (left) and immunosuppressive (right) signature scores for DCs in the three subgroups. Comparisons were performed by two-sided Wilcoxon rank-sum test

Altogether, our results demonstrated that immunosuppressed myeloid phenotypes predominated in solid LUADs, which also suggested that myeloid-directed therapeutic interventions were somehow attractive for counteracting the immunosuppressive TME and eliciting effective antitumor responses in solid LUADs.

Stromal status disparities in different histologic subtypes

Previous researches support that cancer-associated fibroblasts (CAFs) are implicated in cancer initiation, progression, extracellular matrix (ECM) remodeling, and treatment resistance [64, 65]. Here we addressed the discrepancies between solid LUAD- and acinar LUAD-derived fibroblasts in terms of phenotype constitution and functional status. A total of 1507 fibroblasts were re-clustered and divided into seven subsets including five subsets of fibroblasts (Fibro-C1-C5; COL1A1, COL1A2), one subset of pericyte (RGS5) and one subset of myofibroblast (ACTA2, MYH11) (Fig. 6A–C, Additional file 12: Fig. S6A–C). Acinar and solid samples contributed 585(38.87%) and 329(21.8%) fibroblasts, respectively (Fig. S6A). It was of particular interest that the Fibro-C2 subset, whose relative proportion decreased gradually from normal to acinar and then solid samples, was marked by low expression of fibroblast activation protein (FAP) and specifically overexpressed the transcription factor 21 (TCF21) (Fig. 6B, C, Additional file 12: Fig. S6B, C). Previous study demonstrated that TCF21 was a master regulator of CAF status, and its overexpression diminished the ability of CAF to contract collagen gels, promote tumor growth, invasion and chemotherapy resistance [66]. The Fibro-C1 subset, which was highlighted by classic CAF markers (e.g., FAP, PDPN, MMP11) as well as EMT signature genes (e.g., CTHRC1, POSTN, IGFBP3), was enriched in tumors relative to normal samples (Fig. 6B, C, Additional file 12: Fig. S6C). Simultaneously, Fibro-C1 up-regulated the expression of various collagen fibers (e.g.; COL1A1, COL1A2, COL3A1, COL10A1, COL11A1). This corresponded to the result of gene ontology annotation of the specific markers of Fibro-C1, which revealed that these tumor-enriched fibroblasts were possibly actively participated in the processes of collagen fibril organization and ECM remodeling (Additional file 12: Fig. S6D). Remodeling of the ECM by CAFs might result in a mechanical barrier preventing the migration and infiltration of immune effector cells into the tumor parenchyma [67]. This notion led us to further compare the functional differences of fibroblasts between the two subtypes. Analysis of differentially expressed genes disclosed that CTHRC1 and IGFBP5, both of which played important roles in collagen production and fibrotic progression, were significantly upregulated in solid LUAD-derived fibroblasts (Fig. 6D). In addition, the immunomodulatory gene TDO2 was also found to be upregulated in solid LUAD-derived fibroblasts. We further quantified the fibrillar collagens transcriptional score of all fibroblasts based on the expression of main collagens (i.e., collagens I/ III/V, and fibronectin) that were assembled into large mechanically resilient fibers (Additional file 6: Table S6) [28]. Comparisons between different histologic subgroups demonstrated that the fibrillar collagens transcriptional score was significantly higher in tumor tissues than in adjacent normal lung tissues, and more importantly, higher in solid LUADs than in acinar LUADs, suggesting a more extensive interstitial fibrosis in solid LUADs (Fig. 6E). Moreover, solid LUAD-derived fibroblasts had a significantly higher EMT score (Fig. 6F), which could contribute to the transition of malignant cells from an epithelial phenotype to a more migratory mesenchymal phenotype, thereby promoting distant tumor metastasis and treatment failure. By applying ssGSEA to bulk RNA-seq data from the TCGA LUAD cohort, we subsequently quantified and compared the enrichment scores of fibrillar collagen, matrix remodeling and EMT signatures [36] across histologic subtypes, and the results were consistent with our previous findings (Additional file 12: Fig. S6E–G; Additional file 2: Table S2).

Fig. 6
figure 6

Stromal status disparities in different histologic subtypes. A. Cluster-colored UMAP plot of all fibroblasts. B. Boxplot showing cellular fractions of each fibroblast cluster in normal (n = 5), acinar (n = 4) and solid (n = 4) samples. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Line segments of different lengths indicate which two groups were compared. Two-sided Wilcoxon rank-sum test. C. Heatmap of normalized expression of fibroblast marker genes among clusters. D. Volcano plot showing the upregulated genes in solid LUAD-derived fibroblasts. E–F. Violin and box plots of fibrillar collagens (E) and EMT (F) signature scores for fibroblasts in the three subgroups. Comparisons were conducted using two-sided Wilcoxon rank-sum test. G. Cluster-colored UMAP plot of all endothelial cells. H. Boxplot showing cellular fractions of each endothelium cluster in normal (n = 5), acinar (n = 4) and solid (n = 4) samples. Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5 × the interquartile range. Line segments of different lengths indicate which two groups were compared. Two-sided Wilcoxon rank-sum test. I. Heatmap of normalized expression of endothelial marker genes among clusters. J. Differentially enriched KEGG pathways between tumor endothelial cells from solid and acinar samples revealed by GSVA. K. Violin plots showing the expression of selected MHC molecules in endothelial cells

Endothelial cells in the TME not only serve as a physical barrier to tumor cell metastasis but also act as semi-professional antigen-presenting cells involved in the interactions with immune cells (e. g., recruitment and activation of T cells) [68]. Consistent with previous findings, endothelial cells were enriched in normal lungs on account of their hypervascular nature [69]. A total of 1422 endothelial cells, of which 294 (20.7%) and 163 (11.5%) cells were separately from acinar and solid samples, were further re-clustered. The resulting sub-clusters obtained were assigned to known endothelial types in line with the expression of well-established marker genes (Fig. 6G–I and Additional file 13: Fig. S7A–C). When comparing the percentage composition of endothelial categories in different histologic subgroups, we found that tumor endothelial cells (tumor ECs) account for the majority of endothelial population in tumor samples, especially in solid samples. Inversely, lymphatic ECs and capillary ECs were relatively absent in solid samples (Fig. 6H). Notably, tumor ECs were featured by upregulation of angiogenesis-related genes including collagens (e.g., COL4A1, COL4A2), PXDN, SPARC and HSPG2, as well as proliferation markers (e.g., MKI67, TOP2A) (Fig. 6I). These hyperproliferative endothelial cells did not hold the structure and function of normal blood vessels, but might instead exacerbate the hypoxic situation in the microenvironment of solid LUADs [70]. To further elucidate the functional status of tumor endothelial cells, GSVA was conducted and pathway activities were compared between solid and acinar samples. We uncovered that graft versus host disease, allograft rejection, antigen processing and presentation pathways were significantly enriched in acinar LUAD-derived endothelial cells (Fig. 6J). It was documented that increased antigen presentation by endothelial cells assisted in T cell priming and tumor-specific immune activation [71]. Differently, solid LUAD-derived endothelial cells upregulated multiple metabolic pathways such as glycine, serine and threonine metabolism, linoleic acid metabolism. This consisted with previous findings that tumor endothelial cells were transcriptionally reprogrammed to regulate their metabolic functions, which might be associated with the downregulation of their antigen presentation and homing immune cell recruitment functions, and ultimately contributed to tumor escape from immune destruction [71, 72]. In fact, we did find that multiple MHC class I and MHC class II molecules were under-expressed in solid LUAD-derived endothelial cells (Fig. 6K).

Taken together, these encouraging findings enlighten us that therapeutic regimen targeting stromal components may enhance antitumor benefit by facilitating immune-tumor cells interactions in solid LUADs and hold promise to synergize with immunotherapeutic regimens.

Discussion

Researches dedicated to dissecting the heterogenous TME of distinct histologic patterns of LUADs are well underway. Here we performed a comparative analysis of the TME characteristics between solid and acinar LUAD samples primarily from a single-cell perspective, with the relevant results complemented by bulk transcriptomic and proteomic datasets and validated by immunohistochemistry. Our results suggested that the degree of acidity and hypoxia and the tumor metabolic preferences varied between histologic subtypes and might correspondingly impinge on the metabolism and function of immune components. In addition, ubiquitination modifications might also be involved in the progression of histologic patterns. Indeed, a prevalent state of restrained immune effector function was found in the solid histologic subtype, and discrepancies in stromal cell function could also contribute to the specific immune phenotype of solid LUADs.

Metabolic remodeling fuels tumorigenesis and evolution. However, horizontally, metabolic heterogeneity exits amongst tumors as well as inside an individual tumor, and vertically, adaptive metabolic remodeling may arise along with malignant progression [73]. A concomitant concern is whether the progression of histologic patterns is accompanied by the transformation in metabolic profiles of LUADs. Here we did find a gradient of metabolic alterations and relatively specific metabolic preferences between histologic subtypes, these metabolic properties coincided with the malignant potential of the histologic subtypes and might have a direct or indirect impact on intra-tumoral immune function. Energetically, tumor cells from solid LUADs upregulated glycolytic activity, confronting immune cells, which also relied on glycolysis for effector functions, with a scarcer energy source [56, 74]. On substance metabolism, augmented serine metabolism and pentose phosphate metabolic pathways were observed to supply biomolecules for the exuberant cell replication of solid LUADs. Indeed, we observed a positive correlation between the metabolic activities of tumor glycolytic and pentose phosphate pathway and the proportion of exhausted CD8 + T cells through an independent scRNA-seq dataset of LUAD. Interestingly, unlike the utilization of oxidative phosphorylation by M2 macrophages, the predominant metabolic pattern of M1 macrophages is precisely the glycolytic and pentose phosphate pathway [75]. Furthermore, both hypoxia and higher lactate levels exhibited a facilitative effect on macrophage polarization toward the M2 phenotype. Furthermore, both the hypoxia and lactate levels facilitated the polarization of macrophages toward the M2 phenotype [76]. Hence, oxygen and nutrient deprivation, as well as acidification could be crucial contributors to the predominance of T cell exhaustion and M2 macrophage polarization in solid LUADs. Alongside tumor cells, immunosuppressive cell metabolism was also of concern, an example being the upregulation of IDO1 by solid LUAD-derived mregDCs, which attenuated effector T cell responses by depleting tryptophan and producing the immunosuppressive metabolite kynurenine [76].

Immune checkpoint blocker (ICB) therapies are currently held in high esteem, yet their overall response rate for monotherapy in lung cancer is barely satisfactory [77]. Our study shows that solid LUADs exhibit higher exhausted transcriptional score in T cells relative to acinar LUADs, which raises the question of whether ICB can improve the dysfunctional state of T cells in solid LUADs and enhance their immunosurveillance and immunocidal effects on tumors. It is therefore of great interest to specify the association of histologic composition with the response and the survival benefit of perioperative ICB treatment for LUADs. There is promise in targeting tumor metabolic pathways and moderating metabolic competition in the TME to improve the efficacy of immunotherapy [9, 76]. Study shows that tumors with low glycolysis rates are more sensitive to CTLA-4 blockade, while those tumors with higher glycolytic activities, for example, solid LUADs, may benefit more from the combinational treatment with glycolytic inhibitor and CTLA-4 blockade compared to other histologic subtypes of LUADs [78]. In addition, improved local supply of L-arginine to TME synergistically enhances the efficacy of PD-L1 blockers via a T cell-dependent manner [79]. Remarkably, with respect to other histologic subtypes, we find that folate-mediated one-carbon metabolism and its key gene, MTHFD2, are upregulated in tumor cells from the solid subtype. Intriguingly, the main contributor to the one-carbon unit, serine, is also observed to be metabolically enhanced in tumor cells from solid LUADs [80]. The incrementally increased expression of MTHFD2 with histologic progression is further substantiated at both transcriptional and translation levels by diverse experimental technical scales. MTHFD2, a folate cycling metabolizing enzyme, is considered an attractive metabolic target for tumors by virtue of its multiple roles in metabolic reprogramming and immune regulation [45, 81]. Actually, in the scenario of acute myeloid leukemia, specific MTHFD2 inhibitors potently and selectively repress cancer cell replication while protecting non-neoplastic lymphocytes [82]. To date, specific MTHFD2 inhibitors have shown impressive therapeutic efficacy in preclinical models of diverse cancer types, particularly exhibiting promising antitumor activities in lung adenocarcinoma cell lines, either as a single agent or in synergy with pemetrexed [44, 82,83,84].Considering the disparities in tumor metabolism and the immune microenvironment in distinct histologic subtypes of LUADs, metabolic targets represented by MTHFD2 hold promise for developing their application in histologic subtype-directed combinatorial immunotherapy. Notwithstanding, more efforts on exploring and validating the efficacy and clinical application protocols of MTHFD2 inhibitors in lung cancer are warranted. We hold out the prospect of adopting MTHFD2 inhibitor as monotherapy or in combination with immunotherapy for the treatment of tumors with elevated MTHFD2 expression, especially for solid LUADs.

Tavernari et al. demonstrated on the spatial scale that regions of solid LUADs exhibited a geographic signature of immune exclusion, i.e., immune infiltration was significantly declined from the immune-enriched tumor margins toward the tumor core, whereas suppressive immune markers such as FOXP3, TIM3 and CTLA4 were distinctively elevated [10]. Notably, the reasons for the spatial distribution and functional differences of immune cells in the solid histologic region remain elusive. Here we propose the following two potential explanations. Firstly, the potential contribution of differences in the spatial distribution of oxygen and nutrients in the tumor regions of solid LUADs; and secondly, the obstruction by ECM components to the migration and movement of immune cells. In the case of the former, we introduce here the tumor model proposed by Lloyd et al. whereby tumor cores tend to maximize their population density and exhibit static, less proliferative phenotypes, while tumor margins are characterized by aggressive proliferative phenotypes [85]. Intriguingly, this model fits highly with the spatial characteristics of the solid pattern of LUAD identified by Tavernari et al. [10]. The harsh metabolic microenvironment created by vicious competition for limited resources in the tumor core may be detrimental to the survival and functional execution of immune cells. Indeed, Lambrechts et al. also suggested that the degree of hypoxia increases progressively from the tumor margin toward the core, whereas most immune cells are inclined to accumulate at the normoxic tumor margin [69]. In the case of the latter, CAFs and their remodeling of the ECM are key factors in structuring the immune infiltration barrier [86]. Based on comparative analysis of transcriptional profiles of the identified fibroblasts, we find that the fibrillar collagen transcriptional level is significantly higher in solid LUAD-derived fibroblasts [28]. And bulk transcriptome-based analysis further confirmed the elevated fibrillar collagen transcription and extracellular matrix remodeling activities in solid LUADs. In addition, it was noteworthy that solid LUADs are often accompanied by substantial intracellular and extracellular mucus production and secretion [1]. This implies that therapeutic regimens targeting CAFs or local ECM potentially promote immune infiltration into the tumor core of solid LUADs, thereby increasing the inter-contact between immune and tumor compartments.

Conclusions

Collectively, we herein proposed some potential entry points to disrupt the immune exclusion and immunosuppressive phenotype and to potentiate immunotherapeutic efficacy for solid LUADs, yet the realization of these notions requires further investigation and validation at different experimental techniques scales, such as microdissection and spatial omics techniques, as well as tumor models. Furthermore, considering the prospect of possible future applications of histologic subtype-directed LUAD treatment, the development of methods to determine the histologic composition or the presence of certain key components in the tumor prior to treatment is crucial.

Availability of data and materials

All data used for this study are publicly published and available with detailed access links described in the Data resources section. No new algorithms were developed for this manuscript. All code generated for analysis is available from the authors upon request.

Abbreviations

LUAD:

Lung adenocarcinoma

TME:

Tumor microenvironment

TMB:

Tumor mutation burden

ICB:

Immune checkpoint blocker

EMT:

Epithelial-mesenchymal transition

GSVA:

Gene set variation analysis

ssGSEA:

Single sample gene set enrichment analysis

HIF1A:

Hypoxia-inducible factor-1 alpha

LDHA:

Lactate dehydrogenase A

TCGA:

The cancer genome atlas

EAS:

The east Asian ancestry

CPTAC:

The clinical proteomic tumor analysis consortium

MTHFD2:

Methylenetetrahydrofolate dehydrogenase 2

ROS:

Reactive oxygen species

pDCs:

Plasmacytoid dendritic cells

DCs:

Dendritic cells

cDCs:

Conventional dendritic cells

CAFs:

Cancer-associated fibroblasts

ECM:

Extracellular matrix

MHC:

Major histocompatibility complex

FFPE:

Formalin-fixed and paraffin-embedded

IHC:

Immunohistochemistry

References

  1. Travis WD, Brambilla E, Noguchi M, Nicholson AG, Geisinger KR, Yatabe Y, et al. International association for the study of lung cancer/american thoracic society/european respiratory society international multidisciplinary classification of lung adenocarcinoma. J Thorac Oncol O Publ Intern Assoc Study Lung Cancer. 2011;6(2):244–85.

    Google Scholar 

  2. Travis WD, Brambilla E, Nicholson AG, Yatabe Y, Austin JHM, Beasley MB, et al. The 2015 world health organization classification of lung tumors: impact of genetic, clinical and radiologic advances since the 2004 classification. J Thorac Oncol Off Publ Intern Assoc Study of Lung Cancer. 2015;10(9):1243–60.

    Google Scholar 

  3. Tsao MS, Marguet S, Le Teuff G, Lantuejoul S, Shepherd FA, Seymour L, et al. Subtype classification of lung adenocarcinoma predicts benefit from adjuvant chemotherapy in patients undergoing complete resection. J Clin Oncol Off J Am Soc Clin Oncol. 2015;33(30):3439–46.

    Article  CAS  Google Scholar 

  4. Moreira AL, Ocampo PSS, Xia Y, Zhong H, Russell PA, Minami Y, et al. A grading system for invasive pulmonary adenocarcinoma: a proposal from the international association for the study of lung cancer pathology committee. J Thorac Oncol Off Publ Intern Assoc Study Lung Cancer. 2020;15(10):1599–610.

    Google Scholar 

  5. Hou L, Wang T, Chen D, She Y, Deng J, Yang M, et al. Prognostic and predictive value of the newly proposed grading system of invasive pulmonary adenocarcinoma in Chinese patients: a retrospective multicohort study. Mod Pathol. 2022;35(6):749–56.

    Article  PubMed  Google Scholar 

  6. Caso R, Sanchez-Vega F, Tan KS, Mastrogiacomo B, Zhou J, Jones GD, et al. The underlying tumor genomics of predominant histologic subtypes in lung adenocarcinoma. J Thorac Oncol Off Publ Intern Assoc Study Lung Cancer. 2020;15(12):1844–56.

    CAS  Google Scholar 

  7. Dong ZY, Zhang C, Li YF, Su J, Xie Z, Liu SY, et al. Genetic and immune profiles of solid predominant lung adenocarcinoma reveal potential immunotherapeutic strategies. J Thorac oncol Off Publ Intern Assoc Study Lung Cancer. 2018;13(1):85–96.

    CAS  Google Scholar 

  8. Ci B, Yang DM, Cai L, Yang L, Girard L, Fujimoto J, et al. Molecular differences across invasive lung adenocarcinoma morphological subgroups. Transl lung cancer res. 2020;9(4):1029–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Li X, Wenes M, Romero P, Huang SC, Fendt SM, Ho PC. Navigating metabolic pathways to enhance antitumour immunity and immunotherapy. Nat Rev Clin Oncol. 2019;16(7):425–41.

    Article  CAS  PubMed  Google Scholar 

  10. Tavernari D, Battistello E, Dheilly E, Petruzzella AS, Mina M, Sordet-Dessimoz J, et al. Nongenetic evolution drives lung adenocarcinoma spatial heterogeneity and progression. Cancer Discov. 2021;11(6):1490–507.

    Article  CAS  PubMed  Google Scholar 

  11. Nguyen TT, Lee HS, Burt BM, Wu J, Zhang J, Amos CI, et al. A lepidic gene signature predicts patient prognosis and sensitivity to immunotherapy in lung adenocarcinoma. Genome Med. 2022;14(1):5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Chen J, Tan Y, Sun F, Hou L, Zhang C, Ge T, et al. Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer. Genome Biol. 2020;21(1):152.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. He D, Wang D, Lu P, Yang N, Xue Z, Zhu X, et al. Single-cell RNA sequencing reveals heterogeneous tumor and immune cell populations in early-stage lung adenocarcinomas harboring EGFR mutations. Oncogene. 2021;40(2):355–68.

    Article  CAS  PubMed  Google Scholar 

  14. Cancer Genome Atlas Research N. Comprehensive molecular profiling of lung adenocarcinoma. Nature. 2014;511(7511):543–50.

    Article  Google Scholar 

  15. Chen J, Yang H, Teo ASM, Amer LB, Sherbaf FG, Tan CQ, et al. Genomic landscape of lung adenocarcinoma in east Asians. Nat Genet. 2020;52(2):177–86.

    Article  CAS  PubMed  Google Scholar 

  16. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2(5):401–4.

    Article  PubMed  Google Scholar 

  17. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6(269):pl1.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Gillette MA, Satpathy S, Cao S, Dhanasekaran SM, Vasaikar SV, Krug K, et al. Proteogenomic characterization reveals therapeutic vulnerabilities in lung adenocarcinoma. Cell. 2020;182(1):200-25 e35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kim N, Kim HK, Lee K, Hong Y, Cho JH, Choi JW, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun. 2020;11(1):2285.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell. 2019;177(7):1888-902.e21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573-87.e29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Xing X, Yang F, Huang Q, Guo H, Li J, Qiu M, et al. Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing. Sci Adv. 2021. https://doi.org/10.1126/sciadv.abd9738.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Buffa FM, Harris AL, West CM, Miller CJ. Large meta-analysis of multiple cancers reveals a common, compact and highly prognostic hypoxia metagene. Br J Cancer. 2010;102(2):428–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Liu Y, He S, Wang XL, Peng W, Chen QY, Chi DM, et al. Tumour heterogeneity and intercellular networks of nasopharyngeal carcinoma at single cell resolution. Nat Commun. 2021;12(1):741.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Piersma B, Hayward MK, Weaver VM. Fibrosis and cancer: a strained relationship. Biochim Biophys Acta Rev Cancer. 2020;1873(2): 188356.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. 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 Res. 2019;47(W1):W556–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Wu Y, Yang S, Ma J, Chen Z, Song G, Rao D, et al. Spatiotemporal immune landscape of colorectal cancer liver metastasis at single-cell level. Cancer Discov. 2022;12(1):134–53.

    Article  CAS  PubMed  Google Scholar 

  31. DeTomaso D, Jones MG, Subramaniam M, Ashuach T, Ye CJ, Yosef N. Functional interpretation of single cell similarity maps. Nat Commun. 2019;10(1):4376.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9(7):671–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Yang SM, Chen LW, Wang HJ, Chen LR, Lor KL, Chen YC, et al. Extraction of radiomic values from lung adenocarcinoma with near-pure subtypes in the international association for the study of lung cancer/the american thoracic society/the european respiratory society (IASLC/ATS/ERS) classification. Lung cancer. 2018;119:56–63.

    Article  PubMed  Google Scholar 

  34. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013;14:7.

    Article  Google Scholar 

  35. Wang B, Zhao Q, Zhang Y, Liu Z, Zheng Z, Liu S, et al. Targeting hypoxia in the tumor microenvironment: a potential strategy to improve cancer immunotherapy. J Exp Clin Cancer Res. 2021;40(1):24.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. 2021;39(6):845-65.e7.

    Article  CAS  PubMed  Google Scholar 

  37. Lequeux A, Noman MZ, Xiao M, Van Moer K, Hasmim M, Benoit A, et al. Targeting HIF-1 alpha transcriptional activity drives cytotoxic immune effector cells into melanoma and improves combination immunotherapy. Oncogene. 2021;40(28):4725–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Martinez-Reyes I, Chandel NS. Cancer metabolism: looking forward. Nat Rev Cancer. 2021;21(10):669–80.

    Article  CAS  PubMed  Google Scholar 

  39. Du X, Song H, Shen N, Hua R, Yang G. The molecular basis of ubiquitin-conjugating enzymes (E2s) as a potential target for cancer therapy. Int J Mol Sci. 2021;22(7):3440.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Pan YH, Yang M, Liu LP, Wu DC, Li MY, Su SG. UBE2S enhances the ubiquitination of p53 and exerts oncogenic activities in hepatocellular carcinoma. Biochem Biophys Res Commun. 2018;503(2):895–902.

    Article  CAS  PubMed  Google Scholar 

  41. Jung CR, Hwang KS, Yoo J, Cho WK, Kim JM, Kim WH, et al. E2-EPF UCP targets pVHL for degradation and associates with tumor growth and metastasis. Nat Med. 2006;12(7):809–16.

    Article  CAS  PubMed  Google Scholar 

  42. Liu W, Xin H, Eckert DT, Brown JA, Gnarra JR. Hypoxia and cell cycle regulation of the von Hippel-Lindau tumor suppressor. Oncogene. 2011;30(1):21–31.

    Article  CAS  PubMed  Google Scholar 

  43. Sun T, Liu Z, Yang Q. The role of ubiquitination and deubiquitination in cancer metabolism. Mol Cancer. 2020;19(1):146.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Yang C, Zhang J, Liao M, Yang Y, Wang Y, Yuan Y, et al. Folate-mediated one-carbon metabolism: a targeting strategy in cancer therapy. Drug Discov Today. 2021;26(3):817–25.

    Article  CAS  PubMed  Google Scholar 

  45. Shang M, Yang H, Yang R, Chen T, Fu Y, Li Y, et al. The folate cycle enzyme MTHFD2 induces cancer immune evasion through PD-L1 up-regulation. Nat Commun. 2021;12(1):1940.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Yu C, Yang L, Cai M, Zhou F, Xiao S, Li Y, et al. Down-regulation of MTHFD2 inhibits NSCLC progression by suppressing cycle-related genes. J Cell Mol Med. 2020;24(2):1568–77.

    Article  CAS  PubMed  Google Scholar 

  47. Harris IS, DeNicola GM. The complex interplay between antioxidants and ROS in cancer. Trends Cell Biol. 2020;30(6):440–51.

    Article  CAS  PubMed  Google Scholar 

  48. Sahoo BM, Banik BK, Borah P, Jain A. Reactive oxygen species (ROS): key components in cancer therapies. Anticancer Agents Med Chem. 2022;22(2):215–22.

    Article  PubMed  Google Scholar 

  49. Lee CQE, Kerouanton B, Chothani S, Zhang S, Chen Y, Mantri CK, et al. Coding and non-coding roles of MOCCI (C15ORF48) coordinate to regulate host inflammation and immunity. Nat Commun. 2021;12(1):2130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Donnelly RP, Finlay DK. Glucose, glycolysis and lymphocyte responses. Mol Immunol. 2015;68(2 Pt C):513–9.

    Article  CAS  PubMed  Google Scholar 

  51. Caushi JX, Zhang J, Ji Z, Vaghasia A, Zhang B, Hsiue EH, et al. Transcriptional programs of neoantigen-specific TIL in anti-PD-1-treated lung cancers. Nature. 2021;596(7870):126–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Ren X, Zhang L, Zhang Y, Li Z, Siemers N, Zhang Z. Insights gained from single-cell analysis of immune cells in the tumor microenvironment. Annu Rev Immunol. 2021;39:583–609.

    Article  CAS  PubMed  Google Scholar 

  53. Guo X, Zhang Y, Zheng L, Zheng C, Song J, Zhang Q, et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat Med. 2018;24(7):978–85.

    Article  CAS  PubMed  Google Scholar 

  54. van der Leun AM, Thommen DS, Schumacher TN. CD8(+) T cell states in human cancer: insights from single-cell analysis. Nat Rev Cancer. 2020;20(4):218–32.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Geltink RIK, Kyle RL, Pearce EL. Unraveling the complex interplay between T cell metabolism and function. Annu Rev Immunol. 2018;36:461–88.

    Article  CAS  PubMed  Google Scholar 

  56. Xia L, Oyang L, Lin J, Tan S, Han Y, Wu N, et al. The cancer metabolic reprogramming and immune response. Mol Cancer. 2021;20(1):28.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Borst J, Ahrends T, Babala N, Melief CJM, Kastenmuller W. CD4(+) T cell help in cancer immunology and immunotherapy. Nat Rev Immunol. 2018;18(10):635–47.

    Article  CAS  PubMed  Google Scholar 

  58. Mould KJ, Moore CM, McManus SA, McCubbrey AL, McClendon JD, Griesmer CL, et al. Airspace macrophages and monocytes exist in transcriptionally distinct subsets in healthy adults. Am J Respir Crit Care Med. 2021;203(8):946–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Cheng S, Li Z, Gao R, Xing B, Gao Y, Yang Y, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell. 2021;184(3):792-809.e23.

    Article  CAS  PubMed  Google Scholar 

  60. Tie Y, Zheng H, He Z, Yang J, Shao B, Liu L, et al. Targeting folate receptor beta positive tumor-associated macrophages in lung cancer with a folate-modified liposomal complex. Signal Transduct Target Ther. 2020;5(1):6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Chen J, Yao Y, Gong C, Yu F, Su S, Chen J, et al. CCL18 from tumor-associated macrophages promotes breast cancer metastasis via PITPNM3. Cancer Cell. 2011;19(4):541–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Cardoso AP, Pinto ML, Castro F, Costa AM, Marques-Magalhaes A, Canha-Borges A, et al. The immunosuppressive and pro-tumor functions of CCL18 at the tumor microenvironment. Cytokine Growth Factor Rev. 2021;60:107–19.

    Article  CAS  PubMed  Google Scholar 

  63. Jhunjhunwala S, Hammer C, Delamarre L. Antigen presentation in cancer: insights into tumour immunogenicity and immune evasion. Nat Rev Cancer. 2021;21(5):298–312.

    Article  CAS  PubMed  Google Scholar 

  64. Grauel AL, Nguyen B, Ruddy D, Laszewski T, Schwartz S, Chang J, et al. TGFbeta-blockade uncovers stromal plasticity in tumors by revealing the existence of a subset of interferon-licensed fibroblasts. Nat Commun. 2020;11(1):6315.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Chen Y, McAndrews KM, Kalluri R. Clinical and therapeutic relevance of cancer-associated fibroblasts. Nat Rev Clin Oncol. 2021;18(12):792–804.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Hussain A, Voisin V, Poon S, Karamboulas C, Bui NHB, Meens J, et al. Distinct fibroblast functional states drive clinical outcomes in ovarian cancer and are regulated by TCF21. J Exp Med. 2020. https://doi.org/10.1084/jem.20191094.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Biffi G, Tuveson DA. Diversity and biology of cancer-associated fibroblasts. Physiol Rev. 2021;101(1):147–76.

    Article  CAS  PubMed  Google Scholar 

  68. Carman CV, Martinelli R. T Lymphocyte-endothelial interactions: emerging understanding of trafficking and antigen-specific immunity. Front Immunol. 2015;6:603.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Lambrechts D, Wauters E, Boeckx B, Aibar S, Nittner D, Burton O, et al. Phenotype molding of stromal cells in the lung tumor microenvironment. Nat Med. 2018;24(8):1277–89.

    Article  CAS  PubMed  Google Scholar 

  70. Liang P, Ballou B, Lv X, Si W, Bruchez MP, Huang W, et al. Monotherapy and combination therapy using anti-angiogenic nanoagents to fight cancer. Adv Mater. 2021;33(15): e2005155.

    Article  PubMed  Google Scholar 

  71. Huinen ZR, Huijbers EJM, van Beijnum JR, Nowak-Sliwinska P, Griffioen AW. Anti-angiogenic agents - overcoming tumour endothelial cell anergy and improving immunotherapy outcomes. Nat Rev Clin Oncol. 2021;18(8):527–40.

    Article  PubMed  Google Scholar 

  72. Li C, Guo L, Li S, Hua K. Single-cell transcriptomics reveals the landscape of intra-tumoral heterogeneity and transcriptional activities of ECs in CC. Mol Ther Nucleic Acids. 2021;24:682–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science. 2020;368(6487):eaaw5473.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Chang CH, Qiu J, O’Sullivan D, Buck MD, Noguchi T, Curtis JD, et al. Metabolic competition in the tumor microenvironment is a driver of cancer progression. Cell. 2015;162(6):1229–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Vitale I, Manic G, Coussens LM, Kroemer G, Galluzzi L. Macrophages and metabolism in the tumor microenvironment. Cell Metab. 2019;30(1):36–50.

    Article  CAS  PubMed  Google Scholar 

  76. Leone RD, Powell JD. Metabolism of immune cells in cancer. Nat Rev Cancer. 2020;20(9):516–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Martinez-Usatorre A, Kadioglu E, Boivin G, Cianciaruso C, Guichard A, Torchia B, et al. Overcoming microenvironmental resistance to PD-1 blockade in genetically engineered lung cancer models. Sci Transl Med. 2021. https://doi.org/10.1126/scitranslmed.abd1616.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Zappasodi R, Serganova I, Cohen IJ, Maeda M, Shindo M, Senbabaoglu Y, et al. CTLA-4 blockade drives loss of T(reg) stability in glycolysis-low tumours. Nature. 2021;591(7851):652–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Canale FP, Basso C, Antonini G, Perotti M, Li N, Sokolovska A, et al. Metabolic modulation of tumours with engineered bacteria for immunotherapy. Nature. 2021;598(7882):662–6.

    Article  CAS  PubMed  Google Scholar 

  80. Kory N, Wyant GA, Prakash G, Uit de Bos J, Bottanelli F, Pacold ME, et al. SFXN1 is a mitochondrial serine transporter required for one-carbon metabolism. Science. 2018. https://doi.org/10.1126/science.aat9528.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Nilsson R, Nicolaidou V, Koufaris C. Mitochondrial MTHFD isozymes display distinct expression, regulation, and association with cancer. Gene. 2019;716: 144032.

    Article  CAS  PubMed  Google Scholar 

  82. Bonagas N, Gustafsson NMS, Henriksson M, Marttila P, Gustafsson R, Wiita E, et al. Pharmacological targeting of MTHFD2 suppresses acute myeloid leukemia by inducing thymidine depletion and replication stress. Nat Cancer. 2022;3(2):156–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Mo J, Gao Z, Zheng L, Yan M, Xue M, Xu J, et al. Targeting mitochondrial one-carbon enzyme MTHFD2 together with pemetrexed confers therapeutic advantages in lung adenocarcinoma. Cell Death Discov. 2022;8(1):307.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Scaletti ER, Gustafsson Westergren R, Andersson Y, Wiita E, Henriksson M, Homan EJ, et al. The first structure of human mthfd2l and its implications for the development of isoform-selective Inhibitors. ChemMedChem. 2022. https://doi.org/10.1002/cmdc.202200274.

    Article  PubMed  Google Scholar 

  85. Lloyd MC, Cunningham JJ, Bui MM, Gillies RJ, Brown JS, Gatenby RA. Darwinian dynamics of intratumoral heterogeneity: not solely random mutations but also variable environmental selection forces. Can Res. 2016;76(11):3136–44.

    Article  CAS  Google Scholar 

  86. Chen X, Song E. Turning foes to friends: targeting cancer-associated fibroblasts. Nat Rev Drug Discov. 2019;18(2):99–115.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank all the participants in this study.

Funding

This research was funded by the National Science Foundation of China (Grant No.82125001, No.81972172), Clinical Research Plan of SHDC (Grant No. SHDC2020CR2020B), the Shanghai Science and Technology Committee (Grant No. 201409001000), the Shanghai Sailing Program (21YF1438300) and Funding of Shanghai Pulmonary Hospital (Grant No.FKCX1904, No.FKLY20004).

Author information

Authors and Affiliations

Authors

Contributions

PZ and DL conceived the study. JH and YY contributed to the data collection. DL and JH conducted the data analyses and interpreted the results. GJ, LZ, and HY provided administrative support and gave critical advice. LS provided help for statistic tests. LH and SL collected clinical samples and conducted the immunohistochemical analysis. DL and SL contributed to the design and writing of the manuscript. HY, JH, LZ, and PZ participated in the revision of the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Likun Hou, Lele Zhang or Peng Zhang.

Ethics declarations

Ethics approval and consent to participate

The study was conducted in accordance with the principles of the Declaration of Helsinki, and the study protocol was approved by the ethics committee of Shanghai Pulmonary Hospital. Because of the retrospective nature of the study, patient consent for inclusion was waived.

Consent for publication

Not applicable.

Competing interests

The authors declare no potential conflicts of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Characteristics of the LUAD samples included in this study.

Additional file 2: Table S2.

Data of epithelial cells analysis. (1) Differential expression genes between epithelial cells from solid and acinar samples; (2) Differentially enriched HALLMARK signatures in epithelial cells from solid and acinar samples; (3) Differentially enriched KEGG pathways in epithelial cells from solid and acinar samples; (4) Gene list of signatures for ssGSEA in bulk RNA-seq dataset of the TCGA LUAD cohort.

Additional file 3: Table S3.

Markers used for differing between primary immune cell types.

Additional file 4: Table S4.

Data of T/NK cells analysis. (1) Differentially expressed genes between CD8+ T cells from solid and acinar samples; (2) Differentially enriched KEGG pathways between CD8+ T cells from solid and acinar samples; (3) Gene list of signatures for T/NK cells.

Additional file 5: Table S5.

Data of myeloid cells analysis. (1) Gene list of signatures for macrophages; (2) Differentially expressed genes between macrophages from solid and acinar samples; (3) Gene list of signatures for DCs; (4) Differentially expressed genes between mregDCs from solid and acinar samples.

Additional file 6: Table S6.

Gene list of signatures for fibroblasts.

Additional file 7: Fig. S1.

Primary classification of epitheliums, immune cells and stromal cells. AB. Representative hematoxylin–eosin (HE) staining images of LUAD manifesting as acinar (A) and solid (B) growth pattern. The acinar pattern is predominantly glandular, round to oval in shape, with a central lumen surrounded by tumor cells. While the solid pattern consists of cytoplasm-rich polygonal tumor cells forming dense sheets and lacking any other recognizable patterns. The box regions in the upper panel are shown at higher magnification below. Scale bars, 200 μm (top panels) and 50 μm (lower panels). C. UMAP plots depicting all cells labeled as epitheliums, immune cells or stromal cells and split by sample types. D. UMAP plots of canonical markers for labeling general cell types. E. Heatmaps showing large-scale CNVs for individual epitheliums from tumor samples. Each row represented a cell and the columns represented chromosomal regions. Stromal cells were treated as references (top) and large-scale CNVs were observed in tumor cells (bottom).

Additional file 8: Fig. S2.

Characteristics of epithelial cells from solid and acinar samples. AB. Boxplots showing HIF1A and LDHA mRNA expression across LUAD histologic subtypes in the EAS (A) and the CPTAC cohorts (B). Box centerlines, median; box limits, the 25th and 75th percentiles; box whiskers, 1.5× the interquartile range. For all comparisons of molecular expression between histologic subtypes, the statistical significance was determined by two-sided Wilcoxon rank-sum test (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, n.s not significant). CD. UBE2S and UBE2C mRNA expression in the the EAS cohort (C), and protein expression in the CPTAC cohort (D). Comparisons were performed using two-sided Wilcoxon rank-sum test (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, n.s not significant).

Additional file 9: Fig. S3.

Metabolic differences between epithelial cells from solid and acinar samples. A. The average optical density of MTHFD2 immunohistochemical staining in tumor regions from different histologic patterns as semi-quantified by the Image J software. In the box plot, the centre line represents the median, box edges show the 25th and 75th percentiles, and whiskers extend to 1.5× the interquartile range. The statistical significance was determined by two-sided Wilcoxon rank-sum test. B. Kaplan–Meier survival curves showing the prognostic difference between the low and high MTHFD2 expression groups in the TCGA LUAD cohort. C. Correlation between the expression of MTHFD2 and UBE2S in the TCGA LUAD cohort. P-value was determined by Pearson's correlation test. DF. Violin plots showing enrichment scores of one-carbon pool by folate (D), pyrimidine metabolism (E) and galactose metabolism (F) signatures by histologic subtypes in the TCGA LUAD cohort. Global differences were measured by the Kruskal-Wallis test.

Additional file 10: Fig. S4.

Primary immune cell types identification and T/NK cells analysis. AC. UMAP plots of all immune cells colored by major immune types (A), sample origins (B) and histologic subgroups (C). D. UMAP plots of selected canonical markers for identifying major immune types. EF. UMAP plots of T/NK cells colored by sample origins (E) and histologic subgroups (F). G. Dot plots showing the top marker genes for each T/NK cell cluster.H. Volcano plot showing differential expression genes between CD8+ T cells from solid and acinar samples. I. Differentially enriched KEGG pathways between CD4+ T cells from solid and acinar samples revealed by GSVA.

Additional file 11: Fig. S5.

Myeloid cells analysis. A. UMAP plot of annotated myeloid cells. B. UMAP plots of selected canonical markers for annotating myeloid cells. C. UMAP plot of monocytes and macrophages colored by histologic subgroups. D. Dot plots showing the top marker genes for each monocyte/macrophage cluster. E. Gene ontology annotation of marker genes for the Macro-C2 subset. F. Gene ontology annotation of marker genes for the Macro-C3 subset. G–I. Violin plots showing enrichment scores of macrophage and DC traffic (G), M1 phenotype (H) and immune suppression by myeloid cells (I) signatures by histologic subtypes in the TCGA LUAD cohort. Global differences were measured by the Kruskal-Wallis test.

Additional file 12: Fig. S6.

Fibroblasts analysis. A. UMAP plot of fibroblasts colored by histologic subgroups. B. UMAP plots showing the expression of selected canonical marker genes for fibroblasts. C. Dot plots showing the top marker genes for each fibroblast cluster. D. Gene ontology annotation of marker genes for the Fibro-C1 subset. EG. Violin plots showing enrichment scores of fibrillar collagens (E), matrix remodeling (F) and EMT (G) signatures by LUAD histologic subtypes in the TCGA cohort. Global differences were measured by the Kruskal-Wallis test.

Additional file 13: Fig. S7.

Endothelial cells analysis. A. UMAP plot of endothelial cells colored by histologic subgroups. B. UMAP plots showing the expression of selected canonical marker genes for endothelial cells. C. Dot plots showing the top marker genes for each endothelium cluster.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, D., Yu, H., Hu, J. et al. Comparative profiling of single-cell transcriptome reveals heterogeneity of tumor microenvironment between solid and acinar lung adenocarcinoma. J Transl Med 20, 423 (2022). https://doi.org/10.1186/s12967-022-03620-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12967-022-03620-3

Keywords