Integrative single-cell meta-analysis reveals disease-relevant vascular cell states and markers in human atherosclerosis

Abstract

Coronary artery disease (CAD) is characterized by atherosclerotic plaque formation in the arterial wall. CAD progression involves complex interactions and phenotypic plasticity among vascular and immune cell lineages. Single-cell RNA-seq (scRNA-seq) studies have highlighted lineage-specific transcriptomic signatures, but human cell phenotypes remain controversial. Here, we perform an integrated meta-analysis of 22 scRNA-seq libraries to generate a comprehensive map of human atherosclerosis with 118,578 cells. Besides characterizing granular cell-type diversity and communication, we leverage this atlas to provide insights into smooth muscle cell (SMC) modulation. We integrate genome-wide association study data and uncover a critical role for modulated SMC phenotypes in CAD, myocardial infarction, and coronary calcification. Finally, we identify fibromyocyte/fibrochondrogenic SMC markers (LTBP1 and CRTAC1) as proxies of atherosclerosis progression and validate these through omics and spatial imaging analyses. Altogether, we create a unified atlas of human atherosclerosis informing cell state-specific mechanistic and translational studies of cardiovascular diseases.

Keywords: CP: Genomics; atherosclerosis; coronary artery disease; genome-wide association studies; integration analyses; single cell RNA-seq; smooth muscle cells.

Figure 1.. Pipeline devised to build the integrated human scRNA atherosclerosis atlas

Doublets scRNA libraries were identified using scDblFinder while ambient RNA was removed using decontX. Decontaminated matrices were used for downstream filtering of cells based on 1) detected gene number 2) UMI number 3) % reads mapping to mitochondrial genome 4) % reads mapping to hemoglobin genes using Seurat. Libraries were normalized using SCTransform and a subset used for integration benchmarking using four approaches (Canonical Correlation Analysis + Mutual Nearest Neighbors (CCA+MNN), reciprocal PCA (PCA), Harmony and Scanorama . PCA embeddings were used for calculating the following metrics: running time, integration and cell-type local inverse simpson index (iLISI and cLISI) scores, k-Nearest Neighbors Batch Effect Test (kBET) rejection rates, principal component regression coefficients and cluster average silhouette widths (cASW). rPCA was used to harmonize the 22 libraries and level 1 annotations were added using Transfer learning with the Tabula Sapiens Vasculature dataset as reference. Finally, automated and manual curation enabled identification of granular cell subtypes (level 2 annotations). See also Method details, Figure S1 and Table S1.

Figure 2.. Integration of single cell data identifies major cell compartments in human atherosclerosis

(A) UMAP embeddings of 118,578 cells based on reciprocal PCA (rPCA) integration of 22 sequencing libraries. Dot colors depict broad cell lineage annotations (level 1). (B) Dot plot of top five marker genes by major cell lineage compartment. Dot size represents the portion of cells expressing the gene per level 1 compartment. (C) Stacked bar plot showing the distribution of level 1-annotated cells across included studies, arterial beds (coronary/carotids), and lesion status (lesion, non-lesion). (D) UMAP embeddings of level1-annotated cells across lesion status. (E) Heatmap showing ES_μ values for established markers of mural and immune level 1 cell annotations previously identified as differentially expressed (DE) using a Wilcoxon Rank sum test in (B). (F) LD Score Regression applied to specific genes (LDSC-SEG) and MAGMA analyses prioritizing the contribution of level 1-annotated cell type to Cardiovascular and non-Cardiovascular GWAS traits. LDSC-SEG- and MAGMA-based regression analyses were carried out using an expression specificity matrix generated with CELLEX. The black line depicts the FDR significance threshold (FDR < 5% at -log10(P) = 1.301). See also Figures S2, S3 and Tables S2 and S4.

Figure 3.. Atherosclerosis cell subpopulations (level 2) and distribution of myeloid subtypes across disease status

(A) UMAP representation of cell subtypes (level 2 annotations) within the largest level 1 cell compartments (T/NK, Macrophages, Endothelial, Fibroblast). (B) UMAP embeddings of canonical marker genes delineating immune and non-immune cell subtypes. SCTransform-normalized gene expression is indicated by color. (C) UMAP embeddings and bar plot of level 2 myeloid subtypes according to lesion status. Frequencies for each subtype are normalized to the total number of cells in each condition (lesion n=59691; non-lesion n=58887) and shown as percentages. See also Figure S2.

Figure 4.. Characterization of etiologic SMC phenotypes for cardiovascular traits and diseases

(A) UCell enrichment of meta-analyzed SMC murine gene modules (add ref) (Contractile, Lgals3+ transitional, Fibrochondrocytes) and non-SMC-derived fibroblasts in the level 1 SMC compartment. A subset of human fibroblasts and pericytes were used as negative enrichment controls for murine SMC modules. UCell scores were calculated using the Mann-Whitney U statistic. (B) UMAP embeddings of SMC level 2 annotations in addition to pericytes and subset of fibroblasts. Annotations were defined using UCell scores to guide SMC differentiation state in addition to differentially expressed markers from Louvain clusters. (C) Dot plot showing expression of top marker genes after SCTransform normalization for SMC level 2 annotations. Dot size represents the portion of cells expressing the gene. (D) LDSC-SEG and MAGMA analyses prioritizing the contribution of SMC phenotypes, pericytes, and fibroblasts to cardiovascular GWAS traits. Type 2 diabetes and AD were used as negative controls. The black line depicts the FDR significance threshold (FDR < 5% at -log10(P) = 1.301). (E) Meta-analysis UMAP embeddings showing normalized scDRS scores for CAD and immune traits (WBC count and AD) previously shown as highly enriched in level 1 myeloid annotations. Red indicates cells enriched for the above mentioned traits while non-relevant cells are denoted in dark blue. WBC=White Blood Cell; AD=Alzheimer’s disease. See also Figure S4 and Tables S2 and S3.

Figure 5.. Summary of cell crosstalk in human atherosclerosis

(A) Circle plots depicting aggregated cell-cell communication network for level 1-labeled cell compartments leveraging the CellChat human database. Interactions include secreted signaling, ECM-receptor and cell-cell contacts. Interactions were calculated separately across disease status (lesion vs non-lesion). Top 30% of interactions are shown in the plot. (B) Stacked bar plot showing conserved and disease status-specific signaling pathways. Signaling enrichment is based on pathways information flow changes (sum of communication probability among all pairs of cell groups in the inferred network). Pathways in bold denote those that showed statistically significant (P < 0.05) enrichments in each disease condition. (C) Circle plot depicting sources and targets for SPP1 signaling using level 2 annotations for myeloid cells and level 1 SMC annotations. (D) Circle plot depicting SPP1 signaling sources and targets for level 2 myeloid and SMC annotations. (E) Summary dot plot of ligand-receptor interactions for level 2 myeloid and SMC annotations. Myeloid subtypes were defined as signaling sources while SMCs were defined as signaling targets. Width of the edges in the circle plot depicts the weight/strength of the interactions in (A, C and D). ECM = Extracellular matrix. See also Figure S5 and Table S5.

Figure 6.. Pseudotime and TF inference activity for ECM-rich SMC phenotypes

(A) UMAP embeddings showing the contractile-to-modulated SMC pseudotime trajectory calculated with Monocle3. SMC phenotypes for this analysis included contractile, transitional SMCs, fibromyocytes and fibrochondrocytes (FCs). The numbered circle depicts the trajectory root defined as the contractile SMC subset with highest MYH11 expression. (B) Heatmap showing the top 500 SMC differentially-expressed (DE) genes across pseudotime. Bolded genes depict GWAS CAD SMC effector genes differentially expressed across pseudotime. The black arrow indicates gradual increase in pseudotime towards the right. (C-D) Cubic spline interpolation of SCTransform-normalized gene expression as a function of pseudotime. Genes plotted include hits from Monocle 3 and Seurat DE tests (FDR < 0.05). Genes from SMCs to fibromyocytes trajectory: FN1, LGALS3, AEBP1, LTBP1, PDGFRB. From SMCs to FC trajectory: COL1A2, IBSP, CRTAC1, COMP, MMP2. (E) TF activity prediction with VIPER leveraging DoRothEA regulons for SMC phenotypes. Only regulons with high confidence scores (A-C) were used for this analysis. Highly variable TFs were selected for plotting and scale indicates relative predicted activity. TF = Transcription Factor. See also Method details and Figures S3 and S6.

Figure 7.. CRTAC1 and LTBP1 are candidate markers of atherosclerosis progression

(A) UMAP embeddings of CRTAC1 and LTBP1 expression within the subset defined in Figure 4A. IBSP and VCAN are provided as calcification and fibromyocyte markers, respectively. (B) Pearson correlation plot of CRTAC1 versus all genes across contractile SMC and FC clusters. Examples of osteochondrogenic and contractile genes are shown. (C) Pearson correlation plot of LTBP1 versus all genes across contractile SMC and fibromyocyte clusters. Examples of contractile and ECM-related genes upregulated during SMC modulation are shown. (D) Bulk RNA-seq expression of CRTAC1, IBSP, LTBP1 and TCF21 in human coronary arteries from lesion (n=27) and non-lesion samples (n=21). The Y axis represents normalized expression counts (TPMs). P values calculated using a non-parametric Wilcoxon Rank Sum Test. (E) Log-normalized protein expression of CRTAC1, SPP1, LTBP1, and VCAN in non-lesion (n=27) and lesion (n=29) human coronary segments. P values calculated using an unpaired Student’s t-test. Boxplots in (D) and (E) represent the median and the inter-quartile (IQR) range with upper (75%) and lower (25%) quartiles shown and each dot represents a separate individual. (F) Clinical trait enrichment for CRTAC1-containing module in subclinical mammary artery STARNET gene regulatory network. Pearson’s correlation P values (gene-level) were aggregated for each co-expression module using a two-sided Fisher’s exact test. Case/control differential gene expression enrichment was estimated by a hypergeometric test. (G) Representative H&E and LTBP1 IF staining in human atherosclerotic coronary artery segments from normal (control/subclinical, n=4) and stage IV-V lesions (n=4) (Stary classification). Two regions of interest (ROI) per sample stained with LTBP1 are shown on the right, respectively at 10x and 20x magnification, with LTBP1 (red), alpha-smooth muscle actin (a-SMA; green), and DAPI-stained nuclei (blue). Arrows point to LTBP1+ cells. Scale bars = 0.1mm (left) and 30μm (right) in ROIs. (H) Representative H&E and RNAscope staining of LTBP1 and CRTAC1 in human coronary artery segments from normal (control/subclinical, n=2) and stage IV-V lesions (n=4) (Stary classification). Two ROIs per condition are shown on the right at 40x magnification, with LTBP1 (red), LMOD1 (green), CRTAC1 (white), and DAPI-stained nuclei (blue). Pink arrows point to LTBP1+/LMOD1+ cells while white arrows point to CRTAC1+ cells. Scale bars = 0.1mm (left) and 30μm (right) in ROIs in RNAscope. Scale bars = 1mm for H&E images in (G) and (H). See also Figures S4, S7, S8, S9 and Table S6.