Single-Cell Gene-Regulatory Networks of Advanced Symptomatic Atherosclerosis

Abstract

Background: While our understanding of the single-cell gene expression patterns underlying the transformation of vascular cell types during the progression of atherosclerosis is rapidly improving, the clinical and pathophysiological relevance of these changes remains poorly understood.

Methods: Single-cell RNA sequencing data generated with SmartSeq2 (≈8000 genes/cell) in 16 588 single cells isolated during atherosclerosis progression in Ldlr^-/-Apob^100/100 mice with human-like plasma lipoproteins and from humans with asymptomatic and symptomatic carotid plaques was clustered into multiple subtypes. For clinical and pathophysiological context, the advanced-stage and symptomatic subtype clusters were integrated with 135 tissue-specific (atherosclerotic aortic wall, mammary artery, liver, skeletal muscle, and visceral and subcutaneous, fat) gene-regulatory networks (GRNs) inferred from 600 coronary artery disease patients in the STARNET (Stockholm-Tartu Atherosclerosis Reverse Network Engineering Task) study.

Results: Advanced stages of atherosclerosis progression and symptomatic carotid plaques were largely characterized by 3 smooth muscle cells (SMCs), and 3 macrophage subtype clusters with extracellular matrix organization/osteogenic (SMC), and M1-type proinflammatory/Trem2-high lipid-associated (macrophage) phenotypes. Integrative analysis of these 6 clusters with STARNET revealed significant enrichments of 3 arterial wall GRNs: GRN33 (macrophage), GRN39 (SMC), and GRN122 (macrophage) with major contributions to coronary artery disease heritability and strong associations with clinical scores of coronary atherosclerosis severity. The presence and pathophysiological relevance of GRN39 were verified in 5 independent RNAseq data sets obtained from the human coronary and aortic artery, and primary SMCs and by targeting its top-key drivers, FRZB and ALCAM in cultured human coronary artery SMCs.

Conclusions: By identifying and integrating the most gene-rich single-cell subclusters of atherosclerosis to date with a coronary artery disease framework of GRNs, GRN39 was identified and independently validated as being critical for the transformation of contractile SMCs into an osteogenic phenotype promoting advanced, symptomatic atherosclerosis.

Keywords: coronary artery disease; gene expression; lipoportein; macrophages; subcutaneous fat.

Figure 1.

Schematic overview of the overall study design. FACS indicates fluorescence-activated cell sorting; GRN, gene-regulatory network; SMC, smooth muscle cell; and STARNET, Stockholm-Tartu Atherosclerosis Reverse Network Engineering Task.

Figure 2.

Single-cell RNA sequencing (scRNAseq) characterization of the subcellular transformation taking place during atherosclerosis progression in the Ldlr^−/−Apob^100/100 mouse model. scRNAseq data were generated from pairs of aortic arches during atherosclerosis progression obtained from a total of 46 mice at the 10-week baseline (n=8 mice), the 20- to 30-week early atherosclerosis stage (n=18 mice), and the 45- to 60-week advanced atherosclerosis stage (n=20 mice). Aortic arches, dissected from the aortic root to the third rib in anesthetized and euthanized mice, were reduced to single-cell suspensions in a 384-well plate format (1 plate per cell type/mouse). Each well was then incubated with fluorophore-conjugated antibodies specific to the major vascular cell types; anti-PDGFRb/SMCs, anti-CD31/ECs, and anti-CD45/immune cells. Full-length RNA sequencing of the single-cell suspensions was performed by SmartSeq2 using the Illumina Hiseq3000 platform (single read, 50-bp length). A through C, Clusters of individual B-cells, endothelial cells (ECs), macrophages (MPs), myofibroblasts (FBs), smooth muscle cells (SMCs), and T-cells from the atherosclerosis-free baseline (A), the early atherosclerosis stage (B), and the advanced atherosclerosis stage (C) visualized by Uniform Manifold Approximation and Projection (UMAP). Arrows indicate major changes to the cell clusters during the progression of atherosclerosis. D, Progression of atherosclerosis in the chow-fed, Ldlr^−/−Apob^100/100 mice measured as the percentage oil-red-O-stained lesion area in relation to the total surface area in pinned-out aortas. Adopted from Björkegren et al. E, Expression heat map showing the 5 most significant markers of each cell type. Wilcoxon Rank-Sum test was used to establish significantly upregulated markers comparing expression levels in each cell type against all other cell types. The scalebar indicates fold-gene expression levels relative background. F, Immunofluorescence-staining (αSMA for SMCs [red], PECAM1 [platelet endothelial cell adhesion molecule 1] for ECs [green], and CD68 for MPs [blue]) from a Ldlr^−/−Apob^100/100 mouse with average-size aortic root lesion at the advanced stage of atherosclerosis. The arrow points to the atherosclerotic lesion. Scale bar bottom right indicates 75 µm. G, Dot plots and a Venn diagram of cell type–specific genes enriched in the advanced stage of atherosclerosis progression. Cell type–specific gene enrichment was calculated using Wilcoxon Rank-Sum test comparing gene expression levels between baseline/early stage, and advanced stage, atherosclerosis (log₂-fold >0.3, Bonferroni-adjusted P<0.005). Dot plots show top-ranked biological processes of cell-specific and combined advanced stage atherosclerosis genes indicated in the Venn diagram. Gene ratio (x-axis) indicates the number of overlapping genes divided by the total number of genes in the gene ontology (GO) category. GO enrichment −log10 P values were calculated with Fisher exact test. Genes indicate the number of differential expressed genes for cell type overlapping with the GO category.

Figure 3.

Single-cell RNA sequencing (scRNAseq) characterization of the subcellular transformation taking place in human atherosclerotic carotid plaques. scRNAseq data were generated from carotid plaques obtained from a total of 15 asymptomatic (n=7) and symptomatic (n=8) patients during carotid endarterectomy surgery. Carotid plaque cells were dissociated and incubated with fluorophore-conjugated antibodies: anti-PDGFRb (Platelet Derived Growth Factor Receptor Beta) for smooth muscle cells (SMCs), anti-CD31/CD144 for endothelial cells (ECs) and anti-CD45 for immune cells. Full-length RNA sequencing of single-cell suspensions in a 384-well plate format (1 plate per cell type/plaque) was performed by SmartSeq2 using the Illumina Hiseq3000 platform (single read, 50-bp length). A, Clusters of carotid plaque ECs, macrophages (MPs), pericytes (PCs), smooth muscle cells (SMCs), and T-cells visualized with Uniform Manifold Approximation and Projection (UMAP). B, Expression heat map showing top-5 ranked markers for each cell type. Wilcoxon Rank-Sum test was used to establish significantly upregulated markers comparing expression levels in each cell type against all other cell types. Scalebar indicates fold-gene expression levels relative background. C, Bar plots showing the relative contribution of symptomatic and asymptomatic cells of the main cell types in the carotid plaques. Values in parenthesis indicate total number of collected cells for each cell type. D, Dot plots and a Venn diagram showing the number and overlap of cell type–specific genes enriched in carotid plaques isolated from symptomatic patients. Cell-specific gene enrichment was calculated using a Wilcoxon Rank-Sum test comparing gene expression levels between carotid plaques isolated from symptomatic and asymptomatic patients (log₂-fold >0.5; Bonferroni adjusted, P<0.005). Dot plots show top-ranked biological processes of cell-specific and combined symptomatic carotid plaque genes indicated in the Venn diagram. Gene ratio (x-axis) indicates the number of overlapping genes divided by the total number of genes in the gene ontology (GO) category. GO enrichment −log10 P values were calculated with Fisher exact test. Genes indicate the number of differential expressed genes for cell type overlapping with the GO category.

Figure 4.

Single-cell RNA sequencing (scRNAseq) characterization of smooth muscle cell (SMC) subclusters active during atherosclerosis progression in Ldlr^−/−Apob^100/100 mice and in carotid plaques from humans with or without neurological symptoms. A and B, Subclusters (A) and trajectory (B, color coded by pseudotime) visualized with Uniform Manifold Approximation and Projection (UMAP) of SMCs during atherosclerosis progression in Ldlr^−/−Apob^100/100 mice. A, top left, Main SMC cluster in Figure 1A. C, Dot plot showing expression levels of established markers in the seven SMC subclusters active during atherosclerosis progression in mice. Average expression indicates the average fold expression levels of the marker in the indicated subcluster compared with all other clusters as background. % expressed, indicates the percentage of cells in each subcluster where the gene marker is expressed. D, Bar plot showing the relative cell contributions of the seven SMC clusters (x-axis) in relation to progression of atherosclerosis in Ldlr^−/−Apob^100/100 mice (weeks 10 to 60). Values in parenthesis indicate the total number of cells obtained in each of the seven SMC subclusters. E, Subcluster and trajectory of SMCs isolated from the human carotid plaques. A, bottom right, The main SMC cluster in Figure 2A. F, Violin plots of expression levels of the SMC marker of contractility, CNN1 (top) and the fibroblast marker DCN (bottom) in the 8 SMC clusters (x-axis) identified in human carotid plaques. G, Bar plot showing the relative contribution of carotid plaque cells in the 8 SMC subcluster isolated from asymptomatic (n=7) or symptomatic (n=6) carotid plaques. Values in parenthesis indicate the total number of cells obtained in each SMC subclusters. H, Venn diagram of subcluster genes in mSMC6/mSMC7 (D) and hSMC7 (G). Dot plot shows top-ranked biological processes of common mSMC6/mSMC7/hSMC7 genes. Gene ratio (x-axis) indicates the number of overlapping genes divided by the total number of genes in the gene ontology (GO) category. GO enrichment −log10 P values were calculated with Fisher exact test. Genes indicate the number of differential expressed genes for cell type overlapping with the GO category. I, Ingenuity pathway analysis of mSMC6/mSMC7/hSMC7 (human smooth muscle cell subcluster 7) genes in relation to genes previously identified in MP clusters of advanced-stage atherosclerosis.^,

Figure 5.

Single-cell RNA sequencing (scRNAseq) characterization of macrophage (MP) subclusters active during atherosclerosis progression in Ldlr^−/−Apob^100/100 mice and in carotid plaques from humans with or without neurological symptoms. A and B, Subclusters (A) and trajectory (B, color coded per pseudotime) of MPs during atherosclerosis progression in Ldlr^−/−Apob^100/100 mice (n=18). A, top right, the main MP cluster in Figure 1A. C, Dot plot showing expression levels of well-established markers in the 6 MP subclusters active during atherosclerosis progression in mice. Average expression indicates the average fold expression levels of the marker in the indicated subcluster compared with all other clusters as background. % expressed indicates the percentage of cells in each subcluster where the gene marker is expressed. D, Bar plot showing the relative contributions of cells to the 6 MP clusters (x-axis) in relation to progression of atherosclerosis in Ldlr^−/−Apob^100/100 mice (weeks 10–60). Values in parenthesis indicate total number of cells in each MP subcluster. E, Subclusters with cell trajectory visualized of MPs isolated from human carotid plaques (n=15). E, bottom right, The MP cluster shown in Figure 2A. F, Violin plots of expression levels of the monocyte-marker F13A1 and the M2-macrophage marker SPP1 in the 8 MP subclusters (x-axis) of the human carotid plaques. G, Bar plot showing the relative contribution of cells in the 8 MP cluster isolated from asymptomatic (n=7) or symptomatic (n=6) carotid plaques. Values in parenthesis indicate the number of cells isolated in each of the MP subcluster. H, Venn diagram showing of subclusters mMP5/mMP6 (D) and hMP7 (human macrophage subcluster 7; G) genes. Dot plot shows top-ranked biological processes of common mMP5/mMP6/hMP7 genes. Gene ratio (x-axis) indicates the number of overlapping genes divided by the total number of genes in the GO category. GO enrichment −log10 P values were calculated with Fisher exact test. Genes indicate the number of differential expressed genes for cell type overlapping with the GO category. I, Ingenuity pathway analysis of mMP5/mMP6/hMP7 genes in relation to genes previously identified in MP clusters of advanced stage atherosclerosis.^,

Figure 6.

Enrichment of 6 late-stage/symptomatic macrophage (MP) and smooth muscle cell (SMC) subcluster genes in arterial wall gene-regulatory networks (GRNs) of patients with coronary artery disease (CAD). One hundred thirty-five tissue-specific GRNs inferred from genotype and bulk RNAseq data of 2 arterial wall (atherosclerotic aortic root [AOR], and mammary artery [MAM]) and 4 metabolic (liver [LIV], skeletal muscle [SKLM], subcutaneous fat [SF], and visceral abdominal fat [VAF]) tissues obtained from 600 patients with CAD of the STARNET (Stockholm-Tartu Atherosclerosis Reverse Network Engineering Task) study were inferred using blockwise weighted gene co-expression network analysis (WCGNA) and the GENIE3 algorithms, as described. A, Dot plots showing GRNs (x-axis) with gene enrichment (y-axis, −log10 false discovery rates [FDRs]) for genes in the advanced stage and symptomatic SMC (left, n=12 GRNs) and MP (right, n=10 GRNs) subclusters. B, Horizontal bar grafts showing the statistical significances (x-axis, −log10 FDRs) of associations with Synergy Between Percutaneous Coronary Intervention With Taxus and Cardiac Surgery (SYNTAX; left) and Duke (right) scores of the arterial wall–specific GRNs in A. C, Pie charts showing the relative cell-type specificity of GRN39 (top), GRN33 (middle), and GRN122 (bottom) and their key driver genes according to the scRNAseq data (Methods). D, The arterial wall–specific GRN39 (top), GRN33 (middle), and GRN122 (bottom) color-coded according to cell-specificity (C). Bigger size nodes are key drivers. Bar plots show examples of top-key driver expression patterns during atherosclerosis progression in the Ldlr^−/−Apob^100/100 mice, and their corresponding expression pattern in subtype cell clusters according to Uniform Manifold Approximation and Projection (UMAP). Scale bars in red, indicates fold expression levels of the individual key drivers in the indicated subclusters compared with background subclusters of the same cell type. Top/rank, indicates key driver’s hierarchical ranking in the GRN. H², broad sense heritability contribution per GRN (%). E, Radar plot showing statistical significance of key cardiometabolic phenotype associations with GRN39 (top), GRN33 (middle), and GRN122 (bottom) and below each radar plot; abbreviations of the GWAS CAD candidate genes present in each GRN. The significance of GRN-phenotype associations (−log10; P=0–100) was calculated by aggregating GRN gene-level phenotype associations (Pearson correlation coefficients, Student t test) corrected for the total number of STARNET GRNs (n=135) and the number of genes in each GRN using the Benjamini-Hochberg procedure.

Figure 7.

GRN39 replication and evaluation in 5 independent human arterial wall datasets. A, GRN39 replication according to the NetRep algorithm comparing the original STARNET (Stockholm-Tartu Atherosclerosis Reverse Network Engineering Task) data with 2 coronary (University of Virginia/Stanford and Genotype-Tissue Expression project [GTEx]), and 1 aortic root (GTEx) bulk RNAseq datasets. cor.cor measures the concordance of correlation structures; cor.degree measures the concordance of the weighted degrees; cor.contrib measures the concordance of node contributions; avg.cor measures the average magnitude of correlation coefficients and avg.contrib measures the average magnitude of node contributions. Bar shows log⁻¹⁰ P values. B, Bar plots showing the odds ratio (OR) of enrichment of GRN39 for differentially differential expressed genes in coronary artery tissue samples annotated as (1) lesions vs nonlesions, (2) ischemic cardiomyopathy vs controls, or (3) combined lesion and ischemic (Methods). Fisher test was performed on overlapping gene sets, P values shown represent the Fisher test P values, and error bars represent the 95% CIs. C, GRN39 gene cell-type specificity according to enrichment of cell-type determined genes in scATACseq data obtained from 41 human coronary artery tissue samples (Methods). Dot plot showing the odds ratio of enrichment of GRN39 module genes for marker genes determined from integrated snATACseq and scRNAseq profiles obtained from atherosclerotic coronary artery tissues (Methods). D, Bar plots showing the distributions of pairwise spearman correlations of GRN39 genes compared with random background genes in RNAseq data from quiescent and proliferative human vascular smooth muscle cells (SMCs) obtained from 151 donors (Methods). Raw P values were calculated using Kruskal-Wallis test. E, Pie chart showing the proportion of GRN39 genes found to affect SMC-related atherosclerosis-related phenotypes, as described. snATAC indicates single nucleus assay for transposase-accessible chromatin.

Figure 8.

Experimental validation of GRN39 by perturbing its top key drivers FRZB and ALCAM in human coronary artery smooth muscle cells (HCASMCs). The top key drivers of GRN39, FRZB (Frizzled Related Protein), and ALCAM were overexpressed and knocked down, respectively, in HCASMCs that were either untreated or treated with transforming growth factor (TGF)-beta or cholesterol to, respectively, to induce proliferative and synthetic conditions followed by RNAseq. A, Average pair-wise correlation of GRN39 coding genes compared with the average pairwise correlations of 1000 equally sized random groups of background (Bkg) coding genes from the same dataset (ie, mock-treated [Mock], silencing RNA-treated ALCAM [siALCAM] or overexpressed FRZB [oeFRZB] RNAseq datasets). B, Average pair-wise correlation of GRN39 coding genes compared between Mock, siALCAM, and oeFRZB. C, Gene set enrichment analysis (GSEA) of RNAseq data of HCASMCs overexpressing FRZB and downregulating ALCAM. The conditions representing significant enrichment of GRN39 module genes are shown. D, Relative proliferation in the FRZB, compared with the control (CTRL) group, was significantly reduced using a paired 2-tailed t test (P=0.016) and borderline significant using a Wilcoxon 2-tailed matched-pairs signed rank test (P=0.13; n=4, mean and SD are shown). E, Fluorescence intensity indicative of calcification deposit levels in cells transduced with CTRL or FRZB lentivirus following a 14-day inorganic phosphate-induced calcification. Statistical analysis was performed using the Mann-Whitney U test. Both mean and SD are shown (n=6).