Introduction

Classical view of vascular smooth muscle cells as resident cells mediating vasoactive responses and adventitial fibroblasts as structural support cells producing of fibrillar collagen has been appended by single cell transcriptomics studies. Single cell resolution has helped to define the tissue specific heterogeneity within and delineate the distinctions of vascular resident cell populations to other smooth muscle and fibroblast populations ^1–4. Additional layer of cellular complexity and heterogeneity is acquired through the disease process. Vascular smooth muscle cells and adventitial fibroblasts preserve their lineage identity in animal models of systemic and pulmonary vascular disease ^5–10. However, diseased cells consistently display altered expression profiles compared to their healthy counterparts ^4,9,10.

Pulmonary vascular disease is an umbrella term encompassing functional and morphological alterations in lung vasculature, either as a standalone condition or, more often, accompanying chronic lung and heart disease ¹¹. We have previously shown that pulmonary vascular disease and remodeling of pulmonary artery compartment is associated with a skewed cellular communication centered on smooth muscle cells and adventitial fibroblasts ⁴. Moreover, molecular alterations identified in situ and freshly isolated cells are at least partially preserved and propagated in vitro. For example, pulmonary artery smooth muscle cells (PASMC) and pulmonary artery adventitial fibroblasts (PAAF) isolated from remodeled pulmonary arteries display metabolic and proinflammatory reprograming that was linked with the progression of pulmonary vascular remodeling ^12–14. On transcriptional and proteomic level, PASMC display characteristics of phenotypic shift, namely downregulated expression of contractile machinery and upregulation of extracellular matrix components ^4,15. PAAF were shown to direct the behavior of immune cells and serving as a homing place for the immune cells during vascular remodeling ^16–18. Cumulatively, results imply a cell-type specific communication network that is perturbed in the vascular remodeling. It is however unclear to which extent are the changes in expression profiles ultimately mirrored phenotypically.

In the current study, we investigated early passage PASMC and PAAF sourced from the same pulmonary artery of healthy, downsized lung donors or patients with idiopathic form of pulmonary arterial hypertension (IPAH). We have combined integrative omics and network-based analysis with targeted phenotypic screens to identify common and cell-type specific differences. Our aim was to uncover significant phenotypic differences upon pulmonary vascular remodeling and identify a set of key molecular correlates that could serve as surrogate marker defining healthy and diseased cell states in a cell-type specific manner. Our results provide insight into bidirectional cell communication between PASMC and PAAF and the critical role of heterotypic cell-to-cell interaction in cell state maintenance and transition in pulmonary artery compartment.

Results

Distinct omics profiles underlie specialized homeostatic functions of normal PASMC and PAAF

The extent of PASMC and PAAF functionality beyond classical contractile and structural roles is unclear. At first, we focused on defining the normal, healthy cell states by using very early passage (p1) PASMC and PAAF isolated from same (source-matched) pulmonary arteries (PA) of healthy, downsized donors. We performed deep molecular characterization of RNA and protein expression profiles (Figure 1A, Data S1, S2). Principle component analysis confirmed that molecular distinction between two cell types was preserved in early passage on both transcriptomic and proteomic level (Figure 1B). It also revealed higher interpersonal heterogeneity within PASMC samples, while source-matched PAAF represented as a more homogenous group (Figure 1B). We identified a set of significantly enriched genes/proteins in PASMC (Regulator Of G Protein Signaling 5 - RGS5, Insulin Like Growth Factor Binding Protein 5 - IGFBP5) and PAAF (Alcohol Dehydrogenase 1C - ADH1C, Scavenger Receptor Class A Member 5 - SCARA5, Complement Factor D - CFD) (Figure 1C). We previously described RGS5 as cluster-enriched marker for contractile SMC, while SCARA5 and CFD were identified as human pulmonary adventitial fibroblast markers ^4,19. We validated our omics results and mapped location of selected enriched targets by multicolor staining of human lung tissue (Figure 1D, E). RGS5 and IGFBP5 positivity overlapped predominantly with ACTA2 (alpha smooth muscle) stain, while ADH1C, SCARA5 and CFD were observed mostly in adventitial region marked by PDGFRA (Platelet-Derived Growth Factor Receptor Alpha) (Figure 1D, E). Comparative analysis of transcriptomic and proteomic data sets revealed a strong, but not complete level of linear correlation between the gene and protein expression profiles (Figure S1). We therefor decided to use an integrative dataset and analyzed all significantly enriched genes and proteins (-log₁₀(P)>1.3) between both cell types (Figure 1F). A gene set enrichment analysis (GSEA) on gene ontology (GO) terms for biological processes was used to infer functional implications of observed expression profiles. The analysis confirmed polarized and distinct cell-type functional associations, revealing vascular development and cell migration processes for PASMC, and endoplasmatic reticulum protein targeting for PAAF. Additionally, PASMC were associated with processes linked to sterol transport, O-glycosylation, and nerve growth factor signaling (Figure 1F). PAAF on the other hand showed enrichment in mitochondrial metabolic processes and heavy ion detoxification (Figure 1F). Applying an additional cutoff by taking into consideration only strongly expressed genes/proteins, we performed an overrepresentation analysis of GO biological processes (Figure 1G). Among top10 processes, there was a clear separation and functional distinction between PASMC and PAAF, with overlap in angiogenesis and adhesion (Figure 1G). Contrasting the canonical functional roles, it was PASMC, rather than PAAF, that showed enrichment for extracellular matrix (ECM) organization (Figure 1G), while PAAF showed enrichment for heavy ion response, lipid homeostasis and cell differentiation (Figure 1G). The enrichment of lipid-associated processes in PAAF prompted us to validate the observed results phenotypically. We performed a fuel usage assay using Seahorse metabolic analyzer. Real-time monitoring of oxygen consumption under basal condition and upon specific metabolic inhibitors revealed that PASMC could equality utilize glucose/pyruvate and long chain fatty acids as a fuel source for oxidative metabolism (Figure 1H, I). In contrast, upon blocking of fatty acid utilization in Krebs cycle, PAAF had stronger drop in basal respiration rate while PASMC were less sensitive to this maneuver (Figure 1I). PAAF displayed additional significant metabolic differences to PASMC (Figure S2). PAAF had lower oxygen consumption and extracellular acidification rates (Figure S2). Ratio of oxygen consumption to extracellular acidification (OCR/ECAR) revealed that PAAF have relative preference for oxidative phosphorylation compared to PASMC (Figure S2), in line with observed enrichment for mitochondrial respiration processes (Figure 1F). This was further supported by smaller glycolytic capacity (Figure S2) and smaller mitochondrial proton leak (Figure S2). Cumulatively, molecular and phenotypic profiling revealed specialized homeostatic roles of normal, healthy PASMC and PAAF.

Figure 1.

IPAH-defining cell states of PASMC and PAAF

We next set out to delineate expression changes defining the diseased, remodeling state. Omics profiling was performed on source-matched PASMC and PAAF isolated from patients with end stage idiopathic pulmonary arterial hypertension (IPAH) in comparison to donor samples (Figure 1A). PCA analysis of transcriptomic and proteomic profiles revealed a clear sample separation based on cell type (Figure 2A). However, sample separation based on disease state was less obvious (Figure 2A). We used the supervised multivariate Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA) method to model the extent by which cell-type and disease-state parameters have on sample separation. In the first model iteration, we considered either cell type or disease state as input parameter. OPLS-DA modeling found significant differences in the global transcriptomics and proteomics profiles based only on the cell type, but not disease state (Figure 2B), validating PCA analysis observations. For the second model iteration, we took both cell type and disease state simultaneously as input parameters. In this OPLA-DA model, only in the case of transcriptomic dataset there was significant separation between all four groups (Figure 2C). These results imply that remodeling process results in substantial transcriptional changes (cell-state), but that each cell type retains their original cell lineage designation, in line with results from fate mapping and fresh cell single cell transcriptomics ^4,6,20. Gene ontology analysis of selected gene sets from combined transcriptomic and proteomic data (significantly regulated genes/proteins, Figure 2D), showed a cell-type specific footprint of molecular changes upon IPAH (Figure 2E). IPAH-PASMC were characterized by upregulation of nucleotide-sugar and lipid biosynthesis, Golgi vesicle transport and endoplasmatic reticulum (ER) stress processes, while simultaneously downregulating of cytoskeletal components, RNP export and p53 signaling (Figure 2E). IPAH-PAAF upregulated genes and proteins involved in cell proliferation, DNA repair and downregulated expression of NF-kB pathway components, NK-cell activation and transcriptional stress response (Figure 2E). Enrichment for activated sugars and Golgi transport process is associated with increased synthesis and secretion of glycosylated proteins, such as ECM components. We therefore measured production of proteoglycans by IPAH-PASMC in situ. Glycosaminoglycan (GAG)-containing ECM components were visualized using Alcian blue staining and showed enhanced medial staining compared to donor and prominent signal in neointimal PA compartment (Figure 2F). Quantitative measurements using 1,9-dimethylmethylene blue assay on isolated PA confirmed different distribution of GAG-content in IPAH compared to donor PA (Figure 2G). In line with observed enrichment for DNA replication and cell cycle phase regulation (Figure 2E), IPAH-PAAF displayed higher proliferative capacity compared to donor-PAAF, upon both stimulated and unstimulated conditions (Figure 2H). IPAH-PASMC in contrary showed lower proliferative rate compared to donor-PASMC, corresponding to our previous observation of G1 cell cycle accumulation of IPAH-PASMC ⁴ (Figure 2H). The observed phenotypic effect is strongest in very early passage and is progressively lost with passaging (Figure 2I). In vitro results mirror indirect determination of cellular proliferation rates in situ using PCNA staining of human lung tissue. We observed frequent PCNA+ fibroblasts, identified by PDGFRA staining, in the adventitial layer of PA from IPAH patients (Figure 2J).

Figure 2.

Shared IPAH expression disturbances in PASMC and PAAF converge on mitochondrial dysfunction

Distinct remodeling expression profiles and cell behavior, prompted us to search if there are any common molecular perturbations shared between IPAH-PASMC and IPAH-PAAF. Number of co-directionally regulated genes and proteins were limited to 47 (Figure 3A). Further selection based on significance and fold change provided a list of most significantly regulated molecules in IPAH that are shared between PASMC and PAAF (Figure 3B, highlighted corners bottom left and top right). Using this as an input for network analysis based on the STRING database, we built the underlying gene interaction network and identify shared upstream regulating factors (Figure 3C). Only a fraction of commonly regulated genes and proteins were encompassed in the regulatory network (Figure 3C). For the remaining elements, it was not possible to construct an interaction map, implying regulation by independent mechanisms and upstream regulators. We therefore focused further on those molecules that were encompassed in the interaction network. Interestingly, their regulatory nodes converged on the immuno-inflammatory and survival-apoptosis regulators and effectors (Figure 3C). Specifically, a majority of the nodes building the interaction network belonged to a limited set of genes comprising mostly cytokines/chemokines (IL1B, IL2, IL4, IL6, IL10, TNF, CCL2, CXCL8) and inflammatory-stress response signal transduction pathways (JNK, AKT, NF-κB, STAT3, MYC, JUN) (Figure 3C). Additional pathway enrichment analysis based on the BioPlanet database, uncovered regulation of apoptosis as a shared pathway by both cell types in IPAH (Figure 3D). Shared SOD2 (mitochondrial superoxide dismutase 2) downregulation prompted us to investigate the mitochondrial function and redox defenses. We performed tetramethylrhodamine methyl ester (TMRM) staining for the detection of mitochondrial content and potential. Cellular mitochondrial content was lower both in IPAH-PASMC and IPAH-PAAF upon removal of serum and growth factors (starvation) (Figure 3E). However, the remaining mitochondria in IPAH cells were hyperpolarized (Figure 3F). Mitochondrial membrane hyperpolarization is linked with increased reactive oxygen species (ROS) production. We measured cellular ROS production using CellRox Deep Red dye at basal levels and upon PDGF-BB stimulation revealing that only IPAH-PASMC significantly upregulate ROS production (Figure 3G). Although mitochondrial dysfunction seems as a common underlying disturbance between cell types in IPAH, the downstream consequences appear to be cell-type selective.

Figure 3.

PASMC and PAAF acquire divergent disease-cell states

Contrasting the limited overlap in regulated elements between PASMC and PAAF, vast majority of IPAH transcriptomic and proteomic changes (more than a thousand genes/proteins) were uniquely regulated between the cell types, including substantial number of inversely regulated elements (Figure 4A). PASMC possessed more pronounced changes than PAAF as shown by the number of differentially regulated genes and proteins in IPAH (Figure 4B). A STRING based interaction network, overlaid with information about the initial cell-type association under normal condition (preferential enrichment in donor PASMC or PAAF) (Figure 4B), revealed that the majority of down-regulated molecules in IPAH-PASMC came from a set of genes and proteins preferentially expressed by PASMC in healthy state (donor-PASMC) (Figure 4B). This reflects the fact that most downregulation in PASMC takes place in contractile machinery elements (Figure 2E) that are characteristic for smooth muscle cell lineage. In contrast, molecules upregulated in IPAH-PASMC contained targets from both PASMC and PAAF enriched data set under healthy cell states (donors). IPAH-regulated molecules in PAAF on the other hand could be traced back to both PASMC- and PAAF-enriched normal cell state (Figure 4B). Gene ontology analysis of IPAH regulated molecules revealed a clear cell-type distinct fingerprint of molecular changes: PASMC harbored expression changes in biological processes related to ECM, while cell cycle and replication processes were enriched in PAAF (Figure 4C). Looking selectively at top20 differentially regulated genes, a similar pattern was observed: genes that were down-regulated in IPAH-PASMC, such as SHROOM3, BLID or CCK, are normally expressed by donor PASMC, while upregulated genes such as chemokines (CCL7, CCL26) are not restricted to smooth muscles (Figure 4D). For comparison, expression of the same genes in PAAF is not significant (Figure 4D). The same was observed for top20 differentially regulated genes in IPAH PAAF – their expression is not significantly changed in PASMC (Figure 4E). However, osteoprotegerin (TNFRSF11B) and matrix metalloproteinase 1 (MMP1), genes with most pronounced expression change in IPAH PAAF, are enriched in PASMC (Figure 4E). These examples show the importance of distinguishing cell-type specific contributions to gene/protein expression levels. It also implies cell-type specific regulatory mechanisms at play.

Figure 4.

Biological function assignment of top20 regulated elements in PASMC and PAAF entailed a broad spectrum of assigned biological functions, ranging from motility, adhesion, hydrolysis, metabolism, growth, differentiation, and ECM components. We wondered if a gene/protein-level evaluation of significantly changed elements for each cell type between IPAH and donors, could give indication of putative phenotypic changes in addition to pathway enrichment approach. We manually curated significantly changed elements in IPAH PASMC (Figure 4F, outermost circle) and IPAH PAAF (Figure 4F, middle circle) and grouped them in ECM, immune system, metabolism, and cell cycle categories (Figure 4F). Elements that additionally had significantly differing regulations between PASMC and PAAF under IPAH conditions were marked bold (Figure 4F, black bordering and bold font). We also provide reference if a particular element was preferentially expressed in PASMC or PAAF (Figure 4F, innermost circle). Striking feature is that both cell types display perturbations in all four major functional categories, the involved genes/proteins are different. For instance, IPAH PASMC upregulate proinflammatory cytokines (CCL7, CCL26) and ECM components claudin 14 (CLDN14) and metalloproteinase ADAMTS18, while IPAH PAAF downregulate anti-inflammatory NF-kB component RelB and laminin 1 (LAMA1) and matrix metallopeptidase 1 (MMP1) ECM components (Figure 4F). Generally, PASMC and PAAF take a divergent path upon vascular remodeling: while IPAH-PASMC acquire synthetic phenotype with up regulation of chemotactic elements, IPAH-PAAF display hallmarks of fast-cycling cells. Based on these results, we wondered if ECM-based phenotypic assays could uncover cell type and state distinct responses, similarly to results obtained by proliferation and metabolic assays.

ECM components centered signaling elicits cell type and cell state distinct phenotypic responses

Changes in ECM composition are a recurring hallmark of pulmonary vascular remodeling ^21–23 recapitulated also in current omics dataset (Figure 4F). Ligand-receptor interaction analysis of our omics dataset indicated that ligand expressing PAAF might act as sender cell type and PASMC as receiver cells through their integrin receptor repertoire expression (Figure 5A). ITGA1 and ITGB1 subunits expressed in PASMC build a heterodimer complex functioning both as laminin and collagen receptor ²⁴. We speculated that change in laminin-to-collagen ratio in ECM would have functional consequence and influence phenotypic behavior of cells in contact with altered ECM. We measured attachment and motility properties of cells with respect to these two substrates. Both PASMC and PAAF displayed striking phenotypic differences depending on surface coating. Whereas cells on collagen-coated surface had spread-out shape with several protrusions (Figure 5B), laminin coating inhibited cell spreading and cells obtained a more symmetric, round-like appearance (Figure 5B). PASMC, both donor and IPAH, attached less efficiently to laminin compared to collagen-coated surfaces (Figure 5C), but only donor-PAAF displayed comparable behavior (Figure 5C). Although not significant, IPAH-PASMC had tendentially higher attachment on laminin-coated surface compared to donor-PASMC (Figure 5C). This trend was again observed in an additional assay measuring gap closure on collagen- or laminin-coated surface (Figure 5D). At early time points, donor PASMC showed widening gap area before final gap closure on laminin coating (Figure 5D). On collagen coating, potential difference was observed at later time point and showed a reverse trend (Figure 5D), probably due to influence of higher proliferative rate of donor PASMC compared to IPAH PASMC (Figure 2H), rather than acute effects of attachment and migration that have more influence in early time points. On both substrates, PAAF revealed a faster gap closure compared to PASMC (Figure 5D), likely in line with their higher proliferative capacity (Figure 2H).

Figure 5.

In addition to changes in collagen and laminin class of ECM components, expression profile of proteoglycans showed both cell type and cell state distinct pattern with IPAH-PASMC characterized by significant enrichment in this ECM family (Figure 2E-G). Among other functions, proteoglycans serve as regulators of the complement system ²⁵. PAAF display enrichment for complement components such as CFD (Figure 1E) and IPAH-PAAF were identified to drive inflammation through classical complement pathway ²⁶. Two biologically active cleavage products, C3a and C5a, are released in this process and both PASMC and PAAF express corresponding receptors (Figure S3). We therefore tested a possible feedback loop between the activated complement system and enhanced proteoglycan expression in IPAH. We stimulated donor and IPAH cells with C3a and C5a and measured expression levels of versican (VCAN, PASMC enriched) and decorin (DCN, PAAF enriched) proteoglycans ⁴. Both VCAN and DCN showed a trend or had significantly higher expression level in IPAH cells (Figure 5E, F). However, only donor PAAF responded to C3a and C5a stimulation and displayed IPAH-state-like expression levels of DCN (Figure 5F). In summary, these results provide a differential phenotypic readout and link the observed ECM expression changes in a cell type and cell state specific manner, with PAAF serving generally as sender type and PASMC as responders, excluding the PAAF autocrine complement response.

IPAH-PAAF modulate phenotypic hallmarks of PASMC

In order to analyze PASMC-PAAF communication in a more comprehensive manner, we took advantage of our published single cell RNA sequencing dataset from human PA ⁴. We extracted PAAF and PASMC clusters and performed intercellular communication network analysis using CellChat. This revealed a skewed PAAF-to-PASMC signaling in IPAH state towards secreted ligands (Figure 6A, Figure S4, S5). We also performed alternative cell-cell interaction analysis using NicheNet and filtered resulting ligand-receptor pairs for secreted signaling using CellChat database annotation. Merged CellChat and NicheNet analysis results revealed a set of functionally-relevant soluble ligands such as IL1A, POSTN, PDGFA/B, FGF2, BMP2/4, CD40LG perturbed in IPAH state (Figure 6B). We therefore asked if cellular crosstalk is able to modulate phenotypic markers using a co-culture model. As reference, we used source-matched donor-PASMC and PAAF that were directly co-cultured for 24 hours followed by sorting and quantitative PCR analysis (Figure 6C). Expression levels of selected cell-type and cell state markers were determined as a readout of phenotypic modulation. Simultaneously, same donor cells were co-culture with IPAH cells from the other cell type, producing a combinatorial readout to assess PASMC-to-PAAF and PAAF-to-PASMC crosstalk (Figure 6C). Results revealed that co-culture of donor-PASMC with IPAH-PAAF lead to partial loss of healthy state PASMC markers (CNN1, RHOA) without affecting the expression of disease-state markers (VCAN, BGN) (Figure 6D). Expression of neither healthy nor disease state PAAF markers was affected in donor-PAAF co-cultured with IPAH-PASMC (Figure 6D). To get a better understanding on donor-PASMC state shift induced by co-culture with IPAH-PAAF, we performed bulk RNA-Seq on collected samples. As the effect of 24h co-culture was mild compared to sample-driven variability, we performed model correction to account for paired PASMC samples in our experimental setup. The resulting gene set enrichment analysis on significantly regulated genes (Data S3) showed that genes upregulated in PASMC by IPAH-PAAF co-culture mediate inflammatory processes and signaling, apoptosis and proliferation (Figure 6E).

Figure 6.

We next went on to identify possible PAAF ligands that could drive the observed phenotypic change in PASMC. As our proteomics approach sampled only proteins from the cell lysate, we additionally performed protein array based screening of 137 unique soluble factors secreted from donor-compared to IPAH-PAAF (Figure 6F). We selected for factors whose expression levels were changed 10% or more in at least two samples and performed clustering of tested samples and factors (Figure 6G). Similarity matrix revealed one major and six minor clusters of factors showing correlated changes (Figure S6). One set of factors correlated with chitinase 3-like 1 upregulation; other clusters with downregulated EMMPRIN, ENA-78, and VEGF secretion, respectively (Figure S6). Three additional clusters involved upregulated secretion of EGF, HGF, and chemerin (Figure S6). We focused our further investigations on the cluster containing hepatocyte growth factor (HGF). Alongside to HGF upregulation, IPAH-PAAF displayed pentraxin-3 (PTX3) downregulation (Figure S6). Intriguingly, HGF and PTX3 seem to have opposite effect on vascular smooth muscle cells, with PTX3 inhibiting neointima formation ²⁷ and HGF promoting pulmonary vascular remodeling ²⁸. Indeed, stimulation of donor PASMC with PTX3 preserved the expression of the contractile markers (Figure 6H), while HGF led to significantly downregulated expression of the same marker genes (Figure 6H). These results indicate that PAAF secretome has a crucial role in maintaining healthy PASMC cell state, but that PAAF can be a driver or phenotypic changes in PASMC.

Discussion

Effective treatment of pulmonary arterial remodeling has been elusive over the years ²⁹. IPAH is characterized by recalcitrant morphological changes and persistent alterations in cellular function. Studies have shown that PASMC and PAAF retain these disease-associated changes in vitro conditions ^30–32. Our current results using transcriptomic and proteomic profiling in combination with phenotypic screening demonstrate that PASMC and PAAF follow different cell fate trajectories and acquire distinct IPAH-associated cell states, while retaining their original cell type designation. This transition from healthy to diseased state is a cell type-specific process and results with functionally meaningful differences between cell types. Furthermore, we identify that maintenance and cell state transition of PASMC is dependent on external cues, such as ECM components and soluble ligands, provided by neighboring PAAF.

The divergent cell fate trajectories can be traced to initial cell lineage differences that are formed and imprinted in a tissue- and anatomical location-specific manner ^1,2. These anatomical-location based distinctions provide fibroblast and smooth muscle cell lineages with phenotypic potential beyond canonically assigned roles. Expression of three selected markers in PAAF, scavenger receptor class A member 5 (SCARA5), complement factor D (CFD) and alcohol dehydrogenase 1C (ADH1C), provide the best example of PAAF pleiotropic roles beyond mere ECM production and structural support. PAAF, being at the outer barrier of vessels, might act as innate immunity sentinels with adventitia being the main site of inflammatory reaction and homing site for inflammatory cells ³³. CFD activates the alternative pathway of the complement triggering chemotaxis and inflammation ³⁴. SCARA5, a ferritin scavenger receptor that mediates cellular iron uptake, limits the extracellular iron concentration ³⁵. Going beyond putative immunomodulatory functions, we additionally reveal PAAF involvement in lipid metabolism at expression profile and phenotypic level. This functionality perhaps reflects broader ontogenetic potential of PDGFRA+ cells as adipocyte precursors ³⁶: PAAF express ADH1C, enzyme involved in metabolism of lipid peroxidation products and hydroxysteroids ³⁷, PPARG, a key transcriptional regulator of lipid metabolism, and APOE, a lipid handling protein (Supplemental Table 1, 2). All three homeostatic functions, complement regulation, iron handling and lipid metabolism, seem to be co-opted and distorted in PAH pathomechanism ^26,38,39.

Simultaneously, cell fate trajectory of PASMC charts a different course. Firstly, we observed that already healthy PASMC display intrinsically higher heterogeneity compared to PAAF. Enrichment in developmental, morphogenesis, migration, and signal transduction biological processes coupled with higher metabolic rate, both oxidative and glycolytic, provides PASMC with functional flexibility and adaptability in response to various injuries. Expression of regulator of G protein signaling 5 (RGS5) and insulin-like growth factor binding protein 5 (IGFBP5) in healthy PASMC, provides two examples of the plasticity potential imprinted in normal PASMC. RGS5, for instance, as a gatekeeper of G protein signaling and signal transduction modulator, influences and fine-tunes the PASMC response to external clue ⁴⁰. IGFBP5 seems to exert a dual role functioning extracellularly as an IGF pathway regulator and intracellularly as SMC transcriptional regulator⁴¹. Under normal condition, expression of contractile machinery coupled with higher metabolic rate, enables PASMC to fulfill their homeostatic role in providing vasoactive function and structural support for endothelial lined lumen. However, one of the defining vascular SMC characteristic and adaptation to vascular insult, including that of IPAH-PASMC cell state, is downregulated expression of the cytoskeletal contractile elements.

Distinctive cell type characteristics of PASMC and PAAF come into play during disease process of pulmonary vascular remodeling, testified by low overlap in commonly perturbed pathways and associated genes/proteins between PASMC and PAAF in IPAH. PAAF and PASMC follow lineage-distinct cell fate trajectories upon remodeling with IPAH-PAAF showing enrichment in DNA replication and repair processes and accompanying hyperproliferative phenotype. IPAH-PASMC in addition to partial loss of contractile markers alter their metabolic profile to support the increased biosynthetic demand. In line with this notion, we have observed significant downregulation of enolase 3, malate dehydrogenase 1B, pyruvate dehydrogenase kinases, and phosphofructokinase, and up-regulation of glutamine-fructose-6-phosphate transaminase 2 in IPAH-PASMC (Supplemental Table 1, 2). This profile indicated potential glucose shunting into hexosamine pathway ^42,43. Hexosamine biosynthetic pathway is a side branch of glycolysis utilized for generation of aminosugars and subsequent protein and lipid glycosylation, including glycosaminoglycans and PG synthesis⁴³. The mostly enriched ECM class next to classical fibrillar collagens, were PG, in line with previously published reports ^4,21,42,44. Due to the complex structure of PG, they can withstand high pressure, modulate cellular phenotype and thus be a necessary part of the adaptation or maladaptation process ⁴⁵.

Interestingly both cell types share a common disturbance in mitochondria characterized by lower IPAH mitochondrial content and hyperpolarized mitochondrial membrane potential. Nevertheless, ultimate phenotypic consequences of such a change still seems to be divergent for each cell type, as implied by distinct stimulated ROS response. Whether mitochondrial disturbance is one the initial cause of the disease or its consequence, as a common point it offers a potential action point to revert cell states from different lineages towards a normal, homeostatic condition, in line with reports on restoring SOD2 expression as a viable therapeutic target ^18,46.

Going beyond cataloguing molecular and functional changes, we address the underlying communication mechanisms responsible for maintenance and transition between identified cell states. We first dissect the role of ECM composition and reveal how a relative change in laminin to collagen ratio, as a consequence of PAAF cell state shift, alters the spreading, attachment, migration, and expansion of PASMC. These processes could induce a more invasive cellular phenotype and potentially promote neointima formation. One could trace this PASMC susceptibility to ECM composition change to their enrichment in integrin family of receptors expression making them uniquely primed to sense and respond to external cues. However, the same mechanism renders cells susceptible to phenotypic change induced simply by extended vitro culturing, testified by loss of phenotypic characteristics upon prolonged culturing. We also uncover an autocrine loop centered on PAAF and showing how activated complement pathway can drive expression of PAAF-enriched ECM components. We extend our investigations more broadly and uncover a general shift in IPAH-PAAF secretome being sufficient to induce significant phenotypic changes in normal PASMC. Among several deregulated ligands, we identify downregulated expression of PTX3, a protective factor, and upregulation of HGF, pathogenic factor, as PAAF-sourced regulators of PASMC cell state. PTX3 counteracts pro-proliferative action of fibroblast growth factor and its overexpression attenuates neointima formation in systemic arteries²⁷. Conversely, upregulation of HGF/cMet pathway in endothelial cells upon Sox17 loss was associated with PAH promotion ²⁸. We would therefore argue that potential exploratory clinical trial for pulmonary vascular disease could entail a combination therapy consisting of recombinant PTX3 and HGF blocking antibody or c-Met inhibitor. However, the identified factors likely represent just a fraction of ligands expressed by other cell types present in the pulmonary artery compartment. The influence of endothelial and inflammatory cells on PASMC and PAAF phenotype requires further investigation.

Nevertheless, our results provide detailed understanding of functional changes in relation to underlying molecular expression differences upon pulmonary vascular remodeling. Furthermore, we identify novel targets for cell-type specific correction of diseased cell states and potential reverse remodeling strategies.

Acknowledgements

We appreciate excellent technical assistance from Julia Kohlbacher, Hans Peter Ziegler and Elisabeth Blanz. Schematic figure panels for this article are created with BioRender.com. This work was supported by Cardio-Pulmonary Institute EXC 2026 390649896 (GK) and Austrian Science Fund (FWF) grant I 4651-B (SC). For open access purposes, the author has applied a CC BY public copyright license to any author accepted manuscript version arising from this submission.

Author contributions

Conceptualization: SC, HTP, VJP, GK

Methodology: HTP, JW, MB, AW, IM, AM, RL, ASO, NB, SLH

Investigation: SC, VB, JW, EB, AW, IM, BMN, ME, MW

Visualization: HTP, SM, FV Supervision: VJP, GK

Writing—original draft: SC, HTP

Writing—review & editing: KH, AO, KS, LMM, VJP, GK

Declaration of interests

The authors declare no competing interests.

Data and materials availability

Genome-wide expression profiling and RNA-Seq data is deposited to the National Center of Biotechnology Information Gene Expression Omnibus database (accession numbers GSE255669, GSE262069). Human tissue samples can be provided pending ethical approval and completed material transfer agreement.

Experimental Procedures

Human tissue samples

Human lung samples from IPAH patients and downsized non-transplanted donor lungs, serving as healthy control, were obtained from Division of Thoracic Surgery, Medical University of Vienna, Austria. The protocol and tissue usage was approved by the local ethics committee (976/2010; 1417/2022) and patient consent was obtained before lung transplantation. Small size, resistance pulmonary arteries (diameter < 0.5 cm) were dissected from the lungs and used for smooth muscle cell and fibroblast cell isolation or flash-frozen in liquid nitrogen. Briefly, adventitial layer was pulled off the PA and used for outgrowth of PAAF. Remaining PA was cut open, endothelial layer scraped off and remaining media mechanically cut into small size pieces and used for smooth muscle cell outgrowth. Cells were grown in full media (VascuLife SMC or FibroLife S2, LifeLine Cell Technology) and used up to passage four. Respective basal media without added supplements and serum was used for starvation. Expression profiling (transcriptomic and proteomic) and functional cellular assays were performed on cells isolated from the same pulmonary artery. Histological staining of formalin-fixed paraffin-embedded human lung samples (n=6) was done on same control/patient cases that were used for in vitro studies. Frozen PA used for glysoaminoglycan measurement (n=8 for each condition) consisted of an independent control/patient cohort. Clinical data is given in Table S1.

Genome-wide expression profiling

Briefly, whole genome expression profiling was performed on total RNA isolated using RNeasy Mini kit (Peqlab) from PAAF and PASMC in first passage (4 healthy donors and 4 IPAHs). Purified total RNA was amplified and Cy3-labeled using the LIRAK kit (Agilent) following the kit instructions. Per reaction, 200 ng of total RNA was used. The SureTag DNA labelling kit (Agilent, Waldbronn, Germany) was used to Cy5- and Cy3-label the samples and subsequently hybridized to 60-mer oligonucleotide spotted microarray slides (Agilent Human G3 8×60k, Design ID 072363). The following hybridization, washing, and drying steps were performed following the Agilent hybridisation protocol. Thereafter, the slides were scanned at 2 µm/pixel resolution using the InnoScan 900 (Innopsys, Carbonne, France). Image analysis was performed with Mapix 6.5.0 software, and calculated values for all spots were saved as GenePix results files. Stored data were evaluated using the R software (R Development Core Team 2007) and the limma package ⁴⁷ from BioConductor ⁴⁸. Log mean spot signals were taken for further analysis. Data was background corrected using the NormExp procedure on the negative control spots and quantile-normalized before averaging ⁴⁹.

Genes were ranked for differential expression using a moderated t-statistic ⁵⁰. Results were uploaded to the National Center of Biotechnology Information Gene Expression Omnibus database (accession number GSE255669). List of differentially expressed genes is given in Data S1.

Proteomic Analysis

Protein samples from isolated PASMC and PAAF were digested overnight at 37°C with Trypsin/LysC (Promega). Proteolytic digestion was quenched with 1% formic acid. The dried peptides were dissolved in reconstitution buffer (2% acetonitrile + 0.1% formic acid) and equimolar amounts (based on starting total protein concentration) of sample injected into the mass spectrometry instrument. Experiments were performed on the Orbitrap Fusion Tribrid mass spectrometer (Thermo Scientific) coupled with ACQUITY M-Class ultra-performance liquid chromatography (UPLC, Waters Corporation). A flow rate of 450 nL/min was used for this Liquid Chromatography/Mass Spectrometry experiment, where mobile phase A was 0.2% formic acid in water and mobile phase B was 0.2% formic acid in acetonitrile. Analytical columns were pulled using fused silica (I.D. 100 microns) and packed with Magic 1.8-micron 120Å UChrom C18 stationary phase (nanoLCMS Solutions) to a length of ∼25 cm. Peptides were directly injected onto the analytical column using a gradient (2-45% B, followed by a high-B wash) of 80 minutes. The MS was operated in data-dependent fashion using CID (collision induced dissociation) for generating MS/MS spectra, collected in the ion trap. Raw data was processed using Byonic v3.2.0 (ProteinMetrics) to infer protein isoforms using the Uniprot Homo sapiens database. Proteolysis with Trypsin/LysC was assumed to be semi-specific allowing for N-ragged cleavage with up to 2 missed cleavage sites. Precursor mass accuracies were held within 12 ppm and 0.4 Da for MS/MS fragments. Proteins were held to a false discovery rate (FDR) of 1% or lower, using standard target-decoy approaches and only the proteins with >3 spectral counts were selected for further data processing (keratins and KAPs were removed at this stage). List of detected proteins is given in Data S2.

Gene/Protein and pathway-centered functional data exploration and graphical representation

Gene Set Enrichment Analysis (GSEA) and Over-Representation analysis (ORA) of transcriptomics and proteomics data was performed by using GO dataset ^50,51 or BioPlanet ^50,52 and applying the indicated significance and LFC cutoffs. The STRING (Version 11.0) Protein-Protein Interaction Networks/Functional Enrichment Analysis database was used to retrieve data on gene/protein interactions and examine the connections between genes. The gene-list enrichment tool Enrichr ^50,53 was utilized to compare functional results. The open-source software platform Cytoscape (Version 3.7.2) was used for visualizing complex networks and representation of pathway terms with root nodes. Furthermore, RStudio (Version 1.1.383) and the package for circular visualization ⁵⁴ and ggplot package ⁵⁵ were used for creating the circular plots.

Multiplex immunofluorescent and RNA in situ hybridization imaging

Formalin-fixed, paraffin-embedded lung sections were deparaffinized, rehydrated and subjected to heat induced antigen retrieval (pH6). Multiplex immunofluorescence staining was done using Opal kit (Akoya). Successive rounds of primary antibody, detection reagent (Opal Polymer HRP, Akoya or Immpress HRP Polymer, Vector Labs) and fluorescent signal development (Opal dyes, Akoya or CF tyramide dye, Biotium), followed by antibody removal were performed according to manufacturer’s instructions. List of antibodies is given in Table S2. DAPI (Thermo Scientific) was used as a nuclear counterstain at 2.5 µg/mL final concentration. In the case of single molecule RNA in situ hybridization, slides were first against anti-PDGFRA antibody and processed further for RNA in situ hybridization according to manufacturer’s instructions using RNAscope Multiplex Fluorescent V2 assay (ACD). Hybridization probes against human IGFBP5, SCARA5 and CFD were obtained from ACD. Following hybridization signal detection, slides were incubated with primary labelled antibody against alpha smooth muscle actin and counterstained with DAPI. Slides were imaged using SP8 fluorescence confocal microscope (Leica) and apochromatic glycerol immersion objectives 40x (1.25 numerical aperture).

SeaHorse metabolic flux measurements

Isolated PASMC and PAAF (20’000/well) were plated directly into Seahorse XF96 cell culture microplates (Agilent Technologies) and allowed to reach confluency. After 96 and 120 hours the medium was replaced by Standard Oxidation Stress Test assay medium (10 mM Glucose, 1mM Pyruvate, 2mM Glutamine) and plates were incubated in a non-CO2 incubator for 60 minutes prior to the assay. Oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) were measured in a Seahorse XFe96 extracellular flux analyzer (Agilent) using XF Substrate Oxidation Stress Test Kits (XF Long Chain Fatty Acid Oxidation Stress Test Kit, XF Glucose/Pyruvate Oxidation Stress Test Kit). First, basal respiration was established, followed by injection of the relevant pathway inhibitor (4 µM etoxomir, 2 µM UK5099). The response to the inhibitor was monitored, then ATP production was blocked by oligomycin (1.5 µM), followed by the dissipation of proton gradients with FCCP (1 µM). The electron transport chain was inhibited by rotenone/antimycin (0.5 µM) to reveal the non-mitochondrial respiration. OCR/ECAR values are representing the ratio of Mitochondrial Respiration/Glycolysis derived from the very first rate measurements of the Substrate Oxidation Stress Test. The acute response to the inhibitor represents alterations of basal respiration as OCR response to pathway inhibitors (etoxomir, UK5099). The glycolytic capacity was calculated as the difference of maximum ECAR measurement after oligomycin injection and the last ECAR measurement right before oligomycin injection. Non-mitochondrial respiration was defined as the minimum OCR measurement after rotenone/antimycin injection. Proton leak was calculated as minimum OCR measurement after oligomycin injection minus non-mitochondrial respiration. Measurements were performed as septuplets, rate measurement values from each well were normalized to their corresponding protein amounts. Data were analyzed using the Wave Software (Agilent Technologies, Version-2.6.1).

Proteoglycan/Glycosaminoglycan detection

Combination of Alcian blue and Verhoeff’s staining of formalin-fixed, paraffin-embedded lung sections was used for localization of glycosaminoglycans and elastic fibers, respectively. Glycosaminoglycans were quantified in isolated cryo-preserved PA of donor and IPAH patients as described previously ⁵⁶. Briefly, isolated PA were weighed and digested overnight at 60°C in a phosphate-EDTA buffer papain solution. Resulting lysate was mixed with 1,9-dimethylmethylene blue assay and the absorbance of the resulting complex measured spectrophotometrically at 525 nm. Chondroitin sulphate was used for standard curve preparation.

Cellular turnover measurements

Proliferation of PASMC and PAAF (6 healthy donors and 6 IPAHs) was determined by [³H]thymidine incorporation assay (BIOTREND Chemikalien, Germany). 7500 cells per well were seeded in a 96 well plate. The cells were starved in basal medium (without supplement) overnight and next day they were stimulated with 2% FCS or maintained in media without supplements as controls. Following 24h [³H]thymidine incorporation measurement was evaluated in each well. Alternatively, cell turnover in different passages was done by plating 100’000 cells in a flask and allowed to grow for 72 hours. Cells were detached by trypsination and cell number were determined using Neubauer cell counting chamber.

Mitochondrial content, membrane potential and cellular ROS

Cellular reactive oxygen species (ROS) production was based on fluorescent imaging of the ROS indicator CellRox Deep Red (Thermo Scientific). Cells were loaded with 5µM of CellROX dye for 1 hour at 37°C, followed by confocal fluorescence imaging (LSM 510 Meta with a Plan Neofluar 40x/1.3 oil-immersion objective, Carl Zeiss) before and after stimulation with 50 ng/mL PDGF-BB (R&D Systems) for 40 minutes at 37°C. Mitochondrial membrane potential was measured in quench/de-quench mode following staining with 1 µM TMRM dye (Sigma Aldrich). Image analysis was performed in ImageJ (NIH). A binary image of the TMRM signal was created, to measure the cell area covered by mitochondria. Mitochondrial density was then determined by dividing the cell area with the area covered with mitochondria. Mitochondrial membrane potential was calculated by normalizing the TMRM signal to the mitochondrial density.

RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA from PAAF and PASMC was extracted using peqGOLD Total RNA Kit (Peqlab) or innuPREP RNA Mini Kit (IST Innuscreen). The cDNA was reverse transcribed using a qScript cDNA Synthesis kit (QuantaBio). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed using the S’Green qPCR kit (Biozym) on LightCycler 480 (Roche Applied Science). The threshold cycle difference (ΔCT) was calculated using the equation ΔCT = ([CT(B2M)+CT(HMBS)]/2) – CT(Gene of Interest)]. Beta-2-microglobulin (B2M) and hydroxymethylbilane synthase (HMBS) were used as reference genes. The primer list is given in Table S3.

Adhesion and gap closure assays

Adhesion assay was performed on Biocoat collagen I and laminin coated 24-well plates (Corning). Briefly, 50’000 cells were allowed to attach for 40 min at 37°C, followed by PBS wash, formalin fixation and staining with crystal violet. Images were transformed into 8-bit image type and the adhered cells were counted with the plugin in cell counter of the Fiji-Image J. Gap closure assay was perfomed in two chamber silicon inserts with a defined cell-free gap were used (Ibidi). Each chamber seeded with 15’000 cells that were allowed to grow till confluency. Experiment was started by removing the insert and image of the cell-free gap taken (t=0). Subsequent images of the gap were taken 4, 8, 24 and 48 hours after removal of the insert. Analysis of wound area has been performed with Fiji-ImageJ. Briefly, the wound area of each picture was first identified with the MRI Wound healing tool and the area enclosed was then measured.

Cell stimulations and co-culture experiment

PASMC were cultured in 24well or 6well plates up to 90% confluence and serum starved for 24 hours. Cells were stimulated for 24 hours with ligands C3a (R&D Systems, 100 ng/mL), C5a (R&D Systems, 100 ng/mL), pentraxin-3 (R&D Systems, 5 μg/mL), hepatocyte growth factor (Peprotech, 25 ng/mL) in starvation medium. Following stimulation, cells were washed with PBS and lysed for subsequent RNA isolation.

For PAAF-PASMC co-culture experiment, cells were fluorescently labelled with either CellTracker Green CMFDA (Thermo Scientific) or CellTracker Red CMTPX (Thermo Scientific) for 30 min at 37°C. Labelled cells were trypsyized and plated in a 1:1 ratio (PAAF:PASMC) in a 6well plate in fully supplemented smooth muscle cell medium for 24 hours. Next day, cells were detached and single cell suspension processed for two-way fluorescence activated cell sorting. CMFDA and CMTPX viable cells (Sytox Blue dead cell stain negative cells) were sorted using FACSAria IIIu (BD Biosciences) at the flow cytometry core facility (Center for Medical Research, Medical University of Graz). Collected cells were centrifuged and processed for RNA isolation. One part of RNA was used for qRT-PCR, while other RNA was used for bulk RNA sequencing analysis. Briefly, libraries were prepared using NEB Next® Ultra II RNA Library Prep Kit for Illumina (New England Biolabs) with rRNA Depletion kit v2 (Human/Mouse/Rat) (New England Biolabs) and sequenced upon library pooling on NovaSeq 6000 (Illumina). Results were uploaded to the National Center of Biotechnology Information Gene Expression Omnibus database (accession number GSE262069). List of differentially expressed genes is given in Data S3.

Cell-cell communication analysis

Ligand-receptor expression analysis was performed on previously published single cell RNA sequencing dataset of human pulmonary arteries from donors and PAH patients GSE210248 ⁴. For the analysis, the data matrix was transformed into a Seurat object and was processed with the Seurat package (Version 4.3.0) ⁵⁷ in the R studio environment (Version 4.1.2). The 10x files were converted into a Seurat object with CreateSeuratobject. Doublets were filtered out with the scDblFinder tool. Furthermore, all cells with more than 5% mitochondrial DNA (mtDNA) and cells with under 200 and over 2500 RNA counts were filtered out.

Final quality control step included filtering out outliers defined as cells with an RNA-count outside the range of 5 MADs. Data normalization was performed with the NormalizeData function from Seurat, followed by data scaling (ScaleData). The FindVariableFeatures Seurat-function detected the variable features in the data set (selection method=vst, nfeatures=2000). A principal component analysis was run on the normalized and scaled data. The number of PC was set after plotting an elbow plot which showed an elbow between 20 and 30 clusters. The following steps included data integration with harmony(Version 0.1.1) and Uniform manifold Approximation (UMAP) depiction of the data (RunUMAP). The data was clustered with the help of the FindNeighbours- and FindClusters function. The cluster annotation to cell type was done manually based on gene enrichment of canonical cell type marker genes.

For the cell-cell-communication analysis the clustered datasets were normalized so that the number of cells per cell type in each group (PAH or donor) were comparable. Two independent analysis methods were used. Nichenet analysis: A comparative Nichenet analysis was performed according to the Nichnet method (https://github.com/saeyslab/nichenetr, Version 2.0.4) and as described in the corresponding paper ⁵⁸. Donor and PAH niches were defined and smooth muscle cells and fibroblasts were each set as sender and receiver cell type. The analysis was performed twice to obtain results for communication in both directions. CellChat analysis: A decriptive analysis of cell-cell-communication for donor and PAH samples was performed according to the CellChat tool ((https://github.com/sqjin/CellChat, Version 1.6.1) ⁵⁹. The Seurat object was subsetted, and the analysis was performed separately for each subset.

Secretome analysis

Proteome Profiler Human XL Cytokine Array Kit (R&D Systems) and Proteome Profiler Human Adipokine Array Kit (R&D Sytems) were used to measure factors secreted by donor and IPAH PAAF. Briefly, PAAF were grown till 80-90% confluence and after 3-4 days medium was collected, centrifuged at 4°C, 400g for 5 min, aliquoted and stored at −80°C till measurement. Equal volume of cell culture supernatant was used in Proteome Profiler Array Kits and assay performed according to manufacturer’s instructions. Signal intensity values for each spot were normalized to the mean value of the positive control spots for each membrane. Normalized intensity values were then used to calculate pairwise (IPAH-to-donor) relative expressions for each factor above the detection limit. In total, 4 such pairwise comparisons were done. We took into further consideration only those factors whose relative normalized expression levels changed by more than 10% in at least two pairwise comparisons. Thus obtained list of ligands with corresponding normalized relative expressions were hierarchically clustered using Clustergrammer online tool ⁶⁰ to provide a heatmap and correlation matrix visualization.

Statistical analysis

Graphs and statistical calculations were performed using GraphPad Prism (Version 8.0, GraphPad Software) or R (v3.5.3, packages stringr, data.table, readxl, openxlsx, MetaboAnalystR, ggplot2, colorspace, circlize) and Tibco Spotfire (v11.1.0). Probability values of p<0.05 were considered statistically significant. Differences between groups were investigated using Student’s t test, Mann-Whitney or Friedman test. 2-way or 3-way ANOVA was used for interaction analysis. For orthogonal PLS-DA (OPLS-DA) data was filtered to retain only features with less than 30% missing data (proteomics 1552 genes, transcriptomics 23541 genes). OPLS-DA was performed centered and scaled to unit variance with a standard 7-fold cross validation for the classification factor. Model stability was additionally verified with 1000 random label permutations and models with Q2>50% were considered significant.