Network analysis of transcriptomic diversity amongst resident tissue macrophages and dendritic cells in the mouse mononuclear phagocyte system.

The mononuclear phagocyte system (MPS) is a family of cells including progenitors, circulating blood monocytes, resident tissue macrophages, and dendritic cells (DCs) present in every tissue in the body. To test the relationships between markers and transcriptomic diversity in the MPS, we collected from National Center for Biotechnology Information Gene Expression Omnibus (NCBI-GEO) a total of 466 quality RNA sequencing (RNA-seq) data sets generated from mouse MPS cells isolated from bone marrow, blood, and multiple tissues. The primary data were randomly downsized to a depth of 10 million reads and requantified. The resulting data set was clustered using the network analysis tool BioLayout. A sample-to-sample matrix revealed that MPS populations could be separated based upon tissue of origin. Cells identified as classical DC subsets, cDC1s and cDC2s, and lacking Fcgr1 (encoding the protein CD64) were contained within the MPS cluster, no more distinct than other MPS cells. A gene-to-gene correlation matrix identified large generic coexpression clusters associated with MPS maturation and innate immune function. Smaller coexpression gene clusters, including the transcription factors that drive them, showed higher expression within defined isolated cells, including monocytes, macrophages, and DCs isolated from specific tissues. They include a cluster containing Lyve1 that implies a function in endothelial cell (EC) homeostasis, a cluster of transcripts enriched in intestinal macrophages, and a generic lymphoid tissue cDC cluster associated with Ccr7. However, transcripts encoding Adgre1, Itgax, Itgam, Clec9a, Cd163, Mertk, Mrc1, Retnla, and H2-a/e (encoding class II major histocompatibility complex [MHC] proteins) and many other proposed macrophage subset and DC lineage markers each had idiosyncratic expression profiles. Coexpression of immediate early genes (for example, Egr1, Fos, Dusp1) and inflammatory cytokines and chemokines (tumour necrosis factor [Tnf], Il1b, Ccl3/4) indicated that all tissue disaggregation and separation protocols activate MPS cells. Tissue-specific expression clusters indicated that all cell isolation procedures also co-purify other unrelated cell types that may interact with MPS cells in vivo. Comparative analysis of RNA-seq and single-cell RNA-seq (scRNA-seq) data from the same lung cell populations indicated that MPS heterogeneity implied by global cluster analysis may be even greater at a single-cell level. This analysis highlights the power of large data sets to identify the diversity of MPS cellular phenotypes and the limited predictive value of surface markers to define lineages, functions, or subpopulations.

Introduction

The mononuclear phagocyte system (MPS) [1] is a family of cells including progenitors, circulating blood monocytes, resident tissue macrophages, and dendritic cells (DCs) that are present in every tissue in the body [2-5]. Within each tissue, resident macrophages occupy niches or territories with a remarkably regular distribution (reviewed in [5, 6]). The proliferation, differentiation, and survival of most resident macrophage populations depends upon signals from the macrophage-colony–stimulating factor receptor (CSF1R) initiated by one of 2 ligands, CSF1 or IL34 [7, 8]. Based upon detection of macrophage-restricted mRNA, including Csf1r, the relative abundance of resident macrophages in most organs in mice was shown to reach a maximum in the first week of postnatal life and remains stable thereafter during postnatal growth [9]. Lineage-trace studies in the C57BL/6 strain suggest that many macrophage populations established in the mouse embryo are maintained in adults mainly by self-renewal, whereas others are replaced progressively to differing extents by blood monocytes derived from bone marrow progenitors throughout life [10-12]. Most, if not all, tissue macrophage populations can be generated and maintained in the absence of blood monocytes because of the intrinsic homeostatic regulation by circulating CSF1 [13]. The precise details of ontogeny, turnover, and homeostasis of resident macrophages may not be conserved across mouse strains or species [5]. However, regardless of their steady-state turnover, all tissue-resident macrophages, including the microglia of the brain, can also be rapidly replaced by blood monocytes following experimental depletion ([3–6, 13] and references therein).

Within individual tissues, resident macrophages acquire specific adaptations and gene expression profiles [2, 4, 5, 14–16]. These adaptations contribute to survival as well as function and involve inducible expression of transcription factors and their downstream target genes. At least some of these transcription factors act by regulating Csf1r expression. Deletion of a conserved enhancer in the mouse Csf1r gene leads to selective loss of some tissue macrophage populations, whereas others express Csf1r normally and are unaffected [17]. In the mouse embryo, where abundant macrophage populations are engaged with phagocytosis of apoptotic cells [18], the macrophage transcriptome does not differ greatly between organs. Tissue-specific macrophage adaptation occurs mainly in the postnatal period, as the organs themselves exit the proliferative phase and start to acquire adult function [9, 16].

Classical DCs (cDCs) are commonly defined functionally on the basis of a proposed unique ability to present antigen to naïve T cells, a concept that requires a clear distinction between DCs and macrophages [19]. The situation is confused by the widespread use of the term DC to describe any antigen-presenting cell (APC), including cells that are clearly derived from blood monocytes [20]. An attempt at consensus proposed an MPS nomenclature classification based upon ontogeny and secondarily upon location, function, and phenotype [21]. The proposal separates monocyte-derived APCs from cDC subsets: cDC1s, dependent on the transcription factor BATF3, and cDC2s, dependent upon IRF4. Some support for this separation came from analysis of an Ms4a3 reporter transgene, which labelled cells derived from committed granulocyte-macrophage (GM) progenitors and distinguished monocyte-derived cells from tissue DCs [11]. Secondary classification is based upon cell surface markers that are presumed to be linked in some way to ontogeny. The proposed development pathway of these DC subsets from a common myeloid progenitor via a common DC progenitor (CDP) has been reviewed recently [22]. However, it remains unclear as to whether cDCs should be considered part of the MPS and the extent to which they can be separated from other MPS cells based upon surface markers [13].

Even within individual tissues, resident macrophages are extremely heterogeneous [23, 24]. Since the advent of monoclonal antibodies and later development of transgenic reporter genes [25], numerous markers have been identified that segregate the MPS into subpopulations. More recently, mouse tissue macrophage heterogeneity has been analysed using multiparameter flow cytometry and single-cell RNA sequencing (scRNA-seq) [26]. Amongst the recent suggestions, LYVE1 was proposed as a marker of macrophages associated with the vasculature [27], CD64 (Fcgr1 gene) and MERTK as markers that distinguish macrophages from cDCs [28, 29], and CD206 (Mrc1 gene) as a marker of so-called M2 macrophage polarisation [30, 31]. Several surface markers have also been identified that are encoded by genes expressed only in macrophages in specific tissues (for example, Clec4f, Tmem119, Siglecf) [16, 32]. Other markers define macrophages in specific locations within a tissue, for example, CD169 (encoded by Siglec1) in the marginal zone of spleen and haematopoietic islands in bone marrow [33]. In the case of blood monocytes, the subpopulations are clearly a differentiation series in which short-lived LY6Chi ‘classical’ monocytes give rise in a CSF1R-dependent manner [34] to long-lived LY6Clo nonclassical monocytes via an intermediate state [12, 34, 35]. This is likely also the case in tissues such as the liver [32] and intestine [36, 37].

Mechanistically, the association between marker expression and cellular function depends upon coordinated transcriptional regulation. One way to identify coregulated sets of transcripts is to cluster large transcriptomic data sets. This approach was used to create transcriptional atlases in multiple species and identify lineage-specific transcription factors and their target genes [38-42]. It enabled the extraction of a generic tumour-associated macrophage signature from multiple large cancer data sets [43]. Previous meta-analysis of large microarray data sets [38, 39, 42], as well as a reanalysis of data from the ImmGen Consortium [44], indicated a clear separation in the mouse of MPS cells from all other leukocyte lineages but did not support the basic premise that markers can separate macrophages from DCs or define lineages within the MPS.

Over the past 5 years, RNA-seq has supplanted microarrays as an approach to expression profiling. The recent cascade of interest in tissue-specific macrophage adaptation has produced RNA-seq data for MPS cells isolated from most major organs of C57BL/6 mice. There has been no previous effort to integrate this data deluge into a cohesive view of MPS transcriptional diversity and to identify sets of transcripts that are stringently co-ordinately regulated. To enable comparative analysis of data sets from multiple laboratories, we devised an automated informatics pipeline employing random sampling of RNA-seq data to a common depth and quantification using the pseudoaligner Kallisto. Robust transcriptional atlases for the chicken [45] and pig [46] were generated using data sets from numerous divergent sources. The analysis of these merged data sets, as well as large multitissue data sets from sheep, human, and water buffalo [47-49], provided strong support for the principle of guilt by association, namely that genes that contribute to a specific biological function tend to be coregulated. Each of these analyses identified transcripts that were enriched specifically in MPS cells relative to other haematopoietic cells. Using the same basic pipeline as in the chicken and pig projects, we identified a total of 466 RNA-seq libraries generated from isolated macrophage and cDC populations from 24 different studies that sample mouse MPS transcriptional diversity (Table 1). Here, we apply network clustering to this large data set to identify shared and divergent transcriptional adaptation of tissue-resident MPS cells and revisit the relationships between macrophages and DCs.

Table 1

GEO and BioProject accession numbers for samples used in the analysis.

SRA and NCBI accessions and sample descriptions are available in S1 Data.

Accession	BioProject	Reference	Description (Markers Used in FACS Purification)
GSE125691	PRJNA517169	[27]	Interstitial subsets from lung, skin, fat, heart + monocytes, and alveolar macrophages (LYVE1, SIGLECF).
GSE84586	PRJNA330530	[51]	Resident macrophages from heart, kidney, and liver (F4/80, CD11B).
GSE94135	PRJNA369038	[52]	Three interstitial subsets from lung (MERTK, CD64, CD11B, CD11C, CD206, MHCII) + alveolar macrophages.
GSE95859	PRJNA378611	[53]	Brown adipose macrophages (CX3CR1-EGFP).
GSE114434	PRJNA471340	[37]	Monocytes and small intestinal macrophage subsets (CD4, TIM4, CD64).
GSE116094	PRJNA478258	[54]	Kidney-resident and monocyte-derived subpopulations, effect of ischaemia (F4/80, CD64, CD11B, CD11C, MHCII).
GSE122766	PRJNA506249	[55]	Brain microglia, bone marrow-derived brain macrophages (CD45, CD11B, CX3CR1).
GSE123021	PRJNA507265	[56]	Brain microglia, cortex, cerebellum, hippocampus, striatum (TMEM119).
GSE127980	PRJNA525977	[57]	Erythroblastic island macrophages from marrow (EPOR-EGFP, F4/80, VCAM1, SIGLEC1).
GSE135018	PRJNA557178	[58]	Alveolar macrophages and peritoneal macrophages, effect of Bhlhe40/41 mutation (SIGLECF, CD11B, CD11C, F4/80).
GSE128662	PRJNA528430	[32]	Monocyte to KC differentiation series. Effects of Nr1h3 and Smad4 mutations (F4/80, CD11B, LY6C, CLEC4F).
GSE128781	PRJNA529096	[59]	Nonparenchymal brain macrophages, microglia, and peritoneal macrophages (MHCII, CD64, CD11B).
E-MTAB-6977	PRJEB27719	[36]	Macrophage subsets from intestinal lamina propria, serosa, and muscularis (CD64, CX3CR1 lineage trace).
GSE112002	PRJNA438927	[60]	Pancreatic islet and peri-islet macrophage populations. Effect of high-fat diet (F4/80, CD11B, CD11C).
GSE103847	PRJNA407286	[61]	White adipose and sympathetic neuron-associated macrophages, spleen, microglia (CD45, CX3CR1-EGFP, F4/80).
GSE68789	PRJNA283850	[62]	Mucosal and skin LCs and DCs (CD103, CD11B, EPCAM, CD207).
GSE128518	PRJNA527979	[63]	White adipose macrophages, effect of Trem2 mutation (CD11B, F4/80).
GSE107130	PRJNA419127	[64]	Brain microglia developmental time course: male and female. Role of microbiome (CD45, CD11B, F4/80, CD64).
GSE83222	PRJNA325288	[65]	Spleen, intestine, bone marrow macrophages. Effect of engulfment of apoptotic cells (F4/80, CD11B).
GSE95702	PRJNA378162	[66]	Monocyte subsets and bone marrow progenitors. Effect of Cebpb mutation (CD115, CD135, LY6C, CD11B, CD11C).
GSE130201	PRJNA534273	[67]	DCs, LN, and spleen. cDC1s/cDC2s (CD11C, CD64, MHCII, CD103, TBX21).
GSE120012	PRJNA491337	[68]	Cardiac vessel macrophages (MHCII, CCR2, CD64, CD11B).
GSE140919	PRJNA519465	[69]	Monocyte engraftment of colon/ileum (CX3CR1-EGFP, CD115, LY6C, CD64).
GE131751	PRJNA 544681	[70]	Kidney-resident and monocyte-derived macrophage and DCs (F4/80, CD64, CD11B, CD11C, MHCII, CLEC9A lineage trace).

Abbreviations: cDC, classical DC; DC, dendritic cell; EGFP, enhanced green fluorescent protein; FACS, fluorescence activated cell sorting; GEO, Gene Expression Omnibus; KC, Kupffer cell; LC, Langherhans cell; LN, lymph node; NCBI, National Center for Biotechnology Information; SRA, Sequence Read Archive.

Materials and methods

The RNA-seq data sets from within the BioProjects shown in Table 1 were downloaded from the European Nucleotide Archive (ENA). S1 Data contains all the Sequence Read Archive (SRA) and National Center for Biotechnology Information (NCBI) accessions and sample descriptions. Individual BioProjects differ in methods of mRNA isolation, library preparation and sequencing methods, length, depth, and strandedness, but previous analysis in other species [45, 46] indicated that they can still produce comparable expression level estimates. Prior to expression quantification and for the purpose of minimising variation between samples, all libraries were randomly downsampled to 10 million reads, 5 times each, as described previously [45, 46]. The expression levels were then requantified using Kallisto v0.44.0 [50], and the expression level was taken as the median transcripts per million (TPM) across the 5 downsampled replicates. Kallisto quantifies expression at the transcript level as TPM by building an index of k-mers from a set of reference transcripts and then ‘pseudo-aligning’ reads to it, matching k-mers in the reads to k-mers in the index. We built a custom index (k = 21) containing the combined set of 154,627 unique protein-coding transcripts from Ensembl and NCBI RefSeq, representing 24,149 protein-coding genes (Mus musculus annotation GRCm38.p6). Because expression is quantified relative to this index, Kallisto is robust to the presence of spurious k-mers in the reads so that ranked TPM estimates are largely unaffected by fastq preprocessing. For paired-end samples, Kallisto estimates the fragment length from the reads. For single-end samples, fragment length cannot be empirically derived from read mapping and is assumed to follow a truncated Gaussian distribution with user-specified mean and standard deviation. For the single-end libraries, we considered the mean fragment length to be 1.2 × the median read length and the standard deviation to be 0.1 × the mean fragment length. Varying these parameters did not substantially alter the expression profile of each sample.

The selected BioProjects include subsets of resident tissue macrophages defined using surface markers or reporter genes as indicated in Table 1 and separated by FACS, as well as temporal profiles of adaptation from monocytes to tissue macrophages. Several studies involve the analysis of the impact of mutations in specific transcription factors and surface receptors. The focus is on tissue-resident cells. Data sets related to inflammatory macrophages or macrophages stimulated in vitro have been excluded.

The purpose of this analysis was to identify clusters of transcripts that are robustly correlated regardless of the tissue of origin. The outcomes of such an analysis may reveal tissue-specific modules but may equally include modules that are shared by several tissues or specific niches within tissues. For this purpose, the size of the data set and the diversity of transcriptomic space sampled is a major strength.

In an RNA-seq library, the distribution of TPM estimates should comply, to a reasonable approximation, with Zipf’s law, which states that the probability of an observation is inversely proportional to its rank [71, 72]. We confirmed that each of the 466 libraries obeyed the predicted power-law relationship.

Prior to network analysis, transcripts that were not detected at an arbitrary threshold of 10 TPM in at least 1 sample were removed to further minimise stochastic sampling noise intrinsic in RNA-seq data. Given the nature of the samples, this also helps to reduce the low-level representation of transcripts derived from contaminating cells of nonmyeloid origin. Of the 18,175 genes that met this minimum threshold, 11,578 were detected in at least 90% of the RNA-seq data sets, and 6,901 had a median expression >10 TPM across the 466 samples. The TPM estimates for the 18,175 genes quantified in all 466 samples included are provided in S1 Data.

Network analysis was performed using the program BioLayout (http://biolayout.org). Pairwise Pearson correlations (r) were calculated between all samples to produce a sample-to-sample correlation matrix and inversely between all pairs of genes to produce a gene-to-gene correlation matrix. Gene coexpression networks (GCNs) were generated from the matrix, in which nodes represent genes and edges represent correlations between nodes above a defined correlation threshold. Note that BioLayout supports 3-dimensional visualisation of network graphs, and each of the graphs generated in this study can be regenerated using the data provided in Supporting Tables and the freely available software. For each gene-to-gene analysis, the value of r was chosen to generate an optimal network that retains the maximum number of transcripts (nodes) with the minimum number of edges [45]. Networks based on Spearman correlation coefficients (nonparametric) were also calculated to assess the impact of distribution shape on the analysis.

Results and discussion

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.

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.

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.

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.

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).

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.

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).

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.

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.

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.

Table 2

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

Cluster Number	Description	Representative Genes
1	MPS	Acp2, Atp6 subunits, Cd276, Cd53, Cd68, Cd84, Clec5a, Cln5/8, Csf1r, Ddx/Dhx family, Fcgr1, Gpr107/108, Hk3, lysosomal enzymes, Ifngr1/2, Il10ra, Il13ra1, Il6ra, Irak1/2, Jak1/3, Lamp1/2, Lgals8/9, M6pr, P2ry6, P2rx7, Sirpa, Tlr6/7/8, Tnfrsf11a, Cebpg, Creb3, Crebzf, Elf1, Etv5, Fli1, Foxj2, Foxn3, Foxo1, Gabpa, Hdac3/10, Hif1a, Hsf1, Klf3, Maf1, Mafg, Mitf, Nfatc1, Nfx1, Nfyc, Nr1h2, Nr2c1, Nr2f6, Nr3c1, Prdm4, Rela, Smad1/2/4, Sp3, Spi1, Srebf1, Stat6, Tcf3, Tfe3
3	MPS	Abca1/2, Aim2, Akt2/3, Arrb1, Arrb2, Atxn7, Bak1, Cbl, Cd180, Cdk8/10/12/13/19, Csk, Ddi2, Ddx3/6/17/19a/21/23/39b/46, Dhx9/15, Grk2, H6pd, Ly9, Megf8, Mertka, Mpeg1, Naip2/5/6, Nirp1b, Ptprj, Socs4/7, Syk, Taok1/2, Traf7, Tram2, Atf1, Bach1, Bcor, Cebpa, Elf2/4, Erf, Foxk1, Foxk2, Foxo3, Foxo4, Fus, Hsf2, Ikzf1, Maf, Maz, Mef2d, Ncoa3, Ncoa6, Nfat5, Nfatc3, Nfya, Pbx2, Prdm2, Smad5
4	Microglia and brain macrophages	Abi3, Acvr1, Adrb2, Bcl9, Bmp1/2k, Card6, Ccr5, Cd34, Csf3r, Cx3cr1, Cxxc5, Ddx31/43, Entpd1, Fcrls, Fgf13, Gabbr1, Gpr155, Gpr165, Gpr34, Hexb, Itgb3/b5, Lpcat1/2/3, Mrc2, Nckap5l, Olfml3, P2ry12/13, Paqr7, Plexna4, Nanos1, Siglech, Slc1a3/4, Slco2b1, Slc2a5, Sipa1, Tgfbr1, Tmem119, Tmem173, Trem2, Vav1, Vsir, Bhlhb9, Ebf3, Elk3, Ets1, Hivep3, Lefty1, Mef2c, Prox2, Sall1/2/3, Sox4
7	Mitochondria and ribosome	Atp5e/g2/h/j2/l, Cox5b/6a1/6b1, Mrpl family, Nduf family, Rpl and Rps families
9	Cell cycle	Aurka, Aurkb, Birc5, Bub1, Ccna2/b1/b2/e2, Cdk1, Cenpe, Haus family, Kif family, Mcm family, Plk1, Foxm1, Mybl2
10	Lung macrophages	Anxa2, Atxn10, Car4, Cd2, Cd200r4, Cd9, Chil3, Ctsk, Cx3Cl1, Cxcr1, F7, Fabp1, Ffar4, Flt1, Flvcr2, Gal, Htr2c, Igflr1, Il1rn, Lpl, Ly75, Nceh1, P2rx5, Plscr1, Serpine1, Siglecf, Slc6a4, Tmem138, Nlrx1, Pparg, Tcf7l2
12	Liver KCs, peritoneal and splenic red pulp macrophages	Acp5, Adgre4, Apoc1, C6, Cd5l, Cdh5, Clec1b, Clec4f, Fabp7, Fcgr4, Il18bp, Itga9, Kcna2, Lrp5, Ly9, Pecam1, Pira1/2, Ptger1, Ptprj, Scarb1, Scarf1, Sema6d, Siglec1, Siglece, Slc11a1, Slc40a1, Slc1a2, Stab2, Tmem65, Treml4, Trpm2, Vsig4, Elk1, Id3, Nr1h3, Rxra, Smad6, Thrb, Zbtb4
13	CCR7 DCs	Arc, Birc2, Cacnb3, Cblb, Ccl19, Ccl22, Ccr7, Cd1d1, Cd200, Cd40, Cd70, Dpp4, Fas, Icosl, Glipr2, Gpr68, Heatr9, H2-Q6/7/8/9, Il15, Il15ra, Itgb8, Laptm4b, Lrrk1, Slamf1, Socs2, Tank, Tmem19, Tnfrsf4, Tnfrsf9, Traf1, Tyk2, Vsig10, Zc3h12c, Zmynd15, Foxh1, Id2, Ikzf4, Spib, Stat4
15	Monocytes	C3, Camkk2, Ccr2, Cd177, Cd244a, Celsr3, Clec2g, Erbb4, Fgr, Gpr15, Gpr35, Gpr141, Hpse, Il17ra, Itga4, Met, Mmp8, Ms4a4c, Nlrc5, Ptgir, Ptprc, Sell, Sgms2, Slk, Vcan, E2f2, Foxn2, Jarid2, Rara, Rfx2, Stat2b
21	Peritoneal macrophages	Ackr3, Alox15, Arg1, C4a/b, Car6, Cyp26a1, F5, F10, Fgfr1, Fzd1, Icam2, Itga6, Itgam, Jag1, Lbp, Lrg1, Mst1r, Naip1, Nt5e, Padi4, Pycard, Selp, SerpinB2, Slpi, Tgfb2, Thbs1, Wnt2, Gata6, Rarb, Smad3, Sox7, Tox2
22	LYVE1-positive macrophages	Adam9, C3ar1, C5ar1, Cd36, Cfh, Clcn5, Ctsb, Dab2, Egfr, Epor, F13a1, Fcgrt, Frmd6, Gas6, Gpr160, Igfbp4, Lyve1, Mrc1, Nrp1, S1pr1/2, Tlr5, Tmem9, Trf, Trpv4, Etv1, Nfatc2, Tcf4
28	DCs	Adam11, Bcl2a1b/d, Ccr6, Cd7, Clec4a4, Ddr1, Dtx1, Flt3, H2-DMb2, H2-Eb2, H2-Oa/b, Kit, Lta/b, Nlrp10, P2ry10, Siglecg, Sirpb1a, Tnfrsf18, Relb
38	Intestinal macrophages	Adam19, Asb2, Cxcl9, Cxcr4, Dna1l3, Fgl2, Gpr31b, Gpr55, Il10, Il12rb1, Kynu, Mmp9/13/14, Ocstamp, P2rx6, Pgf, Tlr12, Wnt4, Fosb, Hes1, Hic1
41	Immediate early genes	Ccrl2, Dusp1, Mcl1, Tnfaip3, Trib1, Zfp36, Atf3, Egr1, Fos, Ier2/5, Jun, Junb, Jund, Klf2, Klf6, Nfe2l2, Nfkbiz, Tgif1
43	LCs	Cd207, Dkk1, Dpep3, Hapin3, Il1r2, Mfge8, P2rx2, P2rx5, Plek2, Sema7a, Serpind1, Tnfaip2
49	cDC1s	Cd8a, Clec4b2, Clnk, Ctla4, Gcsam, Gpr33, Gpr141b, Gpr171, Ildr1, Itgae, Il12b, P2ry14, Procr, Plekha5, Tlr11, Xcr1, Ncoa7
165	Class II MHC	Cd74, H2-Aa, H2Ab1, H2-DMa/b1, H2-Eb1

Bold text indicates transcription factors. For descriptions and accession numbers of all genes, please see S1 Data.

Abbreviations: cDC, classical DC; DC, dendritic cell; KC, Kupffer cell; LC, Langerhans cell; MHC, major histocompatibility complex; MPS, mononuclear phagocyte system.

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 LYVE1hi macrophages relative to LYVE1lo 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 LYVE1hi 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 LYVE1hi macrophages express Mrc1, but the reciprocal relationship does not hold.

The set of coexpressed genes in Cluster 22 suggests a function for LYVE1hi 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 LYVE1hi macrophage populations relative to LYVE1lo 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 MHCIIlo, 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 MHCIIhi 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).

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.

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 LYVE1hi/MHCIIlo subpopulation. They did not comment upon the reciprocal pattern observed with Tnf and Il1b, which were both more highly expressed in the LYVE1lo 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.

Table 3

Major contaminant clusters.

Cluster Number	Description	Representative Genes
2	General neuronal contamination	Cacna family, Cdh family, Chrn family, Gabrg1/g2, glutamate receptors, etc.
5 and 17	Intestinal epithelial	Multiple SLCs, Cdx1/2, Hox family, Isx, Ihh
8 and 14	Kidney epithelia	Pax8, Cldn4/8, Hnf1b, Hoxb2/7
16	Hepatic parenchymal cells	Alb, C8/9, Cyp2 family, Igfbp1, Serpina1, Nr1l3
18	Pancreatic islets	Ins1, Gcg, Isl1
26	Skin/keratinocytes	Krt4/5/6, Stfn, Pitx1/2
32	Bone-marrow–specific, neutrophil contamination	Elane, Camp, Fcer1a, Gpc1, M6s4a3, Mpo, Prg2/3, S100a8/S100a9, Gata2, Gfi1, Cebpe, Myb
33	Immature erythroid	Hemgn, Klf1
36	Neuronal	Tnfrsf14, Pax6, Sox8
45	Pancreatic acinar cells	Cel, Cpa1, Ctrb1, Pnlip
65	Smooth muscle (intestine muscularis	Acta2, Cnn1, Des, Mylk, Tpm1, Nkx3-2
67	NK cells	Cd3g, Cd160, Gzma/b/c, Il2rb, Itga2, Kirg1, Klra4/7/8/9, Klrc2/3, Ncr1
76	Endothelial	Adgrf5, Clec4g, Ehd3, Flt4, Kdr, Ptprb, Robo4, Tie1, Sox18, Gata4
87	B cells	Blk, Cd19, Cd79a, Cxcr5, Fcer2a, Fcmr, Itk, Lax1, Tnfrsf13c, Mef2b

Abbreviations: NK cell, natural killer cell; SLC, solute carrier.

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).

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.

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/80hi 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.

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.

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 LYVE1hi and LYVE1lo 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 LYVE1hi 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 LYVE1hi cells relative to LYVE1lo 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.

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.

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 LYVE1hi/MHCIIlo and LYVE1lo/MHCIIhi, 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 LYVE1hi 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.

Spearman correlation coefficients.

(A) Spearman correlation coefficients for expression patterns of different housekeeping genes compared with the Pearson correlation coefficients from Fig 1. (B) Spearman correlation coefficients for expression patterns of different MPS genes compared with the Pearson correlation coefficients from Fig 2. MPS, mononuclear phagocyte system.

(PDF)

Click here for additional data file.

Each sphere (node) represents a sample, and lines between them (edges) show Pearson correlations between them of ≥0.85. The network includes 458 samples. (A) Samples coloured by tissue of origin. (B) Samples coloured by cell type. (C) Samples coloured by BioProject. DC, dendritic cell.

(PDF)

Click here for additional data file.

Each sphere (node) represents a sample, and lines between them (edges) show Pearson correlations between them of ≥0.95. The network includes 418 samples. (A) Samples coloured by tissue of origin. (B) Samples coloured by cell type. (C) Samples coloured by BioProject. DC, dendritic cell.

(PDF)

Click here for additional data file.

Each sphere (node) represents a sample, and lines between them (edges) show Spearman correlations between them of ≥0.85. The network includes 443 samples. (A) Samples coloured by tissue of origin. (B) Samples coloured by cell type. (C) Samples coloured by BioProject. DC, dendritic cell.

(PDF)

Click here for additional data file.

Each sphere (node) represents a sample, and lines between them (edges) show Spearman correlations between them of ≥0.9. The network includes 427 samples. (A) Samples coloured by tissue of origin. (B) Samples coloured by cell type. (C) Samples coloured by BioProject. DC, dendritic cell.

(PDF)

Click here for additional data file.

Graph size compared with correlation threshold for the analysis of the mouse macrophage data set.

The chosen correlation threshold of 0.75 resulted in inclusion of 12,775 nodes, making 1,113,125 edges (correlations of ≥0.75) between them.

(PDF)

Click here for additional data file.

Average expression of genes in Cluster 12 during differentiation of monocytes to KCs.

Data from BioProject PRJNA528435. Clec4f-cre Rosa26iDTX mice were treated with DTX to remove mature KCs. Livers were harvested at indicated time points after DTX treatment. Control animals were treated with PBS and harvested at 72 hours. The experiment shows the repopulation of the liver with cells derived from monocytes. DTX, diphtheria toxin; KC, Kupffer cell.

(PDF)

Click here for additional data file.

Graph size compared with correlation threshold for the analysis of the mouse macrophage transcription factor data set.

Red line shows the highest threshold to include all 1,103 nodes (r ≥ 0.28). Black broken lines show the 3 correlation thresholds used in the analysis: r ≥ 0.5 (1,064 nodes), r ≥ 0.6 (949 nodes), and r ≥ 0.7 (714 nodes).

(PDF)

Click here for additional data file.

Excel spreadsheet containing gene expression data for all MPS samples expressed as TPM.

Separate sheet highlights genes of interest encoding surface markers and transcription factors. Analysis includes means, standard deviation, CoV, and Mac:DC expression ratios. CoV, coefficient of variance; DC, dendritic cell; Mac, macrophage; MPS, mononuclear phagocyte system.

(XLSX)

Click here for additional data file.

Cluster lists for the gene-centred network analysis of the complete MPS data set including graphs of average expression profiles.

Separate sheet shows the GO term enrichment scores for each cluster. GO, gene ontology; MPS, mononuclear phagocyte system.

(XLSX)

Click here for additional data file.

Clusters lists for gene-centred network analysis of transcripts encoding transcription factors at 3 r values: 0.5, 0.6, and 0.7.

(XLSX)

Click here for additional data file.

Cluster lists for the expression data for kidney and spleen MPS populations from [119].

MPS, mononuclear phagocyte system.

(XLSX)

Click here for additional data file.

Excel spreadsheet containing expression data for scRNA-seq and total RNA-seq from lung MPS populations from [27].

MPS, mononuclear phagocyte system; RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA-seq.

(XLSX)

Click here for additional data file.

Cluster lists for the GCN analysis of lung MPS scRNA-seq data from [27] including graphs of average expression profiles.

GCN, gene coexpression network; MPS, mononuclear phagocyte system; RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA-seq.

(XLSX)

Click here for additional data file.