Atherosclerotic plaque epigenetic age acceleration is characterized by mesenchymal reprogramming and poor prognosis

Study population

The Athero-Express is an ongoing, prospective biobank study, collecting atherosclerotic plaques from patients undergoing carotid endarterectomy (CEA) in two Dutch tertiary referral hospitals: University Medical Centre Utrecht and St Antonius Hospital, Nieuwegein. Study design, inclusion criteria have been published previously¹⁹. A baseline table for the 485 patients (148 women and 337 men) in this study is provided in Supplemental Dataset 1. Indication for surgery was based on international guidelines for carotid and iliofemoral atherosclerotic disease^20–23, and standardized treatment protocols and operative techniques were applied. The medical ethics committees in both participating centres approved the study. All patients provided written informed consent. A composite cardiovascular endpoint for the outcome analysis of patients who underwent CEA was used. This consisted of: (sudden) cardiovascular death, stroke, myocardial infarction, coronary intervention (coronary artery bypass grafting or percutaneous coronary intervention), peripheral (re)intervention or leg amputation. For patients who reached multiple endpoints during follow-up, only the first manifestation of a cardiovascular event was used for analysis of the composite endpoint.

Plaque histopathology

As described previously^19,24,25, the atherosclerotic plaque was processed directly after surgery and (immune-)histochemical staining was routinely performed on the culprit lesion (segment with the highest plaque burden) for identification of macrophages (CD 68), calcification (haematoxylin-eosin (HE)), smooth muscle cells (alfa actin), collagen (Picro Sirius red (PSR)), intra-plaque hemorrhage (HE, Elastin von Gieson staining), vessel density (CD34) and fat (PSR, HE).

Bioinformatics

We used DNA methylation data from 485 atherosclerotic plaques and 93 overlapping blood samples, for which histological procedures, data collection, and normalization procedures were previously published²⁶. In brief, DNA was extracted from stored plaque segments and stored blood samples of patients using standardized in-house protocols as described before in Van der Laan et al²⁷. DNA purity and concentration were assessed using the Nanodrop 1000 system (Thermo Scientific, Massachusetts, USA). DNA concentrations were equalized at 600 ng, randomized over 96-well plates and bisulfite converted using a cycling protocol, and the EZ-96 DNA methylation kit (Zymo Research, Orange County, USA). Subsequently, DNA methylation was measured on the Infinium HumanMethylation450 Beadchip Array (HM450k, Illumina, San Diego, USA), which was performed at the Erasmus Medical Center Human Genotyping Facility in Rotterdam, the Netherlands. Processing of the sample and array was performed according to the manufacturer’s protocol. We calculated epigenetic age using Horvath’s clock with the R package cgageR (v0.1.0, accessed from: https://github.com/metamaden/cgageR), which uses methylation values at 353 CpGs to generate a continuous value mimicking a person’s epigenetic age.

The dimensionality reduction of 485 patient samples and B-values of the DNA methylation data at 483,731 CpGs was performed using a t-distributed stochastic neighbor embedding with the rtsne package (v0.15). Differential methylation analysis was performed using limma (v3.46.0)³². We called a CpG with a Bonferroni-corrected p-value <0.05 significantly differentially methylated. We used ChIPseeker³³ (v.1.26.2) to annotate the differentially methylated CpGs. Differentially methylated promoter CpGs were enriched using clusterProfiler³⁰.

In 420 of these plaques we were able to collect methylation and RNA sequencing data which was processed accpring to previously published methods²⁸. Briefly, the raw RNAseq read counts were corrected for UMI sampling: then normalized by sample sequencing depth and log-transformed. Differential gene expression analysis was performed using DESeq2²⁹ (v1.30.1). A gene was called differentially expressed if FDR-adjusted p-value < 0.1. Gene enrichment analysis were performed using Clusterprofiler v(3.18.1)³⁰, enrichplot (v1.10.2), and enrichr³¹.

Statistical analysis

All analyses were performed using R (v4.0.4) in R-Studio. Baseline tables were generated using the tableone package (accessed from: https://github.com/kaz-yos/tableone). We calculated the EAA as the residual from a linear model: Epigenetic Age ∼ Chronological Age. This removes the influence of chronological age on epigenetic age, since both are heavily correlated. The same procedure also means that positive EAA values indicate accelerated ageing, whereas negative EAA values indicate decelerated ageing. Linear models were constructed in R adjusting for potential confounders when needed. Cox regressions, plotting, and adjusted cox regressions were performed using rms (v6.2.0), survminer (v0.4.9) and survival (v3.2.10). Code can be found online at www.github.com/

Single-cell RNA-sequencing

We used previously published single-cell RNA-sequencing data of human atherosclerotic plaques, as well as new data from the same cohort³⁴. We were able to distinguish 20 different cell clusters from 46 patients. We used the addModuleScore within seurat³⁵ to calculate module scores for genes higher expressed within accelerated plaques, and created violin plots to plot the obtained module scores over the different clusters (Fig. 4C).

Chronological age and biological plaque age

To determine the epigenetic age of the atherosclerotic plaque, we used DNA methylation data of plaques from 485 patients (148 women and 337 men) who underwent a carotid endarterectomy. A baseline table with the clinical characteristics of the 485 patients is provided in the Supplemental Dataset 1. Chronological age for the patients ranged from 40 to 91 years in females, and 43 to 88 years in males, with a median age of 69 for both sexes. We used the Horvath pan tissue clock since it applies to all nucleated cells including those found in plaques. Epigenetic age measured using the Horvath pan tissue clock¹ ranged from 50.4 years of age for females to 81.3 years of age, while it ranged from 51.0 to 88.8 years of age for males, with a median for females of 65.0 years, and 67.6 years of males (Fig. 1A). Chronological age of the patient and epigenetic age of atherosclerotic plaques were strongly correlated (r = 0.695, p-value = 3.5e-71, Fig. 1B), indicating that the Horvath clock is relatively accurate in predicting age in diseased atherosclerotic tissue. We calculated epigenetic age acceleration (EAA) by taking the residuals of the linear model in which epigenetic age was predicted by chronological age. It is a mathematical fact that the resulting measure of age acceleration is not correlated with age (correlation zero). A negative EAA means that the chronological age is higher than the epigenetic age of the sample under study, which points to decelerated ageing of the atherosclerotic plaque, whereas a positive EAA indicates accelerated ageing of the plaque. EAA ranged from -15.4 years to 12.6 years in females, and from -10.6 years to 14.4 years in males (Fig. 1C). The median EAA significantly differed between females and males by approximately 2 years (females: -2.2, IQR [-4.3, 2.2] years vs. males: 0.3, IQR [-2.9, 3.8] years; p-value: 5.99e-5).

• Download figure
• Open in new tab

Figure 1. Chronological age versus plaque epigenetic age in 485 carotid endarterectomy patients.

A) A boxplot shows the distribution and range of chronological age for the 485 patients in the study, stratified by sex. B) A boxplot shows the distribution and range of epigenetic age determined from DNA of 485 atherosclerotic plaques using the Horvath clock, stratified by sex. C) A scatterplot shows the strong correlation between chronological age (x-axis) and epigenetic age (y-axis), dots are color based on sex. D) A boxplot shows the distribution of range of epigenetic age acceleration (EAA), sex-stratified.

Epigenetic age acceleration is associated with secondary cardiovascular events

Patients within the Athero-Express biobank are followed up for 3 years for secondary events (see Methods). First, we performed a standard Cox-regression on EAA and secondary events (Fig. 2). To identify any risk factors associated with EAA, we categorized our patient population in two groups, those with accelerated ageing (EAA > 0, n = 224) and those with decelerated ageing (EAA < 0, n = 261) (Table 1). Importantly, since EAA is defined as the residual of the linear model specified above, chronological age is not different between the two groups. Notable differences between the groups are the number of males and females (accelerated: 80% male, decelerated: 60% male, p-value < 0.001), their body-mass index (BMI) (accelerated: 27.1 ± 4.2 SD, decelerated: 26.1 ± 3.6 SD, p-value = 0.004), and the presence of diabetes (accelerated: 26.8% diabetic, decelerated: 18.4% diabetic, p-value = 0.0035). No other risk factors (plasma low-density lipoprotein (LDL), plasma high-density lipoprotein (HDL), hypertension, smoking and high sensitivity C-reactive protein (C-RP)) significantly differed between patients with accelerated and decelerated plaque ageing.

• Download figure
• Open in new tab

Figure 2. Epigenetic age acceleration of carotid plaques and follow-up in 485 male and female carotid endarterectomy patients.

A) A line-graph portrays the relationship between the log of the relative hazard and EAA. The line intersects the y-axis at EAA = 1.33e-16, indicated by the dashed purple line. Dashed grey lines indicate the confidence interval. B) Cox-regression survival graphs, adjusted for sex, BMI, and diabetes, are shown for accelerated (purple) and decelerated (green) plaques.

To investigate prognosis in the two EAA-groups we performed a Cox-regression, which showed that being in the group with accelerated plaque ageing was detrimental for secondary outcome (HR = 1.7, p-value = 0.0034). This result did not substantially change when we corrected for sex, BMI, and diabetes (adjusted Cox-graphs Fig. 2B, HR = 1.56, p-value = 0.02). Adjustment by blood EAA (see Methods) in a subset of patients (n = 93) revealed this to be a plaque-specific effect (HR plaque EAA = 1.3, p-value = 0.018; HR blood EAA = 0.1, p-value = 0.81), furthermore blood EAA was not predictive for outcome in a separate model (n = 93, HR = 1.6, p-value = 0.44, Supplemental Fig. 1).

Interestingly, we observed a moderate correlation between plaque EAA and blood-based EAA across the subset of patients for whom both measurements were available (correlation coefficient = 0.44, p-value < 0.001, n = 93). The moderate correlation illustrates that blood is only a poor surrogate for plaque tissue. There was no significant interaction between EAA and sex (p-value = 0.88). As prognosis in these patients is affected by plaque phenotypes⁹, we studied if EAA was specifically higher in plaques with plaque hemorrhage. However, plaque hemorrhage was not related with EAA, and adjustment for plaque hemorrhage of the Cox-regression showed independent association with secondary outcome for both plaque hemorrhage (HR = 1.6, p-value = 0.02) and EAA (HR = 1.7, p-value = 0.007).

Accelerated ageing is independent of plaque characteristics

As plaque accelerated ageing independently predicts cardiovascular secondary outcome, we tested if this coincided with a certain plaque phenotype. The group with accelerated ageing did not present a different profile of semi-quantitative hallmarks of plaque histology: fat, plaque hemorrhage, collagen, calcifications, smooth muscle cell content and macrophage content (see Methods for staining procedures, Supplemental Fig. 2). Plaque phenotype (as shown in Table 1) can be seen as a proxy for plaque cellular composition which can significantly affect epigenetic age calculations². Therefore, we investigated if adjusting the Cox-regression model for secondary cardiovascular outcome by plaque phenotype changed outcome. All analyses remained unchanged and significant, independently of plaque phenotype (HR = 1.6, p-value = 0.021). This suggests that EAA predicts secondary cardiovascular outcome independent of histological plaque composition. This prompted us to investigate the effect of EAA on plaque biology at a molecular level by using genome wide epigenetic signatures (beyond the 353 CpGs on which the EAA clock was built on). In addition, we investigated gene expression patterns in accelerated and decelerated plaques.

Epigenetic age acceleration drives the genome-wide epigenetic state in human plaques

EAA has previously been described as a proxy for the state of the epigenetic maintenance system¹. An accelerated ageing in our case would indicate a that the epigenetic machinery in the plaque has been under an increased biological workload. In contrast, a decelerated plaque ageing would indicate that the epigenetic machinery has not been under a strong epigenetic maintenance influence as would be expected from plaques which are biologically less active¹⁰. We therefore studied DNA-methylation patterns in accelerated and decelerated aged plaques (n=485 male and female plaques). For this, we used dimensionality reduction on the DNA methylation data of 483,731 CpGs. We excluded the Horvath-clock CpGs to exclude ageing-related methylation patterns (Fig. 3A). We observed two large clusters representing male and female plaques in our cohort (Fig. 3A). This is likely driven by the differences in DNA methylation of the sex chromosomes: X inactivation leads to hypermethylation of many X chromosomal markers in females ^11,12. Decelerated and accelerated aged plaque clustered separately within the female, and to some extent, the male plaque groups (Fig. 3A). This suggests that besides sex, EAA is one of the largest drivers of DNA methylation in plaques. Taking male and female plaques together, we found 65,488 differentially methylated CpGs between decelerated and accelerated aged plaques (Supplemental Dataset 2). EAA was associated with hypermethylation of CpGs (45,287 out of 65,488, Fig. 3B). By using gene enrichment analysis on genes marked with differentially methylated CpGs in their promoter (within 1000 bases of the transcription start site), we found enrichment for smooth muscle cell contraction and extracellular matrix organization. We found immune activation and cell adhesion enrichment for hypomethylated promoters of EAA plaques (Fig. 3C). These enrichments were consistent in both females and males when stratifying by sex (Supplemental Fig. 3). We also studied enrichment of atherosclerotic pathways within the Horvath clock CpGs, but no signal was found for immune cell activation and smooth muscle cell pathways. This is in line with previous results that highlighted the processes enriched within the clock CpGs which point to cell death/survival, cellular growth/proliferation, organismal/tissue development, and cancer².

• Download figure
• Open in new tab

Figure 3. Epigenetic age acceleration and DNA methylome of 485 male and female carotid plaques.

A) A tSNE-plot is shown, reducing data of 485 samples based on project 483731 CpGs into 2 dimensions. Shape of the data point indicates sex, color indicates EAA group. We find visual evidence that females cluster into 2 separate groups that can be distinguished by EAA. B) A density plot shows the distribution of log2 fold changes of all significant differentially methylated CpGs between accelerated and decelerated plaques. C) A dot-plot shows the gene ontology Biological Processes enrichment of promoter CpGs significantly hyper- and hypomethylated with accelerated epigenetic ageing corrected for hospital of inclusion and sex. Color indicates significance of the enrichment, whereas size of the dot indicates gene ratio.

Smooth muscle cells and endothelial cells show transcriptional changes in epigenetically accelerated plaques

To understand which molecular processes are active during plaque EAA, we studied differential gene expression between accelerated and decelerated plaques using whole plaque RNAseq. Of the 19,291 genes tested, 3,031 were differentially expressed (FDR < 0.1) between decelerated and accelerated plaques (Supplemental Dataset 3). The majority of the differentially expressed genes were higher expressed in accelerated plaques (2839 out of 3031) as compared to decelerated plaques (Fig. 4A). We found upregulated genes encoding for ribosomal proteins, as well as mitochondrial genes, in accelerated plaques. This points towards a general increase in translation turn over and mitochondrial respiratory activity within the plaque (Fig. 4B). Gene enrichment analysis using other established databases (i.e. Reactome, KEGG) on the same genes also highlighted TGFβ signaling, endothelial to mesenchymal transition and extracellular matrix-related processes (Supplemental Dataset 4) which point towards smooth muscle cells.

• Download figure
• Open in new tab

Figure 4. Epigenetic age acceleration and RNA transcriptome of 420 male and female carotid plaques.

A) An enhanced volcano-plot is shown for the differential expression analysis between accelerated and decelerated plaques, where the x-axis shows the –log10 of the p-value and the y-axis shows the log 2 fold change (accelerated over decelerated ageing). Notable examples are highlighted in the plot by gene symbol. B) An enrichment map is shown for the gene ontology Biological Processes enrichment performed on genes higher expressed with accelerated ageing. Nodes indicate a gene set, whereas edges indicate overlap between the gene sets. The vicinity of nodes indicates similarity between gene sets. Color highlights significance, whereas size of the nodes indicates the number of genes within that node. C) A violin plot shows the module score of genes higher expressed with accelerated ageing from the bulk RNA-sequencing over the single-cell clusters detected from the single-cell RNA-sequencing of the human atherosclerotic plaque. Colors indicate different cell clusters.

Our transcriptomic and epigenetic analyses of EAA suggest that gene regulation and translation is altered in plaque smooth muscle cells. With the knowledge that smooth muscle cells display remarkable plasticity in atherosclerotic disease, we hypothesized that plaque EAA might alter the transcriptional state of these smooth muscle cells. Therefore, we combined our whole plaque RNA dataset (420 patients, 127 women and 293 men) with single cell RNA data to project the plaque EAA genes over 20 single cell clusters of the plaque including 17 cluster of immune cells (macrophages, T-cells, B-cells, natural killer cells and mast cells) and 3 clusters of mesenchymal cells (ACTA2⁺ smooth muscle cells and CD34⁺ endothelial cells) within the plaque (4,984 total cells from 46 patients, 20 women and 26 men)¹³. We generated a module score of the genes that were higher expressed in accelerated compared to decelerated plaques and projected this module score onto the individual cells within the 20 identified plaque cell clusters (Fig. 4C). This module score was significantly higher in ACTA2⁺ smooth muscle cells (mean module score: 0.15) and the CD34⁺ endothelial cell clusters I and II (mean module scores: 0.12 and 0.12, respectively) compared to the rest of the cells clusters within the plaque (mean module score: 0.02; all p-values <0.001) suggesting that mesenchymal cells are most sensitive to accelerated ageing.