Epigenome-Guided Analysis of the Transcriptome of Plaque Macrophages during Atherosclerosis Regression Reveals Activation of the Wnt Signaling Pathway

Summary

Atherosclerosis, a progressive accumulation of lipid-rich plaque within arteries, is an inflammatory disease in which the response of macrophages (a key cell type of the innate immune system) to plasma lipoproteins plays a central role. In humans, the goal of significantly reducing already-established plaque through drug treatments, including statins, remains elusive. In mice, atherosclerosis can be reversed by experimental manipulations that lower circulating lipid levels. A common feature of many regression models is that macrophages transition to a less inflammatory state and emigrate from the plaque. While the molecular regulators that control these responses are largely unknown, we hypothesized that by integrating global measurements of macrophage gene expression in regressing plaques with measurements of the macrophage chromatin landscape, we could identify key molecules that control macrophage responses to the lowering of circulating lipid levels. Our systems biology analysis of plaque macrophages yielded a network in which the Wnt signaling pathway emerged as a candidate upstream regulator. Wnt signaling is known to affect both inflammation and the ability of macrophages to migrate from one location to another, and our targeted validation studies provide evidence that Wnt signaling is increased in plaque macrophages during regression. Our findings both demonstrate the power of a systems approach to uncover candidate regulators of regression and to identify a potential new therapeutic target.

Introduction

Underlying many cardiovascular diseases is atherosclerosis, an arterial accumulation of lipid-laden macrophages (foam cells) with attendant chronic inflammation [1]. Atherosclerosis is endemic in adults [2] and can progress over decades before presenting with overt symptoms [3]. Atherosclerosis is correlated with risk factors such as hyperlipidemia, especially high plasma levels of low-density lipoprotein (LDL) cholesterol [4], [5]. There is significant interest in finding new therapies to promote the regression of atherosclerotic plaques [6]–[8], with particular emphasis on the plaque macrophage as a potential target because of its roles in lipoprotein uptake and vascular inflammation [9]–[12], plaque destabilization [13], [14], and plaque remodeling in response to dietary changes [15].

In mice, plaque regression has been achieved through expression [16]–[18] or inactivation [19], [20] of lipid metabolism-modifying genes; aortic tissue grafting [21]–[23]; and treatment with small-molecule therapeutics [24]–[26] or biologics [27], [28], frequently in conjunction with a diet alteration. We have been studying the role of macrophages during plaque regression in two mouse models of regression, the Reversa mouse [19], [20], [29], [30] and an aortic arch transplant model [21], [22], [31]. Reversa mice are Ldlr^−/− Apob^100/100 [32], [33] and they possess a loxP-flanked allele of the gene encoding microsomal triglyceride transfer protein (Mttp) and a drug-inducible Cre recombinase (Cre) transgene [19], [34]. In the Reversa mouse, hyperlipidemia can be reversed by the transient induction of Cre expression in the liver, which causes recombinant knockout of Mttp, and thus significantly lowers the plasma levels of very low density lipoprotein (VLDL) and LDL [19]. In prior work we have shown that inactivation of Mttp in the Reversa mouse, when combined with a switch from Western diet to chow, leads to plaque regression at the aortic root [20]. In the aortic transplant model, an atherosclerotic aortic arch from an Apoe^−/− mouse [35] is grafted in place of a segment of the abdominal aorta [31] of either a wild-type (WT) or an Apoe^−/− mouse. Grafts within WT recipients exhibit plaque regression, while plaques within Apoe^−/ recipients continue to progress [22]. A common feature of the Reversa and transplant regression models is the substantial reduction in plaque macrophages during regression [20], [23], [29], [30], [36]. It is at present unknown whether there are molecular pathways that may initiate the depletion of plaque macrophages that are shared across regression models. Identifying such core pathways could yield new therapeutic approaches to stimulate regression and/or beneficial remodeling of plaque.

Recently, we studied transcriptome differences between macrophages in regressing versus progressing plaques in the aortic transplant model [36]. We observed patterns of differential expression that suggest that macrophages in regressing plaques up-regulate genes associated with cell motility, down-regulate genes associated with cell adhesion, and up-regulate genes associated with an anti-inflammatory, “M2 macrophage” [37] phenotype. Although our analysis pointed to several candidate molecular regulators, it is unknown whether there are transcription factors that are common to the transplant and Reversa regression models that act as master controllers for the responses of plaque macrophages to lipid lowering.

In the current study we carried out a systems investigation of the Reversa model, with three goals: (i) to identify pathways and gene functions that are associated with plaque macrophage responses to lipid lowering in vivo; (ii) to understand the connections between gene sets identified in our study and gene sets from previous studies of macrophage foam cell formation, aortic wall changes during plaque progression and regression, and atherosclerosis risk alleles in humans; and (iii) to identify transcriptional regulators that are associated with macrophage responses to lipid lowering in vivo. We therefore carried out microarray transcriptome profiling of plaque macrophages isolated from the aortic root of Reversa mice undergoing lipid lowering in vivo and analyzed differentially expressed genes using multiple bioinformatics approaches. Gene functional enrichment analysis (i) detected over-representation of cytoskeletal-binding and Rho GTPase genes among genes that are upregulated in response to lipid lowering, pointing to cytoskeletal reorganization during regression. Comparative analysis (ii) of gene sets from our study with those from previous studies of in vitro models of macrophage foam cell formation revealed a statistically significant overlap. Additionally, we found a significant overlap between human monocyte expression quantitative trait loci (eQTLs) for human orthologs of genes that are differentially expressed during plaque regression, and genetic loci that are associated with risk of atherosclerosis or coronary artery disease. Promoters of genes identified in our study were analyzed for transcription factor binding site over-representation (iii) using a novel chromatin-guided method, REMINISCE (Figure S1 and Methods). REMINISCE leverages macrophage epigenomic and chromatin measurements, such as histone acetylation [38], [39] and DNase I hypersensitive sites [40], in the analysis of the 5′ regulatory regions of differentially expressed genes, including enhancers as well as promoters. In both the Reversa and aortic transplant regression models, we found that the consensus binding site sequence for the T-cell specific, HMG-box factors (TCF) and lymphoid enhancer factors (LEF), together known as the TCF/LEF family of transcription factors, was overrepresented within the 5′ regulatory regions of genes that are upregulated in plaque macrophages during regression.

TCF/LEF transcription factors are activated by nuclear β-catenin (CTNNB1) that accumulates in response to activation of the canonical Wnt signaling pathway [41], [42]. The Wnt pathway controls multiple functions in development, tissue organization, cell proliferation, cell migration [42], and inflammation [43]. In macrophages, canonical Wnt signaling is thought to promote cell motility through a β-catenin-dependent mechanism [44]. In this study we show, for the first time, that activation of canonical Wnt signaling, defined by up-regulation of β-catenin, occurs within macrophage-rich plaque in two mechanistically distinct lipid-lowering models of plaque regression, the Reversa and aortic transplant models.

Results

Identification of pathways and gene functions that are associated with plaque macrophage responses to lipid lowering in vivo

Reversa mice were maintained on Western diet for 16 weeks, which allowed for the development of hypercholesterolemia (903±71 mg/dL total cholesterol [mean ± standard error, SE]) and substantial, macrophage-rich plaques at the aortic root, as determined by immunostaining with the macrophage marker Cluster of Differentiation 68 (CD68). At 16 weeks, a group of animals were sacrificed to measure the baseline plasma total cholesterol levels and CD68+ areas within plaque within aortic root sections. The remaining animals were switched to chow diet and divided into two groups, an experimental group that received four injections of polyinosinic:polycytidylic acid (poly I:C) to induce Mx1:Cre in order to inactivate Mttp^fl/fl, and a control group that was injected with vehicle (saline) on the same schedule. On day seven or 14 after the final injection, animals were sacrificed and plasma cholesterol and plaque CD68+ areas were measured. By day 14, cholesterol levels had decreased by 88% in the Mttp-inactivated animals vs. baseline (Figure 1A and Table 1, P<0.001) and plaque CD68+ area had decreased by 36% in the Mttp-inactivated animals vs. baseline (Figure 1B and Table 2, P<0.05). Additionally, plaque CD68+ area was reduced in Mttp-inactivated vs. vehicle-treated animals at seven and 14 days (P<0.05); at day 14, CD68+ area was 26% smaller in the Mttp-inactivated animals to the vehicle-treated animals (Figure 1C). These results indicated that there is plaque regression that is attributable to inactivation of Mttp. In contrast, comparisons of the saline-treated sample groups to baseline showed no statistically significant differences in CD68+ area that could be attributed to diet switch alone.

In order to focus on the molecular signaling events that initiate the reduction in plaque macrophages, we selected an early time point (day seven) when there was a statistically significant difference in the CD68+ area in the treatment group vs. vehicle (Figure 1, B and C). CD68+ plaque macrophages were isolated from aortic root sections from Mttp-inactivated Reversa, vehicle-treated Reversa, and poly I:C-injected Ldlr^−/− animals (Table S1) by laser capture microdissection (LCM) [45], [46], a procedure that we have previously demonstrated specifically enriches for the mRNAs of macrophage-related genes by more than 30-fold over levels seen in aortic homogenates [45]. Global transcriptome profiles were then obtained from the LCM-derived samples using high-density oligonucleotide microarray hybridization. A total of 213 genes were detected as differentially expressed in plaque macrophages from Mttp-inactivated vs. vehicle-treated animals (93 had higher expression levels in Mttp-inactivated than in control, and 120 had lower levels; Table S2). In order to control for any effect of the poly I:C injections on the macrophage transcriptional profile we compared the transcriptomes of plaque CD68+ cells from saline-treated Reversa mice with those from poly I:C-treated Ldlr^−/− mice. We found that only one of the 213 lipid-responsive genes, Trnt1, was detected as differentially expressed in the Reversa saline vs. Ldlr^−/− poly I:C comparison, and this gene was excluded from further bioinformatic analysis. To further substantiate that the 213 genes in plaque CD68+ cells are responding to lipid lowering rather than directly to poly I:C, the list of 213 genes was compared with a list of genes that are differentially expressed (vs. vehicle) in bone marrow-derived macrophages stimulated in vitro with poly I:C (Table S3) to determine whether the degree of overlap in the two gene lists is greater than would be expected by chance. The analysis revealed that there is not more overlap with in vitro poly I:C responses than would be expected by chance (P>0.1, Fisher's Exact Test).

To gain insight into the biological processes that occur during lipid lowering in vivo, we performed gene functional annotation enrichment analysis. The most significantly overrepresented processes associated with the upregulated genes are related to cell motility, cytoskeletal binding, and Rho GTPases (Figure 2A and Table S4). For the genes that were downregulated in macrophages from lipid-lowered vs. control animals, the gene annotation term that was the most significantly overrepresented was mitochondrion (Figure 2B and Table S4). In order to determine whether lipid lowering in vivo might have increased plaque macrophage apoptosis, multiple gene functional enrichment analysis tools were applied, and in each case, no finding of a statistical enrichment for apoptosis-related gene annotations was detected (see Methods).

Contextualization of lipid lowering-responsive genes with genes identified in previous studies and atherosclerosis risk alleles in humans

We hypothesized that among the plaque macrophage genes that are responsive to lipid lowering in vivo, genes that had been previously reported to be differentially expressed in response to lipid loading in vitro would be statistically overrepresented. We therefore combined differentially expressed gene lists from 21 datasets from 15 studies in the atherosclerosis and macrophage foam cell literature (Table S5), and scored genes based on the number of studies in which they were detected as differentially expressed during foam cell formation or disease progression (similar to the approach that was used in a recent meta-analysis to identify obesity-related genes [47]). We found that 19 out of the 213 differentially expressed genes in our study were also detected among the high-scoring genes in the in vitro lipid-loaded macrophage meta-analysis gene set, corresponding to a 2.9-fold enrichment (P<0.0001) vs. for all mouse genes.

Because in vivo lipid lowering is associated with plaque regression in the Reversa and aortic transplant models (Figure 1 and [19], [20], [23], [36]), we conjectured that the gene expression changes in plaque CD68+ cells in response to lipid lowering would be associated with atheroprotective effects in humans. Because such effects might be subtle for any individual gene, and because “master regulators” of the transcriptional response are often not transcriptionally regulated themselves, we tested our hypothesis by computationally searching for connections between human atherosclerosis/CAD risk-associated single nucleotide polymorphisms (SNPs) and monocyte eQTLs for orthologs of genes that are differentially expressed in plaque CD68+ cells in response to lipid lowering (as had been previously applied to a transcriptome study of whole aortic arches [30]). We assumed that if there are no atheroprotective effects associated with the observed gene expression changes in CD68+ cells, then the fraction of atherosclerosis/CAD risk SNPs within monocyte eQTLs for the 213 lipid lowering-responsive genes should be the same as the fraction of atherosclerosis/CAD SNPs within all monocyte eQTLs. Using a list of human eQTLs from two population genetic studies of peripheral blood monocyte gene expression [48], [49], we identified SNPs associated with mRNA levels (eSNPs) of genes whose mouse orthologs are differentially expressed in CD68+ cells during regression in the two mouse models. We expanded the regression-gene-set-associated eSNPs to include all validated SNPs in linkage disequilibrium (using haplotype information from the 1,000 Genomes Project [50]), obtaining a set of regression-associated monocyte eQTLs. We then assessed the frequency of atherosclerosis/CAD risk SNPs (starting with a list of 761 risk SNPs obtained from the NCBI PheGenI database [51]) within these regression-associated monocyte eQTLs and compared these frequencies to the frequency of atherosclerosis/CAD risk SNPs within monocyte eQTLs in general. We found that the regression-associated eQTLs are substantially more enriched for atherosclerosis/CAD risk SNPs than monocyte eQTLs in general (Figure 2C and Table S6), in both regression models (P<0.05 for each model); overall, a total of 17 atherosclerosis/CAD risk-associated SNPs were found within regression-associated eQTLs (Table S7). Further investigation revealed that one of the Reversa regression-associated eQTLs encompasses rs599839, a SNP in the 1p13.3 locus that has been found to be associated with LDL cholesterol (LDL-C) levels and with risk of CAD and MI [52]–[55], as well as with hepatic mRNA expression of the intracellular lipoprotein sorting receptor Sortilin 1 (SORT1) [56], [57]. In human monocytes [48] and in liver cells [56], variant rs599839 is also an eSNP for mRNA expression of the gene Proline/serine-rich coiled coil 1 (PSRC1), which is one of most strongly upregulated genes in CD68+ cells in Mttp-inactivated vs. vehicle-treated Reversa mice (Table S2).

Identification and targeted validation of transcriptional regulators that are associated with macrophage responses to lipid lowering in vivo

We applied our novel promoter scanning analysis method, REMINISCE (Figure S1), to identify candidate transcription factors that are associated with the transcriptional responses of plaque macrophages to lipid lowering in vivo. REMINISCE combines 17 different types of macrophage chromatin and epigenomic measurements (see Table S8) that were obtained from NCBI GEO and from the Mouse ENCODE project [58] in a machine-learning algorithm (using genome regions that coincide with transcription factor binding sites that are obtained from ChIP-seq as a training dataset of regulatory regions) to map macrophage cis-regulatory regions, and then scans genomic sequences within these regions to identify statistically overrepresented transcription factor binding site sequence patterns (see Methods). In a direct comparison of REMINISCE with non-epigenome-guided motif-scanning of all noncoding sequence within ±1 kbp or ±5 kbp of the transcription start site for gene sets derived from macrophage transcriptional responses to the atherogenic lipoprotein oxidized LDL, REMINISCE detected substantially more binding site motifs above the significance threshold (Figure S2). Applying REMINISCE to the 213 lipid lowering-responsive genes from the Reversa model, we identified 15 transcription factor binding site sequence motifs that are statistically overrepresented within the 5′ regulatory regions of upregulated genes (Figure S3A), and 14 motifs that were associated with the downregulated genes (Figure S3B). Because any experimental model is subject to model-specific artifacts, we also analyzed our previously published data on the transcriptional response in CD68+ cells isolated from regressing plaques in the aortic transplant model using the REMINISCE method. We found that of the 15 motifs that were overrepresented among upregulated genes in the Reversa model, eight were also overrepresented among upregulated genes in the aortic transplant model (Figure 3A), including the motifs for PPARγ (nuclear receptor direct repeat 1, or NR DR1), the TCF/LEF family of transcription factors, and activating protein 1 (AP-1). The transcription factors represented by the 15 motifs that were associated with upregulated genes in the Reversa model are highly interconnected in the human protein-protein interactome [59] (clustering coefficient of 0.31, vs. 0.14 for the whole interactome [60]) (Figure 3B). Similarly to the case with transcription factors for upregulated genes, of the 14 motifs that were overrepresented among downregulated genes in the Reversa model, six were also detected in the aortic transplant model (Figure 3C and Table S9).

Under the hypothesis that the transcription factors that control macrophage responses to lipid lowering in vivo would have binding site sequences within the 5′ regulatory regions of a high proportion of differentially expressed genes, we ranked the motifs from Figure 3 by the fraction of the differentially expressed genes that possess a match for each motif (Table S10). Notably, for the upregulated genes, the motif whose matches were present within the largest percentage of genes, 33%, was the CTTTGA motif that is recognized by TCF/LEF and SOX transcription factor families. The finding of enrichment of TCF/LEF binding site sequences within 5′ regulatory regions of genes that are upregulated in the Reversa mouse was confirmed by additional computational analysis using high-precision motifs based on three-dimensional structural modeling of protein-DNA interactions [61] (Table S9). Because the reduced macrophage content of the plaque in the transplant regression model has been associated with macrophage emigration from the plaque, we also ranked transcription factor binding site motifs based on the proportion of cytoskeletal-binding and RhoGAP domain-containing genes (Table S4) whose 5′ regulatory regions contain matches for each motif. We found that the motif with the highest percentage, 50%, was the CTTTGA motif (corresponding to TCF/LEF and SOX) (Table S10). Next, we used the global mammalian protein interaction network to examine which signaling pathways could connect between lipoproteins or cholesterol and these transcription factors. We found that the Wnt signaling pathway [62], [63], through β-catenin and through TCF/LEF and SOX family members [64], satisfies these criteria (Figure 4).

Gene regulation likely occurs at multiple levels in plaque macrophages. Therefore, we scanned the 3′ untranslated regions (UTRs) of the two sets of differentially expressed genes in plaque CD68+ cells in the Reversa model for microRNA target sequences. We identified eighteen microRNAs whose target sequences were overrepresented within the 3′ UTRs of the upregulated and downregulated genes (six for upregulated genes, and twelve for downregulated) (Figure S4, Table S11). Under the hypothesis that some of these microRNAs might work in concert in response to lipid lowering in vivo, we analyzed molecular pathways for enrichment of the number of pathway-associated genes that possess target sequences for any of the microRNAs in Figure S4A, using a microRNA pathway analysis tool [65]. For the upregulated genes, the two pathways with the strongest enrichment were axon guidance (P<10⁻⁴) and focal adhesion (P<0.001); also enriched was the Wnt signaling pathway (P<0.05) (Table S12).

β-catenin-dependent activation of TCF/LEF is the final step in the canonical Wnt signaling pathway. To investigate the possibility of canonical Wnt pathway activation during plaque regression, we measured the relative mRNA level of β-catenin (Ctnnb1) in CD68+ cells isolated from aortic grafts from WT and Apoe^−/− recipient mice in the aortic transplant regression model. By quantitative PCR (qPCR), we found that the relative level of Ctnnb1 mRNA is approximately 2.2-fold higher in CD68+ cells in regressing plaques than in progressing plaques (Figure 5A) (P<0.05). Additionally, we measured mRNA levels of two known Wnt target genes, Lrp6 (for which loss-of-function has been associated with premature CAD [66]) and Gja1 [67], both of which were indicated by microarray profiling to be upregulated in regressing vs. progressing plaque CD68+ cells in the aortic transplant model [36]. By qPCR, both Lrp6 and Gja1 were found to be upregulated in CD68+ cells in regressing vs. progressing plaques (Figure S5). Given a previous report that β-catenin deficiency impairs macrophage migration [44] and in light of evidence of plaque macrophage emigration during regression in both the Reversa [20] and aortic transplant [22] models, we investigated whether canonical Wnt pathway activation stimulates macrophage migration in vitro. In a scratch-wound cell migration assay [68], [69], primary murine macrophages treated with the canonical Wnt pathway agonist Wnt3a scored more than two-fold higher than vehicle-treated macrophages (Figure S6).

To confirm the prediction that Wnt signaling is activated in plaque macrophages during regression in vivo, we analyzed the expression of β-catenin in plaque within aortic root sections from Reversa mice at day seven following Mttp-inactivation or vehicle treatment, by fluorescence immunohistochemistry. Consistent with the systems biological analyses, the level of β-catenin was significantly greater (approximately 2.6-fold, P<0.001) in the plaques isolated from the Mttp-inactivated vs. vehicle-treated animals (Figure 5B, 5C). Because the TCF/LEF binding site was detected as being overrepresented among genes that are upregulated during plaque regression in the transplant regression model as well, we also analyzed β-catenin expression in plaques in the aortic arch graft of WT and Apoe^−/− recipient animals. As with the Reversa model, we found that β-catenin is expressed at an approximately 2.5-fold higher level (P<0.0001) in regressing vs. progressing plaques (Figure 5D, Figure 5E, Figure S7), with increased abundance in the nuclear area (Figure S8).

Discussion

It is clear that atherosclerosis regression involves the interactions of many cellular pathways and multiple cell types and plaque components. To address this complexity, a systems biology approach employing large-scale molecular profiling is needed to identify core, model-independent mechanisms [8], [10], [70], [71]. To our knowledge, the present work, along with complementary measurements obtained from the aortic transplant regression model [36], represent the first reports of transcriptome measurements specifically of CD68+ macrophages within regressing plaques. In a previous transcriptome study [29] of changes in vessel wall tissue after lipid lowering in vivo in the Reversa mouse, a modest number of genes were detected as differentially expressed in the aortic wall (37 genes in [29], and 42 genes for arteries containing early lesions in [30]). Notably, none of these was detected as differentially expressed in our study of responses of plaque macrophages to lipid lowering, pointing to the need to separately measure the gene expression changes of constituent cell types of the plaque [71], as we have done for macrophages in the present work. While the 213 genes we identified as differentially expressed in plaque CD68+ cells in the Reversa mouse in response to lipid lowering are likely only a subset of the transcriptional changes in plaque macrophages over various time points, the enrichment of atherosclerosis/CAD risk alleles within human monocyte expression QTLs for orthologs of these genes (Figure 2C) provides intriguing support for their collective function in macrophages in the context of atherosclerosis.

Given the central role of macrophages in vascular inflammation [72] and in plaque destabilization in humans [73]–[75], our finding that Mttp-inactivation in the Reversa mouse reduces the macrophage content of the plaque (Figure 1 and [20]) points to the potential for beneficial remodeling. However, it also raises the question of how much of this reduction is driven by three processes (all of which have been reported in different models of regression): emigration from the plaque [22], deceased infiltration of monocytes [17], and increased macrophage cell death [25]. To the extent that accumulation of apoptotic macrophages within shoulder regions of advanced plaque is considered a marker of vulnerable plaques [76], the balance of these mechanisms is potentially clinically significant.

In the aortic transplant regression model, plaque macrophages in WT recipients rapidly emigrate from the graft to regional lymph nodes [22], while remaining macrophages do not exhibit abnormal levels of apoptosis [22]. Measurements of macrophage transcriptional changes during plaque regression versus progression in the aortic transplant model are consistent with an induced emigration hypothesis; macrophages in regressing (vs. progressing) plaques had lower expression levels of genes associated with adhesion and higher expression levels of genes associated with cell motility [36]. In regards to the emigration hypothesis in the Reversa model (which is supported by macrophage trafficking studies in vivo [20]), the plaque macrophage transcriptome measurements in this work provide several key insights. First, among genes that are upregulated in macrophages from Mttp-inactivated versus control animals, genes encoding cytoskeletal-binding proteins and Rho GTPases are overrepresented, consistent with increased motility. Second, among genes that are differentially expressed, no overrepresentation of apoptosis-related genes was found at the time point that we examined. A third insight, involving the Wnt signaling pathway, comes from the transcription factor analysis as described below.

In this study, we used an integrative computational approach, REMINISCE, to identify transcription factors that are likely to mediate macrophage responses during plaque regression. A novel aspect of the approach is that it incorporates macrophage epigenome measurements to guide the search for transcription factor binding sites (resulting in improved sensitivity; Figure S2) in the context of an enrichment analysis. In the analysis, we made use of transcriptome measurements from both the aortic transplant and Reversa regression models as well as macrophage-specific chromatin measurements to detect overrepresentation of known transcription factor binding site sequences upstream of genes that are differentially expressed during plaque regression. To focus on transcription factors that may have a model-independent role in macrophage responses to lipid lowering, we selected factors whose binding site sequences were overrepresented in both regression models. The significant overlap between the sets of transcription factors that were identified from the two models was unexpected due to the distinct differences in genetic backgrounds, plaque locations, and regression time-scales of the models. Of these factors, the association of PPARγ with genes that are upregulated during plaque regression is consistent with one of the principal findings from a previous transcriptome profiling study of whole aortic arches in the Reversa mouse with early lesions [30]. While our previous work investigating the effect of treatment of Reversa mice with a PPARγ agonist did not observe a significant effect of the agonist on the abundance of macrophages in the plaque [20], PPARγ is certainly worthy of future studies to clarify its precise function in plaque regression (especially given conflicting findings in different models and disease stages [20], [77]–[79]).

Among the other transcription factors whose binding site motifs were found to be associated with genes that are upregulated in CD68+ cells in regressing plaques in both regression models, the TCF/LEF and SOX transcription factor families (whose binding site sequences were associated with genes that are upregulated in response to lipid lowering) were particularly notable because of their connection to Wnt signaling, and thus, cell motility. Additionally, pathway enrichment analysis based on genes with target sequences for the microRNAs that were associated with upregulation in plaque CD68+ cells identified axon guidance and Wnt signaling as enriched pathways. Based on the transcription factor binding site analysis, we hypothesized that canonical Wnt signaling may be increased in plaque macrophages during regression. This hypothesis was confirmed by multiple observations including Ctnnb1 mRNA measurements in CD68+ cells, mRNA measurements of Wnt pathway target genes, and immunohistochemical analysis demonstrating increased β-catenin in regressing plaques in the Reversa and aortic transplant mouse models. Although the molecular mechanisms that cause increased Wnt signaling in plaque in response to lipid lowering are not known, and, of course, in spite of our associative findings, it is possible that plaque regression is mediated by Wnt-independent pathways, nonetheless, the observation of increased Wnt signaling in response to lipid lowering is consistent with previous reports that, in cultured cells, depletion of cholesterol increases canonical Wnt signaling [63].

In both the Reversa and aortic transplant models of regression, lipid lowering in vivo leads to decreased macrophage content in plaque (Figure 1B, [20], [23], [28]–[30], [36]), and trafficking studies in both of these models suggest that lipid lowering triggers macrophage emigration from the plaque [20], [22]. Our own functional studies of macrophage responses to Wnt3a stimulation (Figure S6) as well as previous studies of the role of β-catenin in macrophage migration [44], [80], are consistent with the model that increased canonical Wnt signaling – resulting in increased β-catenin levels – stimulates macrophage motility. A transitioning of macrophages from a sessile to a motile state would be expected to facilitate macrophage egress from the plaque [81].

This work is the first to our knowledge to propose that canonical Wnt signaling is involved in atherosclerosis regression in response to lipid lowering in vivo, and adds to several lines of evidence that support that alterations in Wnt signaling play a role in the progression of atherosclerosis. LRP6 genetic variants that impair Wnt/β-catenin signaling have been reported that are associated with increased risk of carotid atherosclerosis [82] and early CAD [66] in humans, suggesting that decreased canonical Wnt/β-catenin signaling may be pro-atherogenic; conversely, by this reasoning, increased Wnt signaling would be expected to have protective, pro-regressive effects. Consistent with this model, the canonical Wnt signaling molecule Wnt3a has been demonstrated to mediate anti-inflammatory effects via β-catenin-dependent suppression of the expression of proinflammatory cytokines [83]. Canonical Wnt signaling is also thought to play an important role in maintaining vascular cellular homeostasis, for example, in orchestrating lineage commitment within pericytes in the vessel wall [84]. Furthermore, Wnt signaling has been shown to potently inhibit lipid accumulation in cultured cells and to regulate plasma cholesterol levels in vivo [85]–[87], as well as (noted before) to promote macrophage migration ([44] and Figure S6).

While the present work provides evidence implicating β-catenin-interacting transcription factors in atherosclerosis regression, it is highly probable that other transcription factors are involved, given the complexity of macrophage transcriptional regulation in response to changing lipoprotein levels [70]. Of particular interest for future studies are the transcription factors that have overrepresented binding site sequences within the upregulated gene group, and that are connected to the TCF/LEF family members in the protein interaction network (Figure 3B). It should also be noted that a potential limitation of the analytical approach used in this work is the partial reliance on proximity to key factors in the global protein interaction network, in prioritizing transcriptional regulators for targeted validation; the use of such a proximity criterion could introduce bias against transcription factors that are less well-studied.

While this work used curated chromatin measurements from primary macrophages to guide the bioinformatic search for transcription factor binding sites, the numbers of cultured cells that are required for a ChIP-seq assay is rapidly being reduced; continued assay improvements may make possible genome-wide location analysis of chromatin factors within specific cellular constituents of plaque, which would likely enable more precise mapping of regulatory elements that mediate transcriptional changes during plaque regression. Equally important, comprehensive datasets of genetic variants and their associations with cardiovascular phenotypes (which may achieve significance only in a meta-analysis), such as those provided by the CARDIoGRAMplusC4D and dbGAP databases, will undoubtedly be beneficial for contextualizing the functions of these regulatory elements in atherosclerosis. Another enhancement in future studies would be to use lineage-tracing techniques in vivo to definitively characterize the extent of cellular phenotypic conversion in the Reversa model at the aortic root. Our previous work, however, has demonstrated that CD68-guided LCM of aortic root plaque yields a highly macrophage-enriched cell population [45]. Furthermore, no transcriptional changes were observed during plaque regression that would suggest a change in the proportion of macrophages in the cells captured by CD68-guided LCM in regressing vs. progressing plaques.

Beyond the enrichment of TCF/LEF transcription factor binding sites, the finding of enrichment of binding sites for the oxidative stress-responsive transcription factor nuclear respiratory factor 2 (NRF-2) within the regulatory regions of the downregulated gene group (Figure S3) is interesting in light of the observed down-regulation of mitochondrial transcription factor A (Tfam, a regulator of mitochondrial biogenesis and a known NRF target gene [88]). Downregulation of NRF-2 activity might reflect decreased oxidative stress in plaque macrophages due to the substantially reduced plasma lipoprotein levels [89].

Another intriguing finding from our work is the higher level of expression of Psrc1 in plaque macrophages in the Reversa regression model. Psrc1 mRNA has been previously been reported to be elevated in macrophage-rich regions of human atheroma [90]. The function of PSRC1 in macrophages is not known, but in other cell types it stabilizes microtubules [91], [92] and functions in cell migration [92]. Given our findings of elevated β-catenin in regressing plaques, it is interesting that a study of PSRC1 function found that PSRC1 stimulates the β-catenin pathway [93]. While strong evidence indicates that the effect of SNP rs599839 on LDL-C is mediated by hepatic expression of Sortilin 1 (SORT1) [56], the link between rs599839 and monocyte expression of PSRC1 seems worthy of further study in the context of plaque regression.

The multi-model analyses of plaque regression we have carried out is likely to ultimately shed light on the process of plaque stabilization in humans, given that atherosclerosis in genetically modified mice recapitulates many features of clinical disease. Large clinical trials have shown that in humans treated with lipid-lowering drugs, plaque area, as assessed by intravascular ultrasound (IVUS), changes very little [94], [95], while atherothrombotic event rates are dramatically decreased [96], [97], presumably due to stabilization resulting from changes in plaque macrophages not detected by IVUS. Indeed, even in mouse models, we have found that plaque size measurements do not always reflect changes in plaque macrophage content because of increased matrix proteins (e.g., [98]). Demonstrating similar phenomena in human plaques awaits improvements in non-invasive imaging. Overall, in human plaques, a drastic reduction in the abundance of macrophages within the plaque would be presumed to be associated with protective remodeling of the atheroma, and given the similarities in the disease process between mice and people, some of the factors and pathways we have identified are likely to be clinically relevant (e.g., Wnt signaling [66]).

In conclusion, this work demonstrates the value of a genomic, systems biology approach for elucidating the cellular pathways and transcriptional regulatory mechanisms that are dysregulated or altered by perturbations in vivo. Although the present study is focused on macrophage responses to lipid lowering in vivo, our approach of integrating LCM-derived transcriptome measurements with epigenomic and chromatin structural information from primary cell models could be applied to other constituents of the plaque (e.g., smooth muscle cells and endothelial cells) and then ultimately combined in a comprehensive regulatory network. More broadly, with today's panoply of publicly available epigenomic and chromatin structural measurements from a wide array of mouse and human cell types and tissues (e.g., through the ENCODE project), we anticipate that integrative approaches such as the one used in this work could be useful in the regulatory analysis of transcriptome measurements from unique, highly specific cell populations in other in vivo contexts.

Materials and Methods

Animal work

This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals (National Institutes of Health), and in accordance with animal use protocols that were approved by the NYU Langone Medical Center Institutional Animal Care and Use Committee (PHS Animal Welfare Assurance Number A3435-01) or the Institute for Systems Biology Institutional Animal Care and Use Committee (PHS Animal Welfare Assurance Number A4355-1). All tissue harvest (except mouse peritoneal macrophage harvest) and all survival surgeries were carried out under ketamine/xylazine-induced anesthesia, and for survival surgeries, animals were carefully monitored post-surgery for signs of pain or distress. For peritoneal macrophage harvest, mice were euthanized by CO₂ asphyxiation following standard operating procedures established by the Institute for Systems Biology IACUC.

The development of the Reversa mouse (JAX Stock Number 004192, Ldlr^−/− Apob^100/100 Mttp^fl/fl Mx1:Cre^+/+ animals on mixed background, in which injection of poly I:C causes silencing of hepatic Mttp expression) has been previously described [19]. In this work, for experimental studies of the Reversa model of plaque regression, four-week-old Reversa mice were weaned onto a Western diet (Catalog no. 100244, 40% kcal from milk fat; Dyets Inc., Bethlehem, PA) for 16 weeks (baseline time point) to establish aortic plaque and then switched to a chow diet and injected intraperitoneally every other day, for a total of four injections, with either 500 µg poly I:C (Sigma) or vehicle (saline). Blood samples were obtained from the retro-orbital plexus at baseline, seven days after the last poly I:C injection, and 14 days post-injection. Animals were sacrificed for arterial tissue collection at either baseline, seven days after the last poly I:C injection, or 14 days post-injection, as indicated in the Results section.

The aortic transplant model of plaque regression has been previously described [21], [23], [31]. In this work, for studies of plaque regression in the aortic transplant model, C57BL/6–Apoe^−/− mice [35] were maintained on Western diet for 16 weeks and then the aortic arches were harvested and grafted into the abdominal aortas of Apoe^−/− and WT mice. The grafts were harvested at three or five days post-surgery, as indicated in the Results section.

For all experimental data reported in this work, animal cohort sizes are as reported in the corresponding figure captions, table captions, or main article text, except for the microarray transcriptome profiling of plaque macrophages in the Reversa model, for which the cohort sizes are described in Table S1.

Cell culture

All cell culture was performed at 37°C and in the presence of 5% CO₂.

Bone marrow-derived macrophages (BMDM) were prepared as follows: L929-conditioned medium was prepared by growing L929 cells to confluency in Dulbecco's Modified Eagle Medium (DMEM)+10% (v/v) fetal bovine serum (FBS)+100 I.U./mL penicillin and 100 µg/mL streptomycin (hereafter, 1% penicillin/streptomycin). The media was then replaced with DMEM containing 2% FBS+1% penicillin/streptomycin. Conditioned media was collected every three days, filtered through a 0.22 µm filter and fresh media containing 2% FBS was added. Bone marrow cells were isolated from the femur and tibia of 8–12 week old mice. After red blood cell lysis (RLB, Sigma), cells were plated in Petri dishes in L929-conditioned medium and incubated for seven days. Cells were then serum-starved in 0.2% FBS overnight and treated with 400 ng/mL Wnt3a (R&D Systems, Minneapolis, MN) for 24 h. RNA was isolated from cells using TRIzol (Life Technologies). cDNA synthesis was performed using Verso cDNA kit (Thermo Scientific) and qPCR analysis was performed as described below (see section qPCR).

Peritoneal macrophages were prepared as follows: In each of three separate replicate experiments, resident peritoneal macrophages were isolated from 2–4 female C57BL/6J mice (ages 8–12 weeks) using intraperitoneal lavage as described previously [70]. Peritoneal exudate cells were pooled, resuspended in peritoneal macrophage medium (Roswell Park Memorial Institute 1640 medium plus L-glutamine plus 10% FBS plus 1% penicillin/streptomycin) and then plated into two tissue culture wells at approximately 10⁶ cells/well. After two hours, wells were triple-washed with warm PBS to select for adherent cells. Cells were incubated for 24 h in the presence of 25 µg/mL (by protein content) CuSO₄-oxidized human LDL (oxLDL, Intracel, Frederick, MD) or vehicle.

Scratch-wound macrophage migration assay

BMDM (prepared from N = 3 mice) were differentiated in six well plates and grown to confluency. Cells were serum-starved in 0.2% FBS overnight and a single scratch wound was made with a plastic p200 tip. Cells were washed with starvation media to remove non-adherent cells and incubated with or without Wnt3a (400 ng/mL) for 24 h. The cells were then fixed in 3.7% formaldehyde, and wound area visualized. The percentages of area covered in 24 h were quantified using NIH ImageJ software. All conditions were done in quadruplicates.

Blood and tissue analysis

Plasma cholesterol levels were determined by a colorimetric enzymatic assay (Wako Diagnostics, Richmond, VA). Mice were sacrificed at the time points indicated below (see Results) and hearts (for aortic root analysis) or aortic arch grafts were harvested and frozen in Optimum Cutting Temperature (OCT) compound and serial-sectioned at a thickness of 6 µm onto positively charged glass slides (Superfrost Plus, Fisher Scientific, Pittsburgh, PA).

For CD68 immunostaining and laser capture microdissection, slides were hematoxylin-stained, cleared with xylenes, air-dried, and foam cells were identified by light microscopy. Every seventh section was immunostained with a CD68-specific antibody (AbD Serotec, MCA 1957) as previously described [23] and used as a guide slide for laser-capture microdissection for the next six serial sections. Sections were imaged at 40× on a Leica DM4000B microscope. In addition, slides that were used for morphometric analysis were counterstained with eosin and the CD68-immunostained areas were quantified by computer-aided morphometric analysis of digitized images (Image-Pro Plus software v5, Media Cybernetics, Silver Spring, MD). For isolation of CD68+ cells from the plaques, laser-capture microdissection (LCM) was performed with the PixCell II instrument (Arcturus Bioscience) as previously described [45], [46]. All LCM procedures were performed under RNase-free conditions. Student's t-test was used for comparisons of CD68+ pixel areas and cholesterol levels for each sample group vs. baseline. A two-way ANOVA was used to test for differences in CD68+ pixel areas and cholesterol levels, between Mttp-inactivated and vehicle-treated sample groups at the two time points (seven and 14 days post-injection). In the linear model for the ANOVA, the post-treatment time point (day seven or day 14) was taken as the first factor, and the treatment (poly I:C or saline) was taken as the second factor.

β-catenin immunohistochemistry (IHC) was carried out using frozen tissue sections from aortic grafts from eight animals (N = 4 animals per genotype group). Sections were fixed in acetone for 10 min, blocked with DAKO blocking buffer (DAKO, Cat. # X0909) and incubated with β-catenin rabbit polyclonal antibody (Santa Cruz, Cat, # sc7199, 1∶200) overnight at 4°C. As a negative control, sections from Apoe^−/− recipient animals were assayed without the primary antibody. Sections were then incubated with FITC-conjugated secondary antibody for one hour, washed and mounted onto slides using VECTASHIELD mounting medium with diamidino-2-phenylindole (DAPI) (Vector Labs, Burlingame, CA). Images at 63× magnification were obtained using a Leica SP5 confocal microscope. Quantification of the mean fluorescent pixel intensity was performed on immunofluorescent images using Image-Pro Plus software v5. Six sections per animal were analyzed. An unpaired, two-tailed Student's t-test was used to test for differences in the numbers of pixels that were positive for β-catenin fluorescence, between the sample groups.

Microarrays

For LCM-derived samples, RNA was isolated from LCM caps using the PicoPure RNA Isolation Kit (Life Technologies) following the manufacturer's instructions. Total RNA was quantified using the NanoDrop spectrophotometer and assayed for quality using the Agilent Bioanalyzer 2100. Single-stranded cDNA was produced in microgram quantities from the RNA template using the Ovation Pico WTA System v2 (NuGEN), with the QIAquick PCR Purification Kit (QIAGEN) used for cDNA purification. The cDNA was analyzed on the Agilent 2100 Bioanalyzer and then fragmented and biotin-labeled using the Encore Biotin Module (NuGEN).

For peritoneal macrophage-derived samples, total RNA was isolated using TRIzol and RNA quality was analyzed with an Agilent 2100 Bioanalyzer. The three independent experiments thus yielded six RNA samples. Sample mRNA was amplified and labeled with the Affymetrix One-Cycle Eukaryotic Target Labeling Assay protocol and reagents. Labeled cDNA was hybridized to the Affymetrix GeneChip Mouse Exon 1.0 ST array using standard protocols and reagents from Affymetrix.

Fragmented, biotin-labeled cDNA was hybridized to the Affymetrix Mouse Exon Array 1.0 ST GeneChip and the GeneChip was washed and stained using the protocol and reagents in the Affymetrix GeneChip Hybridization, Wash, and Stain Kit. The stained GeneChips were imaged using the Affymetrix GeneChip Scanner 3000 and probe intensities were quantify from the scanned images using the Gene Chip Operating Software from Affymetrix. The transcriptome profiling data for this study are available online in the NCBI Gene Expression Omnibus (GEO) database, under accession numbers GSE52482 and GSE58913.

qPCR

LCM-derived samples

Real time quantitative PCR was performed using Taqman probes (Life Technologies) and using the Applied Biosystems 7300 Real time PCR system (Life Technologies). The input for each assay was 50 ng of cDNA amplified from LCM-isolated RNA using the Ovation Pico WTA system V2 (NuGEN). Taqman Universal Master Mix II (Life Technologies) was used for all reactions. C_t values were extracted from the fluorescence-versus-cycle measurements using the instrument software (Applied Biosystems). Relative mRNA levels were quantified from the C_t values using the exp(ΔΔC_t) method [99], using Gapdh as the endogenous normalizer. The specific Taqman assays that were used for the study are listed in Table S13.

Computational meta-analysis

Affymetrix microarray data from the Cho et al. study [100] and the Hägg et al. study [101] (see Table S5 for references) were obtained from NCBI GEO and processed using the Robust Multichip Average method [102] in the software package Bioconductor, to obtain background-adjusted, quantile-normalized, and probeset-summarized log₂ intensities. Human gene probesets were mapped to mouse orthologs using the Bioconductor tools for querying NCBI HomoloGene, and differential expression testing was carried out using a two-way ANOVA with a Benjamini-Hochberg [103] false discovery rate cutoff of 0.15. From all other studies, lists of differentially expressed genes were obtained from data tables from the publications referenced in Table S5, and mapped to mouse orthologs using HomoloGene. For each mouse gene in the Entrez Gene database, the number of independent analyses in which the gene was reported as differentially expressed during foam cell formation or atherosclerosis progression (among the studies in Table S5) was tabulated and used as the gene's meta-analysis score. Genes with a score greater than or equal to two (see Table S2, column L) were counted as high-scoring in the meta-analysis, for statistical enrichment test (see Results). Statistical enrichment of the set of genes with meta-analysis scores of two or greater, within the set of 213 genes that are differentially expressed in plaque CD68+ cells, was tested using Fisher's Exact Test, with the larger sample set consisting of the set of 21,594 mouse genes for which a meta-analysis score was tabulated.

Microarray data analysis

For microarray transcriptome profiling of plaque macrophages in the Reversa model, GeneChip scan images were processed into probe intensities using the Affymetrix GeneChip Operating Software system and then background-adjusted, quantile-normalized, and probeset-summarized using the Robust Multichip Average feature of Affymetrix Power Tools, with transcript-level probeset definitions from the Version 15 of the CustomCDF Project [104] (Ensembl Transcript-derived probeset definitions). Probe intensities were separately processed into exon-level normalized probeset intensities using the Robust Multichip Average Feature of Affymetrix Power Tools (version 1.10.2), and (separately) using exon-level probeset definitions from the CustomCDF project (Ensembl Exon) and from the Affymetrix Mouse Exon Array 1.0 ST NetAffx Annotations version 29. Exon-level detect-above-background calls were made using the two separate sets of (CustomCDF and Affymetrix) exon-level probesets using the Detection Above Background algorithm in Affymetrix Power Tools. For each exon-level probeset, detect-above-background P values were combined for all samples within a sample group using the geometric mean, and then across sample groups by taking the minimum P value, to obtain a single “present” P value for each exon-level probeset. For each Ensembl Transcript corresponding to one of the transcript-level probesets for which at least 75% of the constituent exon-level probesets was interrogated on the array, the transcript was considered to be above background in the sample if at least half of its constituent exons (that were interrogated with probes on the array) were detected above background with a “present” P value of less than or equal to 0.01 [105], [106]. For all probesets that passed the filter, log₂ probeset intensities across all 24 Reversa samples (see Table S1) were analyzed for differential expression between the Mttp-inactivated and control sample groups, as described below.

For comparative transcriptome profiling of plaque macrophages in the transplant regression model, probe intensity values were obtained from a two-color 3′ microarray dataset of CD68+ cells from aortic arch grafts harvested from Apoe^−/− and WT host animals three days post-surgery (N = 14 for Apoe^−/−, N = 18 for WT). This previously published dataset was obtained from the lab of one of the authors (EAF) [36], and was based on the Rosetta/Merck Mouse 23.6K 3.0 A1 microarray (custom-synthesized by Agilent; design files are available at GEO accession number GPL9733). Probe intensities in each of the two color channels (corresponding to sample-derived RNA and a universal mouse RNA reagent (Stratagene) that was used as an internal control) were combined to quantify absolute expression (geometric mean) and normalized expression (ratio of the macrophage sample-derived intensity to the internal control-derived intensity) for each of the samples. Probes were filtered for minimum absolute log₂ intensity (greater than or equal to −5) to remove probable below-background probe intensities. For all probes that passed the filter, log₂ expression ratios across all 32 samples were analyzed for differential expression between the WT and Apoe^−/− sample groups, as described below.

Testing for differential expression was performed using a two-factor model (sex and treatment for the Reversa study, sex and genotype for the transplant study) with empirical Bayes variance estimates, using the Limma package in Bioconductor [107], with a proportion factor of 0.03. Transcripts were selected as differentially expressed between the sample groups if the treatment factor P value was less than 0.01. For each unique gene among the set of differentially expressed transcripts, a representative transcript was selected by selecting (from among the differentially expressed transcripts that are associated with the gene in the Ensembl database) the transcript with the largest number of exons represented on the microarray. Enrichment tests for gene annotation terms were carried out separately for the sets of upregulated and downregulated genes in the Reversa transcriptome dataset, using the gene annotation analysis tool DAVID [108] with a P value cutoff of 0.05. Gene sets were analyzed with three different functional annotation enrichment tools to determine whether gene annotations related to apoptosis or cell death are overrepresented: DAVID, Gene Set Enrichment Analysis [109], and Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA).

Microarray analysis of poly I:C-treated macrophages was carried out using the raw microarray intensity data from a previous study [110] in which primary murine bone marrow-derived macrophages (BMDM) were treated in vitro for two hours with 6 µg/mL poly I:C and then profiled using the Affymetrix Mouse Exon Array 1.0 ST GeneChip. Six CEL files (three from each of the two sample groups, untreated and poly I:C-treated) were obtained from NCBI GEO (accession number GSE7052) and processed as described above (save for the proportion factor used for the empirical Bayes step in LIMMA being 0.1). The statistical cutoff used for differential expression tests was P<0.01. The statistical test for overlap between the poly I:C-responsive genes in BMDM and the genes that are transcriptionally responsive to lipid lowering in vivo was carried out using Fisher's Exact Test (two-sided).

Microarray transcriptome profiling probe-level measurements from mouse peritoneal macrophages were quantile-normalized without background adjustment and then summarized for probeset intensities using the Affymetrix PLIER algorithm [111] and for detection above background using the Affymetrix DABG algorithm [112] using probe-to-probeset mappings from the CustomCDF project [104] version 15, based on the mouse transcripts and mouse exons lists from the Ensembl database release 65. Transcripts were marked as “present” only if at least 2/3 of the constituent exons had a DABG P value of less than 0.01 [105]. Differential expression testing for each Ensembl transcript was carried out using a t-test with empirical Bayes variance estimates, using the Limma package [107] in Bioconductor with 0.1 as the prior expected proportion of differentially expressed genes. P values for differential expression between the sample groups (vehicle and oxLDL) were converted to q-values using the Bioconductor package qvalue [113]. Transcripts were filtered with a false discovery rate cutoff of 0.05 and a minimum absolute fold-change cutoff of 2.0, and filtered to include only transcripts that could be mapped to the RefSeq NM identifier. Transcripts were mapped start site coordinates on the mm9 genome assembly, and noncoding sequences within ±5 kbp of the transcription start sites were retrieved, using the BioMart tool in Ensembl release 65.

Transcription factor binding site motif analysis

In this section, we describe the computational method REMINISCE (Regulatory Element Machine-learning to Infer Networks Instructing Specific Cellular Expression; see Figure S1) that was used to identify transcription factor binding site motifs whose matches are overrepresented within 5′ regulatory regions of differentially expressed genes. Using the BioMart interface to the Ensembl database (release 65), transcription start site locations (on the NCBI M37 mouse genome assembly) were obtained for the sets of differentially expressed transcripts that were detected by microarray in CD68+ cells in plaques in the Reversa or transplant regression models. Noncoding DNA sequences (defined by excluding locations marked as protein coding within at least one exon within the Ensembl database) within ±5 kbp of a transcription start site were analyzed using a Random Forest classifier [114] to identify probable cis-regulatory regions, with a score assigned at 100 bp intervals based on the voting fraction of trees in the classifier. The features used for the classification (Table S8) include genome-wide measurements (binarized peak locations) of histone acetylation in bone marrow-derived macrophages obtained by ChIP-seq [39], “valley scores” in histone acetylation ChIP-seq signal [39], DNase I hypersensitive sites in bone marrow-derived macrophages [39] detected by short-read sequencing, vertebrate PhastCons interspecies sequence conservation scores from the University of California Santa Cruz Genome Bioinformatics database [115], GC content, computational predictions of CpG island locations [116], and basepair distances to the nearest binarized feature locations of each type. The Random Forest classifier for predicting cis-regulatory regions was trained using compendium of transcription factor binding site locations obtained by combining ChIP-seq datasets targeting 19 transcriptional regulators in murine macrophages, that were obtained from the literature (Table S8). Noncoding sequence within regions predicted to have a cis-regulatory function by the Random Forest model (voting fraction greater than 0.5) was scanned for transcription factor binding site motif matches using a combined set of vertebrate binding site position weight matrices from the TRANSFAC Professional 2010 [117] and JASPAR 2010 [118] databases. The scanning was performed using the software tool Clover [119] (Feb. 19, 2010 release) on noncoding sequences within cis-regulatory regions (after masking low-complexity repeats using NCBI DustMasker) proximal to transcription start sites of differentially expressed transcripts as well as proximal to transcription start sites of a set of 1,000 randomly selected genes that were detected above background in plaque macrophages by microarray hybridization (see Microarray Analysis). For the Clover scanning, a pseudocount value of 0.3 [120] was used. Matches for binding site sequence patterns within the promoter sequences for upregulated or downregulated genes were assessed for statistical overrepresentation using the frequencies compiled on the promoter regions for the 1,000 randomly selected macrophage-expressed genes, using the binomial test (P value threshold of 0.01). Subsequent to the motif analysis using TRANSFAC 2010, 5′ cis-regulatory regions for the upregulated Reversa gene set were re-analyzed using a substantially larger set (numbering 56) of TCF/LEF transcription factor family motifs from an updated motif database (TRANSFAC Professional 2013.4), using the same methods.

Transcription factor network analysis

Official gene symbols for genes encoding protein components of transcription factors that were implicated in the transcription factor binding site enrichment analysis (see Figure 3A) were uploaded to the GeneMANIA network analysis web application [59]. The interaction network for the transcription factors was assembled using protein-protein interactions (direct physical interactions only), and visualized using Cytoscape [121].

MicroRNA target sequence enrichment analysis

MicroRNA enrichment analysis was carried out using the CORNA software algorithm [122] using conserved microRNA target predictions from miRBase (v5) and using Ensembl BioMart [123] (release 70) for gene identifier translation. Statistical enrichment was assessed using Fisher's Exact Test (P value threshold of 0.01). Pathway enrichment analysis of genes with target sequences for the microRNAs shown in Figure S4 was carried out using the analysis program DIANA-miRPath [65], with target sequences predictions from the DIANA-microT-CDS v.5.0 algorithm [124].

Human SNP analysis of plaque regression gene sets

A set of 239,165 eSNP-to-gene associations (that collectively connect 157,668 eSNPs to monocyte expression levels of 13,989 genes) was compiled by merging (eSNP,gene) pair lists from two human population genetic studies of monocyte gene expression regulation [48], [49] comprised of samples from 1,490 and 283 individuals, respectively. From these 239,165 pairs, eSNPs that were associated with human genes whose mouse orthologs (obtained via successive mapping using the HomoloGene version 68 and the INPARANOID [125] version 8.0 databases) were detected as differentially expressed in plaque CD68+ cells in the Reversa model of plaque regression (for Reversa, the 213 genes in Table S2) or in the aortic transplant model of plaque regression (from Table S2 of [36], the 549 genes with absolute geometric-mean expression ratio >1.25 between the WT and Apoe^−/− recipient sample groups) were then obtained, resulting in lists of 2,380 (Reversa) and 5,024 (transplant) plaque regression-associated eSNPs, respectively. The three sets of eSNPs – Reversa eSNPs, transplant eSNPs, and the full 157,668 “monocyte eSNPs” – were expanded to include proxy SNPs that are both in linkage disequilibrium (LD) with an eSNP (R>0.9) and within 250 kbp of the eSNP, yielding LD-expanded sets of 6,214, 13,248, and 363,283 SNPs, respectively. The expansion to proxy SNPs was carried out using LD maps from the 1,000 Genomes Project (CEU population group) and using the SNP analysis software SNAP [126]. A set of 761 SNPs that are associated with risk of atherosclerosis or CAD was obtained by searching the NCBI Phenotype-Genotype Integrator database [51] for the MeSH terms “atherosclerosis”, “coronary artery disease”, and “carotid stenosis”. Overlaps were computed for the set of 761 atherosclerosis/CAD SNPs with four different sets of SNPs: the 6,214 LD-expanded eSNPs from the Reversa model (overlap = 6), the 13,248 LD-expanded eSNPs from the aortic transplant model (overlap = 11), the 363,283 LD-expanded monocyte eSNPs (overlap = 134), and the set of all 44,278,189 validated SNPs in dbSNP build 138 (overlap = 761). From these overlaps, two 2×2 contingency tables were constructed to test for association between atherosclerosis/CAD SNPs and regression gene-associated eSNPs in the Reversa and aortic transplant models, respectively, using the LD-expanded set of monocyte eSNPs as a background set. Additionally, two 2×2 contingency tables were constructed to test for association between atherosclerosis/CAD SNPs and regression gene-associated eSNPs in the Reversa and aortic transplant models, respectively, using the set of all validated SNPs as a background set. P values and odds ratios (Figure 2C and Table S6) were obtained for these four contingency tables using Fisher's Exact Test. Confidence intervals (shown in Figure 2C) on the expected overlap of atherosclerosis/CAD SNPs and regression-associated eSNPs were obtained using the quantile function of the hypergeometric distribution.

Analysis software

Computational analyses were carried using scripts in the R statistical computing environment, also using the Bioconductor set of packages for R), and in the Perl programming language. Bar graphs were prepared using Prism (GraphPad Software, La Jolla, CA). All custom software scripts used in the computational analysis are available from the authors upon request. Statistical tests were carried out using R and GraphPad Prism.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21