Network analysis of large-scale ImmGen and Tabula Muris datasets highlights metabolic diversity of tissue mononuclear phagocytes

Abstract

The diversity of mononuclear phagocyte (MNP) subpopulations across tissues is one of the key physiological characteristics of the immune system. Here, we focus on understanding the metabolic variability of MNPs through metabolic network analysis applied to three large-scale transcriptional datasets: we introduce (1) an ImmGen MNP open-source dataset of 337 samples across 26 tissues; (2) a myeloid subset of ImmGen Phase I dataset (202 MNP samples); and (3) a myeloid mouse single-cell RNA sequencing (scRNA-seq) dataset (51,364 cells) assembled based on Tabula Muris Senis. To analyze such large-scale datasets, we develop a network-based computational approach, genes and metabolites (GAM) clustering, for unbiased identification of the key metabolic subnetworks based on transcriptional profiles. We define 9 metabolic subnetworks that encapsulate the metabolic differences within MNP from 38 different tissues. Obtained modules reveal that cholesterol synthesis appears particularly active within the migratory dendritic cells, while glutathione synthesis is essential for cysteinyl leukotriene production by peritoneal and lung macrophages.

Keywords: CP: Immunology; CP: Metabolism; ImmGen; immunometabolism; mononuclear phagocytes; myeloid cells; network analysis; single-cell RNA-seq.

Figure 1.. General overview of ImmGen mononuclear phagocytes open-source (IG MNP OS), ImmGen mononuclear phagocytes phase 1 (IG MNP P1), and myeloid Tabula Muris Senis (mTMS) datasets

(A) Schematic representation of Mus musculus tissues, where samples were derived from (marked with colored dots depending on the dataset). (B) Number of tissues overlapping across all datasets. (C–G) Cell-type distribution across all datasets. Principal-component analysis (PCA) based on 12,000 most expressed genes across all samples colored by the tissue of its origin (E and G) or cell type (D and F). (H and I) Uniform manifold approximation and projection (UMAP) representation of cells colored by the tissue of its origin (I) or its type (H). LN, lymph node.

Figure 2.. TMS single-cell RNA-seq dataset

(A) Dataset preprocessing resulting in myeloid subset derivation. (B–D) UMAP plot with natural clusters (B) and cell types (C) identified based on cell specific markers (D). NP, neutrophil; Mo, monocyte; prog, progenitor; DC, dendritic cell; MF, macrophage; alvMF, alveolar macrophage; MG, microglia; KC, Kupffer cell; pDC, plasmacytoid dendritic cell; migDC, migratory dendritic cell.

Figure 3.. Scheme of analysis approach for multi-sample metabolic network clustering (GAM clustering)

The dataset’s metabolic genes are initially clustered based on a k-medoids algorithm. Averaged gene expression of the obtained clusters is further considered as patterns. For each gene, a score is calculated on the basis of its correlation with each pattern. These scores are superimposed on the KEGG metabolic network. Based on these scores, the most weighted connected subnetwork is found for each parent. After the refinement procedure, metabolic modules as a final version of subnetworks are obtained.

Figure 4.. Metabolic modules as a result of multi-sample metabolic network clustering of all myeloid cells but not inflammatory conditions from ImmGen MNP OS dataset

(A) Heatmap representing samples hierarchically clustered based on averaged gene expression of each of obtained module (from lowest as blue to highest as red). Euclidean distance is used as a clustering metric. YS MF, yolk sac macrophage; EB MF, embryoid body macrophage; alvMF, alveolar macrophage; SPM, small peritoneal macrophage; MG, microglia; MF, macrophage; Mo, monocyte; DC, dendritic cell; pDC, plasmacytoid DC; migDC, migratory DC. (B) Annotation of the obtained modules based on gene enrichment in KEGG and Reactome canonical pathways. Enrichment value is calculated as a percentage of module genes contained in a particular pathway. (C) Radar chart representation of metabolic modules within each metasample. Each individual sample is shown as a gray line, while mean of all samples inside one metasample is shown as a colored line. Nine radii of the radar chart are devoted to the corresponding metabolic modules: 1 and 2: lipid metabolism; 3: FAS pathway; 4: mtFASII pathway; 5: cholesterol synthesis; 6: glycolysis; 7: folate, serine, and nucleotide metabolism; 8: FAO and sphingolipid de novo synthesis; and 9: glycerophospholipid metabolism. Metasamples of EB MFs + alvMFs and alvMFs + SPM cells are shown at one chart as they are extremely close in their metabolic characteristics.

Figure 5.. Cell types shared between IG MNP OS, IG MNP P1, and mTMS datasets have similar patterns of metabolic modules signatures

(A) Population memberships across the datasets: prog, progenitor; SC, stem cell; MLP, multi-lineage progenitor; MF YS, yolk sac macrophage; MF EB, embryoid body macrophage; MF, macrophage; alvMF, alveolar macrophage; SPM, small peritoneal macrophage; MG, microglia; KC, Kupffer cell; Mo, monocyte; pDC, plasmacytoid dendritic cell; DC, dendritic cell; migDC, migratory dendritic cell; NP, neutrophil (Table S3). (B) Enrichment of individual metabolic modules across all datasets obtained during GAM-clustering analysis of IG MNP OS dataset.

Figure 6.. Subnetworks associated with early developmental stages and DCs

(A and E) Heatmaps of module patterns along with the expression of some of its genes or genes related to the same biological subject (from lowest as blue to highest as red). YS MF, yolk sac macrophage; EB MF, embryoid body macrophage; alvMF, alveolar macrophage; SPM, small peritoneal macrophage; MG, microglia; MF, macrophage; Mo, monocyte; DC, dendritic cell; pDC, plasmacytoid DC; migDC, migratory DC. (B, C, and F) Metabolic modules per se where edges of modules are attributed with color according to correlation of its enzyme’s gene expression to this particular module pattern and thickness according to its score. (D) Enrichment of modules genes expression (from lowest as blue to highest as red, transparent dots correspond to treated samples) across all three analyzed datasets: IG MNP OS, IG MNP P1, and mTMS datasets. (G and H) Flow cytometry analysis of DC staining with cholesterol-dependent cytolysin perfringolysin O (PFO) that indicates the level of cholesterol in the cell membrane. (H) Mean fluorescence intensity (MFI) levels of PFO binding in DC subsets. n = 4 mice per group (2 male and 2 female); each dot indicates an independent mouse. Statistics by one-way ANOVA with Dunnett’s multiple comparison test. (I) DC migrations experiment scheme. (J) Total plasma cholesterol levels in control and treated-with-simvastatin animals; n = 5 mice in each group; statistical analysis by unpaired two-tailed t test. (K) Percentage of migrated FITC⁺CD11c⁺ DCs in draining lymph nodes after FITC application in control and treated-with-simvastatin animals; n = 20 mice in each group; statistical analysis by unpaired two-tailed t test. **** p < 0.0001. (H, J, and K) Data shown as mean ± standard error of the mean.

Figure 7.. Subnetworks associated with fatty acid synthesis and degradation

(A) Heatmaps of module patterns along with the expression of some of its genes (from lowest as blue to highest as red). YS MF, yolk sac macrophage; EB MF, embryoid body macrophage; alvMF, alveolar macrophage; SPM, small peritoneal macrophage; MG, microglia; MF, macrophage; Mo, monocyte; DC, dendritic cell; pDC, plasmacytoid DC; migDC, migratory DC. (B) Enrichment of module gene expressions (from lowest as blue to highest as red, transparent dots correspond to treated samples) across all three analyzed datasets: IG MNP OS, IG MNP P1, and mTMS datasets. (C) Metabolic modules per se and corresponding schematic diagrams. Edges of modules are attributed with color according to correlation of its enzyme’s gene expression to this particular module pattern and with thickness according to its score. (D) Schematic representation of metabolic module. (E) Schematic illustrating the design of the experiment with mouse peritoneal (Per) and alveolar (Alv) macrophages (MΦs) treated with BSO for 12 h to inhibit GSH synthesis followed by activation by zymosan for 4 h. (F) Secretion levels of PGE2; n = 3 mice per group; statistics by unpaired two-tailed t test. N.S., non-significant (p > 0.05). (G) Secretion levels of cysteinyl leukotrienes; n = 3 mice per group; statistics by unpaired two-tailed t test. (F and G) Data shown as mean ± standard error of the mean.