Overview

The massive data set in S1 Data presents expression estimates for 18,175 transcripts in 466 RNA-seq data sets from MPS cells in a wide diversity of differentiated states. This provides the framework for a meta-analysis and review of the current state of knowledge of MPS cell differentiation. We first consider the expression profiles of individual transcripts, a network analysis of the relationships between MPS cell populations, and identification of stringently coregulated clusters of transcripts. This analysis leads to detailed consideration of the content of individual clusters enriched in tissue-specific MPS populations and the recognition of likely artefacts associated with isolation of MPS cells. Separate analysis of expression of transcripts encoding transcription factors addresses the control network that regulates differential gene expression in the MPS, whilst analysis of solute carriers and metabolic pathways contributes to the emerging field of immunometabolism. scRNA-seq is a rapidly emerging approach to identification of MPS subpopulations. Based upon comparative analysis, we critically review the validity of conclusions about MPS heterogeneity based upon scRNA-seq and the relationship with coregulated clusters identified by network analysis. In the light of our analysis, we question the validity of cell surface markers in analysis of MPS cell diversity, especially as they relate to the separate identity of DC and macrophage polarisation states.

Expression profiles of individual transcripts

To survey the heterogeneity of the MPS cells (monocytes, macrophages, and DCs), we first considered the expression profiles of selected individual transcripts, including candidate housekeeping genes and surface markers commonly used in studies of this lineage. The choice of appropriate reference genes for quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) quantification of RNA levels is a significant issue in many studies, including macrophage differentiation. For example, Stephens and colleagues [73] proposed the use of a weighted geometric average of the most stably expressed genes for studies of differentiating macrophages and osteoclasts. Fig 1A shows the expression profiles of candidate housekeeping genes (Hprt, Actb, B2m, Gapdh, Ppia) commonly used in qRT-PCR as reference genes. We envisaged that these transcripts would be relatively consistent between data sets and would be correlated with each other. However, each of these transcripts varied by >10-fold among the MPS populations, and in pairwise analysis, they were only weakly correlated (Pearson correlations in Fig 1B and Spearman correlations in S1 Fig). S1 Data also includes the mean, standard deviation, and coefficient of variance (CoV) of each transcript across the entire data set. Only 200 transcripts had a CoV < 0.5. This analysis indicates that pathways normally considered as housekeeping (intermediary metabolism, protein synthesis, endoplasmic reticulum (ER) and Golgi membrane trafficking and secretion, endocytosis, cytoskeleton, etc.) are independently regulated in MPS cells. There are few, if any, identifiable housekeeping genes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Expression of housekeeping genes across MPS cell populations.

The data underlying this Figure can be found in S1 Data. (A) Expression patterns across cells from different tissues. Each column represents a sample. Upper bar along the X axis shows the cell type (black—monocytes and macrophages; red—DCs). Lower bar shows the tissue, coloured as shown in the key. Y axis shows expression level in TPM, calculated using Kallisto. (B) Correlations (Pearson correlation coefficient) between expression patterns of different housekeeping genes. DC, dendritic cell; MPS, mononuclear phagocyte system; TPM, transcripts per million.

https://doi.org/10.1371/journal.pbio.3000859.g001

Fig 2A shows the expression pattern of transcripts encoding surface markers used to separate some of the subpopulations herein: Adgre1 (encoding F4/80), Cd4, Cd74 (Class II MHC), Csf1r (CD115), Cx3cr1, Fcgr1 (CD64), Icam2, Itgax (CD11C), Lyve1, Mertk, Mrc1 (CD206), and Tnfrsf11a (RANK). Fig 2B shows a summary of the Pearson correlations between them; Spearman correlations are shown in S1 Fig. Consistent with studies using Csf1r reporter transgenes [74, 75], Csf1r mRNA was universally expressed in MPS cells, albeit with significant variation in level, being highest in microglia and lowest in cDC1s. Csf1r was correlated (Pearson r > 0.5, Spearman r > 0.65) with Adgre1, Fcgr1, Cx3cr1, Mertk, and Tnfrsf11a, but these transcripts were less correlated with each other. Mrc1 was reported to be correlated with expression of Lyve1 and inversely with Class II MHC transcripts [27, 76]. Across the entire spectrum of macrophage transcriptomes, Mrc1 was correlated with Lyve1 but was more widely expressed (Fig 2A). However, there was no evidence of an inverse correlation between Mrc1 and Cd74 or other Class II MHC-associated transcripts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Expression of cell surface marker genes across MPS populations.

The data underlying this Figure can be found in S1 Data. (A) Expression patterns across cells from different tissues. Each column represents a sample. Upper bar along the X axis shows the cell type (black—monocytes and macrophages; red—DCs). Lower bar shows the tissue, coloured as shown in the key. Y axis shows expression level in TPM, calculated using Kallisto. (B) Correlations (Pearson correlation coefficient) between expression patterns of different MPS genes. DC, dendritic cell; MPS, mononuclear phagocyte system; TPM, transcripts per million.

https://doi.org/10.1371/journal.pbio.3000859.g002

Network analysis of relationships of MPS populations and expressed transcripts

To determine whether any transcripts encoding surface markers were correlated with cellular phenotype, we used the graph-based network analysis tool BioLayout. Fig 3 presents a sample-to-sample correlation matrix generated using the Fruchterman–Rheingold algorithm in BioLayout showing the clear segregation of the MPS populations based on the tissues from which they were isolated (Fig 3A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Sample-to-sample network analysis of gene expression in MPS cell populations.

Each sphere (node) represents a sample, and lines between them (edges) show Pearson correlations between them of ≥0.68 (the maximum value that included all 466 samples). (A) Samples coloured by tissue of origin. (B) Samples coloured by cell type. (C) Samples coloured by BioProject. MPS, mononuclear phagocyte system.

https://doi.org/10.1371/journal.pbio.3000859.g003

Initial iterations of the sample-to-sample analysis identified 3 technical issues that served to validate our approach and data quality. We initially included data from the large ImmGen UL1 project (GSE127267l; GSE124829; see [77]), but this project uses an ultra-low–input RNA-seq pipeline based upon 1,000 sorted cells and the scRNA-seq platform Smartseq2. The initial sample-to-sample analysis revealed a large batch effect for these samples relative to all other samples, and we therefore excluded these data. Secondly, we mistakenly included populations of stellate cells, endothelial cells (ECs), and hepatocytes profiled within a study of the influence of the hepatic niche on Kupffer cell (KC) development [78]. The replicate samples of these non-MPS cells formed 3 entirely separate clusters each distinct from the main MPS cluster, and these were also excluded. Thirdly, the initial analysis indicated that macrophages isolated from the choroid plexus clustered with resident peritoneal macrophages. Consultation with the authors of the study [59] confirmed that this was due to a 3-way error in upload to the SRA, and the samples labelled choroid plexus macrophages were indeed resident peritoneal macrophages; those labelled peritoneal macrophages were actually microglia, and those labelled microglia were choroid plexus macrophages. The discovery of these technical issues highlights the power and validity of the network approach as implemented in BioLayout.

Consistent with previous analysis of microarray data sets, in which all MPS cells, including blood monocytes and DCs, clustered together and were clearly distinct from other haematopoietic and nonhaematopoietic lineages [39, 41, 42, 44], the isolated spleen, lung, and lymph node (LN) DC subpopulations did not form a separate element in the network (red nodes in Fig 3B). Based upon their overall transcriptomic profile, the DCs were no more divergent from other MPS populations than the isolated monocytes and macrophages purified from different tissues were from each other. The apparent relationship to BioProject (Fig 3C) occurs mainly because most studies were focussed on a particular tissue or cell type. There may also be minor impacts from differing methods of extracting and processing RNA and low-depth and single-end libraries compared with high-depth/paired-end libraries, but nonetheless, when different groups had profiled the same cell populations, the profiles were clustered together.

To further test the robustness of the associations between different MPS populations, we repeated the analysis with increasing r values of the Pearson correlation and based upon a Spearman correlation. Fig 4 presents the networks with nodes coloured by cell type at different Pearson and Spearman correlation coefficients. The networks with nodes coloured by tissue type and BioProject are presented in S2–S5 Figs. These networks clearly show that the DCs are more similar to monocytes and macrophages from the same tissue than they are to DCs from another tissue. Creating the network of samples at a Pearson r value of 0.85 clearly separated the microglia from the main network as the most divergent MPS cells (S2 Fig). Interestingly, the lung DCs were also separated from the lymphoid tissue DCs. At an r value of 0.95, the different MPS cells formed separate elements in the network (S3 Fig), but the cDC2 samples were still in the same element as macrophages from several nonlymphoid tissue sources. The close relationship of DCs with other MPS cells was equally evident when based upon Spearman correlation at the r threshold of 0.85 (S4 Fig) and at higher r (S5 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Sample-to-sample 2D network analysis of gene expression in monocyte, macrophage, and DC populations.

Each sphere (node) represents a sample, and lines between them (edges) show Pearson correlations between them. (A) Network laid out at Pearson correlation coefficient of ≥0.85. The network includes 458 samples. (B) Network laid out at Pearson correlation coefficient of ≥0.95. The network includes 418 samples. (C) Network laid out at Spearman correlation coefficient of ≥0.85. The network includes 443 samples. (D) Network laid out at Spearman correlation coefficient of ≥0.9. The network includes 427 samples. The networks with nodes coloured by tissue and BioProject are shown in S2–S5 Figs. DC, dendritic cell.

https://doi.org/10.1371/journal.pbio.3000859.g004

The GCN for the same data set was developed at an optimal Pearson r value of 0.75 chosen based on the graph of network size versus correlation threshold (shown in S6 Fig). Fig 5A shows the whole network, and Fig 5B highlights the tissue-specific clusters and those that contain markers of other cell types, as discussed below. S2 Data summarises the coexpressed gene clusters and the average gene expression profiles of the clusters containing at least 10 nodes (transcripts). The graphs are colour-coded to indicate the tissue origin and cell type as in Fig 1 (samples are listed in the Readme sheet of S2 Data). An additional sheet in S2 Data provides gene ontology (GO) term enrichment of the larger clusters. For ease of visualisation relative to sample information, profiles of surface markers and transcription factors discussed below are provided as an additional sheet in S1 Data. Table 2 provides an overview of the major functional clusters discussed in more detail below. It is beyond the scope of this study to analyse and cite published evidence related to every transcript in detail. In Table 2, individual genes from within the cluster have been included based their candidate role as transcriptional regulators and upon known associations with mononuclear phagocyte biology determined by PubMed search on ‘Genename’ AND ‘macrophage’ or ‘dendritic cell’. On the principal of guilt by association [38–42], there are hundreds of other genes within these clusters that have inferred functions in innate immunity and mononuclear phagocyte biology.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Description of major functional clusters of coexpressed genes in mouse MPS cell samples.

https://doi.org/10.1371/journal.pbio.3000859.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. GCN analysis of gene expression in MPS cell populations.

Each sphere (node) represents a gene, and lines between them (edges) show Pearson correlations between them of ≥0.75. Nodes were grouped into clusters with related expression patterns using the MCL algorithm with an inflation value of 1.7. Lists of genes and expression profiles of clusters are presented in S2 Data. (A) The network generated by the BioLayout analysis. Elements with ≥5 nodes are shown. Nodes are coloured by MCL cluster. Lists of genes and average expression profiles for all clusters are presented in S2 Data. Monocyte and macrophage genes (black ovals), DC genes (red oval). (B) Network showing only major clusters of monocyte and macrophage genes (black ovals), DC genes (red oval), and other cell types. DC, dendritic cell; GCN, gene coexpression network; MCL, Markov clustering algorithm; MPS, mononuclear phagocyte system; NK cell, natural killer cell.

https://doi.org/10.1371/journal.pbio.3000859.g005

Major macrophage-enriched coregulated clusters

At the chosen r threshold of 0.75, the GCN approach using the normalised data from multiple laboratories identified many coregulated clusters of transcripts that are known to be functional in MPS cells based upon prior knowledge inferred from smaller data sets. In the large transcriptional atlas projects across many species and tissues discussed in the introduction, the largest clusters tend to contain housekeeping genes that show relatively little variation amongst tissues. Consistent with the analysis of individual housekeeping genes in Fig 1 and analysis of the variance for the entire data set in S1 Data, there is no such cluster in the MPS data set.

Cluster 1 is a generic MPS cluster that drives the relatively close association between all of the samples, including the different subclasses of DCs, in the sample-to-sample network (Figs 3B and 4) and distinguishes MPS cells from other leukocytes. It includes Csf1r, Fcgr1, Cd68, Sirpa, Tnfrsf11a and the core myeloid transcription factor gene Spi1 alongside many other known MPS-enriched transcription factors [79, 80]. One notable inclusion is the glucocorticoid receptor gene, Nr3c1, which mediates transcriptional activation of a wide range of anti-inflammatory genes in macrophages [81]. As one might expect from the known endocytic and secretory activity of MPS cells, the cluster is enriched for GO terms related to endosome/lysosome and intracellular transport/secretion, which are major constitutive functions of mononuclear phagocytes [38] (S2 Data). Transcripts in Cluster 3 were also expressed widely in MPS cells, but the cluster has a distinct average expression profile. Cluster 3 includes genes encoding several forkhead transcription factors (Foxo3, Foxo4, Foxk1, and Foxk2), the key transcriptional regulators of autophagy [82–84], and Nfat5, which controls macrophage apoptosis [85]. This cluster also contains Mertk, the perforin-like immune effector gene Mpeg1, Aim2 (which encodes a sensor for cytoplasmic DNA [86]), and transcripts for numerous DEAH- and DEAD- box helicases all implicated in DNA sensing in innate immunity [87]. There are also members of the neuronal apoptosis inhibitory protein (NAIP) family of inflammasome regulators (Naip2, 5, 6); reviewed in [88]). We infer that this cluster of transcripts reflects an independently regulated capacity for innate immune recognition of internalised pathogens. Other than Mertk, there is no other plasma membrane marker associated with this set of candidate innate immune effector genes.

Genes in Cluster 4 were strongly expressed in samples from brain and include microglia-enriched markers that are depleted in brains of Csf1r-deficient mice and rats, such as Cx3cr1, Tmem119, P2ry12, and the key transcription factor genes Sall1, Sall2, and Sall3 [17, 89]. This cluster drives the separation of microglia as the most distinctive member of the MPS family. Cluster 9 contains the S phase transcription factor gene Foxm1 and numerous cell-cycle–associated transcripts [90], and the GO term enrichment supports a cell-cycle role. Genes of the cell-cycle cluster were expressed in all isolated MPS populations at various levels, consistent with evidence that they are capable of self-renewal in the steady state [5, 6, 13]. The separation of this cluster indicates that proliferative activity is not tightly linked to any MPS differentiation state or surface marker.

Identification of a capillary-associated expression cluster

Most macrophages and DCs included in this analysis were purified by FACS based upon their expression of specific markers including those shown in Fig 2 (see Table 1). Chakarov and colleagues [27] identified a population of pericapillary cells in the lung that expressed LYVE1 and extended their analysis to FACS-separated cells from fat, heart, and dermis. Their RNA-seq results are included in our data set. Based upon analysis of differentially expressed genes, the authors identified a set of genes with high expression in sorted LYVE1^hi macrophages relative to LYVE1^lo macrophages across the 4 tissues, including Mrc1, Timd4, Cd5l, Fcna, and Vsig4 [27, 76]. The GCN reveals that there is, indeed, a set of transcripts (Cluster 22; see S2 Data) that is strongly correlated with Lyve1 expression across MPS cells from a larger spectrum of tissues. The cluster includes Mrc1 but excludes Timd4, Cd5l, Fcna, and Vsig4, which were associated with distinct tissue-specific clusters (Table 2). The correlation between Lyve1 and Mrc1 is actually lower than the network threshold of 0.75 (r = 0.62, Fig 2B). The 2 genes were included in Cluster 22 because of shared links to other genes. In fact, Mrc1 was only marginally enriched in the purified LYVE1^hi macrophages from fat, heart, lung, and skin [27]. It was highly expressed in isolated MPS cells from adipose, brain, intestine, kidney, and liver that lack Lyve1 mRNA (see S1 Data/selected transcripts and Fig 2). We conclude that most LYVE1^hi macrophages express Mrc1, but the reciprocal relationship does not hold.

The set of coexpressed genes in Cluster 22 suggests a function for LYVE1^hi macrophages in control of endothelial biology and vascular permeability. It includes genes for 2 of the sphingosine-1-phosphate (S1P) receptors (S1pr1 and S1pr2) that have been implicated in many aspects of inflammation, lymphangiogenesis, and angiogenesis [91, 92]; the vanilloid receptor (Trpv4), which controls capillary barrier function and inflammation [93, 94]; and neuropilin 1 (Nrp1), which controls endothelial homeostasis [95]. Cluster 22 also contains the erythropoietin receptor gene (Epor), which was shown to synergise with S1P to promote apoptotic cell clearance by macrophages [96], and the epidermal growth factor (EGF) receptor gene (Egfr), which has also been shown to regulate macrophage function in a range of inflammatory models [97]. Indeed, the coexpressed genes might support the known functional association of macrophages with lymphatic as well as blood vessels [98]. The Lyve1-associated cluster contains genes for 3 candidate transcriptional regulators: Etv1, Nfatc2, and Tcf4. Etv1 expression in macrophages has been implicated in functional polarisation in vitro and the response to altered mitochondrial membrane potential [99]. NFATC2 is required for osteoclast differentiation in vitro [100], but its roles in macrophage differentiation/function have not been explored. Tcf4 encodes a transcription effector of the Wnt/β-catenin pathway, which is implicated in responses to E-cadherin and other effectors in macrophage differentiation [101].

Mrc1 is commonly referred to as a marker for alternative or M2 macrophage polarisation [31]. Another putative marker of M2 polarisation is the somatic growth factor insulin-like growth factor 1 (Igf1 gene) [102]. Igf1 was correlated with Mrc1 (r = 0.67) but did not form part of a coexpression cluster. It was absent from monocytes and DCs but was highly expressed in most resident tissue macrophages (see S1 Data/selected transcripts). Igf1 is CSF1-inducible and of particular interest because of the profound impact of Csf1r mutations in multiple species on postnatal growth and development [8]. Unlike hepatocytes and mesenchymal cells, tissue macrophages did not express transcripts encoding the growth hormone receptor (Ghr), Igf1r, or the Igf1 binding protein genes (Igfals, Igfbp1, 2, 3, 5, 6). The exception is Igfbp4, which was highly expressed in most macrophage populations and did form part of the Lyve1/Mrc1-associated Cluster 22. Interestingly, Igfbp4 knockout in mice mimics impacts of Csf1r deficiency on somatic growth and adipose formation [103, 104].

The intimate association of macrophages with capillaries was evident from the first localisation of the F4/80 antigen [105]. Adgre1 expression was also correlated with Mrc1 (r = 0.64; Fig 2B), but it was more widely expressed than either Mrc1 or Lyve1 and therefore not within Cluster 22. Adgre1 was not enriched in any of the purified LYVE1^hi macrophage populations relative to LYVE1^lo cells from the same tissue [27]. It was high in most isolated tissue macrophages and induced during differentiation of monocytes in situ, as in the liver data set [32] and the intestinal developmental series [36, 37]. F4/80 (now known as ADGRE1) was proposed as a marker of macrophages of embryonic origin [106], but Adgre1 mRNA was equally high in intestinal macrophages, which turn over rapidly from monocytes [107, 108], and in cDC2s. It was also strongly induced during monocyte differentiation to occupy a vacant KC niche [32]. Whatever the association with ontogeny, the pattern is rodent-specific. Adgre1 is a rapidly evolving gene, and the expression pattern also varies across species [109].

Tissue-specific macrophage clusters

Several coexpressed clusters were associated with MPS cells isolated from a single tissue. Aside from the large brain-enriched expression cluster (Cluster 4) that contains many microglia markers, Cluster 10 was lung-enriched and contains the alveolar macrophage marker Siglecf and key transcription factor Pparg [16]. Cluster 12 was shared amongst liver KCs, peritoneal macrophages, and splenic macrophages and includes the transcription factors Id3, Nr1h3, and Smad6 and markers Cd5l, Clec4f, and Vsig4 [16, 32, 37]. Within Cluster 12, we noted the strong coexpression (r = 0.81) between Nr1h3 and Rxra, the gene encoding its promiscuous heterodimerisation partner, which is also implicated in control of KC lipid and iron metabolism [110] and may have independent function in innate immune regulation [111].

The average expression of Cluster 12 increased progressively in the monocyte-KC differentiation series [32] included in this data set (see profile in S7 Fig). Cluster 12 also reveals the regulated and coordinated expression of the thyroid hormone receptor (Thrb gene), likely mediating the many impacts of thyroid hormones in innate immune function [112]. One other potential regulator identified in this cluster is Zbtb4, which encodes an epigenetic regulator with a high affinity for methylated CpG. Zbtb4^−/− mice are viable and fertile but growth retarded compared with littermates [113]. Impacts on myeloid differentiation have not been reported. The transcription factor SPIC is implicated in splenic red pulp macrophage differentiation and iron homeostasis [114, 115]. Although Spic mRNA was highest in red pulp macrophages, KCs, and bone marrow macrophages, it was detected in other macrophage and DC populations and therefore has a unique expression profile. Cluster 21 contains transcripts most highly expressed in resident peritoneal macrophages and includes the genes for the transcription factor Gata6 and the retinoic acid receptor (Rarb), which control peritoneal macrophage survival and adaptation [116, 117]. The data confirm the specific high expression of the enigmatic plasminogen activator inhibitor encoded by Serpinb2 in resident peritoneal macrophages, first described >20 years ago [118] and still seeking a function [119].

Genes in Cluster 15, including the monocyte-specific chemotactic receptor Ccr2, were highly expressed in classical monocytes. Genes in Cluster 43 were expressed specifically in Langerhans cells (LCs). They include the marker Cd207 (langerin) used in the purification of LCs [62] but also expressed at lower levels in many other tissue macrophage populations. This cluster did not include the gene for another LC marker, Epcam [62]. It was highly expressed in LCs but also detected in one set of intestinal macrophage samples, most likely a contamination with epithelial cells (Cluster 5, see below). Epidermal LCs have at times been considered as DC-like because of their migratory and APC properties but are now considered to be specialised resident tissue macrophages [120]. Unlike most classical DCs in lymphoid tissue but in common with nonlymphoid DCs, they are clearly CSF1R-dependent and share with several other macrophage populations dependence on the conserved enhancer in the Csf1r locus [17]. Cluster 43 did not include a transcriptional regulator specific to LCs. In common with several other macrophage populations, LC differentiation is regulated by transforming growth factor β (TGFβ) signalling, involving transcription factors RUNX3 and ID2 [120]. Both transcription factor genes were highly expressed in LCs but also present in several other tissue macrophage populations.

Intestinal macrophage-enriched gene expression profiles, which have not previously been identified, emerge in Cluster 38. Two large separate data sets of intestinal macrophages were included here [36, 37], both likely reflecting a differentiation series of adaptation from blood monocytes to resident intestinal tissue macrophages [5]. In one case, CD4 and TIM4 were used as markers [37]. Surprisingly, despite the fact that TIM4 was used as a marker to sort intestinal macrophages, the Timd4 gene was not part of a significant coexpression cluster; it was highly expressed in the sorted intestinal macrophages but idiosyncratically in several other tissue macrophage populations. Cd4 mRNA expression was shared uniquely with lung, skin, and kidney macrophage subpopulations (see Fig 2A). A third data set tracks the adaptation of transferred blood monocytes to the intestinal niche [69]. Cluster 38 identified Cxcr4 as a candidate intestinal macrophage marker consistent with their continuous derivation from CXCR4⁺ monocytes. The high expression of Wnt4 in lamina propria macrophages was recently confirmed by immunohistochemistry (IHC). Conditional deletion of Wnt4 using Itgax-cre led to dysregulation of immunity against an intestinal parasite [121]. WNT4 is a candidate mediator of the key trophic role of lamina propria macrophages in the intestinal stem cell niche [122]. Fosb, Hes1, and Hic1 encode identified potential transcriptional regulators of intestinal macrophage differentiation and adaptation. HES1 inhibits inflammatory responses in macrophages and contributes to gut homeostasis [123, 124]. FOSB has not previously been implicated in macrophage adaptation to any niche. Unfortunately, we were not able to include data from a microarray analysis of resident colonic macrophages that identified a set of 108 genes >2-fold higher in the colon relative to other macrophage populations in the ImmGen database [125]. However, Cluster 38 confirmed the gut-macrophage–specific expression of several of these transcripts, including Dna1l3, Fgl2, Gpr31b, Hes1, Mmp13, Ocstamp, Pgf, and Tlr12.

There were no unique expression profiles enriched in macrophages isolated from any other major tissues including adipose, brain (nonmicroglia), heart, kidney, pancreas, or skin. The abundant resident macrophages of adipose are especially topical in light of the obesity epidemic. The literature on adipose macrophages focusses on ‘M2-like’ markers [126]. Amongst resident macrophage populations, Apoe and Retnla, both detected in most tissue macrophages and not included in a cluster, were highest in adipose-derived macrophages. RETNLA (also known as RELMα) has been referred to as an adipokine, regulated by food intake and controlling lipid homeostasis [127]. Kumamoto and colleagues [128] claimed that Retnla was coexpressed with Mgl2 (another putative M2 marker) in many mouse tissues, including adipose, and attributed it a role in maintenance of energy balance. The 2 transcripts were not correlated in this larger data set. In fact, Mgl2 was part of a small cluster (Cluster 83) with Cd14. Like Retnla, mRNA for the related lectin, MGL1 (Clec10a gene)—also considered an M2 macrophage marker [126]—was highest in the adipose-associated macrophages but also expressed in macrophages from other tissues, including dura, heart, lung, and skin (Cluster 101).

DC coexpression clusters

It has become a central dogma in immunology that DCs are uniquely adapted to present antigen to naïve T cells. This view has driven the search for surface markers to enable isolation of DCs for immunotherapy and receptors that can mediate selective antigen uptake to enhance immunisation. If the central dogma is correct and currently used DC markers have any validity, network analysis should uncover coexpression clusters associated with antigen uptake, processing, and presentation to T cells. This prediction is not supported by the data. Despite evidence that it is expressed by many resident tissue macrophages (reviewed in [25]), CD11C (encoded by Itgax) is still widely used as a surface marker in mouse DC purification. Ongoing studies of the impacts of conditional mutations using Itgax-cre continue to be interpreted solely in terms of DC specificity (for example, [121, 129, 130]). Consistent with the literature [25], Itgax was expressed in multiple macrophage populations other than DCs (Fig 2A) at levels at least as high as in DCs purified using CD11C as a marker and correlated only with Cd22, Cd274 (encoding programmed cell death 1 ligand 1), Csf2rb, Csf2rb2, solute carrier (Slc)15a3, Tmem132a, and the transcription factor gene Prdm1 (Cluster145). Class II MHC is also often used as a marker to purify DCs, and expression is obviously a prerequisite for antigen presentation to T cells. The ImmGen Consortium compared DCs from multiple sources with various macrophage populations to identify transcripts that distinguish DC from macrophages [28, 29]. Since the macrophages used for comparison were MHCII^lo, the DC signature included class II MHC genes. In our meta-analysis, one small cluster (Cluster 165) contained the transcription factor gene Ciita and Cd74, H2-Aa, H2Ab1, H2-DMa/b1, and H2-Eb1 encoding its targets. The genes in this cluster were clearly highly expressed in many tissue macrophages (see profile for Cd74 in Fig 2A) but regulated independently of any other markers and expressed no higher in cells annotated as DCs than in cells annotated as macrophages from intestine, lung, heart, and kidney. Interestingly, again highlighting the issue with a definition of DCs based upon unique APC function, isolated lung MHCII^hi interstitial macrophages were as active as cDC2s in antigen-presentation assays in vitro [27]. These results are consistent with the sample-to-sample analyses that showed that DCs and monocytes/macrophages from the same tissue are more similar than DCs are to DCs from other tissues (Figs 3 and 4).

The GCN analysis did identify 3 separate DC-associated coexpression clusters that are consistent with current knowledge of putative DC subsets and adaptation in mice [21, 22, 131]. Cluster 13 includes Ccr7 and transcription factor genes Spib and Stat4; Cluster 28 includes Flt3, Kit, and the transcription factor gene Relb; and Cluster 49 includes cDC1 markers Itgae (CD103) and Xcr1. CCR7 is associated with DC migration [132] and the transcript was abundant in both cDC1s and cDC2s isolated from spleen and LN. By contrast, the expression was much lower in isolated lung DCs and in kidney DCs from a separate data set (see below), similar to levels in isolated macrophages from multiple tissues. Several putative DC markers were excluded from DC-specific clusters. The transcription factor gene Batf3, implicated in cDC1 differentiation [133], did not form part of a cluster and was detected in most macrophage populations (consistent with [16]). Similarly, IRF4 has been attributed a specific function in cDC2 differentiation [130]. Irf4 mRNA was more abundant in cDC2s than in cDC1s but was also expressed in monocytes and monocyte-derived macrophage populations. Transcripts encoding NFIL3 and IRF8, which interact in the regulation of cDC1 differentiation [134], were also highly expressed in cDC2s and in monocytes and many tissue macrophages. Although the transcription repressor gene Zbtb46, encoding a putative DC lineage marker [135], was highest in DCs, it was also detectable in most isolated tissue macrophages, notably in kidney and lung. Another putative DC marker gene, Clec9a [70], also clustered independently because of expression in isolated intestine, kidney, liver, and lung macrophages.

Interestingly, tissue macrophages may contribute to homeostatic regulation of cDC differentiation. The transcript encoding the FLT3 ligand (Flt3l) was expressed constitutively to varying degrees in all of the MPS populations studied. Fujita and colleagues [129] showed that FLT3L is cleaved from the cell surface of expressing cells by ADAM10. Conditional deletion of Adam10 using Itgax-cre led to reduced differentiation of cDC2s. Adam10 is also expressed by CD11C⁺ macrophages; it forms part of Cluster 3, containing genes that are low in monocytes and expressed by all resident macrophages at higher levels than in DCs.

Aside from CLEC9A, many other lectin-like receptors have been proposed as DC markers and inferred to have a function in antigen uptake. Fig 6 shows the profiles for the 12 members of the so-called DC immunoreceptor (DCIR) family. The original member of this family, Clec4a2, the likely ortholog of the single CLEC4A gene in humans, encodes a lectin with a broad binding specificity for mannose and fucose [136]. Studies on knockout mice lacking Clec4a2 continue to be based upon the claim that the lectin is mainly expressed by DCs [137], but the global analysis showed that the mRNA is more highly expressed in most isolated macrophage populations. Two of the DC-associated clusters contained other members of the family, Clec4a4 and Clec4b2. Clec4a4 has been attributed a specific role in cDC1 function [138], but it was equally expressed in cDC2s and forms part of Cluster 28. Most of the Clec4 genes in the mouse genome are in a single location on Chromosome 6. They also include macrophage-inducible C-type lectin (Mincle) encoded by Clec4e, which mediates innate immune responses to Candida [139]. The related Clec4f (KC marker) and Cd207 (langerin) are located together in a separate locus on Chromosome 6. Each of the Clec4 genes had a unique expression profile in tissue MPS populations. Analysis of the entire data set reveals that ‘DCIR’ is a misnomer for this family. The DC designation has also been misapplied to other surface receptors, including DC-SIGN (CD209 in humans), DEC205 (Ly75), and DC-HIL (Gpnmb). In mice, there are multiple Cd209 paralogs. Cd209b was highly expressed in marginal zone macrophage populations in spleen and is Csf1r-dependent [89]. These cells have not been successfully isolated by tissue disaggregation. Four members of the CD209 family (Cd209a, d, f, g) were coexpressed in a unique pattern (Cluster 100) together with Cbr2, Ccl24, and Clec10a. Ly75 was detected in both cDC subpopulations but was most highly expressed in lung macrophages (Cluster 10).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Expression of members of the DCIR (CLEC4) family across MPS cell populations.

(A) Expression patterns across cells from different tissues. Each column represents a sample. Upper bar along the X axis shows the cell type (black—monocytes and macrophages; red—DCs). Lower bar shows the tissue, coloured as shown in the key. Y axis shows expression level in TPM, calculated using Kallisto. (B) Correlations (Pearson correlation coefficient) between expression patterns of different Clec4 genes. DC, dendritic cell; DCIR, DC immunoreceptor; MPS, mononuclear phagocyte system; TPM, transcripts per million.

https://doi.org/10.1371/journal.pbio.3000859.g006

CD64 (Fcgr1 gene) was used as an exclusion criterion to remove or separate macrophages from DCs or to enrich macrophages in all of the data sets included herein based upon the earlier studies of the ImmGen Consortium [28]. This exclusion was clearly successful in that all the purified DC data sets have very low Fcgr1 (Fig 2A and S1 Data), but the expression of this gene in macrophage populations was also highly variable. As a simple screen for additional markers that distinguish all CD64⁺ ‘macrophages’ from all CD64⁻ ‘DCs’, we averaged expression across all macrophage and DC samples and compared them (see S1 Data). Amongst the transcripts that were robustly expressed and highly enriched in macrophages to at least the same extent as Fcgr1, those encoding surface markers were also variably expressed amongst macrophage populations. However, we identified 3 transcription factor genes—Cebpb, Mafb, and Klf10—that were apparently excluded from all of the CD64⁻ cDCs. The role of Cebpb in macrophage differentiation is well-recognised [66, 140, 141], and one of the data sets includes progenitors from Cebpb^−/− mice [66]. There is evidence of a negative feedback relationship with Irf8 in monocyte-derived DCs [142]. Cebpb was detected in most tissue macrophages but uniquely excluded from some populations, notably the heart and intestinal muscularis. MAFB has been proposed previously as a lineage marker separating macrophages from DCs [143, 144]. However, the most parsimonious explanation of the data would be that MAFB regulates expression of Fcgr1. The literature on KLF10 is more limited, with evidence that it participates in TGFβ-induced macrophage differentiation [145].

In overview, although our analysis identifies coregulated clusters associated with cells currently defined as DCs, most proposed markers of this population are clearly shared with other MPS cells, and there is no enrichment for any genes that could confer unique APC activity.

Resident macrophage activation during isolation

Cluster 41 contains numerous immediate early genes (IEGs) encoding transcription factors and feedback regulators (for example, Fos, Egr1, and Dusp1), consistent with evidence from scRNA-seq of disaggregated cells that isolation of cells from tissues produces cell activation [59, 146]. In many samples, IEGs were amongst the most highly expressed transcripts. The majority of isolated MPS populations also had high levels of macrophage-specific lipopolysaccharide (LPS)-inducible genes. Cluster 224 contains Ccl2, Ccl7, Ccl12, Cxcl1, and Il6; Cluster 329 includes Il1b and Ptgs2 (also known as Cox2); and Cluster 485 contains Tnf and inducible chemokines Ccl3 and Ccl4. The anti-inflammatory cytokine Il10, which is also LPS-inducible, formed part of the intestinal macrophage cluster (Cluster 38). IL10 is essential to intestinal homeostasis [107], but Il10 mRNA was detected in only 1 of the 3 intestinal macrophage data sets [37] alongside very high expression of IEGs and proinflammatory cytokine genes (for example, Il1b, Tnf). The apparent expression of Fosb in intestinal macrophages discussed above is likely also an artefact because it is undetectable in total intestinal mRNA (see http://biogps.org). Inflammation-associated transcripts were highlighted as evidence of activation in vivo in sensory neuron-associated macrophages [61]. Similarly, Chakarov and colleagues [27] highlighted selective expression of Il10 in interstitial lung macrophages and differential expression in the LYVE1^hi/MHCII^lo subpopulation. They did not comment upon the reciprocal pattern observed with Tnf and Il1b, which were both more highly expressed in the LYVE1^lo macrophages. Both populations of interstitial lung macrophages (and all the samples from other tissues in this BioProject) expressed very high levels of all of the IEG transcripts in Cluster 41. Whereas macrophage-expressed transcripts such as Adgre1 are readily detected in total tissue mRNA and are CSF1R-dependent, inflammatory cytokines and IEG transcripts are not [9, 48]. Accordingly, in each of the RNA-seq data sets we have analysed, the expression of IEGs and inducible cytokines is most likely an artefact of tissue disaggregation and/or cell isolation and separation by FACS.

Interestingly, Acod1, which was massively induced within 1 hour by LPS in mouse macrophages in vitro (see http://biogps.org), was only detected at low levels in a small subset of samples and was not correlated with IEGs or any other inflammatory activation markers. Induction of this gene has been attributed functions in adaptive immunometabolism and accumulation of tricarboxylic acid (TCA) cycle intermediates in activated macrophages [147]. The lack of detection in the isolated macrophages suggests either that induction is specific to recruited inflammatory macrophages or that inducible expression is purely an in vitro phenomenon. The Acod1 expression pattern was correlated only with Il23a (encoding a subunit of the cytokine IL23) at the stringency used here (r ≥ 0.75).

Contamination of macrophage populations with other cell types

Table 3 and Fig 5B highlight other clusters that were tissue-specific and contained markers and transcription factors associated with organ/tissue-specific differentiation, with corresponding enrichment for GO terms associated with specific tissues (S2 Data). There are 3 ways in which mRNA from purified macrophage/DC populations may be contaminated with mRNA from unrelated cells. The most straightforward is poor separation of macrophages from unrelated contaminating cells by FACS for purely technical reasons. A second source derives from active phagocytosis by macrophages of dying (senescent/apoptotic) cells, in which RNA from the engulfed cell may be detected. Finally, there is a phenomenon that arises from the extensive ramification of macrophages and their intimate interactions with other cells. Gray and colleagues [148] found that cells purified from LNs with the surface marker CD169 were in fact lymphocytes coated with blebs of macrophage membrane and cytoplasm. Similarly, Lynch and colleagues [149] found that all methods to isolate KCs for flow cytometry produced significant contamination, with CD31⁺ endothelium tightly adhered to remnants of KC membrane.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Major contaminant clusters.

https://doi.org/10.1371/journal.pbio.3000859.t003

Cluster 2 appears to be a generic ‘rubbish’ cluster, containing transcripts detected at relatively low levels only in specific BioProjects and unrelated to tissue of origin. Other clusters were driven by a single RNA-seq result from within 1 BioProject. These clusters most likely represent technical noise as well as contamination.

Consistent with the proposal from Lynch and colleagues [149], 3 endothelial-associated transcripts—Cdh5, Pecam1, and Stab2—were contained with the KC-enriched cluster (Cluster 12) and apparently increased in expression during KC differentiation. However, other abundant endothelial transcripts were absent. Bonnardel and colleagues [78] generated RNA-seq data from purified liver sinusoidal ECs. We examined the profiles of the most highly expressed EC genes in the macrophage data set. Many of them were detectable in isolated KCs but at much lower levels than Cdh5, Pecam1, and Stab2. They contributed to a separate liver-specific endothelial cluster (Cluster 76). So, whilst there is evidence that ECs contaminate KC preparations, reflecting the close apposition in the sinusoids, Chd5, Pecam1, and Stab2 are likely also genuine KC-expressed transcripts.

The detection of mature red cell transcripts encoding haemoglobins (Hba, Hbb), which are quite abundant in many macrophage populations, most likely reflects ongoing erythrophagocytosis. Macrophages isolated from the intestinal lamina propria in 1 of the 2 large data sets from small intestine [36] were heavily contaminated with markers of intestinal epithelium (Clusters 5 and 17). This might be a separation artefact but could also reflect an active role of macrophages in homeostatic turnover of epithelial cells [150]. Cluster 18 and Cluster 45 were restricted to samples from a study of pancreatic islet and peri-islet macrophage populations [60]. The authors noted the expression of insulin (Ins1) mRNA in their islet macrophage populations and attributed it to an intimate interaction with β-cells. Contamination or β-cell–macrophage fusion was said to be excluded on the basis that β-cell markers such as Pdx1 were not detected. However, many other islet-associated transcripts were abundant and formed part of Cluster 18, notably transcription factors Isl1, Foxa2, Nkx6.1, and Nkx2.2as well as other islet-specific transcripts, Inhba, Chga/b, Iapp, Gipr, and Gcg. Similarly, Cluster 45 was relatively enriched in the peri-islet macrophages and contains transcripts encoding many pancreatic enzymes. Cluster 65 includes Acta2 and other smooth muscle markers that selectively contaminated macrophages isolated from the intestinal muscularis [36].

The bone marrow contains several populations of macrophages [33], including those associated with haematopoietic islands expressing CD169 (Siglec1 gene) and VCAM1. One of the data sets included in the present study profiled the transcriptome of macrophages associated with erythroblastic islands, based upon isolation using an Epor-EGFP reporter gene [57]. A second bone marrow data set separated macrophages based upon their engagement in phagocytosis of blood-borne material [65]. The putative erythroblastic island macrophages did not actually express increased Epor mRNA (although Epor was detected in other macrophage populations as reported recently [96] and fell within Cluster 22). However, in the isolated bone marrow macrophages, Siglec1 was correlated with high levels of both immature neutrophil (Cluster 32) and erythroid-associated (Cluster 33) mRNAs. The separation of these 2 clusters implies that the contamination occurs in distinct macrophage populations, enriched selectively in each preparation and perhaps derived from separate haematopoietic islands [33]. Cluster 32 also contains the myeloid progenitor transcription factor Myb and the GM-progenitor marker Ms4a3. Given the extensive ramification of marrow macrophages and their intimate interactions with progenitors [33], this contamination likely reflects the same isolation artefact reported in LN [148], namely haemopoietic progenitor cells cloaked in macrophage clothing.

There are separate clusters including B-cell and natural killer (NK)-cell–specific markers. The B-cell cluster, Cluster 87, shows the highest average expression in intestine, bone marrow, lung, and spleen and likely reflects close association between macrophages and B cells in lamina propria and germinal centres [36]. The cluster containing NK cell markers, Cluster 67, had the highest average expression in one of the DC preparations. Those DCs came from a study that proposed a further subdivision of cDC2s based upon expression of transcription factors T-bet (Tbx21) and RORγT [67] and separated cDC2s based upon expression of a Tbx21 reporter allele. Tbx21 was detected in all purified splenic cDC preparations presented on http://biogps.org but at much lower levels than in NK cells. NK cells also express Itgax, used in purification of the cDCs. Accordingly, it seems likely that apparent Tbx21 expression in DCs is due to NK cell contamination.

Clustering of transcription factor expression

Most of the coregulated clusters identified above contain genes encoding transcriptional regulators that are known to be essential for tissue-specific adaptation. These represent only a small subset of the transcription factors detected in MPS cells. The r value of 0.75 was selected for the analysis of the whole data set to maximise the number of genes included whilst minimising the number of edges between them (S6 Fig) and aimed at assessing the predictive value of markers including those shown in Fig 2. To test the effect of reducing the stringency, we focussed on annotated transcription factors [151] to reduce the complexity and remove noise. One thousand, one hundred and three transcriptional regulators were detected above the 10 TPM threshold in at least 1 MPS population. The sample-to-sample matrix including all samples formed a single network, including the annotated DCs, as shown in Fig 7. As with the whole data set, increasing the r value resulted in separation of various cell populations, but the DCs remained in a group with monocytes/macrophages. We generated GCNs at 3 different Pearson correlation coefficient thresholds, 0.5, 0.6, and 0.7 (S8 Fig). The results are provided in S3 Data. As the cutoff was reduced, more transcription factor transcripts were included in the network. At the highest stringency r value (≥0.7), the largest cluster includes Spi1 alongside many of the transcription factors identified in the largest generic MPS coexpression clusters above (Clusters 1, 3, and 4). We conclude that the basic shared identity of MPS cells involves coordinated expression of around 100–150 transcription factors. Even at the lowest r value (≥0.5), transcription factor genes identified as specific to particular tissue-specific MPS populations made few additional connections, indicating that local adaptation is dependent on highly correlated and regulated expression of a small cohort of transcription factors. Nevertheless, associations that become evident at lower r values may identify combinatorial interactions in particular cell populations: Mycl, associated with DC fitness [152], was weakly correlated with Irf8 and Zbtb46; Cebpb with Nfil3; and interferon-related transcription factors (Batf2, Irf1/7/9, Stat1/2) were connected at the threshold of 0.5 (S3 Data).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Network analysis of transcription factor gene expression in MPS cell populations.

The sample-to-sample network was generated by BioLayout analysis at r ≥ 0.66, which included all 466 samples. Nodes representing samples are coloured by source tissue (left) and cell type (right). Lists of genes and expression profiles of clusters at different r values are presented in S3 Data. MPS, mononuclear phagocyte system.

https://doi.org/10.1371/journal.pbio.3000859.g007

Expression of solute carriers and metabolism genes in MPS cell populations

The burgeoning field of immunometabolism has focussed on regulation of intermediary metabolism in recruited monocytes and macrophages in various states of activation or polarisation [147]. Amongst emerging concepts is the view that M1 polarisation (classical activation) is associated with aerobic glycolysis and mitochondrial dysfunction, whereas M2 polarisation requires an active tricarboxylic acid cycle [147]. We used the MPS transcriptome to infer likely metabolic adaptations of tissue-resident MPS cells.

Cluster 7 contains mitochondria-associated transcripts and transcripts encoding ribosomal subunits, with variable expression across all samples even from the same tissue, indicating that resident tissue macrophages vary in their dependence upon mitochondrial oxidative phosphorylation irrespective of surface markers or differentiation state.

In many cases, metabolic pathways are regulated at the level of solute transport [147]. There were > 400 members of the SLC family expressed in mononuclear phagocytes above the 10 TPM threshold. Some were more highly expressed in intestine and kidney epithelial cells and clustered with tissue-specific epithelial markers. However, many contributed to macrophage-enriched expression clusters. One such gene, Slco2b1, which encodes an organic anion transporter of unknown function, has been proposed as a marker gene to distinguish macrophages from DC subsets, and the promoter was used in an inducible macrophage depletion strategy [27]. The larger data set herein does not support this dichotomy. Slco2b1 is part of Cluster 4, enriched in microglia and absent from multiple other macrophage populations, in addition to both cDC subsets.

Macrophages depend to varying degrees upon glutamine, glucose, and fatty acids as fuels [153], and glutamine is an important immune regulator [154]. 14 different solute carriers from 4 families have been shown to transport glutamine [155]. Of the genes encoding these carriers, Slc38a1 was widely expressed in MPS cells and did not fall within a cluster, whereas Slc7a5, Slc7a7, Slc7a8, and Slc38a7 were part of distinct macrophage-enriched clusters. Consistent with the importance of glutamine as a fuel for resident MPS cells, transcripts encoding enzymes of glutamine metabolism (Gls, Glud1, Glul, Slc25a11) were also highly expressed and part of Clusters 1 and 3. By contrast, resident MPS cells apparently have very limited expression of glucose transporters. Slc2a1 (encoding glucose transporter GLUT1) was low, highly variable, and idiosyncratic amongst tissues. A myeloid-specific conditional knockout of Slc2a1 confirmed that GLUT1 was the major glucose transporter in macrophages analysed in vitro, but the loss of glucose as a fuel had remarkably little impact on macrophage function [156]. The expression of Slc2a1 in cells isolated from tissues is difficult to interpret because the transporter is induced by hypoxia [157], which might arise during isolation.

In the absence of Slc2a1, macrophages in vitro increased oxidation of fatty acids [156]. The Slc27a1 gene, encoding the fatty acid transporter FATP1—which also contributes to functional regulation in macrophages [158, 159]—was widely expressed in tissue macrophages and, with carnitine acyl transferase genes (Crat, Crot), formed part of Cluster 1. Slc2a5 (found in Cluster 4) encodes a fructose-specific transporter [160] and was expressed primarily in microglia. Slc2a6 is a lysosome-associated glucose transporter that was recently knocked out in the mouse genome [161]. It also has a novel expression profile, being highest in monocytes and cDC2s.

One of the best known functional solute carriers in macrophages is natural resistance associated membrane protein 1 (NRAMP1; Slc11a1 gene), which is associated with genetic resistance to intracellular pathogens. SLC11A1/NRAMP1 is expressed in lysosomes and contributes to pathogen resistance by restricting available iron [162]. The role in iron metabolism is reflected by its presence in Cluster 12 alongside Slc40a1, encoding ferriportin, the macrophage-enriched iron exporter [163]. One other prominent class of solute carriers highly expressed in macrophages (Slc30a6, Slc30a7, Slc30a9, Slc39a3, Slc39a7, and Slc39a9 in Cluster 1 and Slc39a12 in Cluster 4) is involved in transport of zinc, which is a component of antimicrobial defence [164, 165]. Two further zinc transporters, Slc39a2 and Slc39a11, were enriched in lung macrophages (Cluster 10). This lung-macrophage–enriched cluster also contains Slc52a3, encoding a riboflavin transporter, Slc6a4 (sodium- and chloride-dependent sodium symporter), and 2 members of the Slc9 family of sodium–hydrogen exchange (NHE) transporters (Slc9a4 and Slc9a7), which are more traditionally associated with epithelial function [166].

Validation of coexpression clustering with an independent kidney data set

The abundant macrophage populations of the kidney were first described in detail using F4/80 as a marker in situ [167]. There has been considerable debate about the relationships between resident macrophages, monocyte-derived macrophages, and cDC subsets in the kidney [168]. The main cluster analysis did not reveal a separate kidney-resident macrophage-enriched profile. The kidney data set in the preceding analysis included F4/80⁺, CD64⁺ macrophages isolated from control and ischaemic kidneys, further subdivided based upon expression of CD11B and CD11C [54]. Salei and colleagues [70] recently produced RNA-seq data for isolated populations of resident macrophages, monocyte-derived cells, cDC1s, and cDC2s from kidney compared with similar populations from spleen. The primary data were not available for download by our automated pipeline through the ENA at the time we pooled and froze our data set (February 2020). We therefore obtained the processed data directly from the authors and carried out network analysis using the 33 samples and 9,795 genes with normalised expression of at least 10 in at least 1 sample. This analysis served to validate the approach to analysis and the core conclusions using an independent data set.

The macrophages of the kidney are intimately associated with the capillaries [167, 169], but Lyve1 was not detectable in resident macrophages in this data set or in [54]. Published IHC on mouse kidney reveals that LYVE1 is restricted to lymphatic vessels [170]. Fig 8 illustrates the way in which the sample-to-sample matrix revealed relationships between the cell populations with increasing correlation coefficient threshold. Even at the lowest correlation cutoff used in the main atlas (0.75), the splenic red pulp macrophages separated from all kidney and DC samples. As the cutoff was made more stringent, the cDC1s from both spleen and kidney separated, but the resident kidney macrophages, cDC2s, and monocyte-derived macrophages remained closely connected until r ≥ 0.98, when the spleen cDC2s separated from the monocyte-derived macrophages and kidney cDC2s. At r ≥ 0.99, the kidney cDC2s and monocyte-derived macrophages were still not separated, indicating that the expression profiles of these cell types are very similar. Salei and colleagues [70] performed a principal components analysis based upon the 500 most variable transcripts and also identified the close relationship between cDC2s and monocyte-derived cells. Our analysis further emphasises their conclusion that the main axis of difference is between spleen macrophages and all other cells. cDC1s from both tissues were more similar to each other than to the other cells, but spleen cDC2s were only separated from kidney cDC2s and monocyte-derived macrophages at the highest stringencies. We also performed a gene-to-gene analysis on these data. The profiles of kidney myeloid cells other than cDC1s were very similar and differed by only a small number of genes. Consistent with this conclusion, the 2 largest clusters in this analysis (see S4 Data) were shared between all of the isolated populations and contain Spi1 as well as many of the DC-enriched markers identified in the main analysis. However, Ccr7 and many of the genes associated with it in the main data set (Cluster 13, Table 2; for example, Spib, Stat4, Vsig10, Cd200, and Itgb8) were expressed at low levels in isolated kidney DCs as in lung DCs. Cluster 3 of the kidney analysis was specific to splenic red pulp macrophages and contains the known transcriptional regulators Pparg, Spic, and Nr1h3. Transcripts in Cluster 4 were enriched in the resident kidney macrophages compared to both splenic macrophages and other kidney myeloid populations. Interestingly, the resident kidney macrophage cluster includes many genes that are also highly expressed in microglia and depleted in the brain in Csf1r mutant mice and rats, including Cx3cr1, C1qa/b/c, Csf3r, Ctss, Fcrls, Hexb, Laptm5, Tgfbr1, Tmem119, and Trem2 [17, 89]. These were also detected in the isolated kidney macrophages in S1 Data. Both microglia and resident F4/80^hi kidney macrophages are selectively lost in a mouse line with a mutation in a conserved enhancer of the Csf1r locus [17]. RUNX1 regulates the activity of the Csf1r enhancer [171] and has also been implicated in the establishment of microglial cells during development [172]. The Runx1 gene was within this cluster. Csf1r mRNA was expressed at high levels in cells defined as cDC2s, as well as monocyte-derived cells and resident macrophages. All cells expressing a Csf1r-EGFP reporter in the kidney were depleted by treatment with anti-CSF1R antibody [34]. This suggests that despite their expression of FLT3, renal cDC2s are CSF1R-dependent, a conclusion consistent with previous evidence that cells classified as cDC2s in other nonlymphoid tissues are dependent upon CSF1R rather than FLT3 [173]. Cluster 6 of the kidney analysis, including Itgam, was enriched as expected in the selected CD11B⁺ populations from kidney but highly expressed in all of the populations. This cluster includes all of the co-regulated IEGs identified in Cluster 41 in the extended MPS data set above, suggesting that recent monocyte-derived cells may be more susceptible to activation during isolation. In summary, sample-to-sample and gene-to-gene networks on this smaller independent data set are entirely consistent with conclusions from the global MPS analysis that question the basis for the separation of DC from macrophages using surface markers.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Sample-to-sample 2D network analysis of gene expression in macrophage and DC subpopulations from kidney and spleen.

The sample-to-sample network was generated by BioLayout analysis at the indicated Pearson r values, which all included all 33 samples up to r ≥ 0.98. Above r ≥ 0.98, 1 red pulp macrophage sample was lost. Red, kidney DC1; pink, kidney DC2; dark blue, spleen DC1; light blue, spleen DC2; black, kidney macrophages; grey, monocyte-derived macrophages; white, spleen red pulp macrophages. Lists of genes and expression profiles of clusters are presented in S4 Data. DC, dendritic cell.

https://doi.org/10.1371/journal.pbio.3000859.g008

The relationship between single-cell and bulk RNA-seq data

The advent of scRNA-seq has been heralded as a revolution, promising new approaches to classification of myeloid heterogeneity [26, 174, 175]. To determine whether scRNA-seq produces distinct insights about MPS heterogeneity, we wished to assess the relationship between the coexpression networks we have generated from large total RNA-seq data sets and the coexpression modules inferred from scRNA-seq. scRNA-seq is an intrinsically noisy, nonquantitative stochastic sampling of a subset of the most abundant mRNAs in individual cells [176, 177]. Algorithms that support nonlinear dimensional reduction (for example, t-distributed stochastic neighbour embedding [t-SNE] or Uniform Manifold Approximation and Projection [UMAP]) [178] followed by some form of clustering are then used to join together groups of cells in which the members share detectable expression of an arbitrarily defined set of markers. The number of populations defined depends upon the parameters applied, and different approaches do not always give the same answers [178]. There is an implicit assumption in this approach that defined cell types with approximately identical transcriptomes actually exist and that sampling noise can be overcome by analysis of a sufficiently large number of cells. Based upon scRNA-seq analysis of interstitial lung macrophages, Chakarov and colleagues [27] inferred the existence of a subpopulation that selectively expressed Lyve1. They then generated bulk RNA-seq data from separated LYVE1^hi and LYVE1^lo subpopulations. Their data uniquely support a critical comparison of the 2 approaches to transcriptome analysis and the outcomes of our network analysis. For this purpose, the primary scRNA-seq data were downloaded, reanalysed, and expressed as TPM using the Kallisto pseudoaligner, as described in the ‘Materials and methods’. S5 Data contains these reprocessed scRNA-seq data, alongside the bulk RNA-seq data for the lung macrophage subpopulations from the same study, with the level of expression ranked based upon the bulk RNA-seq data for the purified LYVE1^hi interstitial macrophages.

Consistent with Zipf’s law, the power-law distribution of transcript abundance [71, 72], the top 200 expressed transcripts in the bulk RNA-seq data (approximately 1% of the total) contribute around 50% of the total detected transcripts in the scRNA-seq data (S5 Data). The abundant transcripts from bulk RNA-seq that were also detected in scRNA-seq samples include many cell-type–specific surface markers, which explains the ability to use scRNA-seq to discover such markers. These abundant transcripts also include IEGs such as Dusp1, Egr1, Fos, Ier2, and Junb, indicative of the activation that occurs during isolation as discussed above. The inducible cytokines, including Ccl2, Tnf, Il1b, Il6, and Il10, were each detected in a subset of the single cells, most likely also induced during isolation. Of the 200 most highly expressed transcripts identified in the bulk RNA-seq data, only a very small subset (including Actb, Apoe, B2m, Ccl6, Cd74, Ctsb, Fth1, Ftl1, and Lyz2) had nonzero values in all cells in the scRNA-seq output. The average expression of the top 500 transcripts in the single cells was similar to the bulk RNA-seq, but the detected expression level varied over 4 orders of magnitude amongst individual cells. Fcgr1 and Mertk mRNAs, encoding markers used to purify the interstitial macrophages for scRNA-seq, as well as other commonly used markers (Cx3cr1, Itgax) were actually detected in only a small subset of the cells and were not correlated with each other. Both this study and a subsequent study [76] state that Mrc1 and Lyve1 expression is shared by overlapping populations of lung interstitial macrophages. That conclusion is not supported by the data. Even in the bulk RNA-seq data from lung interstitial macrophages, the expression of Mrc1 was only marginally enriched in purified LYVE1^hi cells relative to LYVE1^lo cells (S1 Data). The separation of these 2 markers was evident from the separate study of lung interstitial macrophage populations [52] included in our global MPS analysis and has been discussed above. Consistent with that conclusion, in the scRNA-seq data, the two are not strictly correlated with each other, with Mrc1 being detected in many more cells than Lyve1 (Fig 9A) despite similar absolute levels of expression in the total RNA-seq data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Network analysis of scRNA-seq data.

The reprocessed data and the bulk RNA-seq data for the lung macrophage subpopulations from the same study are available in S5 Data. (A) Expression profiles in single cells for selected genes. Each column represents RNA from a single cell. Y axis shows expression in TPM, calculated using Kallisto. Only the first 6 cells expressed Lyve1 (coloured red). (B) The sample-to-sample network was generated by BioLayout analysis at r ≥ 0.53, which included all 54 single-cell samples. Nodes represent samples; red nodes show the samples with high expression of Lyve1. (C) Gene-to-gene network (r ≥ 0.5), clustered at MCL inflation value of 1.7. Cluster lists and expression profiles are available in S6 Data. MCL, Markov clustering algorithm; RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA-seq; TPM, transcripts per million.

https://doi.org/10.1371/journal.pbio.3000859.g009

To identify whether any robust correlations actually exist in the scRNA-seq data, the top 500 expressed transcripts in the scRNA-seq samples were used for network analysis. The sample-to-sample network (r ≥ 0.53) is shown in Fig 9B, and the gene-to-gene network (r ≥ 0.5) in Fig 9C. The cluster list and average expression profiles are provided in S6 Data. One clear-cut finding is the coexpression of genes involved in APC activity (H2-Aa, H2-Ab1, H2-Eb1, Cd74, and Ctss; Cluster 13 of the scRNA-seq analysis, indicated in Fig 9C), which were effectively present or absent in individual cells. Chakarov and colleagues [27] defined 2 subpopulations as LYVE1^hi/MHCII^lo and LYVE1^lo/MHCII^hi, but only 6 of the scRNA-seq samples expressed Lyve1, and only half of those also expressed detectable MHC II genes (Fig 9A). This is consistent with the lack of any inverse correlation between Lyve1 and Cd74 in Fig 2B. Even at this low r value, known highly expressed markers segregated from each other. Lyve1 forms a cluster with Mgl2, Cd209, and Cd302 (Cluster 7 in the scRNA-seq analysis; Fig 9C). Adgre1 is in a coexpression cluster that includes Lyz2 and Msr1 (Cluster 4 of the scRNA-seq analysis), Csf1r is coexpressed with Mrc1 and Cd163 (Cluster 2), and Lgals3 with Retnla and Fcrls (Cluster 1). The coregulation of MHC-related genes and genes located in the same chromosomal region (for example, C1qa, C1qb, C1qc; Cluster 25 of the scRNA-seq analysis), as well as the relatively uniform detection of genes such as Actb (Fig 9A), suggests that a significant proportion of the all-or-nothing differences in expression between cells in the scRNA-seq data is real. That conclusion may be considered a reflection of the limitations of the technology, but it is actually supported by other evidence. Tan and Krasnow [179] defined subpopulations of interstitial lung macrophages based upon expression of F4/80, Mac-2 (Lgals3), and Class II MHC and tracked the changes in their relative abundance during development. Interestingly, they did not detect LYVE1 on adult lung interstitial macrophages by IHC. Consistent with their data, in the scRNA-seq data examined here, most lung interstitial macrophages expressed high levels of either Adgre1 or Lgals3, but some expressed both or neither.

Across our data set, the selection of MPS populations based upon surface markers was also strongly enriched for the mRNA encoding those markers (for example, CD11C⁺ cells express Itgax). The reciprocal relationship need not be the case. Protein expression at a single-cell level clearly does not vary to the same extent as mRNA because proteins have different rates of translation, turnover, and decay [180]. Markers such as F4/80 and CD11C and transgenes based upon macrophage-enriched promoters such as those of Csf1r and Cx3cr1 do appear to label the large majority of MPS cells in most tissues. The disconnect between scRNA-seq and cell surface markers may partly reflect the nature of transcription. At the single-cell level, transcription occurs in pulses interspersed by periods of inactivity and mRNA decay, which can manifest as random monoallelic transcript expression [181]. If gene expression is genuinely probabilistic at the level of individual loci [180], the assumption of transcriptomic homogeneity in definable cell types upon which scRNA-seq analysis is based is clearly invalid. The number of macrophage subpopulations that can be defined in any scRNA-seq data set becomes a matter of choice and model. As an extreme example, one recent scRNA-seq study identified 25 distinct myeloid cell differentiation ‘states’ in a mouse lung cancer model [182].

A critical view of the validity of markers

The RNA-seq data included as representative of cDC subsets [67] were from cells purified using CD64 as a marker to exclude macrophages. Despite this choice, an unbiased assessment of the sample-to-sample networks in Figs 3B and 4 (based on all genes) and Fig 7 (based on transcription factor genes) would class all of these DCs as part of the same family as monocytes and macrophages. The use of CD64 as a definitive marker distinguishing macrophages from DCs was criticised when it was proposed [44], and it remains untenable. It is actually a curious choice as a marker to define a cell as a macrophage because the encoded protein FCGR1 (CD64) has been implicated functionally in APC activity [183]. Interestingly, a recent study of inflammation in the lung now posits the existence of CD64⁺⁺ cDC2s and suggests that CD26 (encoded by Dpp4) is a more definitive marker [184]. In the separate kidney data set [119], monocyte-derived macrophages were separated from cDC2s based upon CD64 as a marker. Indeed, the cDC2s lacked Fcgr1 expression, but there was no enrichment for Dpp4. As noted above, despite the use of CD64 as a marker, the transcriptomes of the 2 populations (cDC2s and monocyte-derived macrophages) were almost indistinguishable. Like other MPS cells, DC expression profiles may also be influenced by tissue-specific signals. Both the larger data set and the kidney data [70] suggest that there is tissue-specific adaptation of ‘cDC2s’ that may remain more ‘macrophage-like’, which may be a reflection of their dependence upon CSF1R. From our analysis, the proposed separation of DCs from all other members of the MPS based upon APC function, surface markers, transcription factors, or ontogeny [21] remains problematic [19].

The concept of M1/M2 polarisation derives from analysis of classical and alternative activation of recruited monocytes by Th1 (interferon gamma [IFNγ]) and Th2 (IL4/IL13) cytokines [31]. Previous meta-analysis indicated that proposed M2 markers defined by others [31] correlate poorly with each other in isolated inflammatory macrophages and are not conserved across species [30]. The M1/M2 concept was also challenged in a recent comparative analysis of in vitro and in vivo data on macrophage gene expression [185], which concluded that ‘valid in vivo M1/M2 surface markers remain to be discovered’. We would suggest that they do not exist; each proposed marker also has its own unique transcriptional regulation and tissue-specific function in resident macrophages. Aside from proposed M2 markers already mentioned that each have idiosyncratic expression (Mrc1, Retnla, Igf1, Mgl2), Chil3 (also known as Ym1) was highly expressed in lung macrophages (Cluster 10 of the whole data set analysis; S2 Data), Arg1 and Alox15 were restricted to peritoneal macrophages (Cluster 21), and Cd163 was part of a small cluster of 4 transcripts (Cluster 312). Detection of M2 markers on resident macrophages cannot imply that they share any functions with alternatively activated recruited monocytes. Nevertheless, IL4/IL13/STAT6 signalling could contribute to resident MPS cell differentiation. The IL13 receptor (Il13ra) is part of the generic MPS Cluster 1, and Il4ra is also highly and widely expressed. IL4 administration to mice can drive resident tissue macrophage proliferation beyond levels controlled homeostatically by CSF1 [186].

How do transcriptional networks contribute to understanding macrophage heterogeneity in situ?

One concern with analysis of the cells isolated by tissue digestion and evaluated here is whether recovered cells are representative of the tissue populations. Several macrophage populations have resisted isolation, notably those of the marginal zone of spleen and subcapsular sinus of LNs and the abundant macrophage populations detected in skeletal muscle [5], whilst others are clearly fragmented during isolation as discussed above. Given the remarkable ramification of resident MPS cells in situ, it seems unlikely that they are quantitatively recovered intact in any isolation protocol. Subpopulations of isolated cells defined by markers can sometimes be linked to precise location within the tissue. One recent example is the clear separation of distinct myeloid populations in the liver. A unique subcapsular MPS cell population distinct from KC was initially classified as DC-like based upon expression of CD11C and apparent lack of macrophage markers [187], but these cells were subsequently characterised as F4/80⁺, CSF1R-dependent macrophages uniquely expressing CD207 [188]. Another is the apparent location of LYVE1^hi macrophages with capillaries in the lung [27]. On the other hand, it is unclear where the putative long-lived CD4⁺, TIM4⁺ population in the gut [37] is located. In broad overview, macrophages in every organ, detected with Csf1r reporter transgenes that are expressed in all myeloid cells including DCs, have a remarkably regular and uniform spatial distribution. The concept of a macrophage territory [5] or a niche [189, 190] has been proposed. Despite this apparent homogeneity in location and morphology, multicolour immunolocalisation of macrophage surface markers indicates that they are almost infinitely heterogeneous (reviewed in [24]). Some of this variation may be purely stochastic [180]. However, most of the data sets analysed here suggest that monocytes and macrophages in each organ are a differentiation series. We have taken the view that macrophages in tissues have a defined half-life such that some cells survive by chance and continue to change their gene expression [5]. Each macrophage that occupies a new territory, following either infiltration as a monocyte or self-renewal by cell division, starts a life history that involves changes in gene expression and surface markers with time. In that view, many MPS subpopulations may be no more than arbitrary windows within a temporal profile of adaptation.

Conclusion

The transcriptional network analysis confirms that using our unique approach to downsizing and a common quantification pathway, the RNA-seq data from different laboratories can be merged to provide novel insights. The network analysis indicates the power of large data sets to detect sets of coregulated transcripts that define metastable states of MPS adaptation and function. The merged data set we have created provides a resource for the study of MPS biology that extends and complements resources such as ImmGen (http://www.immgen.org). It can be readily expanded to include any new RNA-seq data for comparative analysis, and clusters can be regenerated using BioLayout or the further development of this network approach, Graphia, which are freely available (http://biolayout.org; https://graphia.app). For example, 2 recent studies have profiled isolated macrophage subpopulations from peripheral neurons [191, 192].

Clusters of transcripts that are robustly correlated give clear indications of shared functions and transcriptional regulation. However, our analysis also revealed 2 important artefacts in the study of isolated tissue macrophages: the clear evidence of inflammatory activation during isolation and the extensive contamination of isolated preparations with transcripts derived from other cell types. One recent innovation that may obviate these issues is the use of a so-called MacTrap transgene to isolated macrophage-specific polysome-associate mRNA from tissues without isolating the cells [193].

A discussion review of MPS heterogeneity in 2010 [194] suggested that in order for the field of immunology to advance and communicate, ‘all cells have to be called something’. This Linnaean view continues to drive efforts to classify MPS cells into subsets based upon markers. The analysis we have presented shows that surface markers are poorly associated with each other and have very limited predictive value. Aside from Class II MHC, there are no markers that can be correlated with predicted APC activity. Resident tissue MPS cells, including cells that are currently defined as DCs, belong to a closely related family of cells in which the transcriptomic similarities are much greater than the differences. The cumulative function of the population of MPS cells acting together within each tissue is likely to be more important to homeostasis and immunity than the individual heterogeneity.