Identification of cellular heterogeneity and key signaling pathways associated with vascular remodeling and calcification in young and old primate aortas based on single-cell analysis

Single-cell profile of vascular wall cells

The arterial vascular wall is mainly composed of ECs, SMCs, AFs and ECM [6]. All components of the vessel wall interact with each other to maintain the physiological homeostasis of blood vessels. In the healthy artery, the vessel wall can change its structure and function to adapt extravascular environment or hemodynamic demands [3]. However, with the aging and function losing of vasculature, these adaptive changes do not return to baseline level but instead initiate pathological vascular remodeling, which further leads to the development of cardiovascular disease [3]. These vascular remodeling was obviously observed in senile aortas, represented as increased wall thickness, fibrous cap formation, arterial calcification and fragmentation of the elastic lamina [20].

In order to analyze the cellular heterogeneity and the changes of cellular communication between different cell types of aged monkey vasculature. We analyzed a published scRNA-seq dataset of aortic arches from four young (age, 4–6 years) and four old (age, 18–21 years) male cynomolgus monkeys (GEO accession nr. GSE117715). As we focused on cells forming the wall of blood vessel, we only extracted ECs, SMCs and AFs from cell types that were identified in the original studies. The cells of ECs, SMCs and AFs were merged into an integrated Seurat object. After data preprocessing and dimensional reduction, ECs, SMCs and AFs were separately displayed in the UMAP plot (Figure 1A). The traditional marker genes of ECs (PECAM1), SMCs (ACTA2) and AFs (LUM) were specifically expressed in ECs, SMCs and AFs respectively (Figure 1B), which indicated that the definition of these cell types was correct. To further verify it, we identified the top 5 highly expressed genes in each cell type, and found that those genes were all cell type specific marker genes (Figure 1C). These results further verified the correctness of cell types.

Figure 1. Transcriptomic profile of aortic arches from young and old monkeys. (A) UMAP plots of the scRNA-seq datasets displayed by group and cell type. (B) The expression of classic markers in each cell type. (C) Top 5 marker genes in each cell type. Abbreviation: UMAP: Uniform Manifold Approximation and Projection.

EC1 may be a subpopulation related to arterial remodeling and vascular calcification

ECs, a single cell layer in the internal surface of the vasculature, play a crucial role in sensing hemodynamic changes and maintaining vascular homeostasis. To figure out the effect of aging on ECs of monkey vasculature, we performed differentially expressed genes (DEGs) analysis between old and young aortic arteries, and found 162 upregulated genes and 223 downregulated genes (Figure 2A). Most of these upregulated genes were associated with extracellular matrix formation, such as COL3A1, FN1 and BGN. Moreover, genes related to lipid metabolism and transport, including APOD and APOE, were highly expressed in old aortic arteries. Functional enrichment analysis showed that these upregulated genes were mainly enriched in biological processes of response to oxidative stress, wound healing and protein folding (Figure 2B). These results indicated that the function and structure of ECs may have changed during aging.

Figure 2. Characteristics of ECs and their subpopulations. (A) Volcano plot of DEG (average log fold change > 0.25, P value < 0.01) between young and old group. (B) GOBP analysis of upregulated DEG. (C) UMAP plots of EC subpopulations displayed by cell type. (D) Heatmap showing the top 10 genes of each EC subpopulation (left) and GOBP enriched by upregulated DEG in each subpopulation (right). Abbreviations: DEG: differentially expressed genes; GOBP: gene ontology biological process.

To further analyze the structural and functional heterogeneity of ECs, we performed unsupervised clustering to ECs alone and found out four distinct subpopulations (Figure 2C). We then identified the top 10 highly expressed genes of each subpopulation. Each of the four ECs subsets expressed distinct marker genes, suggesting the presence of discrete subpopulations (Figure 2D). We observed that EC1 highly expressed genes related to extracellular matrix organization and ossification, such as FN1, COL3A1, BMP4 and SPARC (Supplementary Figure 1). The biological processes of EC1 were mainly enriched in blood coagulation, wound healing and response to transforming growth factor beta, which are involved in ECM remodeling and vascular fibrosis (Figure 2D). While, the biological processes of EC2 subpopulation were related to negative regulation of MAPK cascade, negative regulation of proteolysis and response to reactive oxygen species. The upregulated genes in EC3 were enriched in biological processes of response to glucocorticoid/corticosteroid, extrinsic apoptotic signaling pathway and aging (Figure 2D). EC4 was enriched for glycerolipid biosynthetic process and glycerophospholipid metabolic process. These results showed that ECs had structural and functional heterogeneity and EC1 may be a subpopulation related to arterial remodeling and vascular calcification.

The structural and functional heterogeneity of SMCs

The SMCs reside in the medial layer of major blood vessels. They play a vital role in vascular remodeling and calcification [9]. To figure out the effect of aging on SMCs of monkey vasculature, we performed DEG analysis between old and young aortic arteries, and found 573 upregulated genes and 816 downregulated genes (Figure 3A). These upregulated genes in old aortic arteries were mainly enriched in biological processes of actin filament organization, cell−substrate adhesion and cytoplasmic translation (Figure 3B).

Figure 3. Characteristics of SMCs and their subpopulations. (A) Volcano plot of DEG (average log fold change > 0.25, P value < 0.01) between young and old group. (B) GOBP analysis of upregulated DEG. (C) UMAP plots of SMC subpopulations displayed by cell type. (D) Heatmap showing the top 10 genes of each SMC subpopulation (left) and GOBP enriched by upregulated DEG in each subpopulation (right). Abbreviations: DEG: differentially expressed genes; GOBP: gene ontology biological process.

To further analyze the structural and functional heterogeneity of SMCs, we performed unsupervised clustering to SMCs and identified two different subpopulations (Figure 3C). We then identified the top 10 highly expressed genes, and performed functional enrichment analysis for these upregulated genes of each subpopulation (Figure 3D). SMC1 had a contractile phenotype and highly expressed traditional marker genes of SMCs, such as ACTA2 and CNN1 (Supplementary Figure 2). The biological processes enriched by the upregulated genes in SMC1 were involved in muscle tissue development and response to mechanical stimulus, which are required for maintaining normal function of arteries (Figure 3D). While, SMC2 seemed to have a synthetic VSMC phenotype, with losing contractility markers and expressing proteins involved in extracellular matrix organization (Supplementary Figure 2). Besides, SMC2 highly expressed ossification related genes (EGFR, CHRDL1) and aging-related genes (SERPINF1, APOE), and the biological processes of these upregulated genes were enriched in ossification and aging.

Deciphering cell-cell communication in all cell subpopulations of vascular wall

Intercellular communication between different cell types of vascular wall is important for aging process and aging related arterial remodeling [23]. In order to distinguish the changes of cellular communication between different cell subpopulations of aged aortas, using CellChat, we evaluated cell-cell communication patterns in seven subpopulations of aortic wall cells, including four EC subpopulations, two SMC subpopulations and one AF subpopulation. Results of CellChat analysis revealed 2,085 total ligand-receptor interactions in the young group and 2,246 interactions in the old group (Figure 4A). But the total interaction strength of the old group was moderately lower than that of the young group (Figure 4A).

Figure 4. Overview of intercellular communication between young and old group. (A) Total interaction number and strength in young and old group. (B) Heatmap of differential interaction strength in young and old group. (C) Scatter plot of incoming and outgoing interaction strength of each cell population in young and old group. (D) Overall information flow of each signaling pathway in young and old group.

To figure out the signal changes of each subpopulation, we compared the outgoing and incoming signals of each cell population in young and old aortas. We found that the signals from EC1 to SMC1 and SMC2 were increased significantly, while the signals received by AFs were all slightly decreased (Figure 4B). Compared with young group, the incoming signals of SMC1, SMC2 and EC1 were all strengthened, and that of EC2 were weakened obviously (Figure 4C). To identify specific signaling pathways that changed in the aged aortas, we compared the overall information flow for each signaling pathway between young and old group. Results showed that a number of signaling pathways in old group were stronger than that in young group, such as VEGF, BMP, NOTCH and EGF. Moreover, there were four signaling pathways (BSP, CHEMERIN, IGF, WNT) unique to old group (Figure 4D).

The altered signals in subpopulations of aged aortas

To further analyze the alteration of signals in each cell subpopulation, we compared the signal strength of each signaling pathway between young and old group (Figure 5A). Most of the overall signaling pathways in AFs were weakened in old group, such as COLLAGEN and LAMININ (Supplementary Figure 3). However, we found that the VEGF signaling pathway was increased in aged aortas and it mainly sent from AFs to ECs (Supplementary Figure 3). VEGF was a vascular endothelial-derived growth factor, which can regulate multiple biological functions, including vascular permeability, angiogenesis and endochondral ossification [24]. In contrast, a number of signaling pathways were increased in EC1 subpopulation of aged aortas, such as PECAM1, ESAM and CDH5. PECAM1 is closely correlative with cell migration, proliferation, apoptosis, signal transduction and cellular immunity [25]. ESAM is involved in leukocyte extravasation and VEGF-induced vascular permeability [26]. Interestingly, these signaling pathways have similar changes in the aged aortas. The cell-cell communications between EC2 and EC3 were weakened, but the communications between EC1 and EC3 were strengthened in the aged aortas (Supplementary Figure 4). These results indicated that the cellular communications of ECs have changed during aging.

Figure 5. Overall signaling patterns and inferred signaling pathways. (A) Overall signal strength of each signaling pathway in each cell population. (B) Circle plot of PDGF signaling network in young and old group. (C) Circle plot of NOTCH signaling network in young and old group. (D) Circle plot of FN1 signaling network in young and old group. (E) Circle plot of BMP signaling network in young and old group.

Similarly, a number of signaling pathways were increased in SMC subpopulations. For example, the ANGPT, THBS and CSPG4 were overall highly expressed signaling pathways in SMC1 (Supplementary Figure 5). ANGPT is a second endothelial cell specific ligand–receptor signaling system, which is involved in postnatal angiogenesis, vessel remodeling, vascular permeability and inflammation to maintain vascular homoeostasis in adult physiology [27]. Moreover, THBS are evolutionarily conserved, calcium-binding glycoproteins that contribute to wound healing and angiogenesis, vessel wall biology, connective tissue organization, and synaptogenesis [28]. Recent analysis found that THBS signaling pathway was associated with abdominal aortic aneurysm [22], which is a kind of aged related cardiovascular disease. These results indicated that the cellular communications of SMCs have changed during aging.

The altered signals between EC1 and SMC

Adequate cooperation between ECs and SMCs is required to maintain vascular normal function. Disturbances in this balance predispose the vessel to initiate pathological change. As EC1 is a subpopulation related to vascular remodeling and calcification, we analyzed the changed signals that sent from EC1 to SMCs. PDGF is a platelet-derived growth factor, which plays a prominent role in the pathogenesis of cardiovascular disease [29]. We observed that increased PDGF signal was mainly sent from EC1 to SMC1 and AFs (Figure 5B). The NOTCH signaling pathway, which is involved in cell fate determination and tissue homeostasis [30], was also strengthened in the old group (Figure 5C). As we mention above, the aged aortas have phenotypes of vascular remodeling and calcification. We then examined the signaling pathway of FN1, which is associated with vascular remodeling and calcification. We found that the signaling pathway of FN1 was significantly increased from EC1 to SMC1 in the old group (Figure 5D). Bone morphogenetic proteins (BMPs) are a group of ligand proteins that involved in regulating vascular calcification and bone formation [31]. In our study, the BMP signaling pathway and its ligand-receptor pairs were all increased from EC1 to SMC1 in the old group (Figure 5E), implying its key role in vascular aging. These results indicated that these changed signals between EC1 and SMCs may relate to the arterial remodeling and calcification of age aortas.

SMCs are the major cells in the vascular wall and play an important role in arterial remodeling and calcification. To further analyze the signaling changes in SMCs, we compared the incoming and outcoming signals between young and old groups. Results showed that FN1, LAMININ and COLLAGEN were the most changed signals in SMCs (Figure 6A). We noticed that EGF signaling pathway was specifically received by SMCs in the old group. This signal was mainly sent from EC3 to EC2 in the young group, but in the old group, it was changed and mainly sent to EC3, SMC1 and SMC2 (Figure 6B). Besides, the ligand-receptor pairs of EGF signaling pathway were all strengthened from ECs to SMCs in the old group (Figure 6C). Among these ligand-receptor pairs, BTC-EGFR ligand-receptor was increased significantly in the old group. EC1 and EC3 were the main output subpopulations of BTC-EGFR signaling, while the SMCs and AFs were the receiver of that (Figure 6D). To further verify it, we analyzed the gene expression of ligand-receptors of EGF signaling pathway in each subpopulation. We observed that the genes of ligands were highly expressed in ECs, while the receptors were mainly expressed in AFs and SMCs (Figure 6E). Notably, the expression of EGFR, an important receptor of EGF signaling pathway, was significantly higher in EC1, SMC1 and SMC2 of the old group than the young group (Figure 6E). Besides, the average expression of EGFR in old group was also higher than that in young group (Figure 6F). These results indicated that EGF signaling pathway may play an important role in vascular aging.

Figure 6. EGF signaling in young and old group. (A) Signaling changes of SMC1 and SMC2. (B) Circle plot of EGF signaling network in young and old group. (C) The communication probability of all the significant ligand-receptor pairs that contributed to EGF signaling. (D) Circle plot of BTC-EGFR ligand-receptor pairs in young and old group. (E) The expression of genes of EGF signaling pathway in each subpopulation and group. (F) The expression of EGFR in young and old group.

Data acquisition and preprocessing

We downloaded the integrated dataset provided in the NCBI GEO data repository (GEO accession no. GSE117715), which is a merged dataset containing expression matrix of all samples from eight young (4–6 years old) and eight old (18–21 years old) monkeys (Macaca fascicularis) [20]. To avoid the effect of sex-dependent factors on the composition of vascular cells, we only extracted the aortic dataset of 4 young and 4 old male monkeys. Since we focused on studying the changes in the structure and function of blood vessel walls with aging, only vascular component cells (ECs, SMCs and AFs) were selected for subsequent analysis. The R package Seurat (version 4.1.0) was used for subsequent analysis in R (version 4.1.2) environment [54]. To exclude poor-quality cells, cells that gene count per cell <1000 and >6500, and total counts >300000 were discarded. After quality control, data normalization and scaling were performed with Seurat package. The top 3000 variable genes for all samples were used for principal component analysis (PCA).

Clustering of scRNA-seq data

The cluster of ECs and SMCs were performed with Seurat respectively. First, the extracted data of ECs or SMCs was normalized and scaled according to the Seurat manual. Similarly, the top 3000 variable genes were used for PCA, and uniform manifold approximation and projection (UMAP) was used for nonlinear dimensional reduction. The top 20 principal components were then used to cluster cells at the resolution of 0.3 for ECs and SMCs.

Differential gene analysis and functional annotation

Cell type-specific and subcluster-specific markers of ECs or SMCs were analyzed using the “FindAllMarkers” function of Seurat with the Wilcox test. Differentially expressed genes (DEGs) between indicated cell types or sub-clusters were found with “FindMarkers” function, and genes with adjusted p values < 0.01 and |logFC| > 0.25 were considered to be DEGs [55]. Gene set enrichment analysis of the DEGs was performed using the R package clusterProfiler (version 4.2.2). GO terms with corrected P < 0.05 were considered significantly enriched.

Analysis of cell–cell communication

Cell–cell communication between each cell subpopulation of vascular wall was performed using CellChat (version 1.1.3) with default settings [21]. Briefly, cells of ECs, SMCs and AFs were merged into an integrated Seurat object, and then subjected to CellChat pipeline to analyze. The ligand-receptor interaction database CellChatDB in human was used for following analysis, and the “computeCommunProb” function from CellChat was used to compute the communication probability. The signaling pathways were visualized with Hierarchy plot, Circle plot or Chord diagram, and the communication patterns of target cells were visualized with dot plots.

Statistics

Statistical analysis was performed with default method in Seurat package. As we previously described [56–59], the P values for gene expression analysis in different groups were estimated using Wilcoxon Rank Sum test. As for cell communication analysis, statistical analysis was also performed with default method in CellChat package, and default method of “triMean” was used for computing the average gene expression per cell group. P-value < 0.05 was used for determining significant interaction.

Data availability statement

The data that support the findings of this study are available in GEO database at (https://www.ncbi.nlm.nih.gov/geo/), reference number (GSE117715). Further inquiries can be directed to the corresponding authors.