Research Paper Volume 15, Issue 20 pp 11571—11587

Metabolic and senescence characteristics associated with the immune microenvironment in non-small cell lung cancer: insights from single-cell RNA sequencing

Hongliang Liao1, *, , Zihao Wan2, *, , Yaqin Liang3, *, , Lin Kang4, , Renping Wan1, ,

  • 1 Department of Thoracic Surgery, The Yuebei People’s Hospital of Shaoguan, Shaoguan, Guangdong 512025, China
  • 2 College of Physical Education and Health, Guangzhou University of Chinese Medicine, Guangzhou, Guangdong, China
  • 3 Department of Nursing Medical College, Shaoguan University, Shaoguan, Guangdong 512005, China
  • 4 Department of Gynaecology and Obstetrics, The Qujiang District Maternal and Child Health Care Hospital, Shaoguan, Guangdong, China
* Equal contribution

Received: June 12, 2023       Accepted: October 6, 2023       Published: October 26, 2023      

https://doi.org/10.18632/aging.205146
How to Cite

Copyright: © 2023 Liao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Non-small lung cancer (NSCLC) has been defined as a highly life-threatening heterogeneous disease, with high mortality and occurrence. Recent research has indicated that tumor-infiltrating lymphocytes play a key determinant role in cancer progression. Emerging single-cell RNA sequencing (also termed scRNA-seq) has been extensively applied to depict the baseline landscape of the cell composition and function phenotype in the tumor environment (TME). Herein, we dissected the cell types in NSCLC samples (including tissue and blood) and identified three types of cell marker genes including cancer cells, T cells, and macrophages by integrating two NSCLC-associated scRNA-seq datasets in GEO. Survival analysis indicated that 17 marker genes were related to tumor prognosis. Function annotation was used to scrutinize the molecular mechanism of these marker genes in different cells. Besides, we investigated the developmental trajectory and T cell receptor repertoire diversity of tumor-infiltrating T cells. Our analysis will help further understand the complexity of cell components and the heterogeneity of TME in NSCLC.

Introduction

Non-small cell lung cancer (NSCLC) is the leading life-threatening cancer worldwide and accounts for approximately 85% of lung cancers [1, 2]. NSCLC tumors generally harbor extensive genomic variations, with broad alternation load associating with better response to checkpoint blockade therapies, albeit with some exceptions [3, 4]. With the rapid development of surgical and chemotherapeutical therapies, overall survival has improved considerably; however, clinical outcomes in advanced NSCLC patients remain unfavorable [5]. The past several decades have witnessed a revolution in cancer treatment by moving away from target tumors typically (e.g., radiotherapy, and chemotherapy), although with great advancement in improving the prognosis of NSCLC, toward antibody-based immunotherapy that nudges immune responses against tumors [6, 7].

Recent advances in single-cell sequencing technology have provided a powerful tool for basic research and hold significant potential for deciphering biological systems with unprecedented resolution [8]. Compared with traditional bulk sequencing, the most unique feature of single-cell sequencing lies in the analysis of intercellular heterogeneity and the requirement of obtaining single cells in good status [9, 10]. Sequencing the genome of individual cells, encompassing transcriptomic and epigenomic assays, is promising to reveal somatic mutations and comprehensively characterize the diversity of tumor cell types [11]. The lung tumor microenvironment, comprised of tumor bulk plus supporting cells, plays a critical role in the onset, progression, and malignant transformation of NSCLC and is tightly correlated to therapeutic outcomes [12]. Single-cell DNA sequencing can delineate the genomic variations within a single cell, but it fails to identify salient expression differences in heterogeneous cells [13, 14]. However, the generation of transcriptional profiles of single cells by single-cell RNA sequencing (scRNA-seq) enables the investigation of intercellular heterogeneity on a transcriptome-wide and single nucleotide and level [15]. scRNA-seq of NSCLC samples can disclose the intrinsic gene regulatory mechanisms that determine cellular properties and reveal the developmental and evolutionary relationships of NSCLC cell populations. In addition, it contributes to dissect the correlation between NSCLC heterogeneity, signaling pathways, drug resistance and microenvironment shaping, which is of great importance for NSCLC treatment.

In the current study, we investigated the cell types and cell features genes in TME of NSCLC by integrating the two scRNA datasets of NSCLC from GEO. Also, we studied the association between cell marker genes and tumor survival. Of note, we analyzed the single-cell trajectory features and characteristic differences of T cell receptors (TCR) in NSCLC. Our research highlighted the decisive function of the cell component of TME and emphasized the presence of TME heterogeneity in tumor progression.

Materials and Methods

Data collection and processing

The original scRNA-seq gene expression matrix information related to NSCLC was downloaded from the GEO database. We selected GSE162498 with 17 samples (13 tumor samples and 4 blood samples) and GSE117570 with 8 samples (4 tumor samples and 4 paracancerous samples) as our metadata. The single cell TCR-seq (T cell receptor) annotation file was acquired from GSE162498 with 10 samples (6 tumor samples, 2 para-tumor normal samples, and 4 blood samples). The workflow of data processing was performed according to previous research [16]. Generally, the cells that meet the requirements are screened out according to the following standard: min. features = 100, min. cells = 10, nFeature RNA > 100 and nFeature RNA < 5000, percent.mito < 20, nCount RNA > 10. There was a batch effect between GSE117570 and GSE162498. Therefore, we removed the batch effect by using FindAnchor and IntegrateData functions in the Seurat package (V4) in R software [17]. By integrating GSE117570 and GSE162498 data, a total of 23886 genes in 153,595 cells were identified.

Cell cluster analysis

Principal component analysis (PCA) is a linear method used to reduce data complex dimensionality. PCA analysis is to condense the information of a large number of genes in the data into a few variables (representing the main effects in the sample). Then, just 2–3 variables (named PC1, PC2, and PC3) can represent most of the information contained in tens of thousands of genes. Then the difference in expression between cells is reflected in the differences in the values of PC1 and PC2 [18]. t-distributed Stochastic Neighbor Embedding (t-SNE) algorithm is a machine learning algorithm for dimensionality reduction, which is used for exploring high-dimensional data [19]. Uniform Manifold Approximation and Projection (UMAP), is an emerging dimensionality reduction algorithm, which can retain the features of the original data to the maximum extent and reduce the feature dimension [20]. We normalize the scRNA-seq data by using the Seurat package. After the dimension of PCA, a total of 15w cells were clustered into 29 clusters by tSNE/UMAP analysis. Among that, the number of cells in the tumor sample was the highest (113448 cells) and normal was the lowest (6432 cells).

Identification of cell marker genes

The specific gene list expressed in each cluster was collected by using the FindAllMarkers function in the Seurat package. Herein, the cell type in each cluster was annotated by the combination of the singleR (3.17) method and manual annotation (Supplementary Table 1). The marker gene annotated manually comes from this article [21] and the immune cell marker website (http://www.cellsignal.com/pathways/immune-cell-markers-human). Also, these cluster-specific genes identified by FindAllMarkers were regarded as cell type-specific genes of annotated cells in the matching cluster. In each cell type, the differentially expressed genes (DEGs) among different groups (tumor, normal, and blood samples) were acquired by using the limma and DESeq method with the threshold of p-value < 0.05. The overlapping genes between cell-specific genes and DEGs were defined as the matching cell marker genes.

Survival analysis and gene expression pattern of marker genes

The influence of cell marker genes on NSCLC prognosis was assessed based on a Cox model, which was created to estimate the survival risk of patients under different conditions. Kaplan-Meier curve analysis was performed to compare the survival differences in patients with different marker gene expression contexts by using the survival (3.5.5) package in R.

Function enrichment and ssGSEA analysis

GO and KEGG enrichment analysis of cell marker genes was used to identify the major biological function and involved signaling pathways of each type of cell. In addition, we also conducted the gene set enrichment analysis to further explore the predominant signaling axis implanted by each type of cell according to previous research. ssGSEA was termed as a non-parametric and unsupervised approach, which was used to calculate the normalized enrichment score of T cell subsets in our metadata with reference to the marker genes of 17 T cell subsets.

Hallmark of cancer and immune signature analysis

We downloaded the hallmark gene set and immunologic signature gene set from the MSigDB database (http://www.gsea-msigdb.org/gsea/msigdb/). The GSEA method was used to analyze the predominant molecular signals of cell marker genes in tumor-associated pathways and immune-related pathways.

Trajectory analysis

We arranged T cells according to pseudotime along a trajectory with the orderCells function by using the monocle3 method. Then we conducted a data dimensionality analysis of cell clustering by using the DDRTree approach. After that, the trajectory image of T cells was visualized by UMAP. Here we defined the naïve T cell as the root node. Likewise, the trajectories of myeloid cells are analyzed with monocle3, and the node at the bifurcation is selected as the root node.

TCR analysis

The T-cell receptor sequencing (TCR-seq) technique was applied to comprehensively and rapidly detect the TCR diversity of antigen recognition decisive surface molecules. When the immune response occurs, the gene rearrangement in the VDJ region of TCR formed many different immune receptors, that is, the immune library. The main purpose of TCR-seq data analysis was to count the occurrence frequency of genes in different regions, namely geneUsage. Here we analyzed the distribution pattern of TCR by using the immunarch package in R based on GSE16249. The results were visualized by the vis function in ggplot2 (3.4.3).

Availability of data and materials

The datasets used and/or analyzed in the present study are available from the corresponding author on reasonable request. R code script was shown in Supplementary File 1.

Results

The cell distribution pattern in NSCLC

We identified a total of 23886 genes in 153,595 cells by integrating the GSE117570 and GSE162498 data with Seurat algorism. According to tSNE/UMAP analysis, these cells were mainly enriched into 29 distinct clusters (Figure 1A). Also, the cell distribution pattern in three sample types (tumor tissues, corresponding para-tumor normal tissues, and blood samples) was visualized. Most cells were mainly augmented in tumor tissue samples while the corresponding para-tumor normal tissue samples gained the lowest cell abundance (Figure 1B, 1C). The overall cell could be annotated into 12 cell types (Figure 1D). The cellular compositions, including cancer cells, T cells, B cells, endothelial cells, and so on (a total of 12 types of cells) in each sample type differed from other types. We found that blood samples had the highest abundance of T cells with the highest percentage of CD4+ T cells. Tumor tissue carried most cellular components such as cancer cells, T cells, B cells, and endothelial cells, among which macrophage had the highest percentages, following CD4+ and CD8+ T cells (Figure 1E). Furthermore, we found that CD4+ and CD8+ T cells were evenly distributed in tumor samples compared with the normal group (Figure 1F). Consistently, the macrophage was the most abundant of three myeloid cells (macrophage, monocyte, and M2-macrophage) in tumor tissues compared to the normal tissue (Figure 1G). These results suggested the active anti-cancer response in tumor tissues.

The cell distribution pattern in NSCLC. (A) The cell cluster analysis of scRNA data. (B, C) The cell distribution in three sample types (tumor tissues, para-tumor normal tissues, and blood samples). (D) Cell annotation results of UMAP cluster analysis. (E) Percentages of clustered total cells in three types of samples including blood, normal tissue and tumor tissue group. (F) The cell percentages of clustered T cells in three types of samples. (G) The cell percentages of clustered macrophage subpopulations in three types of samples.

Figure 1. The cell distribution pattern in NSCLC. (A) The cell cluster analysis of scRNA data. (B, C) The cell distribution in three sample types (tumor tissues, para-tumor normal tissues, and blood samples). (D) Cell annotation results of UMAP cluster analysis. (E) Percentages of clustered total cells in three types of samples including blood, normal tissue and tumor tissue group. (F) The cell percentages of clustered T cells in three types of samples. (G) The cell percentages of clustered macrophage subpopulations in three types of samples.

Gene expression features in cell clusters

The heat map of the top 10 DEGs in each cluster was shown in Figure 2A. Interestingly, the top 10 DEGs identified by cluster 0 were distributed at cluster 0 other than resident clusters (Figure 2B), validating the efficacy of our cell cluster analysis. Observing the abundance composition of cluster 0, we obtained DEGs in cluster 0. GO and KEGG analysis of DEGs in cluster 0 was displayed in Figure 2C, 2D. They were associated with T cell activation and autoimmune disease pathways.

Gene expression features in cell clusters. (A) Heatmap of the top 10 DEGs in 29 clusters. (B) Expression of top 10 DEGs of cluster 0 in 29 clusters. (C, D) GO and KEGG analysis of DEGs in cluster 0.

Figure 2. Gene expression features in cell clusters. (A) Heatmap of the top 10 DEGs in 29 clusters. (B) Expression of top 10 DEGs of cluster 0 in 29 clusters. (C, D) GO and KEGG analysis of DEGs in cluster 0.

Identification of prognostic associated marker genes (PAMGs) of cells

We obtained the DEGs between three samples of each cell type. The top 10 DEGs in were shown in Figure 3A. The overlapping genes between DEGs and cell-specific gene lists were used as marker genes for matching cells. 285 marker genes were identified in the cancer cell, which was consistent with the tumor tissue feature. 142 marker genes were identified in T cells and 225 marker genes in macrophages. Here, we acquired cancer cells, T cells, and macrophage marker genes. Furthermore, we investigated the influence of these cell marker genes on NSCLC survival based on TCGA-LC clinical data. Results indicated that 17 marker genes were found to have significant effects on tumor survival (3 representative genes, KRT6A, NASPA, and ADM Figure 3B3D).

Survival analysis of cell marker genes. (A) DEGs between different groups. (B–D) Survival analysis of 3 representative marker genes (KRT6A, NASPA, and ADM).

Figure 3. Survival analysis of cell marker genes. (A) DEGs between different groups. (BD) Survival analysis of 3 representative marker genes (KRT6A, NASPA, and ADM).

Expression pattern of PAMGs in different types of cells

Considering these PAMGs in tumor prognosis and development, we also explored their expression pattern in our training set (GSE117570 and GSE162498). Gene expression patterns of 17 PAMGs in three sample types were shown in Figure 4A. The heatmap of 17 PAMGs in different types of cells was displayed in Figure 4B. Furthermore, the box plots of gene expression of 3 representative genes, KRT6A, NASPA, and ADM were also displayed in Figure 4C. KRT6A was mainly expressed in the cancer cell. We could conclude that most PAMGs were up-regulated in cancer cells and epithelial cells. These results demonstrated that cancer cells play a crucial role in cancer progression while other cells such as T cells were not directly engaged in tumor development.

Expression pattern of PAMGs in different types of cells. (A) Expression of 17 PAMGs in three sample types. (B) Heatmap of 17 PAMGs in different types of cells. (C) Box plot of 3 representative genes expression (KRT6A, NASPA, and ADM) among tumor, blood, and normal tissues.

Figure 4. Expression pattern of PAMGs in different types of cells. (A) Expression of 17 PAMGs in three sample types. (B) Heatmap of 17 PAMGs in different types of cells. (C) Box plot of 3 representative genes expression (KRT6A, NASPA, and ADM) among tumor, blood, and normal tissues.

Enrichment analysis of marker genes

The GO and KEGG enrichment analysis of cell marker genes in different cell types were performed to examine the functional difference. The most significant pathway in T cells was “response to interferon-gamma”, which was found to be related to antigen recognition (Figure 5A). Recent studies have revealed that lung cancer–specific tumor-infiltrating CD8+ T cells had a distinctive differentiation trajectory, which acquired effector and exhausted phenotypes. This unconventional T cell dysfunction explained the ICB resistance for some NSCLC patients [22]. Of course, IFN-γ has a pivotal influence on ICB efficacy, evidenced by that IFN-γ-related mRNA profile exerted a predictive influence on PD-1 blockade response in a clinical trial setting [23]. On the contrary, the most significant pathway enriched in macrophage was neutrophil activation (Figure 5B), which is related to innate immune activation. Increasing data has demonstrated that there were two lineages of macrophages: local tissue-resident macrophages and those derived from monocytes in NSCLC [24]. Early-stage tumors were influenced by tissue-resident macrophages which enhance tumor cell invasiveness and create a protective shield by triggering a powerful regulatory T cell response [24]. Interestingly, as the tumor expands, monocyte-derived macrophages rise in number, pushing tissue-resident macrophages to periphery of the TME [24]. These data highlighted the potential in focusing on tissue-resident macrophages for early NSCLC interventions. These results indicated that the gene expression pattern of maker genes was strictly governed by the cell functional phenotypes.

Enrichment analysis of marker genes. (A, B) GO and KEGG enrichment analysis of cell marker genes.

Figure 5. Enrichment analysis of marker genes. (A, B) GO and KEGG enrichment analysis of cell marker genes.

Hallmark and immune signature analysis of marker genes

We downloaded the hallmark and immunologic signature and gene set from the MSigDB database to scrutinize the function of our marker genes in tumor-associated and immune-related pathways. By KEGG analysis in tumor signal pathways, we found that B cell maker genes were mainly enriched in estrogen response, interferon-gamma response, androgen response, TNF-α signaling via NK-κB, and complement response (Figure 6A). The immune signature pathway of B cell maker genes was shown in Figure 6B. GSEA results illustrated that the cancer cell marker gene was significantly enriched in xenobiotic metabolism (Figure 6C), suggesting the hyperactive metabolism response of cancer cells. Also, GSEA of immune signature enrichment showed that the macrophage marker gene was mainly engaged in influenza induced by PBMC inactivation (Figure 6D), suggesting macrophages play a key role in innate immunity. This was consistent with our previous results (Figure 4B).

Hallmark and immune signature analysis of marker genes. (A, B) KEGG analysis of B cell marker genes in tumor-associated signal pathways and immune signature pathway. (C, D) GSEA results of cancer cell marker and macrophage marker genes.

Figure 6. Hallmark and immune signature analysis of marker genes. (A, B) KEGG analysis of B cell marker genes in tumor-associated signal pathways and immune signature pathway. (C, D) GSEA results of cancer cell marker and macrophage marker genes.

Identification of T cell subpopulations

Considering the great impact of T cells on tumor progression, we focused on the T cells subpopulation with distinct phenotypes. Here we collected the maker of 17 T cell subpopulations with different phenotype functions (https://www.cellsignal.com/pathways/immune-cell-markers-human). The abundance of various T cell subpopulations in three NSCLC samples was shown in Figure 7A. Most T cells were up-regulated in tumor samples (including tissue and blood) compared with the normal group. Furthermore, we found that the distribution pattern of the 17 T cell subpopulation in CD8+ T cells was similar to CD4+ T cell, with the CD8+ T cell and naïve T cell highest infiltration (Figure 7B7E).

T cell subpopulation analysis. (A–C) The abundance of 17 T cell subpopulations in three NSCLC samples. Comparisons between two groups were calculated based on Student’s t test. (D, E) The abundance of CD4+ and CD8+ T cell subpopulations in 17 T cell subpopulations.

Figure 7. T cell subpopulation analysis. (AC) The abundance of 17 T cell subpopulations in three NSCLC samples. Comparisons between two groups were calculated based on Student’s t test. (D, E) The abundance of CD4+ and CD8+ T cell subpopulations in 17 T cell subpopulations.

Analysis of T cell development trajectory

Trajectory analysis can reshape the dynamic evaluation process of cell development and maturity by focusing on the different stage cells with distinct gene signatures. Here we dissected the trajectory development process of T cells and myeloid cells, respectively by performing correlation analysis between marker genes and developmental trajectories of T cell and myeloid cell via Monocle3. These significant genes were selected for the following trajectory analysis via Monocle3. The developmental trajectory of T cells was shown in Figure 8A, 8B. There was heterogeneity in the differentiation pathway of T cells with distinct phenotypes. Conversely, both M2-macrophage and monocyte evaluation toward macrophage (Figure 8C, 8D). In addition, we identified 57 marker genes that were found to exert vital function in the developmental trajectories of T cells (Figure 8E). We also investigated the expression pattern of their expression in T cells (Figure 8F). Most marker genes were down-regulated in naïve T cells and they were expressed in CD8+ and CD4+ T cells.

T cell trajectory analysis. (A) The evolution pathway of T cells. (B) The pseudo-time sequence of evolution of T cells. (C) The evolution pathway of macrophage. (D) The pseudo-time sequence of evolution of macrophages. (E) Significant marker gene lists involved T cell development. (F) Gene expression pattern of the significant marker gene in T cell.

Figure 8. T cell trajectory analysis. (A) The evolution pathway of T cells. (B) The pseudo-time sequence of evolution of T cells. (C) The evolution pathway of macrophage. (D) The pseudo-time sequence of evolution of macrophages. (E) Significant marker gene lists involved T cell development. (F) Gene expression pattern of the significant marker gene in T cell.

TCR repertoire analysis

TCR belongs to the immunoglobulin superfamily and is a characteristic mark on the surface of T cells. TCR is a heterodimer composed of two different peptide chains, which consisted of TCR1 (γ and δ chains) and TCR2 (α and β chains). The random rearrangement of genes produces TCR polymorphism, forming a huge TCR library. The polymorphism of the TCR indicates the richness of T cell types or subpopulations. Here we downloaded TCR data of 8 NSCLC samples (blood and tissue sample) from GSE162499, including 162,839 TCRs of 42,163 cells. Results suggested that TCR distribution in the NSCLC blood sample was similar to the disposal pattern in NSCLC tumor tissue (Figure 9A, 9B). Most TCRs were identified as T cells rather than cancer cells or B cells, among which CD4+ and CD8+ T cells had the highest percentages while the CD4- and CD8- (double negative, DN T cell) carried the lowest TCRs, suggesting the maturity and function complexity of T cells (Figure 9C, 9D). In addition, we observed that the CD4+ T cell in tissues was higher than that in blood, demonstrating the great potential of tumor infiltration of CD4+ T cells. Notably, there was an additional TCR repertoire of Treg cell subpopulation in tissue, implying that Treg cells engaged in tumor microenvironment reshaping.

TCR repertoire analysis. (A) The distribution of TCR in different cell types. (B) The cell percentages of TCR in tumor samples. (C, D) TCR colony numbers and length in different types cell of tumor samples.

Figure 9. TCR repertoire analysis. (A) The distribution of TCR in different cell types. (B) The cell percentages of TCR in tumor samples. (C, D) TCR colony numbers and length in different types cell of tumor samples.

Discussion

The tumor ecosystem consists of multiple cell components including cancer cells, tumor-infiltrating lymphocytes, stromal cells, endothelial cells, and other supporting ingredients. The cross-talk between cell partners and their inter-communications of themselves created a dynamic and complex local microenvironment manipulating the overall tumor progression. The tumor ecosystem was also termed a tumor microenvironment (TME). More importantly, a plethora of research has manifested that plastic TME with heterogeneity is recognized as a predominant hallmark in cancer immunotherapy. For example, according to a recent study, a platelet protein, TLT-1 was found to recognize the specific CD3ε region of CD8+ T cells, impeding the anti-tumor function of CD8+ T cells by regulating the NF-κB pathway [25]. Enhancing the CD137 expression of T cells boosted the cytotoxic effector function of T cells by activating CD8+ T cells, improving the tumor OS [26]. In addition, anti- PD-1 immunotherapy manipulated CD4+ T cell chemotaxis response by reshaping CD11b+ neutrophil degranulation activity, which was considered a clinical indicator of successful pancreatic ductal adenocarcinoma (PDAC) immunotherapy [26]. Anti-cancer immunotherapy has made a great breakthrough in recent years. Of note is the need to consider the heterogeneity of tumor-infiltrating immune cells. Here, we investigated the overall cell components in the TME of NSCLC. By integrating two NSCLC-associated scRNA-seq datasets (including tumor tissue, para-tumor normal tissue, and blood samples), we identified 29 distinct cell clusters and most cell clusters were mainly enriched in the tumor tissue group. These results indicated that tumor tissue was a complex cell ecosystem, implying the existence of heterogeneity. Because of the high proportion of cluster 0, we obtained significant cluster-specific genes in cluster 0. They were chiefly involved in T cell activation and autoimmune disease-associated pathways.

Furthermore, 285 cancer cell marker genes, 142 T cell marker genes, and 225 macrophage marker genes were identified by taking the intersection between DEGs (tumor vs. normal) and immune cell marker genesets. 17 marker genes were associated with NSCLC survival. KRT6A has been reported to be a negative regulator in NSCLC survival [27], which was in agreement with our analysis. An isolated study based on transcriptome data also validated the efficacy of ADM as an independent prognosis factor [28]. Here we exemplified that ADM was a key clinical factor in evaluating NSCLC survival at the single-cell level. Enrichment analysis indicated that T and B cell marker genes principally engaged in response to interferon-gamma (IFN-γ), which was indispensable to the cancer-killing of CAR-T cell therapy. A recent study suggested that the knockout of genes in the IFN γ receptor signal pathway (IFNGR1, JAK1, or JAK2) endowed solid tumor cells with resistance to CAR-T cells by reducing the overall binding time and affinity of CAR-T cells with the cancer cells [29]. Macrophage marker genes were associated with neutrophil activation. Within the tumor background, neutrophils derive bone marrow-derived suppressor cell (MDSCs) subsets, which suppressed the anti-tumor response of T cells by regulating PD-L1 expression to promote angiogenesis [30, 31]. The cancer cell marker gene was significantly enriched in xenobiotic metabolism. The metabolic preferences acquired by cancer cells were a favorable factor in the nutritional requirements of cancer cell growth and proliferation, validating the profound influence of cell metabolism reprogramming on tumorigenesis [32]. The unique metabolic characteristic of cancer cells endowed them with the ability to obtain necessary nutrients from a nutrient-deficient environment, which was vital to maintain viability and generating new biomass [33].

As previously described, T cells subpopulations with distinct phenotypes varied in anti-tumor activity. Herein, 17 T cell subpopulations were enrolled in our analysis. naïve T (DN T cell) cell was identified as the highest infiltration T cell subpopulations in tumor tissue and blood samples. It is widely accepted that CD3+CD4-CD8- DNT cells are a distinct subset of T cells. DNTs carried an effective cancer-killing potential in NSCLC models. The expanding DNTs in vivo acquired an anti-tumor phenotype with the NK cell marker NKG2D, DNAM-1, which mediated the cytotoxicity response targeting cancer cells by secreting cytotoxic cytokines such as IFNγ and soluble TRAIL (sTRAIL) [34]. Emerging evidence validated that TCRαβ+ DNTs and TCRγδ+ DNTs displayed equivalent anti-tumor activity in a series of preclinical models. There was the DNT cell subpopulation in tumor-infiltrating lymphocytes derived NSCLC samples, together with increased expression of PD-1. Interestingly, the DNT cells extracted from healthy donors curtailed late-stage lung cancer progression through which anti-PD-1 therapy contributed to DNT cell infiltration into primary tumor biomass [35]. These results, combined with our findings, manifested that targeting DNT cells could be a novel platform in immunotherapy.

The cell trajectory analysis based on the correlation between cell marker genes and trajectories-associated genes via Monocle3 indicated the temporal heterogeneity of NSCLC-infiltrating T cells. We found that 57 marker genes were associated with the T cells maturation process. Most marker genes were down-regulated in naïve T cells (considered as immature T cells) and expressed in mature T cells such as CD8+ and CD4+ T cells. These results demonstrated that 57 marker genes were important driving forces in T cell development. Consistent with this point, our TCR analysis indicated that the TCR numbers of early-stage T cells were lower than the mature-stage T cells, including CD4+ and CD8+ T cells. These data showed that the dynamic changes of the TCR repertoire library could be a sensitive indicator of the T cell development process. In addition, our data revealed an additional TCR repertoire of Treg cell subpopulation in tumor tissue.

From a clinical perspective, understanding the behavior and molecular signatures of Treg cells in the TME is paramount. While some studies have demonstrated the immune-suppressive attributes of Treg cells in the TME, which aids tumor growth and immune evasion [36, 37]. However, another study discovered Treg cells restrained tumor growth and increased response to immunotherapy by up taking lactic acid [38]. This dichotomy suggests a nuanced metabolic balance within Treg cells. Clinically, manipulating this balance could offer a dual therapeutic approach: impairing tumor growth and promoting immune response by reshaping the metabolic landscape within the TME.

Given this data, it becomes clear that understanding the intricacies of T cell trajectory and the TCR repertoire can provide critical insights for enhancing patient care. Harnessing this point can guide the development of personalized therapeutic strategies for NSCLC, thereby improving patient outcomes and response rates to treatments.

Conclusion

Taken together, our data revealed the heterogeneity of TME, the T with different phenotypes cell of TME provides a theoretical basis for variable immunotherapy response.

Author Contributions

All the authors contributed substantially to the work presented in this article. Hongliang Liao, Conceptualization, Data curation, Methodology, Validation, Visualization, Writing - original draft; Zihao Wan, Yaqin Liang, Data curation, Visualization, Writing – review and editing; Lin Kang, Validation, Writing – review and editing; Renping Wan, Supervision, Project administration. The corresponding author had full access to all of the data and the final responsibility for the decision to submit this article for publication.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

No funding was provided for this study.

References

  • 1. Lutschg JH. Lung cancer. N Engl J Med. 2009; 360:87–8. https://doi.org/10.1056/NEJMc082208 [PubMed]
  • 2. Herbst RS, Morgensztern D, Boshoff C. The biology and management of non-small cell lung cancer. Nature. 2018; 553:446–54. https://doi.org/10.1038/nature25183 [PubMed]
  • 3. Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, Lee W, Yuan J, Wong P, Ho TS, Miller ML, Rekhtman N, Moreira AL, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015; 348:124–8. https://doi.org/10.1126/science.aaa1348 [PubMed]
  • 4. Topalian SL, Taube JM, Anders RA, Pardoll DM. Mechanism-driven biomarkers to guide immune checkpoint blockade in cancer therapy. Nat Rev Cancer. 2016; 16:275–87. https://doi.org/10.1038/nrc.2016.36 [PubMed]
  • 5. Forde PM, Chaft JE, Smith KN, Anagnostou V, Cottrell TR, Hellmann MD, Zahurak M, Yang SC, Jones DR, Broderick S, Battafarano RJ, Velez MJ, Rekhtman N, et al. Neoadjuvant PD-1 Blockade in Resectable Lung Cancer. N Engl J Med. 2018; 378:1976–86. https://doi.org/10.1056/NEJMoa1716078 [PubMed]
  • 6. Tumeh PC, Harview CL, Yearley JH, Shintaku IP, Taylor EJ, Robert L, Chmielowski B, Spasic M, Henry G, Ciobanu V, West AN, Carmona M, Kivork C, et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature. 2014; 515:568–71. https://doi.org/10.1038/nature13954 [PubMed]
  • 7. Levy BP, Signorovitch JE, Yang H, Patterson-Lomba O, Xiang CQ, Parisi M. Effectiveness of first-line treatments in metastatic squamous non-small-cell lung cancer. Curr Oncol. 2019; 26:e300–8. https://doi.org/10.3747/co.26.4485 [PubMed]
  • 8. Grün D, van Oudenaarden A. Design and Analysis of Single-Cell Sequencing Experiments. Cell. 2015; 163:799–810. https://doi.org/10.1016/j.cell.2015.10.039 [PubMed]
  • 9. Kang K, Wang X, Meng C, He L, Sang X, Zheng Y, Xu H. The application of single-cell sequencing technology in the diagnosis and treatment of hepatocellular carcinoma. Ann Transl Med. 2019; 7:790. https://doi.org/10.21037/atm.2019.11.116 [PubMed]
  • 10. Lei Y, Tang R, Xu J, Wang W, Zhang B, Liu J, Yu X, Shi S. Applications of single-cell sequencing in cancer research: progress and perspectives. J Hematol Oncol. 2021; 14:91. https://doi.org/10.1186/s13045-021-01105-2 [PubMed]
  • 11. Armand EJ, Li J, Xie F, Luo C, Mukamel EA. Single-Cell Sequencing of Brain Cell Transcriptomes and Epigenomes. Neuron. 2021; 109:11–26. https://doi.org/10.1016/j.neuron.2020.12.010 [PubMed]
  • 12. Zhang J, Song C, Tian Y, Yang X. Single-Cell RNA Sequencing in Lung Cancer: Revealing Phenotype Shaping of Stromal Cells in the Microenvironment. Front Immunol. 2022; 12:802080. https://doi.org/10.3389/fimmu.2021.802080 [PubMed]
  • 13. Evrony GD, Hinch AG, Luo C. Applications of Single-Cell DNA Sequencing. Annu Rev Genomics Hum Genet. 2021; 22:171–97. https://doi.org/10.1146/annurev-genom-111320-090436 [PubMed]
  • 14. Li X, Wang CY. From bulk, single-cell to spatial RNA sequencing. Int J Oral Sci. 2021; 13:36. https://doi.org/10.1038/s41368-021-00146-0 [PubMed]
  • 15. Eling N, Richard AC, Richardson S, Marioni JC, Vallejos CA. Correcting the Mean-Variance Dependency for Differential Variability Testing Using Single-Cell RNA Sequencing Data. Cell Syst. 2018; 7:284–94.e12. https://doi.org/10.1016/j.cels.2018.06.011 [PubMed]
  • 16. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, Rodman C, Luo CL, Mroz EA, Emerick KS, Deschler DG, Varvares MA, Mylvaganam R, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017; 171:1611–24.e24. https://doi.org/10.1016/j.cell.2017.10.044 [PubMed]
  • 17. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, Lee MJ, Wilk AJ, Darby C, Zager M, Hoffman P, Stoeckius M, Papalexi E, et al. Integrated analysis of multimodal single-cell data. Cell. 2021; 184:3573–87.e29. https://doi.org/10.1016/j.cell.2021.04.048 [PubMed]
  • 18. Zhao L, Li Z, Vong JSL, Chen X, Lai HM, Yan LYC, Huang J, Sy SKH, Tian X, Huang Y, Chan HYE, So HC, Ng WL, et al. Pharmacologically reversible zonation-dependent endothelial cell transcriptomic changes with neurodegenerative disease associations in the aged brain. Nat Commun. 2020; 11:4413. https://doi.org/10.1038/s41467-020-18249-3 [PubMed]
  • 19. Sugimura R, Jha DK, Han A, Soria-Valles C, da Rocha EL, Lu YF, Goettel JA, Serrao E, Rowe RG, Malleshaiah M, Wong I, Sousa P, Zhu TN, et al. Haematopoietic stem and progenitor cells from human pluripotent stem cells. Nature. 2017; 545:432–8. https://doi.org/10.1038/nature22370 [PubMed]
  • 20. Cao J, Spielmann M, Qiu X, Huang X, Ibrahim DM, Hill AJ, Zhang F, Mundlos S, Christiansen L, Steemers FJ, Trapnell C, Shendure J. The single-cell transcriptional landscape of mammalian organogenesis. Nature. 2019; 566:496–502. https://doi.org/10.1038/s41586-019-0969-x [PubMed]
  • 21. Wu F, Fan J, He Y, Xiong A, Yu J, Li Y, Zhang Y, Zhao W, Zhou F, Li W, Zhang J, Zhang X, Qiao M, et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun. 2021; 12:2540. https://doi.org/10.1038/s41467-021-22801-0 [PubMed]
  • 22. Horton BL, Morgan DM, Momin N, Zagorulya M, Torres-Mejia E, Bhandarkar V, Wittrup KD, Love JC, Spranger S. Lack of CD8+ T cell effector differentiation during priming mediates checkpoint blockade resistance in non-small cell lung cancer. Sci Immunol. 2021; 6:eabi8800. https://doi.org/10.1126/sciimmunol.abi8800 [PubMed]
  • 23. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, Albright A, Cheng JD, Kang SP, Shankaran V, Piha-Paul SA, Yearley J, Seiwert TY, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017; 127:2930–40. https://doi.org/10.1172/JCI91190 [PubMed]
  • 24. Casanova-Acebes M, Dalla E, Leader AM, LeBerichel J, Nikolic J, Morales BM, Brown M, Chang C, Troncoso L, Chen ST, Sastre-Perona A, Park MD, Tabachnikova A, et al. Tissue-resident macrophages provide a pro-tumorigenic niche to early NSCLC cells. Nature. 2021; 595:578–84. https://doi.org/10.1038/s41586-021-03651-8 [PubMed]
  • 25. Tyagi T, Jain K, Yarovinsky TO, Chiorazzi M, Du J, Castro C, Griffin J, Korde A, Martin KA, Takyar SS, Flavell RA, Patel AA, Hwa J. Platelet-derived TLT-1 promotes tumor progression by suppressing CD8+ T cells. J Exp Med. 2023; 220:e20212218. https://doi.org/10.1084/jem.20212218 [PubMed]
  • 26. Li K, Tandurella JA, Gai J, Zhu Q, Lim SJ, Thomas DL 2nd, Xia T, Mo G, Mitchell JT, Montagne J, Lyman M, Danilova LV, Zimmerman JW, et al. Multi-omic analyses of changes in the tumor microenvironment of pancreatic adenocarcinoma following neoadjuvant treatment with anti-PD-1 therapy. Cancer Cell. 2022; 40:1374–91.e7. https://doi.org/10.1016/j.ccell.2022.10.001 [PubMed]
  • 27. Zhou J, Jiang G, Xu E, Zhou J, Liu L, Yang Q. Identification of SRXN1 and KRT6A as Key Genes in Smoking-Related Non-Small-Cell Lung Cancer Through Bioinformatics and Functional Analyses. Front Oncol. 2022; 11:810301. https://doi.org/10.3389/fonc.2021.810301 [PubMed]
  • 28. Liu W, Dong R, Gao S, Shan X, Li M, Yu Z, Sun L. Assessing the Prognostic Capability of Immune-Related Gene Scoring Systems in Lung Adenocarcinoma. J Oncol. 2022; 2022:2151396. https://doi.org/10.1155/2022/2151396 [PubMed]
  • 29. Larson RC, Kann MC, Bailey SR, Haradhvala NJ, Llopis PM, Bouffard AA, Scarfó I, Leick MB, Grauwet K, Berger TR, Stewart K, Anekal PV, Jan M, et al. CAR T cell killing requires the IFNγR pathway in solid but not liquid tumours. Nature. 2022; 604:563–70. https://doi.org/10.1038/s41586-022-04585-5 [PubMed]
  • 30. Matsushima K, Yang D, Oppenheim JJ. Interleukin-8: An evolving chemokine. Cytokine. 2022; 153:155828. https://doi.org/10.1016/j.cyto.2022.155828 [PubMed]
  • 31. Yang J, Jin L, Kim HS, Tian F, Yi Z, Bedi K, Ljungman M, Pasca di Magliano M, Crawford H, Shi J. KDM6A Loss Recruits Tumor-Associated Neutrophils and Promotes Neutrophil Extracellular Trap Formation in Pancreatic Cancer. Cancer Res. 2022; 82:4247–60. https://doi.org/10.1158/0008-5472.CAN-22-0968 [PubMed]
  • 32. Yang M, Liu K, Chen P, Zhu H, Wang J, Huang J. Bromodomain-containing protein 4 (BRD4) as an epigenetic regulator of fatty acid metabolism genes and ferroptosis. Cell Death Dis. 2022; 13:912. https://doi.org/10.1038/s41419-022-05344-0 [PubMed]
  • 33. Stine ZE, Schug ZT, Salvino JM, Dang CV. Targeting cancer metabolism in the era of precision oncology. Nat Rev Drug Discov. 2022; 21:141–62. https://doi.org/10.1038/s41573-021-00339-6 [PubMed]
  • 34. Yao J, Ly D, Dervovic D, Fang L, Lee JB, Kang H, Wang YH, Pham NA, Pan H, Tsao MS, Zhang L. Human double negative T cells target lung cancer via ligand-dependent mechanisms that can be enhanced by IL-15. J Immunother Cancer. 2019; 7:17. https://doi.org/10.1186/s40425-019-0507-2 [PubMed]
  • 35. Fang L, Ly D, Wang SS, Lee JB, Kang H, Xu H, Yao J, Tsao MS, Liu W, Zhang L. Targeting late-stage non-small cell lung cancer with a combination of DNT cellular therapy and PD-1 checkpoint blockade. J Exp Clin Cancer Res. 2019; 38:123. https://doi.org/10.1186/s13046-019-1126-y [PubMed]
  • 36. Nishikawa H, Koyama S. Mechanisms of regulatory T cell infiltration in tumors: implications for innovative immune precision therapies. J Immunother Cancer. 2021; 9:e002591. https://doi.org/10.1136/jitc-2021-002591 [PubMed]
  • 37. Li C, Jiang P, Wei S, Xu X, Wang J. Regulatory T cells in tumor microenvironment: new mechanisms, potential therapeutic strategies and future prospects. Mol Cancer. 2020; 19:116. https://doi.org/10.1186/s12943-020-01234-1 [PubMed]
  • 38. Watson MJ, Vignali PDA, Mullett SJ, Overacre-Delgoffe AE, Peralta RM, Grebinoski S, Menk AV, Rittenhouse NL, DePeaux K, Whetstone RD, Vignali DAA, Hand TW, Poholek AC, et al. Metabolic support of tumour-infiltrating regulatory T cells by lactic acid. Nature. 2021; 591:645–51. https://doi.org/10.1038/s41586-020-03045-2 [PubMed]