Comparative transcriptomic analysis of rabbit interscapular brown adipose tissue whitening under physiological conditions

ABSTRACT Interscapular brown adipose tissue (iBAT) of both rabbits and humans exhibits a similar whitening phenomenon under physiological conditions. However, a detailed characterization of iBAT whitening in them is still lacking. Here, we chose rabbits as a model to gain a better understanding of the molecular signature changes during the whitening process of iBAT by transcriptomic analysis of rabbit iBAT at day 1, day 14, 1 month and 4 months after birth. We applied non-invasive MRI imaging to monitor the whitening process and correlated these changes with analysis of morphological, histological and molecular features. Principal component analysis (PCA) of differentially expressed genes delineated three major phases for the whitening process as Brown, Transition and Whitened BAT phases. RNA-sequencing data revealed that whitening of iBAT was an orchestrated process where multiple types of cells and tissues participated in a variety of physiological processes including neovascularization, formation of new nervous networks and immune regulation. Several key metabolic and signalling pathways contributed to whitening of iBAT, and immune cells and immune regulation appeared to play an overarching role.


Introduction
Obesity is a global pandemic and is associated with many metabolic diseases, such as type 2 diabetes, atherosclerosis, cardiovascular diseases, hypertension, fatty liver and some cancers [1]. Notably, China may already have the most obese people in the world according to Chinese criteria [2]. In order to improve quality of life and relieve medical discomfort as well as economic burden, effective ways to prevent and combat obesity for the healthy and sustainable development of human beings are urgently needed.
Activation of brown adipose tissue (BAT) and/or recruitment of beige tissue have been considered as a promising therapeutic strategy to counteract obesity and type 2 diabetes in recent years [3][4][5]. There is general consensus that two major adipose tissues exist in mammals: BAT and white adipose tissue (WAT). WAT is a major lipid reservoir for excess energy in the form of triglycerides and also acts as an endocrine organ. Beige adipose tissue has been recognized recently as being derived from WAT but behaves metabolically like BAT. Functionally, BAT and beige adipose tissue differ from WAT due to their capacity for uncoupling mitochondrial ATP synthesis from the electron gradient potential to generate heat through the mitochondrial inner membrane protein, uncoupling protein 1 (UCP1) [6]. This process occurs in non-shivering thermogenesis (NST) and is believed to be important for the maintenance of normal body temperature in mammals, such as rodents, hibernating animals and newborn infants [7,8]. Unlike rodents where BAT is present functionally throughout their lifespan, human beings and other animals such as ruminants and rabbits [9][10][11] are believed to possess functional BAT for only a short period of time after birth whereafter it is replaced by WAT [12]. Thus, historically it has been widely accepted that adults possess little or no metabolically active BAT. However, metabolically active BAT in adult humans was identified in 2009 by positron emission tomography/computed tomography with the glucose analog, F18-fluorodeoxyglucose (FDG-PET /CT) [13][14][15]. These findings reignited interest in brown fat biology and propelled the idea that BAT may be a promising target to augment energy expenditure and reduce obesity.
Adult human BAT is mainly present in the cervical, supraclavicular, and paravertebral areas, and can be activated by cold exposure [13][14][15][16]. Decreased BAT content is associated with accumulation of body fat with age [17]. BAT can also be recruited through repeated stimulation by cold exposure or capsinoid intake, leading to decreased total body fat content in the absence of body weight change [16].
In human infancy, BAT exists principally in the interscapular area, namely iBAT [8,18]. iBAT content is persistent within the first decade but gradually dissipates where it is barely detectable between 30 and 80 years of age [19]. Instead, a fat pad with a phenotype similar to WAT emerges to take its place [20], suggesting that brown adipocytes within iBAT change into white adipocytes with age in humans. A similar whitening phenomenon occurs in rabbits where iBAT is developmentally reprogrammed to a 'WAT-like' phenotype after birth [12,21]. However, a detailed molecular description together with a non-invasive live imaging system to monitor the dynamic changes during the process of iBAT whitening is lacking.
Here, we utilized a rabbit animal model, magnetic resonance imaging (MRI) technology, histology and whole-genome RNA-Seq analysis to correlate MRI data with histological and molecular changes to gain insights into the process of iBAT whitening.

BAT converts to white-like adipose tissue in rabbits
To gain a better understanding of the whitening process of BAT, we systematically and non-invasively examined the developmental changes of adipose tissues in rabbits using MRI. Three rabbits for each group at 1 day, 14 days, 1 month, 2 months, 3 months and 4 months after birth were selected (Figure 1a).
The fat fraction of interscapular, dorsal and cervical adipose tissues increased with age (Figure 1b, 1c, S1A and S1B). Visual inspection after anatomical dissection revealed that the colour of iBATs and dorsal adipose tissues (dBATs) were brown in 14 days old rabbits, and iBATs were darker than dBATs (Figure 1d). However, by 4 months, adipose tissues at the same sites became white in appearance, indicating loss of iBAT content (Figure 1d). This visual change was confirmed by haematoxylin and eosin staining showing that the adipocytes in iBAT and dBAT reprogramed phenotypically from multilocular brown adipocytes to unilocular white-like adipocytes (Figure 1e and S1C). In contrast, adipocytes within the inguinal adipose tissue (ingWAT) depot appeared to be unilocular initially but became larger during the whitening process (Figure 1e).
Corresponding analyses of UCP1 protein expression in iBAT and dBAT showed increased expression from day 1 to 14, followed by a gradual decrease and disappearance at 2, 3 and 4 months of age (Figure 1f and S1D). By comparison, UCP1 protein expression in ingWAT was not detected at any time point by western analysis (Figure 1f). However, multilocular UCP1 positive cells in rabbit ingWAT are detected when immunohistochemical analysis was conducted, which is similar to that of the inguinal adipose tissue of young mice [22,23]. The ingWAT has multilocular UCP1 positive cells on day 1 rabbits and these cells essentially disappeared by day 14 (Figure 1g). Overall, these results show that iBAT and dBAT gradually convert with age to a white-like adipose tissue phenotype in rabbits

The transcriptional profile of iBAT shifts towards that of WAT during whitening
To delineate the regulatory mechanisms underlying the whitening process of rabbit BAT, we performed wholegenome transcriptomic analysis of rabbit iBAT at different ages. We selected iBAT at day 1, day 14, 1 month and 4 months for RNA sequencing analysis based on UCP1 protein expression content as determined by western blot analysis (Figure 1f). The ingWAT of 4-month-old rabbits was used as a control for iBAT after whitening. Fifteen samples in total were sequenced and 110.05 Gb of clean data were obtained with each sample greater than 6.05 Gb. Average clean reads for iBAT_D1, iBAT_D14, iBAT_M1, iBAT_M4 and ingWAT_M4, 47 (Table S1).
To compare the time-dependent transcriptional changes of rabbit iBAT globally, we performed PCA of RNA sequencing data obtained from rabbit  50 µm. (f) UCP1 protein levels of rabbit iBAT and ingWAT at different ages. (g) UCP1 staining of iBAT_D14, ingWAT_D1 and ingWAT_D14 of rabbits. Scale bars, 20 µm. Data represent mean ± SEM, **P < 0.01, ***P < 0.001, ****P < 0.0001. iBAT_D1, iBAT_D14, iBAT_M1, iBAT_M4 and ingWAT_M4 (Figure 2a). Gene expression analyses indicated that these samples could be divided into three different groups. The first group comprising iBAT_D1 and iBAT_D14 had remarkably similar transcriptomic signatures. The second group includes iBAT_M1 and the third group comprises iBAT_M4 and ingWAT_M4 (Figure 2a). These three groups clustered as distinct phases during the whitening process of iBAT. Thus, iBAT_D1 and iBAT_D14 represent the Brown Phase, while iBAT_M4 represents a Whitened BAT Phase similar to ingWAT_M4. Since iBAT_M1 has a transcriptional profile between Brown Phase and Whitened BAT Phase (Figure 2a), it was defined as a Transition Phase.
Although iBAT_M4 and ingWAT_M4 displayed closely related transcriptomes, they could nevertheless be subdivided into two different groups by PCA ( Figure 2a). This analysis indicates that Whitened iBAT and ingWAT are not identical and that fundamental differences may exist between them. Consistent with these results, global correlation analysis of adipose tissue transcriptomes during iBAT whitening revealed that iBAT_D1, iBAT_D14 and iBAT_M1 were tightly correlated ( Figure 2b and Table S2). iBAT_M4 was clearly correlated with ingWAT_M4, but its correlation was weaker than that between iBAT_M4 and iBAT_M1 ( Figure 2b and Table S2). In contrast, ingWAT_M4 displayed a poor correlation with iBAT_D1, iBAT_D14 and iBAT_M1 ( Figure 2b and Table S2). Taken together, these results indicate that rabbit iBAT undergoes a profound transcriptomic shift during whitening.

Identification of new marker genes for BAT and WAT
We performed a differential gene expression analysis on samples from each group and identified differentially expressed genes (DEGs) between samples from each group ( Figure S2A, Table S3 and Table S4). These data show that the number of DEGs in iBAT increased with age during the whitening process ( Figure S2A). In addition, the number of DEGs between iBAT and ingWAT decreased with age although iBAT_D1 vs ingWAT_M4 and iBAT_D14 vs ingWAT_M4 had a comparable number of DEGs ( Figure S2A). The changes in the number of DEGs also suggest that the expression profile of iBAT becomes closer to that of ingWAT during the process of whitening.
Numerous marker genes have been proposed to differentiate among brown and beige and white adipocytes or brown and white adipose tissues in mice and humans [24][25][26][27][28]. We first examined whether the expression patterns of respective marker genes in rabbits fit those in mice and humans. As shown in Figure 2c, conventional thermogenic marker genes such as PRDM16, UCP1, PPARGC1A, and DIO2 were highly expressed in the Brown phase. Their expression decreased in Transition and Whitened phases. By comparison, white adipocyte marker genes such as RETN, LEP, AGT, SLC16A12, HOXC8, HOXC9 and ADIPOQ exhibited an opposite pattern with increased expression during whitening. Thus, the expression patterns of the above genes in rabbits were consistent with those in mice and humans.
Nevertheless, some BAT and WAT marker genes exhibited unique rabbit expression patterns. For example, mitochondrial tumour suppressor 1 (MTUS1) is highly enriched in UCP1-positive human adipocytes and identified as a molecular marker in addition to UCP1 for assessing thermogenic adipocyte content within human adipose tissues [27]. However, in rabbits, MTUS1 displayed higher expression in the Whitened phase as compared to the Brown and Transition phases ( Figure 2c). Likewise, some other genes which are brown adipocyte markers in mice and humans, such as CIDEA, SLAC29A1 and LHX8, were enriched in the Whitened phase and ingWAT in rabbits ( Figure 2c).
Interestingly, TBX1 and TMEM26, which represent beige adipocyte marker genes in mice, exhibited different expression patterns in rabbits. Specifically, TBX1 was expressed at higher levels in ingWAT as compared to iBAT. By comparison, TMEM26 displayed a relatively low expression level in ingWAT and in iBAT at day 1, comparable expression levels in Transition and Whitened phases, and a higher expression level in iBAT at day 14 ( Figure 2c). The different expression pattern may reflect rabbit specific physiology. Overall, our findings indicate that the majority of marker genes found for BAT and WAT in mice and human also apply to rabbits.
We then determined whether additional BAT and WAT marker genes could be identified from our RNA data and therefore compared DEGs between iBAT_D1 and ingWAT_M4 in rabbits to those of humans BAT and WAT and mouse brown and white adipocytes as reported in BATLAS [29]. By comparing DEGs from all three species, a total of 188 BAT and 102 WAT marker genes were identified using │log 2 FC│≥ 0.5 and Padjust ≤ 0.05 as screening criteria (Figure 2d, 2e, S3A and S3B). This finding suggests the existence of a panel of evolutionarily conserved genes that are fundamental for brown and white adipocyte function. In addition, 100 BAT and 128 WAT marker genes were identified that are shared only by humans and rabbits (Figure 2d, 2e, S3C and S3D). By definition, BAT marker genes should represent properties characteristic of brown adipose tissue. Consequently, these genes should express at the highest levels in the Brown phase, at a lower level in Transition phase, and at the same or even lower level in the Whitened phase and ingWAT_M4. In contrast, WAT marker genes should express at the highest level in ingWAT_M4 and/or iBAT_M4, and at the lowest level in the Brown phase. As can be seen in Figure S3A-S3D, the expression patterns of some putative BAT and WAT marker genes do not align with the above criteria, and thus visual inspection and manual adjustment were put in place to eliminate these outliers. After manual elimination, only 80 BAT and 99 WAT marker genes were shared among humans, rabbits and mouse (Table 1). Also, 23 BAT ( Figure 2f) and 126 WAT marker genes were only shared between humans and rabbits ( Table 1).
Closer examination of the additional 126 WAT marker genes revealed that only 12 genes, i.e. RNF112, LMX1B, FMOD, TMEM119, GDF10, ADRA2A, HAND2, IRX1, HOXB6, HOXB8, UNC5C and FBLN7, exhibited an expression pattern characteristic of WAT. These genes were highly expressed in ingWAT_M4 with comparable expression across the three phases of iBAT whitening (figure 2f). Therefore, these genes could serve as bona fide new WAT marker genes that may play an instrumental role in maintaining WAT characteristics. The remaining 114 genes may play a role in promoting the whitening process since their expression levels gradually increased during iBAT whitening and eventually reached levels comparable to those in ingWAT. Similarly, the well-characterized brown fat marker genes, such as PRDM16 and PPARGC1A, as well as several potential new brown fat marker genes, such as ADRB3, CKMT1B and CKMT2, were present in the 23 BAT marker genes. Therefore, the remaining 18 genes could serve as new potential marker genes for BAT (figure 2f). The function of these newly identified marker genes for BAT or WAT, however, requires further investigation and characterization in terms of their roles in the browning and whitening of adipose tissue as well as adipose biology.
We next carried out a Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analysis of higher expression genes in the Brown, Transition and Whitened BAT phases. KEGG pathway enrichment analysis revealed that genes highly expressed in the Brown Phase were significantly enriched in metabolic pathways related to mitochondrial function ( Figure 3b and Table S5). Genes within these pathways were down-regulated during the whitening process ( Figure 3c). DEGs with higher expression in the Transition Phase showed that only the Pertussis pathway was significantly enriched when setting Padjust ≤ 0.05. With a more relaxed threshold P ≤ 0.05, 28 pathways including the nod-like receptor signalling pathway, peroxisome function and autophagy were significantly enriched ( Figure S4A and Table S5). Finally, 47 KEGG pathways were significantly enriched in the Whitened BAT Phase such as protein processing within the endoplasmic reticulum, endocytosis, lysosome, ubiquitin mediated proteolysis, N-glycan biosynthesis, AMPK and the sphingolipid signalling pathway ( Figure S4B and Table S5).
Mitochondrial function is closely tied to mitochondrial content. Adipose tissue mitochondrial homoeostasis is tightly regulated by a balance between mitochondrial biogenesis and degradation [30]. ESRRA and PGC1α (PPARGC1A), two major regulators of mitochondrial biogenesis, were down-regulated during the whitening process (Figure 3d), suggesting mitochondrial biogenesis is decreased. On the other hand, 74 autophagy-related genes were up-regulated ( Figure S4C), 34 of them were mitophagy-related ( Figure S4D), implying mitochondrial degradation is likely increased during the whitening process. Consistent with these results, mtDNA copy number decreased sharply from day 1 to 14 after birth, followed by a gradual decline and eventual plateau to the same level as ingWAT at 4 months of age (Figure 3e), indicating mitochondrial content is diminished during rabbit iBAT whitening.

Expression of genes involved in angiogenesis and innervation is augmented during iBAT whitening
The STEM program was used to classify the 14,899 DEGs into 26 major possible model profiles based on temporal gene expression patterns (Figure 4a). Seven significantly different (P < 0.05) gene expression patterns were identified and defined as Profile 15,16,25,8,24,22 and 13 (Figure 4a). A gene set was created for KEGG and Gene Ontology (GO) pathway enrichment analysis utilizing these seven Profiles. KEGG pathway enrichment analyses indicated that 96 pathways were significantly enriched, such as lysosome activation, axon guidance and regeneration, synapses of dopaminergic, glutamatergic and cholinergic nervous systems, platelet activation and the VEGF signalling pathway (Figure 4b and Table S6).
GO enrichment analysis also revealed that 1276 pathways, such as blood vessel development and morphogenesis, angiogenesis, axonogenesis, positive regulation of axon extension and axon guidance ( Figure 4c and Table S7), were significantly enriched. DEGs associated with blood vessel development and morphogenesis ( Figure S5A), platelet activation ( Figure S5B), the VEGF signalling pathway ( Figure S5C), angiogenesis ( Figure S5D) and axon guidance ( Figure S5E) were more abundant in the Transition and Whitened BAT Phases as compared to the Brown Phase. Indeed, the expression of factors that promote vasculogenesis was up-regulated ( Figure 4d). On the other hand, the expression of axon guidance cues SEMA3A, which induces RHOA translation [31], pan-neuronal marker genes, as well as glutamate neuronal marker genes, was also up-regulated during the whitening process ( Figure 4d). These findings indicate an essential role for angiogenesis and neo-vascularization during the whitening process.

Expression of immune-related genes and infiltration of immune cells are increased during iBAT whitening
Our findings showed that 95% of DEGs were upregulated during iBAT whitening. We believed that the triggering forces for the whitening process of rabbit iBAT must come from the biological processes driven by these upregulated genes. We therefore reasoned that if these up-regulated genes control and co-ordinate the whitening process, then their expression would be relatively stable during the Brown phase (iBAT_D1 and iBAT_D14) and increased during the iBAT_D14 to Transition phase (iBAT_M1). As shown in Figure 4a, the gene expression patterns of profile 15 and profile 16 met the above criteria.
GO analysis uncovered 326 and 474 significantly enriched pathways for Profile 15 and Profile 16, respectively (Table S8 and Table S9). For Profile 15, most of the pathways were closely associated with immune  processes such as regulation of cell proliferation for lymphocytes, mononuclear cells and leukocytes, regulation of innate and adaptive immune response ( Figure 5a and Table S8). For Profile 16, enriched pathways include those for histone modification, chromatin organization, cell cycle processes, protein ubiquitination, covalent chromatin modification and regulation of chromatin modification ( Figure S6A and Table S9). Many of these pathways are intimately linked to remodelling of chromatin and epigenetic regulation.
For further verification, GO analysis was performed using DEGs between iBAT_D1 and iBAT_D14, iBAT_D14 and iBAT_M1 as well as iBAT_M1 and iBAT_M4. Up-regulated DEGs between iBAT_D1 and iBAT_D14 were significantly enriched in metabolic pathways related to regulation of T cell, lymphocyte and leukocyte differentiation and activation, and regulation of brown fat cell differentiation ( Figure S6B and Table S10). As expected, these pathways were similar to those for Profile 15.
Between iBAT_D14 and iBAT_M1, up-regulated DEGs were also significantly enriched in metabolic pathways related to immune processes, including regulation of innate and adaptive immune responses, and regulation of cell activation comprising leukocyte, myeloid leukocyte, lymphocyte, macrophage and B cell activation (Figure 5b and Table S10). Monocyte marker gene CD14, dendritic cell marker gene CD40, NK cell marker gene CD59, macrophage marker gene CD68, M1 macrophage marker genes CD80 and CD86, M2 macrophage marker genes IL10, CD163 and ARG1, were all up-regulated from the Brown to Transition Phases (Figure 5c). These findings suggest an increased number of immune cells by either proliferation or infiltration during iBAT whitening.
The infiltration and passing through tissues of immune cells are coordinated by a large number of chemokines and their receptors [32]. Also, some chemokines and chemokine receptors such as CXCL8, CCL24, CXCL11, CCL14 and CXCR1, are only present in humans and not in mice. We therefore examined the expression of these genes during iBAT whitening. It is intriguing that the expression levels of the chemokines, CXCL8, CCL24, CXCL11, and the chemokine receptor, CXCR1, were much higher in the Transition phase as compared to the Brown and Whitened BAT phases (Figure 5c). These findings further suggest that the infiltration of immune cells likely increases during iBAT whitening.
The up-regulated DEGs between iBAT_M1 and iBAT_M4 were significantly enriched in metabolic pathways associated with lipogenesis, principally including acylglycerol, neutral lipid, triglyceride and sterol biosynthetic processes (Figure 5d and Table  S10). Similar processes were reported for the whitening of murine BAT in response to thermoneutrality, where lipogenesis was prevalent during iBAT whitening [33].
In contrast, down-regulated DEGs between iBAT_D1 and iBAT_D14 were significantly enriched in processes involving oxidoreductase activation as well as fatty acid and lipid metabolism ( Figure S6C and Table S10). The down-regulated DEGs between iBAT_D14 and iBAT_M1 were significantly enriched in metabolic pathways involved in mitochondrial function ( Figure S6D and Table S10). Unexpectedly, downregulated DEGs between iBAT_M1 and iBAT_M4 were significantly enriched not only in processes comprising the respiratory electron transport chain, regulation of cold-induced thermogenesis and adaptive thermogenesis, but also in immune response and immune-related processes (Figure 5e and Table S10).
Subsequently, we examined the expression of genes involved in immune regulation between iBAT_D14 and iBAT_M1 as well as between iBAT_M1 and iBAT_M4. Several immune-related genes were initially upregulated between iBAT_D14 and iBAT_M1 and then down-regulated between iBAT_M1 and iBAT_M4 (figfigure 5f and Table 2). Taken together, these results indicate that the up-regulation of immune-related genes precedes the up-regulation of lipid synthesis relevant genes and that immune cells are likely the primary determinant driving iBAT whitening.

Identification of potential transcription factors that promote rabbit iBAT whitening
ChIP-X Enrichment Analysis 3 (ChEA3) was used to identify potential transcription factors (TFs) that drive iBAT whitening (https://amp.pharm.mssm.edu/ChEA3, Keenan et al., 2019). We first looked for TFs that control up-regulated and down-regulated DEGs between iBAT_D1 and iBAT_M4, iBAT_D1 and iBAT_14, iBAT_D14 and iBAT_M1 as well as iBAT_M1 and iBAT_M4. The top 10 TFs between each group are listed in Table 1. During iBAT whitening, transcription factors CHCHD3 and ESRRA, which regulate mitochondrial function showed higher expression in the Brown phase as compared to the Transition and Whitened BAT phases (Figure 6a). This finding indicates that mitochondrial function diminishes in the later phases of iBAT whitening.

Discussion
In this study, we established a detection technique to accurately correlate non-invasive MRI measurement with morphological, histological and molecular events during physiological iBAT whitening in rabbits. To the best of our knowledge, MRI imaging coupled with transcriptomic analysis of iBAT whitening under physiological conditions has yet to be reported.
Adipose tissue is a highly plastic and dynamic tissue. In response to physiological stimulation, significant changes can be induced in its metabolism, structure, and phenotype [35]. WAT can acquire a BAT-like cellular and molecular program in response to various stimuli, such as cold exposure, β-adrenergic receptor agonists or genetic alterations, in a process termed browning [3,36,37]. BAT is present throughout the life of mice and provides resistance to diet-induced obesity through UCP1 [38,39]. Nevertheless, BAT whitening leading to obesity and insulin resistance can also be induced in mice by various interventions, such as high fat diet regimens, thermoneutrality, and genetic manipulation [28,[40][41][42][43]. Rabbits have a relatively wide thermoneutral zone of 15-25°C [44,45]; however, the cold environment does not prevent rabbit BAT whitening [46], suggesting that thermoneutrality is not the main reason driving iBAT whitening.
Decreased mitochondrial biogenesis and activation of mitochondrial autophagy are also important mechanisms driving BAT whitening in mice [30,47,48]. Under physiological conditions, BAT whitening in rabbits is spontaneous. Our findings suggest that this process involves a decrease in mitochondrial biogenesis and increased autophagy, which may lead to an overall decreased mitochondrial content and number. Recently, BAT-specific deletion of the TFEB gene, a master regulator of lysosomal biogenesis and autophagy, attenuated BAT whitening in mice at thermoneutrality [49]. These results support the idea that inhibiting autophagy and/or promoting mitochondrial biogenesis of brown adipocytes counteracts the process of BAT whiteness [30,48,50].
In mice, vascular rarefaction resulting from chronic high fat feeding regimens is a significant causal factor for BAT whitening. Chronic high fat diet causes BAT whitening, and Vegfa expression was significantly down-regulated at 12, 16 and 20 weeks on a high fat diet compared with that of mice on a normal chow diet [51]. Targeted deletion of Vegfa in adipose tissue resulted in BAT whitening even in mice fed with a normal chow diet [42]. Denervation of BAT is also reported to cause its whitening [28]. Nevertheless, our RNA seq analysis shows that VEGFA expression was in fact upregulated during iBAT whitening in rabbits. KEGG pathway and GO enrichment analysis also found that expression of genes involved in both angiogenesis and innervation was increased during rabbit iBAT whitening, indicating that angiogenesis and innervation are likely boosted, which suggests that the rabbit iBAT whitening process differs from that of mice.
Angiogenesis is also closely related to adipogenesis [52] where the vascular system transports nutrients, oxygen, growth factors, cytokines and hormones required for adipocyte function, growth and survival; the vascular system also controls changes in the microenvironment for adipose tissue [53]. Innervation of sympathetic nerve fibres is also essential for BAT thermogenesis; however, glycinergic stimulation or vagal afferent activation inhibits BAT sympathetic nerve activity and thermogenesis in rats [54,55]. BAT and WAT depots in mice have abundant vascular and nerve distribution as shown by 3D volume fluorescenceimaging [56,57]. Therefore, this imaging technique could be used to determine the distribution of blood vessels and nerves as well as the types of synapses involved during rabbit iBAT whitening in future studies. An adipose tissue is an extraordinarily heterogeneous organ with mature adipocytes constituting less than 50% of the adipose cell fraction, while the rest is classified as the stromal vascular fraction comprising many different cell types, including various immune cells, endothelial and vascular cells, fibroblasts, and adipose precursor cells [58][59][60]. Furthermore, mature adipose cells and various other cell types in adipose tissue together orchestrate adipose tissue development and function [61,62]. Our detailed transcriptomic analysis shows that up-regulation of genes involved in immune-related processes precedes the up-regulated expression of genes regulating lipid synthesis, and upregulated expression of marker genes in a variety of immune cells, such as macrophages, monocytes, dendritic cells and NK cells during rabbit iBAT whitening. These findings indicate that immune cells may participate in and orchestrate the whitening process. However, adipocytes themselves as well as adiposederived stromal cells could also express immunerelated genes. Increased IRF7 expression in adiposederived stromal cells of ageing mice leads to impaired mitochondrial function and BCAA metabolism [63]. Therefore, we cannot rule out that non-immune cells increase the expression of immune-related genes during rabbit iBAT whitening.
The immune system gradually matures during infancy [64] with maturation of immune competence after 1 year in humans, and 8 weeks in rabbits [65]. Therefore, the timing of iBAT whitening appears to be synchronized with the maturation of the immune system, further emphasizing the importance of immune cells in this process. It has been demonstrated that senescent T cells produce and release IFN-γ to inhibit differentiation of preadipocytes from brown adipocytes, thereby contributing to age-induced BAT whitening in mice [66]. Single-cell RNA sequencing (scRNA-seq) and single-nucleus RNA-sequencing (snRNA-seq) technologies have greatly advanced our understanding of the cellular complexity and plasticity of adipose tissue, and can be utilized to identify all major cell types in adipose tissue [67]. Detailed cellular atlases of human and mouse subcutaneous and visceral WAT at singlecell resolution have recently been generated, and subpopulations of various cells, such as adipocytes, adipose stem and progenitor, vascular and immune cells, as well as potential cell-cell interactions have been identified [68]. Using snRNA-seq, it was found that a population of CYP2E1+/ALDH1A1+ adipocytes secrete acetate by paracrine means to suppress adjacent adipocyte thermogenesis [69]. These results indicate that cellular crosstalk between adipocytes themselves, as well as with various other cell types such as immune cells in the adipose microenvironment influences adipose function. Future work may employ scRNA-seq to further pinpoint which specific immune cells play a central role in the whitening process of rabbit iBAT. As a matter of fact, while this manuscript was being prepared, a report was published in Cell Reports where scRNA-seq analysis was used to reveal that brown adipocyte progenitor cells expressing FSTL1 contribute to rabbit iBAT whitening, and their subsequent results elegantly demonstrated that removal of the Fstl1 or Fstl1 + progenitors results in BAT paucity in mice [46]. Despite the fact that the iBAT in their article refers to the dBAT depot in our study, we noticed an increased proportion of immune cells such as B cells, T cells and macrophages during the whitening process in rabbits [46], which is consistent with our conclusion that the immune cells and immune regulation may orchestrate the physiological whitening of iBAT in rabbits.
Transcriptional and epigenetic regulation are two of several mechanisms that regulate gene expression to ensure coordinated cellular behaviour and fate determination. In the current study, the discovery of TFs that regulate iBAT whitening reinforces an important role of immune-related genes. In addition, SREBF1 and SREBF2, which encode SREBP-1a, SREBP-1c, and SREBP-2, respectively, as well as MLXIPL (also known as ChREBP), are major TFs that drive adipocyte de novo lipogenesis [70]. Their expression was upregulated during iBAT whitening, suggesting enhanced de novo lipogenesis. Recent studies also show that BAT whitening is enhanced by BAT-specific ChREBP-β overexpression and prevented in ChREBP-deficient mice [33,43]. Whether identified TFs, such as AR, TWIST2, and the IRF family promote BAT whitening remain open research avenues. Moreover, epigenetic regulations such as DNA methylation and histone modification alter epigenetic effector activity, thereby altering the chromatin landscape and gene regulation, ultimately affecting adipocyte function [28,71,72]. Our results also show that epigenetic regulations, such as chromatin and histone modification, are enhanced during iBAT whitening.

Limitations
Current research focuses on bioinformatic analysis to discover potential regulatory mechanisms promoting iBAT whitening in rabbits, and our findings require further experimental validation. Since it is currently difficult to conduct genetic manipulations in rabbits, especially for construction of rabbit models with tissuespecific gene deletions, some verification studies may need to be carried out in mice.

Conclusion
In summary, we demonstrate that the phenomenon of iBAT whitening occurs in rabbits using MRI imaging together with anatomical and histological observations. Thus, brown adipocytes in BAT change from small multi-locular fat droplets to large white adipocytes-like single-locular fat droplets. The present study provides valuable insights into gene regulation during BAT whitening. A detailed transcriptome analysis was performed to understand the underlying regulatory mechanisms of iBAT whitening in rabbits under physiological conditions. We found that most of the marker genes used to distinguish BAT and WAT in humans and mice can be applied to a rabbit model. New marker genes were also identified and used to discern BAT and WAT. We also observed that iBAT whitening was associated with an up-regulation of approximately 95% of genes, suggesting it is an active process where extensive cellular remodelling or reprogramming takes place. Transcriptomic analyses highlighted that blood vessels, innervation and immune cells within iBAT change dramatically during whitening. Our findings are consistent with a mechanism whereby immune cells and immune regulation orchestrate the whitening process and increase lipid biosynthesis through transcription and epigenetic regulation in adipocytes. The genes and pathways associated with the iBAT whitening revealed here are candidates for subsequent validation studies and could serve as potentially new targets for the prevention and treatment of obesity by reactivating the Whitened BAT.

Animals
Male New England rabbits of different age states were purchased from Guangzhou Bai Yun District Long Gui Xing Ke Animal Farm. Rabbits were housed at a temperature of 22-23°C with free access to food and water. According to the literature [12] and the results of our MRI pilot experiment, rabbit samples were divided into \six groups by the age: 1 day (D1), 14 days (D14), 1 month (M1), 2 months (M2), 3 months (M3) and 4 months (M4). Each group had 3-6 rabbits. iBATs, dBATs and ingWATs were collected, then frozen in liquid nitrogen and stored at −80°C until analysis. All rabbit experiments were approved by the Animal Care and Use Committee of the Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences.

Haematoxylin-eosin staining and immunohistochemistry
Brown and white adipose tissues were fixed in 4% formaldehyde overnight at room temperature and embedded in paraffin before sectioning, then cut into 5 μm section with a microtome. Slides were deparaffinized, rehydrated and stained with haematoxylin and eosin and photographed under bright-field microscopy. Alternatively, sections were incubated with UCP1 antibody (MAB6158) followed by DAP chromogenic reaction. Imaging was performed on a slice scanner (Pannoramic MIDI).
RNA-seq transcriptome libraries were prepared using a TruSeqTM RNA sample preparation kit from Illumina (San Diego, CA) with 1 μg of total RNA. mRNA was isolated according to poly A selection method using oligo(dT) beads and then fragmented by fragmentation buffer. Double-stranded cDNA was then synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Synthesized cDNA was subjected to end-repair, phosphorylation and 'A' base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantification by TBS380, paired-end RNA-seq sequencing libraries were sequenced with the Illumina NovaSeq 6000 sequencer (2 × 150 bp read length).
Data were analysed on the free online platform of Majorbio Cloud Platform (www.majorbio.com). Methods of analysis are briefly described below.

Differential expression analysis
The expression level of each transcript was calculated according to the transcripts per million reads (TPM) method. RSEM [74] (http://deweylab.biostat.wisc.edu/ rsem/) was used to quantify gene abundance. Differential expression analysis was performed using DESeq2 [75] where Q value ≤ 0.05, and DEGs with Q value ≤ 0.05 and |log2FC| > 1 were considered to be significantly different expressed genes.

Identification of novel BAT and WAT marker genes
To identify potentially new marker genes for BAT and WAT, we compared the DEGs between iBAT_D1 and ingWAT_M4 in rabbits with those of human BAT and WAT and mouse brown and white adipocytes as reported in BATLAS [29]. A │log2FC│ ≥ 0.5 and Padjust ≤ 0.05 was used as discriminating criteria to filter possible BAT and WAT marker genes. The candidate genes were further examined by inspection and comparison to iBAT_D14, iBAT_M1 and iBAT_M4 before manual elimination of outliers was carried out.

Gene Ontology and Kyoto Encyclopaedia of Genes and Genomes pathway enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and genomes (KEGG) pathway enrichment analysis were performed to identify which DEGs were significantly enriched in GO and KEGG terms with a Bonferroni-corrected P-value (Padjust) ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tan ghaibao/Goatools) and KOBAS (http://kobas.cbi.pku. edu.cn/home.do) [76].

Time course differential expression analysis
Time course differential expression analysis was performed according to maSigPro (http://www.bioconduc tor.org/packages/release/bioc/html/maSigPro.html), which can identify genes with different changes in the whole time-series node samples [77]. The Hclust clustering algorithm was adopted to obtain six clustering numbers.

Temporal expression trend analysis
The Short Time-series Expression Miner (STEM) program [78] was used to identify temporal expression profiles. The STEM clustering algorithm was used and all parameters set to the default value. The results show that the time-series expression pattern with colour conforms to a significant change trend, while the timeseries expression pattern without colour is a statistically non-significant change trend. Profiles with the same colour belong to the same cluster.

Statistical analysis
Data are expressed as mean ± SEM. GraphPad Prism 8.0 was used for graphing and statistical analysis and one-way ANOVA with multiple comparison adjustment by Tukey's test was used. P < 0.05 was considered statistically significant.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This study was supported in part by an International Partnership Program of the Chinese Academy of Sciences.