Dissecting the polygenic basis of atherosclerosis via disease-associated cell state signatures

Abstract

Coronary artery disease (CAD) is a pandemic disease where up to half of the risk is explained by genetic factors. Advanced insights into the genetic basis of CAD require deeper understanding of the contributions of different cell types, molecular pathways, and genes to disease heritability. Here, we investigate the biological diversity of atherosclerosis-associated cell states and interrogate their contribution to the genetic risk of CAD by using single-cell and bulk RNA sequencing (RNA-seq) of mouse and human lesions. We identified 12 disease-associated cell states that we characterized further by gene set functional profiling, ligand-receptor prediction, and transcription factor inference. Importantly, Vcam1+ smooth muscle cell state genes contributed most to SNP-based heritability of CAD. In line with this, genetic variants near smooth muscle cell state genes and regulatory elements explained the largest fraction of CAD-risk variance between individuals. Using this information for variant prioritization, we derived a hybrid polygenic risk score (PRS) that demonstrated improved performance over a classical PRS. Our results provide insights into the biological mechanisms associated with CAD risk, which could make a promising contribution to precision medicine and tailored therapeutic interventions in the future.

Keywords: GWAS; atherosclerosis; cell state; coronary artery disease; genetics; genome-wide association study; polygenic risk score; scRNA-seq; single cell.

Graphical abstract

Figure 1

Clustering and identification of atherosclerosis-associated cell states (A) Schematic overview of the experimental setup. (B and C) (B) UMAP projection of the scRNA-seq profiles represented as eleven manually annotated clusters. (C) UMAP regional occupancy analysis demonstrating relative changes in cell density comparing atherosclerotic vascular wall to healthy controls. Atherosclerosis-associated cell states are revealed by increased local abundance of cells in regions of the UMAP plot (log2FC, log2 fold change). (D) UMAP plot depicting the 12 disease-associated cell states and the selected top marker genes. (E) Relative changes in the cell state proportions during different stages of atherosclerosis shown for each of the three biological replicates. Diamond represents the average of the three replicates.

Figure 2

Atherosclerosis-associated cell state signatures are activated in mouse and human lesions based on bulk RNA-seq (A–C) (A) Cell population mapping using the scRNA-seq cell state space to plot the differential cell abundance estimated from the mouse bulk RNA-seq data. The gene set activity scores for each cell state were further investigated in the (B) Tampere Vascular Study representing 68 advanced atherosclerotic plaques (15 aortic, 29 carotid, and 24 femoral plaques) and 28 controls (left internal thoracic artery) and (C) the Maastricht Pathology Tissue Collection representing atherosclerotic carotid artery segments from 13 early intimal thickening/xanthoma lesions and from 16 advanced fibrous cap atheroma lesions.

Figure 3

Characterization of atherosclerosis-associated cell states and key biological pathways (A) Marker gene counts for the 12 most abundant disease-associated cell states. (B) Common markers between disease-increased cell states. (C) UpSet plot showing the gene overlaps between cell state signatures and the gene ontologies enriched (log10 of adjusted p value) in the intersections. (D) Single-cell gene set enrichment scoring for selected biological processes. Enrichment score is shown from the lowest 5% (Q5) to highest 95% (Q95).

Figure 4

Modeling intercellular communication between cell states Ligand–receptor–target gene analysis was carried out with NicheNet. (A) Top 10 prioritized upstream ligands for cell state signature gene sets. (B) Ligand target gene networks presented for Vcam1+ SMCs, Col2a1+ SMCs, and Lrg1+ ECs. (C) Expression of the prioritized ligands by the different cell states and types. Row normalized gene expression (TPM = transcripts per million) is shown. (D) Schematic of predicted ligand-mediated signaling between Vcam1+ and Col2a1+ SMC states involving autocrine and paracrine signaling.

Figure 5

Prediction of cell-state-specific transcription factor activities (A–D) (A) The most abundant atherosclerosis-associated cell states were selected for SCENIC analysis along with disease-unperturbed cells of the same cell type. Differentially active gene regulatory networks identified for (B) smooth muscle cell (SMC), (C) endothelial cell (EC), and (D) macrophage (MP) cell states. The predicted regulon activity and transcription factor gene expression (row normalized TPM) are shown. (E) Selected examples of regulon activities and transcription factor gene expression plotted on UMAP.

Figure 6

The contribution of cell states to CAD heritability (A) Overlap of the CAD GWAS candidate causal gene lists from nine different sources (different colors) with the cell state markers. NA indicates no overlapping genes. (B) Results from LD score regression (LDSC) applied to cell state marker genes (using 100 kb flanking regions) to partition CAD heritability within the genome.

Figure 7

Pathway- and cell-state-specific polygenic risk scores shed light into the genetic basis of CAD (A) Cell-state-specific PRS was constructed with (A) gene coordinates (−35 kb upstream to 10 kb downstream) using PRSet. To obtain the empirical p value, random SNP sets containing the same number of post-clumping SNPs were selected from background regions of the genome, selected from all genic regions. (B) Explained variance of each pathway-specific PRS to polygenic risk of CAD calculated for the gene sets listed in Figure 3C. (C) Cell-state-specific PRS analysis constructed with plaque scATAC-seq peak coordinates that were found within ±500 kb of the TSS. (D) Proportion of variance of CAD explained by PRS in genome-wide analysis. The values represent PRS calculated for all cell-type-specific scATAC-seq peaks at different p value thresholds, which are compared to the classical genome-wide clumping and thresholding PRS.