This study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Union Hospital, Tongji Medical College, Huazhong University of Science and Technology (No. 2019/S877). Four HS and four initially diagnosed COPD patients were enrolled, and all participants provided written informed consent. COPD patients were diagnosed according to the guidelines of the global initiative for chronic obstructive lung disease (GOLD, 2018). All COPD patients were current smokers and had not received systemic treatment, including anticholinergics and glucocorticoids, within three months prior to the study. The HS cohort consisted of current smokers with a smoking history of more than 20 pack-years and normal lung function. Current smokers were defined as individuals who had smoked within the last three months prior to bronchoscopy. Participants were excluded if they had a history of malignant tumors, unstable cardiac disease, allergic and autoimmune diseases, or other acute or chronic pulmonary diseases such as asthma, bronchiectasis, fibrosis, pneumonia, or sarcoidosis. Supplementary Table S1 summarizes the demographic and clinical characteristics of the participants.

Paired BALF and peripheral blood samples were collected from each participant. BALF was obtained from the middle lobe of the right lung of the participants via bronchoscopy at Union Hospital, Tongji Medical College, Huazhong University of Science and Technology. The BALF was diluted with phosphate buffered saline (PBS) to a final volume of 50 mL and filtered through a 100 μm nylon cell strainer (ThermoFisher Science). After centrifugation at 300 g for 10 minutes, the supernatant was removed. Then, cells were stained with CD45 antibodies, and CD45+ immune cells were isolated through flow sorting for downstream 10x scRNA-seq, scTCR-seq, and flow cytometry analyses. Peripheral blood was drawn on the same day as bronchoscopy. Peripheral blood mononuclear cells were isolated using Ficoll-Hypaque gradient centrifugation (Pharmacia, Uppsala, Sweden) and resuspended in PBS for subsequent 10x scRNA-seq, scTCR-seq, and flow cytometry analyses.

Following the manufacturer’s instructions (10x Genomics, Pleasanton, CA), Chromium Single Cell 5’ Library & Gel Bead Kit (PN-1000006) was utilized for single-cell capture and library preparation. In brief, the cell suspension, barcoded gel beads, and partitioning oil were introduced into the 10x Genomics Chromium Kit to form single-cell Gel Beads-in-Emulsion (GEMs). Captured cells were lysed, and transcripts were barcoded through reverse transcription within individual GEMs. The cDNA, along with the corresponding cell barcodes, was then amplified via PCR. Libraries for scRNA-seq were constructed using the 5’ Library Kits (PN-1000006), and libraries for scTCR-seq were prepared using the V(D)J Enrichment Kits for Human T Cells (PN-10000005). Each sample was processed independently without a hashing cell. Sequencing was performed on an Illumina NovaSeq 6000 platform.

The 10x Chromium scRNA-seq data were processed using the CellRanger toolkit (v.6.0.0) for alignment, barcode assignment, and unique molecular identifier (UMI) counting (using the GRCh38 human reference genome). The Seurat package (v.4.1.1) (17) in R (v.4.3.1) was used to analyze the filtered gene expression matrices. Quality control (QC) was applied based on three metrics: (1) the number of detected genes was between 200 and 6,000; (2) the percentage of mitochondrial genes was below 10%; and (3) doublets were identified and removed using the R package DoubletFinder (v.2.0.3) (18). After QC filtering, a total of 62,738 single cells were retained for subsequent analysis.

Data were normalized and scaled using the Seurat functions NormalizeData and ScaleData, and highly variable genes (HVGs) were identified using FindVariableGenes. Principal component analysis (PCA) was performed using the top 2,000 HVGs, and the top 20 principal components (PCs) were selected as the features in the PCA space. To enable joint analysis across samples, we applied the Harmony algorithm (19) for batch effect correction. The integration efficacy was validated by comparative visualization of cellular distributions before and after correction ( Supplementary Figures S1A, B ). Clustering was performed using the Seurat functions FindNeighbors and FindClusters, and dimensionality reduction was visualized using Uniform Manifold Approximation and Projection (UMAP) via the RunUMAP function (20). Cluster markers were identified using the FindAllMarkers function, and clusters were annotated using canonical and signature markers. Clusters expressing two or more sets of canonical markers of cell types were classified as doublets or undefined cells and excluded from further analysis.

To characterize cell types in detail, 1–3 rounds of clustering were performed on BALF and PBMC separately. The first round of clustering (resolution = 0.5) identified three cell types: T/NK cells, myeloid cells, and B cells. The second round of clustering was performed on T/NK cells and myeloid cells to further characterize subsets. Then, we performed the third round of clustering on CD4+ T and CD8+ T cells. The second and third rounds of clustering were performed using the top 20–30 PCs and the resolution ranging from 0.2 to 0.6.

We identified differentially expressed genes (DEGs) based on the Wilcoxon rank-sum test using the Seurat function FindMarkers. Unless noted otherwise, we selected the genes with P value < 0.05 and min.pct > 0.1 as significant DEGs. Pathway enrichment analysis was performed using the enricher and GSEA functions in the clusterProfiler package (v.4.2.2) (21). Gene identifiers were mapped using the R package org.Hs.eg.db, and pathway terms were obtained from the msigdbr package (v.7.5.1), including Kyoto Encyclopedia of Genes and Genomes, Gene Ontology (GO), Hallmark, and Reactome databases. Pathways with P value < 0.05 were considered significantly enriched.

Gene set variation analysis (GSVA) was performed using the GSVA package (v.1.42.0) (22) to predict the functional states of CD8+ T cell and macrophage subtypes. Gene sets were obtained from the msigdbr package (v.7.5.1), and heatmaps displayed the mean expression level of each cell cluster.

The AUCell package (v.1.16.0) was used to score individual cells for signature gene sets derived from previous literature (23–28). Cytotoxicity- and phagocytosis-related gene sets were obtained from the cytolysis and phagocytosis pathway terms in the GO Biological Process category. Details of the gene sets are provided in Supplementary Table S2 .

Spearman’s rank correlation was used to evaluate the similarity across CD8+ T cell clusters. The Spearman’s rank correlation was calculated using the mean expression values of the top 2,000 HVGs of each cell cluster.

Cell pseudotime trajectories were inferred by Monocle2 using the monocle package (v.2.26.0) (29). Marker genes for each cluster were identified using the differentialGeneTest function, and the top 2,000 genes with the lowest q-values were used for pseudotime ordering via the reduceDimension and orderCells functions. Dimensionality reduction and visualization were performed using DDRTree and plot_cell_trajectory functions. After the cell trajectories were constructed, DEGs along the pseudotime were detected using the differentialGeneTest function. Branch-dependent gene expression patterns were identified using branched expression analysis modeling (BEAM). A specialized heatmap generated by the visCluster function in the ClusterGVis package (v.0.1.1) visualized the top 2,000 branch-dependent genes. Additionally, trajectory analysis was also performed using diffusion maps via the Destiny package (v.3.8.1) (30), and connectivity between cell clusters was assessed using the partition-based graph abstraction (PAGA) algorithm in Scanpy (31).

Raw data of scTCR-seq from 16 samples were processed using Cell Ranger (v.6.0.0) against the GRCh38 human VDJ reference genome. TCR analysis was performed using the immunarch (v.0.9.0) and scRepertoire (v.1.12.0) packages. We used alpha-beta T cells to quantify unique clonotypes scaled to the total number of clonotypes recovered based on gene sequencing. The frequency of clonotypes was categorized by the number of cell counts and classified as single, small, medium, large, and hyperexpanded clonotypes. We calculated the Shannon index to measure clonal diversity and analyzed the clonal homeostasis and clonal overlap of T cells.

CellChat v2 package (32) was used to explore cell-cell communication networks via ligand-receptor interactions. The CellChatDB was set as “Secreted Signaling”, and separate cellchat objects were generated for HS and COPD groups. Interaction comparisons were performed between COPD and HS groups using the compareInteractions and RankNet functions, and the upregulated/downregulated signaling pathways were identified.

The expression of surface markers and intracellular molecules of cells was determined using flow cytometry. Cells were stained with fluorochrome-conjugated antibodies, which were purchased from BD Biosciences or Biolegend. Immune cells were surface-stained with fluorochrome-conjugated antibodies. The samples were incubated with antibodies for 15 min at 4°C. Cells were resuspended in PBS and washed at 400 g for 6 min. After fixation and permeabilization (eBioscience), intracellular proteins were labeled with the corresponding mAbs conjugated with fluorescent molecules, according to the manufacturer’s instructions. Flow cytometry was performed on a BD LSRFortessa X-20 and analyzed with FlowJo V10 software.

All statistical analyses were implemented via R (v.4.3.1) or Graphpad prism 9. For continuous variables, t-tests or Wilcoxon tests were used to compare differences between two groups, while Kruskal-Wallis tests were performed for multiple groups. Correlations between variables were estimated with Pearson or Spearman correlation analysis. Statistical significance was defined as P value < 0.05.

To explore the immunological characteristics of COPD, we performed single-cell transcriptomic analysis of immune cells in BALF and PBMC from HS and COPD patients ( Figure 1A ). Following QC and filtering, we obtained raw data for 62,738 cells for downstream analysis, including 18,697 cells from BALF and 44,041 cells from PBMC ( Supplementary Table S3 ). Clustering analysis was performed on BALF and PBMC samples, respectively. Firstly, initial clustering characterized three same major cellular compartments in BALF and PBMC, encompassing T/NK cells (CD3D, CD4, CD8A, and NKG7), myeloid cells (CD14, CD68, FCER1A, and CD1C), and B cells (JCHAIN, MS4A1, CD79A, and CD19) ( Figure 1B , Supplementary Figures S1C, D , Supplementary Table S4 ). Then, T/NK cells showed the identical cell types between BALF and PBMC, including CD4+ T, CD8+ T, γδ T, MAIT, NKT, NK, and proliferating T (proli.T) cells ( Figure 1B ). However, myeloid cells had distinct cell types in different sample types, with macrophages and dendritic cells (DCs) in BALF but monocytes, DCs, and megakaryocytes in PBMC ( Figure 1B ), demonstrating the accuracy of our data analysis. B/plasma cells only constituted a minor fraction of BALF and PBMC ( Figure 1C ). Importantly, all cell types were derived from multiple patients ( Figure 1C ), confirming that cells were clustered according to immune characteristics rather than patient specificity.

Next, we compared the proportions and features of T/NK cell types between different anatomic compartments. Among T/NK cells, CD8+ T and MAIT cells showed higher percentages in BALF, while γδ T and NKT cells exhibited higher proportions in PBMC ( Figure 1D ). Flow cytometry data demonstrated a lower CD4/CD8 ratio in BALF than PBMC ( Supplementary Figures S2A, B ), indicating the predominance of CD8+ T cells in the airway. Further analysis of CD62L and CD45RA expression by flow cytometry revealed distinct T cell subset distributions ( Supplementary Figure S2C ). Compared to PBMC, CD4+ T cells showed reduced naive T (Tn) cell frequencies but increased effector memory T (Tem) cell proportions in BALF, and CD8+ T cells displayed lower proportions of central memory T (Tcm) and Temra cells but higher percentages of Tem cells in BALF ( Supplementary Figure S2C ). Then, we performed gene set score analysis and DEG analysis for the integrated T/NK cells from BALF and PBMC. Both CD4+ T and CD8+ T cells demonstrated lower naïve scores and reduced expression of Tn cell markers (e.g., SELL, CCR7, and LEF1) in BALF compared to PBMC ( Figures 1E, F , Supplementary Table S5 ). Moreover, they upregulated the expression of tissue-resident markers (e.g., CD69, ITGAE, and CXCR3) in BALF, which was corroborated by flow cytometry demonstrating elevated expression of CD69, CD103, CXCR3, and CCR6 in CD4+ T and CD8+ T cells from BALF ( Figure 1F , Supplementary Figure S2D ). These findings highlight their identity as tissue-resident memory T cells (Trms) in the airway compartment. Importantly, CD4+ T and CD8+ T cells showed higher exhaustion signatures in BALF than PBMC, which was confirmed by elevated expression of PD-1 in BALF ( Figure 1E , Supplementary Figure S2E ). CD8+ T, γδ T, NKT, and NK cells were identified as the predominant cytotoxic lymphocyte populations ( Figures 1G, H ). Notably, these cells exhibited lower cytotoxicity scores and reduced expression of cytotoxic molecules (GZMB and perforin) in BALF compared to PBMC ( Figure 1I , Supplementary Figure S3 ). Pathway enrichment analysis revealed that T/NK cells in BALF upregulated the genes associated with metabolic pathways (hypoxia, glycolysis, and fatty acid metabolism) and stress responses (cellular response to starvation, response to virus, IFN-α/β signaling, IFN-γ response, and response to endoplasmic reticulum stress) compared to their PBMC counterparts ( Figure 1J ), indicating the different immune microenvironment between the airway and peripheral blood.

In conclusion, our analysis reveals compartment-specific immune status in the airway, especially Trms with heightened exhaustion and reduced cytotoxicity, likely driven by local microenvironmental cues.

Compared with HS, PBMC CD8+ T cells showed higher cytotoxicity scores in COPD, whereas BALF CD8+ T cells exhibited lower cytotoxicity scores in COPD ( Figure 2A ). CD8+ T cells from BALF had a lower proportion and higher exhaustion scores in COPD than HS, while no significant differences were observed in PBMC ( Figures 2B, C ). Thus, we further focused on CD8+ T cells in the airway. As expected, BALF CD8+ T cells upregulated PDCD1 expression but downregulated expression of effector function-associated genes (e.g., GZMB, IFNG, and TNF) in COPD compared to HS ( Supplementary Figure S4A ). Flow cytometry analysis further corroborated the elevated PD-1 expression and concomitant reduction in GZMB and perforin levels in CD8+ T cells from COPD airways ( Supplementary Figures S4B, C ). Furthermore, compared to HS, CD8+ T cells from COPD airways downregulated the functional pathways, including cell-cell adhesion, T cell migration, defense response, IL-2 production, and TCR signaling; conversely, they upregulated the T cell exhaustion-associated pathways, such as influenza infection, PD-1 signaling, IFN-γ signaling, and response of EIF2AK4 (GCN2) to amino acid deficiency ( Figure 2D ).

On deep phenotyping, CD8+ T cells in PBMC were divided into a Tn cluster (CD8Tn), a Tcm cluster (CD8Tcm), two Tem clusters (CD8Tem_early and CD8Tem), and three effector T cell (Teff) clusters (CD8Teff_c1, CD8Teff_c2, and CD8Teff_ZNF683) ( Figure 2E , Supplementary Figure S5A ). In BALF, CD8+ T cells were subdivided into six clusters ( Figure 2F ). Concretely, CD8T_IL7R and CD8T_KLRG1 clusters were Tcm (IL7R, LEF1, and CCR7) and Tem/Temra (KLRG1, GZMK, EOMES, CST7, and S1PR5) cells, respectively ( Supplementary Figure S5B ). CD8Tex_PDCD1 cells highly expressed exhaustion markers (PDCD1, HAVCR2, TIGIT, and CTLA4) and late activation markers (e.g., HLA-DRA, HLA-DRB1, and HLA-DRB5) ( Supplementary Figure S5B ), consistent with previous studies showing the co-expression of HLA-DR alongside the characteristic exhaustion phenotype of CD8+ T cells (33, 34). In comparison, CD8Trm_IFNG, CD8Trm_FOS, and CD8Trm_LAG3 clusters exhibited high expression of tissue-resident markers (CD69 and ITGAE) and high activities of Trm differentiation-related pathways, including TGF-β signaling, IL-15 signaling, and integrin signaling ( Supplementary Figures S5B, C ). CD8Trm_IFNG cluster highly expressed cytokines and chemokines (e.g., IFNG, TNF, XCL1, and XCL2), while CD8Trm_FOS cluster showed high expression of genes associated with TCR signaling (e.g., FOS, JUN, and NR4A1) ( Supplementary Figure S5B ). CD8Trm_LAG3 cells were characterized by high expression of exhaustion markers (LAG3, PDCD1, HAVCR2, and CTLA4) and IFN-stimulated genes (ISGs; e.g., IFI44, ISG15, and ISG20) ( Supplementary Figure S5B ), reflecting a terminal exhaustion state driven by chronic IFN-I stimulation (35). By examining signature gene sets, we observed distinct functional status for each CD8+ T cell subset, with the highest cytotoxicity scores of CD8Teff, the highest inflammation scores of CD8T_KLRG1, the highest exhaustion scores of CD8Tex_PDCD1, and the highest dysfunction and immunosuppression scores of CD8Trm_LAG3 ( Figure 2G ). These findings reveal two different exhausted CD8+ T (CD8+ Tex) cell subsets (CD8Tex_PDCD1 and CD8Trm_LAG3) in the airway.

Notably, CD8T_KLRG1 and CD8Tex_PDCD1 clusters were preferentially enriched in COPD airways compared to HS ( Figure 2H ). Both CD8Tex_PDCD1 and CD8Trm_LAG3 clusters showed higher exhaustion and inflammation but lower cytotoxicity scores, while CD8Teff clusters exhibited higher cytotoxicity scores in COPD ( Supplementary Figure S5D ). We further analyzed the TCR clonotypes of all CD8+ T cell subsets from BALF and PBMC. The distinct patterns of clonal expansion indicated heterogeneous TCR repertoires between HS and COPD ( Figures 2I, J ). Compared to HS, COPD patients showed reduced TCR diversity of CD8+ T cells, as evidenced by the decreasing trend of Shannon indices in both BALF and PBMC ( Figure 2K ). Hyperexpanded TCR clonotypes accounted for 38.8% of CD8+ T cells in COPD blood compared to 22.1% in HS blood ( Figure 2L , Supplementary Table S6 ). 6.6% of CD8+ T cells were hyperexpanded in COPD airways, whereas no hyperexpanded TCR clonotypes were detected in HS airways ( Figure 2L , Supplementary Table S6 ). CD8Teff cells exhibited predominant hyperexpanded clonotypes, with greater clonal expansion in COPD compared to HS ( Figure 2L , Supplementary Table S6 ). Within COPD airways, hyperexpanded TCR clonotypes were predominantly enriched in CD8Tex_PDCD1 cells ( Figure 2L , Supplementary Table S6 ), suggesting that the hyperexpansion of a minority of clonotypes may contribute to their increased overall abundance.

Together, these data indicate increased exhaustion, impaired cytotoxicity, and reduced TCR diversity of CD8+ T cells in COPD airways.

Next, we combined gene expression and TCR data to construct the potential developmental trajectories for all CD8+ T cells. On a global scale, CD8+ T cells could differentiate from Tn cells to either Teff (Path 1, effector trajectory) or Tex cells (Path 2, exhaustion trajectory). Monocle2 analysis showed that CD8Tn cells were at the beginning of the trajectory, whereas CD8Teff_ZNF683 cells and CD8Trm_LAG3 cells were at the terminal state of Path 1 and Path 2, respectively ( Figures 3A, B , Supplementary Table S7 ). Such developmental trajectories for Teff and Tex cells were further supported by diffusion map analysis ( Figure 3C ). Cytotoxicity scores were gradually increased along the pseudotime in Path 1, while exhaustion scores were progressively increased in Path 2 ( Figure 3D ), further confirming the differentiation process of Teff cells and Tex cells, respectively. Next, PAGA analysis showed that CD8Tcm and CD8Tem_early clusters in PBMC exhibited high connectivity with CD8T_IL7R and CD8T_KLRG1 clusters in BALF, respectively ( Figure 3E ). Meanwhile, correlation analysis revealed the high similarity between CD8Tcm and CD8T_IL7R as well as the high similarity between CD8Tem_early and CD8T_KLRG1 ( Figure 3F ). Combined with the pseudotime data ( Figure 3B ), these results suggest that peripheral blood CD8Tcm and CD8Tem_early cells may migrate into the airway and subsequently differentiate into CD8T_IL7R and CD8T_KLRG1 cells. Then, TCR analysis showed that CD8T_KLRG1 cells exhibited a high degree of clonal overlap with CD8Tex_PDCD1 cells ( Figure 3G ). Combined with the pseudotime data, these findings suggest that CD8Tex_PDCD1 cells may develop from CD8T_KLRG1 cells; similarly, CD8T_IL7R cells may differentiate into CD8Trms and ultimately terminate in CD8Trm_LAG3 cells ( Figures 3B, G ). Hence, there are two potential exhaustion paths for airway CD8+ T cells: one via Tem cells to CD8Tex_PDCD1 cells and the other via Trms to CD8Trm_LAG3 cells.

Subsequently, we investigated the transcriptional changes associated with the trajectory branch point using BEAM analysis and observed that branch-dependent genes were categorized into five clusters ( Figure 3H ). Cluster 1 genes (e.g., CX3CR1, PRF1, and GZMB) showed progressive upregulation along Path 1 pseudotime, related to cytotoxic T cell differentiation, whereas Cluster 2 genes (e.g., PDCD1, LAG3, CTLA4, and HAVCR2) increased along Path 2, associated with T cell exhaustion ( Figure 3H ). In contrast to the enrichment of T cell activation, proliferation, cytotoxicity, and migration pathways in Cluster 1 genes, Cluster 2 genes were involved in response to virus, oxidative stress, amino acid deficiency, and cholesterol homeostasis ( Figure 3I ), which can trigger endoplasmic reticulum (ER) stress and unfolded protein response (UPR) (36, 37). Furthermore, both CD8Tex_PDCD1 and CD8Trm_LAG3 cells in COPD airways upregulated the pathways associated with protein synthesis and processing, such as infectious disease, ribosome biogenesis, translation, and SRP-dependent cotranslational protein targeting to membrane ( Figure 3J ). This may increase the protein-folding burden on the ER, potentially leading to ER stress. As expected, the UPR-related pathway and genes (ERN1 and XBP1) were upregulated in CD8+ Tex cells from COPD compared to HS ( Figures 3K, L ).

Taken together, these results identify two distinct CD8Tex cell clusters originating from different developmental paths; ER stress may be a potential driver of CD8+ T cell exhaustion and dysfunction in COPD airways.

We identified four clusters of CD4+ T cells in BALF ( Figure 4A ). CD4T_LEF1 and CD4T_GZMK clusters were associated with central memory (IL7R, TCF7, CCR7, and LEF1) and effector memory (GZMK, GZMA, EOMES, and NKG7), respectively ( Figure 4B ). CD4T_FOXP3 cluster exhibited high expression of FOXP3, IL2RA, CTLA4, and TIGIT, known markers of regulatory T cells (Tregs), and displayed significant anti-inflammatory and regulatory signatures ( Figures 4B, C ). CD4T_XCL1 cluster represented CD4+ Trms, strongly expressing CD69, ITGAE, and CXCR6 ( Figure 4B ). In PBMC, CD4+ T cells were classified into a Tn cluster (CD4T_CCR7), a memory T (Tm) cluster (CD4T_GPR183), a Temra cluster (CD4T_GZMK), a Treg cluster (CD4T_FOXP3), and an ISG+ T cluster (CD4T_IFI44L) ( Figures 4D, E ).

Compared to HS, Tregs (CD4T_FOXP3 cells) exhibited a reduced proportion in COPD airways, with only a decreasing trend in COPD peripheral blood ( Figure 4F ). Next, we focused on the airway Tregs. They showed lower regulatory and anti-inflammatory scores but higher pro-inflammatory scores in COPD than in HS ( Figure 4G ). Besides, they displayed elevated expression of IFN-responsive genes (e.g., IFI44L, IFI6, IFIT5, MX1, ISG15, and OAS1) but decreased expression of several genes critical for Treg development, survival, and function, such as JUN, JUNB, CTLA4, CD69, NR4A2, KLF2, IL2RA, and TGFB1 ( Figure 4H ). Additionally, airway Tregs upregulated influenza infection, IFN-α/β signaling, and IFN-γ signaling pathways, but downregulated IL-2 production, TGF-β, TNF-α, Toll-like receptor (TLR), TCR, and MAPK signaling pathways in COPD ( Figure 4I ).

CD4+ Trms are a heterogeneous population, including T-helper (Th) 1 Trms, Th2 Trms, and T follicular helper (Tfh)-like resident helper cells (38). In this study, the CD4+ Trm (CD4T_XCL1) cluster was polyfunctional, co-expressing cytotoxic molecules (GZMA), chemokines (CCL4, CCL5, XCL1, and XCL2), and Th1/Th2-related cytokines (IL2, IFNG, TNF, IL4, and IL13) ( Figure 4B ). Compared to HS, CD4+ Trms exhibited elevated Th2 signatures and genes (IL4, IL5, and IL13) in COPD ( Figures 4J, K ). Polyfunctional CD4+ Trms are also indispensable in protecting the lungs against respiratory pathogens. For instance, they constitutively express high transcript levels of cytotoxic mediators, such as IFN-γ (39). However, CD4+ Trms showed lower Th1 signatures and expression of IFNG, TNF, GZMB, and PRF1 in COPD than HS ( Figures 4J, K ).

Collectively, these results reveal the impaired anti-inflammatory function of Tregs and imbalanced Th1/Th2 response of CD4+ Trms in COPD airways.

Myeloid cells were divided into four macrophage clusters and four DC clusters in BALF ( Figure 5A ). Concretely, Macro_FCGR3A cells strongly expressed MARCO, MRC1, and PPARG ( Figure 5B ), representing tissue-resident AMs (TRAMs) (40). They exhibited high activities of fatty acid beta oxidation, oxidative phosphorylation, and glycolysis pathways ( Figure 5C ), which were essential for maintaining long-lived, self-replenishing tissue-resident macrophages (41, 42). As the predominant macrophage subset, Macro_SPP1 cluster highly expressed lysosomal genes (e.g., CD63, NPC2, CTSB, and LGMN) and lipid metabolism-related genes (e.g., FABP5, SPP1, GPNMB, ABCA1, and CD9), showing the high activities of phospholipid transport, cholesterol homeostasis, collagen metabolism, and extracellular matrix (ECM) disassembly pathways ( Supplementary Figure S6A , Figures 5C, D ). Macro_SPP1 cluster deviated from the classical M1/M2 paradigm, as it simultaneously expressed both pro- and anti-inflammatory genes (e.g., TNF, CCL2, IL10, IL1RN, and LGALS3) and concurrently exhibited high pro- and anti-inflammatory scores ( Figures 5D, E ). Interestingly, SPP1+ macrophages have been described in numerous diseases, such as fatty liver (43), COVID-19 (44), dystrophic muscle (45), and cancers (46), indicating a conserved functional phenotype across diseases. Macro_FCN1 cluster highly expressed CD14, VCAN, FCN1, and S100A12, indicating monocyte-like macrophages with high responses to chronic inflammation and pathogens ( Figures 5B, C ). Macro_CCL18 cluster, marked by CCL18, FABP4, and C1QB, was implicated in acetylcholine metabolism and catecholamine secretion pathways ( Figures 5B, C ). In addition to the three traditional types, cDC1 (DC_CLEC9A), cDC2 (DC_CLEC10A), and pDC (DC_LILRA4), we also identified a non-classical DC type (DC_LAMP3) in BALF, characterized by high expression of CCR7 and LAMP3 ( Figures 5A, B ). In PBMC, we identified a classical monocyte cluster (Mono_CD14), a non-classical monocyte cluster (Mono_FCGR3A), a megakaryocyte cluster, and three traditional DC clusters mentioned above ( Figure 5F , Supplementary Figure S6B ).

Next, we integrated all macrophages and monocytes to investigate their developmental trajectories ( Supplementary Figure S6C ). Firstly, Monocle2 analysis revealed two distinct trajectory branches, both starting from Mono_CD14 cells and terminating at Mono_FCGR3A cells (Cell Fate 1) and Macro_SPP1 cells (Cell Fate 2), respectively ( Figure 5G , Supplementary Figure S6D ; Supplementary Table S8 ). Diffusion map analysis corroborated these trajectory patterns ( Supplementary Figure S6E ). Furthermore, we assessed the likelihood of connections among all subsets using the connectivity matrix of the PAGA network. Macro_SPP1 cells showed the strongest PAGA connectivity with Macro_FCN1 cells, and Macro_FCN1 cells exhibited higher connectivity with Mono_CD14 than Mono_FCGR3A cells ( Figures 5H, I ). These findings suggest that Macro_SPP1 cells are monocyte-derived AMs (MoAMs) originating from classical monocytes. Subsequently, the branch-dependent genes identified by BEAM analysis were divided into five clusters ( Figure 5J ). Notably, Cluster 2 genes (e.g., MARCO, SPP1, C1QC, TREM2, GPNMB, FABP5, APOE, IL1B, IL10, CXCL8, and CCL2) showed a gradual increase during the development of Macro_SPP1 cells ( Figure 5J ). Pathway enrichment analysis revealed that Cluster 2 genes were associated with hypoxia, response to reactive oxygen species (ROS), low-density lipoprotein (LDL) clearance, lysosome, fatty acid metabolism, and lipid storage pathways ( Figure 5K ). This aligns with a previous study showing that increased ROS levels promote the formation of lipid-laden macrophages in the lungs (47).

Interestingly, Macro_SPP1 cells were preferentially enriched in COPD, while Macro_FCGR3A cells were in HS ( Figure 6A ). Among all macrophage subsets, Macro_SPP1 cells exhibited the largest numbers of DEGs between COPD and HS ( Figure 6B , Supplementary Table S9 ), highlighting their high plasticity and immunoreactivity. Importantly, Macro_SPP1 cells in COPD airways might undergo lipid metabolic reprogramming. They upregulated genes involved in cholesterol and lipid transport and metabolism (e.g., PPARG, CD36, FABP5, LIPA, SREBF1, SOAT1, ABCA1, APOC1, NR1H3, and APOE) and pathways of long-chain fatty acid transport, VLDL particle clearance, lysosome, reverse cholesterol transport, and regulation of cholesterol esterification ( Figures 6C, D ). Moreover, the mTORC1 and PPAR signaling pathways, which modulated lipid metabolism and promoted anti-inflammatory phenotypes, were upregulated, while pro-inflammatory pathways such as Notch, NF-κB, and MAPK signaling were downregulated in COPD ( Figure 6E ). A recent study has demonstrated the induction of the monocyte-derived macrophage population (marked by Spp1, Gpnmb, Fabp5, Cd9, and Arg1) in lungs by Notch2 blockade (48). Functionally, Macro_SPP1 cells in COPD exhibited higher anti-inflammatory scores but lower phagocytosis scores, with elevated expression of immunosuppressive genes (e.g., IL1RN, CCL18, CD274, IL4I1, and IDO1) but decreased expression of efferocytosis-related genes (e.g., LRP1, PECAM1, CD44, SIGLEC1, ICAM1, and FCGR3A) ( Figures 6C, F ). Additionally, Macro_SPP1 cells upregulated proteolysis, collagen degradation, and degradation of the ECM pathways in COPD ( Figure 6D ). Consistently, they showed increased expression of protease genes (e.g., ADAM9, CTSB, CTSL, CTSS, and MMP9) but decreased expression of the antiprotease gene (CST3) in COPD ( Figure 6C ).

In summary, these results reveal that Macro_SPP1 cells undergo lipid metabolic reprogramming, exhibiting an anti-inflammatory phenotype, reduced phagocytosis, and protease-antiprotease imbalance in COPD airways.

We performed CellChat analysis to explore potential intercellular communication patterns among airway immune cells ( Figures 7A, B ). Overall, we detected 25 significant secreted signaling pathways mediating interactions across all immune cell types in the airway ( Figure 7C ). Among these, six pathways (SPP1, PLAU, CCL, ANNEXIN, IL16, and GALECTIN) were upregulated in COPD, while 13 pathways (LT, PARs, IFN-II, FASLG, TNF, GAS, RESISTIN, BAFF, TGFb, APRIL, LIGHT, GRN, and VISFATIN) were upregulated in HS ( Figure 7C ). Although the majority of ligand-receptor (L-R) pairs mediating cell-cell interactions were shared between COPD and HS, distinct L-R pairs were also identified in each group, such as SPP1-CD44, SPP1-(ITGA5+ITGB1), and SPP1-(ITGA4+ITGB1) in COPD while LTA-TNFRSF1B, LTA-TNFRSF1A, LTA-TNFRSF14, and LTA-(LTB+LTBR) in HS ( Supplementary Figure S7 ).

Notably, macrophages harbored the largest numbers of ligands and receptors, actively interacting with all cell types in both HS and COPD ( Figure 7D ). In total, macrophages in COPD received and sent more signals than HS ( Figure 7E ). Among the signals received by macrophages, we observed that MIF-(CD74+CD44), CCL3-CCR1, CCL3-CCR5, CCL4-CCR5, CCL5-CCR1, CCL5-CCR5, ANXA1-FPR1, and ANXA1-FPR3 were upregulated in COPD compared with HS ( Figure 7F ). These L-R pairs were involved in MIF, CCL, and ANNEXIN signaling pathways. The transcriptomic analysis further revealed that macrophages, CD4+ T, CD8+ T, γδ T, NK, and proli.T cells exhibited higher expression of the ligand gene MIF, and macrophages showed elevated expression of receptor genes CCR1, CCR5, FPR1, and FPR3 in COPD compared to HS ( Figure 7G ). More importantly, among the signals sent by macrophages, SPP1-CD44, SPP1-(ITGA4+ITGB1), LGALS9-P4HB, LGALS9-CD45, and LGALS9-CD44 were upregulated in COPD compared to HS ( Figure 7H ). These L-R pairs were involved in immunosuppressive SPP1 and GALECTIN signaling pathways, by which macrophages interacted with T cells. Notably, the upregulation of these L-R pairs might be attributed to the significantly increased expression of ligand genes (SPP1 and LGALS9) in macrophages from COPD ( Figure 7G ). Especially, SPP1 and LGALS9 were mainly upregulated in the Macro_SPP1 subset from COPD compared with HS ( Figure 7I ).

Thus, these findings suggest macrophages (particularly Macro_SPP1 cells) as key regulators of T cell dysfunction and identify SPP1 and GALECTIN signaling pathways as potential therapeutic targets for restoring immune homeostasis in COPD airways.