A common framework of monocyte-derived macrophage activation

Macrophages in inflammatory conditions co-exist in diverse functional states

We aimed to contextualize macrophage heterogeneity across tissues in diverse inflammatory conditions, specifically defining common aspects of macrophage activation during infection and inflammation (fig. 1A). For this purpose we built a reference dataset (fig. 1A–B) based on 2 inflammatory conditions representing either a classical inflammatory response during bacterial infection using Listeria monocytogenes (L. mono), or a type-2 immune response during Heligmosomoides polygyrus (H. poly) helminth infection (fig. S1). We reasoned that these settings, which induce multi-cellular and systemic responses, provide a broad spectrum of macrophage activation improving the completeness of our reference dataset. We performed scRNA-seq on all stromal vascular fraction cells from mesenteric fat (fig. S2A), adjacent to the site of H. poly infection, and popliteal fat (fig. S2B), which is directly invaded by L. mono following footpad injection. We observed a rich infiltration of CD45⁺ cells expressing archetypical cytokines associated with type-2 and type-1 immunity (fig. S1). We then examined macrophage gene expression markers in these datasets (fig. S2A–B) and extracted, balanced (500 cells per condition), integrated and re-clustered identified macrophages. Gene expression was distinct within each resulting cluster (fig. S2C & data file S1) and associated with specific biological processes (fig. 1C & table S1).

We observed that clusters 0 and 3 were enriched for genes associated with macrophage alternative activation (e.g. Cd36, Clec10a, Mrc1), antigen presentation (e.g. H2-Aa, H2-Eb1, H2-Ab1) and the complement cascade (e.g. C1qc, C1qb) (fig. 1C, S2C & data file S1). Cluster 2 was enriched for genes involved in extracellular matrix (ECM) receptor-interactions (e.g. Cd44, Sdc1, Fn1) and cytoskeleton regulation (e.g. Pfn1, Actg1, Tmsb4x). Cluster 4 displayed high expression of genes participating in antigen presentation (e.g. H2-Oa, H2-DMb2, Cd74). Clusters 5, 6 and 7 were enriched for genes associated with the phagosome (e.g. Fcgr1, Ncf4, Fcgr3) and oxidative stress (e.g. Prdx5, Txn1, Gsr), with cluster 6 in particular enriched for innate immune response genes (e.g. Ifitm3, Fcgr1, Isg20). ECM organization genes (e.g. Col1a1, Col3a1, Ddr2) were highest in Cluster 8, while Cluster 9 displayed high expression of cell cycle associated genes (e.g. Cks1b, H2afx, Cks2). Cluster 1, which occupied the center of the UMAP, had no distinctly regulated genes (fig. S2C & data file S1) and therefore no pathway assignment (fig. 1C). Thus, our analysis shows that macrophages within a tissue simultaneously specialize into multiple functional categories, corroborating findings in other studies where this diversity has been reported (21, 23, 26).

Macrophages in inflammatory conditions are arranged along activation paths

Cluster 1 could not be associated with a distinct function as it displayed no up-regulated genes using the thresholds we established (average log fold change > 1, percent cells expressing gene > 0.4 and an adjusted p value < 0.01). We reasoned that cluster 1 could represent an intermediate state of activation, suggesting that rather than different populations of macrophages, our analysis captured activation paths being followed by infiltrating monocyte-derived macrophages. To test this hypothesis, we used Slingshot (28) to generate a model of macrophage activation in the tissue (fig. 1D–E). To establish a starting point for this analysis, we calculated a monocyte gene set score (fig. 1D & table S1) (1), then selected cluster 4 (fig. 1B–D, black arrow), as it displayed a monocyte-signature as well as increased Major histocompatibility complex class II (MHC-II) expression, in line with reported monocyte-to-macrophage transitions (29).

Our analysis identified 4 activation paths that we labeled as “phagocytic”, “oxidative stress”, “inflammatory” and “remodeling” according to the enriched pathway at the end point clusters of each lineage (fig. 1E). Our analysis also revealed three common stages of macrophage activation, which we named according to their relative position in the pseudotime progression as initial-, early- and intermediate-stage (clusters 4, 2, 1 respectively). The remaining clusters were renamed as either late.P#- or final.P#-stage (fig. 1F); for example, late.P1-stage corresponded to cells in the cluster in activation path 1 (P1) that is between the intermediate- and final-stage. We complemented our approach with two alternative strategies of data integration, specifically the Seurat (30) implementation of Harmony (31) (fig. S2D) and the scanpy (32) implementation of BBKNN (33) followed by a Partition-based graph abstraction (PAGA) (34) to infer activation trajectories (fig. S2E). Both approaches broadly recapitulated our initial results. Initial-stage cells were closely connected to those in early-stage, which then transitioned through an intermediate-stage that split into cells in the “phagocytic” (P1), “inflammatory” (P3) and “remodeling” (P4) paths. Notably, our clustering and labeling results were not greatly altered by alternative data processing approaches (fig. S2D–E), with the overall conservation of global structure and relationship between the proposed stages supporting further exploration.

Next, we investigated how these activation-stages were distributed across biological conditions (fig. 1G). We observed that the distribution in macrophages from naive mice into different functional stages in different fat deposits was comparable, and diverse, perhaps indicating an active process of macrophage activation as monocytes replenish these sites in steady state (fig. 1G). We noted that L. mono infection favored the “oxidative stress” path, while H. poly infection favored the “phagocytic” path when compared to each other (fig. 1G). However, all activation- stages were present in each condition, underscoring the limitations of bulk phenotyping to capture macrophage activation ex vivo. Moreover, the changes in diversity induced by infection were notably different between the inflammatory conditions studied. While H. poly infection resulted in increases in numbers of cells within less well-represented activation-stages, L. mono infection resulted in 2 activation-stages becoming dominant. These differences could be explained both by the specific immune responses tailored to the pathogens involved, but also by the time at which samples were collected (table S2; day 1 p.i, for the L. mono dataset, day 14 p.i. for the H. poly dataset).

In summary, our model predicts that macrophages in an inflamed tissue progress through activation-stages with distinct transcriptional profiles and that the balance of this progression is influenced by the type of immune response that dominates the microenvironment. Finally, our model suggests that tissue can become dominated by relatively few activation-stages, yet macrophages committed to multiple paths co-exist.

Macrophages transition through a RELMɑ-expressing activation-stage

We complemented our differential gene expression analysis (fig. S2C) by modeling gene expression as a function of pseudotime (fig. 2A), reasoning that transcript levels would be regulated as macrophages progress along activation paths (fig. 1E–F). We included in this analysis only the top 2000 most variable genes in the cells of each path, and ranked the resulting models based on the p value of the association of pseudotime and gene expression (fig. S3A). We found that collectively, the expression of 828 genes (p value < 1×10⁹) could be modeled in this way and we show the top most significant association for each pathway (fig. 2B). We observed that the regulation of typical macrophage activation markers corresponded to defined activation paths (fig. 2C & S3B). For instance, cells in the “phagocytic” path exhibited a steep and continuous increase in the alternative activation marker Mrc1 and the mitochondrial metabolism gene mt-Co1 as a function of pseudotime (fig. 2C). Conversely, expression of inflammatory Il1b was only retained at high levels in the “inflammatory” and “remodeling” paths (fig. 2C).

Exploring the relationship between RELMɑ encoding gene Retnla and pseudotime, we observed that our model predicted a wave of expression early during activation (fig. 2D). Indeed, our data revealed that Retnla was expressed both in cells committed to the “phagocytic” path (P1) (fig. 2D), which also expressed other markers of alternative macrophage activation (fig. 2C and data file S1), but also in the early-stage shared by all paths (fig. 2D). Moreover, we identified Retnla and Ear2 co-expression as unique to early-stage cells (fig. 2D). We reasoned that if our model was correct, then monocytes would express Retnla shortly after entering a tissue regardless of the inflammatory condition. To test this, we isolated bone marrow monocytes from CD45.1⁺ mice (fig. S3C), and adoptively transferred them into the peritoneal cavity of CD45.2⁺ naive hosts (fig. 2E). We evaluated the abundance of Retnla mRNA and RELMɑ protein in CD45.1⁺ macrophages recovered on days 2, 4 and 8 days after adoptive transfer (fig. 2F & S3D). Confirming our predictions, monocyte-derived macrophages expressed RELMɑ 4 days post-adoptive transfer and ~80% of the cells were positive for this molecule after 8 days (fig. 2F). RELMɑ induction occurred in the absence of IL-4Rɑ stimulation, as we did not detect changes in RELMɑ expression in host tissue-resident macrophages in the peritoneal cavity and IL-4Rɑ^−/− monocytes displayed similar RELMɑ expression compared to IL-4Rɑ^+/+ cells (fig. 2G). Finally, mature CD45.1⁺ peritoneal macrophages did not express RELMɑ after transfer into CD45.2⁺ naive recipients (fig. 2G). Thus, our experimental data support the view that monocytes differentiating into macrophages early after entry into a tissue begin to express RELMɑ independently of IL-4 signaling.

Intriguingly, macrophage specific RELMɑ deficiency leads to poor monocyte/macrophage infiltration of the lung in a model of fibrosis (35), which suggests that RELMɑ might be required for macrophage tissue engraftment. To test this, we isolated monocytes from the bone marrow of RELMɑ^−/− animals (36) then adoptively transferred them into the peritoneal cavity of CD45.1⁺ hosts and evaluated their ability to engraft compared to wild type macrophages (fig. 2H). We recovered fewer RELMɑ^−/− cells 8 days after transfer (fig. 2H). Furthermore, the proportion of CD115⁺F4/80⁺ macrophages was significantly reduced in RELMɑ^−/− vs wild type controls. These results indicate that RELMɑ serves as a marker of macrophage tissue engraftment, and may also function to maintain the identity of the infiltrating cells, as CD115 (CSF1R) is a critical macrophage growth factor receptor. These findings support our predicted model of macrophage activation and explain evidence of historical Retnla expression in tissue-resident macrophages (37). Moreover, RELMɑ expression has been used to identify cells of a distinct tissue-restricted phenotype (21) or in an immature state of differentiation (8). Our data is broadly in agreement with these past reports and indicates that historical RELMɑ expression should be evident in tissues where circulating monocytes replace tissue-resident macrophages. Moreover, our data suggest that transit via a RELMɑ⁺ stage is a common feature of all macrophages and not restricted to a single tissue or macrophage subset.

Macrophage activation-stages are conserved across tissues and inflammatory conditions

Our findings in the adipose tissue datasets indicate that observed heterogeneity in macrophage activation can be accounted for by macrophages transiting through defined activation paths after they enter a tissue. To evaluate these ideas further, we used our adipose tissue dataset as a reference to interrogate multiple other situations of tissue inflammation (table S2). We took advantage of a recent data transfer implementation that enables the construction of large cell atlases (38), assigning probability-scored labels and performing data imputation (i.e. inferring and assigning missing gene expression values). This strategy is highly accurate for label transfer from a reference dataset (39), yet can over-correct batch effects while sacrificing biological variability (40). We tested this approach on datasets which contained a mixture of macrophages and other CD45⁺ cells, using the adipose tissue as a reference and extensively optimizing label transfer/imputation parameters, reasoning that only macrophages should have high probability scores resulting from transfer of activation-stage labels (fig. S4A–H). Using a probability threshold of 80% (or 0.8), almost exclusively macrophages were assigned an activation-stage label in datasets with abundant (fig. S4A–D) and scarce macrophage populations (fig. S4E–H). Thus, we proceeded to interrogate multiple tissues and inflammatory conditions using this approach.

We retrieved publicly available scRNA-seq datasets spanning 9 tissues and 13 inflammatory conditions plus healthy controls, including infections, injuries, cancer and dietary interventions (fig. 3 & table S2). We retained only macrophages by filtering cells based on established markers (table S1), harmonized the data within each tissue to remove batch effects and applied the transfer process, assigning cells to a previously defined activation-stage, with a label probability (fig. 3) and imputed gene expression data (fig. S4I–Q). Imputed data was used to cluster and calculate a UMAP for each dataset (fig. S5). Where a single label dominated a cluster, and a portion of the cells within the cluster had a high label probability score, the corresponding label was assigned to the entire cluster (e.g. fig. S5B, clusters 3, 4, 6 & 7). Macrophages that did not meet this criteria were marked as “not classified”. We identified most of the activation-stages defined in our reference with a reasonable proportion of cells with a high label transfer probability in all interrogated datasets (fig. 3).

The distribution of activation-stage labels was different in each studied tissue and then modified by the corresponding inflammatory condition (fig. 3). This could reflect the influence of the tissue micro-environment in shaping the immune response, as well as the way in which the immune response is tailored to a specific insult. For example, in the large intestine lamina propria (fig. 3A), infiltrating initial-stage cells were abundant in steady state (fig. 3A, bottom right) in line with the reported turn-over of macrophages in this tissue (9). However, after 12 weeks of high fat diet (HFD) the proportion of initial-stage macrophages diminished, being replaced by cells in the final.P2-stage of the “inflammatory” path (fig. 3A), which is consistent with HFD-induced gut inflammation (20), and TNF production by intestinal macrophages (41). This correspondence between our labeling strategy and established biology could be seen in all datasets. In sciatic nerve injury (21) a wave of “inflammatory” path cells (final.P3) could be seen 1 day post-wounding (dpw) which receded by day 5, when “phagocytic” path final.P1-stage cells took over (fig. 3B). In breast tumors (22), the macrophage landscape was dominated by “phagocytic” path cells (fig. 3C), which displayed markers of alternatively activated cells. The same appeared true in regressing atherosclerotic plaque lesions (23) (fig. 3D) and in liver fibrosis (24) (fig. 3E), while fungal infection in the lung (25) resulted in an increase in “inflammatory” path cells (fig. 3F). In infarcted heart (26) and retinal damage (27), an expansion of “phagocytic” path cells was also evident (fig. 3G–H), although the diversity of activation-stages in each tissue was strikingly different, with few identified stages in the retina both in steady state and after light-induced neurodegeneration (fig. 3H). In contrast, skeletal muscle macrophages (18) displayed diverse activation-stages, with chronic parasitic infection having a modest effect on the stage distribution in this tissue (fig. 3I), although increased “inflammatory” path macrophages were apparent (Final.P3). Finally, we observed early-stage cells, co-expressing Retnla and Ear2, in nearly all analyzed datasets (fig. 3A–D, G & I), underscoring how this activation step appears to be common to macrophages in most tissues.

Despite demonstrable utility of our labeling approach, there remained a number of cells with no label assignment (i.e. “not classified”). We developed 2 potential explanations for the abundance of these cells in the studied datasets. First, our approach hinges on stringently identifying anchor pairs between the data, obtaining a high score and/or having dominant labels in the clusters. Consequently, we are less efficient at identifying transcriptional profiles of cells in between defined activation-stages. For this reason, many of the “not classified” cells in our analysis could be seen in between labeled clusters in the UMAP, and would often share low probability labels for flanking clusters with a more defined signature (e.g. fig. S5D, cluster 2 flanked by 1, 4, 7 & 9). Likely for this reason, intermediate-stage cells were relatively rare in query datasets, as these were the least defined transitional state that we identified in our reference data. Second, tissue-resident macrophages have a distinct transcriptional profile (4, 5), and thus might not easily relate to activated macrophages originating from circulating monocytes in inflammatory settings. Indeed, we observed more unclassified cells in tissues where monocyte infiltration is rare (Retina - fig. 3H) or where specialized macrophage subsets are common (alveolar macrophages in the Lung - fig. 3F). In fact, the distinct cluster of “not classified” cells present on the left of the Lung UMAP (fig. 3F) was enriched for alveolar macrophage markers, thus explaining the striking bimodal label distribution in this dataset (fig. 3F). By contrast, the label assignment in the atherosclerotic plaque dataset was nearly global (fig. 3D) likely as only circulating monocytes-derived macrophages were studied (23).

We explored the issue of macrophage embryonic origin and tissue immune privilege in more detail by studying a dataset where microglia were recovered at different developmental steps from naive mice (fig. 4A and table S2) (42). In line with our expectations, we observed poor label probability distributions for all investigated ages (fig. 4A). Interestingly, the probability threshold was never surpassed and indeed these distributions were skewed progressively towards 0 as the age of the investigated animals increased.

Taken together, our label transfer analysis shows that macrophages across tissues and inflammatory conditions share common transcriptional profiles that correspond to definable activation paths. Our analysis also suggests that embryonically seeded and highly specialized tissue resident macrophages do not respond to inflammatory conditions in a way analogous to that of infiltrating monocytes, with the latter encapsulating most of the functional diversity found in all the studied datasets. Consequently, our model best captures the activation of monocyte-derived macrophages that replace tissue resident cells during inflammation. Finally, well-established paradigms of macrophage biology are reinforced by the functional stages we defined, making these labels a potential tool to probe deeper into the functional specialization of macrophages during inflammation.

Exploiting the predictive nature of the proposed macrophage activation model

We decided to explore in greater detail potential biological insights that might be obtained from our activation model. For this purpose we turned to the atherosclerotic plaque and breast tumor datasets (fig. 4), both of which included interventions that ameliorated disease progression (table S2) (22, 23). Dietary and pharmacological intervention (table S2) were reported to induce regression of atherosclerotic plaque lesions (23), which are dominated by “phagocytic” path macrophages (fig. 4B). Strikingly, there was a shift between late.P1- and final.P1-stage cells in regressing lesions, with a sizable decrease in the proportion of late.P1-stage cells (fig. 4C). We investigated which genes were altered in expression within cells in this activation-stage between progressing and regressing lesions. We performed this analysis on the original unimputed expression data to enable our label predictions to orient the analysis without affecting the underlying measurements. We compared these regulated genes with those associated with the “phagocytic” path based on our pseudotime analysis (fig. 2). Interestingly, expression of all but one of the genes in this activation-stage increased in regressing lesions (fig. 4D). Moreover, all but one of these genes tended to be expressed more strongly as cells progressed from the late.P1- (fig. 4E - dashed lines) to the final.P1-stage. Collectively, this data suggest that the intervention causing lesions to regress, induced the accelerated transit of macrophages along the “phagocytic” path, as indicated by the increased expression of genes associated with this trajectory in late.P1-stage cells, which concomitantly decrease in abundance.

Next, we examined the breast cancer dataset, where macrophage-specific Dab2 deletion was reported to dampen tumor progression (fig. 4F–J & table S2) (22). We observed that Dab2 was most highly expressed in “phagocytic” path late.P1- and final.P1-stage macrophages (fig. 4F–G), with the former being more abundant. In fact, we observed an increase in late.P1-stage Dab2 deficient macrophages (fig. 4H), leading us to examine differentially expressed genes between these cells and their WT counterparts. Our results show that from 15 regulated genes also included in the pseudotime analysis, 11 were differentially down-regulated between these two groups (fig. 4I). Of these 15 genes, 13 were highly expressed at this stage of activation in our reference data (fig. 4J - dashed lines). The downregulation of these path-associated genes and the accumulation of late.P1-stage cells could indicate that the absence of Dab2 stalls the progression of macrophages in the “phagocytic” path.

Monocytes enter wounds populating the functional diversity predicted by the proposed macrophage activation model

Our model predicts that incoming monocytes are able to assume the phenotype of existing macrophages in the tissue and populate all functional stages described. Indeed, in the atherosclerotic plaque dataset, all sequenced cells were derived from circulating precursors (23). To formally evaluate this, we traced the influx of monocytes into wounded skin (fig. 5A–B & S6). Red fluorescent monocytes (tdRFP⁺) were administered i.v. 2 or 12 dpw, and all wound macrophages were harvested 4 and 14 dpw (fig. 5B), using index sorting to retain fluorescence values from barcoded cells for further analysis. We mapped, clustered (fig. S6A) and labeled the sequenced cells (fig. 5A) as described above, applying similar thresholds and evaluating label probability distributions for the entire dataset (fig. 5C) and across clusters (fig. S6A–B). Wounded skin exhibited most previously defined activation-stages, again demonstrating the robustness of our activation model.

We observed a global label probability distribution skewed towards a value of 1 (fig. 5C), indicating a good agreement with our reference data. Moreover, label distribution across conditions was consistent with expectations based on established literature (43), with an early wave of inflammatory cells at 4 dpw (fig. 5D) and a later increase in “phagocytic” path cells at 14 dpw (fig. 5D). Transferred monocytes were detected at both 4 and 14 dpw (fig. S6C), although only when given on day 2 (fig. S6D). Infiltrating monocytes were distributed across all detected clusters and assigned stage labels (fig. 5E–G) mirroring closely the distribution for all macrophages sequenced (fig. 5D & G). As predicted by our model, fluorescent cells were assigned to the initial-stage and “inflammatory” path 4 dpw (fig. 5G), emphasizing the transitory nature of the early inflammatory wave which occurs during tissue repair. By contrast, transferred monocytes mapped preferentially to the early-stage and “phagocytic” path 14 dpw (fig. 5G). The fact that transferred cells only mapped to the initial-stage on day 4, and that by day 14 they had assumed other identities, reflects the number of days since these cells accessed the wound and supported our model, showing monocyte flux through distinct activation-stages. More early-stage macrophages co-expressing Retnla and Ear2 were detected at 14 dpw (fig. 5G–H) compared to day 4, despite the presence of “inflammatory” path cells at this time point that, according to our model, should have transited through the early-stage. Previously, 2 days were insufficient to observe Retnla expression in infiltrating cells (fig. 2F), and increased proportion of early-stage cells at 14 dpw (fig. 5G) aligns with our earlier findings that after 8 days a higher percentage of infiltrating cells would fall in this stage (fig. 2F). To reconcile these observations, we propose that the relative speed at which macrophages traverse defined activation paths is influenced by the inflammatory conditions at the site of immunological insult, so that cells may accelerate their passage through identified stages to better adapt to the required immune response.

Next, we took advantage of the indexed nature of the skin wound dataset to validate our labeling strategy with common macrophage phenotypic markers. For this purpose, we scaled the fluorescent signal of CD301b, CD45, F4/80, CD11b, MHCII and Ly6C, as well as the side and forward scatter parameters, across all cells and calculated a UMAP from this flow cytometry data, retaining the stage labeling based on the transcriptional profile of the cells (fig. 5I). We then identified clusters using k-means (fig. 5J & S6E). Strikingly, we observed that late.P1-, final.P1- and final.P3-stage cells each dominated a cluster (fig. 5L - clusters 4, 5 and 6 respectively), while early-stage cells localized predominantly to cluster 2 (fig. 5J&L). As expected, transferred monocytes were evenly distributed across these clusters (fig. 5K & S6F), further emphasizing the ability of these cells to differentiate into all identified functional stages. Finally, we show that these flow cytometry-based clusters are associated with significantly different protein expression levels (fig. 5M - p value < 0.001). Collectively, these data show that monocytes entering a wound respond to this environment by becoming activated in accordance with the model proposed in this study and provide surface protein validation of these stages.

Macrophage activation-stages in tissues have distinct transcriptional markers

Our approach relies on the quality of anchor pairs identified across datasets to transfer labels and impute gene expression data. In this manner, we identified all activation stages in different proportions across all the datasets studied (fig. 6A). As this approach by necessity transforms the original expression data, we sought to determine its robustness by interrogating labeled cells directly (fig. 6B). For this purpose we took all high label probability cells (>80%) from 10 query datasets, and combined these with a randomly sampled portion (n= 500) of macrophages from our reference, retaining only the label assignment and original uncorrected gene expression data (fig. 6B). Once extracted, these macrophages (n = 2843) were integrated across tissues, without giving priority to any dataset, then clustered and visualized as a UMAP (fig. 6C).

Our expectation for this analysis was that the transcriptional tissue signature would not obscure the activation-stage label. That is, our model would predict that tissues would not define the resulting cell clustering, but rather that the activation-stage of these cells would be sufficient to group them in this unsupervised analysis. Critically, this was the outcome we observed (fig. 6D). Cells with identical activation labels clustered together, regardless of the tissue of origin or the inflammatory condition, demonstrating that our initial approach was valid. We next performed a complementary trajectory inference analysis on this data using Monocle (44), which organized cells into two main branches of macrophage activation, chiefly one directed to the phagocytic path (P1) and one towards the inflammatory pathways (P2 and P3) (fig. S7A). Although this analysis is unable to resolve the differences between P2 and P3, the overall architecture of the projection is reminiscent of our reference data. Indeed, cells in the initial-, early- and cycling-stage clustered near each other, while intermediate-stage macrophages separated the bottom “phagocytic” and “remodeling” paths from the “oxidative stress” and “inflammatory” paths at the top of the UMAP (fig. 6D). This distribution of clusters adds weight to our proposed activation trajectories and emphasizes the relative relationship between the activation-stages defined.

Having all macrophages clustered in this manner, allowed for the extraction of tissue-independent transcriptional markers for all activation-stages. We found upregulated genes (fig. 6E & data file S2) corresponding to each identified activation-stage, which largely aligned with the genes we originally associated with each label (fig. S2C & data file S1). Moreover, we extracted from this result potential cell surface markers (data file S3) and selected those with commercially available antibodies (fig. 6F & S7B). We reasoned that combining these markers would enable us to identify macrophages at different activation-stages. We injected the peritoneal cavity of mice with 10⁶ CFU of L. monocytogenes and tracked macrophage recruitment 1, 4 and 8 days post-infection (fig. 6G). The resulting fluorescence intensities for the macrophage populations were scaled using a linear regression model and used to cluster cells and calculate a PCA and UMAP for visualization (fig. 6H). These scaled fluorescence values were relative to the cells retained in the analysis, so rather than positive or negative, we used high (hi) or low (lo) to score protein expression. This analysis was able to capture a high degree of heterogeneity in the population of macrophages in the peritoneal cavity (fig. 6H & S7C), with cells originating from different time points after infection clustering separately. Next, we assigned our model’s activation-stage labels to the identified clusters based on the distinct staining patterns we observed in the data (fig. S7C). For some markers, we found no correlation between transcript and protein levels (notably Nrp1 - CD304, Fcgr4 - CD16–2, Cd14 - CD14, Tnf - TNF), however, we were reasonably confident of the annotation based on the remaining markers and identified at least three proteins that are distinctly expressed in each of the identified activation-stages (fig. 6H - right & fig. S7C).

Next we used the assigned labels to explore the population dynamics across our experimental peritonitis model (fig. 6I). Our annotation was consistent with the population dynamics that our model would predict and are expected of an inflammatory process (fig. 6I). Indeed, we observed a wave of initial-, early- and late.P2-stage macrophages replacing the tissue resident cells shortly after infection, followed by a progressive appearance of cells in the intermediate-stage, and lastly an accumulation of final.P2-, final.P3-, late.P1- and final.P1-stage cells (fig. 6I). Notably, RELMɑ expression was greatly enriched in the early-stage macrophages, which were most abundant 4 days post-infection (fig. 6H). These findings provide further support to our proposed model and the identified markers may be used by researchers to explore macrophage heterogeneity in other contexts.

Transcriptional network analysis reveals information hubs and upstream regulators

Earlier we hypothesized that changes in the regulation of genes associated with defined activation paths was likely to stall or promote macrophage activation. To explore this final aspect of the data we returned to the pseudotime-regulated genes (fig. 2), built then filtered a network based on known protein-protein interactions and gene expression regulation along identified paths (fig. S8A–B). Next, we clustered the resulting network, performed GO term enrichment analysis in each cluster (fig. S8C) and labeled the network according to the most enriched term in the gene set. The resulting network (fig. 7A, 242 genes) showed 10 clusters associated with distinct functions. These gene sets were not engaged similarly by all activation paths, as revealed by subsetting the network to include only genes significantly associated with pseudotime in each path (fig. 7B).

We extracted three types of information from this transcriptional network. First, we examined highly connected nodes, which articulated the network by connecting 2 or more clusters and were overrepresented in the paths connecting pairs of nodes in the network (i.e. high betweenness). Examining genes meeting these criteria (fig. 7C) highlighted both some well-known macrophage regulators (e.g. Lyz2, Csf1r), but also genes whose function has not been widely studied in the context of macrophage activation (e.g. Gngt2, Srgn), thus warranting further exploration. Second, we identified transcription factors (TF) upstream of the transcriptional network clusters (fig. 7D), of which six (Rel, Maf, Fos, Jun, Spi1, Jund) were themselves regulated dynamically along the activation paths. Interestingly, some TF had opposing behaviors (Rel vs. Maf), while others behaved similarly in all paths except one, where they suddenly veered in opposing directions (Spi1 vs. Fos/Jun). Finally, we took the top 3 GO enriched terms in each cluster (fig. 7E & S8C) and all enriched KEGG pathways detected (fig. 7F) and calculated gene set scores for each of these (table S1). We then estimated the variance of each gene set score within the activation paths and represented these data as a heatmap (fig. 7E–F), finding a distinct profile for each path for these functions. In this manner, we highlight the relative regulation of several processes of interest in the activation paths we defined, as a guide to researchers wishing to explore these aspects of macrophage function in more detail.

In summary, we employed a predictive model of label transfer to encompass macrophage activation irrespective of tissue or inflammatory condition. We demonstrate that this model is robust, aligning with well-established paradigms of macrophage function, while providing original avenues for investigation. We provide surface and global gene expression profiles for these activation-stages to aid in their identification in future studies. We have prepared an online tool (https://t.jh.edu/macrophage-framework) to aid in exploring the data contained in this study. Our results emphasize the conservation and relative homogeneity of macrophage activation across tissues, transcending macrophage tissue residence, while still allowing for activation diversity.

Study design

The objective of this study was to integrate the transcriptional profile of macrophages retrieved from diverse inflammatory conditions and in steady state. We used single cell RNA sequencing technologies and flow cytometry to study cells ex vivo. Tissues were obtained from animals infected with several pathogens, wounded or fed experimental diets. We adoptively transferred monocytes and phenotyped these and other macrophages after specified times. Collected data were analyzed with bioinformatic and computational tools with the goal of developing a unified framework for monocyte-derived macrophage activation. Groups sizes and endpoints were determined based on experimental purposes using previous experience or estimated based on preliminary data and expected power of 95%. No outliers were detected in the data. Age and sex-matched mice were assigned randomly to experimental groups. This study was not blinded. Sampling and experimental replicates are indicated in each figure legend.

Mouse Models

C57BL/6J (RRID: IMSR_JAX:000664), B6.129P2-Lyz2^{tm1(cre)lfo/J} (RRID: MGI:5014089), CD45.1 congenic (RRID: IMSR_JAX:002014) and C57BL/6N (RRID: MGI:5882504) mouse strains were purchased from The Jackson Laboratory. IL-4Rɑ^−/− mice were generated from B6.129P2-Lyz2^{tm1(cre)lfo/J} and B6-Il4ra^tm(loxp) (65). Mice were maintained at the Max Planck Institute for Immunobiology and Epigenetics or at The Bloomberg-Kimmel Institute for Cancer Immunotherapy at Johns Hopkins. C57BL/6JRj (RRID: MGI:2670020) and R26^LSL-tdRFP animals (66) mice were maintained at the Experimental Center of the Medical Faculty, TU Dresden. Bone marrow from RELMɑ deficient animals was kindly donated by Dr. Meera Nair from University of California Riverside. Animal care was undertaken in accordance with Institutional Animal Use and Care Guidelines with approval by the corresponding animal care committee. All animals used for tissue harvest or experimental procedures were aged between 6–9 weeks at the start of the experiment. Animals were humanely sacrificed by carbon dioxide asphyxiation followed by cervical dislocation and tissue dissection. Mice were bred under specific pathogen free standards.

Experimental infections and diets

L3 infectious stage Heligmosomoides polygyrus (H. poly) larvae were kindly provided by Dr. Joseph Urban Jr, USDA, ARS, Beltsville Human Nutrition Research Center, Diet Genomics and Immunology Laboratory, Beltsville, USA. To induce H. poly infection mice were gavaged with 200 L3 infectious stage larvae in PBS. Mice were left for 13 days before being sacrificed. A wild type strain of Listeria monocytogenes (L. mono) was used for infections. Mice were infected subcutaneously on the footpad with a sublethal dose of 1 × 10⁶ colony-forming units (CFU). Mice were left for 1 day before being sacrificed. Alternatively, mice were injected intraperitoneally 1 × 10⁶ colony-forming units (CFU) and left for 1, 4 or 8 days before being sacrificed. Obesity was induced by ad libitum feeding of C57/BL6 mice for 12 week with irradiated high fat diet (Rodent Diet 60% kcal from fat, Research Diets, Inc., Cat #D12492). Control diet (chow) containing 24% protein, 47.5% carbohydrate, and 4.9% fat, was given to age and sex matched animals as a control group. Cells were isolated from tissues obtained from sacrificed animals (see Supplementary Methods). Single live cells were purified via fluorescent activated cell sorting, excluding dead cells labeled with LIVE/DEAD dyes (Thermo) and doublets based on side and forward light scatter. Lamina propria cells were further sorted based on CD45 expression (BioLegend, clone: 30-F11).

Monocyte/Macrophage isolation and transfer

Monocytes were isolated from bone marrow using with MojoSort^™ Mouse Monocyte Isolation Kit (Biolegend) following the manufacturer’s instructions or via fluorescence activated cells sorting using fluorochrome-conjugated monoclonal antibodies: CD11b (Biolegend, clone: M1/70), F4/80 (Biozol, clone: BM8), SiglecF (BD Horizon, clone: E50–2440), Ly6G (BioLegend, clone: 1A8), CD11c (BioLegend, clone: N418), MHC-II (BioLegend, clone: M5/114.15.2), Ly6C (BioLegend, clone: HK1.4). Macrophages were isolated from peritoneal lavage via magnetic sorting using biotinylated monoclonal antibody against TIM4 (Miltenyi Biotec, clone: REA999), and anti-biotin Microbeads (Miltenyi Biotec, Cat# 130–090-485) following the manufacturer’s instructions. Isolated cells were stained with Cell Trace Violet (Life Technologies), assessed for purity via flow cytometry and then 0.5–1×10⁶ cells/mouse were transferred via intraperitoneal injection to littermates randomly assigned to experimental groups. RFP⁺ monocytes were isolated from B6.RFP mice via Immuno-magnetic depletion of whole bone marrow by incubating samples with biotinylated antibodies against CD3 (eBioscience, clone: 145–2C11), CD4 (eBioscience, clone: GK1.5), CD8 (eBioscience, clone: 53–6.7), CD45R (Biolegend, clone: RA3–6B2), CD19 (eBioscience, clone: eBio1D3), NK1.1 (eBioscience, clone: PK136), Ter119 (Biolegend), CD49b (Biolegend, clone: DX5), Ly6G (Biolegend, clone: 1A8) and CD117 (Biolgend, clone: 2B8). Anti-Biotin microbeads (Miltenyi) were then used according to the manufacturer’s protocol. Enrichment was validated by staining purified bone marrow monocytes with antibodies against CD115 (Biolegend, clone: AFS98), Ly6C (BD Bioscience, clone: AL-21).

Skin wounding model

Wounding, preparation of wound tissue and cell isolation was performed as previously described (67). 3×10⁶ RFP⁺ monocytes were adoptively transferred into each previously wounded C57/BL6 recipient via intravenous injection, either 2 or 12 days after injury. Wounds were excised 4 or 14 days after injury. Recovered cells were stained for flow cytometry and cell sorting using CD11b (eBioscience, clone: M1/70), F4/80 (eBioscience, clone: BM8 or AbD Serotec, clone: CI:A3–1), MHC-II (eBioscience, clone: M5/114.15.2), Ly6C (BD, clone: AL-21), CD45 (eBioscience, clone: 30-F11) and CD301b (Biolegend, clone: URA-1). Dead cells were excluded by DAPI staining (50 ng/mL - ThermoFisher scientific).

Flow cytometry

Cells were stained for flow cytometric analysis using TIM4 (BioLegend, clone: F31–5G3), CD45.1 (BioLegend, clone: A20), CD115 (BioLegend, clone: AFS98) and F4/80 (BioLegend, clone: BM8). Detection of intracellular Retnla mRNA and RELMα protein (Peprotech, Cat# 500-P214) was performed using PrimeFlow RNA Assays (Thermo scientific) following the manufacturer’s instructions. Data from stained cells were collected using LSR Fortessa flow cytometers (BDBiosciences) with FACSDiva v. 9.0 and data were processed using FlowJo v. 10.6. For spectral flow cytometry, cells were stained with specified antibodies (see Supplementary Methods) and collected using Cytek^® Aurora (Cytek) spectral flow cytometer at the Sidney Kimmel Comprehensive Cancer Center High Parameter Flow Core. Data were processed using CytoExploreR v.1.1.0 (68), followed by UMAP calculation and fluorescence data scaling with a linear regression model in R v. 4.0.0 (69) using Seurat v.4. Data clustering was performed using the Leiden algorithm (70). Stage assignment was performed manually.

Single cell sequencing and analysis:

Detailed strategy is included in the supplementary methods. Briefly, scRNA-seq of cells from adipose tissue and lamina propria was performed using a 10X Genomics Chromium Controller. Single cells were processed with GemCode Single Cell Platform using GemCode Gel Beads, Chip and Library Kits (v2) following the manufacturer’s protocol targeting 5000 cells sequenced on HiSeq 3000 (Illumina) with 50000 reads per cell. scRNA-seq of wounded skin was performed as described previously (71, 72) with Illumina 50 bp paired-end sequencing on a Novaseq6000 aiming for 0.5 mio reads per cell. Samples were demultiplexed and aligned using Cell Ranger (10X genomics), gsnap v.2018–07-04 (73) and featureCounts v. 1.6.2 (74) to genome build GRCm38 release 87. Read count matrices were processed, analyzed and visualized using Seurat v.4 (38) in R v. 4.0.0 (69) ad scanpy v.1.8 (32) in Python v.3 (75), complemented by SCTransform (30), Clustree (76), Harmony (31), BBKNN (33), PAGA (34), slingshot v. 1.6.1 (28) and Monocle v.3 (44) plus visualized with a Uniform Manifold Approximation and Projection (UMAP) (77) as a dimensionality reduction approach. Pathway enrichment and protein interaction network was determined using goseq v. 1.40.0 (78), StringDB (79), igraph v. 1.2.5 (80) and ggraph v. 2.0.2 (81). Network clustering was performed using the cluster_louvain function (82). Trajectory dependent gene regulation was calculated on top 2000 most variable genes using gam v. 1.20 (83). Query datasets available from public repositories were retrieved (table S2), macrophages extracted based on a “Macrophage score” (table S1) and processed as described above. Flow cytometry data associated with individually barcoded cells was used to calculate a UMAP projection using uwot (84), which was then clustered with kmeans v. 4.0.0 (69), setting k based on the total within cluster sum of squares. Transcription factor enrichment analysis was performed with RcisTarget v. 1.8.0 (85, 86) and visualized using wordcloud2 v 0.2.1 (87).

Quantification and statistical analysis

Statistical analysis was performed in R v. 4.0.0 (69), using functions from the base stats package to calculate a single factor anova with aov followed by Tukey Honest Significant Differences to determine statistically significant differences between means (fig. 2G–H) with a Bonferroni correction or only a single factor anova to establish significant clustering effects (fig. 5M).Tests were two-sided and alpha = 0.05. Ex vivo results are represented as dots for individual mice.