Research Article

Specification of tissue-resident macrophages during organogenesis

See allHide authors and affiliations

Science  09 Sep 2016:
Vol. 353, Issue 6304, aaf4238
DOI: 10.1126/science.aaf4238

Structured Abstract

INTRODUCTION

Embryonic development and tissue homeostasis depend on cooperation between specialized cell types. Resident macrophages are professional phagocytes that survey their surroundings; eliminate unfit cells, microorganisms, and metabolic waste; and produce a large range of bioactive molecules and growth factors. Resident macrophages also serve tissue-specific purposes: For example, microglia in the central nervous system support neuronal circuit development, Kupffer cells scavenge blood particles and dying red blood cells in the liver, and alveolar macrophages uptake surfactant and remove airborne pollutants and microbes from the airways. Resident macrophage diversity in adult mice is reflected in tissue-specific gene expression profiles, which may be due to responses to specific cues from their microenvironment, different developmental processes, and the contribution of distinct progenitors cell types. Altogether, the mechanisms responsible for the generation of tissue-resident macrophage diversity remain unclear.

RATIONALE

Tissue-resident macrophages originate, at least in part, from mesodermal erythro-myeloid progenitors (EMPs) from the yolk sac, which invade the embryo proper at the onset of organogenesis. These tissue-resident macrophages are also self-maintained in postnatal tissues, independently of definitive hematopoietic stem cells (HSCs) in a steady state. We therefore hypothesized that resident macrophages represent a founding cell type within most organ anlagen. In this model, the generation of macrophage diversity, as observed in the tissues of postnatal mice, may be integral to organogenesis.

RESULTS

To test this hypothesis and explore the molecular basis of macrophage diversity in mammals, we performed a spatiotemporal analysis of macrophage development in mice, from embryonic day 9 (E9) to 3 weeks after birth. Unbiased single-cell RNA sequencing (RNA-seq) analysis of CD45+ cells, combined with RNA-seq analyses of sorted cell populations, genetic fate mapping, and in situ analyses, revealed that EMPs give rise to a population of premacrophages (pMacs) that colonize the whole embryo from E9.5, as they acquire a core macrophage differentiation program that includes pattern recognition, scavengers, and cytokine receptors. The chemokine receptor Cx3cr1 is up-regulated in pMacs and is important for embryo colonization, which is delayed in Cx3cr1-deficient embryos. Fate mapping of pMacs using a Tnfrsf11a–Cre reporter labels homogeneously fetal and adult tissue-resident macrophages but not HSCs and their progeny. Transcriptional regulators that identify postnatal tissue-resident macrophages in the brain, liver, kidney, skin, and lung were specifically up-regulated immediately after colonization. These dynamic changes mark the onset of diversification into adult macrophages. We identified Id3 as a Kupffer cell–specific transcriptional regulator. Deletion of Id3 in pMacs resulted in Kupffer cell deficiency but did not affect development of microglia and kidney macrophages.

CONCLUSION

Our study shows that EMP-derived precursors colonize embryonic tissues and simultaneously acquire a full core macrophage program. This is followed by their diversification into tissue-specific macrophages during organogenesis, likely via the expression of distinct sets of transcriptional regulators. These results indicate that differentiation of tissue-resident macrophages is an integral part of organogenesis and identify a spatiotemporal molecular road map for the generation of macrophage diversity in vivo. Our findings provide a conceptual framework to analyze and understand the consequence(s) of genetic variation for macrophage contribution to development, homeostasis, and disease pathogenesis in different tissues and will support efforts to differentiate specialized macrophages in vitro.

Specification of tissue-resident macrophages.

Erythro-myeloid progenitors (EMPs) from the yolk sac colonize the fetal liver and give rise to macrophage precursors (pMacs) that acquire a core macrophage transcriptional program and colonize the embryo from E9.5 in a Cx3cr1-dependent manner (green arrows). Specification of F4/80+ resident macrophages (brown arrows), starting from E10.25, is initiated by the expression of tissue-specific transcriptional regulators. Id3 (red) is important for Kupffer cell development. Transcription factors noted in blue have been shown to be important for the differentiation or the maintenance of the corresponding macrophage subsets. MΦ, macrophage.

Abstract

Tissue-resident macrophages support embryonic development and tissue homeostasis and repair. The mechanisms that control their differentiation remain unclear. We report here that erythro-myeloid progenitors in mice generate premacrophages (pMacs) that simultaneously colonize the whole embryo from embryonic day 9.5 in a chemokine-receptor–dependent manner. The core macrophage program initiated in pMacs is rapidly diversified as expression of transcriptional regulators becomes tissue-specific in early macrophages. This process appears essential for macrophage specification and maintenance, as inactivation of Id3 impairs the development of liver macrophages and results in selective Kupffer cell deficiency in adults. We propose that macrophage differentiation is an integral part of organogenesis, as colonization of organ anlagen by pMacs is followed by their specification into tissue macrophages, hereby generating the macrophage diversity observed in postnatal tissues.

Tissue-resident macrophages are a diverse family of cells found in most organs. They include brain microglia, liver Kupffer cells, lung alveolar macrophages, and epidermal Langerhans cells. In mice, tissue-resident macrophages share an embryonic origin and differentiate, at least in part, from yolk sac (YS) erythro-myeloid progenitors (EMPs) (1, 2). Tissue-resident macrophages are also self-maintained in adult tissues, independently of hematopoietic stem cells (HSCs) under steady-state conditions (36). However, the mechanisms responsible for the generation of macrophage diversity observed in adult mice remain unclear. It has been proposed that resident macrophage diversity reflects their exposure to specialized tissue environments (710) or the contribution of distinct embryonic or fetal progenitors to distinct subsets (2, 1113). The preferential expression of transcription factors in macrophage subsets was also noted (7) and appears functionally important. Several such cases—including Gata6 for large peritoneal macrophages (9, 10, 14), Runx3 for Langerhans cells (15), Nr1h3 for splenic marginal zone macrophages (16), SpiC for splenic red pulp macrophages (17), and Pparg for alveolar macrophages (18)—have been functionally validated by knockout mice. To better understand how macrophage diversity is generated, we performed a molecular and spatiotemporal analysis of macrophage development in mice.

Colonization of developing tissues by EMP-derived macrophage precursors

Erythro-myeloid progenitors (Csf1r+ Kit+ CD45low AA4.1+) are first detected in the YS at embryonic day 8.5 (E8.5) (2, 19) and subsequently colonize the fetal liver (1, 2). Previous fate-mapping analysis of EMP differentiation indicated that their progeny lose Kit expression and increase CD45 expression as they invade the embryo before acquiring F4/80 expression to give rise to fetal and postnatal tissue-resident macrophages (2) (Fig. 1A). To explore the spatiotemporal and molecular determinants of macrophage differentiation and diversification, we first performed whole-transcriptome sequencing of sorted CD45low Kit+ EMPs, CD45+ Kit Lin (Ter119, Gr1, F4/80) cells and F4/80+ macrophages from embryonic and postnatal tissues up to 3 weeks after birth (Fig. 1A and fig. S1, A and B). We identified genes that were significantly up-regulated between EMPs and CD45+ Kit Lin cells [adjusted P ≤ 0.05, DESeq2 (20), Benjamini Hochberg (BH) correction, table S1]. Subsequent summarization and visualization of these genes via scorecard analysis (21) (Fig. 1B) indicated that the signature of CD45+ Kit Lin cells was also present in F4/80+ macrophages across tissues and over the entire time-course analysis in the embryo and postnatal mice in the kidney, liver, and brain. Epidermal Langerhans cell and lung alveolar macrophage signatures were modified after birth (Fig. 1B). A second scorecard analysis of genes up-regulated between EMPs and early (E10.25 to E10.5) F4/80+ macrophages identified a signature that was already detectable in CD45+ Kit Lin cells and conserved in later macrophages across tissues (fig. S1C and table S1). Unsupervised principal component analysis (PCA) showed a distinct grouping of EMPs, CD45+ Kit Lin cells, and macrophages, regardless of their tissue of origin—i.e., YS, liver, head, or caudal region (fig. S1D). Morphologically, CD45+ Kit Lin cells from the YS, fetal liver, head, and caudal embryo displayed a similar morphology, as did macrophages from the same tissues (Fig. 1C). CD45+ Kit Lin cells resembled EMPs, albeit with the presence of occasional phagocytic vacuoles, whereas phagocytic features become prominent in F4/80+ cells (Fig. 1C).

Fig. 1 A core macrophage program is initiated simultaneously in pMacs in all tissues.

(A) Summary of surface phenotype used for EMPs, pMacs, and macrophages (MΦ). (B) Scorecard visualization of differentially up-regulated genes (DESeq2 Wald test, adjusted P < 0.05, BH correction) in pMacs (E9.5 and E10.25) compared with EMPs, on the basis of whole-transcriptome sequencing of indicated populations. The table shows the relative enrichment of differentially up-regulated genes in pMacs across cell types and tissues (y axis) and developmental time points (x axis, from E9 to P21). Each sample represents the mean of at least two biological replicates and two technical replicates, except for E14.5 liver macrophages and P8 and P21 lung macrophages, which consist of one biological replicate and two technical replicates. See table S1, fig. S1, and methods for details of the scorecard. YS, yolk sac; FL, fetal liver; epi., epidermis. (C) May-Grünwald-Giemsa–stained cytospin preparations of sorted EMPs, pMacs, and early macrophages from YS, head, limbs, and fetal liver at E10.25 and E12.5. Images are representative of n = 3 independent experiments. Scale bars, 10 μm. (D) t-Distributed stochastic neighbor embedding (tSNE) plot of single-cell RNA-seq data showing distribution of CD45low/+ cells from E10.25 embryos into three clusters (see also fig. S2). Cluster distribution based on DBScan is overlaid onto the graph. In total, we analyzed 408 single cells from n = 2 independent experiments with four to six embryos per experiment. (E) Superimposition of EMP-, pMac-, or macrophage-specific signatures defined by the bulk RNA-seq on the tSNE plot shown in (D). (F) tSNE plot as in (D), overlaid with the relative expression values for Kit and Maf.

These data suggest that a macrophage differentiation program was initiated simultaneously in the whole embryo in CD45+ Kit Lin cells; these macrophage precursors will subsequently be referred to as premacrophages (pMacs). To further test this hypothesis, we performed independent, unbiased whole-transcriptome single-cell RNA sequencing (scRNA-seq) of CD45low/+ cells purified from the whole embryo at E10.25 [30 to 34 somite pairs (sp)]. Nonlinear dimensionality reduction in combination with unsupervised clustering of cells indicated that these cells are best described by three major clusters (Fig. 1D and fig. S2). Overlay of EMP, pMac, and macrophage signatures from differentially expressed genes in bulk RNA-seq analysis (table S2) indicated superimposition on clusters 1, 2, and 3, respectively (Fig. 1E). Intermediate differentiation states in EMPs, pMacs, and macrophages were clearly apparent, which suggests a gradual differentiation path from EMPs to macrophages via pMacs (Fig. 1E), consistent with the scorecard analysis.

Core macrophage transcriptional program

Analysis of genes differentially expressed in EMPs, pMacs, and early macrophages by scRNA-seq (Figs. 1F and 2A and table S2) and bulk RNA-seq (Fig. 2B and table S1) confirmed that a core macrophage transcriptional program was initiated in pMacs. As Kit, Gata1, and Gata2 expression was lost, pMacs up-regulated expression of Csf1r; the transcription factors Maf, Batf3, Pparg, Irf8, and Zeb2; the chemokine receptor Cx3cr1; cytokine receptors; complement and complement receptors; pattern-recognition receptors; phagocytic receptors; Fc gamma receptors; the inhibitory receptor Sirpa; MerTK; cathepsins; Aif1 (Iba1); Emr1 (F4/80); and Grn (Granulin). Expression of selected genes was confirmed at the protein level by fate mapping in EMP-derived cells. Fate-mapping of EMPs was performed by pulse labeling of Csf1rMeriCreMer; Rosa26LSL-YFP E8.5 embryos, with 4-hydroxy-tamoxifen (OH-TAM) (2, 5). Short-lived OH-TAM allows transient nuclear translocation of the estrogen receptor–Cre recombinase fusion protein (MeriCreMer) in cells expressing the Csfr1MeriCreMer transgene and deletion of a floxed stop cassette (LSL) in the Rosa26LSL-YFP allele, resulting in stable yellow fluorescent protein (YFP) expression in targeted cells and their progeny (2, 5). Expression of cytokine receptors for interleukin (IL)–4, IL-13, interferon-γ, and tumor necrosis factor–α; phagocytic and activating receptors Mrc1 (CD206), Trem2, and Dectin-1 (Clec7a); Fc-gamma receptors (Fcgr1, Fcgr2/3, Fcgr4); and Iba1 and Grn was confirmed by flow cytometry and immunofluorescence in situ on EMP-derived pMacs and early macrophages from the YS, as well as from the head, caudal region, limbs, and liver of the embryo proper; fetal macrophages; and adult tissue macrophages (Figs. 3, A and B; 4A; and 5A and figs. S3 and S4). Protein expression was first detected in pMacs and increased as pMacs differentiated into F4/80+ macrophages (Fig. 3, A and B, and fig. S3). Of note, pMacs and macrophages represented 70 to 90% of EMP-derived cells in the head and caudal embryo, whereas 80 to 90% of YFP+ cells in the E10.25 fetal liver represented progenitors (Fig. 3B and fig. S4).

Fig. 2 Differentially expressed genes during differentiation from EMP to macrophage.

(A) Kinetic diagrams show the pseudotemporal behavioral of 408 single CD45low/+ cells from scRNA-seq at E10.25, as aligned by Monocle 2, a tool used to analyze differentially expressed genes during an ongoing biological process such as differentiation. See fig. S2G for progression of single cells over developmental pseudotime. Cells within cluster 1 (Fig. 1D) are found farther left on the pseudotime scale, cluster 2 cells are in the middle, and cells within cluster 3 are farther to the right. (Upper panel) Developmental pseudotime diagram (q < 0.05) showing down-regulation of EMP-specific genes [differentially expressed compared with the macrophage and pMac cluster, P < 0.05, fold change (FC) > 1.4] over the differentiation path from EMP to pMacs and macrophages. (Lower panel) Similar plot depicting macrophage-specific genes significantly regulated over pseudotime (q < 0.05) and differentially expressed compared with the EMP and pMac cluster (P < 0.05, FC > 1.4). (B) Heat-map representation of selected genes differentially regulated between EMPs versus pMacs and EMPs versus early macrophages in bulk RNA-seq analysis. Black boxes were drawn around those samples used for differential expression analysis. See also tables S1 and S2.

Fig. 3 EMP-derived pMacs colonize the embryo to generate macrophages.

(A) Flow cytometry analysis of E10.25 Csf1rMeriCreMer; Rosa26LSL-YFP (OH-TAM at E8.5) tissues showing expression of Il4ra, Il13ra1, CD16.2, CD64, Ifngr, Tnfr2, Tim4, and CD206 on YFP+ Kit+ progenitors, pMacs, and macrophages. MFI, mean fluorescence intensity. Data are representative of n = 4 independent experiments with four to six embryos per marker. See also fig. S3A. (B) Quantification of immunostainings on cryosections from E10.25 Csf1rMeriCreMer; Rosa26LSL-YFP embryos, pulse-labeled with OH-TAM at E8.5 with antibodies against YFP, Iba1, and additional candidate genes (CD16/32, Dectin-1, Trem2, F4/80, CD206 or Granulin). Orange bars represent the percentage of YFP+ cells expressing Iba1 and the indicated candidate; blue bars denote the percentage of YFP+ cells expressing only Iba1. Red and green bars represent the percentage of YFP+ cells negative for Iba1 and expressing (red) or negative for (green) the indicated candidate gene. n = 2 to 4 embryos and two sections per embryo per marker. See also fig. S4. (C) tSNE plot (as in Fig. 1D) overlaid with the relative expression values for Tnfrsf11a and Cx3cr1. (D) Genetic labeling efficiency in Tnfrsf11aCre+; Rosa26LSL-YFP mice measures efficiency of Cre-mediated recombination of the floxed stop sequence at the Rosa26LSL-YFP locus and, therefore, transcription of Tnfrsf11a in the indicated cell types and their progenitors. Cells analyzed were pMacs and F4/80+ macrophages in the YS and whole embryo at E10.25; fetal liver HSCs [long-term (LT) (LinKit+Sca1+CD150+CD48), short-term (ST) (LinKit+Sca1+CD150CD48) and multipotent progenitor (MPP) (LinKit+Sca1+CD150CD48+)] and tissue macrophages at E14.5 and 6 weeks; and blood leukocytes [B cells (CD19+), T cells (CD19Ly6GCD115CD3+), natural killer (NK) cells (CD19Ly6GCD115CD3NKp46+], neutrophils (CD19Ly6G+) and Ly6Chi monocytes (CD19CD115+Ly6GLy6Chi), and tissue CD11bhigh myeloid cells from 6-week-old mice. Circles represent individual mice from n = 4 independent experiments. See also fig. S5.

Fig. 4 Tissue colonization by pMacs is dependent on Cx3cr1.

(A) Expression of GFP and Dectin-1 in Cx3cr1gfp/+ mice during development (E8.5 to 10.5) in Kit+ cells (CD45low, Kit+), pMacs, and macrophages. sp, somite pairs. Data are representative of n = 9 independent experiments. Biological replicates have been aggregated per cell type, time point, and tissue. (B) Flow cytometry analysis in Cx3cr1+/− and Cx3cr1−/− of pMacs and macrophages from the YS, head, and caudal region at E9.5 (upper panel) and the liver, YS, head, and limbs at E10.5 (lower panel). Circles represent individual mice. Data are pooled from n = 6 independent experiments. P were calculated using Student’s t test.

Fig. 5 Early specification of tissue-resident macrophages.

(A) Flow cytometry analysis of Csf1rMeriCreMer; Rosa26LSL-YFP (OH-TAM at E8.5) liver, brain, lung, and skin F4/80+ cells from postnatal mice (4 weeks old) showing expression of Il4ra, Il13ra1, Tnfr2, Ifngr, Dectin-1, CD64, Tim4, and CD206. Black dashed lines depict expression by all macrophages (CD45+CD11blowF4/80high); green tinted lines show expression by the subsets of YFP+ macrophages. Gray histograms show the fluorescence intensity of the FMO (fluorescence minus one) controls for all macrophages (gray lines) and YFP+ macrophages (gray tinted lines). Histograms are representative of five mice and n = 2 independent experiments. (B and C) Scorecard analysis of all differentially up-regulated genes in postnatal macrophages from bulk RNA-seq. The scorecards show the relative enrichment of each set of up-regulated genes across each cell type (y axis) and developmental time point (x axis). See Materials and methods for details of the score card. Numbers for each population indicate differentially up-regulated transcripts in postnatal (P2 to P21) brain, liver, kidney, epidermis, or lung macrophages when comparing one population versus the others. See also table S3. (D) Heat map representing all differentially up-regulated transcriptional regulators from bulk RNA-seq (twofold change, adj. P < 0.05, BH correction) between postnatal macrophages from the brain, liver, kidney, skin, and lung, as well as their relative expression in tissue macrophages from E10.25 to P21. LC, Langerhans cells; AM alveolar macrophages.

Detection of the early expression of the cytokine receptor Tnfrsf11a in pMacs by RNA-seq and scRNA-seq analyses (Figs. 2B and 3C) predicted that pMacs and their progeny can be genetically labeled in Tnfrsf11aCre (22); Rosa26LSL-YFP mice. The Tnfrsf11aCre transgene allows expression of the Cre recombinase in cells expressing Tnfrsf11a and deletion of the floxed stop cassette from the Rosa26LSL-YFP allele. We observed YFP labeling by flow cytometry in ~80% of pMacs and early macrophages from the YS and embryo of Tnfrsf11aCre; Rosa26LSL-YFP mice (Fig. 3D and fig. S5). Fetal macrophages in all tissues also expressed YFP at comparable levels at E10.25 and E14.5 (Fig. 3D). In addition, ~80% of brain, lung, epidermis, kidney, and liver macrophages from 6-week-old mice expressed YFP (Fig. 3D). Therefore, the whole resident macrophage lineage is labeled in Tnfrsf11aCre; Rosa26LSL-YFP mice, although Tnfrsf11a expression itself is lost in postnatal Langerhans cells and alveolar macrophages (Fig. 2B). Moreover, we noted that YFP expression was observed in only ~15% of fetal HSCs, adult HSCs, and HSC-derived cells in the blood and tissues of adult mice (Fig. 3D and fig. S5). Tnfrsf11aCre; Rosa26LSL-YFP mice thus represent an efficient and relatively specific model for genetic labeling of tissue-resident macrophages in fetuses and adult mice.

Early Cx3cr1 expression by pMacs (Figs. 2B and 3C) is consistent with previous reports showing Cx3cr1 expression in macrophage precursors (23) and suggests that this chemokine receptor may be involved in colonization of embryo tissues by pMacs. Time-course analysis of green fluorescent protein (GFP) expression in Cx3cr1gfp/+ mice (where GFP expression denotes transcription of Cx3cr1) from E8.5 (19 to 21 sp) to E10.5 (38 to 39 sp) indicated that GFP expression, and therefore Cx3cr1 expression, is not detected in Kit+ progenitors but is up-regulated in Kitlow Dectin-1+ pMac and is highest in F4/80+ macrophages that appear at E10.25 (Fig. 4A). Next, we found that colonization of the head and caudal tissues is delayed in Cx3cr1-deficient embryos, as pMac and macrophage numbers are decreased in the head, caudal region, and limbs of E9.5 and E10.5 embryos in comparison with Cx3cr1+/− littermates while they accumulate in the YS and fetal liver (Fig. 4B). Nevertheless, tissue macrophage numbers even out in the consecutive days of embryonic development in most tissues, with the exception of the kidney, where a 50% lower number in resident macrophages is still observed in adult mice, in line with previous research (24) (fig. S6).

A progressive enrichment of pMac- and macrophage-specific gene expression signatures was observed in gene set enrichment analysis (GSEA) of bulk RNA-seq data (fig. S7, A and B). Moreover, the core macrophage signature identified in adult mice by the Immunological Genome Consortium (25) was already enriched in the genes up-regulated in pMacs and early macrophages compared with EMPs (fig. S7, C and D). In silico investigation of transcription factor binding sites identified using chromatin immunoprecipitation sequencing (ChIP-seq) in the proximity of up-regulated genes [transcription start site (TSS) ± 20 kb] in pMacs and macrophages further supports the proposition that pMacs undergo a coordinated macrophage differentiation program: A Locus Overlap Analysis (LOLA) (26) yielded a statistically significant association (adjusted P ≤ 0.001, Benjamini Yekutieli correction) with binding sites for Spi1, Egr1, Irf1, Irf8, Maf, Jun, Stat1, Stat3, Stat5b, Stat6, Rela, and Relb (fig. S8A). These factors are expressed in our data set (fig. S8B and tables S1 and S5), and their binding sites align at the same loci in enhancers and superenhancers associated with differentially regulated genes such as thrombospondin1 (Thbs1) (27), Cx3cr1 (8), F4/80 (Emr1), and Mrc1 but not with control genes such as Gata1 and MyoD (fig. S8C). Comparable results, with a higher statistical significance, were found for genes differentially up-regulated in early macrophages when compared with EMPs (fig. S8A).

Together, these results characterize, at the cellular and molecular level, the EMP-derived macrophage precursors (pMacs) that acquire a core macrophage transcriptional program as they colonize the head and caudal embryo from E9.5 in a Cx3cr1-dependent manner to give rise to tissue-resident macrophages. These data are consistent with our previous demonstration that the vast majority of resident macrophages in these tissues originate from a progenitor that does not express the stem cell and endothelial marker Tie2 after E10.5 (2). However, our results do not exclude the possibility that later precursors may also contribute to the resident macrophage pool.

Early specification of tissue-resident macrophages

Following the acquisition of this core program, F4/80+ macrophages soon display heterogeneity between different tissues (Figs. 2B and 5A). For example, expression of Timd4 is lost at the transcriptional and protein level in microglia, alveolar macrophages, and Langerhans cells but is maintained in Kupffer cells throughout development and postnatally. Expression of the mannose receptor (CD206) is lost in microglia and Langerhans cells. Dectin-1 and CD64 expression is maintained in all subsets except Langerhans cells (Figs. 2B and 5A and fig. S3). To systematically investigate the kinetics and molecular determinants of macrophage diversification, we first characterized tissue-specific signatures of genes differentially up-regulated in postnatal microglia, kidney macrophages, Langerhans cells, alveolar macrophages, and Kupffer cells in the bulk RNA-seq data set. Unsupervised clustering analysis suggested that macrophages in different tissues undergo characteristic differentiation trajectories (fig. S9A). Supervised analyses identified lists of genes differentially up-regulated in each cell population (table S3 and fig. S9, B and C). The signatures of adult tissue-resident macrophages of the brain, lung, and liver previously defined by independent research (7, 25, 28) were progressively enriched in our developing resident macrophage populations (fig. S10, A to D). Finally, scorecard (21) and GSEA analysis of differentially up-regulated genes for each postnatal macrophage population in all tissues and over time from E9 to P21 (Fig. 5, B and C, and fig. S10C) revealed that the tissue-specific signature of postnatal microglia, Kupffer cells, and kidney macrophages could be traced back to fetal macrophages as early as E12.5 (Fig. 5B). The signatures of Langerhans cells and alveolar macrophages reflected important postnatal changes in gene expression (Fig. 5C), noted in previous studies (29), which may reflect their anatomical location at epithelial barriers.

Tissue-specific sets of transcriptional regulators define macrophage diversity

A heat-map visualization of all transcriptional regulators present in the postnatal tissue-specific signatures (twofold change, adj. P < 0.05, BH correction, Fig. 5D) confirmed the tissue-specific expression of the transcription factors Sall1 and Sall3 in microglia (7), Nr1h3 (Lxra) in Kupffer cells (30), Pparg in lung alveolar macrophages (18, 31), and Runx3 and Ahr in Langerhans cells (15, 32) (Fig. 5D) and identified additional tissue-specific transcriptional regulators such as Id1 and Id3 in Kupffer cells (Fig. 5D). In addition, this analysis showed that many of these transcriptional regulators start to be differentially expressed in tissue macrophages as early as E10.25—for example, in Sall1 and Sall3 in head macrophages; Nr1h3, Id1, and Id3 in the liver; and Ahr in limb macrophages (Fig. 5D). Some genes, like Id1 and Sall3, are expressed by progenitors and pMacs before their expression becomes restricted to macrophages in the liver and the head, respectively (fig. S10E), whereas expression of other genes, such as Id3 and Sall1, is low in progenitors and up-regulated in pMacs (fig. S10E). In situ immunofluorescence confirmed expression of Id3 by E10.25 macrophages in the liver and head and of Id1 in liver macrophages (Fig. 6A). Together these data suggest that the transcriptional programs of tissue-specific resident macrophages start to be established early on, as soon as macrophages or pMacs are present in tissues, and identify a number of “candidate” tissue-specific transcriptional regulators.

Fig. 6 Tissue-specific macrophage signatures are not detected in pMacs.

(A) Immunostaining with antibodies against Id3 (red, upper panel) or Id1 (red, lower panel), F4/80 (cyan), and YFP (green) on cryosections from E10.25 Csf1rMeriCreMer; Rosa26LSL-YFP embryos (OH-TAM at E8.5), analyzed by confocal microscopy. Nuclei are counterstained with 4′,6-diamidino-2-phenylindole (DAPI) (white). Scale bars, 2 μm. Images are representative of n = 3 independent experiments. (B) tSNE plots of scRNA-seq data from CD45low/+ cells from E10.25 embryos (as in Fig. 1), showing coexpression of Id1, Id3, and Sall3. See also fig. S11. (C) PCA plots of scRNA-seq data of cells from cluster 2 (pMacs, defined in Fig. 1D) with superimposed fetal tissue macrophage-specific signatures. See also fig. S11, table S2, and Materials and methods.

We thus investigated whether specification of tissue macrophages takes place in F4/80+ macrophages or at the level of their pMac or EMP precursors. We plotted transcriptional coexpression of Id1, Id3, and Sall3 onto the t-distributed stochastic neighbor embedding (tSNE) representation of single CD45low/+ cells from our scRNA-seq data set (Figs. 1D and 6B and fig. S11). Coexpression of Id1 and Id3 was found in pMacs and macrophages. However, Id1 and Sall3 were coexpressed in EMPs and pMacs, and Id3 and Sall3 were coexpressed in pMacs, suggesting that cells at the pMac state are not completely committed to exclusive expression of tissue-specific transcription factors. These data confirmed that Id1 and Sall3 are expressed by progenitors and pMacs (fig. S10E), although their expression is ultimately lost by tissue macrophages outside the liver and the brain, respectively (Figs. 5D and 6B and fig. S10E). When common and tissue-specific E14.5 to E18.5 macrophage signatures (table S2) were superimposed on the pMac population, we found that the common macrophage signature was expressed in pMacs with a gradual enrichment within this population (Fig. 6C). However, we did not observe pMac subsets expressing tissue-specific signatures (Fig. 6C). The lack of tissue specificity within pMacs was confirmed using an alternative bioinformatic approach based on analysis of multimodal expression followed by hierarchical clustering of genes with subsequent analysis of enrichment of tissue signatures within these clusters in pMacs (fig. S11E).

Taken together, these data suggest that expression of Id1 and Id3 or Sall3 does not specify pMac subsets precommitted to give rise to either Kupffer cells or microglia. More generally, diversification of tissue macrophages appears to take place after pMacs have colonized tissues and differentiated into F4/80+ macrophages.

The transcriptional regulator Id3 is essential for Kupffer cell development

To functionally validate the role of early tissue-specific transcriptional regulators on the development and specification of tissue macrophages—and, thereby, to validate the importance of our findings—we studied the role of Id3 expression in Kupffer cell differentiation. E10.25 Id3-deficient embryos had normal or increased numbers of pMacs and early macrophages in the YS, but macrophages were reduced in numbers in the embryo proper (liver and head) in comparison with littermate controls (Fig. 7A). The further development of liver macrophages was severely impaired in E14.5 and E18.5 Id3-deficient embryos, as determined by flow cytometry and histology, and 4-week-old Id3-deficient mice still presented with a marked Kupffer cell deficiency, whereas development of microglia and kidney macrophages appeared normal (Fig. 7, B, C, and E, and fig. S12C). The role of Id3 in Kupffer cells appears to be cell-autonomous, as targeted deletion of an Id3 floxed allele in pMacs (Tnfrsf11aCre+; Id3f/f) recapitulated the phenotype of the Id3-deficient mice in embryo and postnatal mice (Fig. 7, C and D). Expression of Id3 in postnatal Kupffer cells was confirmed by quantitative reverse transcription polymerase chain reaction (qRT-PCR) (fig. S12A) and by immunofluorescence in Kupffer cells from Csf1rMeriCreMer; Rosa26LSL-YFP mice pulse-labeled with OH-TAM at E8.5 (Fig. 7F and fig. S12B). In contrast, fate-mapped microglia do not express Id3 (fig. S12B), in line with our RNA-seq data (Fig. 5D). Of note, the partial Kupffer cell deficiency observed in Id3-deficient and Tnfrsf11aCre+; Id3f/f embryos and adults was not associated with an abnormal liver lobular architecture or vasculature (Fig. 7, C and D, and fig. S12D). Kupffer cell proliferation in the steady state was not affected by Id3 deficiency (Fig. 7G), but RNA-seq analysis of Id3−/− and Id3+/− adult Kupffer cells indicated the up-regulation (>threefold) of Id1 expression, and gene ontology (GO)–term analysis evidenced that Id3−/− cells overexpressed genes involved in the control of cell death and cytokine responses and down-regulated genes involved in metabolic processes (Fig. 7H; fig. S12, E and F; and table S6). These data suggest that Id3 is important for the development and maintenance of Kupffer cells in the liver but is dispensable for other macrophage subsets. As Id1 and Id3 are coexpressed, it will be interesting to investigate whether up-regulation of Id1 partially compensates for Id3 deficiency (33, 34).

Fig. 7 Id3 is important for Kupffer cell development.

(A) Flow cytometry analysis of the percentage of pMacs and macrophages in the YS and liver from E10.25 Id3−/− and Id3+/− embryos. Circles represent individual mice from n = 3 independent experiments. (B) Flow cytometry analysis of F4/80+ macrophages in the liver, brain, and kidney from E14.5 and E18.5 Id3−/− and Id3+/− mice. Cell suspensions were prepared, stained, and acquired from whole organs, and the number of live cells per tissue was directly obtained from flow cytometry standard (FCS) files (see Materials and methods). Circles represent individual mice from n = 4 independent experiments. (C and D) Immunostaining with antibodies against CD31, which labels endothelial cells, and F4/80 on liver cryosections from E14.5 and 2-week-old (2 w.o) Id3−/−, Id3+/−, Tnfrsf11aCre-;Id3f/+, Tnfrsf11aCre+;Id3f/+, and Tnfrsf11aCre+; Id3f/f littermates. The figures display is representative of isovolume-rendered images from n = 4 independent experiments. Circles in bar graphs represent quantification of F4/80+ cells per square millimeter in individual fields of view. See also fig. S12. (E) Flow cytometry analysis of F4/80+ macrophages in the liver, brain, and kidney from 4-week-old Id3−/− and Id3+/− mice. Circles represent individual mice from n = 2 independent experiments. (F) Immunostaining with antibodies against ID3 (red), F4/80 (cyan) and YFP (green) on cryosections from the livers of 4-week-old Csf1rMeriCreMer; Rosa26LSL-YFP (OH-TAM at E8.5) mice. Nuclei are counterstained with DAPI (white). Scale bars, 5 μm. Images are representative of three independent experiments. (G) Immunostaining with antibodies against F4/80 (green) phospho-histone 3 (PHis3, red) and DAPI (gray) on liver cryosections from 4-week-old Id3−/− and Id3+/− mice. Scale bars, 10 μm. Individual dots in the bar graph represent mean percentages of F4/80+ cells labeled with PHis3 from five sampling areas (830 μm2) per liver from n = 5 mice per genotype. (H) Scatter-plot comparison of gene expression of 3-week-old Id3−/− and Id3+/− Kupffer cells performed on bulk RNA-seq data. Both axes (in log2 scale) represent normalized gene expression values (average value from three Id3+/− and two Id3−/− replicates). Red circles mark the threefold cut-off in both directions in gene expression level. Top GO terms for genes enriched in either Id3+/− or Id3−/− are indicated. P values were calculated using Student’s t test. See also fig. S12.

Conclusions

We show here that EMPs rapidly differentiate into a population of cells that we call pMacs, because they simultaneously colonize the whole embryo from E9.5 in a Cx3cr1-dependent manner while differentiating into macrophages. pMacs do not yet have macrophage morphology but are in the process of establishing a full core macrophage differentiation program that includes cytokine receptors, phagocytic and pattern recognition receptors, and complement. Starting from E10.5, almost immediately after colonization of embryonic tissues, tissue-specific expression of transcriptional regulators in F4/80+ macrophages initiates their specification into adult-type resident macrophages (Fig. 8). Our data, therefore, suggest that a broad or core macrophage program is progressively restricted or refined to a tissue-specific one in response to the absence and presence of tissue-specific cues, rather than the alternative possibility that committed subsets of early macrophages, pMacs, or even EMPs choose their future tissue of residence. The present results are consistent with an EMP origin of tissue-resident macrophages (2), although they do not exclude the additional contribution of primitive precursors (23), fetal HSCs, or a second wave of EMPs (11, 12), which would adopt a similar differentiation program when entering tissues. However, the transcriptional data set described here may lead to reinterpretation of fate-mapping models used to characterize the contribution of such later waves. For example, the progressive YFP expression by fetal macrophages in S100a4Cre; Rosa26LSL-YFP mice (12) is compatible with our present findings and does not require the existence of a second wave of precursors because expression of S100a4 and several other S100a family members was found to be part of the pMac transcriptional profile (fig. S13).

Fig. 8 Graphic summary of the establishment of the core macrophage program and subsequent specification.

Cx3cr1 is expressed in pMacs and is important for colonization of the embryo. Id3 is a liver macrophage–specific gene, which is essential for Kupffer cell development.

We propose that establishment of a core macrophage differentiation program in EMP-derived pMacs, as they colonize the embryo, is followed by the initiation of macrophage specification via the acquisition of tissue-specific transcriptional regulators, such as Id3 in Kupffer cells. This process is initiated simultaneously in the whole embryo during the first 2 days of organogenesis, from E8.5 to E10.5. Differentiation of resident macrophages is thus a developmental process and an integral part of organogenesis, independent of postnatal changes in the environment in the kidney, liver, and brain but not in the skin or lung. These results also identify a developmental window where the molecular mechanisms of macrophage specification can best be studied, tools to selectively label resident macrophages, and sets of tissue-specific transcriptional regulators expressed by developing and adult macrophages that may control their differentiation, maintenance, and function. Additionally, our results elucidate a molecular road map that will support efforts to differentiate specialized macrophages—microglia, Kupffer cells, kidney macrophages, alveolar macrophages, and Langerhans cells—in vitro from multipotent progenitors. Finally, our work provides a framework to analyze and understand the consequence(s) of genetic variation for macrophage contribution to disease pathogenesis in different tissues.

Materials and methods

Mice

Csf1rMeriCreMer, Csf1riCre, Rosa26YFP-LSL, Cx3cr1gfp/+ and Tnfrsf11aCre mice (2, 5, 22, 35) and Id3−/−; Id1fl/fl and Id3f/f mice—kindly provided by R. Benezra (34, 36, 37)—were maintained under SPF conditions. Animal procedures were performed in adherence to our project license issued by the United Kingdom Home Office under the Animals (Scientific Procedures) Act 1986 and by the Institutional Review Board (IACUC 15-04-006) from MSKCC. Genotyping was performed according to protocols described previously for Csf1riCre (38) Csf1rMeriCreMer (5), Id3−/−; Id1fl/fl (36), and Tnfrsf11aCre (22) mice. Genotyping of Cx3cr1gfp/+ was performed as recommended by Jackson Laboratory, Bar Harbor, Maine with the following primers: Cx3cr1 F 5′-CCC AGA CAC TCG TTG TCC TT-3′, Cx3cr1 R 5′-GTC TTC ACG TTC GGT CTG GT and Cx3cr1 R mut 5′CTC CCC CTG AAC CTG AAA C-3′. Cre recombination in Csf1rMeriCreMer; Rosa26LSL-YFP embryos was induced by single injection at E8.5 of 75 mg per g (body weight) of 4-hydroxytamoxifen (OH-TAM, Sigma) into pregnant females as described (2). The OH-TAM was supplemented with 37.5 mg per gram (body weight) progesterone (Sigma) to counteract the mixed estrogen agonist effects of tamoxifen, which can result in fetal abortions. Embryonic development was estimated as previously considering the day of vaginal plug formation as 0.5 days post-coitum (dpc), and staged by developmental criteria. E9: 20 to 25 sp, E9.5: 26 to 29 sp, E10.25: 30 to 35 sp, E10.5: 36 to 44 sp.

Preparation of cell suspensions and cell sorting

Pregnant females were killed by cervical dislocation or by exposure to CO2. Embryos ranging from E9 to E18.5 were removed from the uterus, washed in 4°C phosphate-buffered saline (PBS, Invitrogen) and dissected under a Leica M80 microscope. Yolk sacs (YS) were harvested from embryos between E9 and E10.5. To obtain single-cell suspensions for FACS sorting, tissues were included in cold PBS and mechanically disrupted under a 100um filter. Postnatal tissues were collected following the same procedure. For collection of Langerhans cells, epidermal sheets (from E18.5 to P21) were separated from the dermis after incubation for 45 min at RT in 4.8 mg/ml of dispase (Invitrogen), 3% fetal calf serum (FCS, Invitrogen) and 1 μM of flavopiridol (Sigma). The epidermis was further digested for 30 min at RT in PBS containing 2 mg/ml of collagenase D (Roche), 200 U/ml DNase I (Sigma), 4.8 mg/ml of dispase (Invitrogen), 3% FCS (Invitrogen) and 1uM of flavopiridol (Sigma) followed by mechanical disruption under a 100um filter. For blood phenotyping, mice were anaesthetized and blood was collected by cardiac puncture. For flow cytometry experiments, organs were incubated in PBS containing 1 mg/ml collagenase D (Roche), 100 U/ml DNase I (Sigma), 2.4 mg/ml of dispase (Invitrogen) and 3% FCS (Invitrogen) at 37°C for 30 min prior to mechanical disruption. Epidermal sheets were obtained as previously described (2). For embryonic tissue incubation time at 37°C was reduced to 20 min. Cell suspensions were centrifuged at 320g for 7 min, resuspended in FACS buffer (PBS, 0.5% BSA and 2 mm EDTA) containing purified anti-CD16/32 (FccRIII/II) (1:100 dilution) and incubated for 15 min at 4°C. Samples were immunostained with antibodies mixes for 30 min at 4°C, and analyzed by flow cytometry using a LSR Fortessa X-20 (BD-Bioscience). Single live cells were gated on the basis of dead cell exclusion (Hoechst 33258), side (SSC-A) and forward scatter (FSC-A) gating, and doublet exclusion using forward scatter width (FSC-W) against FSC-A. The full list of antibodies used can be found in table S4. Cell numbers per organ or per gram of tissue were calculated as follow. In E9.5 to E12.5 embryonic tissue cell suspensions were prepared, stained, and acquired from whole organs, and the number of live cells per tissue was directly obtained from FCS files. In E14.5 to E18.5 embryos, cells suspensions were prepared as above, but half were stained and acquired. In adult mice, organs were weighted, cell suspensions were prepared from 100 to 500 mg of tissues, and the number of live single F4/80+ cells per gram of tissue was determined using a cell counter (GUAVA easyCyte HT).

Cell sorting was performed using an Aria II BD cell sorter. Cell suspensions were prepared as above, except for the gating of single live cells which was performed using dead cell exclusion with DAPI. For bulk sequencing, EMPs were identified after gating on Kit+CD45lo cells based on AA4.1 expression. pMacs were identified after gating on KitCD45+ based on CD11b expression and no expression of F4/80. Additionally, Gr1+ or Ter119+ cells were excluded from the F4/80CD11b+ gate. Macrophage populations were identified after gating on CD45 based on expression of F4/80 and CD11b. 100 cells for each sample were directly sorted into a 96 well plate (Biorad) in 4 μl of H2O containing 0.2% of triton TXT (Sigma) and 0.8 U/μl of RNAse inhibitor (Clontech), and processed as indicated below (Generation and analysis of “bulk” transcriptomes from candidate EMP, pMac, and macrophage populations from E9 to P21). For single-cell sequencing, all CD45low/+ single cells from a E10.25 (30 to 34 sp) embryo proper were sorted into 384-well plates filled with 2 μl lysis buffer (Triton-X 0.2%) (Sigma) in molecular biology grade H2O (Sigma) supplemented with 0.4 U/μl protein-based RNase inhibitor (Takara) and barcoded poly-T primers, and processed as described below (Generation and analysis of single-cell transcriptomes).

Cytology

Cells were collected into FCS and centrifuged (800 rpm, 10 min, low acceleration) onto Superfrost slides (Thermo Scientific) using a Cytospin 3 (Thermo Shandon). Slides were air-dried for at least 30 min, and fixed for 5 min in methanol, stained in 50% May-Grünwald solution for 5min, 14% Giemsa for 15 min, washed with Sorensons buffered distilled water (pH 6.8) for 5 min and rinsed with Sorensons buffered distilled water (pH 6.8). Slides were air-dried and mounted with Entellan New (Merck) and representative pictures were taken using an Axio Lab.A1 microscope (Zeiss) under a N-Achroplan 100x/01.25 objective.

Immunofluorescence, imaging, analysis, and illustrations

E10.25 Embryos were fixed for 4 hours in 4% formaldehyde (Sigma) under agitation and >E10.25 embryos were fixed overnight. After fixation, embryos were incubated overnight in 30% sucrose and embedded in OCT compound (Sakura Finetek). Cryoblocks were cut at a thickness of 10-12 μm and then blocked with PBS containing 5% normal goat serum (Invitrogen); 1% BSA (w/v); 0.3% Triton X-100 for 1 hour at room temperature. Samples were incubated overnight at 4°C with rat anti-mouse F4/80 (1:200, Biorad), rabbit anti-mouse Iba1 (1:200; Wako), chicken anti-GFP for YFP detection (1:500, Invitrogen), rat anti-mouse Dectin-1 (1:200, Biorad), anti-Granulin (1:200, abcam), rat anti-mouse CD206 (1:200, Biorad), Armenian hamster anti-mouse PECAM-1 (1:300, Thermo Scientific), rabbit anti-mouse/human ID3 (1:500, Biocheck), rabbit anti-mouse ID1 (1:500, Biocheck), rat anti-mouse CD16/CD32 (BD Biosciences), goat anti-mouse Trem2 (1:200, Abcam), armenian hamster anti-mouse CD119 (1:200, Clone 2E2, eBioscience), armenian hamster anti-mouse CD120b (1:200, clone 55R-286, Biolegend). Secondary antibodies used were anti-rabbit Cy3 (1:500; Invitrogen), anti-chicken Alexa Fluor 488 (1:500; Invitrogen), anti-rat Alexa Fluor 555 (1:500; Invitrogen), anti-rat Alexa Fluor 488 (1:500; Invitrogen), anti-rat Alexa Fluor 647 (1:500; Invitrogen), anti-goat Alexa Fluor 568 and anti-armenian hamster Dylight 649 (Jackson ImmunoResearch Laboratories). Samples were then mounted with Fluoromount mounting medium with DAPI (eBiosciences) and visualized using a LSM880 Zeiss microscope with 20x/0.5 (dry), performing a tile scan and Z-stack on whole embryos or tissues. Image analysis and cell quantification was performed using Imaris (Bitplane) software. To determine the CD31+ area in Id3+/− and Id3−/− liver sections, maximum intensity Z-projections pictures were converted into binary images and the CD31+ area was measured using Image J (NIH, Bethesda, MD, USA) (39). Results were normalized per mm2 of tissue. Illustrations were created by adapting templates from Servier Medical Arts (www.servier.com/Powerpoint-image-bank, licensed under a Creative Commons Attribution 3.0 Unported License).

Generation and analysis of Kupffer cell transcriptome in Id3−/− and littermate controls

For qRT-PCR experiments, cells were sorted as described above. Hepatocytes were enriched by centrifugation of the whole liver cell suspension at 50 g for 3 min (Sorvall Legend XTR centrifuge). Supernatant was taken for further staining of macrophages (CD45+, CD11blo, F4/80+). Hepatocytes were sorted using the FSC-A and SSC-A gate with subsequent exclusion of doublets and CD45+ cells. Cells were sorted into RNA lysis buffer and RNA extraction was performed as per manufacturers protocol (Macherey-Nagel). cDNA was synthesized using the QuantiTect Reverse Transcription Kit (Qiagen) as per manufacturers protocol. qRT-PCR was performed on a QuantStudio 6 Flex using TaqMan Fast Advanced Master Mix (Applied Biosystems) and TaqMan probes for Id3 (Mm00492575_m1), Nr1h3 (Mm00443451_m1), and GAPDH (Mm99999915_g1) (Life Technology). For RNA-seq analysis of Id3+/− and Id3−/− Kupffer cells, cells were sorted into Trizol as described above. RNA from sorted cells was extracted using RNeasy mini kit (Qiagen) according to instructions provided by the manufacturer. After ribogreen quantification and quality control of Agilent BioAnalyzer, 400 pg of total RNA underwent amplification (12 cycles) using the SMART-seq V4 (Clonetech) ultra low input RNA kit for sequencing. 10 ng of amplified cDNA was used to prepare Illumina hiseq libraries with the Kapa DNA library preparation chemistry (Kapa Biosystems) using 8 cycles of PCR. Samples were barcoded and run on a Hiseq 2500 1T in a 50 bp/50 bp Paired end run, using the TruSeq SBS Kit v3 (Illumina). An average of 54 million paired reads were generated per sample and the percent of mRNA bases was closed to 77% on average. The output data (FASTQ files) were mapped to the target genome using the rnaStar aligner (40) that maps reads genomically and resolves reads across splice junctions. 2 pass mapping method was used, outlined in Engstrom et al. (41) in which the reads are mapped twice. The first mapping pass uses a list of known annotated junctions from Ensemble. Novel junctions found in the first pass are then added to the known junctions and a second mapping pass is done (on the second pass the RemoveNoncanoncial flag is used). After mapping we post process the output SAM files using the PICARD tools to: add read groups, AddOrReplaceReadGroups which in additional sorts the file and coverts it to the compressed BAM format. We then compute the expression count matrix from the mapped reads using HTSeq and one of several possible gene model databases. The raw count matrix generated by HTSeq are then processed using the R/Bioconductor package DESeq, which was used to both normalize the full data set and analyze differential expression between sample groups. Gene Ontology analysis was performed using the GO analysis function in GeneSpring GX 13.0 (Agilent), with the P value calculated using a hypergeometric test with Benjamini–Yekutieli correction. For that, genes with a fold change difference of ±2 between Id3+/− and Id3−/− Kupffer cells were selected. Significantly regulated genes (t test P < 0.05; FDR < 0.05) from this selection were grouped into GO terms.

Generation and analysis of “bulk” transcriptomes from candidate EMP, pMac, and macrophage populations from E9 to P21

RNA isolation and library construction

cDNA synthesis and enrichment was performed following the Smart-seq2 protocol as described (42, 43). ERCC spike-in RNA (Ambion) was added to the lysis buffer in a final dilution of 1:1,000,000. After the cDNA was synthesized and amplified from single cells using 18 cycles, quantitative PCR was performed with GoTaq-PCR master mix (Promega) on a C1000 Touch Thermal Cycler qPCR instrument (Bio-Rad) to test for house keeping gene expression. Library preparation was conducted on 1ng of cDNA using the Nextera XT library preparation kit (Illumina) as described (43). Sequencing was performed by the Biomedical Sequencing Facility at CeMM using a 50 bp single-read setup on the Illumina HiSeq 2500 platform.

RNA-seq analysis

We first trimmed off sequencing adapter from the reads generated, and then aligned the reads using Bowtie v 1.1.1 (44) (parameters: -q -p 6 -a -m 100 --minins 0 --maxins 5000 --fr --sam --chunkmbs 200) to the cDNA reference transcriptome (mm10 cDNA sequences from Ensembl). For genome browser track visualization, we generated a second alignment with Tophat2 (45) (v 2.0.13, parameters: --b2-L 15 --library-type fr-unstranded --mate-inner-dist 150 --max-multihits 100 --no-coverage-search --GTF) against the reference genome (mm10). Next, we removed duplicate reads before quantifying transcript levels with BitSeq (46) (v 1.12.0). The raw transcript counts were loaded into R and processed further. At this stage, we removed samples with a substandard quality, that is, all samples that had less than 2 million reads, less than 33% of reads aligned, or less than 5,000 transcripts detected (with ≥ 25 reads). 20 out of 178 data sets failed these criteria (11.2%). Of the remaining data sets, we merged technical replicates creating 93 unique biological samples. We took forward only transcripts that were detected reliably (≥50 reads) in at least four samples (n = 37,521). For statistical analysis, we used raw read counts as input for DESeq2 (20), factoring in the flowcell identifier as a covariate to reduce the effect of technical variation. To identify genes of particular interest to the development of tissue-resident macrophages over time and in different tissues, we performed two types of comparisons: (a) Cell type-specific: Pairwise comparisons (Wald test) between EMPs, pMacs, and macrophages (up to E10.5) independent of their tissue of origin (treating time and tissue as covariates). (b) Tissue-specific: Pairwise comparisons (Wald test) between macrophages from all tissues stratified and stage (post-natal). We considered genes with an FDR-corrected P ≤ 0.05 as differentially expressed. For visualization and illustration purposes (PCA, supplementary tables, heatmaps), we used values adjusted by variance stabilizing transformation from DESeq2 in which batch-effects had been corrected for with ComBat (47).

Evaluation of lists of differentially expressed genes from the RNA-seq data

We bioinformatically investigated the differentially expressed transcripts (table S5) from our statistical comparisons in several ways: To identify transcription factors specifically regulating each set of transcripts, we used LOLA (26) (v 0.99.4) together with its core database of ChIP-seq binding peaks from CODEX (48) to identify enrichment of experimentally-derived transcription factor binding locations in a window around the promoter regions (TSS ± 20 kb) of differentially expressed transcripts. We corrected for multiple testing using the Benjamini and Yekutieli method. To visualize and summarize the expression patterns of many genes (lists of differentially expressed genes) in many different conditions (different tissues at different time points) and across replicates, we sought to use an adaptation of lineage scorecards (21). Briefly, we considered each list of differentially up-regulated genes as a set of marker genes and determined the relative enrichment of these marker sets in each individual condition (tissue by cell type by time) in comparison to all other data sets using a modified version of parametric gene set enrichment analysis in R (package: PGSEA). We also used GSEA to test for the relative overrepresentation of gene signatures in sorted gene lists. To this end, we extracted lists of genes sorted by logarithmic fold change between the mean expression levels in any one tissue Mac, pMac, or EMP sample stratified by stage (early = E9 to E10.5, fetal = E12.5 to E18.5, postnatal = P2 to P21) compared to all other samples at the same stage. Additionally, we extracted lists of all genes we found differentially up-regulated in any tissue-specific signature, cell type–specific signature, or in lists of genes extracted from (7, 25, 28). External gene lists were translated to Ensembl Transcript IDs using g:profiler (49). Both sets of data were loaded into and analyzed using the GSEA tool (50) and the results read and summarized using the metaGSEA R library. We also incorporated ChIP-seq binding profiles (bigWig) for factors identified in the LOLA analysis from the CODEX database. Heatmaps were generated using GeneSpring GX 13.0 (Agilent). All other analyses and plotting was performed in R.

Generation and analysis of single-cell transcriptomes

For single-cell RNA sequencing the MARS-Seq approach described by Jaitin et al. (51) was applied using the Biomek FXP system (Beckman Coulter). In brief, single cells were sorted into each well of 384-well plates filled with 2 μl lysis buffer [Triton-X 0.2% (Sigma) in molecular biology grade H2O (Sigma) supplemented with 0.4 U/μl protein-based RNase inhibitor (Takara)] and barcoded poly-T primers [400 nM (IDT); for barcode details see (51)]. Samples were pre-incubated (3 min at 80°), reverse transcriptase mix (10 mM DTT (Invitrogen), 4 mM dNTPs (NEB), 2.5 U/μl Superscript III RT enzyme (Invitrogen) in 50 mM Tris-HCl (pH 8.3; Sigma), 75 mM KCl (Sigma), 3 mM MgCl2 (Sigma), ERCC RNA Spike-In mix (Life Technologies) at 1:80 × 107 dilution per cell) was added to each well and mRNA was reverse transcribed to cDNA (2 min at 42°C, 50 min at 50°C, 5 min at 85°C). Excess primers were digested [Exo I (NEB); 37°C for 30 min then 10 min at 80°C] followed by a 1.2 x SPRI bead (Beckman Coulter) cleanup, samples were pooled and second strands synthesized (second strand synthesis kit (NEB); 2.5 hours at 16°C) followed by a 1.4 x SPRI bead (Beckman Coulter) cleanup. Samples were linearly amplified by T7-promoter guided in vitro transcription (T7 High Yield RNA polymerase IVT kit (NEB); 37°C for 12 hours). DNA templates were digested [Turbo DNase I (Ambion); 15 min at 37°C] followed by a 1.2 x SPRI bead (Beckman Coulter) cleanup and the RNA was fragmented [Zn2+ RNA fragmentation solution (Ambion); 1.5 min at 70°C] followed by a 2 x SPRI bead (Beckman Coulter) cleanup. Barcoded ssDNA adapters [IDT; for barcode details see (51) were ligated to the fragmented RNA (9.5% DMSO (Sigma), 1 mM ATP, 20% PEG8000 and 1 U/μl T4 RNA ligase I (NEB) in 50 mM Tris HCl pH 7.5 (Sigma), 10 mM MgCl2 and 1 mM DTT; 22°C for 2 hours) and a second RT reaction (Affinity Script RT buffer, 10 mM DTT, 4 mM dNTP, 2.5 U/μl Affinity Script RT enzyme (Agilent); 2 min at 42°C, 45 min at 50°C, 5 min at 85°C) was performed followed by a 1.5 x SPRI bead (Beckman Coulter) cleanup. Final libraries were generated by subsequent nested PCR reaction (0.5 μM of each Illumina primer [IDT; for primers details see (51) ] and KAPA HiFi HotStart ready mix (Kapa Biosystems) for 15 cycles according to manufacturer’s protocol followed by a 0.7 x SPRI bead (Beckman Coulter) cleanup. Library quantity and quality were assessed using the Agilent 2200 Tapestation system and libraries were subjected to next generation sequencing using an Illumina HiSeq1500 instrument (PE with no index; read 1: 61 bases (3 bases random nucleotides, 4 bases pool barcode, 53 bases specific sequence), read 2: 13 bases [6 bases cell barcode, 6 bases unique molecular identifier)].

Preprocessing, quality assessment, and control of single-cell transcriptome data

From sequenced data, pool barcodes, cell specific tags and Unique molecular identifiers (UMI) were extracted (576 cells sequenced). Subsequently, sequencing reads with ambiguous plate/cell-specific tags or UMI sequence with low quality (Phred < 27) and reads which map to E. coli were eliminated using Bowtie with parameters “-M 1 -t --best --chunkmbs 64 –strata.” Next, fastq files were demultiplexed using the fastx_barcode_splitter from the fastx_toolkit and R1 reads (after trimming of pool barcode sequences) were mapped to the mouse mm10 and ERCC pseudo genome assembly using Bowtie “-m 1 -t --best --chunkmbs 64 –strata.”

Valid reads were then counted using unique molecular identifiers if they mapped to the exon based gene model derived from Ensembl’s biomart, mm10. Following this, a gene expression matrix was generated containing the number of unique UMIs associated with valid reads for every cell and every gene. Additionally, UMI sequencing errors were corrected for and filtered as described in (51).

Filtering single-cell transcriptome data

In order to avoid biases introduced by low quality data we performed the following filtering of single-cell data. Removal of cells with a ratio of mitochondrial versus endogenous genes exceeding 0.15 and cells with less than 320 molecule counts or less than 150 unique genes were removed from the analysis. Prior to analysis expression tables were filtered for mitochondrial, ribosomal and predicted genes to reduce noise.

Analysis of MARS-seq single-cell transcriptome data

Analysis of the normalized and filtered single-cell gene expression data (8657 genes across 408 single-cell transcriptomes used in the final expression table) was done using several functions of the SEURAT single-cell analysis package (52) and Monocle 2 (53). First, highly variable genes were determined as genes exceeding the dispersion threshold of 0.75. To infer the structure of the gene expression data a PCA was performed on the basis of highly variable genes. Following tSNE DBScan clustering was performed to identify clusters of cells. The optimal number of clusters was identified by calculating several cluster indices by the NbClust R package (54). Relative expression for a cell was calculated as gene expression of a gene/gene set in relation to the total molecule counts in this cell. To identify clusters within the MARS-seq data the relative gene expression profiles of cell type specific gene signatures were overlaid to the tSNE plots. Pseudotime analysis was performed by the Monocle 2 algorithm by genes exceeding the average expression cutoff of 1 while having an empirical dispersion higher than 1. To analyze expression of single genes, relative expression was overlaid onto the tSNE plots or visualized as dot plots.

Generation of cell signatures for analysis of single-cell data

In order to unambiguously identify cell state specific genes for EMPs, pMacs and early macrophages we generated exclusive gene signatures. Here, differentially expressed (DE) genes were identified by a one-way ANOVA model [|FC| > 1.4, FDR-adjusted P < 0.05 (55)] between EMPs (E10.25 YS and fetal liver), pMacs (E10.25 YS, fetal liver, head, and limbs) and early macrophages (E10.25 and E10.5 YS, fetal liver, head, and limbs). The nonoverlapping DE genes between these three contrasts were chosen as exclusive gene signatures for further analysis.

Generation of tissue signatures for analysis of single-cell data

To identify genes that are up-regulated in tissue macrophages in relation to early macrophages a one-way ANOVA model [|FC| > 1.5, FDR-adjusted P < 0.05 (31)] was calculated. Up-regulated nonoverlapping DE genes between early macrophages (E10.25 and E10.5, fetal liver, head, and limbs/skin) versus brain macrophages (E14.5 and E18.5 brain), liver macrophages (E14.5 and E18.5 liver) or limb/skin macrophages (E14.5 and E18.5 limb/skin) were used as tissue macrophage specific gene signatures. Furthermore, a common early macrophage signature was defined as being up-regulated in early macrophages (E10.25 and E10.5 fetal liver, head, and limbs) versus all other late tissue macrophages (E14.5 and E18.5 liver, head, and limbs/skin). We assessed enrichment of these signatures in scRNA-seq data from pMacs by calculating an relative enrichment score for each signature in pMacs (molecule count of signature genes/(total molecule count × number of signature genes).

Enrichment of tissue signatures in scRNA-seq data

To assess enrichment of tissue macrophage-specific signatures or differentiation signatures in scRNA-seq data from pMacs, we identified genes with multimodal expression using Hartigan’s Dip test statistic for unimodality. Next, multimodal genes were grouped by hierarchical clustering and enrichment of signatures within the clusters tested using hypergeometric testing with FDR correction (Benjamini-Hochberg).

Supplementary Materials

www.sciencemag.org/content/353/6304/aaf4238/suppl/DC1

Supplementary Text

Figs. S1 to S13

Tables S1 to S6

References and Notes

  1. Acknowledgments: We thank R. Benezra (MSKCC, New York) for the Id3−/− and Id3f/f strains. We thank the Biomedical Sequencing Facility at CeMM and the MSKCC Integrated Genomics and Bioinformatics Cores for assistance with next-generation sequencing. The data presented in this manuscript are tabulated in the main paper and the supplementary materials. Sequencing data sets described in this work have been submitted to the Gene Expression Omnibus (GEO) repository (GEO accession number GSE81774). Additionally, a genome browser track hub and additional supplementary data are available at http://macrophage-development.computational-epigenetics.org/. This work was supported by the National Cancer Institute of the U.S. NIH (grant P30CA008748) and by investigator awards from the Wellcome Trust (WT101853MA) and the European Research Council (2010-StG-261299) to F.G. E.M. is supported by a European Molecular Biology Organization long-term fellowship (ALTF 530-2015). F.H. is supported by a Deutsche Forschungsgemeinschaft (DFG) postdoctoral fellowship (HA 7723/1-1). J.K. is supported by a Doktorand/inn/enprogramm der Österreichischen Akademie der Wissenschaften (DOC) Fellowship of the Austrian Academy of Sciences. C.B. is supported by a New Frontiers Group award from the Austrian Academy of Sciences. M.B. and J.L.S. are members of the Excellence Cluster ImmunoSensation. M.B. is supported by a DFG grant (BE 4427/3-1). J.L.S. is supported by DFG grants SFB 704, SFB 645, INST 217/575-1, INST 217/576-1, and INST 217/577-1.
View Abstract

Navigate This Article