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

We report the first systems biology investigation of regulators controlling arterial plaque macrophage transcriptional changes in response to lipid lowering in vivo in two distinct mouse models of atherosclerosis regression. Transcriptome measurements from plaque macrophages from the Reversa mouse were integrated with measurements from an aortic transplant-based mouse model of plaque regression. Functional relevance of the genes detected as differentially expressed in plaque macrophages in response to lipid lowering in vivo was assessed through analysis of gene functional annotations, overlap with in vitro foam cell studies, and overlap of associated eQTLs with human atherosclerosis/CAD risk SNPs. To identify transcription factors that control plaque macrophage responses to lipid lowering in vivo, we used an integrative strategy – leveraging macrophage epigenomic measurements – to detect enrichment of transcription factor binding sites upstream of genes that are differentially expressed in plaque macrophages during regression. The integrated analysis uncovered eight transcription factor binding site elements that were statistically overrepresented within the 5′ regulatory regions of genes that were upregulated in plaque macrophages in the Reversa model under maximal regression conditions and within the 5′ regulatory regions of genes that were upregulated in the aortic transplant model during regression. Of these, the TCF/LEF binding site was present in promoters of upregulated genes related to cell motility, suggesting that the canonical Wnt signaling pathway may be activated in plaque macrophages during regression. We validated this network-based prediction by demonstrating that β-catenin expression is higher in regressing (vs. control group) plaques in both regression models, and we further demonstrated that stimulation of canonical Wnt signaling increases macrophage migration in vitro. These results suggest involvement of canonical Wnt signaling in macrophage emigration from the plaque during lipid lowering-induced regression, and they illustrate the discovery potential of an epigenome-guided, systems approach to understanding atherosclerosis regression.


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 lowdensity lipoprotein (LDL) cholesterol [4,5]. There is significant interest in finding new therapies to promote the regression of atherosclerotic plaques [6][7][8], with particular emphasis on the plaque macrophage as a potential target because of its roles in lipoprotein uptake and vascular inflammation [9][10][11][12], plaque destabilization [13,14], and plaque remodeling in response to dietary changes [15].
In mice, plaque regression has been achieved through expression [16][17][18] or inactivation [19,20] of lipid metabolismmodifying genes; aortic tissue grafting [21][22][23]; and treatment with small-molecule therapeutics [24][25][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 2/2 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 2/2 mouse [35] is grafted in place of a segment of the abdominal aorta [31] of either a wild-type (WT) or an Apoe 2/2 mouse. Grafts within WT recipients exhibit plaque regression, while plaques within Apoe 2/ 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 59 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 59 regulatory regions of genes that are upregulated in plaque macrophages during regression.
TCF/LEF transcription factors are activated by nuclear bcatenin (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 b-catenin-dependent mechanism [44]. In this study we show, for the first time, that activation of canonical Wnt signaling, defined by up-regulation of b-catenin, occurs within macrophage-rich plaque in two mechanistically distinct lipidlowering 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 (903671 mg/dL total cholesterol [mean 6 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

Author 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 alreadyestablished 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.
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 Mttpinactivated vs. vehicle-treated animals at seven and 14 days (P, 0.05); at day 14, CD68+ area was 26% smaller in the Mttpinactivated 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 2/2 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 2/2 mice. We found that only one of the 213 lipidresponsive genes, Trnt1, was detected as differentially expressed in the Reversa saline vs. Ldlr 2/2 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). Mttp inactivation contributes to a reduction in macrophage content of the plaque. Bars are mean 6 SE. CD68+ area differences between vehicle and Mttp-inactivated are: at day seven, 13%; at day 14, 26%. The CD68+ area reduction in Mttp-inactivated vs. vehicle groups, at both days seven and 14, is significant with P,0.05 (1, two-way ANOVA). For all bars, N = 7 except day seven, where N = 9. The difference in CD68 areas between day 14 (or day seven) saline-treated groups and baseline were not statistically significant. (C) Mttp inactivation reduces CD68 cellular content within aortic root plaque. Micrographs are representative aortic root sections from day 14 animals (from vehicle and Mttp-inactivated groups, as indicated) that have been counterstained with hematoxylin (blue) and immunostained for CD68 (red). doi:10.1371/journal.pgen.1004828.g001 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 regressiongene-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 regressionassociated 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][53][54][55], as well as with hepatic mRNA expression of the intracellular lipoprotein sorting receptor Sortilin Reversa mice were fed a Western diet for 16 weeks. One group of animals was sacrificed at 16 weeks as the baseline group. The rest of the mice were switched to chow diet and assigned to two groups, vehicle and Mttp-inactivated; animals within these two groups received injections of saline or poly I:C, respectively (see Methods). Animals were sacrificed on day seven or day 14 after the injections, and plasma total cholesterol levels were measured (see Methods Reversa mice were fed a Western diet for 16 weeks. One group of animals was sacrificed at 16 weeks as the baseline group. The rest of the mice were switched to chow diet and assigned to two groups, vehicle and Mttp-inactivated; animals within these two groups received injections of saline or poly I:C, respectively (see Methods). Animals were sacrificed on day seven or day 14 after the injections, and aortic roots were sectioned, immunostained for CD68, and the CD68+ areas were measured (see Methods 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 Mttpinactivated 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 nonepigenome-guided motif-scanning of all noncoding sequence within 61 kbp or 65 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 59 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 PPARc (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 59 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 59 regulatory regions of genes that are upregulated in the Reversa mouse was confirmed by additional computational  Table S2). Edges connect genes to the statistically enriched annotation categories with which they are associated (see Methods; enrichment P values are given in Table S4). (A) A statistically significant cluster of upregulated genes is associated (by gene annotations) with the cytoskeleton. (B) A statistically significant cluster of downregulated genes is associated with the mitochondrion. (C) Human monocyte expression QTLs for genes that are differentially expressed in plaque CD68+ cells in regression are enriched for SNPs associated with risk of atherosclerosis or CAD, versus the number of atherosclerosis/CAD SNPs for monocyte eQTLs overall. The bars marked ''expected'' show the range of numbers of SNPs associated with risk of atherosclerosis or CAD that would be expected to occur by chance within the eQTLs of the set of 213 lipid lowering-responsive genes, based on the overall percentage of SNPs that are associated with risk of atherosclerosis/CAD (''random'') and based on the percentage of SNPs within all monocyte eQTLs that are associated with atherosclerosis/CAD (''monocytes''). Error bars indicate .90% confidence intervals for the number of atherosclerosis/ CAD SNPs that would be expected by chance (see Methods). (*) P,0.05, Fisher's Exact Test (contingency table data are in Table S6). doi:10.1371/journal.pgen.1004828.g002 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 59 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 b-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 39 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 39 UTRs of the upregulated and downregulated genes (six for upregulated genes, and twelve for downregulated) ( Figure S4, Table S11). Under the Transcription factor binding site motifs whose matches are statistically overrepresented within 59 regulatory regions of genes that are upregulated in plaque macrophages in response to lipid lowering in vivo. Bars indicate the ratio of the number matches of the indicated motifs per kbp, within an upregulated gene set, to the number of motif matches per kbp for randomly selected sets of genes expressed in CD68+ cells. Bar color indicates the regression model: white, Reversa model; gray, aortic transplant model. Only motifs for which the ratio exceeds 1.5 (with a statistical significance of P, 0.01) in both models are shown (complete list in Table S9). Motifs are labeled by the transcription factor family name with which they are associated (motif IDs are provided in Table S9). Each motif is represented by a four-color motif logo. (B) The transcription factors that are associated with upregulated genes are highly interconnected at the level of protein-protein interaction network. Blue hexagons represent transcription factors or transcription factor families, and edges denote protein-protein interactions between them. (C) Transcription factor binding site motifs whose matches are statistically overrepresented within 59 regions of genes that are downregulated in plaque macrophages in response to lipid lowering in vivo. ''NR (HS)'', nuclear receptor half-site. doi:10.1371/journal.pgen.1004828.g003 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 24 ) and focal adhesion (P, 0.001); also enriched was the Wnt signaling pathway (P,0.05) ( Table S12).
b-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 b-catenin (Ctnnb1) in CD68+ cells isolated from aortic grafts from WT and Apoe 2/2 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 b-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 b-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 b-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 b-catenin expression in plaques in the aortic arch graft of WT and Apoe 2/2 recipient animals. As with the Reversa model, we found that b-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). . Predicted molecular interaction network for upregulated genes that possess CTTTGA elements connects cholesterol to motility. Red ellipses represent genes whose expression levels are elevated in plaque macrophages from Mttp-inactivated vs. vehicle-treated animals (darker color indicates a greater differential expression ratio; see Table S2). Yellow hexagons and blue edges represent the upstream signaling pathway based on the literature (the blunt-ended edge indicates inhibition of signaling), with green rectangles representing lipids (cholesterol: membrane cholesterol; lipoprotein: apoB-containing lipoprotein). Blue diamonds represent transcription factor binding site motifs. Gray edges connect the transcription factors to the downstream genes that possess corresponding binding site sequence matches. doi:10.1371/journal.pgen.1004828.g004

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][74][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 PPARc 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 PPARc agonist did not observe a significant effect of the agonist on the abundance of macrophages in the plaque [20], PPARc 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][78][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 b-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][29][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 b-catenin in macrophage migration [44,80], are consistent with the model that increased canonical Wnt signaling -resulting in increased bcatenin 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/b-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/bcatenin 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 b-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][86][87], as well as (noted before) to promote macrophage migration ( [44] and Figure S6).
While the present work provides evidence implicating bcatenin-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 CARDIoGRAM-plusC4D 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 b-catenin in regressing plaques, it is interesting that a study of PSRC1 function found that PSRC1 stimulates the b-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 lipidlowering drugs, plaque area, as assessed by intravascular ultrasound (IVUS), changes very little [94,95], while athero-thrombotic 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 noninvasive 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.

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/xylazineinduced 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 2 asphyxiation following standard operating procedures established by the Institute for Systems Biology IACUC.
The development of the Reversa mouse (JAX Stock Number 004192, Ldlr 2/2 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 mg 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 2/2 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 2/2 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 37uC and in the presence of 5% CO 2 .
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 mg/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 mm 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 6 cells/well. After two hours, wells were triplewashed with warm PBS to select for adherent cells. Cells were incubated for 24 h in the presence of 25 mg/mL (by protein content) CuSO 4 -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 mm 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 406 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.
b-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 b-catenin rabbit polyclonal antibody (Santa Cruz, Cat, # sc7199, 1:200) overnight at 4uC. As a negative control, sections from Apoe 2/2 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 636 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 b-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(DDC 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 2 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 metaanalysis 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-normal-ized, 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. Exonlevel 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 2 probeset intensities across all 24 Reversa samples (see Table S1) were analyzed for differential expression between the Mttpinactivated 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 39 microarray dataset of CD68+ cells from aortic arch grafts harvested from Apoe 2/2 and WT host animals three days post-surgery (N = 14 for Apoe 2/2 , 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 samplederived intensity to the internal control-derived intensity) for each of the samples. Probes were filtered for minimum absolute log 2 intensity (greater than or equal to 25) to remove probable belowbackground probe intensities. For all probes that passed the filter, log 2 expression ratios across all 32 samples were analyzed for differential expression between the WT and Apoe 2/2 sample groups, as described below.
Testing for differential expression was performed using a twofactor 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 marrowderived macrophages (BMDM) were treated in vitro for two hours with 6 mg/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 foldchange 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 65 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 59 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 65 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, 59 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 geometricmean expression ratio .1.25 between the WT and Apoe 2/2 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 LDexpanded 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 262 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 262 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 regressionassociated 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. Figure S1 Diagram of REMINISCE method for detecting enrichments of transcription factor binding sites for specific transcription factors, within predicted regulatory regions. Cell type-specific epigenomic and chromatin information (A) are combined with cell type-generic genomic information (B) within an ensemble decision tree classifier, Random Forest (C). The classifier is trained to use the input features (A,B) to predict the locations of cis-regulatory regions, identified by combining transcription factor location datasets (D) from the specific cell type (macrophages). The classifier predicts cis-regulatory regions (E) upstream of genes that are differentially expressed in plaque macrophages in response to lipid lowering in vivo, and these regions are scanned to identify matches for transcription factor binding site motifs from a precompiled library of vertebrate motifs (F). For each transcription factor, frequencies per bp of predicted cis-regulatory region sequence for its binding sites are tabulated and tested for enrichment (G) above a background frequency from cis-regulatory regions for a list of genes that are expressed above background level in murine macrophages. See Methods section for details. (TIFF) Figure S2 Sensitivity for in silico detection of enrichments of transcription factor binding sites is compared between the epigenome-guided method, REMINISCE, and approaches in which all noncoding sequence within 61 kbp or 65 kbp of transcription start sites are analyzed. The three analysis methods were applied to a transcriptome profiling study of macrophage foam cells in which mouse resident peritoneal macrophages (N = 3 independent replicate experiments; see Methods) were cultured with oxLDL or vehicle for 24 h. High-confidence sets of differentially expressed genes were identified (false discovery rate cutoff of 0.05 and minimum absolute fold-change of 2.0) as upregulated (94 genes) or downregulated (342 genes), and a background set of 2,000 genes that were randomly selected from the set of genes that were detected above background in at least one sample group (vehicle or oxLDL). The 59 regulatory regions for the genes were scanned for matches to transcription factor binding site motifs using the epigenome-guided method, REMI-NISCE, as described in the Methods section (''REMINISCE method'', white bars) and by analyzing all noncoding sequence within 61 kbp of transcription start sites (''1 kbp classic'', gray bars) or within 65 kbp of transcription start sites (''5 kbp classic'', charcoal bars). For all three methods, the number of matches for a binding site motif in the biological gene sets (genes up-or downregulated by oxLDL) per kbp of sequence analyzed, was compared to the number that would be expected by chance for genes that are expressed above background in resident peritoneal macrophages. Bar group labels represent transcription factor binding site motifs. Bar lengths represent the ratio of the number of motif matches (per kbp of sequence analyzed) for the indicated oxLDL-response gene set (left, upregulated in oxLDL; right, downregulated in oxLDL) to the number of motif matches (per kbp of sequence analyzed) for the background set of genes. For each bar group (motif), a missing bar of a particular shade indicates that the motif's binding sites were not detected as enriched (P,0.01) using that method. The results indicate that the REMINISCE method detected substantially more motifs, and at higher enrichment ratios, than the two methods based on analysis of all noncoding sequence. (TIFF) Figure S3 Motif scanning enrichment analysis for promoters of genes that are differentially expressed in plaque macrophages during plaque regression, in the Reversa and/or in the transplant regression models. Bars represent the ratio of the number of transcription factor (TF) binding site motif matches (per kbp of cisregulatory sequence analyzed) for the indicated differentially expressed gene set, to the number of binding site motif matches per kbp for randomly generated sets of genes that are expressed above background in plaque CD68+ cells. Bar shading indicates the regression model from which the gene set was derived (white, Reversa model, gray, transplant regression model). A missing gray bar indicates that binding sites for the indicated TF binding site motif were not detected as significantly higher than the number expected by chance, for the transplant regression model. TF motifs shown on this figure were selected based on a significance threshold (see Methods) for enrichment for gene sets derived from the Reversa mouse, and thus, every TF shown here has a white bar. (A) TF binding site motifs for genes that are upregulated in CD68+ cells from regressing animals vs. control animals. (B) TF binding site motifs for genes that are downregulated in CD68+ cells from regressing animals vs. control animals. (TIFF) Figure S4 microRNA target sequence enrichment analysis identifies many microRNAs whose target sequences are overrep-   Figure 2C