Article https://doi.org/10.1038/s41467-025-60295-2 MACanalyzeR scRNAseq analysis tool reveals PPARγHIGH/GDF15HIGH lipid-associated macrophages facilitate thermogenic expansion in BAT Andrea Ninni1,2, Fabio Zaccaria1,2, Luca Verteramo1,2, Francesca Sciarretta1, Loreana Sanches Silveira 3, José Cesar Rosa-Neto3, Simone Carotti4, Lorenzo Nevi 4, Paolo Grumati 5,6, Satish Patel7, Giulia Carrera8,9, Alessandro Sgambato10,11, Donatella Lucchetti10,11, Filomena Colella10, Ilenia Severi12, Martina Senzacqua12, Antonio Giordano12,13,14, Sergio Bernardini1, Claudia Di Biagio 1, Flavia Tortolici1, Giuseppe Rizzo 15, Clement Cochain15,16, Valerio Chiurchiù 8,17, Stoyan Ivanov 18, Beiyan Zhou 19, Jesse W. Williams 20,21, David B. Savage7, Katia Aquilano 1 & Daniele Lettieri-Barbato 1 Macrophages are key regulators of adipose tissue plasticity. Obesity impairs brown adipose tissue (BAT) function in humans, yet macrophage-mediated mechanisms remain elusive. Here, we introduce MACanalyzeR, a single-cell RNA sequencing (scRNAseq) tool designed for comprehensive monocyte/ macrophage metabolic profiling. Applying MACanalyzeR to BAT from obese male murine models (db/db and HFD-fed mice), we identify lipid-associated macrophages (LAMs) with foamy characteristics. Unlike db/db BAT LAMs, those in HFD BAT correlate with thermogenic gene expression and PPAR sig- naling activation. A distinct PpargHIGH LAM subcluster progressively accumu- lates in thermogenically active BAT. Macrophage-specific Pparg depletion disrupts BAT thermogenesis, inducing a white-like phenotype and metabolic dysfunctions. Mechanistically, PpargHIGH LAMs secrete GDF15, a key regulator of BAT identity and lipid metabolism under high-energy demand. Our study establishes MACanalyzeR as a powerful tool for immunometabolic interroga- tion and identifies PpargHIGH LAMs as critical mediators of BAT homeostasis. Single-cell RNA sequencing (scRNAseq) is amethodused to study gene expression at the single-cell level. This is in contrast to bulk RNA sequencing, which provides information only on the average expres- sion of the predominantly represented cell population in a hetero- geneous cellular system such as that of a tissue. Compelling evidence highlights the importance of brownadipose tissue (BAT)macrophages in governing adaptive thermogenesis1. Recent findings have reported that macrophages participate in BAT physiology at different levels, including release of pro-thermogenic cytokines2, controlling the mitochondrial quality by transmitophagy3,4, or participation in cell profiling adipocyte precursors5. Of note, most trials studied macro- phage dynamics in a context in which BAT was thermogenically acti- vated through exposure of mice to low temperatures. Although some authors have suggested an increase in proinflammatory macrophages Received: 21 July 2023 Accepted: 13 May 2025 Check for updates A full list of affiliations appears at the end of the paper. e-mail: daniele.lettieri.barbato@uniroma2.it Nature Communications | (2025) 16:5063 1 12 34 56 78 9 0 () :,; 12 34 56 78 9 0 () :,; http://orcid.org/0000-0003-4429-0153 http://orcid.org/0000-0003-4429-0153 http://orcid.org/0000-0003-4429-0153 http://orcid.org/0000-0003-4429-0153 http://orcid.org/0000-0003-4429-0153 http://orcid.org/0000-0003-0924-2551 http://orcid.org/0000-0003-0924-2551 http://orcid.org/0000-0003-0924-2551 http://orcid.org/0000-0003-0924-2551 http://orcid.org/0000-0003-0924-2551 http://orcid.org/0000-0002-9942-9389 http://orcid.org/0000-0002-9942-9389 http://orcid.org/0000-0002-9942-9389 http://orcid.org/0000-0002-9942-9389 http://orcid.org/0000-0002-9942-9389 http://orcid.org/0000-0001-9096-3166 http://orcid.org/0000-0001-9096-3166 http://orcid.org/0000-0001-9096-3166 http://orcid.org/0000-0001-9096-3166 http://orcid.org/0000-0001-9096-3166 http://orcid.org/0009-0008-5366-0738 http://orcid.org/0009-0008-5366-0738 http://orcid.org/0009-0008-5366-0738 http://orcid.org/0009-0008-5366-0738 http://orcid.org/0009-0008-5366-0738 http://orcid.org/0000-0002-0316-7311 http://orcid.org/0000-0002-0316-7311 http://orcid.org/0000-0002-0316-7311 http://orcid.org/0000-0002-0316-7311 http://orcid.org/0000-0002-0316-7311 http://orcid.org/0000-0002-0527-2297 http://orcid.org/0000-0002-0527-2297 http://orcid.org/0000-0002-0527-2297 http://orcid.org/0000-0002-0527-2297 http://orcid.org/0000-0002-0527-2297 http://orcid.org/0000-0003-0638-2428 http://orcid.org/0000-0003-0638-2428 http://orcid.org/0000-0003-0638-2428 http://orcid.org/0000-0003-0638-2428 http://orcid.org/0000-0003-0638-2428 http://orcid.org/0000-0003-3815-0891 http://orcid.org/0000-0003-3815-0891 http://orcid.org/0000-0003-3815-0891 http://orcid.org/0000-0003-3815-0891 http://orcid.org/0000-0003-3815-0891 http://orcid.org/0000-0002-5905-9870 http://orcid.org/0000-0002-5905-9870 http://orcid.org/0000-0002-5905-9870 http://orcid.org/0000-0002-5905-9870 http://orcid.org/0000-0002-5905-9870 http://orcid.org/0000-0002-1172-814X http://orcid.org/0000-0002-1172-814X http://orcid.org/0000-0002-1172-814X http://orcid.org/0000-0002-1172-814X http://orcid.org/0000-0002-1172-814X http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-60295-2&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-60295-2&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-60295-2&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-60295-2&domain=pdf mailto:daniele.lettieri.barbato@uniroma2.it www.nature.com/naturecommunications in theBATofmousemodels ofobesity,mostworks have confirmed the role of macrophages in impairing the functionality of white adipose tissue (WAT)6. Through the use of scRNAseq methods, however, it emerges that a subclass of macrophages, commonly annotated as lipid-associatedmacrophages (LAMs), play a protective role in limiting the metabolic degeneration observed in both genetic and dietary models of obesity7,8. It has also been reported that LAMs develop a foamy phenotype and increase metabolic activity and lysosomal mass9,10, such as some authors have shown an increased capacity of LAMs to secrete proteins by virtue of their senescent phenotype11,12. Taken together, these data suggest that there is still no common consensus on the role of BAT macrophages in conditions of obesity. To comprehensively investigate monocyte/macrophage dynam- ics in BAT during obesity, we developed MACanalyzeR, a novel com- putational framework for in-depth single-cell RNA sequencing (scRNAseq) analysis. MACanalyzeR leverages a multi-layered approach, encompassing three distinct modules to comprehensively characterize these crucial immune cells: FoamSpotteR for identifying foamy-like features, MacPolarizeR for assessing polarization status, and PathAnalyzeR for dissecting metabolic phenotypes. By applying MACanalyzeR to scRNAseq data from BAT of obese models, we aimed to unravel the metabolic and phenotypic heterogeneity of BAT mac- rophages and gain deeper insights into their role in obesity-associated BAT responses. In line with the data reported for WAT, in BAT of obese mice we observed a significant increase in the cluster pertaining to LAM which showed a significant enrichment in the metabolic pathways related to glycolysis and fatty acid metabolism. The detailed analysis of pheno- typic andmetabolic parameters reveals that BAT LAMs in diet-induced thermogenesis (HFD model) exhibit unique characteristics related to healing phenotype and PPAR signaling pathway activation. Accumu- lation of PpargHIGH LAMswas alsoobserved in thermogenically stressed BAT. Notably, mice with conditional PPARγ deletion in the myeloid lineage (MACPPARGKO mice) exhibited a loss of BAT thermogenic iden- tity, along with associated metabolic abnormalities. Computational predictions, validated in irradiated HFD-fed mice transplanted with bone marrow from Gdf15 knockout donors (BMTGDF15KO), revealed that PpargHIGH LAMs release GDF15 that supports themaintenance of brown fat identity and its lipogenic characteristics in HFD-fed mice. In con- clusion, we can suggestMACanalyzeR as a versatile computational tool to dissect the macrophage responses from data obtained by scRNA- seq. Our predictions identify PpargHIGH LAMs as a metabolically active macrophage subclass that support BAT identity and promotes meta- bolic health through GDF15 release. Results Lipid-associated macrophages accumulate in BAT of murine models of obesity While many studies have discussed the homeostatic role of macro- phages in BAT during physiological conditions like cold exposure, the behavior of macrophages in conditions deactivating BAT such as obesity remains largely unexplored. To address this gap, we aimed to investigate the phenotypic response of BAT in a genetic model and dietarymodel of obesity such asmicewith type 2 diabetes (db/db) and mice fed with a high-fat diet (HFD), respectively. We initially analyzed bulk RNA sequencing (RNAseq) data from the BAT of db/db and HFD mice, focusing on the expression levels of canonical brown orwhite fat genes. Notably, the BAT of db/dbmice exhibited a reduced identity of brown adipocytes and a more white-like phenotype compared to HFD mice, which preserved their brown fat identity (Fig. 1A). Proteomics data further confirmed these results, showing a loss of thermogenic markers in the BAT of db/dbmice, while BAT fromHFDmice displayed increased levels of thermogenic proteins, including UCP1 (Fig. 1B). To strengthen the hypothesis thatmacrophagesmay play a role in the loss or maintenance of thermogenic function and brown adipocyte identity, we computed our pathway score for the Gene Ontology (GO) pathways (Materials and Methods) related to thermogenic function and macrophages activation for each sample. Subsequently, we con- structed two linearmodels basedon these scores (Fig. 1C), which led to the identification of an inverse relation between macrophage activity and thermogenic function in db/db mice. Interestingly, the HFD mice showed a positive correlation between these two scores, suggesting that under HFD macrophages sustain BAT thermogenesis. Building on these results, we gained valuable insights into the dynamics and transcriptomic changes of immune cells, especially macrophages, in the BAT of db/db and HFD mice. To achieve this, we analyzed single-cell RNA sequencing (scRNAseq) data collected from stromal vascular fraction (SVF) isolated from BAT in both db/db and HFD-induced obese mice (Fig. 1D). After rigorous quality filtering, the analysis led to the identification of a total of 14,715 cells ( ~4900 cells sampled per condition: WT, db/db and HFD). Based on differential gene expression and well-known literature markers, we identified 6 distinct clusters that were annotated as Fibro-Adipose Progenitors (FAPs), Adipose Stromal Cells (ASCs), Endothelial Cells, T Cells and monocytes/macrophages (MonoMacs) (Fig. 1E and Supplemental Fig. 1A). Of note, among all immune cells the MonoMacs cluster is notably the largest in cell number in BAT SVF (Fig. 1F). The population of MonoMacs massively increased in BAT from both db/db and HFD mice when compared with controls (Supple- mental Fig. 1B). Importantly, a recent study highlights that scRNAseq experiments may exhibit biases towards specific cell types13. Conse- quently, while we acknowledge the potential presence of technology- mediated biases in our data, the strong reproducibility observed across the two obese conditions and other studies7 indicates that our compositional estimates remain largely unaffected by such technical factors. In order to givemore insight onMonoMacs dynamics in obese BAT, we sub-clustered and annotated these cells (n = 5578) using well- known biomarkers of adipose tissue macrophages from the literature (Supplemental Fig. 1C). This analysis led to identification of four dis- tinct subclusters (Fig. 1G), including one composed of monocytes expressing Ccr2, Cd44, Cx3cr1, one containing Lyve1HIGH and Mrc1HIGH macrophages denoted as perivascular-like macrophages (PVMs), and two with a very similar expression pattern, which can be attributed to the previously reported lipid-associatedmacrophages (LAMs). Among these two latter groups, one is characterized by genes involved in cellular proliferation (Pola1, Kif11, Kif15), leading us to designate these macrophages as proliferatingmacrophages (proliferating). In line with the massive accumulation of LAMs in expanding WAT and BAT7,14, a significant accumulation in LAMs was also observed in BAT of db/db and HFD mice (Fig. 1H, I). Consistent with scRNAseq data, BAT mac- rophages were labeled with TREM2, a well-known marker of LAM. We observed a significant emergence in the number of TREM2+ macro- phages in the BAT of HFD-fed mice (Fig. 1J and Supplemental Fig. 1D). MACanalyzeR identifies foamy-like macrophages in obese BAT LAMs share key characteristics with foam macrophages (fMACs), which are well-known for their role in atherosclerotic plaques. Like fMACs, LAMs modulate inflammation and release important factors that suggest promoting cell proliferation and differentiation, thereby supporting the regeneration process15. However, a clear description of these foamy-like features in LAMs resident in adipose tissue is lacking. Furthermore, while the LAM with foaming feature within obese WAT has been documented, its occurrence in BAT remains unexplored. To thoroughly decipher the foamy-like features in BAT macrophages, we developed FoamSpotteR, a module based on a machine learning classifier that utilizes cellular transcriptional profiles at a single-cell level (Fig. 2A). To train this model, we analyzed scRNAseq data col- lected from aortic CD45+ (Ptprc) cells isolated from healthy mice and atherosclerotic mice from which we extracted macrophages and annotated them as fMAC- and fMAC+ based on specific markers found Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 2 www.nature.com/naturecommunications in literature16,17 (Fig. 2B and Supplemental Fig. 1E–H). For model training,we utilized 75%of the total cells,while the remaining 25%were used as a test set. The model demonstrated high efficiency and accu- racy in discriminating fMAC+ in single-cell data. Next, we utilized the trained FoamSpotteR model to analyze macrophages in BAT and observed a notable presence of fMACs in both db/db and HFD mice (Fig. 2C, D). To determine the emergence of fMACs within a specific macrophage subcluster, we calculated the proportion of fMACs in each subpopulation. As expected, we observed that the occurrence of fMACs predominantly corresponds to LAMs and proliferating Fig. 1 | scRNA analysis reveals LAM accumulation in BAT of obese mouse models. A Heatmap of BAT and WAT markers in bulk RNAseq data of HFD and db/dbBAT.BHeatmap of thermogenic proteins in proteomics data of db/db (upper panel) and HFD (lower panel) BAT. C Linear models fitted on scores computed on GO terms Macrophage Activation pathway and Adaptive Thermogenesis for BAT db/db data (upper panel) and BAT HFD data (lower panel). P-value by two-side ANOVA test. D Schematic figure delineating the scRNAseq experimental design. E UMAP representing cell clusters identified by scRNAseq of SVFs (n = 14,715) iso- lated from BAT of wild type (WT), db/db and high fat diet (HFD) mice. Clustering identified 5major clusters: Fibro-Adipose Progenitors (FAP), Adipose Stromal Cells (ASC), Endothelial Cell, T cells and Monocytes/Macrophages (MonoMacs). F Barplot representing the cell proportion across the 3 conditions. G UMAP of MonoMacs sub-clustering (n = 5578) identified by scRNAseq analysis. Clustering identified 4 major clusters: Monocytes, Perivascular-like Macrophages (PVM), Lipid-AssociatedMacrophages (LAM) and Proliferating Macrophages. (H) UMAP of BAT MonoMacs splitted by condition. I Bar plot of MonoMacs subcluster abun- dance across the 3 conditions. J Multiplex immunofluorescence stained for Adi- ponectin (orange), UCP1 (yellow), TREM2 (green) and CD9 (grey) in paraffin- embedded tissue section from low fat diet (LFD) and HFD BAT. Nuclei stained with DAPI (blue). Scale bar, 20μm. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 3 www.nature.com/naturecommunications macrophages, constituting 57% and 27% respectively. To validate the predictions generated by FoamSpotteR and clarify the differences between fMACs and LAMs, we compared our results with those from a publicly available dataset from bulk RNA seq (GSE116239), which includes foamy and non-foamy macrophages isolated from athero- sclerotic aortas via lipid probe-based flow cytometry sorting17. We conducted a comparative analysis by examining the differentially expressed genes (DEGs) of fMACs from the bulk RNAseq and those identified through our FoamSpotteR module in the previously anno- tated LAM subcluster (fMAC+ LAMs vs. fMAC− LAMs). Our investigation revealed 593 upregulated genes in fMACs frombulkRNAseq and 193 in fMAC+ LAMs from our scRNAseq, with 73 genes shared between the two datasets. Notably, the commonly upregulated pathways are associated with lipidmetabolism and localization (Fig. 2E). In contrast, pathways uniquely upregulated in fMACs from atherosclerotic plaque predominantly involve cholesterol metabolism (Fig. 2F). Interestingly, the pathways uniquely upregulated in fMACs from BAT are linked to the generation of metabolite precursors, energy production, and vesicle organization (Fig. 2G). This suggests that macrophages in BAT develop a foaming-like phenotype contributing to tissue homeostasis through mechanisms that differ from canonical atherosclerotic foam macrophages. BAT LAMs show a different polarization status between db/db and HFD mice Adipose tissue macrophages (ATMs) with an M1-like phenotype are well-established contributors to impaired insulin sensitivity in both mice and humans4,9, the polarization status of fMACs in BAT of obese mouse models remain unexplored. To investigate this, we developed MacPolarizeR, a machine learning-based module that classifies mac- rophages according to their activation state.MacPolarizeR operates by clustering cells according to the expression of genes associated with M1 and M2 polarization, selected from a public dataset (GSE129253) including bulk RNAseq data of M1 and M2 macrophages polarized in vitro. Following this training process, MacPolarizeR can identify three distinctmacrophage populations: Inflammatory (M1-like), Transitional (M0-like), and Healing (M2-like). In addition to this classification, MacPolarizeR calculates an M1 and M2 polarization score for each cell based on the expression of the previously selected genes. This score system facilitates a more comprehensive and nuanced characteriza- tion of the macrophage polarization state. The calculated scores can be visualized on a two-dimensional plot, with M1 score on the X-axis and M2 score on the Y-axis (Fig. 3A). This graphical representation provides a clearer interpretation of macrophage behavior, indicating whether subpopulations lean towards a healing phenotype (upper Fig. 2 | FoamSpotteR identifies foamy-like macrophages in obese BAT. A Graphical representation of the FoamSpotteR module workflow. B UMAP of macrophages isolated from atheromatous plaque (GSE97310 and GSE116240). Clustering identified twomajor clusters: foammacrophages (fMAC+) and non-foam macrophages (fMAC-). C, D UMAP of BAT macrophages colored by FoamSpotteR result (C) and splitted by sample (D). Values above 0.5 of FoamDEX identify foamy- like macrophages. E–G Gene Set Enrichment for GO Biological Process of shared upregulated genes in bulk RNAseq of foam macrophages and FoamSpotteR pre- dicted fMAC+ (E), uniquely upregulated in foam macrophages (F) and uniquely upregulated in fMAC+ predicted cells (G) (p-value <0.05; log2FC > 1). P-value by BH test. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 4 www.nature.com/naturecommunications left), a transitional state (upper right and lower left), or an inflamma- tory state (lower right). To validate the accuracy of our polarization classification, we appliedMacPolarizeR to a scRNA sequencing dataset derived from in vitro polarizedmacrophages (GSE117176). Specifically, bonemarrowderivedmacrophages (M0phenotype)were treatedwith IFN-γ and LPS to induce anM1 phenotype, and IL-4 plus IL-13 to induce an M2 phenotype. The analysis of this dataset gave rise to the identi- fication of 3 distinct clusters based on the forced polarization status that we used as input for MacPolarizeR (Supplemental Fig. 2A). Although the identification of inflammatory macrophages was very precise, there was an overlap between transitional and healing mac- rophages related to the similarity between these two subclasses, as the original authors pointed out (Li et al. 2019) (Supplemental Fig. 2B). Noteworthy, the polarization score analysis highlights a less pro- nounced overlap between M0 and M2, confirming that the prediction made by MacPolarizeR perfectly aligns with the polarization observed in vitro (Supplemental Fig. 2C). For further validation, we applied MacPolarizeR to analyze a scRNAseq dataset derived from theWAT of mice at one day and one-month post-sepsis induction (PRJNA626597), a representative model system for studying in vivo inflammation (Supplemental Fig. 2D). This dataset represents a perfect fit for our module, as it includes two different inflammatory status: acute (one day post sepsis) and chronic (one month post sepsis). After the scRNAseq analysis, the MonoMacs cluster was extracted from the dataset and used as input for MacPolarizeR (Supplemental Fig. 2E). As expected, after one day of sepsis the macrophage component shifted towards an inflammatory phenotype (M1). In contrast, after onemonth of sepsis, the macrophages acquired a completely different pheno- type, predominantly residing in the M2 quadrant of the MacPolarizeR plot, indicative of healing macrophages. In contrast, macrophages in control WAT exhibited an intermediate localization on the graph, suggesting an inflammation-suppressed phenotype (Supplemental Fig. 2F–G). The use ofMacPolarizeR not only provided insights into the activation states of macrophages but also demonstrated its effective- ness as a tool for offering a more granular and clearer depiction of specific phenomena. By using MacPolarizeR, we were able to trace a sequential activation of macrophages from early to late time points following a specific stimulation. This progression began with a predominantly-inflammatory contribution during the acute phase of the pathology, which then transitioned to a healing state. Given the high accuracy achieved, we applied the MacPolarizeR module to the single-cell profiles of BAT macrophages to analyze the differences between db/db and HFD mice. Consistent with the increased macro- phage component in BAT of the two obese models, our clustering analysis revealed a significant rise in transitional macrophages, a phenomenon that was almost absent in BAT of control mice. Specifi- cally, we observed that the increase in transitional macrophage com- ponents coincided with the accumulation of LAMs in both db/db and HFD mice. In contrast, monocytes exhibited inflammatory behavior, while PVMs primarily displayed healing properties (Fig. 3B–D). To refine polarization characteristics, we conducted a scoring analysis to assessM1orM2polarization. Interestingly, LAMsofHFDmice revealed the acquisition of a healing phenotype in contrast to the LAM of db/db mice, which exhibited a transitional phenotype (Fig. 3E and Supple- mental Fig. 2H). The healing polarization of macrophages in HFDmice suggests a distinct response of BAT compared to db/db mice. BAT LAMs from HFD mice showed an enhanced fatty acid metabolism Given the growing recognition of macrophage metabolism in reg- ulating polarization and tissue function18, we aimed to thoroughly investigate the metabolic signatures of LAMs in relation to thermo- genic versus deactivated BAT. To this end we developed a novel scoring system that compares the expression levels of pathway-related genes under various conditions. This system incorporates cell-specific gene expression, allowing for a more precise analysis of metabolic changes. This approach departs from traditional techniques like over- representation analysis, which neglects cellular expression data, and gene set enrichment analysis (GSEA), which canbe influencedby genes with low or no expression. Our method calculates a pathway activity score by averaging the weighted gene expression of all genes within a pathway across a specific cell type. This score, applicable at both the condition and cell cluster level, provides a clear indication of pathway upregulation (score >1) or downregulation (score <1). Next, we ana- lyzed scRNAseq data from the BAT of db/db and HFDmice to compare the enrichment of KEGG pathways in LAMs. PathAnalyzeR revealed greater enrichment of pathways associated with the pentose phos- phate pathway, glycolysis, and the TCA cycle in both obesity condi- tions compared to the control (Fig. 4A, B). These results underscore LAMs as BAT-resident immune cells with the highest level ofmetabolic Fig. 3 | MacPolarizeR shows a less inflammatory phenotype of LAMs in obese BAT. A Graphical cartoon of MacPolarizeR workflow. B, C Circular bar plot representing MacPolarizeR’s classification across conditions (B) and macrophage subcluster (C). D UMAP of BAT macrophages colored by MacPolarizeR result and splitted by conditions. E 2D-plot of MacPolarizeR M1 and M2 score across conditions. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 5 www.nature.com/naturecommunications activation. However, a key difference emerged that BAT LAM from db/db mice exhibited a prominent mitochondrial oxidative phos- phorylation (OXPHOS) metabolism, while HFD LAM primarily dis- played increased fatty acid (FA) metabolism and lysosome upregulation. To validate these findings, we generated heatmaps for each condition depicting genes belonging to these three pathways, which confirmed the results obtained from PathAnalyzeR (Fig. 4C). PathAnalyzeR enables single-cell pathway scoring, allowing for the comparison of metabolic profiles across conditions or cell types. Using this tool, we found that db/dbmacrophages have high oxidative metabolism, while both db/db and HFD macrophages show increased FA metabolism compared to controls (Fig. 4D and Supplemental Fig. 2I). Additionally, db/db and HFD MonoMacs exhibit higher lyso- somal and mitochondrial components than controls (Fig. 4E and Supplemental Fig. 2I). Next, CD45+/F4.80+ macrophages isolated from the BAT of HFD mice were labeled with MitoTracker and Lysotracker Red. Consistent with the computational predictions, a significant increase in both mitochondrial and lysosomal mass was observed (Supplemental Fig. 3A). To clarify the type of FAmetabolism alteration operating in the two conditions, we computed pathway enrichment scores for specific GO terms related to lipidmetabolism (i.e., catabolic vs biosynthetic). Notably, while therewere no significant differences in FA catabolism between the two conditions, our analysis revealed an upregulation of FA biosynthetic processes in HFD macrophages (Fig. 4F). To investigate the cell type-specific regulation of OXPHOS and fatty acid biosynthesis pathways, we plotted pathway scores for these pathways, separating data by both condition (db/db and HFD) and cell type (Fig. 4G). Our analysis revealed that both conditions primarily upregulated these pathways in LAMs, albeit with distinct patterns. db/db LAMs exhibited a stronger orientation towards OXPHOS pathways, while HFD LAMs showed greater engagement in fatty acid biosynthesis. The observed differential metabolic responses of LAMs in db/db and HFD mice may be attributed to the distinct obesity-inducing mechanisms of these models. This divergence in macrophage metabolic phenotypes could be linked to the observed preservation of brown adipocyte identity in BAT of HFD mice, which feature is absent in db/db mice. Thermogenically-competent LAMs showed an enhanced PPAR pathway activation To explore the transcriptional trajectories that committed macro- phages towards different inflammatory andmetabolic phenotypes, we performed velocity analysis on BAT MonoMacs from both db/db and HFDmodels. The resulting vectors were projected onto the respective UMAP plots (Fig. 5A, B). As expected, the velocity projections reveal that LAMs originate from the differentiation of monocytes and PVMs. A key distinction between the two models is observed in monocyte velocity: while the HFD model exhibits clear directional information for monocyte vectors, the db/db model lacks this guidance. This highlights distinct differentiation patterns of BAT macrophages Fig. 4 | PathAnalyzeR describes metabolic differences in BAT LAMs of db/db and HFD mice. A, B Heatmaps of PathAnalyzeR analysis using KEGG metabolic pathways related to roles or components unique to macrophages. Plotted by condition (A) and by macrophage subcluster (B). C Heatmap representation of gene expression of LipidMetabolism, OXPHOS and Lysosomepathways. Plotted by condition (upper panel) and by macrophage subclusters (lower panel). D, E PathAnalyzeR 2D density plots splitted by condition using OXPHOS and Fatty Acid metabolism KEGG pathways (D) Lysosome and Mitochondria GO terms. F 2D density plots of GO Fatty Acid Biosynthetic and Catabolic processes splitted by condition. G KEGG Fatty Acid Biosynthetic Process and Oxidative Phosphorylation 2D plots splitted both by condition and cell type. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 6 www.nature.com/naturecommunications between the two obesity models. Specifically, the pronounced healing and fatty acid biosynthetic behavior of LAMs in the HFD model may result from a unique molecular differentiation mechanism in mono- cytes, which is absent in the db/dbmodel. This suggests that different molecular processes drive the observed behaviors in each model. To investigate this, and since the differences between the db/db and HFD models were observed in the monocyte cluster, we focused the gene ranking analysis specifically on this subpopulation. This analysis identified 481 genes in db/db and 345 in HFD, with 90 overlapping genes. To elucidate the molecular basis of the distinct polarizations and metabolic phenotypes, we identified the top 10 hub genes among the 345 unique genes differentially expressed in the HFD group (Fig. 5C). Notably, this set included Pparg, a well-known gene reg- ulating the anti-inflammatory response and fatty acid biosynthesis19,20. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 7 www.nature.com/naturecommunications Interestingly, Pparg expression was significantly increased in BAT LAMs of HFDmice. The hypothesis that the different behaviorsmay be connected to the PPARγ signaling pathway is further confirmed by the PathAnalyzeR module, as an upregulation of Pparg was detected in HFD LAM (Fig. 5D). In addition, macrophages (CD45+/F4.80+) isolated from the BAT of HFD mice displayed elevated PPARγ protein levels compared to those from LFD mice (Supplemental Fig. 3B). The collected results led us to postulate that LAMs with high Pparg expression (PpargHIGH) take part in maintaining the BAT identity during HFD. To demonstrate this, we took advantage of mice with conditional PPARγ deletion in the myeloid lineage (MACPPARGKO mice). The mice were fed for 10 weeks with a high fat diet. The analysis of weight gain between LFD and HFD groups demonstrated a significant increase in weight for the obese model compared to the control (Fig. 5E). Subsequent comparison of the two HFD models, WT and MACPPARGKO, revealed a significant trend towards greater weight gain in KO mice compared to the WT. Additionally, a more in-depth exam- ination of BAT weight revealed a significant increase in KO mice compared to WT (Fig. 5F). To further characterize the thermogenic profile of BAT, we per- formed a comparative proteomic analysis betweenMACPPARGKO andWT mice. This analysis revealed a significant disruption in the thermogenic program of BAT in MACPPARGKO mice challenged with a HFD (Fig. 5G). Corroborating this functional impairment, histological examination demonstrated amarked shift towards a white-like AT phenotype in the BAT of MACPPARGKO mice (Fig. 5H). This aberrant BAT identity in MACPPARGKO mice was associated with a systemic metabolic deteriora- tion, manifested by impaired fasting glucose tolerance and elevated circulating cholesterol levels (Fig. 5I, J). Based on these findings, we asked if transplanting BAT LAMs is effective to enhanceBATmetabolic capacity including thermogenesis and lipid metabolism. To test this, macrophages (CD45+/F4.80+ cells) isolated from the BAT of LFD and HFDmice were transplanted (105 cells) into the intrascapular region of adult mice. To assess whether the transplantation of BAT PpargHIGH macrophages effectively enhances thermogenic and metabolic phe- notypes, the transplanted mice were exposed to cold (4 °C) for 24 h. Consistent with the role of LAM in sustaining metabolic capacity in BAT, our results show that mice transplanted with BAT HFD LAMs display an increased UCP1 protein level (Fig. 5K) and a stronger expression of genes involved in thermogenesis and lipid metabolism compared to BAT from control mice (Fig. 5L). These findings collec- tively underscore the pivotal role of PpargHIGH LAMs in BAT thermo- genesis and maintaining metabolic homeostasis in response to a high- fat dietary and cold exposure challenge. PpargHIGH LAMsmaintain BAT identity through GDF15 signalling To elucidate the mechanisms by which PpargHIGH LAMs contribute to the molecular reorganization of BAT during HFD, we analyzed ligand- receptor interactions using scRNAseq data from db/db and HFD mice, employing the CellChat tool21. This comparative analysis revealed four unique communication pathways specifically enriched in HFD mice (Fig. 6A and Supplemental Fig. 3C).We prioritized GDF15 as a potential key mediator of macrophages in maintaining BAT identity under HFD conditions. This aligns with existing data showing that a specific macrophage subpopulation regulates regenerative inflammation through GDF1522. Remarkably, comparative analysis between myeloid and non-myeloid cells in BAT revealed that Gdf15 was mainly up- regulated in HFD LAMs (Fig. 6B, C). To investigate the role of PPARγ in the development of LAMsignatures andGdf15 expression,we analyzed publicly available RNA-seq data (GSE111105) from BMDMs of WT and PpargKO mice, both treated with rosiglitazone, a selective PPARγ ago- nist (Fig. 6D). As shown, rosiglitazone induced Gdf15 expression alongside canonical LAM genes. In contrast, PPARγ ablation led to a significant downregulation of LAM genes and a reduction in Gdf15 expression. This was further confirmed by Western Blot analysis (Fig. 6E), which supports the role of Pparg in regulating LAM features, including Gdf15 expression. To obtain more insight about the mechanism through which PPARγ drives GDF15 induction specifically in BAT LAMs, we isolated LAM from BAT of HFD mice and employed a well-established PPARγ inhibitor, GW9662, a selective antagonist that effectively prevents PPARγ activation. Consistent with the findings in BMDM (Fig. 6D, E), macrophages (CD45+/F4.80+) isolated from the BAT of HFD mice dis- played elevated GDF15 protein levels compared to those from control mice (Supplemental Fig. 3D). Notably, treatment with GW9662 sig- nificantly reduced both GDF15 protein levels and CD36 (Fig. 6F). These results emphasize the crucial role of PPARγ in regulating LAM sig- nature, including the production of GDF15 inmacrophages. It has been reported that brown adipocytes express GDF1523. To demonstrate that GDF15 is primarily expressed by BAT LAMs and not brown adipocytes inHFDmice,we conducted confocalmicroscopy analyses to assess the co-localization between GDF15 and CD206 (a LAM marker) in tissue sections. Our confocal microscopy data revealed a co-localization between CD206 and GDF15 in the BAT of HFD mice, confirming that GDF15+ macrophages accumulate in this tissue under metabolic stress (Fig. 6G). In addition, numerous CD206+/GDF15+ macrophages were observed in BATcrown-like structures (Fig. 6H), further supporting the hypothesis that GDF15-expressing macrophages are recruited from metabolically challenged tissues, where they contribute to BAT homeostasis. Furthermore, morphometric analysis confirms the increaseof LAMsexpressingGDF15 inHFDBAT (Supplemental Fig. 3E). This evidence was further supported by high dimensional flow cytometry to identify the levels of both PPARγ and GDF15 in BAT LAM cells. Accordingly, LAMs (CD45+/F4.80+/TREM2+/CD9+) identified from the BAT of HFD mice showed higher expression levels of both PPARγ and GDF15 compared to controls (Fig. 6I, J and Supplemental Fig. 3F, G), These findings were corroborated by unchanged levels of these markers in non-LAMs (CD45+/F4.80+/Trem2-/CD9-) (Supple- mental Fig. 3H) and were consistent with the significantly higher release of GDF15 observed in their culture media (Fig. 6K). Recent findings reported that GDF15 is crucial for cold-induced thermogenesis in BATandmaintains subcutaneous adipose tissue lipid metabolic signature5,24. To confirm the involvement of PpargHIGH/ GDF15HIGH LAM in maintaining BAT identity and regulating Fig. 5 | PPARγ LysMCremice show an increase of body and BATweight and loss of thermogenic proteins. A, B scRNAseq velocities projected onto the UMAP in db/db MonoMacs (A) and HFD MonoMacs (B). C HubGenes of the PPI network of HFD velocity genes. D Expression of Pparg (upper panel) and PPAR signaling pathway (lower panel) splitted by conditions. Box plots show the median as center line, 25% and 75% percentiles as box limits and whiskers extending to the largest and smallest values. E Body weight analysis. Body weight time course during 10 weeks of LFD and HFD of WT and PPARγ LysMCre mice (left panel) and barplot reporting body weight after 10 weeks (right panel). F Barplot reporting BAT weight after 10 weeks of diet. G Heatmap of thermogenesis-related proteins in proteomics data of WT and PPARγ LysMCre mice in LFD and HFD (left panel) and barplot of UCP1 protein quantification (right panel). H H&E staining of BAT slices of HFD WT and HFD PPARγ LysMCre mice. Scale bar 250 μm. I Glucosemia after 10 weeks of LFD and HFD in WT and PPARγ LysMCre mice. J Blood choles- terol levels after 10 weeks of LFD and HFD in WT and PPARγ LysMCre mice. E–G, I, J Data are reported as mean ± SD. P-value by one-side ANOVA test. (n = 4 mice/ group) (K) Representative immunoblot of UCP1 from BAT extracted from mice exposed to cold and transplanted with LFD (Mac LFD) and HFD (Mac HFD) BAT macrophages (left panel) and densitometry analysis (right panel). Tubulin staining was used as loading control (n = 3 mice/group). L Thermogenesis and lipid metabolism genes relative expression in BAT of mice exposed to cold and transplanted with LFD and HFD BAT macrophages. Data are reported as mean ± SD. P-value by unpaired student’s t test. (n = 3 mice/group). Source data are provided in the Source Data file. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 8 www.nature.com/naturecommunications thermogenesis, we analyzed scRNAseq data of immune cells isolated from cold-exposed BAT (GSE207706). Focusing on the MonoMacs component, we observed a significant increase in LAMs after cold exposure (Fig. 6L). Analysis using MACanalyzeR revealed a foamy macrophage phenotype in cold-exposed BAT, similar to that observed in obese BAT (Supplemental Fig. 3I). Importantly, PPAR signaling and Pparg/Gdf15 expression were increased inMonoMacs of cold-exposed mice, aligning with the thermogenic behavior observed in HFD mice (Fig. 6M, N). We also performed a comparative approach analyzing the serum level of GDF15 from db/db, HFD and cold-exposed mice. As expected, we observed the highest GDF15 quantification in HFD and cold-exposed mice serum (Fig. 6O). In contrast, no differences in GDF15 levels were detected in the sera of db/dbmice. Since PPARγwas also deleted in WAT macrophages, we aimed to rule out the Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 9 www.nature.com/naturecommunications contribution of WAT LAMs to the observed metabolic changes. To do this, weperformed a comparative analysis of scRNAseqdata fromboth BAT and WAT (GSE182233) of HFD mice, focusing specifically on the LAM subcluster. This analysis revealed significantly higher activation of the PPAR signaling pathway, along with increased expression of Pparg and Gdf15, primarily in BAT LAMs compared to WAT LAMs (Supplemental Fig. 3J). These findings suggest that BAT macrophages are more inclined to promote a lipogenic phenotype and enhance thermogenic capacity in BAT, in contrast to WAT macrophages. Building on these findings, we explored whether macrophage- derivedGDF15 contributes tomaintaining BAT identity duringHFD. To determine this, adult mice were irradiated with two 5.5Gy doses of radiation. Two hours after irradiation, 10 million bone marrow cells from eithermale GDF15 KO (BMTGDF15KO) orWTmice were injected into the tail vein. Although no significant changes were observed in body weight gain or BAT mass following HFD25, a notable reduction in the expression of thermogenic and lipid elongation genes was detected in the BATof BMTGDF15KOmice (Fig. 6P). These findings suggest thatGDF15 derived from PpargHIGH LAMs may play a role in modulating the ther- mogenic phenotype and lipid metabolism of BAT during HFD. Discussion Macrophages are highly adaptable cells that play a critical role in maintaining the homeostasis of various tissues. While growing evi- dence highlights their involvement in the thermogenic activation of BAT, there is limited research that thoroughly investigates their molecular and metabolic responses within BAT during obesity. The lack of clear evidence regarding the role of MonoMacs in BAT during obesity prompted us to investigate their molecular and metabolic characteristics in mousemodels of obesity, such as the db/db and HFD mice. To better understand the specific features of MonoMacs in BAT, we developed MACanalyzeR, a comprehensive 3-in-1 tool designed to analyze scRNAseq data. Using a machine learning approach, MACa- nalyzeR returns information relating to the degree of foamy-like characteristics, polarization status and metabolic projections of mac- rophages. Consistent with the observations from WAT depots7, LAMs that accumulate in BAT show foamy-like phenotype26. BAT LAMs, characterized by their foamy appearance, exhibited a marked enrich- ment of genes associated with energy production and lysosomal reorganization, setting them apart from the foamy macrophages observed in atherosclerotic plaques. In contrast to db/db mice, we observed that in HFD-fed mice, LAMs are associated with the expres- sion of thermogenic genes and help maintain an energy-dissipating phenotype in BAT. This observation was consistent with the protective role of LAMs, as depletion of TREM2+ LAMs has been shown to exacerbate white adipose tissue (WAT) dysfunction and metabolic abnormalities7. Remarkably, TREM2+ LAMs exert a protective role also in other metabolically-stressed tissues such as the heart27,28. In con- trast, in a db/db mouse model, LAMs were shown to limit the thermogenic identity of mature brown adipocytes26, suggesting that distinct immunometabolic signatures of LAMs may underlie the vary- ing tissue responses. Accordingly, MACanalyzeR revealed distinct metabolic differences between BAT LAMs from db/db and HFD mice. Specifically, BAT LAMs from HFD mice displayed a less inflammatory phenotype, characterized by enhanced lipid metabolism and lysoso- mal activation, resembling the molecular phenotype identified in adi- pose tissue macrophages of humans and obese mice29. LAMs are believed to play an active role in the uptake, storage, and hydrolysis of lipids from extracellular vesicles derived from adipocytes30. Addi- tionally, through the process of lipophagy, macrophages release fatty acids, which can activate the nuclear receptor PPARγ31,32. Consistent with the higher expression of genes involved in fatty acid metabolism and lysosomal function in BAT LAMs from HFD mice, a significant upregulation of genes associated with the PPAR signaling pathway was also observed. Interestingly, several authors reported a regenerative potential of PpargHIGH macrophages33, which was well described in the transcriptional trajectories of macrophages analyzed in regenerating skeletal muscle22,33,34. Although the role of LAM PpargHIGH in adipose tissue physiology is still unclear, disruption of PPAR signaling pathway in myeloid cells impairs alternative macrophage activation, thereby predisposing these animals to development of diet-induced obesity, insulin resistance, and glucose intolerance35. To investigate the role of PpargHIGH macrophages in BAT expansion, we utilized MACPPARGKO mice fed a HFD. Compared to WTmice, MACPPARGKO mice exhibited a failure to sustain the thermogenic phenotype of BAT. Through a PPARγ- mediated mechanism, macrophages were shown to release GDF15, which actively participates in tissue regeneration22. Recent studies have highlighted GDF15’s critical role in regulating energy balance by promoting BAT thermogenesis and preserving the lipid metabolic profile of subcutaneous adipose tissue23,24. Notably, both PPARγ and GDF15 emerge as key features of LAMs identified in thermogenically activated BAT, highlighting their potential importance inorchestrating immunometabolic functions in this context. By transplanting bone marrow from BMTGDF15KO mice into irradiated WT mice to regenerate the hematopoietic lineage under HFD conditions, we demonstrated that BAT macrophages from BMTGDF15KO mice exhibited a lower ther- mogenic potential compared to controls. Consistentwith this findings, recombinant GDF15 reduced obesity and improves glycaemia in rodents onHFDby suppressing appetite and counteracting reductions in energy expenditure, resulting in greater weight loss and mitigation of NAFLD compared to caloric restriction alone36. Although, we identified several limitations, which are described schematically. MACanalyzeR revealed interesting predictions on the different responses of BAT macrophages, but the study lacked certain validation steps. Specifically, the observed differences between BAT LAMs in db/db andHFDmodels could be influenced by treatment time. A time-course experiment would help better understand the dynamics of BAT LAMs during different stages of obesity. Additionally, there Fig. 6 | GDF15 signaling regulates thermogenesis and brown adipose identity. A Chord diagrams of cell-cell communication analysis of MonoMacs versus Fibro- blasts (upper panel) and Venn graph representing the common and uncommon interactions calculated (lower panel). B, C Expression of Gdf15 projected onto the UMAP plot (B) and as a violin plot through conditions (C). D Heatmap of LAM signature and Gdf15 gene expression in WT and PpargKO BMDM treated with Rosiglitazone. E Representative immunoblot of GDF15 and CD36 of BMDM treated with Rosiglitazone for 4, 8, and 16 h. Ponceau staining was used as loading control. F Representative immunoblot of GDF15 and CD36 of BMDM treated with rosigli- tazone and GW9662. Ponceau staining and Vinculin was used as loading control. G, H Confocal microscopy images of BAT from HFD mice showing CD206 (red), GDF15 (green) and DAPI (blue) stainings in macrophages infiltrating BAT (G) and forming crown-like structures (CLS) (H). (E–H replicates: n = 3). I Gating strategy and (J) high-dimensional flow cytometry analysis of PPARγ and GDF15 in LAM of BAT LFD and HFD mice. K Elisa quantification of GDF15 on CD45+/F4.80+ cells isolated from BAT LFD and HFDmice. J, K Data are reported as mean ± SD. P-value by unpaired student’s t test. (n = 3 mice/group). L UMAP representing cell clusters identified by scRNAseq analysis ofMonoMacs (n = 3631) isolated fromBAT of room temperature (ROOM) and 4 °C cold exposed (COLD)mice (GSE207706). Clustering identified four major clusters: Perivascular-like Macrophages (PVM), Lipid- Associated Macrophages (LAM), Classical Monocytes and Activated Monocytes. M ExpressionofPparg (left panel)Gdf15 (right panel) projectedonto theUMAPplot. N Expression of Ppar signaling projected onto theUMAP plot.O Elisa quantification of GDF15 on serum derived from LFD, HFD, db/dbmice and mice exposed to room temperature and 4 °C cold. Data are reported as mean ± SD. P-value by unpaired student’s t test. (n = 3 mice/group). P Thermogenesis and lipid elongation genes relative expression in BAT of BMTGDF15KO mice. Data are reported as mean± SD. P- value by unpaired student’s t test. (n = 7 mice/group). Source data are provided in the Source Data file. Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 10 www.nature.com/naturecommunications were no in vitro experiments to determine whether the metabolism of LAMs is higher than that of non-LAMs.While the PathAnalyzeRmodule predicted increased metabolism in BAT macrophages and partial validation was achieved via flow cytometry, more robust real-time metabolic measurements would strengthen the study’s findings. Incorporating assays such as oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) would provide a more compre- hensive assessment of metabolic signatures and better validate the predicted changes in BAT macrophages. Further investigation with these techniques is necessary to reinforce the conclusions drawn from the current analysis. Furthermore, although BAT from PPARγ-deficient mice on HFD showed impaired thermogenic capacity, exposing these mice to cold conditions would offer a clearer understanding of how PPARγ in LAMs specifically contributes to thermogenesis. This remains a crucial gap in our knowledge and warrants further study. In summary, our study underscores the critical role of macro- phages, particularly LAMs, in orchestrating BAT thermogenesis and metabolic homeostasis under conditions of metabolic stress. By employing MACanalyzeR, we unveiled distinct metabolic and mole- cular adaptations in BAT LAMs during obesity, revealing their dichot- omous roles in HFD and db/db mouse models. Notably, our findings highlight the importance of lipid metabolism, lysosomal activation, and PPAR signaling in shaping the protective and regenerative potential of BAT LAMs. Furthermore, the identification of GDF15 as a pivotal regulator of BAT thermogenesis and systemic energy balance introduces new therapeutic possibilities for combating obesity and its metabolic sequelae. Together, these insights advance our under- standing of macrophage biology in BAT and offer promising avenues for targeting macrophage-mediated pathways to enhance metabolic resilience in obesity and other conditions of metabolic stress. Methods Mouse models This research was regulated under the Animals (Scientific Procedures) Act 1986 Amendment Regulations 2012 following ethical review by the University of Cambridge Animal Welfare and Ethical Review Body (AWERB) for GDF15 fl/fl LysM-Cremice, the Animal Care Committee of the Institute of Biomedical Sciences approved all the experimental protocols (University of São Paulo, Brazil, Protocol CEUA N° 2951040723) for PPARy fl/fl LysM-Cremice andUniversity of RomeTor Vergata Animal Care Italian Ministry of Health Committee (protocol 256/2022-PR) for HFD mice. All animal experiments were performed using C57BL/6mice aged 4–6 weeks, maintained in ventilated cages with group housing (4 per cage) on a 12 h light/12 h dark cycle, in a temperature-controlled (22 ± 2 °C) facility, with ad libitum access to food andwater. During the experimental protocol, all mice were fed either ad libitum or fasted as stated prior to some tests. Mice were sacrificed by cervical dislocation following a 16-hour fasting period. To limit hormonal variability, all experiments were conducted on male mice. Starting at the age of 6weeks,micewere fed either a control chow (R105-25, Safe Diets: Kcal/Kg: 3,45; 22,4% Kcal from proteins; 64,3% kcal from carbohydrates; 13,3% kcal from lipids) or a 60% high fat diet (D12492i, Research Diets: Kcal/gm 5,24; 20% Kcal from proteins; 20% kcal from carbohydrates; 60% kcal from lipids) for a period of 10–14 weeks. Myeloid-specific deletions of Gdf15 were generated using two separate strategies: (1) GDF15 fl/fl (GDF15 Tm1c) mice were bred to transgenic Lyz2Cre/þ mice carrying Cre recombinase under the con- trol of the LysMpromoter (kindly provided by T.Vidal-Puig, Institute of Metabolic Sciences, Univ of Cambridge). This line wasmaintained on a mixed C57BL6N/J background. Genotypes were confirmed by PCR for the presence of Cre and for the detection of the floxed GDF15 allele. GDF15fl/fl and GDF15 fl/fl LysM-Cre were used in this study. (2) C57BL/ 6N wild type mice (4–6 weeks of age) were purchased from Charles River, Italy and were allowed to acclimatize for 1 week before being irradiated. Mice received two separate doses of 5.5 Gy of radiation using Caesium 60 source with a 4 h gap in between doses. 1–2 h post- irradiation, donor bone marrow cells (10 million/mouse) from either male GDF15KOTm1a orwild type littermateswere injected into the tail veins of the irradiated mice. The mice were then housed under stan- dard conditions for onemonth,monitored andweighed regularly until 12 weeks of age prior to experimentation. Myeloid-specific deletions of PPARγ were produced by crossing C57BL/6N PPARy Lox/P mice (The Jackson Laboratory, stock no. 004584) with LysMCre mice (The Jackson Laboratory, stock no. 004781). The genotypes identificationwas performedusing primers and protocols indicated by the manufacturer. Mice that did not express the Cre recombinase enzyme were used as control littermates. The primers used for genotyping the mice were: PPAR-γ lox (wild-type with 200bp band and mutant with band at 230bp): sense oIMR1934 (5′-TGTAATG- GAAGGGCAAAAGG-3′) antisense oIMR1935 (5′-TGGCTTCCAGTGCA- TAAGTT-3′); CRE enzyme (internal control with product at 324bp and mutant with band at 225bp:22 sense (IMR7338) CTAGGCCACA- GAATTGAAAGATCT (antisense oIMR7339) GTAGGTGGAAATTCTAG- CATCATCC; Sense ADIPOQ-CRE 225BP ACGCGTTAATGGCTAATCGC antisense GGCGCTAAGGATGACTCTGG. ADIPOQ-CRE 202BP: anti-sense CGAACGCACTGATTTCGACCAACCAGCGTTTTCGTCTGC. Type 2 diabetic C57B6.BKS(D)-Leprdb/J (db/db) mice were pur- chased from Jackson laboratories and sacrificed at 8 weeks of age. C57BL/6 J were used as controls. Biochemical parameters and tissue harvesting were conducted 16 h after fasting. Two or threemice were kept in each cage and both, water and the diet, were supplied ad libitum. For BAT macrophage transplantation 10-weeks C57BL/6 J HFD and ND-fed male mice were fasted overnight and CD45+/F4.80+ macrophages were isolated from BAT using magnetic cell sorting. Following isolation, 105 cells were injected into the intrascapular region of recipient 6 weeks-old male mice. Two hours after the transplantation, themice were exposed to cold for 20 h. At the end of the cold exposure, the BAT was harvested and frozen for subsequent analysis. RNAseq data analysis All the RNAseq data were processed using the same pipeline with the same parameters. The bulk RNAseq data were quality-checked using FastQC v0.11.9 and subsequently aligned to the gencode mouse reference genome (mm10) using Hisat2 v2.2.137 with default para- meters. The number of reads for all RefSeq genes were counted using featureCounts v2.0.338 using the multi-mapping option. The resulting count matrix was analyzed in R using the DESeq2 package v1.40.239 and Differential Expression was determined using a cutoff significance of padj< 0.05. TheGene Set ExpressionAnalysis was made with clusterProfiler v4.8.3. All the heatmaps were plotted using ComplexHeatmap40,41. The linear models were calculated using the lm function in R base v4.4.0. The pathway scores were computed on bulk RNAseq data and were used to calculate the linear models. Single-cell RNAseq data analysis All the scRNAseqdatawere processed using the same pipelinewith the same parameters. 10x Genomics CellRanger v7.1.042 was used to per- form sample demultiplexing, filtering and UMI counting, using the reference transcriptome of Mus musculus (mm10) with default para- meters. Low quality cells were removed by filtering on the percentage of mitochondrial reads ( <10%), number of UMIs ( >250 and <30000) and number of detected features ( >500). The scRNAseq analysis of good quality cells was performed in R v4.4.0 using Seurat Package v5.0.343. After normalization, we identified the top 2000 variable genes for each library and the libraries were integrated using the CCA method. For all the integrated objects, we Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 11 www.nature.com/naturecommunications performed lineardimensional reduction (PCA), cell clustering anddata visualization using UMAP and t-SNE. Differentially expressed genes that define each cluster were identified using aWilcoxonRank Sum test in Seurat and, accompanied by well-known literature gene markers, the DEGs were used to anno- tate each cell cluster. Bar plots representing cell proportions were plotted using dittoSeq package44 and the Heatmaps representing gene markers expression were plotted using ComplexHeatmap45. Velocity analysis was conducted using two distinct steps. The computation of velocities employedVelocyto v.0.17.1746 with BAM files provided as input to generate unspliced and spliced abundance counts in loom format. For the inference of gene-specific RNA velocities, we employed scVelo v0.2.5 to explore cell trajectories. The loomfileswere processed using a likelihood-based dynamical model, and subse- quently, the velocity vectors were projected onto the UMAP projection. Velocity ranked genes were identified in each cluster using the scvelo.tl.rank_velocity_genes function. This function calculates, on velocity expression, the genes that show dynamics that is tran- scriptomically regulated differently compared to all other clusters. The velocity ranked genes from the monocytes cluster were used to query the STRING database47 to obtain a protein-protein interaction network. STRING is a database of protein interactions that are derived from a variety of sources, including experimental data, literature curation, and textmining. Theprotein-protein interactionnetworkwas then used to identify hub genes using the cytoHubba48 plugin for Cytoscape49. CytoHubba is a collection of algorithms for identifying important nodes in complex networks. The Degrees algorithm was used, which identifies hub genes based on their degree, or the number of connections they have to other genes in the network. FoamSpotteR module FoamSpotteR module is based on a Random Forest (RF) classifier model. RF is an ensemblemachine learning technique based on a set of classification trees that are trained onpartly independent data subsets. For classification problems, the output is obtained at runtime by majority voting of the trees50. From CD45+ cell dataset of atheromatous plaque (GSE), we sub- clustered monocytes and macrophages, resulting in a total of 4976 cells. We annotated these cells as fMAC+ and fMAC- based on well- established literature markers, yielding 3914 cells as fMAC− and 1062 cells as fMAC+. For model training, we utilized 75% of the cells, while the remaining 25% were reserved as test set. To reduce the dimen- sionality of the data, we employed the Boruta algorithm51 for feature selection. Boruta is a feature selection algorithm in machine learning that identifies and selects the most important features from a dataset by comparing them to random noise features, helping to improve model performance by focusing on relevant input variables. The top 20 genes were selected based on their importance and subsequently used to train the model (Supplementary Data). To create and train the RF, we used the R package randomForest (https://CRAN.R-project.org/doc/Rnews/), which implements the ori- ginal Fortran code in the R environment. The RF was tuned with 500 trees using 4 randomly selected predictors at each split. Furthermore, the FoamSpotteR module provided an additional feature, the Foam- DEX Score, which quantifies the probability of a cell developing foam- like features. Values upper 0.5 refer to fMAC. The final RF showed a high level of accuracy, i.e., K = 0.9879. MacPolarizeR module MacPolarizeR is a module designed to categorize macrophages based on their polarization status. MacPolarizeR operates by clustering cells according to the expression of genes associated in M1 and M2 activa- tion status. To identify genes implicated in the inflammatory process, we interrogated bulk RNAseq of peritoneal macrophages stimulated with IFN-γ/LPS or IL4 to obtain M1 or M2-like macrophages, respec- tively (GSE129253). Through a differential gene expression analysis between M1 andM2 cells, we selected the top 75 genes upregulated in pro-inflammatory (M1-like) macrophages and the top 75 genes upre- gulated in anti-inflammatory (M2-like) macrophages, representing the markers of inflammatory and anti-inflammatory macrophages. In order to minimize the clustering error, MacPolarizeR was designed as a method based on the k-means algorithm, a point-based clustering that starts with the cluster centers initially placed at arbi- trary positions and proceeds by moving at each step the cluster centers52. The k-means algorithm was set with 500 iterations and 3 centers corresponding to the classes of macrophages to be clustered: Inflammatory, Transitional and Healing. Additionally, MacPolarizeR calculates an M1 and M2 polarization score for each cell based on the expression of the previously selected genes. The calculated scores can be visualized on a two-dimensional plot, with M1 score on the X-axis and M2 score on the Y-axis (Fig. 3A). This graphical representation provides a clearer interpretation of macrophage behavior, indicating whether subpopulations lean towards a healing phenotype (second quadrant), a transitional state (first and third quadrants), or an inflammatory state (fourth quadrant). All plots weremade by using the ggplot2 package (https://ggplot2.tidyverse.org). Pathway score and PathAnalyzeR module To quantify pathway expression in single-cell data, we developed a novel pathway score that represents the weighted average of the expression of genes in the pathways. This calculation can be per- formed for both individual cells and cell clusters. To allow for com- parison across multiple conditions, this calculation considers only normalized gene expressions. We excluded pathways with at least 5 genes expressed in the dataset and pathways with at least 3 genes expressed in high value (<0.001) in j conditions. For the i-th pathway gene, we first calculated its mean expression level across cell of the j-th condition: Mi, j = Pnj k = 1gi, k nk ð1Þ In which nk is the number of cells in the j-th condition, gi,k is the expression level of the i-th gene in the k-th cell in the condition. The relative expression level of the i-th gene in the j-th condition was then defined as the ratio of Mi,j to its average over all cell types. For condition: ri, j = Mi, j 1 N PN j Mi, j ð2Þ For single cell: ri, j = gi, j 1 N PN j gi, j ð3Þ Where ri,j quantifies the relative expression level of gene i in condition j compared to the average expression level of this gene in all conditions. The weight for each gene is defined as the variance of gene expression: wi = 1 σi ð4Þ The pathway activity score for the t-th pathway and the j-th con- dition was then defined as the weighted average of ri,j over all genes Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 12 https://CRAN.R-project.org/doc/Rnews/ https://ggplot2.tidyverse.org www.nature.com/naturecommunications included in this pathway: Pt, j = Pm i Wi × ri, j Pm i Wi ð5Þ The final score is a positive real number. Scores greater than 1 indicate upregulation of the pathway, while scores less than 1 indicate downregulation of the pathway. qPCR expression analysis Total RNA was extracted using TRI Reagent (Sigma-Aldrich). RNA (3mg) was retro-transcribed by usingM-MLV (Promega, Madison,WI). qPCR was performed in triplicate by using validated qPCR primers (BLAST), Applied Biosystems Power SYBR Green Master Mix, and the QuantStudio3 Real-Time PCR System (Thermo Fisher, Whal- tam, MA, USA). mRNA levels were normalized to beta-Actin mRNA. The reaction was performed according to the manufacturer’s protocol using QuantStudio 3 Real-Time PCR System. Data were ana- lyzed following the 2-DDCt method. The primers used for reverse transcription quantitative PCR are supplied in Supplementary Data Files. Histochemical analysis Formalin-fixed paraffin-embedded BAT explants were cut into 3mm sections and stained with hematoxylin and eosin (H&E) prior to microscope analysis. After antigen retrieval with citrate buffer (pH 6.0), sections were incubated at room temperature with the following primary antibody anti-TREM2. Negative control was obtained by omitting primary antibodies. Immunohistochemical reactions were visualized by DAB as the chromogen from MACH 1 Universal HRP- Polymer Detection (Biocare Medical, Concord, MA, USA). Primary cells isolation Bone marrow (BM) was extracted from the limbs of 8-week-old male mice by perfusion with PBS and 1% P/S. Bone marrow derived macro- phages (BMDM) were plated at a density of 33,000 cells/mL in alpha- MEM supplemented with 10% FBS, 1% P/S, and 1% GlutaMAX in a humidified incubator containing 5% CO2 at 37 C. Macrophage differ- entiation was induced by adding M-CSF (20ng of cells/mL) in the culture medium for 8 days. BMDMs were treated with rosiglitazone 10μM (Sigma Aldrich) for 4, 8 and 16 h. For PPARγ inhibition, BMDMs were treated with GW9662 10μM for 16 h. After the treatment cells were collected for future analysis. The SVF was isolated from BAT according to the manufactures protocol of Miltenyi. Tissues were minced finely with scissors and incubated in high-glucose DMEM containing enzymes from adipose tissue dissociation kit (Miltenyi). Dissociation was carried on using gentleMACS Octo Dissociator with Heaters. The cell suspension was filtered through a 100 µm nylon mesh. Cells were collected in a 50mL conical tube and centrifuged at 300 × g for 10min at 4 °C. The super- natants and floating adipocytes were aspirated, and the SVF pellet was resuspended in 1mL of ACK RBC lysis buffer (Gibco) and incubated at RT for 2min. The RBC lysis reaction was quenched by adding 10mL of cold wash media (high-glucose DMEM containing 5% heat-inactivated FBS, L-glutamine, and penicillin/streptomycin). The cells were cen- trifuged at 300 × g for 10min at 4 °C. SVF fraction was used to isolate CD45+/F4.80+ cells. Pellet was resuspendedwith 90ul of PBSbuffer (PBS, 0.5%ofBSA, 2mMof EDTA) and 10 ul of anti-CD45 microbeads (Miltenyi) and incubated for 15min at 4 °C. After the incubation, 2ml of PBS buffer were added and then centrifuges at 300 g for 5min at 4 °C. Pellet was resuspended with 0.5mlof PBSbuffer and thenpassed throughMACS column (Miltenyi), on MACS magnetic separator (Miltenyi). CD45+ cells were retained on the magnetic separator and the negative passed through the column. Immediately CD45+ cells were flushed out by firmly pushing the plunger into the column. Cellswere then centrifuged 300 g for 5min at 4 °C. Pellet was resuspended in 90ul of PBS buffer and all the proce- dure was repeated using anti-F4/80 microbeads (Miltenyi). CD45+/ F4.80+ cells obtained were resuspended with DMEM HG containing 10% of FBS, 1% of glutamine and 1% of Pen/strep. F4.80+ were seeded at 20.000 cells/cm2 of density. For PPARγ inhibition, CD45+/F4.80+ cells isolated fromBATHFDmicewere treatedwithGW9662 10μMfor 16 h. After the treatment, cells were collected for future analysis. Serum and media analysis GDF15 levels were quantified in the serum of mice 16 h after fasting. Additionally, GDF15 levels were measured in the media collected from CD45+/F4.80+ cells isolated from LFD and HFD mice, following the manufacturer’s protocol (Fine test EM0087). Immunoblotting Cells were lysed in RIPA buffer (50mMTris-HCl, pH 8.0, 150mMNaCl, 12mM deoxycholic acid, 0.5% Nonidet P-40, and protease and phos- phatase inhibitors). Then, 10–20mg of proteins were loaded on sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS- PAGE) and subjected to immunoblotting. Nitrocellulose membranes were incubated with primary antibodies at a dilution of 1:1000. The membranes were then incubated with the appropriate horseradish peroxidase-conjugated secondary antibodies. Immunoreactive bands were detected using a FluorChem FC3 System (Protein-Simple, San Jose, CA, USA) after incubating the membranes with ECL Select Wes- tern Blotting Detection Reagent (GE Healthcare, Pittsburgh, PA, USA). Densitometric analyses of the immunoreactive bands were performed using FluorChem FC3 Analysis software. The antibodies used for immunoblotting are as follows: GDF15 (Santa Cruz, Cod: sc-515675), CD36 (Invitrogen Cod: PA5-27236), UCP1 (Abcam, Cod: ab10983). Representative immunoblots of at least n = 3 independent experi- ments are shown. Untargeted proteomics Tissues were lysed in an ice-cold RIPA buffer. Samples preparation was done using the in StageTip (iST) method employing the iST Kit from Preomics (cat: P.0.00027). We performed a biological quadruplicate (technical replicate n = 1) for each experimental condition. Samples were separated by HPLC in a single run (without pre-fractionations) and analyzed by LC-MS/ MS. Instruments for LC-MS/MS analysis con- sisted of a NanoLC 1200 coupled via a nano-electrospray ionization source to the quadrupole-based Q Exactive HF benchtop mass spec- trometer (ThermoFisher)53. Peptide separation was carried out according to their hydrophobicity on a home-made column, 75mm ID, 8 Um tip, 400mm bed packed with Reprosil-PUR, C18-AQ, 1.9mm particle size, 120 Angstrompore size (NewObjective, Inc., cat. PF7508- 250H363), using a binary buffer system consisting of solution A: 0.1% formic acid and B: 80% acetonitrile, 0.1% formic acid. Total flow rate: 300nL/min. LC gradient, after sample loading, 7–32% buffer B in 45min, followed by 32-45% buffer B in 5min, 45–95% buffer B in 3min, 95–5% buffer B in 5min. MS data acquisition was performed in DIA (Data Independent Acquisition) mode using 32 variable windows covering a mass range of 300–1650m/z. The resolution was set to 60,000 forMSI and 30,000 forMS2. TheAGCwas 3e6 inbothMS1 and MS2,with amaximum injection timeof 60ms inMS1 and 54ms inMS2. NCE were set 25%, 27.5%, 30%. All acquiredRAWfileswereprocessedusing Spectronaut software (Biognosys). For protein assignment, spectra were correlated with Uniprot Mouse database (UP00000589_2024). Searches were per- formed with tryptic specifications and default settings for mass tol- erances for both MS and MS/MS spectra. The other parameters were set as follow: Fixed modifications: Carbamidomethyl (C) Variable modifications: Oxidation, Acetyl (N-term) Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 13 www.nature.com/naturecommunications Digestion: Trypsin, Lys-C Min. peptide length = 7Da Max peptide length = 470Da False discovery rate for proteins and peptide-spectrum = 1% The Perseus software (1.6.2.3) was used to logarithmize, group and filter the protein abundance. Anova/T-test analysis were per- formed setting FDR =0.05. Proteins with Log2 Difference>1 and q-value < 0.05 were considered significantly enriched. Confocal microscopy and morphometric analysis For immunofluorescence and multiple-labeling experiments, iBAT depots form 3 mice per experimental group were analysed. Paraffin sections (3 µm thick) were dewaxed, washed in 0.1% PBS-Tween and incubated with an antigen retrieval solution (Histo-VT- One, Nacalai Tesque, Kyoto, Japan) for 40min at 70 °C. After a thorough rinse in 0.1% PBS-Tween, they were incubatedwith blocking solution (Blocking OneHisto, Nacalai Tesque) at RT for 40min.Next, theywere incubated with a mixture of two primary antibodies diluted in 1% bovine serum albumin (BSA)-PBS (Rabbit polyclonal anti-Mannose Receptor CD206, ab64693, Abcam, Cambridge, United Kingdom and Mouse mono- clonal anti-GDF15, sc-515675, Santa Cruz Biotechnology, CA, USA) and incubated overnight at 4 °C. The next day sections were washed twice with 0.1% PBS-Tween and incubated in a cocktail of fluorophore-linked secondary antibodies (Anti-rabbit IgG Alexa Fluor-555, A31572 and Anti-mouse IgG Alexa Fluor-488, A21202, both from Invitrogen, Thermo Fisher Scientific, MA, USA) diluted 1:400 in 1% BSA-PBS for 30min at RT. Sections were subsequently washed twicewith 0.1% PBS- Tween, stained with TO-PRO3 (T3605, Invitrogen), mounted on stan- dard glass slides, air dried and coverslipped using Vectashield mounting medium (Vector Laboratories, Newark, California, USA). Sections were viewed under a motorized DM6000 microscope (Leica Microsystems, Wetzlar, Germany) at different magnifications. Fluor- escence was detected with a Leica TCS-SL spectral confocal micro- scope equipped with an Argon and He/Ne mixed gas laser. Fluorophores were excited with the 488, 543 and 649 nm lines and imaged separately. Images (1024 × 1024 pixels) were obtained sequentially from three channels using a confocal pinhole of 1.1200 and stored as TIFF files. The brightness and contrast of the final images were adjusted using Photoshop 6 (Adobe Systems). For morphometric analysis, the percentage of CD206-positive interstitial macrophages also positive for GDF15 were calculated from double-stained representative sections from 3 mice per experimental group. From each section, 10 non-overlapping images were randomly selected using a x63 objective. Multiplex immunofluorescence staining The 5μmsections of FFPEmousebrownadipose tissuewere stainedby multiplex immunohistochemistry (mIHC) using the Opal™ 7-Color IHC Kit (PerkinElmer, Waltham, USA, Cat. No. NEL820001KT) on the Leica BOND RX automated immunostainer (Leica Microsystems, Milton Keynes, UK). Thepanelwas generated using four primary antibodies to identify cell populations positive for UCP1, ADIPONECTIN, TREM2 and CD9. Cell nuclei were counterstained with Spectral DAPI. Detailed information on the antibodies, including clone, dilution, Opal paring and retrieval buffers, is provided in the followingpanel. Anti-rabbit and anti-mouse secondary antibodies are included in the Opal™ 7-Colour IHC Kit. Instead, a goat secondary antibody from Abcam (goat HRP, Abcam AB6741, 1:2000 dilution) was used for the goat anti-TREM2 primary. Multispectral imageswere acquired using the PhenoImagerHT2.0 imaging system (Akoya Biosciences) at 20× magnification. Spectral unmixingwasperformed automatically by the PhenoImager andwhole slide unmixed QPTIFFs were generated. Antibody panel available as supplementary table. Flow cytometry analysis Single cell suspension obtained from BAT of LFD and HFD mice was incubated for 15min at 4 °C with a cocktail mix containing anti-mouse CD45, CD11b, F4.80, CD9, Trem2 and Acqua Zombie live dead dye. After staining, cells were washed in cold running buffer (PBS 1X, 2.5% FBS, 2.5% EDTA) and incubated for 20min at RT in True Nuclear Fix solution (Biolegend). Cells were washed in cold running buffer and incubated at RT with Perm solution (Biolegend) containing primary antibodies against PPARg (rabbit) and GDF15 (mouse). After 30min, cells were washed in PBS and incubated at RT with Perm solution (Biolegend) containing their respective secondary antibodies (anti- rabbit Alexa405 and anti-mouse Alexa488). Macrophages were iden- tified as CD45+/CD11b+/F4.80+ cells and inside this gate, LAMs were identified as CD9+ Trem2+ cells. Intracellular expression of Pparg and Gdf15 was evaluated in LAMs. Samples were acquired on a 13-color Cytoflex Flow Cytometer (Beckman Coulter) and for each analysis at least 0.5 × 106 live cells were acquired by gating on viable cells using Acqua Zombie (1:300) and analyzed by FlowJo Software (TreeStar, v10). Mitochondrial and lysosomal mass was also evaluated in CD45+ F4.80+ cells of BAT tissues following staining respectively with Mito- Tracker red or LysoTracker red and analyzing cells by means of flow cytometry. Quantification and statistical analysis Data were expressed as the mean± SD unless otherwise stated. The exact numbers of replicates are given in each figure legend. A two- tailed unpaired Student’s t-test was performed to assess the statistical significance between two groups. Analysis of variance (ANOVA) fol- lowed by Dunnett’s (comparisons relative to controls) or Tukey’s (multiple comparisons among groups) post hoc tests was used to compare threeormoregroups. Venndiagramswereconstructed using Draw Venn Diagram or Venny 2.1.0 software. Statistical analyses were performed using GraphPad Prism 9 (GraphPad Software Inc., San Diego, CA, USA). In all cases, a p-value or FDR less than 0.05 was set as the minimum significance threshold. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Data availability All data supporting the findings described in this manuscript are available in the article and in the Supplementary Information and from the corresponding author upon request. Source data are also provided with this paper. The sequencing data were obtained from publicly available sources. The main analysis was performed on GSE232278, scRNAseq of BAT db/db and HFD. The bulk RNAseq analysis was per- formed on GSE112740 BAT HFD and GSE232276 BAT db/db. The FoamSpotteR random forest was trained onGSE97310 andGSE116240, scRNAseq of CD45+ cells of atherosclerotic aorta. GSE116239 bulk RNAseq of foamy and non-foamy macrophages isolated from athero- sclerotic aortas via lipid probe-based flow cytometry. The polarization genes were selected using GSE129253, bulk RNAseq of in-vitro polar- ized macrophages. Validation of MacPolarizeR employed GSE117176, scRNAseq of in-vitro polarized macrophages, and PRJNA626597, scRNAseq of septic adipose tissue. Validation of the role of Pparg in macrophages was performed using GSE111105, RNAseq of in vitro PpargKOmacrophages treatedwithRosi and IL4. Comparative analysis was performed using GSE182233 scRNAseq of HFD WAT. Proteomics data of BAT of PPARy fl/fl LysM-Cre mice are available in the Supple- mentary Information. Source data are provided with this paper. Code availability The R package code and full tutorials are available on GitHub at andreeedna/MACanalyzeR [https://github.com/andreeedna/MACanal Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 14 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE232278 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE112740 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE232276 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE97310 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116240 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116239 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129253 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117176 https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA626597 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111105 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE182233 https://github.com/andreeedna/MACanalyzeR www.nature.com/naturecommunications yzeR]. The code is available on GitHub at andreeedna/MACanaly- zeR_Methods [https://github.com/andreeedna/MACanalyzeR_Methods]. References 1. Villarroya, F., Cereijo, R., Villarroya, J., Gavalda-Navarro, A. & Giralt, M. Toward an Understanding of How Immune Cells Control Brown and Beige Adipobiology. Cell Metab. 27, 954–961 (2018). 2. Cereijo, R. et al. CXCL14, a Brown Adipokine that Mediates Brown- Fat-to-Macrophage Communication in Thermogenic Adaptation. Cell Metab. 28, 750–763.e756 (2018). 3. Aquilano, K., Zhou, B., Brestoff, J. R. & Lettieri-Barbato, D. Multi- faceted mitochondrial quality control in brown adipose tissue. Trends Cell Biol. 33, 517–529 (2023). 4. Rosina, M. et al. Ejection of damaged mitochondria and their removal by macrophages ensure efficient thermogenesis in brown adipose tissue. Cell Metab. 34, 533–548.e512 (2022). 5. Burl, R. B., Rondini, E. A., Wei, H., Pique-Regi, R. & Granneman, J. G. Deconstructing cold-inducedbrown adipocyte neogenesis inmice. Elife 11, https://doi.org/10.7554/eLife.80167 (2022). 6. Kriebs, A. Accumulationofmacrophages in adipose tissue.Nat.Rev. Endocrinol. 17, 4 (2021). 7. Jaitin, D. A. et al. Lipid-Associated Macrophages Control Metabolic Homeostasis in a Trem2-DependentManner.Cell 178, 686–698.e614 (2019). 8. Liu, C. et al. TREM2 regulates obesity-induced insulin resistance via adipose tissue remodeling inmiceof high-fat feeding. J. Transl.Med 17, 300 (2019). 9. Coats, B. R. et al. Metabolically Activated Adipose Tissue Macro- phages Perform Detrimental and Beneficial Functions during Diet- Induced Obesity. Cell Rep. 20, 3149–3161 (2017). 10. Xu, X. et al. Obesity activates a program of lysosomal-dependent lipid metabolism in adipose tissue macrophages independently of classic activation. Cell Metab. 18, 816–830 (2013). 11. Rabhi, N. et al. Obesity-induced senescent macrophages activate a fibrotic transcriptional program in adipocyte progenitors. Life Sci. Alliance 5, https://doi.org/10.26508/lsa.202101286 (2022). 12. Sawaki, D. et al. Osteopontin promotes age-related adipose tissue remodeling through senescence-associated macrophage dys- function. JCI Insight 8, https://doi.org/10.1172/jci.insight.145811 (2023). 13. Denisenko, E. et al. Systematic assessment of tissue dissociation and storage biases in single-cell and single-nucleus RNA-seq workflows. Genome Biol. 21, 130 (2020). 14. Duerre, D. J. & Galmozzi, A. Deconstructing Adipose Tissue Het- erogeneity One Cell at a Time. Front Endocrinol. (Lausanne) 13, 847291 (2022). 15. Fernandez, D. M. & Giannarelli, C. Immune cell profiling in athero- sclerosis: role in research and precisionmedicine.Nat. Rev. Cardiol. 19, 43–58 (2022). 16. Cochain, C. et al. Single-Cell RNA-Seq Reveals the Transcriptional Landscape and Heterogeneity of Aortic Macrophages in Murine Atherosclerosis. Circ. Res 122, 1661–1674 (2018). 17. Kim, K. et al. TranscriptomeAnalysis RevealsNonfoamyRather Than Foamy Plaque Macrophages Are Proinflammatory in Athero- sclerotic Murine Models. Circ. Res 123, 1127–1142 (2018). 18. Wculek, S. K., Dunphy, G., Heras-Murillo, I., Mastrangelo, A. & Sancho, D. Metabolism of tissue macrophages in homeostasis and pathology. Cell Mol. Immunol. 19, 384–408 (2022). 19. Charo, I. F. Macrophage polarization and insulin resistance: PPAR- gamma in control. Cell Metab. 6, 96–98 (2007). 20. Zhang, L. & Chawla, A. Role of PPARgamma in macrophage biology and atherosclerosis. Trends Endocrinol. Metab. 15, 500–505 (2004). 21. Jin, S. et al. Inference and analysis of cell-cell communication using CellChat. Nat. Commun. 12, 1088 (2021). 22. Patsalos, A. et al. A growth factor-expressing macrophage sub- population orchestrates regenerative inflammation via GDF-15. J. Exp. Med. 219, https://doi.org/10.1084/jem.20210420 (2022). 23. Jena, J. et al. GDF15 is required for cold-induced thermogenesis and contributes to improved systemic metabolic health following loss ofOPA1 in brown adipocytes. Elife 12, https://doi.org/10.7554/eLife. 86452 (2023). 24. Igual-Gil, C. et al. GDF15 is required for maintaining subcutaneous adipose tissue lipidmetabolic signature.Sci. Rep. 14, 26989 (2024). 25. Patel, S. et al. Combined genetic deletion of GDF15 and FGF21 has modest effects on body weight, hepatic steatosis and insulin resistance in high fat fed mice. Mol. Metab. 65, 101589 (2022). 26. Sciarretta, F. et al. Lipid-associated macrophages reshape BAT cell identity in obesity. Cell Rep. 43, 114447 (2024). 27. Piollet, M. et al. TREM2 protects from atherosclerosis by limiting necrotic core formation.Nat. Cardiovasc. Res. 3, https://doi.org/10. 1038/s44161-024-00429-9 (2024). 28. Zhang, K. et al. TREM2(hi) resident macrophages protect the septic heart by maintaining cardiomyocyte homeostasis. Nat. Metab. 5, 129–146 (2023). 29. Haidar, M. et al. Targeting lipophagy in macrophages improves repair in multiple sclerosis. Autophagy 18, 2697–2710 (2022). 30. Flaherty, S. E. et al. A lipase-independent pathway of lipid release and immune modulation by adipocytes. Science 363, 989–993 (2019). 31. Itoh, T. et al. Structural basis for the activation of PPARgamma by oxidized fatty acids. Nat. Struct. Mol. Biol. 15, 924–931 (2008). 32. Odegaard, J. I. et al. Alternative M2 activation of Kupffer cells by PPARdelta ameliorates obesity-induced insulin resistance. Cell Metab. 7, 496–507 (2008). 33. Varga, T. et al. Macrophage PPARgamma, a Lipid Activated Tran- scription Factor Controls the Growth Factor GDF3 and Skeletal Muscle Regeneration. Immunity 45, 1038–1051 (2016). 34. Hevener, A. L. et al. Macrophage PPAR gamma is required for nor- mal skeletal muscle and hepatic insulin sensitivity and full anti- diabetic effects of thiazolidinediones. J. Clin. Invest 117, 1658–1669 (2007). 35. Odegaard, J. I. et al. Macrophage-specific PPARgamma controls alternative activation and improves insulin resistance. Nature 447, 1116–1120 (2007). 36. Wang, D. et al. GDF15 promotes weight loss by enhancing energy expenditure in muscle. Nature 619, 143–150 (2023). 37. Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph- based genome alignment and genotyping with HISAT2 and HISAT- genotype. Nat. Biotechnol. 37, 907–915 (2019). 38. Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic fea- tures. Bioinformatics 30, 923–930 (2014). 39. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold changeanddispersion for RNA-seqdatawithDESeq2.GenomeBiol. 15, 550 (2014). 40. Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov. (Camb.) 2, 100141 (2021). 41. Yu, G., Wang, L. G., Han, Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287 (2012). 42. Zheng, G. X. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017). 43. Hao, Y. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. 42, 293–304 (2024). 44. Bunis, D. G., Andrews, J., Fragiadakis, G. K., Burt, T. D. & Sirota, M. dittoSeq: universal user-friendly single-cell and bulk RNA sequen- cing visualization toolkit. Bioinformatics 36, 5535–5536 (2021). 45. Gu, Z. Complex heatmap visualization. iMeta 1, e43 (2022). Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 15 https://github.com/andreeedna/MACanalyzeR https://github.com/andreeedna/MACanalyzeR_Methods https://doi.org/10.7554/eLife.80167 https://doi.org/10.26508/lsa.202101286 https://doi.org/10.1172/jci.insight.145811 https://doi.org/10.1084/jem.20210420 https://doi.org/10.7554/eLife.86452 https://doi.org/10.7554/eLife.86452 https://doi.org/10.1038/s44161-024-00429-9 https://doi.org/10.1038/s44161-024-00429-9 www.nature.com/naturecommunications 46. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018). 47. Szklarczyk, D. et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res 49, D605–D612 (2021). 48. Chin, C. H. et al. cytoHubba: identifying hub objects and sub- networks from complex interactome. BMC Syst. Biol. 8, S11 (2014). 49. Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13, 2498–2504 (2003). 50. Breiman, L. Random Forests. Mach. Learn. 45, 5–32 (2001). 51. Kursa, M. B. & Rudnicki, W. R. Feature Selection with the Boruta Package. J. Stat. Softw. 36, 1–13 (2010). 52. Kawamura, A. et al. Correction to: Comparison of the usability of an automatic sleep staging program via portable 1-channel electro- encephalograph and manual sleep staging with traditional poly- somnography. Sleep. Biol. Rhythms 21, 125–126 (2023). 53. Michalski, A. et al. Mass spectrometry-based proteomics using Q Exactive, a high-performance benchtop quadrupole Orbitrap mass spectrometer. Mol. Cell Proteom. 10, M111 011015 (2011). Acknowledgements This work was mainly supported by the European Union - Next Gen- eration EU, call “PRIN 2022 (Prot. 2022Z9HYJH)” and “PRIN 2022 PNRR (Prot. P2022LS3FF)” to D.L.B. This work was partially supported by the MUR-PNRR M4C21.3 PE6 project PE00000019 Heal Italia to K.A. This work was also supported by Grants 2019/09679-2 and 2022/06073-9, São Paulo Research Foundation (FAPESP) to J.C.R.-N. D.B.S. is sup- ported by the Wellcome Trust (WT 219417), the MRC (MRC Reference: MR/X00970X/1), and The National Institute for Health Research (NIHR) Cambridge Biomedical Research Centre and NIHR Rare Disease Translational Research Collaboration. Figures 1D, 2A, 5B, K, 6P and S3(A, J) were created in BioRender. Lettieri barbato, D. (2025) https:// BioRender.com/wy626mu“. Author contributions D.L.-B., A.N. and F.Z. conceived the study, designed the analyses and interpreted the results. A.N. developed the RpackageMACanalyzeR. A.N., F.Z. and L.V. tested MACanalyzeR package. A.N., F.Z. and L.V. performed bioinformatic analyses. F.S. processed samples for the experimental part and contributed to formal analysis. L.S.S and J.C.R.-N. generated and performed experiment on PPARγ LysMCre mice. S.P. and D.B.S. per- formed experiment onWT andGdf15KO bonemarrow transplantedmice. S.C., L.N., and S.B. performed immunohistochemistry analyses. A.S., D.L., F.C. and S.B. performed multiplex immunofluorescence staining. I.S., M.S. and A.G. performed confocal microscopy and morphometric analysis. C.D.B. and F.T. contributed to formal analysis. P.G. performed proteomics analysis. G.C. and V.C. performed and analyzed high- dimensional flow cytometry. A.N., F.Z., F.S., K.A. and D-L.-B. wrote the manuscript and edited the figures. B.Z, J.W.W, S.I., G.R, C.C, V.C. con- tributed to the text. D.L.-B. funded the study. Competing interests The authors declare no competing interest. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-60295-2. Correspondence and requests for materials should be addressed to Daniele Lettieri-Barbato. Peer review information Nature Communications thanks the anon- ymous reviewer(s) for their contribution to thepeer reviewof thiswork. A peer review file is available. Reprints and permissions information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jur- isdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. Youdonot havepermissionunder this licence toshare adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by-nc-nd/4.0/. © The Author(s) 2025 1Department of Biology, University of Rome Tor Vergata, Rome, Italy. 2PhD Program in Evolutionary Biology and Ecology, Department of Biology, University of Rome Tor Vergata, Rome, Italy. 3Immunometabolism Research Group, Department of Cell Biology and Development, Institute of Biomedical Sciences, University of São Paulo (ICB1-USP), São Paulo, Brazil. 4Microscopic and Ultrastructural Anatomy Research Unit, Department of Medicine and Surgery, Università Campus Bio-Medico di Roma, Rome, Italy. 5Telethon Institute of Genetics and Medicine, Pozzuoli, Italy. 6Department of Clinical Medicine and Surgery, University Federico II, Naples, Italy. 7Metabolic Research Laboratories, Wellcome Trust-Medical Research Council Institute of Metabolic Science, University of Cambridge, Cambridge, UK. 8Laboratory of Resolution of Neuroinflammation, IRCCS Santa Lucia Foundation, Rome, Italy. 9Department of Systems Medicine, University of Rome Tor Vergata, Rome, Italy. 10Multiplex Spatial Profiling Facility, Fondazione Policlinico Universitario ‘Agostino Gemelli’ IRCCS, Rome, Italy. 11Department of Translational Medicine and Surgery, Università Cattolica del Sacro Cuore, Rome, Italy. 12Department of Experimental and Clinical Medicine, Marche Polytechnic University, Ancona, Italy. 13IRCSS INRCA, Ancona, Italy. 14Center of Obesity, Marche Polytechnic University-United Hospitals, Ancona, Italy. 15Institute of Experimental Biomedicine, University Hospital Würzburg, D16, Würzburg, Germany. 16Paris Cardiovascular Research Center, Université Paris Cité, INSERMU970, Paris, France. 17Institute of Translational Pharmacology, National Research Council, Rome, Italy. 18Université Côte d’Azur, CNRS, Nice, France. 19Department of Immunology, School of Medicine, University of Connecticut, Farmington, CT, USA. 20Center for Immunology, University of Minnesota, Minneapolis, USA. 21Department of Integrative Biology and Physiology, University of Minnesota, Minneapolis, USA. e-mail: daniele.lettieri.barbato@uniroma2.it Article https://doi.org/10.1038/s41467-025-60295-2 Nature Communications | (2025) 16:5063 16 https://BioRender.com/wy626mu https://BioRender.com/wy626mu https://doi.org/10.1038/s41467-025-60295-2 http://www.nature.com/reprints http://creativecommons.org/licenses/by-nc-nd/4.0/ http://creativecommons.org/licenses/by-nc-nd/4.0/ mailto:daniele.lettieri.barbato@uniroma2.it www.nature.com/naturecommunications MACanalyzeR scRNAseq analysis tool reveals PPARγHIGH/GDF15HIGH lipid-associated macrophages facilitate thermogenic expansion in BAT Results Lipid-associated macrophages accumulate in BAT of murine models of obesity MACanalyzeR identifies foamy-like macrophages in obese BAT BAT LAMs show a different polarization status between db/db and HFD mice BAT LAMs from HFD mice showed an enhanced fatty acid metabolism Thermogenically-competent LAMs showed an enhanced PPAR pathway activation PpargHIGH LAMs maintain BAT identity through GDF15 signalling Discussion Methods Mouse models RNAseq data analysis Single-cell RNAseq data analysis FoamSpotteR module MacPolarizeR module Pathway score and PathAnalyzeR module qPCR expression analysis Histochemical analysis Primary cells isolation Serum and media analysis Immunoblotting Untargeted proteomics Confocal microscopy and morphometric analysis Multiplex immunofluorescence staining Flow cytometry analysis Quantification and statistical analysis Reporting summary Data availability Code availability References Acknowledgements Author contributions Competing interests Additional information