Single-cell analysis of fate-mapped macrophages reveals heterogeneity, including stem-like properties, during atherosclerosis progression and regression

Atherosclerosis underlies coronary artery disease, a leading cause of death in the world. Atherosclerosis already begins its progression in childhood (1), making plaque reversal an important clinical goal during adulthood. There is increasing recognition that plaque progression represents a chronic condition that results from a failure to resolve inflammation (2). The central inflammatory cell in the plaque is the macrophage (3). Though macrophages in tissues can originate from resident macrophages that are seeded in the tissues during embryonic development (e.g., refs. 4, 5), the bulk of plaque macrophages is most likely derived from blood monocytes recruited during disease progression (3). Macrophage proliferation has also been identified as a feature of plaques (6), although the origin of these proliferative cells is unclear.

The two subsets of blood monocytes in mice are often defined by the expression of chemokine receptors: Ccr2⁺Cx3cr1⁺(Ly6C^hi) for classical monocytes and Ccr2^–Cx3cr1⁺⁺(Ly6C^lo) for patrolling nonclassical monocytes, and they have distinct migratory and inflammatory properties (7). Ly6C^hi classical monocytes utilize CCR2 and CX3CR1 to enter atherosclerotic lesions in Apoe^–/– mice (8, 9) and are thought to become classically activated, or M1, macrophages under most inflammatory conditions (9–11). However, alternatively activated M2 macrophages can also be derived from Ly6C^hi CCR2-dependent monocytes during helminth infection (12), in allergic inflammation (13), and, as noted below, in regressing atherosclerotic plaques (14). Hence, as newly emigrating Ly6C^hi monocytes are exposed to different environmental stimuli in the tissues, they will respond to the signals that result in different activation states.

Based on histochemical markers, the majority of macrophages in both mouse and human progressing plaques resemble the activated classical M1 phenotypic state. We have established a number of different mouse models to find that plaque regression is characterized not only by reduced classically activated M1 macrophages, but also by the enrichment of cells expressing markers of alternatively activated (M2 or M[IL-4]) macrophages (3, 15, 16). Alternatively activated M2 macrophages have been shown to participate in resolving inflammation and repairing tissue damage, consistent with features of plaque regression.

This type of macrophage can be derived from tissue-resident macrophages or macrophages derived from classical (Ly6C^hi) or nonclassical patrolling (Ly6C^lo) monocytes. We recently demonstrated that plaque regression is driven by the CCR2-dependent recruitment of macrophages derived from inflammatory Ly6C^hi monocytes that adopt features of the M2 state in a STAT6-dependent manner (14). This suggests that in both progressing and regressing plaques, classically and alternatively activated macrophages are both derived from inflammatory Ly6C^hi monocytes. The full scope of different macrophage activation states after transition from monocytes, however, is only just being revealed by single-cell analysis during plaque progression (17, 18) and, notably, is still unknown for plaque regression. Also, the traditional definition of M1 and M2 macrophage activation states often represents polar extremes that are based on in vitro activation conditions with high concentrations of stimuli and on a small number of markers. Thus, the typical conditions of studies in vitro probably do not reflect the more complex in vivo physiological state in a number of key ways, further contributing to the incomplete understanding of monocyte-to-macrophage maturation process in inflammatory conditions, with the process likely to be tissue specific (19).

To improve the understanding of the origins and fates of macrophages in atherosclerotic plaques undergoing dynamic changes, we have combined single-cell RNA-Seq with genetic fate mapping of myeloid cells derived from CX3CR1⁺ precursors for application in a mouse model in which plaques form and then are induced to regress. This not only greatly increases the resolution of detail over what is afforded by the limited number of markers typically used to study macrophage phenotypes, but also allows extensive characterizations in the in vivo setting. As we will describe, in atherosclerotic plaques there is a spectrum of macrophage activation states with greater complexity than the traditional M1/M2 definitions, with progressing plaques containing more discernible macrophage activation states than during regression. We also found a population of proliferating cells, remarkably, with monocyte markers and stem cell–like signatures, that may represent a new self-renewing source of macrophages in both progressing and regressing plaques.

Fate mapping the conversions of plaque macrophages derived from CX3CR1⁺ precursors during atherosclerosis progression and regression. All blood monocytes that migrate into atherosclerotic plaques express CX3CR1 (20, 21); hence, we first examined the fate of these monocytes during atherosclerosis progression by generating BM chimeras of Ldlr^–/– mice reconstituted with BM from Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice, which were then fed an atherogenic Western diet (WD). We took this approach because we previously utilized this tamoxifen-inducible (TAM-inducible) Cre recombinase (CreER) system under the control of the Cx3cr1 promoter to fate map monocyte-derived macrophages without adoptive transfer in a schistosomiasis model (5). TAM treatment irreversibly and genetically labels CX3CR1⁺ cells and causes them to express tdTomato. Thus, the Ldlr^–/–:Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} BM chimeras were treated with 2 doses of TAM at 14 and 15 weeks of WD, and the aortic root plaques were examined after 18 total weeks of WD feeding, which resulted in advanced plaques (Supplemental Figure 1A; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.124574DS1).

As shown in Figure 1A, newly recruited CX3CR1-EYFP⁺ but TdTomato^– cells were mostly observed in an abluminal, subendothelial location. In contrast, TdTomato⁺EYFP⁺ cells were observed further inward, toward the lipid core. Both populations were found in the adventitia, with somewhat more being TdTomato⁺EYFP⁺ (Figure 1A). Additionally, we observed considerable heterogeneity in shape (including elongated cells, small and foamy macrophages) among the TdTomato⁺ cells derived from CX3CR1⁺ precursors (Supplemental Figure 1B).

Figure 1

Fate mapping the conversion of CX3CR1⁺ cells into plaque macrophages. (A) Representative confocal images of aortic roots stained with CD68 (green), EYFP (yellow), and TdTomato (red) in BM chimeras of Ldlr^–/– mice reconstituted with BM from Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice (n = 3) gavaged with tamoxifen (TAM) at 14 and 15 weeks after feeding on a Western diet (WD) to label cells derived from CX3CR1⁺ monocytes. Scale bars: 50 μm. (B–F) Analysis of aortic arches by flow cytometry of Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice injected with AAV-PCSK9 and fed WD for 18 weeks before TAM treatment. Progression group mice (n = 4) were then kept on WD, while regression group mice were switched to chow and treated with ApoB-ASO for 2 weeks (n = 4). (B) Representative bright-field images and quantification of lesion areas of aortic roots. Scale bars: 50 um. (C) Density plot of t-distributed stochastic neighbor embedding (t-SNE) analysis of CD11b⁺TdTomato⁺ cells from aortic arches of mice in progression and regression groups subjected to aortic digestion. The aortic arches were analyzed by flow cytometry for expression of PD-L2, CD301, EYFP, F4/80, and MHCII markers (n = 8). (D) Heatmaps of geometric mean fluorescence, (E) quadrant plots, and (F) geometric MFI of PD-L2 and CD301 expression on CD11b⁺TdTomato⁺ cells from progression and regression groups in atherosclerotic Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice with AAV-mPCSK9 induction. Progression, n = 4–9; regression, n = 4–8. Statistical significance was calculated using Student’s t test, and data are presented as mean ± SEM (B and F). *P < 0.05, ***P < 0.001, ****P < 0.0001.

To analyze the phenotype cells derived from CX3CR1⁺ precursors during atherosclerosis progression and regression without generating BM chimeras in Ldlr^–/– mice, we utilized a variant of our recently reported model of atherosclerosis regression. Plaque progression is initiated by injecting into Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice an adeno-associated viral vector expressing a gain-of-function mutant of protein convertase subtilisin/kexin type 9 (AAVmPCSK9), which results in LDL receptor deficiency and hypercholesterolemia. After plaques form, to initiate regression, plasma lipid levels are lowered by using an antisense oligonucleotide (ASO) to apolipoprotein B (ApoB) (22, 23), which reduces LDL production. AAVmPCSK9-treated Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice were fed WD for 18 weeks, then treated with TAM. To compare atherosclerosis progression with regression, mice in the progression groups were kept on WD for an additional 2 weeks, whereas mice in the regression groups were switched to chow and treated with ApoB-ASO for 2 weeks to lower plasma LDL levels (Supplemental Figure 2A).

We first confirmed that plaque areas were significantly decreased in the regression compared with the progression group (Figure 1B). We then performed FACS analysis to characterize TdTomato⁺ cells from aortas (Supplemental Figure 2B). While the overall number of TdTomato⁺ cells from aortas of mice undergoing progression and regression was not significantly different (Supplemental Figure 2C), there were phenotypic differences between the fate-mapped cells from the progression and regression groups, as visualized by t-distributed stochastic neighbor embedding (t-SNE) (Figure 1C). Notably, TdTomato⁺ cells from the regression group expressed higher levels of PD-L2 and CD301 (Figure 1, D–F), which are both markers of alternatively activated M2 macrophages (12, 24).

Single-cell RNA-Seq analysis of aortic cells derived from CX3CR1⁺ precursors in atherosclerosis progression and regression. To further define the molecular features associated with progression and regression in cells derived from the CX3CR1⁺ monocyte lineage, we FACS purified TdTomato⁺ CD11b⁺ cells from the aortic arches in the progression and regression groups. This allowed us to focus on cells belonging to the myeloid lineage by excluding T cells, NK cells, B cells, eosinophils, and neutrophils in the dump gate (Supplemental Figure 2B). We obtained between 1000 and 6400 TdTomato⁺CD11b⁺ cells from each mouse and combined 4 samples for a total of 10,000–12,000 cells in each group. Single-cell RNA-Seq was performed on the 10x Genomics platform, which resulted in transcriptome profiles (after quality control filtering) of 3157 cells in the progression group and 2198 cells for the regression group. After data normalization, Louvain clustering of the aggregated data combined from the progression and regression groups identified 11 cell clusters, which were visualized using multicore t-SNE (Figure 2, A and B). Importantly, cells from all of the clusters expressed canonical myeloid cell, monocyte, and/or macrophage markers, including Csf1r, Cd14, Adgre1 (or F4/80), and Cd68 (Supplemental Figure 3A), indicating that we successfully excluded other hematopoietic cells that can express CX3CR1. The majority of the Louvain clusters are contiguous, which suggests a spectrum of related macrophage activation states, apart from the more distinct clusters 0, 3, 5, and 7 (Figure 2A) (discussed below).

Figure 2

Heterogeneity of plaque macrophages derived from CX3CR1⁺ monocyte precursors in atherosclerosis progression and regression. Mice were treated as described in Figure 1. (A and B) We used Louvain clustering and multicore t-SNE to visualize 5355 CD11b⁺TdTomato⁺ single cells isolated from aortic arches of Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} mice from the progression and regression groups (points, n = 4 mice per group; colored by [A] sub-cluster and [B] experimental group). (C) Cluster composition by percentages of experimental group (red, progression; blue, regression) in total sorted CD11b⁺TdTomato⁺ cells. (D) Cell type signatures are shown in heatmap in the relative expression level. Row-wise Z score of ln(X + 1), where X denotes transcript count per cell after normalization (mean = 0, SD = 1); color scale of genes (rows) across cell clusters (columns) is shown.

Cluster 1 included the largest number of cells from both the progression and regression groups, at a relatively similar frequency, indicating a common feature of atherosclerosis (Figure 2C). Folate receptor β (Folr2) is the most differentially expressed gene in cluster 1 (herein referred to as “Folr2^hi macrophages”; Figure 2D, Figure 3, and Supplemental Figure 3B). Expression of Folr2 was previously found to be increased in human atherosclerotic plaques, but not in normal arteries (25), and also to be expressed in a population of macrophages termed in a recent single-cell experiment as “resident-like macrophages” (18). Of special interest, the transcriptional profile of the cells in cluster 4 is dominated (Figure 2D, Figure 3, and Supplemental Figure 3B) by the expression of chemokines and cytokines (Ccl4, Cxcl2, Ccl3, Ccl2, Tnf, Cxcl1, Cxcl10, Ccl5), and appears similar to a population described as “inflammatory macrophages” by Cochain et al. (18) (herein referred to as “chemokine^hi macrophages”). Clusters 6 and 9 are closely related to each other, with each having similar frequencies of cells in the progression and regression groups (Figure 2C). They share expression of Cxcl16, Atox1, CD72, Glipr1, but cluster 9 cells also express at higher levels of Trem2, Cd9, Lgals3, Spp1, Aldoa (Figure 2D, Figure 3, and Supplemental Figure 3B). Hence, cluster 9 is similar to the population of “TREM2^hi macrophages” described in ref. 18 (herein referred to as “Trem2^hi macrophages”). Although clusters 6 and 9 are closely related, there are also distinct differences, for in cluster 6, the AA467197 gene (encoding normal mucosa of esophagus-specific gene 1 protein [NMES1]) is highly expressed (herein referred to as “NMES1^hi macrophages”).

Figure 3

Classifications and gene lists of clusters of plaque macrophages derived from CX3CR1⁺ monocyte precursors in atherosclerosis progression and regression. Ranking of top 20 genes significantly overexpressed in each Louvain cluster, determined by statistical testing of one versus the rest with overestimated-variance t test. Mean of row-wise Z-score in the indicated cluster (y axis) is compared with bulk mean of Z-score for the ranking (x axis).

These results indicate that the 3 macrophage populations described by Cochain et al (resident-like, inflammatory, and Trem2^hi) are all derived from Cx3cr1⁺ monocyte precursors, but they represent only a subset of the macrophage populations we isolated from the aortas. Additionally, since these macrophage populations are present in both progressing and regressing plaques, they represent general inflammatory features of atherosclerosis.

Transcriptional profile of cells more abundant during progression. We next focused on the cells in the distinct clusters 0 (1%), 3 (5%), 5 (9%), and 10 (14%) that contained CX3CR1⁺ cells mostly from the progression group (Figure 2C). Cluster 5 in particular expressed a strong type 1 IFN signature, including IFN regulatory factor 7 (IRF7) and IFN-stimulated gene 15 (ISG15), as well as myeloid cell nuclear differentiation antigen (MNDA) and IFN-induced transmembrane protein 3 (IFITM3) (herein referred to as “IFN signature^hi macrophages”) (Figure 2D, Figure 3, and Supplemental Figure 4). Notably, type I IFN and myeloid type I IFN signaling has been shown to accelerate atherosclerosis progression in Ldlr^–/– mice by promoting chemokine-dependent leukocyte recruitment (26). IFN signature^hi macrophages also exhibit increased expression of Ly6e and Ly6a, which encode for the Sca-1 antigen. With an anti–Sca-1 antibody, by FACS we can identify a distinct pattern of Sca-1 expression in a small subset of cells that are TdTomato⁺ (Supplemental Figure 6, A and B). Additionally, we verified that CD9 (enriched in Trem2^hi macrophages) was highly expressed in a cluster of TdTomato⁺ cells (Supplemental Figure 6C). The CD9-expressing cells were clearly distinct from cells that express high MHC class II (MHCII) levels (Supplemental Figure 6D), representing the cluster of CD74^hiMHCII^hi macrophages. Hence, we could verify several distinct cell clusters from TdTomato⁺ gated cells by antibody staining, as anticipated from the single-cell RNA-Seq results.

In cluster 0, expression of the endonuclease DNase1l3 is of interest and may also be associated with an IFN response (herein referred to as “DNase1l3^hi macrophages”). However, genes expressed in cluster 3 were surprising, because Retnla and Ear2 are associated with exposure to type 2 cytokines such as IL-4 (herein referred to as “Retnla^hiEar2^hi macrophages”) (Figure 2D, Figure 3, and Supplemental Figure 4).

Hence, there is considerable heterogeneity of macrophage activation states under progression conditions, including the M1 features of a IFN type 1 signature, as well as the M2 features of an IL-4 signature. Our previous hypothesis that plaque progression is dominated by inflammatory M1 macrophages while plaque regression is dominated by M2 macrophages based on immunostaining and laser-capture plaque analyses (3, 15, 16) is likely too simplistic, as there was a clearly distinct IL-4 activated macrophage population (Retnla^hiEar2^hi macrophages) in progressing plaques that was actually absent from regressing plaques. These data are more consistent with the concept of a spectrum of macrophage activation states (27–29) than a strict dichotomy between M1 and M2 macrophages.

Use of diffusion pseudotime (DPT) analysis can measure transitions in gene expression between cells, reconstruct cell developmental progression, and identify cell branching decisions and differentiation at a single cell level (30). Since CX3CR1 is downregulated as monocytes differentiate into macrophages, we chose from the cell population with highest average expression of CX3CR1 (cluster 7) (see below) and then within this population, we chose cells with the highest individual CX3CR1 expression to be the “root” cell for DPT analysis to predict monocyte development (Figure 4, A and B). This approach identified Retnla^hiEar2^hi macrophages as having a trajectory separate from the main group comprising the other cell clusters (Figure 4, A and B), and these Retnla^hiEar2^hi macrophages were predominantly from progressing plaques (Figure 2, C and D, Figure 3, and Supplemental Figure 4). When we compared these data with our previous microarray data set (31), we noted that Retnla was also reduced in expression when comparing plaque regression with progression (Figure 4C). Since Retnla is highly expressed in macrophages exposed to IL-4, this initially surprising finding most likely reflects that compared with the regressing plaque, myeloid cells in the progressing plaque may be exposed to higher concentrations of both type 1 and type 2 cytokines.

Figure 4

Diffusion pseudotime and principal component analysis identification of genes associated with atherosclerosis progression and regression. (A) Diffusion pseudotime (DPT) analysis identified a cellular branching point only present in the progression group (blue) for the Retnla^hiEar2^hi cell cluster 3 (red). (B) The population with the highest expression of CX3CR1 was selected as the “root” cells to perform the DPT analysis to predict differentiation of CX3CR1⁺ cells. (C) Merged gene expression log₂ (polyI:C/saline) values from a previous atherosclerosis regression/progression study in Reversa mice (31) compared with gene-level log₂(regression/progression) values from supervised analysis of the single-cell data for 27 differentially expressed genes, which confirms that Retnla expression (red box) is most negatively associated with regression in both data sets. (D) Principal component analysis (PCA) reveals a smooth transition of cells from progression and regression groups along the PC1 axis. (E) Heatmap of the genes with the greatest loading factors for the PC1 axis.

It was striking that 3 of the distinct clusters of activated macrophages (DNase1l3^hi macrophages, Retnla^hiEar2^hi macrophages, IFN signature^hi macrophages; Figure 2A) were highly enriched in cells from progressing plaques (Figure 2C), indicating that during progression, macrophages become more differentiated and activated than during plaque regression. This suggests that the environment of a regressing plaque for CX3CR1⁺ cells may be less complex. These clusters were not noted by the previous single-cell reports (17, 18), likely because we focused on cells derived from CX3CR1⁺ precursors and sampled approximately 5 times more cells.

Transcriptional profile of cells more abundant during regression. For the regression specific clusters, cluster 2 represented only a small fraction (<1%) of total cells, but was particularly interesting because of the high expression of B cell–associated genes, including the early B cell factor 1 (Ebf1) and B cell antigen receptor complex-associated protein α chain (Cd79a) (herein referred to as “Ebf1^hiCd79a^hi macrophages”). Cluster 8 is an important contributor to regression (31%) and less so to progression (5%) (Figure 2C). Genes differentially expressed in this cluster include stabilin-1 (Stab1) and selenoprotein-1 (Sepp1), which are found in human and murine atherosclerotic plaques, with the former proposed to enhance efferocytosis and the latter to be antiinflammatory, both putative M2-associated functions (32). The heat shock protein (HSP) genes Hspa1a and Hspa1b were also differentially expressed in cluster 8 (Figure 2, C and D, Figure 3, and Supplemental Figure 4), indicating a protective role for HSPs in atherosclerosis regression (herein referred to as “HSP^hi macrophages”). A protective role is also suggested by the reduced circulating HSP70-1 (Hspa1a) and HSP70-2 (Hspalb) associated with atherosclerosis progression and heart failure (33, 34).

Genes associated with cells from progression and regression states. During visualization by principal component analysis (PCA), we noticed a smooth transition between cells from the progression and regression groups along the principal component 1 (PC1) axis (Figure 4D). Genes with the largest-magnitude loading factors for PC1 (e.g., CD74 and MHCII molecules) were associated with cells in cluster 10 (Figure 2D, Figure 3, and Figure 4E), which contained cells more abundant in progression than in regression (herein referred to as “Cd74^hiMHCII^hi macrophages”).

Since we observed a smooth transition between cells from the progression and regression samples along the PC1 axis, we used the PC1 score as a “pseudotime” measure to determine the distributions of cells from the 2 sample groups. While overlapping, the distributions were clearly shifted (Figure 5A). We screened for genes whose cell-specific expression levels were correlated or anti-correlated with the cell-specific pseudotime values; this unsupervised analysis yielded 42 genes whose increased expression levels were associated with regression and 7 genes whose expression levels were associated with progression (Figure 5B). Among the 42 genes associated with regression, the analysis highlighted increased gene expression of Cxcr4 (encoding CXC chemokine receptor type 4), Wfdc17 (encoding WAP four-disulfide core domain 17), Serpinb6a (encoding serine peptidase inhibitor, clade B, member 6a), Grn (encoding granulin), Ctsb (encoding cathepsin B), and Ctsd (encoding cathepsin D), whereas increased expression of the genes Cd74 (encoding cluster of differentiation 74), H2-ab1, H2-eb1, H2-aa (encoding MHCII molecules H2-Ab1, H2-Eb1, and H2-Aa), and Malat1 (metastasis-associated lung adenocarcinoma transcript 1–long noncoding RNA) are associated with progression.

Figure 5

Pseudotime scores and related genes associated with atherosclerosis progression and regression. (A) We used PC1 scores as a pseudotime measure and show shifted distributions of progression and regression cells. (B) Heatmap of genes identified to be significantly correlated with the pseudotime score. (C) We identified 53 genes significantly correlated with the pseudotime score that are negatively associated (r = –0.347, P = 0.01103) with merged gene expression log₂ (polyI:C/saline) values from a previous atherosclerosis regression/progression study (38).

In terms of known expression or functions of these genes in atherosclerosis, increased CD74 expression is observed in human atherosclerotic plaques and Malat1 is associated with atherosclerosis-related inflammation and lipid metabolism, although a distinct role for Malat1 in atherosclerosis has not been described (35, 36). Ctsb and Ctsd (37) are associated with regression and have not been linked previously to M2 macrophage features (Figure 4E and Figure 5B).

We also compared the gene-level pseudotime values with regression/progression expression ratios for differentially expressed genes in plaque CD68⁺ cells measured in our previous study of an aortic transplant–based model of plaque regression (38). We found significant anticorrelation (Figure 5C) (P = 0.01, r = –0.35), as expected, since progressing cells have higher mean pseudotime than regressing cells (Figure 5A). Hence, these are features that are consistent between our previous studies and the current single-cell analysis.

Identification of a stem cell–like proliferative CX3CR1⁺ cell cluster during atherosclerosis. As noted above, the population of cells in cluster 7 retained more CX3CR1 expression (Figure 6A). These cells also had a transcriptional profile very distinct from those of the other macrophage populations (Figure 3 and Supplemental Figure 4). This population of cells was enriched for cell cycle genes (Figure 6B), similar to recent findings in cancer stem cells (39), indicating that they are highly proliferative (herein referred to as “stem-like macrophages”). Gene ontology analysis confirmed that “cell division” is an important pathway for this population of cells (Supplemental Figure 7). These cells also express monocyte and macrophage genes such as Cd14, Adgre1, Csf1r, and Cd68 (Supplemental Figure 3A). We performed Ki-67 staining on plaque sections and confirmed that there were indeed CX3CR1-YFP⁺TdTomato⁺ cells expressing Ki-67 (Figure 6C).

Figure 6

Identification of a proliferative “stem-like” cell cluster that retains CX3CR1 expression. (A) Violin plot and heatmap showing higher expression of CX3CR1 in cluster 7. (B) Violin plot and heatmap showing higher expression of cell cycle genes in cluster 7. (C) Representative immunofluorescence images of aortic root staining with Ki-67 (purple), EYFP (yellow), and TdTomato (red) from the Cx3cr1^{CreERT2-IRES-YFP/+}Rosa26^{fl-tdTomato/+} atherosclerotic mice in the progression group. Scale bars: 50 μm. (D) 138 genes that were expressed more than one read in 90% of cells in cluster 7 were compared with macrophages, monocytes, myeloid progenitors, stem cells, and DCs extracted from the ImmGen database. Cells from cluster 7 were randomly grouped into 4 subgroups (shown in black), and the median expression of each gene is shown. Heatmap analysis illustrates the clustering of similar cell types and genes based on normalized gene expression profiles.

We identified 138 genes that were expressed in these stem-like macrophages and compared the expression profile to cells in the ImmGen database. The expression profile of stem-like macrophages hierarchically clustered with myeloid progenitors (MLPs) and stem cells, indicating that they may have self-renewing properties consistent with the expression of cell cycle genes (Figure 6D and Supplemental Figure 5).

To confirm that the stem-like macrophages can be identified in an independent biological experiment, we next merged our data set to the recently published data set from Kim et al., in which the authors single-cell-profiled sorted CD45⁺ cells isolated from whole aortas of Ldlr^–/– mice fed a WD for 12 weeks (40). We performed MergeSeurat, aligned the data sets with canonical correlation analysis (CCA), and identified 11 clusters by Louvain clustering (Figure 7A). We then classified the cell populations with Single Cell Recognition (SingleR), a novel method for unbiased cell type recognition (41), and identified specific cell clusters as macrophages, monocytes, DCs, and innate lymphoid cells (ILCs) (Figure 7B). Most of the cells in the merged data set were classified as macrophages by SingleR, and we identified a distinct macrophage cluster (cluster 6) that expressed the highest level of cell cycle genes and thus exhibited stem-like properties (Figure 7, B–E). Notably, this cell cluster was present in a similar proportion (3.4%) in total sequenced cells (Figure 7D) from each data set (i.e., both progression and regression samples and Kim et al.’s data set), demonstrating that the abundance of this cell cluster is consistent in the different atherosclerosis models and settings. We also confirmed that most of the cells in this cluster maintained high expression of CX3CR1 (Figure 7F).

Figure 7

Validation of the presence of a “stem-like” cell cluster from an independent study of atherosclerotic mice. A recently published single-cell data set from Kim et al. (ref. 40; “Kim”) from sorted CD45⁺ cells isolated from whole aortas of Ldlr^–/– mice fed WD for 12 weeks was combined with our progression and regression data sets (Lin Progression, Lin Regression) with Seurat merging and aligned with canonical correlation analysis (CCA). (A) t-SNE visualization of the 8906 single cells from both data sets after Louvain clustering and colored by sub-cluster.(B) SingleR method was used for unbiased cell classifications of each sub-cluster against the ImmGen database and colored and labeled accordingly on the t-SNE plot. (C) t-SNE plot colored based on experimental group and data sets, showing that cluster 6 includes cells from all 3 experiments. (D) Cluster composition by cell numbers of experimental groups and data sets (red, Kim Ldlr^–/–; green, Lin Progression; blue, Lin Regression) in total single cells. (E) Violin plot and heatmap showing the highest expression of cell cycle (“Stemness”) genes in cluster 6. (F) Violin plot and heatmap showing the expression of CX3CR1 in cluster 6.

Hence, we identify a population of proliferating cells with a stem cell–like signature derived from CX3CR1⁺ precursors that are present in both atherosclerosis progression and regression models, as well as in an independent Ldlr^–/– mouse atherosclerosis experiment (40). This finding extends the current view that tissue-resident macrophages can proliferate in plaques (42).

The above results represent the introduction of single-cell RNA-Seq combined with genetic fate mapping to track the cellular states during the differentiation of CX3CR1⁺ cells into macrophages in atherosclerotic plaques during their progression and regression. In addition to confirming the considerable heterogeneity of macrophages in plaques reported by others (e.g., refs. 43–46), our approach has led to several observations: (i) atherosclerosis progression is associated with differentiation into more distinct macrophage states than during regression; (ii) the spectrum of macrophage activation states in both progressing and regressing plaques has greater complexity than the traditional definition of the M1 and M2 polarization phenotypes; and (iii) there is a cluster of proliferating CX3CR1⁺ monocytes with a stem cell–like signature in both progressing and regressing plaques.

Our first observation is supported in the simplest sense by identification of 4 cell clusters (DNase1l3^hi, Retnla^hiEar2^hi, IFN signature^hi, Cd74^hiMHCII^hi macrophages) enriched from progressing plaques versus 2 cell clusters (Ebf1^hiCd79a^hi, HSP^hi macrophages) enriched from regressing plaques (Figure 2C). In addition, pseudotime analysis suggested additional lineage branching for cells enriched in progression (Retnla^hiEar2^hi macrophage) rather than regression. We also found that 3 of the cell clusters enriched in progressing plaques (DNase1l3^hi, Retnla^hiEar2^hi, IFN signature^hi macrophages) were more transcriptionally distinct (i.e., more separated from the main macrophage cluster/cloud) from the other macrophage populations. This indicates that during progression, macrophages may become more differentiated and activated than during plaque regression. Notably, these distinct clusters were not reported by the previous single-cell reports, but embedded in our data are the 3 macrophage populations described by Cochain et al.: resident-like, inflammatory, and Trem2^hi macrophages (Folr2^hi, chemokine^hi, and Trem2^hi macrophages, respectively in our data). Hence, these 3 macrophage populations noted previously are derived from CX3CR1⁺ monocyte precursors but only represent a fraction of different macrophage activation states present in atherosclerotic aortas, which we were able to achieve by sampling approximately 5 times more cells and focusing our analysis only on cells derived from CX3CR1⁺ precursors. In studies combining cytometry by time of flight (CyTOF) (mass cytometry) with single-cell RNA-Seq, Winkels et al. found evidence to suggest that leukocyte heterogeneity is lower in healthy (vs. atherosclerotic) aortas (17). Our finding that macrophages derived from CX3CR1⁺ cells are less complex in regressing plaques is consistent with this suggestion, as during regression, the plaque microenvironment is improved by the reversal of hyperlipidemia. Alternatively, macrophages in progression could have more distinct activation states because of stronger environmental stimuli.

The greater activation of macrophages derived from CX3CR1⁺ cells in progressing plaques is consistent with atherosclerosis being an inflammatory disease. Even more interesting, perhaps, are the data concerning the spectrum of macrophage activation states in progressing and regressing plaques. While macrophages characteristic of the M1 (classically activated) and M2 (alternatively activated, inflammation resolving) polarization extremes have been found by immunohistochemical markers in human and murine plaques (47), only recently have dynamic changes in marker expression been described. As we first reported in transcriptome analyses of macrophages (CD68⁺ cells) from progressing and regressing plaques, there was an obvious enrichment in cells expressing arginase I and other markers associated with the M2 state (38). The origin of these cells was from circulating Ly6C^hi monocytes, and if STAT6-dependent M2-associated polarization was prevented, M2-marker expression enrichment and, importantly, plaque inflammation resolution, were blocked (14).

Thus, the simple scenario based on those studies is that when the environment of an atherosclerotic plaque became healthier, newly recruited monocytes skewed away from the M1 toward the M2 direction of polarization, resolved the inflammation, and promoted beneficial remodeling of the artery. The present results, using a more incisive approach than previously possible, indicate that there is a diverse spectrum of macrophage activation states in atherosclerosis progression and regression, and the previous association of M1 macrophages only with progression and of M2 macrophages only with regression may be overly simplistic. While macrophages with M2 features (cell surface expression of PD-L2 and CD301 and increased MRC1 expression) were more abundant during regression (Figure 1, D–F, Figure 2, and Figure 3), a distinct cluster of macrophages (Retnla^hiEar2^hi macrophages) expressing M2 signature genes was present in progression, and not regression (Figures 2 and 3). We should note that a broad distribution of M1- and M2-associated gene expression was also found in the macrophage populations classified by Cochain et al., with many displaying heterogeneous patterns within and overlaps between populations (18).

Turning to our third observation, monocytes are known to constitutively enter non-lymphoid organs and recirculate to lymph nodes in the steady state without differentiation to macrophages or DCs (48). These Ly6C^hiMHCII⁺ monocytes may play a role in antigen presentation or surveillance of inflamed or steady-state tissue (48–50). In the present studies, we identified a distinct cluster of cells (stem-like macrophages) with proliferative features that maintain CX3CR1 expression as well as markers of monocytes and macrophages, while sharing transcriptional profiles with stem cells. Macrophage proliferation in progressing plaques in Apoe^–/– mice and other models has been described (ref. 42; reviewed in ref. 6). The present results are also consistent with our finding in regressing plaques of a small population of cells positive for the proliferation marker Ki-67 (14).

The source of the proliferating macrophages could be either tissue-resident macrophages, inflammatory monocyte–derived macrophages (12) or even the proliferative CX3CR1⁺ monocytes that we have now identified in both progressing and regressing plaques (Figures 2 and 6). However, this cell cluster is relatively rare (Figure 2C), and the roles or fate of these proliferative tissue-resident monocytes remains to be elucidated. Nonetheless, uncovering this population raises the possibility that there are proliferating stem cell–like monocytes that can self-renew within inflamed tissues, serving as an additional reservoir for tissue macrophages that adopt different activation states depending on the microenvironments of atherosclerosis (Supplemental Figure 8).

In conclusion, we have significantly extended studies of atherosclerosis by a combination of a state-of-the art mouse model, fate mapping, transcriptomics, and bioinformatics. In addition to analyzing many more cells in progressing plaques compared with other studies (17, 18), we have obtained new data on regressing plaques. As a result, we have identified multiple new cellular phenotypes, with distinct molecular features, that are enriched in populations of macrophages in atherosclerosis progression and regression. We also found a number of macrophage populations with features that are shared between these disease states, which highlights the possibility that therapies based on data from progressing plaques alone may have unintended adverse consequences for the regression process. Further progress will be needed, therefore, to achieve the goal of identifying the targets with the best likelihood of therapeutic success by avoiding this pitfall.

See the Supplemental Methods for a detailed explanation of all experimental procedures.

Study approval. All animal procedures were approved by the NYU School of Medicine IACUC.

Statistics. All data are presented as mean ± SEM. Two-tailed unpaired Student’s t test was used calculated with Prism 7.0 from GraphPad. P < 0.05 was considered significant.

Data availability. The full data set is available in the NCBI’s Gene Expression Omnibus (GEO GSE123587).