Temporal trajectories of human brown and white adipocyte progenitors at single cell resolution

Nagendra Palani The Novo Nordisk Foundation Center for Basic Metabolic Research Pascal Timshel The Novo Nordisk Foundation Center for Basic Metabolic Research Pytrik Folkertsma The Novo Nordisk Foundation Center for Basic Metabolic Research Alexander Grønning The Novo Nordisk Foundation Center for Basic Metabolic Research Tora Henriksen Centre for Physical Activity Research Lone Peijs The Novo Nordisk Foundation Center for Basic Metabolic Research Verena Jensen Steno Diabetes Center Copenhagen Karolina Hvid Centre for Physical Activity Research Wenfei Sun Laboratory of Translational Nutrition Biology, ETH Zurich Naja Jespersen Rigshospitalet https://orcid.org/0000-0002-3488-1657 Christian Wolfrum ETH Zürich https://orcid.org/0000-0002-3862-6805 Tune Pers University of Copenhagen https://orcid.org/0000-0003-0207-4831 Søren Nielsen University Hospital of Copenhagen https://orcid.org/0000-0002-3119-750X Camilla Scheele (  cs@sund.ku.dk ) The Novo Nordisk Foundation Center for Basic Metabolic Research https://orcid.org/0000-0002-60559709


Introduction
Body fat distribution and quality are determinants of metabolic health. Multiple studies have revealed that abdominal obesity is strongly associated with cardiovascular disease and insulin resistance, whereas accumulation of fat in the lower gynoid regions has a lower lipid turnover and is associated with metabolic health 1 . In contrast to the common societal view that all excessive fat tissue is unhealthy, a functional adipose tissue allows for safe storage of excess calories. The essential role of safe lipid storage in a healthy adipose tissue is underscored by conditions like lipodystrophy, which is characterized by atypical accumulation of lipids and development of insulin resistance 2 . Adipocytes are lipid-storing and energy-transforming cells. Subtypes of adipocytes are linked to anatomical depots, displaying nichede ned functions. For example, the subcutaneous depot hosts white adipocytes with tremendous capacity for plasticity, in terms of energy access dependent lipogenesis and lipolysis plasticity 3 . In contrast, the supraclavicular deep neck depot contains brown adipocytes, which are heat producing with specialized uncoupling mitochondria that are turned on in response to cold induced sympathetic signals 4 . Adipocytes originate from mesenchymal stem cells that reside in multiple tissues including adipose tissue, skeletal muscle, and bone marrow 5 . Interestingly, studies suggest that brown adipocytes share some developmental similarities with the mitochondria-rich skeletal muscle cells 6 . Other studies have shown a close relation between osteocytes and adipocytes, which were found to share a transcriptional network 7 . It was found that this common network must be suppressed for cells to differentiate into adipocytes instead of osteocytes 7 . Intriguingly, it has been reported that progenitor cells derived from white adipose tissue of the brown bear, which exhibit an enormous seasonal plasticity, spontaneously differentiated into osteocytes in vitro 8 , further emphasizing the developmental link between these two cell types.

Results
Human WAT-and BAT-derived progenitors share cell fates during early differentiation To investigate adipose depot differences, we cultured 14 different adipose progenitor cell strains derived from supraclavicular or perirenal BAT or from visceral or subcutaneous WAT of adult humans 19,20 ( Figure 1a). We generated droplet-scRNA-seq data from these adipose progenitor cell strains. We used the R package Seurat for pre-processing the data, quality control, regression of cell cycle effects, sample alignment and differential expression analyses. We obtained 56,371 high-quality cells after quality control. For this initial analysis, cells were harvested at a proliferating, sub-con uent state ("T1") ( Figure   1a). For each depot, we included cells from at least three different individuals, matched for age, sex and body mass index (BMI) (Table S1). However, cell cycle effects explained most of the variation in the data ( Figure S1). We therefore regressed out cell cycle effects and kept both datasets for further analysis. We performed principal component analysis (PCA) and the rst 12 principal components were used for t-SNE visualization of the data without cell cycle effects. The t-SNE visualization of the data without cell cycle effects regressed out shows that cells tend to group to their own sample and to the phase of the cell cycle they are in ( Figure S1). To overcome possible batch effects, we aligned our samples with each other using canonical correlation analysis and dynamic time warping. A t-SNE visualization of the aligned data shows all samples aligned in one big cluster without any visible structure ( Figure S2). The alignment almost completely captures cell cycle effects ( Figure S2). A t-SNE visualization of the aligned data without cell cycle effects shows that the structure of the alignment does not change much when regressing out cell cycle effects ( Figure S2). These analyses suggest that the unsupervised clustering of multidepot-derived adipocyte progenitors was based on sample-speci c signals and cell cycle effects. In conclusion, our data suggest that undifferentiated, proliferating adipocyte progenitors have similar expression patterns regardless of BAT or WAT origin (Figure 1b). The high degree of similarity between human brown and white adipocytes at this developmental stage was surprising given our previous ndings that these cell-types are functionally different when fully differentiated 19,20 . However, despite lack of depot-dependent clustering, we decided to explore whether depot-dependent differences in single gene expression occurred. We thus performed differential expression (DE) analysis between the depots. We rst performed DE tests for each depot against the contrary two depots. We found 181 genes for supraclavicular, 129 for perirenal, 116 for subcutaneous and 108 for visceral (Bonferroni corrected P < 0.05, absolute average log fold change > 0.25). However, after inspecting these genes, we observed that most genes were found because of sample-speci c signals. Therefore, we also performed DE tests for each sample against the contrary two depots and only kept the intersecting genes for all samples belonging to the same depot. This resulted in a smaller list of differentially expressed genes per depot: 6 genes for supraclavicular, 19 for perirenal, 11 for subcutaneous and 26 for visceral (Bonferroni corrected P < 0.05, absolute average log fold change > 0.25) (Table S12). Because we used a low threshold for average log fold change, not all genes were equally important. For the brown adipocyte progenitor samples, 3 genes stood out. The rst gene, TM4SF1, was higher expressed in both the perirenal and supraclavicular depots (average log FC 0.87 and 0.34 respectively) and was found for all perirenal samples and 2 of the supraclavicular samples ( Figure S3). LY6K was also higher expressed in the perirenal and supraclavicular samples (average log FC of 0.33 and 0.38 respectively) and was found for 2 perirenal samples and 2 supraclavicular samples ( Figure S3). Last, HOXB7 was found as a negative marker gene for all supraclavicular samples (average log FC -0.42, Figure S3). Furthermore, we found two genes that seemed to discriminate the visceral samples well. BARX1 was found to be higher expressed in a small part of the cells in the visceral samples (average log FC 0.25, Figure S3), and LINC01116 was found as a negative gene for all visceral samples (average log FC -0.32, Figure S3). We next explored the data using Monocle and Velocyto, which underscored the similarities between the adipogenic progenitor cells, assigning cell cycle effects as the major driver of differences ( Figure S4). We observed one small 'common cluster', distinct from the main cluster and with cells from all depots (Figure 1b). This cluster was negative for CD29, a surface marker previously proposed to be a predictive marker for adipogenic precursor cells with thermogenic potential 21 . This observation was intriguing as the main cluster consisted of adipocyte progenitors from all four BAT and WAT depots. We further characterized the two clusters by using scmap 22 to map our data set onto an existing catalogue of cell-types in murine tissues 23 .
In this analysis, the main, CD29-positive cluster was enriched for adipogenic precursor markers, whereas the CD29-negative cluster was enriched for epithelial cell markers ( Figure S5). Most genes distinguishing between brown and white adipocytes, including UCP1, are not expressed until later in differentiation 24 . Thus, to discover transcriptional branching avenues between adipocyte progenitors derived from supraclavicular, perirenal, subcutaneous and visceral adipose depots, we harvested cells from all four depots during four additional time points T2-T5 during differentiation (Figure 1a, Table S1) and performed single-cell transcriptomics, obtaining 23,428 high-quality cells ( Figure S6). At T3, cells were two days post-con uence and a differentiation initiation cocktail was added (please see methods for components). Cells were subsequently harvested at T4 and T5, representing three and six days following addition of the differentiation cocktail, respectively. It should be noted that the full differentiation protocol into mature adipocytes is 12 days from the addition of the differentiation media 25 . However, changes in cell morphology are initiated shortly after adding the cell differentiation media and at T5 (corresponding to 6 days after adding the differentiation media), some accumulation of small lipid droplets has started.
The analysis resulted in clearly separated cell clusters both by time point and by depot within T4 and T5, whereas a PCA plot of the data shows that cell identities from T1-T3 are overlapping (Figure 1c, Figure   S6). To assess common or distinct developmental trajectories across depots, we used Monocle 26 to order cells in 'pseudotime' -a quantitative measure of progress through a biological process (Figure 1d). Notably, we applied a data-driven approach that did not include prior information on adipose depot origin of the cells. The trajectory topology was robust to changes in the input and parameter settings of the Monocle trajectory algorithm ( Figure S7). The trajectory analysis revealed bifurcating cell fates of adipose progenitor cells from the four different depots of human BAT and WAT. Cells from T1-T3 formed a progenitor (P) branch. In line with our initial analysis of proliferating precursor cells, cells from the earlier time points (T1-T3) did not separate in pseudotime and cells from all depots contributed equally to the (P) branch ( Figure 1d). Later in pseudotime, cells separated into two branches: upper branch (U branch) and lower branch (L branch) containing cells after induction of differentiation (T4 and T5) ( Figure 1d). Following branching, we observed a depot-dependent asymmetry in cell distribution where the U branch was dominated by cells from the brown fat depots (>60%) while cells from the white fat depots (>60%) were overrepresented in the L branch. Importantly, all four depots contributed to both upper and lower cell branches, suggesting a striking similarity also at this later developmental stage, ( Figure   1d). When t-SNE plots of cell branch identity was overlaid with time point and depot label, the branch separation occurring at both T4 and T5 was clear ( Figure 1e). Interestingly, a subpopulation of cells at T5, dominated by cells derived from the supraclavicular depot but including cells from all depots, was assigned to the P branch (Table S3). These cells were only sporadically expressed at T4 and clustered in the late part of the P branch in pseudotime, suggesting that dedifferentiation had occurred ( Figure S7).
The mechanism of dedifferentiation of adipocytes has been previously reported 27 and might re ect the ability of adipocytes to interconvert 28 . Supporting that dedifferentiation had occurred, DE genes between the T5 P branch cells and the T1-T3 P branch cells revealed residual high expression of both U branch and L branch speci c markers, suggesting unbiased contribution of dedifferentiated cells and gradual loss of branch speci c markers. These ndings also emphasize that neither of the cells in the U branch or L branch were undifferentiated or dedifferentiated progenitor cells, but rather represented two separate cell fates present among differentiating adipocytes. Importantly, we were able to con rm the differentiation trajectories predicted by Monocle pseudotime using an independent method called Velocyto 29 ( Figure S8).
The identi ed cell fates are also present in human BAT in vivo We next assessed if the subpopulations we observed in vitro also exist in human adipose tissue. We compared our data set with a published snRNA-seq dataset derived from human BAT; hereafter called the "BAT in vivo dataset". As BAT contains multiple cell types, we rst aimed to de ne a subpopulation within the dataset containing developing adipocytes, which could be compared with our single-cell dataset of human adipocytes. Therefore, we utilized the expression of two transcription factors important in brown adipogenesis, PPARG and EBF2, which were expressed in developing brown adipocytes in all three branches (progenitor, lower and upper), while highest expression was observed in the metabolic subpopulation ( Figure 2a). We labeled the EBF2 and PPARG expressing cells in cell atlases of the BAT in vivo dataset (Figure 2b, top) and labelled the sum of these cells (Figure 2b). We next performed a Seurat data integration using the in vitro adipocyte dataset containing P, U and L branches as reference to predict cell labels for the BAT in vivo dataset. Convincingly, we identi ed populations of all three branches (P, U, L) in the human BAT dataset among the cells that were predicted to be developing adipocytes based on PPARG and EBF2 gene expression ( Figure 2c). Interestingly, in vivo, EBF2 expressing cells were mainly (but not exclusively) clustering with the L branch cells, whereas in vitro we observed higher EBF2 expression in the U branch cells, raising the possibility of a switch between cell fates during adipogenesis. In support of the recapitulation of the P, U and L cell branches in the human BAT dataset, we validated this analysis using the Harmony algorithm 30 , which con rmed our conclusion ( Figure 2d).

Adipocyte progenitors split into structural and metabolic cell fates
We next explored the characteristics of the two identi ed cell-types by using BEAM 31 , a bioinformatic approach to identify branch-dependent genes (Table S4). We identi ed six gene clusters of branchdependent genes with distinct kinetic expression pro les (Figure 3a). The U branch was for example characterized by ADIPOQ encoding Adiponectin, a major regulator of lipid metabolism and insulin sensitivity 32 , and mitochondrial uncoupling protein 2 (UCP2) (Figure 3b). UCP2 is structurally in family with the brown fat speci c, thermogenic gene UCP1. The L branch was characterized by genes encoding locally secreted proteins. This included DCN, encoding Decorin, an extracellular matrix protein, and APOD, encoding Apolipoprotein D, which is also present in the extracellular region and is a component of highdensity lipoprotein (Figure 3b). Decorin has previously been reported to counteract adipogenesis 33 . To validate and visually examine the separation of the two cell fates within single cell cultures, we next performed FISH analyses, using RNAscope probes. As predicted, we observed a clear separation between our main markers, ADIPOQ and DCN, as well as between ADIPOQ and APOD expressing cells, whereas ADIPOQ and UCP2 expressing cells were overlapping as expected (Figure 3c). Following the con rmation that UCP2 was in the U branch, we assessed its co-localization with mitochondria in human brown adipocytes harvested at T5. We observed that cells highly abundant for the mitochondrial marker mtHSP70 were also positive for UCP2 mRNA (Figure 3d). Whereas UCP2 has not been demonstrated to rescue thermogenic function in UCP1 knockout mice 34 , the protein is located in the mitochondrial inner membrane, allowing for Ca 2+ to enter the mitochondria. The role of UCP2 in human BAT is to date unexplored and several examples of species differences in thermogenic regulation exist 35,36 . We here show that it is the highest expressed UCP halfway through full maturation of human adipocytes ( Figure  3d). To investigate if any of the cell fates demonstrated higher expression of thermogenic gene signature, we next grouped the cells by each branch, divided them into pseudotime deciles, and utilized the BATLAS 37 tool to predict brown and white adipocyte content ( Figure 3d, Table S6). BATLAS is a computational prediction model for identifying brown versus white fat signatures in samples of unspeci ed content. Strikingly, we found that the U branch cells from all four depots, i.e. both BAT and WAT, were predicted to gradually develop into brown adipocytes in a pseudotime-dependent manner, while L the branch cells developed in the opposite direction ( Figure 3d). This observation raised the idea that adipogenic cells exhibit a multipotent state halfway through full maturation, and that further thermogenic differentiation proceeds dependent on the adipose microenvironment.
This interpretation was consistent with a gene ontology analysis including all BEAM clusters that de ned the U and L branches. Multiple genes encoding proteins localizing to the mitochondria and being part of lipid metabolism and the respiratory chain accumulated in the U branch ( Figure 3e, Table S5). On the other hand, the L branch was clearly associated with an accumulation of processes related to extracellular matrix formation and developmental processes, including cellular adhesion ( Figure 3e, Table  S5). To re ect their predicted functions, we named the U branch cells "Metabolic cells" and the L branch cells "Structural cells". .

Transcriptional network de ning the metabolic and the structural cell types
To assess what de ned the metabolic and the structural cell fates during differentiation, we next explored the landscape of transcription factors expression in pseudotime, identifying TFs with increasing branchspeci c expression over pseudotime as drivers of branch development ( Figure 4a, Table S7). The metabolic cells were de ned by a well-established adipogenic transcription factor program 38 . This included Peroxisome proliferator-activated receptor gamma (PPAR-γ) and CCAAT/enhancer-binding protein alpha (C/EBPα), transcriptional activators that have been shown to synergistically activate adipogenesis by controlling genes important for adipocyte metabolism, including CIDEC, PLIN1, CD36 and ACSL1 38 all of which were also de ning the metabolic cells with higher expression compared to the structural cells. Additional transcription factors selective for the metabolic branch were SREBF1 and NR1H3, stimulators of lipogenesis 39,40 . Intriguingly, the structural cells induced a transcriptional program driving osteogenic proliferation and differentiation. This program included SNAI2 and JUNB, both promoting osteoblast maturation and forming a transcription factor network that counteracts adipogenesis 7 . Hence, the structural cells represent a cellular subpopulation that differentiates in a separate direction despite the addition of an adipogenic cocktail. These cells might be more involved in creating an adipose "skeleton", supporting the continuous adaptation of adipocyte size and structure by adjusting the extracellular matrix. To validate and visualize the separation between the selective transcription factors for the metabolic and structural cells at T5 of human brown adipocytes we performed a FISH analysis using RNAscope probes. We found a clear separation between several transcription factors including the metabolic transcription factors SREBF1 and NR1H3 versus the structural PRRX1, and the metabolic transcription factors PPARG and CEBPA versus the structural SNAI2 ( Figure 4a). We co-stained for additional branching transcription factors, demonstrating a similar cell population separation, with the exception of JUNB and PPARG that seemed to co-localize ( Figure S9). By taking advantage of our pseudotime-ordered dataset and applying a novel tool called Scellnetor 41 we were able to identify transcriptional networks involving a subset of the metabolic branch selective transcription factors (Figure 4b). Scellnetor is designed to identify gene networks enriched in differentiation trajectories 41 . This network analysis identi ed additional genes in the metabolic transcription factor networks including PPARGC1A, PPARGCIB, CIDEA and CKMT1B, further emphasizing the metabolic nature of this branch 42 . We next assessed whether there were any differences in the transcription factor dynamics between brown and white adipocytes. We grouped cells derived from supraclavicular and perirenal brown adipose depot into the "BAT" group and the cells derived from subcutaneous and visceral white adipose into the "WAT" group, and compared the pseudotime windows at which transcription factors diverged in expression between the metabolic and structural branches (Figure 4c). For example, we found that EBF2, a transcription factor well described for promoting brown adipogenesis 43 de ned the metabolic branch in brown adipocytes only, and rather early in pseudotime. However, most differences between brown and white adipocytes became clear later in pseudotime, for example RXRG, which also de ned the metabolic branch in brown adipocytes only. RXRG acts as a cotranscription factor with PPARG to promote thermogenic gene transcription 42 (Figure 4c). Among the genes de ning the structural branch were the above-mentioned SNAI2 and JUNB, driving this branch in both brown and white adipocytes. Interestingly, in late pseudotime, we observed NFIA, a structural branchselective gene in brown adipocytes only, previously reported to promote the brown fat differentiation program 44 (Figure 4c). By prede ning the cell origin from either brown or white adipose depots before sorting into pseudotime, we were thus able to identify branch-speci c transcription factors within brown or white adipocytes (Figure 4d). Some of these brown or white selective transcription factors are previously described in bulk data 24,42,44 , whereas several are novel and might prove powerful as directors of the brown and white adipocyte differentiation programs, respectively.
The Metabolic cell branch associates with traits for waist/hip ratio To investigate the relevance of our cultured cells in human health and disease, we used the computational tool CELLECT 45 , to genetically prioritize our pseudo-temporal ordered data for relevant human traits. Brie y, CELLECT quanti es the association between the expression patterns of cell-types and the genetic components of human complex traits identi ed by genome-wide association studies (GWAS). CELLECT tests for enrichment of genetic association signal for single nucleotide polymorphisms (SNPs) proximal to genes speci cally expressed in each cell-type. For this analysis, we delineated cell populations by stratifying cells by pseudotime decile for each developmental branch, generating 30 strata of cells. These cells from each strata were then used as input to CELLECT to estimate expression speci city of genes in any of these groups. We performed CELLECT analysis for 39 GWAS traits (Table S8). We identi ed signi cant enrichment of genetic association signal for cell populations in the U branch for fat distribution (as assessed by waist-to-hip ratio (WHR) and WHR adjusted for BMI (WHRadjBMI)), and lipid levels (as assessed by low-density lipid levels) (Figure 5b, Table  S9). Interestingly, the P and L branches displayed little or no enrichment. These results indicate that the branch point of adipogenesis, and the metabolic gene signature identi ed in this study is likely to be relevant for fat distribution and lipid metabolism. As a control of the model, we assessed the genetic enrichment for height and BMI GWAS and found no enrichment, consistent with the expectation that adipogenic development does not impact human height. Neither did we nd enrichment for any of the remaining GWAS traits analyzed, underscoring the strong impact of inducing the metabolic cell fate for determining fat distribution in humans. To understand the function and crosstalk between these two cell types, we predicted their cell-type speci c secretomes ordered by pseudotime. We utilized the published human secretome 46 as a scaffold to predict a secreted product from the branch speci c genes. We performed separate analyses for cells derived from BAT (supraclavicular and perirenal) and WAT (subcutaneous and visceral). In the metabolic cells, expression of several humanin-like peptides encoded by MTRNR2L1, MTRNR2L8, MTRNR2L10 and MTRNR2L12 caught our attention. These peptides are transcribed from the nucleus as homologues to the mitochondrial microprotein, humanin. They are transcribed from different chromosomes and include humanin-like 1, 8, 10 and 12 (Figure 5c). Humanin and the humanin-like peptides are 24 amino acid long peptides. Humanin associates with cardiometabolic function 47 and is regulated by IGF-1. Interestingly, IGF-1 is expressed from the structural branch raising the idea of a direct crosstalk between the branches. The structural cells encoded multiple extracellular matrix factors annotated as secreted. Among many others, these included several types of collagens and bronectin, proteins that are highly abundant in the development of obesity-induced brosis 48 (Figure 5d). The secreting capacity of the structural cells were underscored by the outnumbering of the secreting factors de ning the metabolic cells, with around double the amount of genes encoding secreted products in the structural cell compared to the metabolic cells (Figure 5e).

Discussion
We here investigate the rst phase of adipogenesis at single cell resolution across cells derived from four different human BAT and WAT depots. We use single cell trajectory analysis, creating a continuous gene signature using pseudotime analysis. We demonstrate that two cell fates arise from the same progenitor cells. This split occurs during the onset of differentiation and divides the cells into a metabolic and a structural cell type. Halfway through full differentiation, the metabolic cell type expresses a gene signature re ecting multiple metabolic pathways and a brown fat phenotype, regardless of BAT or WAT origin. The metabolic cellular gene signature is associated with genetic variants associated with body for fat distribution. The structural cells express developmental and extracellular matrix genes, many of which are upregulated during obesity-induced adipose tissue brosis.
Our ndings suggest a competing balance between two cell fates with complementing roles in the adipose tissue architecture and function. Importantly, previous studies have proposed that adipocyte progenitors can take on either a brogenic or an adipogenic path 49 . In the current study, we nd that these cell fates correspond to a structural and a metabolic cell type present in differentiating adipocytes from four different white and brown adipose depots of adult humans. These cell fates are thus closely related and interestingly their branching transcription factor signatures aligns with a previously described transcription factor network shared between adipocytes and osteoblasts 7 . This network controls a switch between differentiation into either cell type 7 . The authors demonstrated that suppression of the osteoblast transcription factor network allowed the cells to differentiate into adipocytes, and this was possible upon the induction of PPARG 7 . In the light of the data from Rauch and co-authors, our current results strongly suggest that the structural and the metabolic cell type arise from the same progenitor cells and that the cell fate is determined by asymmetric regulation of an adipocyte/osteoblast transcription factor network. These ndings raise the hypothesis of a closer interaction between these two cell types within the adipose tissue than previously anticipated, arguing for that this interaction should be taken into account when targeting obesity and metabolic disease. Interestingly, a speci c cell type that actively inhibits adipogenesis have been described in previous studies at single cell resolution 14 .
Others have demonstrated a subtype of adipocyte progenitors spatially associated with brotic structures and macrophages 12 . Importantly, the extracellular matrix provides a niche for the adipocytes and has been reported to affect the adipocytes to such an extent that cells stripped and swapped from the extracellular matrix of visceral versus subcutaneous adipocytes resulted in a phenotype corresponding to the origin of the extracellular matrix 33 . Moreover, a comparative secretomics study mapped a distinct extracellular matrix proteome between human brown and white adipocytes 50 , further supporting a role of extracellular matrix in de ning adipocyte characteristics. These observations support the idea of two contradicting cell fates, a structural-brogenic and an metabolic-adipogenic fate, with a balance between these two cell fates that can be shifted by external factors in the adipose microenvironment.
The metabolic cell type displayed a brown adipocyte signature, an intriguing nding, as this cell type was also common among the cells originating from WAT and BAT. Importantly, this nding can provide an explanation of the literature during the last decade describing browning of WAT, a process on differentiation of thermogenic adipocytes within WAT 51 . The phenomenon of thermogenic differentiation in WAT was rst described through induction of PPARG during in vitro differentiation of murine epidydimal adipocytes by adding pharmacological doses of rosiglitazone 52 . At the same time, PPARG is a well-described adipogenic transcription factor also essential for differentiation of white adipocytes 53 .
Taken together with our current ndings, these data suggest that the thermogenic signature is present at a premature state in all adipocytes, but restricted in the later maturation of white adipocytes, perhaps through WAT-speci c interactions between the co-differentiating metabolic and structural cells. Our data also brings perspective to previous ndings that beige adipocytes within WAT initially arise from environmentally primed progenitors, giving rise to de novo beige adipogenesis of adipocytes that subsequently can interconvert between a dormant and an active state 54,55 . Taken together, these observations underscore the thermal plasticity of adipose tissue and the capacity to adapt to environmental temperatures.
As we here show that the metabolic cells are present in all adipose depots and associate with common genetic variants associated with waist-to-hip ratio adjusted for BMI, it is relevant to discuss how metabolic/structural cell fate proportion could affect hip/waist ratio and metabolic health. If the metabolic cell gene signature fails to induce differentiation of a proper metabolic cell phenotype due to the associating SNPs, this would affect the metabolic cells across all BAT and WAT adipose depots in the body. Given the brown adipocyte signature of the metabolic cells, a reduced proportion would likely result in less thermogenic capacity in the BAT depots. In the subcutaneous WAT depots, it would likely result in reduced capacity for safe lipid storage, which in turn could result in an increased ectopic fat deposition in the viscera. Taken together, we nd that the metabolic gene signature associates with genetic variations for fat distribution, and this might in turn provide a common explanation to associations between unfavorable metabolic health and hip/waist ratio on one hand 1 , and BAT activity and cardiometabolic health on the other hand 56 . In support of this notion, we previously observed that UCP1 expression in human BAT directly correlates with waist/hip ratio 57 .
In conclusion, we here provide a seamless differentiation map of progenitor cells from four human BAT and WAT depots, revealing two main cell fates across all depots. Our data challenge the current conception of the origins of adipogenic and brogenic adipose resident cells and provides a missing piece of the puzzle for understanding browning of WAT. From a larger perspective, our data set the base for targeting obesity and its related diseases at a single cell speci c level. Declarations Foundation (www.cbmr.ku.dk) (Grant number NNF18CC0034900). The Centre for Physical Activity Research (CFAS) is supported by a grant from TrygFonden. During the study period the Centre for

Competing interests
The authors declare no competing interests.

Human samples
Human adipogenic progenitor cells were isolated from the stromal vascular fraction of the biopsies on the day they were obtained (surgery or biopsy), from the following four adipose regions: 1) visceral adipose tissue (obtained during gallbladder surgery); 2) perirenal adipose tissue (obtained during nephrectomy surgery); 3) abdominal subcutaneous adipose tissue (obtained with the Bergström needle biopsy method) and 4) supraclavicular adipose tissue (obtained during surgery in patients with suspected cancer in the neck area). Isolated cells were expanded and frozen in liquid nitrogen in a proliferative state as previously described 25

Cell culturing
Biopsies were collected in DMEM/F12 (Gibco) with 1% penicillin/streptomycin (Life Technologies) and tubes were kept on ice during transport. A detailed protocol for isolation and culturing of human adipocyte progenitors has been previously contributed 25 . Brie y, Biopsies were digested with 10 mg collagenase II (C6885-1G, Sigma) and 100 mg BSA (A8806-5G, Sigma) in 10 ml DMEM/F12 for 20 min at 37°C while gently shaken. Following digestion, the suspension was ltered, and cells were washed with DMEM/F12, resuspended in DMEM/F12, 1% PS, 10% fetal bovine serum (FBS) (Life Technologies) and seeded in a 25 cm 2 culture ask. Media was changed the day following isolation and then every second day until cells were 80% con uent; at this point, cultures were split into a 10 cm dish (passage 0). Cells were expanded by splitting 1:3. For the single-cell experiment, cells were seeded in 6 well plates in proliferation media consisting of DMEM/F12, 10% FBS, 1% PS and 1 nM Fibroblast growth factor-acidic (FGF-1) (ImmunoTools). Cells were grown at 37°C in an atmosphere of 5% CO 2 and the media was changed every second day. Adipocyte differentiation was induced two days after adipocyte progenitor cultures were 100% con uent by removal of FGF-1 from the media and addition of a differentiation cocktail (see details below) consisting of DMEM/F12 containing 1% PS, 0.1 μM dexamethasone (Sigma-Aldrich), 100 nM insulin (Actrapid, Novo Nordisk or Humulin, Eli Lilly), 200 nM rosiglitazone (Sigma-Aldrich), 540 μM isobutylmethylxanthine (Sigma-Aldrich), 2 nM T3 (Sigma-Aldrich) and 10 μg/ml transferrin (Sigma-Aldrich). After three days of differentiation, isobutylmethylxanthine was removed from the cell culture media and cells were differentiated for additional three days with the remaining differentiation compounds. For the 10x single-cell sorting, cells were loosened by adding 2 ml of TrypLex and placed in incubator for 3 min. Detachment of cells was con rmed by microscopy and 3 ml of Proliferation Media (PM) was added to the cells to inactivate trypsin. 190μl of cell-suspension was then transferred to an eppendorf-tube and mixed with 10µl of Solution 13 AO-DAPI, and then counted on a nucleocounter NC-3000. Cells were counted as described above, and 8000 cells/donor were pooled in an Eppendorf tube. The pool of cells was centrifuged for 7 min at 700g and resuspended in 80µl PBS with FFA-free BSA.

Single-cell library preparation and sequencing
Single-cell cDNA libraries were generated using Chromium Single Cell platform and 3' v2 Reagent Kit according to manufacturer's protocol (10x Genomics, USA). Single-cell libraries were sequenced on a NextSeq 500 platform to obtain 100 and 32-bp paired end reads using the following read length: read 1,

Proliferating adipocyte analysis
We used the R package Seurat 61 for pre-processing the data, quality control, regression of cell cycle effects, sample alignment and differential expression analyses. We performed quality control on our data to lter out low quality cells and genes, and we preprocessed the data to the format required for further analysis. As all proliferating adipocyte progenitors could be expected to have similar mitochondrial content, we ltered out cells where the mitochondrial gene expression was higher than 8%, as deviating high mitochondrial gene expression indicates stressed cells. Cells with less than 200 genes or more than 9,000 genes were also removed, as well as cells with more than 120,000 UMI's to remove possible doublets. The ltered data was log normalized and scaled, and the number of UMI's and percentage of mitochondrial genes were subsequently regressed out of the data. PCA was performed on the data and the rst 15 PC's were used for clustering and t-SNE visualization. Each cell was then scored for cell cycle phase. Please nd all details on the computational analyses at https://github.com/scheelelab/10xadipocyte-analysis Seurat We used the R package Seurat 61 (version 2.3.4) for preprocessing the data, quality control and differential expression analyses. As input to Seurat, we used the digital gene expression matrix output from the 10x Genomics analysis pipeline CellRanger. Cells with less than 200 genes and genes expressed in less than 3 cells were ltered out. The percentage of mitochondrial gene expression was calculated for each cell.
However, as mitochondrial gene expression increase during adipogenesis, we did not exclude cells with high mitochondrial expression in the developing adipocyte progenitor data set. We performed principal component analysis (PCA) to compute principal components (PCs) needed for clustering and data visualization. PC's were computed on the highly variable genes identi ed on log normalized and scaled data. Clusters were identi ed using a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm, which takes as input the signi cant number of PC's. The data was visualized using Seurat's implementation of t-distributed stochastic neighbor embedding 62 (t-SNE).

Monocle
We used the R tool Monocle 63 (version 2.8.0) to construct the cell developmental trajectory of the preprocessed Seurat object. Feature selection for trajectory construction was performed as follows: rst the dataset was split into two subsets, one containing all cells from T1, T2 and T3 and one containing all cells from T4 and T5. Both subsets were then clustered using Seurat's default clustering algorithm with a resolution of 1.5. Differential expression tests were performed for every cluster against the rest of the cells in the subset (using the negative binomial test, ltering on absolute logFC > 0.25). The union of the resulting gene list (2,464 genes) was used as input feature list for building the Monocle trajectory (DDRTree algorithm, max_components = 2). Monocle orders cells by pseudotime along the trajectories.
(Pseudotime is an abstract unit of progress: it is the distance between a cell and the start of the trajectory, measured along the shortest path. The trajectories total length is de ned in terms of the total amount of transcriptional change that a cell undergoes as it moves from the starting state to the end state.) We used Monocle's Branched Expression Analysis Modeling 64 (BEAM) to identify branch dependent genes. The genes in the resulting list were ltered on q-value < 0.05 (8,647 genes remaining) and subsequently ltered on absolute average logFC > 0.3 between the U branch and L branch (413 genes remaining). To subset transcription factors in the BEAM analysis, we retrieved the gene type, GO term name and GO term de nition for every gene in our dataset using Ensembl Biomart (version 96). From this set of annotated genes, we created two gene sets: a transcription factor gene set by selecting genes annotated with the 'transcription factor' GO term, and a non-coding gene set by ltering out genes annotated with the gene type 'protein_coding'.

Velocyto
We ran the Velocyto 65 command line tool for every sample with a genome annotation le and expressed repeat annotation le (reference genome hg19). We used the Python library for further downstream analysis. First the loom les of every sample were aggregated into one. Cells that were not present in the nal Seurat analysis were discarded and the metadata from the Seurat analysis was added to the remaining cells. We further discarded cells with extremely low unspliced detection, keeping 23,309 cells for the nal analysis. Genes were ltered by ranking the spliced genes based on a coe cient of variation vs. mean t, using the top 3000 to perform a PCA. Both the spliced and the unspliced gene expression matrices were subsequently normalized by size. Using the rst 15 principal components, the data was kNN-smoothed (using the default value of k; 0.025 x nCells). The standard implementation of Velocyto was used with default parameters for tting gene models, predicting velocity, extrapolating and plotting.

Gene set enrichment analyses
Gene set enrichment analyses were performed using the R package gPro leR 66 . GPro leR takes as input a list with gene symbols and returns a table with terms associated with those genes. We ltered the output to only contain Gene Ontology (GO) terms. To generate the gure with GO terms we performed gene set enrichment analyses on every cluster of branch-dependent genes identi ed using BEAM (6 clusters, 413 genes in total, absolute log fold-change > 0.3 between U branch and L branch). The resulting GO terms and their p-values were used as input for REViGO 67 , a web tool to summarize lists of GO terms by removing redundant terms. To visualize the summarized GO terms associated with branch-dependent genes, we used the R package GOplot 68 . GOplot calculates a z-score for each GO term indicating if the term is more likely to be decreased (negative value) or increased (positive value): , where up and down are the number of up-and down-regulated genes, respectively, counting genes with log fold-change > 0 between the U and L branch as up-regulated genes.

BATLAS
We used the webtool BATLAS 37 to predict the percentage of brown adipocyte content in our data. We grouped cells from each Monocle developmental branch (progenitor, lower and upper) by pseudotime decile, to generate 30 groups of cells. The average expression for every gene in each decile was subsequently calculated on the normalized data in non-log space. The resulting matrix was used as input for BATLAS.
Branch-speci c gene expression analysis in pseudotime.
Lists of human transcription factors 69 and of the human secretome 46 were used as input for identifying branch-speci c gene expression over pseudotime. For each gene, smoothed gene expression across cells along pseudotime was determined, and the pseudotime point at which expression increases in one branch compared to the other was identi ed. To be considered as a valid result, we required that genes that were identi ed as diverging in expression between branches maintained the trend of divergence until end of pseudotime.

Genetic prioritization analysis
We used CELLECT [CELL-type Expression-speci c integration for Complex Traits] 70 to genetically prioritize pseudo-temporal ordered groups of cells. Speci cally, we grouped cells from each Monocle developmental branch (progenitor, lower and upper) by pseudotime decile, to generate 30 groups of cells.
We used CELLEX [CELL-type EXpression-speci city] 70 to calculate expression speci city of these cell groups. Brie y, CELLEX normalizes the expression data using a common transcript count (assuming 10,000 transcripts per cell) and apply log-transformation. Next expression speci city likelihood (ES μ ) is computed for groups of cells. We used CELLECT with S-LDSC 71 as the genetic prioritization model. We ran CELLECT with default parameters (100 kb window size around each gene, correcting for baseline v1.1 and 'all genes' annotations, see Timshel2019 70 for details). We performed CELLECT analysis for 39 GWAS traits listed in Supplementary Table le, sheet 8.

RNA Fluorescent in Situ Hybridization
In vitro differentiated human adipocytes derived from either subcutaneous or supraclavicular deep neck adipose depots were xed at day 4 and 6 of differentiation with 10% neutral buffered formalin (Thermo Fisher Scienti c) for 30 minutes, dehydrated and stored in 100% ethanol until the staining procedure. In situ hybridization was performed using RNAscope Fluorescent Multiplex Plex v2 Reagent Kit and RNAscope manual assay probes designed and produced by Advanced Cell Diagnostics (ACD). Nuclei were stained with NucBlue ReadyProbes (Thermo Fisher Scienti c). RNA targets were visualized using an EVOS imaging system (Thermo Fisher Scienti c). The RNA targets were hybridized with RNAscope probes (table S11) and then labelled with the Opal 570, Opal 620 and Opal 690 (Akoya Biosciences). The uorescent signals were detected with an RFP, Cy5 and Texas red light cubes.

Validation of L and U branch cell-fates in human single-nucleus tissue
We obtained human deep neck BAT single-nucleus RNA-seq data 72 . The dataset comprises 50,284 cells from a total of 16 individuals. Cells were pro led using 10x Chromium V3. We applied standard Seurat work ow (as described previously) to QC and cluster the cells. We used Seurat v3 73 to integrate and transfer cell labels from our developmental adipocyte progenitors (T1-T5) to the Wolfrum group dataset.
Brie y, this approach consists of identifying 'cell anchors' between datasets that can then be used to harmonize the datasets and transfer information from one dataset to another. We used the functions FindTransferAnchors(reference= adipocyte progenitors_T1-T5, query=wolfrum_group_BAT) and TransferData() with default parameters to identify anchors and transfer P, U and L cell labels to Wolfrum group data. We retained only high con dent transferred labels with a transfer score above 0.7 (the transfer score ranges between 0-1, where increasing values indicate the con dence of the transfer) and labeled all other cells as 'non-matching cells'.

Data repository and code availability
All processed RNA sequencing data is deposited in ArrayExpress: E-MTAB-XXX and E-MTAB-XXX.

Supporting Information
Supplementary Tables S1-S12 are not available with this version  Comparing the single cell data set with human in vivo data. a) Feature plots highlighting the gene expression of adipogenic markers EBF2 and PPARG in blue, demonstrating the temporal expression of these markers in our in vitro dataset. Violin plots in the bottom demonstrate branch-speci c expression.
b) Cell atlases of brown adipose tissue single-nuclei RNA sequencing from 16 humans (the "BAT in vivo dataset") 72 Selective expression of adipogenic markers EBF2 and PPARG is highlighted in blue in the feature plots on top. The sum of this expression was used to estimate an adipocyte population, highlighted in red in the cell atlas in the bottom. c) The in vitro adipocyte data set containing P, U and L branches was used as reference for Seurat data integration and label transfer to predict cell labels for the BAT in vivo dataset. Gray-colored cells did not match any of the cells in the developing adipocytes data.
On top: Feature plots of main markers for the U (APIPOQ) and L (DCN) branches, respectively. d) Validation of Seurat data integration using Harmony. The gure shows cells projected in the integrated Harmony data embedding 30 . We used default parameters and recommended settings for the Harmony analysis.