Single-cell characterization of the immune heterogeneity of pulmonary hypertension identifies novel targets for immunotherapy

Abstract

Background

Pulmonary arterial hypertension (PAH) is a critical cardiopulmonary vascular disorder marked by the progressive elevation of pulmonary artery pressure, increased pulmonary vascular resistance, and eventual right heart failure. Research has shown that various immune cells play a significant role in the pathogenesis of PAH, both in patients diagnosed with the condition and in experimental models of PAH. Cell–cell communication is important for PAH progression and therapies, while the global cell landscape of intercellular signaling has not been elucidated.

Methods

We performed single-cell RNA sequencing on NCBI Gene Expression Omnibus (GEO) databases GSE169471, GSE 210248, GSE228643 and GSE244781, and analyzed lung tissue samples across healthy controls and PAH patients. In total, approximately 124,561 cells were analyzed and a total 34 clusters were identified. We integrated the sequencing results of multiple samples and used an enhanced single-cell sequencing workflow to overcome the limitations of a single study.

Results

In this study, we elucidated the functional characteristics and potential regulatory interactions of several cell subpopulations that have not been previously documented in similar research. We constructed a comprehensive landscape of cell communications at the single-cell resolution, which is expected to significantly advance the development of personalized diagnostic and therapeutic strategies for PAH. We demonstrated the transcriptomic features of different cell types in PAH patients. We presented an in-depth analysis of T cell subsets, myeloid cell heterogeneity and a comprehensive analysis of SMCs and FBs subsets in PAH. T cell heterogeneity and functional dynamics were exhibited in PAH, which suggests that targeting cytotoxic regulation may be a potential therapeutic strategy. Significant changes and potential functions of myeloid cell subsets in PAH patients and we especially focused on GPNMB⁺ macrophages. In addition, CellChat and NicheNet analyses reveal altered intercellular communication and dys-regulated signaling pathways in PAH progression. The enhanced MIF and IL-1 signaling suggests that the induced inflammatory response in PAH is greatly driven.

Conclusions

We systematically explored the immune heterogeneity and population and target cells in PAH, which may be valuable for developing new and precise therapies.

Introduction

Pulmonary arterial hypertension (PAH) is a life-threatening pulmonary vascular disease characterized by elevated mean pulmonary pressures (mPAP > 20 mmHg), ultimately leading to right heart failure [1]. Significant advancements in PAH management have been made over the past 25 years; however, existing therapies primarily provide symptomatic relief by addressing the imbalance of vasoactive factors, and the disease continues to be life-limiting [2]. The pathogenesis of PAH is complex, involving intricate interactions among various cell types, including smooth muscle cells, endothelial cells, fibroblasts, epithelial cells, and immune cells, which collectively contribute to its heterogeneity. The complex pathobiology of PAH highlights the need for a comprehensive analysis of intercellular communication and its implications for treatment strategies [3, 4].

Single-cell technologies have proven ideal for identifying specific pathogenic cell populations with hundreds of antigens, cytokines, chemokines and transcription factors, which is key for new therapy [5]. The previous microarray analysis mainly focused on whole-genome analysis. Although there have been single-cell sequencing studies, they are limited by small sample sizes from a single center [6, 7]. Here, in our current study, it is novel to utilize a variety of state-of-the-art analytical tools and integrate multiple scRNA-seq datasets to identify the differences of functional variations, immune dynamics, transcriptional regulators, and communication atlas between pulmonary tissues from health control and PAH patients.

Hypoxia, inflammation, immune reactions, and epigenetic modifications all exert significant contributory roles in the development of PAH [8]. Research has shown that various immune cells are implicated in the pathogenesis of PAH in both patients and experimental models. Notably, macrophages, which are the primary inflammatory cells infiltrating PAH lesions, are instrumental in exacerbating pulmonary vascular remodeling associated with PAH. Macrophages typically polarize into either the classical M1 or alternative M2 phenotypes, and they contribute to the progression of PAH through the secretion of diverse chemokines and growth factors, such as CX3CR1 and PDGF [9]. Recently, it has been reported that GPNMB is a type I transmembrane glyco protein with an extracellular segment that can also be cleaved, forming a soluble GPNMB (sGPNMB) [10]. Extracellular RGD domain of GPNMB can bind to receptors (eg, integrin) and subsequently activate downstream signaling pathways. GPNMB has been found to be expressed in multiple tissues, including macrophages anvasc, and involved in multiple functions (e.g., pulmonary fibrosis and inflammation) [11,12,13]. GPNMB was also found to be highly expressed in macrophages and exhibited an anti-inflammatory effect [12]. Interestingly, the functional characterization of GPNMB has been recently reported in vascular [14, 15]. However, more studies are warranted to investigate GPNMB⁺ macrophages in PAH progression.

Through single-cell analysis, we aimed to identify the differences between health control and PAH patients to clarify their distinct pathogeneses and propose directions for prevention and treatment. Pro-inflammatory and dysfunctional GPNMB⁺ macrophages were associated with PAH, which could be expected to be potential therapeutic targets for precision medicine.

Methods

In this study, a comprehensive scRNA-seq analysis was employed to investigate the molecular and cellular mechanisms of PAH. We performed clustering and trajectory analyses to reveal cell-type composition, dynamic changes, and developmental trajectories. Enrichment and landscape analyses were performed to identify key biological pathways, TF activity, and regulatory networks driving PAH cell behavior. Cell–cell interaction analysis using CellChat elucidated cell-to-cell communication and major signaling pathways. Finally, we used scMLnet and multiichenet to construct multilayer signaling networks and infer ligand-receptor-target interactions, providing insights into the regulatory dynamics between cells and within cells. Each analytical step aims to address specific questions about cellular heterogeneity, communication, and disease-specific molecular mechanisms.

Single-cell RNA-seq data acquisition and preprocessing

We utilized scRNA-seq data to investigate cellular heterogeneity and intercellular communication in PAH. This approach enabled us to identify cell clusters, differential gene expression, and transcriptional changes between the control and PAH groups. Publicly available single-cell RNA-seq dataset GSE254617, GSE210248, GSE233189, GSE228643 were downloaded from GEO. The “read10X” function in edgeR package was used to read the single-cell RNA-seq data into R, and “CreateSeuratObject” function in Seurat (v4.3.0) was employed to create Seurat objects for preprocessing and normalization [16]. Quality control and downstream analysis were applied to the single-cell RNA-seq data using Seurat. Reserve cells with unique molecular identifiers (UMI) and gene count (ranging from 200 to 10,000 and > 500), and mitochondrial genes (< 20%). The functions “NormalizeData” and “FindVariableFeatures” were used to normalize the counts matrix of the Seurat object for each sample and obtain the top 2000 highly mutated genes, respectively. The “merge” function integrated individual samples, which was then normalized using the “ScaleData” based on the 2000 highly variable genes. Principal component analysis (PCA) conducted in 2000 highly variable genes using the “RunPCA” function. Harmony, a method bases on PCA, was employed to embed the transcriptome expression profile of the samples into a low-dimensional space, effectively removing batch effects through iterative processing. The first 10 principal components, along with resolution 0.5, were utilized in the “FindClusters” function to identify cell clusters. The identified clusters were visualized using the t-SNE method. Differentially expressed genes across clusters were calculated and identified with the “FindAllMarkers” function. Canonical cell markers for major cell types were obtained from the CellMarker 2.0 [17]. The steps were iterated for each major cell type to further identify sub-clusters and annotate them as cell subtypes in each major cell type.

Trajectory analysis

We utilized Monocle2 (v2.22.0) to investigate the single-cell trajectories and cell state transitions of different cell subtypes [18]. The matrix of UMI counts was taken as input to Monocle and created an object with “newCellDataSet” function. To construct the trajectory, obtain the high variable genes for each Monocle object using the “dispersionTable” function, filter with the mean_expression > = 0.1&dispersion_empirical > = 1*dispersion_fit as the reference. All high variable genes of cell subtypes were selected to order cells using the “orderCells” function. Dimensionality reduction was performed using the “DDRTree” function, and the developmental trajectories were visualized using the “plot_cell_trajectory” under various situations.

Enrichment analysis

Enrichment analysis was conducted to uncover the biological functions and pathways of cell subtypes, providing insights into the mechanisms driving cellular behaviors in PAH. To explore the potential functions and mechanisms of cell subtypes, enrichment analysis was performed with KEGG pathway and GO term set using the ‘clusterProfiler’ (v4.7.1.3) in R (v4.2) [19]. The “compareCluster” function was employed to compare the enrichment results across different clusters.

Single-cell regulatory network inference and clustering (SCENIC) analysis

SCENIC analysis was conducted to infer regulatory networks and identify active transcription factors across major cell types. This step was crucial for understanding transcriptional regulation and identifying key drivers of cellular behavior in PAH. SCENIC analysis was performed on quality-controlled cells of each major cell type using pySCENIC [20, 21]. The counts matrix files for each cell type were extracted and exported, and then converted into a loom file format for SCENIC analysis. “GENIE3” was used to detect the relationship between transcription factors and candidate target genes. The regulatory network between transcription factors (TFs) and genes was constructed with a search region limited to a 10 k distance from the transcription start sites (TSS) or 500 bp upstream of the TSS; “RcisTarget” was employed to validate the regulatory network by performing motif-TF enrichment analysis, selecting upstreram regulators, and identifying highly enriched motif modules while excluding indirect targets lacking motif support. The AUCell algorithm was used to score the regulon activity for each cell, producing a binary regulon activity matrix that determines which cells active regulons is.

Cell–cell interaction analysis

To understand the roles of dominant signaling regulators in PAH, cell–cell interaction analysis was conducted to explore communication between cell subtypes and identify key signaling pathways. The CellChat package (v1.6.0) in R was used for the inference and analysis of intercellular communication [22]. CellChatDB.human, ligand-receptor interaction databases for human, was used to evaluate the major signaling inputs and outputs between cell subtypes.

The upregulated ligand-receptor interaction was calculated with function “netVisual_bubble”. “netAnalysis_computeCentrality” function was used to discover dominant senders, receivers, mediators, and influencers in the intercellular communication networks through the centrality metrics from graph theory.

Construction of multilayer inter- and intra-cellular signaling networks

To elucidate the complex signaling mechanisms underlying PAH, we employed two complementary methods: scMLnet and MultiNicheNet. These approaches aimed to construct multilayer signaling networks and explore intercellular communication, providing insights into stage-specific regulatory dynamics and cell-type-specific signaling pathways. We employed scMLnet to construct multilayer inter- and intra-cellular signaling networks from scRNA-seq datasets [23]. Subsampled and annotated scRNA-seq data were formatted into input files containing cell barcodes and clusters. Using curated ligand-receptor, receptor-transcription factor (TF), and TF-target interaction libraries, we applied the RunMLnet function to infer signaling networks for all cell type pairs, with thresholds of pval = 0.05 and logFC = 0.15. The DrawMLnet function visualized these networks, highlighting stage-specific intercellular communication and regulatory dynamics. Python 3.13 was configured for graphical output, and all interactions were validated against curated databases.

We also employed MultiNicheNet to infer and analyze intercellular communication in pulmonary arterial hypertension (PAH) using single-cell RNA sequencing (scRNA-seq) data [24]. To ensure data quality, we filtered genes using thresholds of min_cells = 1, fraction_cutoff = 0.01, and min_sample_prop = 0.1, retaining genes expressed in at least one cell, with a minimum fraction of 1%, and present in at least 10% of samples. Ligand activity was calculated, and ligand-target gene interactions were predicted using thresholds of logFC_threshold = 0.50 and p_val_threshold = 0.05. The MultiNicheNet pipeline integrated curated ligand-receptor and ligand-target interaction databases to infer ligand activity, predict downstream target genes, and identify significant ligand-receptor-target pathways. This approach provided comprehensive insights into.

Correlation analysis

Spearman correlation is a rank correlation coefficient that captures the direction and strength of the relationship between two variables. The “rcorr” function from the ‘Hmisc’ package was employed to conduct Spearman correlation analysis. Specifically, the top ten highly expressed genes (with fold change > 2 and adjusted p-value < 0.01) in the Seurat cell clusters across all samples were analyzed to determine the correlation strength between different populations and each cell cluster, and the “corrplot” function from the ‘corrplot’ package (v0.92) was used to visualize.

Animal studies

The institutional Animal Research Committee approved all animal procedures, which were performed according to National Institutes of Health guidelines. ARRIVE reporting guidelines were used. Male mice (purchased from Shanghai SLAC Laboratory Animal) 7–8 weeks of age were housed in hypoxia (10% O₂) for six weeks, or in room air for six weeks. Pentobarbital sodium was dissolved in normal saline and administrated intraperitoneally at 50 mg/kg for anaesthesia to carry out following studies. The lungs were perfused with PBS, the left lung fixed, and sections embedded in paraffin for immunohistochemistry and immunofluorescence. The right lung was snap-frozen in liquid N₂ and kept at −80 °C.

Vascular remodeling analysis

After paraffin embedding and sectioning, the lung slides (4 μm thickness) were stained with hematoxylin and eosin (H.E.) and Masson for morphological analysis to assess pulmonary arterial wall thickness.

Multicolor fluorescence imaging

Formalin-fixed, paraffin-embedded lung sections underwent dewaxing and heat-induced antigen retrieval at pH 6. Multiplex immunofluorescence staining was conducted utilizing the Opal kit (Akoya). This process involved sequential application of primary antibodies (GPNMB, CD44, α-SMA, and F4/80), detection reagents (either Opal Polymer HRP from Akoya or Immpress HRP Polymer from Vector Laboratories), and fluorescence signal development using Opal dyes from Akoya or AlexaFluor TSA from Thermo Fisher Scientific, followed by antibody removal, all in accordance with the manufacturer's protocols. DAPI (Thermo Fisher Scientific) was employed as a nuclear counterstain.

Statistical analysis

Data processing and statistical analyses were performed using R software version 4.1.3 and GraphPad Prism version 9.0 (USA). Data are presented as the mean ± standard error of the mean (SEM). Differences in continuous variables between two groups were assessed using either the Wilcoxon test or the t-test, while differences among multiple groups were evaluated using a one-way analysis of variance (ANOVA). All p-values were calculated using a two-tailed approach, with statistical significance set at p < 0.05.

Results

Integrated scRNA‐seq analysis quantified the diversity of major cell populations in pulmonary hypertension samples

To generate a comprehensive map of immune landscape of human PAH at single-cell resolution, four scRNA-seq datasets (GSE254617, GSE210248, GSE233189, and GSE228643) were analyzed by an integrated bioinformatics method (Fig. 1A and Supplementary Material 7). Results of clustering of single-cell transcriptome data for different cell types were presented (Fig. 1B, Supplementary Material 8, Figure S1A, Figure S1B and Figure S1C). In the UMAP, different colors indicate different cell types, including smooth muscle cells and fibroblasts (SMC_Fib), Endothelial cells, Epithelial cells, B cells, T cells, and Myeloid cells. The proportional changes of certain cell types in the disease state are visualized. This analysis provides a basis for understanding the potential role of different cells in the pathogenesis of PAH (Fig. 1C, Figure S1D and Figure S1E). Then, we showed the results of differential gene enrichment analysis for each large population of cell types (Fig. 1D, Table S1 and Table 1). Then, Fig. 1E illustrated the relative proportions of each cell type in PAH patients, highlighting notable changes in cellular composition associated with PAH pathology. These data demonstrate the transcriptomic features of different cell types in the lung tissues of PAH patients.

Fig. 1

scRNA-seq analysis of PAH lung tissue reveals cellular heterogeneity and key gene expression differences between PAH patients and healthy controls. A Schematic representation of the scRNA-seq analysis process. Lung tissue samples were collected from PAH patients and healthy controls, and multiple publicly available datasets (GSE254617, GSE210248, GSE233189, and GSE228643) were integrated for analysis. Single cells from these samples were sequenced to explore gene expression patterns across various cell types, enabling a detailed investigation of cellular and molecular changes associated with PAH. B UMAP plot showing clusters of different cell types identified in the lung tissue samples, with each color representing a unique cell type: SMC, Fib, endothelial cells, epithelial cells, B cells, T cells, and myeloid cells. C UMAP plots (bottom) show the distribution of cell types and their respective proportions in PAH and control groups, highlighting the shifts in cell type composition within diseased tissue. D Heatmap depicting differentially expressed genes across cell types, with pathways significantly enriched in PAH patients annotated on the right. Enrichment terms include pathways related to cell migration, inflammation, vascular remodeling, and immune response, suggesting critical molecular processes that may be driving PAH pathogenesis. The bar plots on the right indicate the number of differentially expressed genes associated with each pathway, grouped by cell type. E Cell proportion analysis: The table illustrates the relative proportions of each cell type in PAH and Con, highlighting notable changes in cellular composition associated with PAH pathology

Table 1 Pathway enrichment results for the total population

Comprehensive single-cell transcriptomic analysis reveals T cell heterogeneity and functional dynamics in PAH

Next, we focused on the change of T cell subsets during PAH development. In Fig. 2A, Figure S2A and Figure S2B, we found that the subset classification of T cells using UMAP dimensionality reduction technology, labeling multiple subsets. These cells include CD8⁺ Effector T cells (CD8⁺ Effector), CD4⁺CCR6⁺ T cells (CL_CD4 + CCR6⁺), memory T cells (Activated T cells), CD4⁺ and CD8⁺ Naive T cells (CD4⁺ Naive and CD8⁺) Naive). Meanwhile, the changes in the number and proportion of T cell subsets in the disease group (including PAH and related disease subtypes) and the Con were presented in Fig. 2B and C. These results revealed the remodeling characteristics of the T cell population, and further supported the close link between PAH and immune activation and dysfunction of T cell subsets. Then, we showed the UMAP distribution of T cells in different disease states (Con, IPAH, PHHF, and SSc-PH) and the curves of T cell numbers in each group as a function of disease progression (Fig. 2D). Pseudo-time series analysis of T cell subsets by the two analysis methods of Monocle2 and Monocle3, revealed the dynamic change trajectories of different T cell subsets in the time axis (Fig. 2E and Figure S2C). In the pseudo-time series analysis, different colors represent different T cell subsets, including CD4⁺ Naive T cells, CD8⁺ Naive T cells, Activated T cells, CD8⁺ Effector T cells, and so on. We found T cells undergo a significant dynamic differentiation during the pathological process of PAH, which is manifested as a gradual transition from Naive T cells to Activated T cells and Effector T cells. The significantly differentially expressed genes in each subgroup were visualized by heat map, and the key pathways of gene function enrichment analysis (GO) were enumerated (Fig. 2F, Table 2 and Figure S2D). SCENIC analysis in Fig. 2G revealed the regulatory network of specific TFs and their Activity changes (AUC, Activity Units of Regulons) in T cell subsets. CellChat analysis of T-cell subtypes revealed intercellular communication network and we observed the altered expression levels of GPNMB, PIEZO1, and PDCD1 in control and PAH groups (Figure S2E and Figure S2F). This analysis provides important insight into the immune mechanisms underlying PAH and suggests that targeting cytotoxic regulation may be a potential therapeutic strategy.

Fig. 2

Comprehensive single-cell transcriptomic analysis reveals T cell heterogeneity and functional dynamics in PAH. A UMAP plot showing the clustering and annotation of T cell subsets in PAH and healthy control samples, including CD8 + Effector, CL_CD4 + CCR6 + , Activated T, CD4 + Naïve, CD8 + Naïve, and NK cells. Each cluster represents a distinct functional subset, highlighting the heterogeneity within the T cell population. B and C Bar and pie charts displaying the proportions and counts of different T cell subsets in the PAH group and healthy controls. Activated T cells and CD8 + Effector T cells show significant increases in the PAH group, while Naïve T cells decrease, indicating enhanced immune activation in PAH. D UMAP plots displaying T cell distribution across different disease states (Con, IPAH, PHHF, SSc-PH), accompanied by line charts showing the proportion of T cell subsets in each condition. Disease-specific shifts in T cell composition are observed, with notable changes in IPAH and SSc-PH. E Monocle2 and Monocle3 pseudotime trajectories of T cell subsets, showing dynamic transitions from Naïve T cells to Activated T cells and CD8 + Effector T cells. Violin plots illustrate gene expression trends along the pseudotime axis, reflecting the functional evolution of T cells in PAH. F Heatmap of differentially expressed genes across T cell subsets, with enriched pathways annotated, including leukocyte adhesion, inflammatory response, and immune activation. Key pathways suggest that T cells contribute to inflammation, tissue remodeling, and disease progression in PAH. G Heatmap of TF activity across T cell subsets, highlighting TFs with altered activity in PAH, such as NFATC1, HIF1A, BATF, STAT1, and FOXP3. These TFs regulate immune activation, metabolic adaptation, and inflammatory processes

Table 2 Pathway enrichment results for T cells

Comprehensive single-cell transcriptomic analysis reveals myeloid cell heterogeneity and functional dynamics in PAH

As exhibited in Fig. 3A, Figure S3A, Figure S3B and Figure S3C, we showed and annotated different subsets of myeloid cells using UMAP dimensionality reduction, including Monocyte, LSAMP⁺ Macro, MARCO⁺ Macro, APOE⁺ Macro, DC1, DC2 and Mast. The number and relative proportions of different myeloid cell subsets in the PAH patient and Con groups were listed in Fig. 3B and C. The differentiation trajectory of myeloid cells from Monocyte to macrophage (such as APOE⁺ Macro and MARCO⁺ Macro) and dendritic cells (DC1 and DC2) by Monocle pseudo-timing analysis was presented (Fig. 3D and Figure S3D). In addition, line plots in Fig. 3E showed the temporal expression trends of GPNMB in specific subsets (APOE⁺ Macro) between PAH and control groups. GPNMB showed increased expression along pseudotime in PAH, reflecting its crucial involvement in immune activation and tissue remodeling. As shown in Fig. 3F and Figure S3E, we found that GPNMB was increased in different subsets of myeloid cells from PAH groups. The results of differential gene enrichment analysis for myeloid cell subsets were indicated in Fig. 3G and Table 3. The transcription factor activities of myeloid cell subsets using SCENIC analysis was performed and the heat map shows the differences in regulatory activities of different transcription factors such as SPI1, IRF8, and HIF1A between PAH and controls in Fig. 3H. For another, Figure S3F manifested the annotation of cell clusters with specific cell type names, clarifying their biological identities. These reveals the significant changes and potential functions of myeloid cell subsets in PAH patients.

Fig. 3

Comprehensive single-cell transcriptomic analysis reveals myeloid cell heterogeneity and functional dynamics in PAH. A UMAP plot showing the clustering and annotation of myeloid cell subsets, including Monocyte, APOE⁺ Macro, MARCO⁺ Macro, LSAMP⁺ Macro, DC1, DC2, and Mast cells. Each subset represents a functionally distinct population of myeloid cells in PAH and healthy controls. B and C Balloon and bar plots displaying the proportions and counts of myeloid cell subsets in PAH (Case) and healthy controls (Con). APOE⁺ Macro and MARCO⁺ Macro are significantly enriched in PAH, while Monocyte and DC1 proportions decrease, reflecting shifts in myeloid cell composition during disease progression. D Pseudotime trajectories revealing the dynamic differentiation of Monocyte into APOE⁺ Macro, MARCO⁺ Macro, and DC subsets. Bubble plots show the temporal expression patterns of key genes (e.g., APOE, TREM2, and GPNMB), indicating functional evolution during disease progression. E Line plots showing the temporal expression trends of GPNMB in specific subsets (APOE⁺ Macro) between PAH and control groups. Both genes show increased expression along pseudotime in PAH, reflecting their involvement in immune activation and tissue remodeling. F UMAP analysis illustrated the expression levels of GPNMB across myeloid cells in PAH and control group. G Heatmap of differentially expressed genes across myeloid cell subsets, with functional pathway enrichment annotations. APOE⁺ Macro and MARCO⁺ Macro are enriched in lipid metabolism and immune regulation pathways, while DC1 and DC2 are enriched in antiviral response and anti-inflammatory pathways, indicating subset-specific functions in PAH. H Heatmap of transcription factor (TF) activity (AUC scores) in myeloid cell subsets, highlighting key TFs such as SPI1, IRF8, and HIF1A. SPI1 drives inflammation in Monocyte and macrophage subsets, IRF8 regulates antiviral responses in DC subsets, and HIF1A mediates metabolic adaptation in APOE⁺ Macro and MARCO⁺ Macro

Table 3 Pathway enrichment results for myeloid cells

Comprehensive analysis of SMC and FB subsets in PAH

Subsequently, we showed the results of subgroup division of smooth muscle cells and fibroblasts. Multiple subsets are annotated by UMAP dimensionality reduction visualization (Fig. 4A and Figure S4A). These included Alveolar FB, Pi16⁺ FB, Contractile SMC-1, contractile SMC-2, and Synthetic SMC and Myofibroblast SMC. Each subpopulation represents a specific functional state or differentiation pathway, revealing the heterogeneity of these cells in PAH. Figure 4B included bar graphs and grouped stacked bar graphs showing the changes in proportions of different subgroups in PAH patients and Con. The ratio of Contractile SMC-1 and Contractile SMC-2 was higher in healthy controls, whereas Synthetic SMC and Myofibroblast SMC were significantly elevated in PAH patients. Figure 4C and Figure S4B illustrated the differential gene functions of SMC and FB subsets were determined by GO enrichment analysis in Fig. 4C and Figure S4B. Then, in Fig. 4D, we showed the marker gene expression of each subgroup, with bubble size indicating the expression ratio and color indicating the expression intensity. The key marker genes included ACTA2, TAGLN, and VCAN. ACTA2 and TAGLN were highly expressed in Contractile SMC. VCAN is highly expressed in Synthetic SMC, suggesting that these cells may contribute to vascular remodeling through ECM organization. Next, Fig. 4E and Figure S4E showed the expression distribution of GPNMB in different cell subsets and compares the differences in the PAH and Con. The activity of key TFs in different SMC and FB subpopulations were analyzed using SCENIC analysis in Fig. 4F. In addition, Fig. 4G demonstrated the quasi-temporal expression dynamics of GPNMB in three SMC subsets. The graphs show the expression trends of these genes in PAH.

Fig. 4

Comprehensive analysis of SMC and FB subsets in PAH. A UMAP plot showing the clustering and annotation of SMC and FB subsets, including Alveolar FB, Pi16⁺ FB, Contractile SMC-1, Contractile SMC-2, Synthetic SMC, and Myofibroblast SMC. Each subset represents distinct functional states, highlighting the heterogeneity in PAH and healthy controls. B Bar and stacked bar plots illustrating the proportions of each subset in PAH and Con. Synthetic SMC and Myofibroblast SMC show significant increases in PAH, while Contractile SMC populations are reduced. C Heatmap showing the enriched pathways for each subset based on GO and KEGG analysis. Synthetic SMC is enriched in ECM organization and angiogenesis pathways, while Myofibroblast SMC is enriched in smooth muscle cell proliferation and inflammatory signaling pathways. D Bubble plot displaying the expression patterns of key marker genes (e.g., ACTA2, TAGLN, VCAN) in different subsets. ACTA2 and TAGLN are highly expressed in Contractile SMC, while VCAN is upregulated in Synthetic SMC, reflecting functional specialization. E Violin plots showing the expression levels of GPNMB across subsets and between PAH and control groups. GPNMB is highly expressed in Synthetic SMC and Myofibroblast SMC, indicating their roles in ECM remodeling and mechanosensation. F Heatmap illustrating transcription factor activity (AUC scores) in SMC and FB subsets. SRF is active in Contractile SMC, regulating contraction-related genes, while TWIST1 and HIF1A are active in Synthetic SMC and Myofibroblast SMC, promoting ECM remodeling and adaptation to hypoxia. G Line plots showing the temporal expression patterns of GPNMB across three SMC subsets (Contractile SMC-1, Contractile SMC-2, Synthetic SMC). Both genes show increased expression in Synthetic SMC at later pseudotime stages, reflecting their roles in pathological progression

Cell-chat and niche-net analyzes reveal altered intercellular communication and signaling pathways in PAH

Clearly, the maintenance of cellular phenotypes involves coordinated actions with many regulatory factors, including the regulation of internal cellular TFs and external cell-to-cell communication. The intercellular communication network was mapped by CellChat analysis, demonstrating the interaction patterns and strength changes between different cell types (Fig. 5A). MIF and IL-1 signaling pathways were significantly enhanced, while some immunosuppressive pathways are weakened in patients with PAH. The heatmap demonstrated the strength of the interaction between different cell types as source and target (Fig. 5B, Figure S5A). Signal pattern analysis showed the overall signal intensity in the healthy control group and the PAH group, respectively (Fig. 5C and D). By presenting the complex ligand-receptor-TF-target gene network, the figure uncovers key mechanisms of cell communication in PAH, with a particular emphasis on the critical role of CD44 in this process (Fig. 5E). The expression of MIF, CD74, CXCL4, and CD44 was compared across cell types in control and PAH groups, revealing their critical roles in PAH pathology. CD44 was upregulated in synthetic SMCs and endothelial cells, highlighting its role in cell adhesion and migration. These results underscore the coordinated roles of these molecules in inflammation and vascular remodeling in PAH (Fig. 5F).

Fig. 5

CellChat analyses reveal altered intercellular communication and signaling pathways in PAH. A Cell–cell communication network highlights enhanced interactions in PAH, particularly involving macrophages and Synthetic SMC. Information flow comparison reveals upregulated inflammatory pathways (e.g., MIF, IL-1) and downregulated immunosuppressive pathways in PAH. Relative Information Flow of Various Molecules in Con and PAH Groups. B Overall signaling patterns show dominant inflammatory signals in PAH compared to controls. C Pathway analysis highlights upregulated MIF signaling. D Heatmap comparing MIF signaling distribution between control and PAH groups. E Ligand-Receptor-Transcription Factor-Target Gene Network. F Comparison of MIF, CD74, CXCL4, and CD44 expression across cell types in Con and PAH Groups

PAH mouse model was built to validate the expression pattern of GPNMB in macrophages

Consistent with the original literature, we used the hypoxia model to induce PAH (Fig. 6A). As shown in Fig. 6B and C, H&E and Masson staining were used to assess the vascular remodeling of mice lung tissues. To elucidate the relationship between SMCs and inflammatory infiltration in PAH, GPNMB-F4/80 and CD44-α-SMA immunofluorescent staining (Fig. 6D) were performed, which indicated the potential cell chat between GPNMB⁺ macrophages and CD44⁺ SMCs.

Fig. 6

In vivo hypoxia-PAH model was established to validate the cell chat between GPNMB⁺ macrophages and CD44⁺ SMCs. A Schematic diagram of the animal experimental design of PAH. B Representative images of H&E staining among groups. Scale bar = 20 μm. C Representative images of Masson staining among groups. Scale bar = 20 μm. D Representative images of GPNMB (purple), F4/80 (pink), α-SMA (red) and CD44 (pink) mmunofluorescent staining. Scale bar = 20 μm. The significance of difference between the means of groups was evaluated using one-way ANOVA; n = 6 per group

Discussion

PAH is a life-threatening condition characterized by elevated blood flow within the pulmonary circulation [25]. This pathology affects the entire vascular layer and is marked by a progressive increase in pulmonary vascular resistance and the development of lesions in the small pulmonary vascular plexus. At present, there is no definitive cure for PAH. Available molecular targeted therapies, including endothelin receptor antagonists, phosphodiesterase inhibitors, calcium channel blockers, and cyclic prostanoids, offer symptomatic relief but do not provide a complete resolution of the disease. Therefore, the contradictory results of the previous extensive anti-inflammatory treatments require a more detailed characterization of the complex immune microenvironment of PAH [26, 27]. In our study, by integrating four human scRNA-seq data- sets and utilizing state-of-the-art analysis tools to avoid the limitations of a single study, we exhibited that the PAH microenvironment is in fact more heterogeneous. Briefly, we conducted a bioinformatics analysis of samples in two groups of health control or PAH patients. We elucidate the functional characteristics and potential regulatory interactions of several cell subpopulations that have not been previously documented in similar research, with a particular focus on GPNMB⁺ macrophages. A comprehensive landscape of cell communications was constructed a at the single-cell resolution, which might significantly advance the development of personalized diagnostic and therapeutic strategies for PAH.

Studies have revealed changes in the composition of T cell subsets, differential expression of key genes, and dynamic changes in disease states in PAH [28]. For example, Tregs play a crucial role in preventing autoimmunity, and there is a growing body of evidence indicating their impaired function within the inflammatory environment of pulmonary arterial hypertension (PAH) [29]. Dysfunctional Treg activity is significantly associated with an increased susceptibility to PAH in both animal models and human patients. In our investigation, we systematically presented an in-depth analysis of T cell subsets in PAH patients by single-cell transcriptomes. We indicated the important roles of specific T cell subsets including activated T and CD8⁺ Effector in PAH. We reported T cell subsets in PAH are involved in disease progression by regulating inflammatory response, enhancing cell adhesion and migration, and metabolic adaptation, which provides potential targets for exploring therapeutic strategies for immune abnormalities. Our analysis provides important insight into the immune mechanisms underlying PAH and suggests that targeting cytotoxic regulation may be a potential therapeutic strategy.

Macrophages are highly versatile and heterogeneous cells present in all bodily tissues. In response to specific conditions, such as tissue injury, additional macrophage subtypes are recruited to the tissue from the monocyte reservoir [30]. This reveals the significant changes and potential functions of myeloid cell subsets in PAH patients: The increased proportion of APOE⁺ Macro and MARCO⁺ Macro and the high expression of key genes such as GPNMB, APOE and TREM2 suggest that these subsets may be involved in the pathological process of PAH by regulating lipid metabolism, inflammatory response and tissue repair. GPNMB is encoded by the GPNMB gene located at locus 7p15 and has 2 isoforms [31]. After transcription, GPNMB is directed to the cell membrane, in which most of its length is located in the extracellular domain a. Several cell types are reported to express GPNMB: these include phagocytes (dendritic cells and macrophages), osteoclasts and melanocytes [32, 33]. Various reports have reported the co-function of GPNMB and macrophages. Pulmonary low-grade inflammation results in the development of ApoE-dependent monocyte-derived alveolar macrophages (ApoE⁺CD11b⁺ AMs), which express high levels of GPNMB [34]. Meanwhile, GPNMB is closely associated with inflammatory responses in multiple diseases. Some scholars studying acute kidney injury found that GPNMB can facilitate transformation of macrophages into M2 type to balance the inflammatory response by M1 type, thereby exerting an anti-inflammatory effect [35]. Our data found GPNMB was increased in different subsets of myeloid cells from PAH groups. Previous study has reported that macrophage-derived GPNMB trapped by fibrotic extracellular matrix can induce pulmonary fibrosis [13]. GPNMB-high macrophages contribute to PN-MES transition and impede T cell activation in glioblastoma [36]. Also, anti-GPNMB antibody has been assessed in phase I/II trials in several cancers, making GPNMB an attractive target in cancer treatment [37]. However, the significance of GPNMB or GPNMB⁺ macrophages to PAH development needs more investigation.

Additionally, we found that the heterogeneity and dynamic functional changes of smooth muscle cells and fibroblasts in PAH. The proportion of synthetic SMC and myofibroblasts is significantly increased in PAH, and the functional changes drive vascular remodeling and inflammatory response. GPNMB are highly expressed in Synthetic SMC and Myofibroblast SMC, which may drive the vascular remodeling and fibrosis process of PAH through ECM remodeling and mechanical signal perception. The abnormal expression of two immunomodulatory genes, PDCD1 and CD274, suggests that smooth muscle cells acquire immune escape ability in the chronic inflammatory microenvironment of PAH and may inhibit T cell function to adapt to the pathological environment. Synthetic SMC not only participates in ECM remodeling in PAH, but also becomes the core of inflammatory regulation by enhancing MIF-CD44 and IL-1 signaling. The weakening of immunosuppressive signals including CXCL12 and other anti-inflammatory signals aggravates the inflammatory microenvironment, which provides a new explanation for the pathological progression of PAH.

Several limitations of this study are subsequent: (1) Due to sample limitations, we just revealed a general disease-related immune heterogeneity, which failed to consider heterogeneity of the diversity of risk factors, degree of disease progression and variability in treatment methods of PAH. (2) The loss of spatial information during cell dissociation results in inadequate prediction of cell communication. (3) Future single cell multin\omics studies with more elaborately designed experiments will further elucidate the pathogenesis of PAH. (4) The function of macrophages release signals for other immune cells recruitment before entering PAH microenvironment remains uninvestigated. (5) Interference of GPNMB in animal models is warranted to investigate their direct effects on PAH immunosuppressive phenotype.

Conclusions

In conclusion, we revealed the immune landscape in PAH by consulting integrating multiple datasets in an unbiased manner. These data may provide insights into the function and modulation of PAH pathogenesis and progression, which could prove valuable immunotherapeutic strategies specifically targeting PAH.