Dynamic Profiling of the Tumor Microenvironment in Lung Squamous Cell Carcinoma by Multidimensional Transcriptomics

Data Collection

We obtained scRNA-seq data from 53 tumor and adjacent normal tissue samples from 17 untreated male patients with LSCC through three databases: the Gene Expression Omnibus (GEO) database (GSE194070), the National Genomics Data Center (https://ngdc.cncb.ac.cn/cancerscem), and the lung cancer database of West China Hospital, Sichuan University (http://lungcancer.chenlulab.com). These datasets included tumor samples categorized by TNM stages I, II, III, and IV (5, 14, 6, and 10 tissues, respectively), along with 18 adjacent normal tissue samples. The ST data used in the analysis were obtained from BioStudies (https://www.ebi.ac.uk/biostudies/) with accession number E-MTAB-13530. Bulk RNA-seq data comes from the TCGA database (https://portal.gdc.cancer.gov/).

Single-cell RNA-seq Data Preprocessing

We processed scRNA-seq data in R (version 4.4.1) with Seurat

⁹

(version 5.1.0). Quality control removed cells with <200 detected genes or >20% mitochondrial expression, and genes present in <3 cells. Datasets were integrated and batch-corrected using Harmony

¹⁰

(version 1.2.1). After normalization, we identified 2,000 highly variable genes (HVGs), computed S and G2/M scores, and regressed cell-cycle effects (CellCycleScoring/ScaleData). PCA on scaled HVGs yielded 50 principal components (PCs) used to construct a shared nearest-neighbor (SNN) graph (FindNeighbors) and perform clustering (FindClusters). Uniform Manifold Approximation and Projection (UMAP) was used for dimensionality reduction and visualization.

Cell Type Annotation

Cell clusters were annotated based on canonical marker genes

¹¹

, ncluding eight distinct cell types: T/NK cells (CD3D), B cells (CD79A), macrophages (CD14), dendritic cells (CD1C), neutrophils (CSF3R), mast cells (TPSAB1), epithelial cells (EPCAM), and fibroblasts (COL1A1).

Identification of Cancer Cells and Cell Stemness Analysis

To identify cancer cells, we used epithelial cells from normal control samples as reference cells and employed Inference of Copy Number Variations (inferCNV)

¹²

(version 1.20.0) to infer copy number variations (CNVs) in epithelial cells from tumor samples. CopyKAT

¹³

was used to automatically identify diploid cells as normal cells. For each non-diploid cell, breakpoints of CNVs were identified, and segments were determined using a Markov Chain Monte Carlo (MCMC) method. Normal and tumor cells were distinguished based on their distinct gene expression profiles. Finally, the results from inferCNV and CopyKAT were combined to define the tumor cell populations within the epithelial cells. Furthermore, cell stemness and differentiation potential was inferred using the CytoTRACE

¹⁴

(version 0.3.3) software package.

Construction and validation of prognostic risk model

We selected genes that were highly expressed in SPP1+ macrophages in stage IV compared to stage I. Using the "glmnet" package

¹⁵

, we performed a least absolute shrinkage and selection operator (LASSO) Cox proportional hazards regression analysis on genes that overlapped with those in the TCGA-LUSC dataset. We compiled a list of genes with nonzero coefficients, and the resulting risk model was carefully constructed by linearly summing the products of the genes and their corresponding risk coefficients. Univariate Cox regression analysis was performed on these genes in conjunction with prognostic information from the TCGA-LUSC dataset to validate their prognostic association. Multivariate Cox regression analysis was performed using age, TNM stage, and risk score to demonstrate the independent predictive power of the model. Patients were stratified into low-risk or high-risk groups based on the median risk score threshold. To systematically validate the prognostic performance of the model, we calculated the area under the curve (AUC) using the "timeROC" package

¹⁶

. Survival analysis was also performed based on the Kaplan-Meier method. The predictive robustness of the model was rigorously validated using AUC calculations from other GEO datasets.

Gene Enrichment Analysis

Gene Ontology (GO) term enrichment analysis was performed using the clusterProfiler R package

¹⁷

(version 4.12.6). For each cell cluster, marker genes with a fold change > 1.5 and adjusted P value < 0.01 were annotated based on the biological process category of GO terms. GO terms with an adjusted P value < 0.05 were considered statistically significant.

Survival Analysis

To evaluate the prognostic value of SPP1 expression, patients were stratified into high- and low-expression groups based on the optimal cutoff value determined by the surv_cutpoint() function from the survminer R package

¹⁸

(version 0.4.9). This function identifies the cutpoint that best separates patients according to overall survival using maximally selected rank statistics. Patients were then categorized using the surv_categorize() function, and survival differences between groups were assessed using Kaplan–Meier survival analysis and the log-rank test. Overall survival time was measured in months.

SCENIC Analysis

Single Cell Regulatory Network Inference and Clustering (SCENIC)

¹⁹

(version 1.3.1) analysis was performed in R to infer transcription factor (TF) regulatory networks, following the standard workflow. Input was the normalized gene expression matrix from Seurat. Regulon activity was quantified using AUC scores to reflect TF activity across cells.

Pathway Activity Scoring via GSVA

To assess hallmark pathway activity across macrophage subclusters, we performed GSVA with the GSVA R package

²⁰

(version 1.52.3). Hallmark gene sets for Homo sapiens were obtained from MSigDB via msigdbr

²¹

(version 7.5.1) for Homo sapiens. Seurat’s AverageExpression() generated subcluster-level expression matrices by averaging gene expression per gene, which served as GSVA input. Enrichment scores were visualized as heatmaps using pheatmap with row-wise Z-score normalization to highlight relative activity differences.

CellChat Analysis

Cell–cell communication was analyzed with CellChat

²²

(version 1.6.1). A CellChat object was created from the merged Seurat object, using CellChatDB’s Secreted Signaling and Cell–Cell Contact references. We computed communication probabilities with computeCommunProb, inferring and aggregating intercellular networks under default settings. Interaction strength was visualized to show aggregated networks and each cell type’s outgoing signals.

Cell type deconvolution and annotation

We used the RCTD

²³

(Deconvolution from Spatial Transcriptome Data) method to deconvolute and annotate cell types from ST data. RCTD is an innovative computational framework designed to leverage spatial transcriptome data and scRNA-seq datasets. In this analysis, we used scRNA-seq data from LSCC as a reference and inferred the cell abundance of each point in the ST data.

Spatial niche and cell distance analysis

To comprehensively understand the spatial niches and cell-to-cell interactions of cells within a sample, we used the MISTy (Multiview InterCellular Spatial Modeling Framework) approach. MISTy analyzes cell-to-cell interactions through different spatial environments, enabling the identification and understanding of the formation and changes of spatial niches. In this study, we used the MISTyR

²⁴

(version 1.12.0) package to analyze spatial niches in ST data based on the spatial distances between cells.

Single-cell landscape across different pathological stages of LSCC

The study workflow is shown in Figure 1A. After quality control, we retained 169,408 cells from 53 lung tissue samples and partitioned them into eight major lineages, based on canonical marker genes expression (Figure 1B, C): T/NK cells, B cells, macrophages, dendritic cells (DCs), neutrophils, mast cells, epithelial cells, and fibroblasts. Compared to adjacent normal lung, tumors exhibited lower infiltration of T cells and NK cells, enrichment of B cells, and reduced infiltration of macrophages and DCs (Supplementary Figure 1A). Interestingly, the abundances of inferred cell populations, for instance epithelial cells, T/NK cells, macrophages etc. in the TME of LSCC also varied markedly with advancing TNM stage (Figure 1D). Collectively, these shifts in cellular composition suggest stage-dependent remodeling of the LSCC tumor microenvironment that may contribute to disease progression.

Epithelial Cell Heterogeneity During LSCC Progression

To futher dissect the epithelial heterogeneity in LSCC, we applied inferCNV to 29,782 epithelial cells using matched normal epithelial cells as the reference. As expected, both the inferCNV heatmap (Figure 2A) and CNV scoring (Supplementary Figure 2A) revealed stage-dependent CNV patterns. Then, integration of CopyKAT results (Figure 2B) with inferCNV outputs enabled robust separation of malignant from normal epithelial cells. Futher, Reclustering of malignant cells classified into seven subclusters (Figure 2C). CytoTRACE analysis demonstrated stemness heterogeneity arcoss TNM stage (Figure 2D), with T_C1 exhibiting the highest stemness potential pattern (Figure 2E). Differential expression analysis highlighted up-regulated genes of T_C1, such as IDO1, POLR2J3.1, FUT3, CLCA3P, TMPRSS2, and LRG1 (Figure 2F). Consistently, Hallmark GSVA showed T_C1 to be preferentially enriched for hypoxia, EMT, TGF-β, Hedgehog, and KRAS signaling relative to other epithelial clusters (Supplementary Figure 2B). Together, these findings strongly indicate that heterogeneity of tumor cells at different TNM stages and T_C1 represents a stem-like epithelial subpopulation within LSCC tumors niche.

Characteristics of immune cell infiltration during the progression of LSCC

We firstly analyzed the changes in T/NK cells, B cells and macrophages at different stages of LSCC. We divided T/NK cells into fourteen subclusters based on the expression of canonical marker genes, including CD4_Tn_LTB/CD28 (naive cells), CD4_Treg_TIGIT/CTLA4 (regulatory T cells), CD8_EM_GZMA (effector memory T cells), CD8_TRM_ITGAE (tissue resident memory T cells), NK_CD16hi_GNLY/NKG7/KLRF1 (NK cells with high CD16 expression), NK_CD56bright_CD7 (NK cells with bright CD56 expression), NK_CD11d_CD247 (CD11d-positive NK cells), NKT_S1PR1 (Natural Killer T cells), gdT_KLRC2 (Gamma delta T cells), and MAIT_IL7R (Mucosal-Associated Invariant T cells) (Figure 3A-B). Compared with other stages, the proportion of CD4_Tn cells in stage IV LSCC decreased significantly, while the proportion of CD4_Treg cells increased significantly (Figure 3C). CD4_Treg cells are the main immunosuppressive cell population in tumors. In line with previous study, higher proportion of CD4_Treg cells often indicates a "cold" state of TME, and the dysfunction or excessive presence of CD4_Treg cells subgroup will significantly weaken the effect of immunotherapy

²⁵

. Our results are consistent with previous studies.

Previous study showed infiltration of B cells was associate with prognosis of lung cancer patients

²⁶

. In the presented study, B cells are annotated into three subclusters, including Bn (naive B cells), Bm (memory B cells) and plasma cells (Figure 3D-E). Interestingly, the proportion of plasma cells was higher in stage I and IV tumor niche, while the proportion of memory B cells was higher in stage II and III tumor tissues (Figure 3F). The above results indicated B cell subtypes-mediated immunity pathways activation might be associated with LSCC prognosis.

We investigated the distribution patterns and cellular states of macrophage subtypes across different pathological stages of LSCC. In total, 17,889 macrophages were categorized into six subclusters, including MT2A⁺_mac, FOLR2⁺_mac, MARCO⁺_mac, CCL4⁺_mac, SPP1⁺_mac, and CXCL9⁺_mac (Figure 3G). The annotation of these macrophage subtypes was supported by the expression profiles of classical macrophage markers as well as subtype-specific signature genes (Figure 3H). We next compared the distribution of these macrophage subtypes across pathological stages I-IV of LSCC. Notably, the proportion of SPP1⁺ mac progressively increased with advancing TNM stage (Figure 3I, Table 1), suggesting a potential association between SPP1⁺ mac and LSCC progression. Previous studies have shown abundance of macrophages significantly increased in the mid‐to‐late‐stage and higher levels of macrophage infiltration are associated with poor prognosis, which supports our hypothesis

²⁷

. Taken together, these results indicated that SPP1+ mac subpopulation potentially play a central role in mediating cellular crosstalk network regulating LSCC prognosis within the immune microenvironment.

Prognostic model constructed based on SPP1+ macrophages

To test whether stage IV–specific alterations in SPP1⁺ macrophages are linked to LSCC prognosis, we constructed a polygenic risk model, named SPP1⁺Mac Differential Gene Score (SPP1⁺Mac-DGS) based on the differentially expressed genes (DEGs) between stage IV and stage I SPP1⁺ macrophages. Genes with avg_log₂FC > 0 and p < 0.05 were considered up-regulated gene cocktails (Figure 4A). Using the male subset of TCGA-LUSC as the training cohort, 516 up-regulated genes as the input; 16 absent from TCGA were excluded, leaving 500 candidates for LASSO Cox modeling (Figure 4B, C). The final model comprised four genes, including TBC1D1, GNG5, ID1, and MAX; with the formula: SPP1+Mac-DGS =0.0473×TBC1D1Exp + 0.0816×GNG5Exp + 0.0065×ID1Exp − 0.1092×MAXExp. Univariate Cox analyses supported the prognostic relevance of these genes (Figure 4D), and multivariable Cox analysis confirmed SPP1⁺Mac-DGS as an independent predictor (Figure 4E).

Patients stratified by the median score showed significantly different overall survival (OS), with longer OS in the low-risk group (log-rank p < 0.0001; Figure 4F). In the training cohort, time-dependent ROC curves yielded AUCs of 0.602, 0.683, and 0.675 for 1-, 3-, and 5-year OS, respectively (Figure 4G); the external validation cohort (GSE30219) showed comparable performance with AUCs of 0.625, 0.655, and 0.682 (Figure 4H). We referenced the AUC values of prognostic models in previous literature. Although our model's AUC value was not as high, it had similar accuracy to those in previous literature

²⁸

. Collectively, these results indicate that SPP1⁺Mac-DGS provides stable, modest prognostic discrimination in LSCC.

Functional Characterization of SPP1+ Macrophages

To delineate subtype-specific programs, we compared DEGs between SPP1⁺ macrophages and other macrophage subsets (Figure 5A). Genes enriched in SPP1⁺ macrophages, such as S100A10, TPT1, and PTGES3 have been implicated in immunosuppression and tumor invasion/metastasis, whereas genes preferentially expressed in other subsets (e.g., PTP4A3, SPDYA, HP, CCL17) encompass cell-cycle control and T-cell activating functions and immunomodulatory activities. These patterns suggest that multiple macrophage lineages may influence tumor progression via distinct mechanisms. The GO enrichment analysis further highlighted biological differences between SPP1⁺ macrophages and other subtypes (Figure 5B–C). SPP1⁺ macrophages showed up-regulation of antiviral response and metabolic reprogramming pathways accompanied by down-regulation of T-cell co-stimulation, antigen-receptor signaling, and leukocyte adhesion. This transcriptomic signature indicates a deviation from classical pro-inflammatory macrophage states, consistent with reduced immunostimulatory capacity and an immunosuppressive phenotype. Such features align with the known tissue-remodeling and tumor-supportive roles of SPP1-high cells.

In the TCGA-LUSC cohort, SPP1 expression was significantly higher in tumors than in normal lung (Figure 5D), and high SPP1 expression level was associated with worse overall survival (Figure 5E), supporting the link between SPP1⁺ macrophage programs and adverse prognosis.

Pathway activity profiling revealed subtype-specific signaling (Figure 5F). Notably, SPP1⁺ macrophages exhibited the highest activity in MYC targets, angiogenesis, oxidative phosphorylation, DNA repair, and EMT, underscoring their pro-tumor features and suggesting contributions to progression through metabolic remodeling, immune suppression, and promotion of tumor cell plasticity. Finally, SCENIC analysis identified possible regulators of SPP1⁺ macrophages, including MAX, CLAF1, HDAC2, ELF1, and ETV6 (Figure 5G).

The SPP1⁺ macrophages-fibroblasts interaction marks an edge-localized, immune-excluding niche in LSCC

To further decipher the signaling features of ligand-receptor interactions within the immune microenvironment at distinct TNM stages. CellChat revealed that patterns of incoming and outgoing interaction strengths varied across LSCC stages; notably, macrophages consistently exhibited higher outgoing strength than most other lineages (Figure 6A), indicating a dominant signal-sender role in the TME. Within the SPP1 signaling pathway, the macrophage-fibroblast axis was among the most pronounced interactions (Figure 6B). Across stages, CellChat predicted macrophage-derived SPP1 engaging integrin heterodimers on fibroblasts, including ITGA8+ITGB1 (Figure 6C). Violin plots of SPP1, ITGA8, and ITGB1 expression across cell types further demonstrated the cellular sources and targets of these interactions (Figure 6D).

To futher dissect SPP1⁺ macrophages-fibroblasts interactions spatial feature, we analyzed a stage IB LSCC spatial transcriptomic dataset. RCTD deconvolution mapped the locations and abundances of SPP1⁺ macrophages, fibroblasts, and tumor cells (Figure 7A), revealing that SPP1⁺ macrophages and fibroblasts were enriched at the tumor edge. MISTy further identified proximity-based associations (Figure 7B, C), and spatial distribution maps showed colocalization in the intra view, consistent with the CellChat-inferred communication. This organization mirrors reports in NSCLC of a marginalized, immune-exclusive TME in which SPP1⁺ TAMs and stromal-active cancer-associated fibroblasts（CAFs） (e.g., FAP⁺/POSTN⁺) assemble at the tumor boundary and are associated with restricted T-cell entry

²⁹

.

In sum, these multi-modal data support an SPP1-centered axis at the tumor–stroma interface in LSCC: macrophages act as dominant senders, signaling through fibroblast integrin receptors to establish an immune-exclusive marginal niche that may impede lymphocyte infiltration.