Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing

Abstract

Lung adenocarcinomas (LUAD) that radiologically display as subsolid nodules (SSNs) exhibit more indolent biological behavior than solid LUAD. The transcriptomic features and tumor microenvironment (TME) of SSN remain poorly understood. Here, we performed single-cell RNA sequencing analyses of 16 SSN samples, 6 adjacent normal lung tissues (nLung), and 9 primary LUAD with lymph node metastasis (mLUAD). Approximately 0.6 billion unique transcripts were obtained from 118,293 cells. We found that cytotoxic natural killer/T cells were dominant in the TME of SSN, and malignant cells in SSN undergo a strong metabolic reprogram and immune stress. In SSN, the subtype composition of endothelial cells was similar to that in mLUAD, while the subtype distribution of fibroblasts was more like that in nLung. Our study provides single-cell transcriptomic profiling of SSN and their TME. This resource provides deeper insight into the indolent nature of SSN and will be helpful in advancing lung cancer immunotherapy.

Fig. 1. Overview of TME in normal lungs and lung tumors.

(A) Workflow showing the scRNA-seq experimental design and initial data exploration. (B) Cellular populations identified. The UMAP projection of 118,293 single cells from nLung (n = 6), SSN (n = 16), and mLUAD (n = 9) samples shows the formation of 10 main clusters with label names. Each dot corresponds to a single cell, colored according to cell type. (C) Canonical cell markers were used to label clusters by cell identity as represented in the UMAP plot. (D) Average proportion of six main types of CD45⁺ immune cells among nLung, SSN, and mLUAD samples. (E) Percentages of the six types CD45⁺ immune cells among three groups. Y axis: Average percent of samples across the three groups. Groups are shown in different colors. Each bar plot represents one cell cluster. Error bars represent ± SEM for normal and tumor samples. Colored dots represent different samples. All differences with P < 0.05 are indicated; two-sided unpaired Wilcoxon rank sum test was used for analysis. (F) Seven-plex staining panel showing the cellular components of nLung, SSN, and mLUAD tissues.

Fig. 2. Identification and characterization of malignant cells in SSN.

(A) Clustering of 1997 epithelial cells from nLung (n = 6). Each dot corresponds to a single cell, colored according to cell type. (B) Canonical cell markers were used to label epithelial subtypes as represented in the UMAP plot. (C) Sample distribution in each cluster. Each bar corresponds to one cell type cluster, colored according to the samples. (D) Heatmap showing large-scale CNVs for individual cells (rows) from one SSN sample (SSN27) with WES paired data. Nonmalignant cells were treated as references (top), and large-scale CNVs were observed in malignant cells (middle). The CNVs of the sample were validated by WES analysis (bottom). The color shows the log₂ CNV ratio. Red: amplifications; blue: deletions. (E) UMAP projection of 9281 malignant cells from SSN (n = 16) and mLUAD (n = 9). Each dot corresponds to a single cell, colored according to the samples. (F) Top 15 up-regulated hallmark pathways in malignant cells. Top: mLUAD versus SSN. Bottom: SSN versus nLung. (G) Heatmap showing differences in metabolic pathways scored per cell by GSVA between normal epithelial cells in nLung and malignant cells in SSN and mLUAD. (H) Heatmap depicting pairwise correlations of intratumoral programs derived from mLUAD (top) and SSN (bottom). Coherent expression programs are identified and labeled.

Fig. 3. Cytotoxic dominant T and NK cells in SSN.

(A) UMAP projection of 57,301 T and NK cells, showing the composition of 12 main subtypes. (B) UMAP projection of 35,185 T and NK cells derived from SSN. (C) Canonical cell markers were used to identify T/NK cell subtypes. (D) Heatmap of functional gene sets in T and NK clusters. T_reg, regulatory T cell. (E) Cumulative distribution function showing the distribution of naïve (left), cytotoxic (middle), and exhausted (right) state scores in each T/NK subtype. A rightward shift of the curve indicates increased state scores. (F) Average proportion of each subtype between nLung, SSN, and mLUAD. (G) Percentages of each T/NK cell subtype among nLung, SSN, and mLUAD. Y axis: Average percent of samples across the three groups. Groups are shown in different colors. Each bar plot represents one cell cluster. Error bars represent ± SEM for normal and tumor samples. Colored dots represent different samples. All differences with P < 0.05 are indicated; two-sided unpaired Wilcoxon rank sum test was used for analysis. (H) Kaplan-Meier plot showing that patients with LUAD in the TCGA dataset with high expression of CD8-C5 cluster markers have shorter overall survival. The high and low groups are divided by the 75% quantile value of the mean expression of the above gene set. (I) Development trajectory of CD8⁺ T cells inferred by diffusion map, colored by cell subtype and expression of example genes. (J) As in (E), but for “cytotoxic/exhausted score” defined as the average expression level of cytotoxic genes divided by the average expression level of exhausted genes to measure the functional state of CD8⁺ T cells in nLung, SSN, and mLUAD. P value was calculated by two-sided unpaired Kruskal-Wallis rank sum test.

Fig. 4. Detailed characterization of myeloid cells.

(A) UMAP projection of 18,380 myeloid cells, showing the composition of 17 main subtypes. (B) UMAP projection of 6655 myeloid cells derived from SSN. (C) Canonical cell markers were used to identity myeloid cell subtypes. (D) Heatmap of marker gene expression in myeloid clusters. (E) Average proportion of each myeloid subtype among nLung, SSN, and mLUAD. (F) Percentages of each myeloid cell subtype among nLung, SSN, and mLUAD. Y axis: Average percent of samples across the three groups. Groups are shown in different colors. Each bar plot represents one cell cluster. Error bars represent ± SEM for normal and tumor samples. Colored dots represent different samples. All differences with P < 0.05 are indicated; two-sided unpaired Wilcoxon rank sum test was used for analysis. (G) Heatmap showing the markers of different DC subtypes. (H) Violin plots showing the expression of IDO1 in DC subtypes, split by sample origin. P values were calculated by differential expression test (DE test) using a pseudo-bulk method with Benjamini-Hochberg–corrected value. NA: P values cannot be calculated because of low expression. FDR, false discovery rate; NA, not applicable. (I) Violin plots showing the expression of M1, M2, and TAM markers in macrophage subtypes. (J) Violin plots showing the expression of example Macro-C7 markers involved in glycolysis and angiogenesis. (K) Kaplan-Meier plot showing that patients with LUAD in TCGA dataset with high expression of Macro-C7 cluster markers have shorter overall survival. The high and low groups were divided by the 75% quantile value of the mean expression level of the Macro-C7 gene set.

Fig. 5. Distinct EC and fibroblast subtype distribution in SSN.

(A) UMAP projection of 3381 ECs, showing the composition of six main subtypes. (B) Heatmap of marker gene expression in endothelial clusters. (C) Average proportion of each subtype between nLung, SSN, and mLUAD. (D) Percentages of each EC subtype among nLung, SSN, and mLUAD. Y axis: Average percent of samples across the three groups. Groups are shown in different colors. Each bar plot represents one subtype. Error bars represent ± SEM for normal and tumor samples. Colored dots represent different samples. All differences with P < 0.05 are indicated; two-sided unpaired Wilcoxon rank sum test was used for analysis. (E) Differentially expressed pathways are scored per cell by GSVA between six endothelial subtypes. The relative activity scores were obtained from a linear model by limma and sorted by pathway activity in Endo-C5 cells. (F) UMAP projection of 2257 fibroblasts, showing the composition of five main subtypes. (G) Heatmap of marker gene expression in fibroblast clusters. (H) Average proportion of each fibroblast subtype among nLung, SSN, and mLUAD. (I) Percentages of each fibroblast subtype in nLung, SSN, and mLUAD. Y axis: Average percent of samples across the three groups. Groups are shown in different colors. Each bar plot represents one subtype. Error bars represent ± SEM. Colored dots represent different samples. All differences with P < 0.05 are indicated; two-sided unpaired Wilcoxon rank sum test was used for analysis. (J) Violin plots showing the expression of selected marker genes of different fibroblast subtypes. (K) Differentially expressed pathways are scored per cell by GSVA between five fibroblast subtypes. The relative activity scores were obtained from a linear model by limma and sorted by pathway activity in Fibro-C3 cells.

Fig. 6. Intercellular interactions in normal lungs and lung tumors.

(A) Circos plot showing the intercellular interactions among different cell types in nLung, SSN, and mLUAD. The strings are directional and represent interactions determined on the basis of expression of a ligand by one cell type and expression of a corresponding receptor by another cell type. The thickness of each string corresponds to the amount of different interaction pairs, colored according to cell type. (B) Dot plot showing the expression level and percentage of selected genes in different cell types among nLung, SSN, and mLUAD. (C) Violin plots showing the expression of CXCL12 and CX3CL1 in different EC subtypes, split by sample sources. (D) Dot plot showing the mean expression level and percentage of selected interaction pairs involved in EMT, lymphocyte homing, and angiogenesis. The expression of each gene was considered separately for each sample source.