Different gene expression profiles were obtained from the GEO database of the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/geo/). All RNA sequencing data were listed as below: Brown and white adipose tissues: GSE44138 (cold 1d), GSE118849 (cold 3d), GSE148361(cold 4d); Brown adipose tissues: GSE100924 (cold 6h), GSE133050 (cold 6h), GSE164936 (cold 6h), GSE149124 (cold 6h), GSE119964 (cold 6h), GSE119452 (cold 6h/1d/2d), GSE207705 (cold 6h/1d/2d/3d/4d/5d), GSE181123 (cold 6h), GSE147392 (cold 6h), GSE110055 (cold 6h), GSE135391(cold 1d), GSE70734 (cold 3d), GSE148361 (cold 4d), GSE86338 (cold 7d), GSE144186 (cold 7d), GSE178720 (cold 10d), GSE51080 (cold 10d); White adipose tissues: GSE140259 (cold 2d), GSE179385 (cold 3d), GSE148361(cold 4d), GSE164219 (cold 7d), GSE145498 (cold 7d), GSE13432 (cold 7d/35d), GSE110420 (cold 10d), GSE51080 (cold 10d), GSE183000 (cold 14d).

Two methods were used for evaluating the infiltration of immune cells in expression profiles. One was the Microenvironment Cell Population counter (mMCP-counter). Based on highly specific transcriptomic markers, 16 immune and stromal murine cell populations were quantified to evaluate their infiltration for the control and cold stimulations using the R package “mMCP-counter” (18). Another method was CIBERSORT (19), which estimated the relative proportion of each relevant cell type using an input matrix of reference gene expression signatures. The proportion of 25 different immune cells between cold and control groups was compared by two-sided statistical tests, where the significant level was set at p < 0.05. The Corrplot package was used for the correlation analysis of infiltrating immune cells.

ScRNA-seq data of GSE207706 was downloaded from the GEO database of the NCBI. The gene expression profiles were paired and normalized based on R software (version 6.0.2). Briefly, sequencing data were analyzed using the “Seurat” package following the provided protocol (20). The “SingleR” package and PanglaoDB database were used for cell type annotation (21, 22). Cells were represented with Distributed Stochastic Neighbor Embedding (TSNE) plots (23).

Male C57Bl/6J mice (8-10 weeks) were either exposed to 4°C for designed periods or at room temperature (RT, 22°C) with free access to water and food under controlled conditions of light and darkness (12h/12h). All animal experiments were approved by the Animal Care and Use Committee of Central South University. Adipose tissue samples were collected for RT-qPCR, flow cytometry or immunofluorescence. The in vivo experiments were repeated two times independently.

SVFs of BAT and SAT were isolated from mice either at RT or exposed to 4°C (Cold) (1d/3d/7d) as previously reported (24). Briefly, BAT and SAT were collected, sliced into pieces with scissors, and digested in Collagenase II buffer for 30 min at 37°C. Then the suspension was filtered and centrifuged for separation of SVFs from adipocyte fractions. The SVFs were incubated with red cell lysing reagent for 5 min, centrifuged at 500 g for 10 min at 4°C, and washed with PBS. Then total RNA was extracted using the Magzol reagent (Magen, #R4801-01) and synthesized into complementary DNA using the Reverse transcription kit (Accurate Biology, #AG11728). Real-time PCR was performed with gene-specific primers and mixture (Yeason, #11202ES03) on Applied Biosystems™ 7900HT Fast Real-Time PCR. Primers and gene information are shown in Table S1 .

SVF was isolated from male mice as mentioned above. Single-cell suspension was divided into two parts. One part was stained with zombie-APC-Cy7 dye (1:200, Biolegend, #423105) for 10 min and then washed with cell staining buffer (BD Pharmingen#554657). After incubation with FcX blocking (1:100, Biolegend, #156604) for 10 min, cells were stained with specific antibodies. For the analysis of monocytes, suspended cells were stained with CD45-PE (1:100, Biolegend, #147711), Ly6C-AF488 (1:100, Biolegend, #128021), Ly6G-BV510(1:100, Biolegend, #127633) and CD11B-Pacific blue (1:100, Biolegend, #101223). For the analysis of M1 and M2 macrophages, suspended cells were stained with CD45-PerCP-Cyanine5.5 (1:100, Biolegend, #103131), CD11B-Pacific blue (1:100, Biolegend, #101223), F4/80-PE/Cy7 (1:100, Biolegend, #123114), CD206-PE (1:100, Biolegend, #141705), and CD11C-Percp (1:100, Biolegend, #117325). For the analysis of NK and NKT cells, suspended cells were stained with CD45-PE (1:100, Biolegend, #147711), CD3-APC (1:100, Biolegend, #100235), and NK1.1-PE/Cy7 (1:100, Biolegend, #156513). For the analysis of endothelial cells, suspended cells were stained with CD45-PerCP-Cyanine5.5 (1:100, Biolegend, #103131) and CD31-APC (1:100, Biolegend, #102409). The other suspension was incubated with FcX blocking (1:100, Biolegend, #156604) for 10 min, and then the cells were stained with specific antibodies. For the analysis of CD4 and CD8 subpopulations, suspended cells were stained with CD45-PerCP-Cyanine5.5 (1:100, Biolegend, #103131), CD3-FITC (1:100, Biolegend, #100203), CD4-Pacific blue (1:100, Biolegend, #100427), CD8-AF700 (1:100, Biolegend, #100729), CD44-APC (1:100, Biolegend, #103011), and CD62L-APC-Cy7 (1:100, Biolegend, #104427). After washes with a cell staining buffer, the cells were finally incubated with Propidium Iodide Solution-PE (1:100, Biolegend, #421301) before fluorescence-activated cell sorting. Stained cells were sorted by flow cytometer (Cytek Northern Lights 3000 and Agilent Novocyte Advanteon) and data were analyzed with FlowJo software version 10.9.0.

Fresh adipose tissues were obtained, fixed, and embedded in paraffin. After antigen repair and hydrogen peroxide blocking, the slides were blocked with a solution containing 3% BSA before incubation with primary antibodies at 4°C overnight. The slides were then treated with horseradish peroxidase (HRP)-conjugated anti-rabbit (Servicebio, GB23303, 1:500) or Alexa Fluor 488-conjugated anti-rabbit (Servicebio, GB25303, 1:400) secondary antibodies for 50 min at room temperature. The primary antibodies included CD44 (Servicebio, GB112054, 1:3000), CD62L (Bioss, bs-1036R, 1:1000), CD31 (Servicebio, GB113151, 1:1000), TH (Servicebio, GB11181, 1:1000), CD11B (Service bio, GB11058, 1:3000), CD115 (Servicebio, GB11581, 1:1000), and NK1.1 (Abcam, AB289542, 1:100). All the fluorescence pictures were captured by a fluorescence microscope (Nikon Eclipse C1) or laser scanning confocal microscope (Zeiss LSM 780).

Statistical analysis was performed using GraphPad Prism 10.1.1. The results were expressed as mean ± SEM. The p < 0.05 obtained from two-sided tests implied statistically significant differences.

To determine changes in different immune cells in the BAT or SAT after different periods of cold exposure, we downloaded distinct datasets of BAT and SAT from the GEO database of the NCBI. The dataset analysis strategy is shown in Figure 1 . After batch normalization, all datasets were analyzed with two immune algorithms. Significantly different immune cells (and several stromal cells) are listed in Table 1 .

To screen the potential immune cells in thermogenic adipose tissues of mice under different periods of cold exposure, we adopted two well-recognized immune cell screening algorithms. One was the Cibersort algorithm, which measured the proportion of 25 types of immune cells. The other was the mMCP-counter algorithm, which analyzed the abundance of 16 types of murine immune and stromal cells. Consistent with the thermogenic activation of brown fat in acute cold (25, 26), most datasets showed changes started as early as a cold duration of 6h in BAT. These algorithm analyses revealed that beige adipogenesis might occur during longer periods of cold exposure, as 35d long-term cold acclimation resulted in more significant changes in immune cells in SAT ( Table 1 ). Next, we focused on several major types of immune cells or stromal cells and pooled the same cold duration datasets to reanalyze their proportions or abundances over time using Cibersort and mMCP-counter methods.

The Cibersort algorithm analysis suggested that in BAT, the proportions of naïve CD4⁺ T cells and naïve CD8⁺ T cells increased throughout most periods of cold exposure compared with that at RT ( Figures 2A, B ). Memory CD4⁺ T cells and memory CD8⁺ T cells displayed similar dynamic change trends ( Figures 2C, D ). In SAT, the proportion of memory CD8⁺T cells, M1 and M2 macrophages, and Treg cells significantly changed in a cold duration of 4d compared to other cold periods ( Figures 2D-G ). Besides, alterations of these immune cells (such as M2 macrophages, Treg cells, and γδT cells) occurred later in SAT than those in BAT during different cold periods, and the cellular dynamics of the two types of adipose tissues were similar during the cold periods ( Figures 2F-H ). Moreover, the proportion of NK resting cells under all cold exposure periods was lower than that at RT in SAT ( Figure 2I ). Due to these changes in immune cells, the corrplot package on R studio was further used to analyze the relationship among these immune cells to explore their possible synergistic or antagonistic interactions. In BAT, these immune cells were less closely related ( Figure 2J ). In SAT, there was an inverse correlation between M2 macrophages and naïve CD8⁺T cells ( Figure 2K ).

The mMCP-counter algorithm analysis disclosed that the abundance of T cells, CD8⁺T cells, B-derived cells, and NK cells was lower than that at RT or TN in all periods of cold exposure ( Figures 3A-D ), while the abundance of monocytes increased at all cold periods in BAT ( Figure 3E ). Intriguingly, during most cold periods (from 2d to 7d), the changes in monocytes/macrophages in BAT and SAT were opposite ( Figure 3F ), whereas the abundance changes of fibroblasts in these two types of adipose tissues were similar ( Figure 3G ). In addition, endothelial cells and vessels showed similar patterns in both BAT and SAT ( Figures 3H, I ). In SAT, however, except for monocytes and NK cells, the abundance of other immune cells or stromal cells (such as T cells, CD8⁺T cells, B-derived cells, fibroblasts, endothelial cells, and vessels) was significantly decreased at 3 d and then gradually increased in longer periods of cold exposure ( Figures 3A-C, G-I ). Finally, the relationship between immune cells in BAT and SAT was evaluated. In BAT, they were not correlated ( Figure 3J ), but in SAT, CD8⁺T cells and B-derived cells were highly positively correlated, consistent with the abundance changes in BAT and SAT. In contrast, monocytes and T cells were negatively correlated ( Figure 3K ).

To validate bulk RNA sequencing results, we collected ScRNA-seq datasets (GSE207706) of specific immune cells in BAT upon 4-d cold stimulation. We identified ten clusters within lineage-positive cells ( Figure 4A ) and presented the top 10 marker genes for each cluster in Supplementary Figure S1A . Consistent with those mMCP-counter algorithm results, we found that the proportions of T cells and B cells decreased within lineage-positive cells, while the proportion of monocytes increased ( Figure 4B ). Additionally, there was a significant increase in the proportion of other immune cells such as macrophages, dendritic cells, and erythrocytes, while the proportion of NKT cells, neutrophils, endothelial cells and fibroblasts showed a decrease. Subsequently, we further analyzed the T cell cluster within lineage-positive cells and classified it into ten subclusters ( Figure 4C ), presenting their top markers in Supplementary Figure S1B ; among which naïve CD4⁺T cells increased upon 4d cold exposure, consistent with the above immune cell analysis by Cibersort ( Figure 4D ).

To further verify the dynamic changes in immune cells in adipose tissues, we reared the male mice at RT or in cold chambers for 1d, 3d, and 7d, and then collected SVF of BAT and SAT for flow cytometry, immunofluorescence, and RT-qPCR ( Figure 5A ). The expression of Ucp1, a thermogenic marker gene, in BAT and SAT was upregulated in all cold periods, suggesting the effectiveness of cold exposure ( Figure 5B ). Then flow cytometry was performed to analyze the proportion changes of immune cells. As expected in BAT, the proportion of monocytes was increased in cold periods ( Figure 5C and Supplementary Figure S2A ), whereas the proportion of NK cells ( Figure 5D and Supplementary Figure S2A ) and mRNA expression of NK1.1 (an NK cell marker) ( Figure 5E ) were decreased under all periods of cold stimulation. It is well known that sympathetic neurons play a vital role in the thermogenesis of BAT by releasing norepinephrine, while tyrosine hydroxylase (TH) is an essential marker of sympathetic innervation (27). Vascular endothelial cells provide nutrition and oxygen to brown and beige adipocytes to sustain thermogenesis (28). Given the well-recognized crosstalk of two types of cells with thermogenic adipocytes (28), we further explored the possible localization relationship between these cells and immune cells in BAT. Interestingly, the immunofluorescence intensity of CD31 (an endothelial cell marker) or TH was highly correlated with that of NK1.1 ( Figures 5F, G, J and Supplementary Figures S2B, C ). The proximity between these two proteins probably suggested a robust functional interaction between these two types of cells, whereas monocytes (marked by CD11B⁺CD115⁺) were less correlated with either endothelial cells or sympathetic neurons ( Figures 5H-J and Supplementary Figures S2D, E ) in BAT. Furthermore, the analysis of immunofluorescence showed an increase in the proportion of monocytes, particularly during 3 days of cold exposure ( Figure 5K and Figure S2F ), while the proportion of NK cells decreased across all periods of cold exposure. ( Figure 5L and Supplementary Figure S3A ).

Also in SAT, the proportion of monocytes was highly upregulated under cold exposure ( Supplementary Figure S3B ), nevertheless the proportion of NK cells was increased in 1 and 3 days of cold exposure ( Supplementary Figure S3C, D ). Previous studies have shown that iNKT, a subset of NKT cells, can stimulate the production of fibroblast growth factor 21 (FGF21) to regulate the thermogenic browning of white fat (29). Therefore, we also examined changes in NKT cells and observed a decreased proportion in longer periods of cold exposure (3 and 7 days) in SAT ( Supplementary Figures S3C, D ).

Moreover, flow cytometry manifested that the proportion of naïve T (both CD4⁺ and CD8⁺) cells and memory T (both CD4⁺ and CD8⁺) cells increased in all periods of cold stimulation ( Figure 6A and Supplementary Figure S4A ). However, these types of cells were not such close to CD31⁺ endothelial cells or TH⁺ sympathetic neurons ( Figures 6B-D ). The immunofluorescence results implied that the proportion of naïve T cells was increased during 1, 3, or 7d cold exposure, while there was an increased trend in the proportion of memory T cells after 3 days of cold exposure in BAT ( Figure 6E and Supplementary Figure S4B ). Additionally, unlike BAT, alterations in naïve and memory T cells were mainly observed after 7 days of cold stimulation in SAT ( Supplementary Figure S4C ).

As a predominant subpopulation of immune cells, M1 and M2 macrophages changed dynamically, as their proportion both upregulated in BAT in 1d acute cold and gradually downregulated in 3-7d cold durations ( Figure 7A and Supplementary Figure S5A ). In SAT, the proportion of M1 macrophages was reduced in all cold periods, while the proportion of M2 macrophages was enhanced in 3d of cold stimulation ( Figure 7B ). Previous studies have shown the important role of M2 macrophages in thermogenesis. Notably, 6h-cold in BAT and 2d-cold in SAT both promoted the proportion of M2 macrophages (30, 31). Here, our study revealed the dynamic changes of M2 macrophages in different cold durations. Bioinformatics analysis demonstrated that the proportion of endothelial cells, a part of the immune microenvironment, differed from that of other stromal cells or immune cells. Thus, we detected such cells by flow cytometry and RT-qPCR and found that they sustained high levels during the longer period (3-7d) of cold stimulation in BAT ( Figures 7C, D and Supplementary Figure S5B ). However, these changes in SAT occurred in the initial 1d cold challenge ( Figures 7E, F ). Besides, ILC2 cells are known for their involvement in thermogenic regulation in SAT under 3 days of cold exposure (32). Their relative markers in both BAT and SAT were further studied. Our results showed that in SAT, the changes in markers IL5, IL13, and IL33 were consistent with those in ILC2 cells under 3-d cold stimulation, but not with those observed in BAT, indicating potential distinct regulations between these two types of adipose tissues ( Supplementary Figures S5C, D ). To summarize, all novel immune cells discovered in BAT and SAT in our studies are shown in Figure 8 and Supplementary Figure S5E .