Intersecting single-cell transcriptomics and genome-wide association studies identifies crucial cell populations and candidate genes for atherosclerosis.

AIMS: Genome-wide association studies (GWASs) have discovered hundreds of common genetic variants for atherosclerotic disease and cardiovascular risk factors. The translation of susceptibility loci into biological mechanisms and targets for drug discovery remains challenging. Intersecting genetic and gene expression data has led to the identification of candidate genes. However, previously studied tissues are often non-diseased and heterogeneous in cell composition, hindering accurate candidate prioritization. Therefore, we analysed single-cell transcriptomics from atherosclerotic plaques for cell-type-specific expression to identify atherosclerosis-associated candidate gene-cell pairs. METHODS AND
RESULTS: We applied gene-based analyses using GWAS summary statistics from 46 atherosclerotic and cardiovascular disease, risk factors, and other traits. We then intersected these candidates with single-cell RNA sequencing (scRNA-seq) data to identify genes specific for individual cell (sub)populations in atherosclerotic plaques. The coronary artery disease (CAD) loci demonstrated a prominent signal in plaque smooth muscle cells (SMCs) (SKI, KANK2, and SORT1) P-adj. = 0.0012, and endothelial cells (ECs) (SLC44A1, ATP2B1) P-adj. = 0.0011. Finally, we used liver-derived scRNA-seq data and showed hepatocyte-specific enrichment of genes involved in serum lipid levels.
CONCLUSION: We discovered novel and known gene-cell pairs pointing to new biological mechanisms of atherosclerotic disease. We highlight that loci associated with CAD reveal prominent association levels in mainly plaque SMC and EC populations. We present an intuitive single-cell transcriptomics-driven workflow rooted in human large-scale genetic studies to identify putative candidate genes and affected cells associated with cardiovascular traits. Collectively, our workflow allows for the identification of cell-specific targets relevant for atherosclerosis and can be universally applied to other complex genetic diseases and traits.

Genome-wide association studies (GWAS) identified a large number of genomic loci associated with atherosclerotic disease. The translation of these results into drug development and faster diagnostics remains challenging. With our approach, we cross-reference the GWAS findings for atherosclerotic disease with single-cell RNA sequencing data of disease-relevant tissue and bring the GWAS findings closer to the functional and mechanistic studies.

Introduction

Since the first publication in 2005 on age-related macular degeneration, genome-wide association studies (GWAS) soared, and to date, the GWAS Catalogue registered over 5000 publications covering thousands of genetic associations. These efforts also identified hundreds of common genetic risk loci for cardiovascular diseases and traits, including coronary artery disease (CAD), ischaemic stroke subtypes, coronary artery calcification, and carotid intima and media thickness. Despite tremendous socio-economic and medical progress, cardiovascular disease remains the leading cause of death, highlighting the need for new biomarkers and pharmaceutical targets. Previous studies convincingly show that GWAS offer a solid foundation for innovative drug development as they agnostically identify genetic relationships with diseases and point to potential causal mechanisms.

Identifying the causal genes and underlying pathological mechanisms that drive the signal in genetic risk loci is not straightforward. One possible approach to narrow down the lists of candidate genes is to utilize gene expression data, by identifying expression quantitative trait loci (eQTLs) that colocalize with genetic risk loci with expression changes in relevant tissues. Alternatively, transcriptome-wide association studies estimate genetic correlation between tissue-specific gene expression and complex traits to predict causal relations between traits. However, these powerful methods take advantage of whole-tissues collected in large biobanks and the variable cell composition confounds the resulting signal. Because tissues are composed of many cell types and tissue expression is averaged across all those cells, the signal from less prevalent but important cell populations is obscured. In order to test potential therapeutic targets, one requires suitable functional assays. This in turn requires an idea of the function and cellular expression of targets, which is complicated by tissue heterogeneity. Single-cell transcriptomics enables measurement of cell-specific expression, cell-to-cell expression variability, transcriptional noise, or temporal dynamics. Intersecting GWAS data with human transcriptomic single-cell data would increase the resolution both in terms of candidate gene and more importantly, it would identify the cell type most relevant for the given disease.

We hypothesized that genetic drivers for atherosclerotic disease can be found in plaque tissue. We aimed to project GWAS signals for 46 traits directly into a cell-type-specific expression using single-cell RNA sequencing (scRNA-seq) data derived from atherosclerotic plaques. To this end, we devised a workflow that utilized gene-based testing to identify candidate genes from GWAS summary statistics and calculated an enrichment score for each candidate list using expression profiles from scRNA-seq driven cell populations. This approach allows us to find the enrichment in expression patterns that point to specific, and potentially causal, cell populations. We report that loci associated with CAD reveal prominent association in mainly plaque smooth muscle cells (SMCs) and endothelial cells (ECs).

Methods

Single-cell sequencing: carotid artery plaques

Atherosclerotic lesions were collected from 12 female and 26 male patients undergoing a carotid endarterectomy procedure. All pathological tissue was included in the Athero-Express Biobank Study biobank (AE, www.atheroexpress.nl) at the University Medical Centre Utrecht (UMCU). Plaque processing for single-cell sequencing is described elsewhere. After sequencing, data were processed in an R 3.5 and 4.0 environment using Seurat version 3.2.2. Mitochondrial genes were excluded and doublets were omitted by gating for unique reads per cell (between 400 and 10 000) and total reads per cell (between 700 and 30 000). Data were corrected for sequencing batches using the function SCTransform. Clusters were created with 30 principal components at resolution 0.8. Multiple iterations of clustering were used to determine optimal clustering parameters. Population identities were assigned by evaluating gene expression per individual cell clusters. Sub-populations of SMCs were determined by isolating the SMC and EC populations from the complete population of cells mentioned above and re-assigned clusters. Clusters were created with 10 principal components at resolution 1.4. New subpopulation identities were assigned by evaluating gene expression per individual cell clusters. The newly defined SMC subpopulation identities were projected onto the SMC populations in the complete scRNA-seq dataset of atherosclerotic lesions for downstream processing.

Genome-wide association studies summary statistics

Genome-wide association studies summary statistics were leveraged from publicly available datasets. Mining these datasets uncovered a plethora of challenges. Not all GWAS summary statistics are available for download, can be requested by the authors, and/or only provide the lead single nucleotide nucleotide polymorphisms (SNPs). Other issues included ambiguous reports of sample size, genome build, or missing information such as standard errors and effect sizes. We were able to partly resolve these issues by using a custom pipeline, which automates tasks such as genome build liftover and SNP alignment and provides a uniform file format (Supplementary material online, ). This significantly decreased workload and improved downstream applications.

Genome-wide association studies summary statistics were selected by means of relevance, ‘latest and the largest’ and publicly available, easily accessible and/or open-source data was given priority. Incomplete summary statistics were omitted during this phase (Supplementary material online, ). A custom Python pipeline (see Data availability) was used to pre-process GWAS summary statistics and primary quality control. This primary control step was used to uniformly transform GWAS summary statistics because of a lack of a standardized format. In short, our custom pipeline was built to automate effect/other allele alignment with 1000 Genomes phase 3 release with a maximum frequency distance of 25%, maximum minor allele frequency (MAF) of 45% for ambivalent variants and lifting genomic positions to hg19 where needed. In preparation for GWAS alignment, dbSNP153 (GCF1405) was referenced to translate MEGASTROKE rs-ids to genetic positions and to augment AF, ASD, NICM, and TAGC (GWAS abbreviations are listed in ) with allele frequency information. Standard errors for HF and NICM were recovered as follows:

Where beta is the effect size, P is the P-value of association which is quantile normalized using the qnorm-function in R. SE is the standard error calculated by dividing the beta by Z. Effect sizes were calculated from odds ratios as beta = log(OR) for BIP, insomnia, IBD, and MDD before entering our custom pipeline.

The final selection was grouped into one of three categories: atherosclerotic disease and other cardiovascular disease (12), risk factors (14), and other (20). Data were further processed with MAGMA (version 1.07) for genome-wide analysis, annotation and characterization of significant hits via the FUMA web platform (version 1.3.6, accessed June 2020, https://fuma.ctglab.nl/). We applied genome-wide significance (P < 5 × 10−8) and linkage disequilibrium (r2 > 0.05) filtering for clumping of independent loci within 1000 kb of the lead variant, and included only variants with MAF >1%. For SVS, logOnset, EvrSmk and SVD no SNP associations below the genome-wide significance threshold we found. Therefore lead variant discovery was performed with a relaxed threshold of 5 × 10−6.

Selecting cell population-specific genes

Differentially expressed genes (DEGs) were selected via the Wilcoxon sum rank test. A gene was considered a DEG if (i) 10% of cells within a cluster expressed the gene with log2 fold change of at least 0.25; (ii) the gene passed the Bonferroni adjusted significance threshold of P < 0.05 for this test. Differentially expressed genes were determined via a ‘one cluster vs. one cluster’ and ‘one cluster vs. remaining clusters’ (Supplementary material online, ). This method results in genes specific for a cluster compared to every other single cluster and genes specific for a cluster compared to the whole single-cell reference set.

Gene score calculation

Fully processed summary statistics and cell population-specific DEG were used for calculating gene scores. The background gene set (n = 20 112 genes) was defined as all the genes that are reported and annotated by MAGMA, regardless of their respective P-value or Z score with duplicates excluded. We calculated the random gene score per cell-population by overlapping the genes in this background set with the cell-population-specific DEG. For every anaysis, the background genes are the same per GWAS. Z-score threshold for the P-values from MAGMA was set at 3 and we included a maximum of 150 candidate genes meeting this threshold (Supplementary material online, ). If there were more than 150 genes meeting the threshold, we only used the first 150. Calculations were made as follows: (i) per cell population, determine the overlap of the candidate genes of a GWAS with the cell population-specific DEG. (ii) Score each of these genes according to average cell population expression: assign one point for expression between 1 and 2, and two points for expression >2 and 0 points for expression <1. Only genes that are cell population-specific receive a score. Assign the points to the genes from the random candidate gene list in the same way (Supplementary material online, ). (iii) Compare the scores from the candidate genes with the scores of the random genes via one-tailed Wilcoxon Rank sum test in R and adjust the results per GWAS for multiple testing using the Benjamini–Hochberg method. Adjusted P-values smaller than 0.1 (10% false discovery rate) were considered significant (Supplementary material online, ). Enrichment was defined as the fold change of the mean of the gene scores divided over the mean of the gene scores in the random set. All scripts are available on GitHub (see Data availability).

Results

Single-cell RNA-seq of atherosclerotic plaques reveals 18 cell populations

We have expanded our previously published scRNA-seq datasets of human atherosclerotic plaque cells in a cohort which now includes 38 patients and identified 18 cell populations amongst 6191 cells. These expanded datasets allowed us to map the genetic loci not only to main cell types (e.g. SMCs) but also to more specific cell subtypes (e.g. synthetic, contractile or transitioning SMCs)—increasing resolution compared to the previous dataset. By combining known and recently published sets of markers we identified clusters of lymphocytes, macrophages, mast cells, ECs, and SMCs (). This data is available for interactive exploration via the online tool Plaqview accessible on www.plaqview.com. Notably, the expression of selected, widely studied GWAS targets for CAD (EDN1, HDAC9, NOS3, VEGFA, FN1, APOE, APOB, SMAD3, PCSK9, and LDLR) clearly demonstrates cell-specific patterns (). This underscores the need for a systematic analysis using large-scale GWAS summary statistics and motivated us to devise an intuitive two-stage approach.

Finding suitable therapeutic targets from genome-wide association studies with single-cell RNA sequencing. (A) UMAP visualization of 6191 cells divided into 18 cell populations derived from 38 human carotid artery plaques. (B) Heatmap of frequently studied coronary artery disease genome-wide association studies hits: EDN1, HDAC9, NOS3, VEGFA, FN1, APOE, APOB, SMAD3, PCSK9, and LDLR in atherosclerotic tissue shown in cells depicted in cell populations shown in A.

Gene-based analyses of atherosclerotic diseases and cell population-specific genes

First, we aimed to identify trait-associated genes using GWAS summary statistics by aggregating per-variant statistics and to calculate an empirically derived P-value for a given gene. Thus, gene-trait associations would be agnostically derived and solely based on genetic signals. To do this, we collected and partitioned GWAS summary statistics for 46 traits across 3 categories: atherosclerotic disease and other cardiovascular disease (12), risk factors (14), and other (20) (Supplementary material online, and ). From these GWASs, 40 were performed on population with European ancestry, and the remaining 6 had mixed or other ancestries. Next, to identify genes associated with individual traits we performed gene-based testing using MAGMA per individual GWAS. Genes mapped by gene-based testing are referred to as candidate genes (Supplementary material online, ).

Second, we aimed for an unbiased selection of genes that would best represent each cell population in a given dataset (Supplementary material online, ). Per tissue, these genes were selected through differential gene expression and include genes that are overexpressed in specific cell populations compared to all other cells (one vs. all) or compared to any individual cell cluster (one vs. one). This approach eliminates ubiquitously expressed genes and genes with low expression without the need of a predefined list (Supplementary material online, ) and identifies a similar number of cell-specific genes in each population. Cell population-specific genes will be further referred to as differentially expressed genes (DEGs). For cell populations in the atherosclerotic plaque, the average number of DEG was 461, with the lowest 307 for B Cells II and highest 678 DEGs for T Cells III.

Plaque cell populations show specific enrichment of genome-wide association studies candidate genes

To test our hypothesis that the genetic drivers for atherosclerotic disease can be found in plaque tissue, we applied our workflow on cardiovascular candidate genes and the DEG sets from individual cell types present in the atherosclerotic plaque (Supplementary material online, , , ). The signals we observed were independent of the total number of DEGs in each cell population (Supplementary material online, ), and did not show a bias towards GWAS sample size (Supplementary material online, ) or trait heritability (Supplementary material online, ). We found cell-specific enrichment for different GWAS in the atherosclerotic disease category, which was mainly observed in ECs and SMCs. Stroke candidate genes were enriched in EC populations, Regulatory T Cells, and B Cells II ().

Smooth muscle cells are predominantly associated with coronary artery disease genes

By highlighting the overlap between plaque-specific DEGs and candidate genes for CAD, we observed a significant enrichment in SMCs (P-adj. = 0.0012, 3.1 fold enrichment) and ECI (P-adj. = 0.0011, 2.7 fold enrichment) populations (). All P-values and fold enrichment values are listed in Supplementary material online, . This uncovered well-studied GWAS loci such as APOE, COL4A1, and COL4A2. The remainder of the overlap is driven by KANK2 and SORT1 in SMCs, SLC44A1 and ATP2B1 in ECI (). KANK2 belongs to the KANK protein family which plays a role in cytoskeletal function. SKI expression was present in both cell populations (). This proto-oncogene inhibits TGF-beta signalling. Interestingly, while using the single-cell transcriptomics enabled us to pinpoint specific cell (sub)populations, these potential targets exhibited ubiquitous expression in bulk transcriptomics data from GTEx tissue samples (). Notably, CD14+CD68+ Macrophages III showed an enrichment of 4.3 for CAD with genes APOE, SKI, and HNRNPUL1 but did not meet the significance threshold (P-adj. = 0.19).

Genome-wide association studies candidate genes show cell-specific enrichment in atherosclerotic plaque cell populations. (A) Plot representing overlaps of genome-wide association studies candidate genes on cardiovascular disease with genes specifically expressed in cell populations from human atherosclerotic plaque. Colour intensity indicates enrichment over random data and size indicates adjusted P-value (Benjamini–Hochberg). (B) Genes for coronary artery disease underlying the overlaps in smooth muscle cell and ECI populations. Only significant overlaps are plotted. (C) Dotplot depicting expression of genes from B in the respective cell populations. (D) Average of normalized expression (zero mean across samples) of genes from coronary artery disease in smooth muscle cells from A in 30 tissues from GTEx.

To verify our results we used an independent scRNA-seq dataset of carotid artery plaques by Pan et al. (Supplementary material online, ). The cell population labels broadly correspond between datasets (Supplementary material online, ). The patterns for enrichment are similar between datasets (Supplementary material online, ). We found that CAD, coronary artery calcification, and plaque presence candidate genes were enriched in the fibrochondrocyte cell population (), which corresponds to ACTA2+ SMCs cells in our dataset. For CAD, the enriched was driven by genes CDKN2A, COL4A1, COL4A2, APOC1, APOE, and CXCL12. Similar to the genes found in our datasets in SMCs and ECs (Supplementary material online, ).

Genome-wide association studies related to stroke show substrate in endothelial cells. (A) Plot representing overlaps of genome-wide association studies candidate genes on stroke with genes specifically expressed in cell populations from human atherosclerotic plaque. Colour intensity indicates enrichment over random data and size indicates adjusted P-value (Benjamini–Hochberg). (B) Heatmap depicting expression of genes from any stroke from A in the respective cell populations. Only significant overlaps are plotted, symbols represent the cell populations with the significant overlaps. (C) Average of normalized expression (zero mean across samples) of genes from any stroke in from B in 30 tissues from GTEx. (D) Genes for any ischaemic stroke underlying the overlaps in ECI and ECII populations. Only significant overlaps are plotted. (E) Genes for large artery stroke underlying the overlaps in T Reg and B cells II populations. Only significant overlaps are plotted.

Endothelial cells are enriched for ischaemic stroke-associated genes

Atherosclerotic disease is one of the major culprits underlying stroke, with carotid plaques being the main driver of large artery stroke. We subsequently investigated overlap of GWAS-associated genes for stroke subtypes any stroke (AS) and any ischaemic stroke (IS) and found they were predominantly enriched in ECI (AS: P-adj. = 0.011, 2.7 fold enrichment; IS: P-adj = 0.077, 2.5 fold enrichment) and ECII (AS: P-adj = 1.9 × 10−5, 3.9 fold enrichment; IS: P-adj = 0.023, 3.2 fold enrichment) subtypes (). Overlapping genes between AS and cell populations include ESAM, LMNA, and SLC44A2 (). Soluble ESAM (EC adhesion molecule) levels have been correlated to myocardial infarction and heart failure.Lmna deficient mouse models develop cardiovascular disease at an accelerated pace due to premature ageing. Of note, expression of these genes was not specific to vascular derived tissue of GTEx when using bulk transcriptomics data, suggesting that our findings are cell population and tissue specific (). No significant enrichment is found for IS in SMCs, but the cell population-specific signal in ECI and ECII is largely different from AS ().

Large artery stroke (LAS) is a subtype of IS. Strikingly, we see enrichment in CD3+ Regulatory T Cells and CD79A+ B Cells II (plasma cells) (P-adj = 0.026 and P-adj = 0.021, respectively) populations for LAS (). Suggesting the involvement of the adaptive immune system on the development of the LAS phenotype.

Transitional smooth muscle cells enriched for atherosclerotic disease candidate genes

Plaque ACTA2+ SMCs can be further divided into three sub-populations with a distinct phenotype; Synthetic SMCs, contractile SMCs, and a smaller cell population consisting of SMCs with an intermediate phenotype (transitional SMCs) (, Supplementary material online, and ). To explore if we can narrow down the signal from SMCs into a specific SMC subtype, we subjected this subset of plaque cells to our workflow (Supplementary material online, ). This revealed that the majority of the signal can be found in transitional SMCs. This population is predominantly enriched for CAD (P-adj. = 0.026, 2.4 fold enrichment), carotid intima-media thickness (cIMT) (P-adj. = 0.027, 2.4 fold enrichment), and plaque presence (P-adj. = 0.081, 2.3 fold enrichment) (). Intima-media thickness is a commonly used tool to measure atherosclerosis. Signal for LAS was present in transitional and synthetic SMCs with unique hits for each cell population ().

Figure 4

Analysis of three sub-populations of atherosclerotic smooth muscle cells. (A) UMAP visualisation of three sub-cell populations of smooth muscle cells in atherosclerotic plaque. (B) Dotplot depicting expression of synthetic and contractile smooth muscle cell markers in the three sub-populations. (C) Plot representing overlaps of genome-wide association studies candidate genes on cardiovascular disease and cardiometabolic traits with genes specifically expressed in the subset of smooth muscle cells from human atherosclerotic plaque. Colour intensity indicates enrichment over random data and size indicates adjusted P-value (Benjamini–Hochberg). (D) Plot depicting the number of unique and intersecting candidate genes from Large artery stroke enriched in transitional smooth muscle cells and synthetic smooth muscle cells.

To test the robustness of these results, we have repeated the analysis with the other public dataset (Pan et al.). For CAD, the overlap between transitional SMCs and corresponding ‘Fibrochondrocytes’ in Pan et al. (Supplementary material online, ) showed significant enrichment and again identified a largely overlapping set of commonly studied GWAS loci COL4A1, COL4A2, APOC1, APOE, and CXCL12.

Functional assessment of smooth muscle cell-associated candidates

To demonstrate the potential role of newly identified gene-cell pairs in atherosclerosis-relevant process we selected three genes enriched in the SMC population of plaques. SKI, KANK2, and EDNRA are CAD candidate genes and are currently understudied in relation to atherosclerotic disease in this cell population. Making them possible interesting targets for further functional studies. To investigate the role of these genes in SMCs, we functionally characterized SMCs from the ascending aortas in 151 healthy and diverse donors for cellular migration, proliferation, and calcification.KANK2 gene expression is correlated with proliferation under different stimuli (P-adj. < 0.1) and was negatively correlated with calcification in quiescent cells [Calcification Osteo: P-adj = 0.0058; Calcification (Pi): P-adj 0.069] (, Supplementary material online, ). Under all stimuli, the response was positive, except upon TGFβ1 exposure. SKI expression positively correlates with calcification in proliferative cells (P-adj = 0.069) but has a negative response towards TGFβ1-induced proliferation (P-adj = 0.069, R −0.21) (, Supplementary material online, ). EDNRA was negatively correlated with proliferation response to IL1β (P-adj. = 0.038) and calcification (P-adj = 0.062) in non-proliferative cells and showed a positive correlation in proliferative cells for migration (P-adj = 0.069) (, Supplementary material online, ).

Smooth muscle cell associated candidates influence migration and proliferation in vitro. (A) Correlation of KANK2, SKIi, and EDNRA expression non-proliferative and proliferative cells. (B) Plot depicting the Hi-C interactions around the main SNP associated with KANK2 (top) and regional association plot (bottom) in Hi-C data derived from IMR90 fibroblasts. AUC, area under curve; DHS, DNaseI hypersensitivity sites; Pi, inorganic phosphate; TAD, topologically associating domain.

Lipid-related risk factors are enriched in liver hepatocyte populations

Causal mechanisms for atherosclerosis have been linked to lipid metabolic processes active in the liver and blood glucose levels. This offered us the opportunity to assess the robustness of our approach across tissues and diseases. We acquired published scRNA-seq datasets of human liver cells and human pancreatic Langerhans islet cells to identify possible cell- and disease-specific targets.

Since one of the major liver functions is lipid metabolism, we expected strong signals for candidate genes involved in triglycerides levels. In line with our expectations, these GWASs indeed show the strongest signal in hepatocytes #3, #4, and #5 driven by several apolipoprotein genes (Supplementary material online, , ).

Langerhans islet cells overall showed weak enrichments for all atherosclerosis-related GWASs (Supplementary material online, ). Signal for CAD in endothelial and mesenchymal populations overlap with hits in comparable cell populations in plaque (Supplementary material online, ). We could detect significant overlap of type 2 diabetes-associated genes with beta cells (P-adj. = 0.001) driven by INS, WFS1, and FTO amongst others (Supplementary material online, ).

Altogether, we show that candidate genes for atherosclerotic disease are significantly enriched in EC and SMC populations of atherosclerotic plaque, notably in transitional SMCs. Stroke-related genes were enriched in ECs and lipid loci were found enriched in hepatocytes. Pancreatic beta cells were enriched for body mass index (BMI) candidate genes. Interestingly, only limited signals were observed in the cells of the adaptive immune system suggesting the genetic component for atherosclerosis is rooted in the structural cells of the vessel. We selected three candidate genes enriched in SMCs for functional testing in cells of ascending aortas for cellular calcification, proliferation and migration that showed that KANK2, SKI, and EDNRA expression correlated with calcification, migration, and proliferation in vascular SMCs.

Discussion

Post-GWAS analyses aiming to identify candidate genes for translation into clinical care are not straightforward. The thousands of variants associated with complex traits, including cardiovascular diseases, have made it abundantly clear that many genes contribute to disease risk and progression. For instance, GWAS for CAD have identified 163 loci linking to hundreds of candidate genes. Further efforts like chromatin interaction mapping linked almost 300 additional gene targets for atherosclerotic disease via distal chromatin interactions. This adds another layer to candidate gene prioritization. Selecting causal GWAS loci for functional testing can be deduced with vast knowledge of the trait or disease, but this remains speculative. There are numerous ways to accelerate the path to translational research. The intersection of loci found amongst different GWAS related to the same trait or disease can highlight some of the primary genes involved but, the cellular resolution is still lacking.,

Overlapping GWAS hits with pathways or within gene networks offers more context to study the genes. However, this is highly dependent on knowledge of the respective pathways and networks, and locating the relevant key players in the network is not indisputable.

Colocalization of tissue-specific eQTLs,, with GWAS data provides another alternative that can be used to estimate if there are overlapping causal variants in both datasets. For example, a recent study on type 2 diabetes leveraged publicly available eQTL data and a separate pancreatic dataset to find targets for glucose- and insulin-related trait loci. In concordance to our data, trait-associated genes are expressed in pancreatic islet tissue. However, cell population-specific signals are lacking and eQTL signals can be confounded by variability in tissue cell composition. Currently, the first efforts are being made to study eQTLs in a single-cell space, with single-cell quantification of gene expression in disease relevant tissue. With this method, it is possible to detangle the effect of SNPs on the single-cell and patient-specific level and construct personalized gene regulatory networks.,

Another straightforward approach to identify the causative genes is overlapping the expression profiles of different tissues and organs to GWAS candidate genes. Utilizing bulk RNA sequencing data applicable to the GWAS helps to find targets that are directly acting in the affected tissue. This powerful approach has one major drawback: the data are composed of multiple cell populations and this signal is difficult to deconvolute into individual cell populations. It is especially limiting when looking for cells that are present in almost all tissues—like ECs or SMCs. We show here that these bulk RNA tissues approaches can be misleading due to lack of specific expression patterns for relevant GWAS loci in specific tissue, like vascular wall tissue (). On the other hand, the overlap with single-cell transcriptomics can detect significant overlaps in distinct cell types such as ECs.

Here, we present a possible approach to the prioritization problem. We projected GWAS loci directly into single-cell transcriptomics datasets derived from disease-relevant tissue. This allowed us to identify both known and novel gene targets together with their associated cell populations. Notably, the CAD-associated gene SKI is significantly enriched in SMC and EC populations (). SKI inhibits TGF-beta signalling, which plays a major role in atherosclerotic disease progression by controlling processes such as cell proliferation and matrix formation. Consequently, SKI increases lesion size in LDR−/− mice. Another identified gene—KANK2 plays a role in cytoskeletal function by reducing the cell’s ability to migrate by activating TALIN. This gene is specifically expressed in the SMC population of plaques. Possible impairment of cell mobility by KANK2 could lead to inefficient SMC migration towards the cap resulting in an unstable plaque phenotype. Interestingly, our functional analysis of KANK2 showed a positive association between gene expression and cell migration in quiescent and proliferative cells of SMCs from ascending aortas (R = 0.15 and 0.2, respectively) (Supplementary material online, )—suggesting more complex mechanisms. Follow-up experiments creating gene-deficient primary cell lines derived from atherosclerotic plaques to study the molecular mechanistic of those genes can provide more information of its role in atherosclerosis. Notably, the CAD-associated locus in the vicinity of KANK2 also contains LDLR gene (, bottom). Even though LDLR involvement in atherosclerosis is widely accepted, we were able to show that KANK2 is another strong candidate in this region as it also shows a strong correlation with SMC proliferation (). This also demonstrates that our approach can identify genes that are located more further from the lead SNP but are positioned within same broader regulatory domain ()—as based on Hi-C data from IMR90 fibroblast cells (, top).

Currently, 32 stroke loci have been implicated by GWAS. From the new main stroke loci found by Malik et al., only SLC44A2 was significantly enriched in scRNA-seq populations. However, we identified other potential targets in SMCs, ECs and in B cells, and regulatory T cell populations (). One of the targets, ESAM, was not yet directly linked to atherosclerosis but is known to enhance vascular permeability in selective tissues and could therefore contribute to the influx of cholesterol and inflammatory cells in the arteries contributing to the stroke phenotype. Another target, Gs-coupled receptor calcitonin receptor-like receptor (CALCRL) is expressed in both ECI and ECII populations. Activation of this receptor leads to an anti-inflammatory response in ECs. Subsequently, CALCRL deficient mice present with more advanced atherosclerotic lesions. A recent study by Örd et al. researched the enrichment of CAD, IS, and seven other cardiometabolic GWAS signals in scATAC-seq peaks of carotid artery plaques. Similarly, they report enrichment of CAD and stroke loci in SMCs and ECs and total cholesterol loci in macrophages.

Remarkable is the absence of significant enrichment in the cells of the adaptive immune system of the plaque for CAD. Out of six T-cell populations, four exhibited the enrichment of 0.00 and CD3+CD8A+ T cells II and II showed enrichment of 1.3 and 0.37, respectively, indicating that there is no or very little overlap between CAD candidate genes and cell population-specific DEGs. This suggests that the genetic component of atherosclerotic disease is not rooted in the adaptive immune system but rather in the structural cells of the vessel and influences the integrity of the vascular wall or plaque stability.

Our workflow can be generalized to other single-cell transcriptomics datasets and complex traits. We applied it to scRNA-seq data from tissues with expected causal involvement in atherosclerotic disease: human liver cells and human pancreatic Langerhans Islet cells. Genome-wide association studies studying lipid levels were associated with hepatocyte subsets but in particular subsets Hep3, Hep4, and Hep5 (Supplementary material online, ). It suggests that these risk factors for atherosclerosis have a stronger biological foundation in lipid metabolism in the liver rather than local lipid metabolism in the plaque. Langerhans Islet cells show less prominent signals for our selected GWAS, suggesting that the genetic mechanisms involved in atherosclerosis do not have a local component in the pancreas (Supplementary material online, ).

Little to no signal can be observed for GWASs that have no relation to the cells composing the atherosclerotic plaque, such as height, cancers and congnition-related traits such as bipolar disorder and neuroticism (Supplementary material online, ). Yet, these GWASs have ample candidate genes and considerable GWAS sample size—suggesting the specificity of our approach. (Supplementary material online, , ). Similarly, no signal for BMI is found in plaque cell populations, intuitively confirming that the mechanisms underlying obesity cannot be accounted for in plaque tissue.

In highly heterogeneous tissues, such as atherosclerotic plaque, multiple genes can be shared amongst the same class of cell populations (e.g. all ECs), whilst they are lacking in others. To resolve this problem and to keep such genes as cell type-specific, we implemented an approach that can unambiguously uncover cell-specific genes without the need of a predefined list of ubiquitously expressed genes. Although this method is inclusive, the balance between the sufficient number of DEGs and their specificity can be adapted depending on the scRNA-seq dataset in question.

Our approach, however, is subject to several limitations. In our workflow, a bias may be introduced by low-powered GWAS, or by threshold introduced in the selection of candidate genes. The depth and quality of the scRNA sequencing is a another limiting factor. The amount and kind of RNA that can be detected and favours the most abundant genes with highest mRNA counts per cell. As a result certain gene classes, like transcription factors or signalling molecules are more difficult to detect. At the same time, cell types like neutrophils and foam cells in atherosclerotic plaque are difficult to capture by single-cell transcriptomics and might be underrepresented or completely missed. Digestion protocols, tissue quality, and sample handling also affect the quality and/or cell composition of scRNA-seq data. Finally, our workflow does not utilize the existing information about eQTLs, non-coding DNA regulatory units and chromatin organization. Integration of these datasets has potential to increase the accuracy and number of identified target genes. However, we aimed to remain agnostic with respect to identifying candidate genes based on genetic data, and thus designed our experiment using GWAS data only. For instance, integration with Hi-C data would be biased towards a specific cell type in a specific tissue, and favour the selection of genes expression in that given tissue and thus bias the cell type enrichment analysis. Our approach moves atherosclerotic candidate genes forward based on the gene-based genetic association and the differential expression plaques.

In conclusion, we systematically projected GWAS summary statistics for 46 traits directly into single-cell transcriptomic profiles across 18 distinct cell populations derived from plaques. We devised a two-step workflow that (i) leverages large-scale human genetic data to prioritize trait-associated genes, and (ii) calculates a cell-type-specific enrichment score. We uncovered 11 potential cell population-specific targets for CAD and provided functional evidence for KANK2, SKI, and EDNRA in SMCs. Additionally, we verified the robustness and transportability of our approach by confirming the tissue- and trait-specific enrichment of circulating lipid genes in liver and glycaemic disorder genes in pancreatic tissue. The principal strength of our framework is the cell-specific resolution it offers for disease-associated genes mapped to GWAS loci, thus putting the emphasis on prioritization of candidates (instead of adding to an ever-growing candidate-list) and accelerating the path to functional testing.

Lead author biography

Lotte Slenders earned a BS in Biology and a MS in Molecular and Cellular Life Sciences from Utrecht University in the Netherlands. In 2018, she became a PhD candidate in the lab of Gerard Pasterkamp at the Utrecht University Medical Centre. During her PhD, she further developed as a data-scientist with a biology background. Her main research topic revolves around the omics data of human atherosclerotic plaques from the AtheroExpress biobank study. She is using single-cell RNAseq data to study complex genetic loci and molecular mechanisms of cell plasticity such as endoMT.

Supplementary material

Supplementary material is available at European Heart Journal Open online.

Click here for additional data file.