Transcriptome analysis of browning human neck area adipocytes reveals strong influence by the FTO intronic SNP variant in addition to tissue and PPARγ specific regulation

Background Brown adipocytes, abundant in deep-neck (DN) area in humans, are thermogenic with anti-obesity potential. Although analyses of tissue samples and cell clones provided important findings, their molecular features are not fully understood. FTO pro-obesity rs1421085 T-to-C SNP shifts differentiation program towards white adipocytes in subcutaneous fat. Methods Human adipose-derived stromal cells were obtained from subcutaneous neck (SC) and DN fat of 9 donors, of which 3-3 carried risk-free (T/T), heterozygous or obesity-risk (C/C) FTO genotypes. They were differentiated to white and brown (long-term PPARγ stimulation) adipocytes, then global RNA sequencing was performed and differentially expressed genes (DEGs) were compared. Findings DN and SC progenitors had similar adipocyte differentiation potential but differed in DEGs. DN adipocytes displayed higher browning features according to ProFAT or BATLAS scores and characteristic DEG patterns revealing associated pathways which were highly expressed (thermogenesis, interferon, cytokine, retinoic acid, with UCP1 and BMP4 as prominent network stabilizers) or downregulated (particularly extracellular matrix remodeling) compared to SC ones. Part of DEGs in either DN or SC browning was PPARγ-dependent determining enrichment for lipid metabolism and induction of some browning genes such as PCK1, CPT1B, and PLIN5. Presence of the FTO obesity-risk allele suppressed the expression of mitochondrial and thermogenesis genes with a striking resemblance between affected pathways and those appearing in ProFAT and BATLAS, underlining the importance of metabolic and mitochondrial pathways in thermogenesis. Interpretation Among overlapping regulatory influences which determine browning and thermogenic potential of neck adipocytes FTO genetic background has a so far not recognized prominence. Funding European Regional Development Fund [GINOP-2.3.2-15-2016-00006] Research in Context section Evidence before this study Browning of white adipose tissue shifts adipocytes from energy storage white to energy expenditure beige types. The balance between the two adipocyte populations in white adipose tissue is highly determined by noncoding variants of the Fat mass and obesity-associated (FTO) locus which has the strongest association with obesity. The rs1421085 FTO risk allele results in a loss of ARID5B repression of IRX3 and IRX5 which promotes excess white adipocyte formation. Recent studies have revealed the presence of brown adipose tissues at several anatomical sites in humans including the deep-neck (DN). Tissue biopsies and brown adipocytes differentiated from immortalized clonal cell of DN origin have displayed high thermogenic potential which has been partially characterized in molecular and regulatory terms. Based on differential gene expression profiles of white and beige or brown adipocytes, ProFAT and BATLAS gene lists as tools were developed to estimate browning and thermogenic potential of tissue samples and cell populations. Added value of this study Ex vivo differentiated, DN derived stem cell populations, compared to those derived from subcutaneous neck fat, keep their higher browning potential displayed by phenotypic, UCP1 content and ProFAT as well as BATLAS scores. It has been revealed that characteristic gene expression profile and associated pathways of brown adipocytes are determined by partially overlapping effects of tissue site specific commitments of the stem cells, PPARγ stimulation and the FTO status of donors. The presence of FTO rs1421085 risk alleles have a strong influence, manifested during differentiation, on browning resulting in compromised expression of metabolic and mitochondrial genes as well as pathways which are decisive in thermogenesis. Furthermore, several molecular pathways could be newly linked to browning, such as highly expressed interferon, cytokine and retinoic acid mediated processes and downregulation of extracellular matrix components and modelling. Implications of all the available evidence Regulatory elements and pathways responsible for determining brown adipocyte and adipose tissue differentiation at specific anatomical sites are more diverse than originally thought. More knowledge of the FTO allele driven regulatory pathways in brown adipose tissue would provide possible pharmaceutical targets in obese patients, particularly those carrying the risk alleles.


Introduction
Brown and beige adipocytes play a major role in maintaining the constant core body temperature of hibernating, small and newborn animals, as well as in humans without shivering 1,2 . Their heat production is mainly mediated by UCP1, a mitochondrial carrier protein, which uncouples ATP synthesis from the respiratory chain activity 1,3 . These adipocytes also conduct effective UCP1independent thermogenic mechanisms, such as the ATP consuming futile cycle of creatine metabolism [4][5][6][7] , calcium redistribution between the endoplasmic reticulum and the cytosol 8 or through release of an enzyme, encoded by PM20D1 gene, which produces acylated amino acids that force uncoupling 9 . The stimulation of these processes leads to increased energy expenditure that can ameliorate the energy balance during obesity and type 2 diabetes mellitus 10,11 . Therefore, we need to learn in detail how regulatory networks drive the differentiation and thermogenesis of brown and beige adipocytes, especially in humans.
In rodents, brown adipose tissue (BAT) contains classical brown adipocytes which derive from myogenic precursors, accumulate numerous small lipid droplets in a multilocular arrangement and convert glucose and fatty acids into heat mostly by the action of the constitutively expressed UCP1 1,2,12 . Beige adipocytes with similar morphologic features were also described in mice.
These cells arise from mesenchymal precursors and have a common developmental origin with white adipocytes, exist in distinct thermogenic fat depots and can be induced by cold and subsequent adrenergic stimulation 2,13 . In adult humans specific adipose depots are enriched in brown adipocytes; these expanse to 1-1.5% of total body mass and are mostly found in the perirenal, deep-neck (DN) and paravertebral regions 14 . It is still unrevealed whether these thermogenic fat depots represent the classical brown or the beige type of adipocytes by origin, even after recent intense studies of DN tissue which could be compared to paired subcutaneous (SC) fat samples [15][16][17][18] . First, it was reported that supraclavicular human adipocytes possess a classical brown character 16 , then data derived from tissue biopsies or clonal cell lines lead to the conclusion that they rather resemble to mouse beige adipocytes 17,19 . In parallel, multiple studies using different approaches proposed specific marker genes for these browned adipocytes that are either connected them to the classical brown or to the beige lineage [16][17][18][19] . For simplicity, in this paper the term brown is used to cover both classical and beige thermogenic adipocytes in humans. 6 Human brown adipocytes differentiated in distinct adipose tissues express thermogenic genes at moderate levels under unstimulated conditions 2,14 . The ratio of brown and white adipocytes is partially determined during the early differentiation of mesenchymal progenitors into adipocyte subtypes which is strongly influenced by genetic predisposition 20,21 . This can be quantified in a given tissue or cell culture sample by determining BATLAS score based on the expression of 98 brown and 21 white-specific genes, which were selected by a combined analysis of gene expression signatures in mouse brown, beige and white adipocytes and human tissue samples 22 .
Similarly, browning probability can be estimated by a recently developed computational tool, ProFAT 23 .
A recent genome-wide association study of body fat distribution identified 98 independent adiposity loci which could affect the appearance of thermogenic fat 24 . In a detailed study by Claussnitzer et al. it has been described that single nucleotide polymorphism (SNP) rs1421085 underlies the genetic association between fat mass and obesity-associated (FTO) locus; obesity and the presence of the C risk-allele of the FTO locus results in a cell autonomous, IRX3 and IRX5 dependent shift in the gene expression programs toward white instead of brown adipocyte with lipid-storage and decreased in thermogenesis. When the T/T risk-free genotype is carried at the rs1421085 position, the ARIDB5 repressor is able to bind to this site not allowing IRX3 and 5 expression which allows adipocyte precursors to be committed for brown differentiation 20 . IRX5 knockout mice had reduced fat mass and were protected from diet-induced fat accumulation. In addition, IRX5 knockdown increased the mitochondrial respiration in isolated murine adipocytes 25 . It was also shown that IRX3 can promote browning of white adipocytes as it directly binds to the UCP1 promoter and its rare variants are associated with human obesity risk 26 .
To learn more about human adipocyte browning and attempting to through more light on unresolved or controversial issues in the regulation of thermogenesis, we decided to study neck adipocyte populations derived from primary human adipose-derived stromal cells (hASCs) instead of whole tissue samples with different cell types present or single cell clones influenced by immortalization protocols. We screened and compared global gene expression patterns by RNA sequencing of hASC-derived white and brown (in response to sustained PPARγ stimulation) differentiated adipocytes. The hASCs were isolated from paired DN and SC adipose tissue samples of nine donors, 3 of each FTO rs1421085 genotype: T/T-risk-free, T/C-7 heterozygous, C/C-obesity-risk. We found that both brown and beige markers, including UCP1, CKMT1A/B 4,15 , CIDEA 27 , or PM20D1 6 were upregulated in DN as compared to SC adipocytes, indicating a large potential of thermogenesis in the deep depots. Novel pathways and biological processes could be linked to browning regulation comparing patterns of upregulated genes. On the other hand, dozens of genes (such as CIDEA, CITED1 28 or PM20D1 thermogenic markers) were upregulated in response to the brown differentiation as compared to white irrespective of the anatomical origin of the hASCs. The gene expression pattern of brown adipocytes was determined to a greater extent by the anatomical origin of the hASCs from which they had been differentiated than the differentiation protocol. Surprisingly, the expression of metabolic, mitochondrial and thermogenic genes was strikingly compromised by the presence of FTO rs1421085 genotypes in the progenitors. Our results suggest that cultivated hASCs from distinct locations of the human neck still keep their differing propensity for adipocyte browning which is strongly influenced by the presence of the obesity-risk alleles at the FTO intronic locus. 8

Ethics statement and obtained samples, isolation and differentiation of hASCs
hASCs were isolated SC and DN adipose tissue of volunteers (BMI<29.9) aged 35-75 years who underwent a planned surgical treatment. During thyroid surgeries a pair of DN and SC adipose tissue samples was obtained to rule out inter-individual variations. Patients with known diabetes, malignant tumor or with abnormal thyroid hormone levels at the time of surgery were excluded 29 .
Tissue collection was complied with the guidelines of the Helsinki Declaration and was approved by the Medical Research Council of Hungary (20571-2/2017/EKU) followed by the EU Member States' Directive 2004/23/EC on presumed consent practice for tissue collection. hASCs were isolated and cultivated; white and brown adipocytes were differentiated from hASCs according to described protocols 6,7,29,30 . Briefly, both white and brown differentiations were induced by hormonal cocktails that contain apo-transferrin (Sigma-Aldrich cat# T2252), insulin (Sigma-Aldrich cat# I9278), T3 (Sigma-Aldrich cat# T6397), dexamethasone (Sigma-Aldrich cat# D1756), and IBMX (Sigma-Aldrich cat# I5879). Later, dexamethasone and IBMX were omitted from both media. In the white cocktail hydrocortisone (Sigma-Aldrich cat# H0888) was constantly present, while the brown contained insulin at 40x higher concentration than the white. However, the major difference between the two protocols was the time interval and concentration of the administered rosiglitazone (Cayman Chemicals cat# 71740). While the white regimen contained 2 µM rosiglitazone on the first four days of the two weeks long process, the differentiating brown adipocytes were treated with the drug at 500 nM concentration between the 4th and 14th days 31,32 . The absence of mycoplasma was checked by polymerase chain reaction (PCR) analysis (PCR Mycoplasma Test Kit I/C, PromoCell cat# PK-CA91) 7,29 .

Flow cytometry
To investigate the phenotype of the undifferentiated hASCs, a multiparametric analysis of surface antigen expression was performed by three-color flow cytometry using fluorochrome-conjugated antibodies with isotype matching controls. See references 7 and 33 for further details about the analysis.

RNA and DNA isolation, genotyping
RNA and DNA preparation was carried out as described previously 6,7,29,30 . Rs1421085 SNP was genotyped by qPCR using TaqMan SNP Genotyping assay (Thermo Scientific cat# 4351379) according to the Manufacturer's instructions 7 .

RNA-Sequencing
To obtain global transcriptome data, high throughput mRNA sequencing analysis was performed on Illumina sequencing platform. Total RNA sample quality was checked on Agilent BioAnalyzer using Eukaryotic Total RNA Nano Kit according to the Manufacturer's protocol.
Samples with RNA integrity number (RIN) value >7 were accepted for the library preparation process. Libraries were prepared from total RNA using NEBNext® Ultra™ II RNA Library Prep for Illumina (New England BioLabs, Ipswich, MA, USA) according to the manufacturer's protocol. Briefly, poly-A RNAs were captured by oligo-dT conjugated magnetic beads then the mRNAs were eluted and fragmented at 94 0 C degree. First strand cDNA was generated by random priming reverse transcription and after second strand synthesis step double stranded cDNA was generated. After repairing ends, A-tailing and adapter ligation steps adapter ligated fragments were amplified in enrichment PCR and finally sequencing libraries were generated. Sequencing run were executed on Illumina NextSeq500 instrument using single-end 75 cycles sequencing.
After sequencing the reads were aligned to the GRCh38 reference genome (with EnsEMBL 95 annotation) using STAR aligner (version 2.7.0a) 34 . To quantify our reads to genes, featureCounts was used (version: 1.6.3) 35 .
Subsequent gene expression analyses were performed in R (version 3.5.2). Genes with low expression values or with outlier value were removed from further analysis. Briefly, after removing the highest read value of each gene we filtered out genes with reads was less than 50 considering raw reads count. Then, to further remove outlier genes we calculated Cook's distance and filtered out genes with Cook's distance higher than 1. After filtering out the low-expressed and outlier ones, the expression profile of 22362 transcripts could be examined. PCA analysis did not show any batch effect, considering sequencing date, the donor origin, and the donor sex or tissue origin. Differential expression analysis was performed using DESeq2 algorithm (version 1.22.2). When Tissue origin and differentiation protocol based differential expression was determined donor origin was controlled to decrease the effect of biological variance on thermogenic capacity 36 . However, when comparison was based on the FTO obesity-risk allele presence, we did not control donor origin. Significantly differentially expressed genes (DEGs) were defined based on adjusted P values <0.05 and log2 fold change threshold >0.85.
Hierarchical cluster analyses and heat-map visualization was performed on the Morpheus web tool (https://software.broadinstitute.org/morpheus/) using Pearson correlation of rows and columns and complete linkage based on calculated z-score of DESeq normalized data after log2 transformation. The z-score was calculated in two ways: to eliminate the donor background we

Antibodies and immunoblotting
Lysis of differentiated adipocytes, denaturation, SDS-PAGE and blotting were performed as described previously. For overnight, membranes were probed at 4ºC with primary antibodies: temperature. Immunoreactive proteins were visualized by Immobilion western chemiluminescence substrate (Merck-Millipore cat# WBKLS0500) 7 .
Immunofluorescence staining, quantification of browning hASCs were plated and differentiated on Ibidi eight-well μ-slides; vital and immunofluorescence staining was carried out as described previously 6,7,29,30 Figure 1A). After the cells received a single bolus dose of cell permeable dibutyril-cAMP, mimicking adrenergic stimulation, we found that adipocytes which were differentiated from DN precursors were more capable to induce their respiration than the adipocytes of SC origin. A similar trend was detected when the differentiated adipocytes were treated with oligomycin that blocks the activity of ATP synthase and provides information on proton leak respiration. In parallel, basal and cAMP-stimulated extracellular acidification rates of DN adipocytes were greater in a significant degree than those of the SC ones (Supplementary Figure 1B). We could also estimate, by applying the creatine analog β-GPA, the contribution of UCP1-independent creatine futile cycle to oxygen consumption 4 and found it was more pronounced in DN adipocytes (Supplementary Figure 1C).
As a further characterization of the hASCs, surface antigen analysis was performed. There was no significant difference in the presence of hematopoietic/monocyte, fibroblast, endothelial and integrin cell surface markers between the SC and DN precursors ( Figure 1A and Supplementary   Table 2) were expressed significantly higher (except STAT3) in all mature adipocyte samples as compared to progenitors and were not significantly different between SC and DN adipocytes or samples differentiated with white and brown protocols. Hierarchical cluster analysis of these genes confirmed that the samples clustered into two main groups according to their differentiation status but not by tissue origin or the differentiation protocol used, suggesting that these factors had no effect on the general adipocyte differentiation efficiency of the progenitor cells ( Figure 2B).

Differentiated adipocytes from DN progenitors display higher browning-related geneexpression features
Next, we investigated morphology and differential transcriptomic characteristics of adipocytes differentiated from hASCs of the nine donors. Laser-scanning cytometric data clearly showed 15 that the DN adipocyte population, differentiated by either the white or the brown protocol, had more brown cells characterized by the appearance of UCP1 and small lipid droplets ( Figure 2C).
To compare tissue origin and differentiation protocol-dependent transcriptional changes related to browning, recently developed computational tools ProFAT and BATLAS were applied. When the ProFAT browning gene-set 23 was used to score browning probability, browning markers were  Figures 3A and B). Relating the DEGs in mature adipocytes to how they were expressed in preadipocytes, four distinct groups could be defined. Figure 3A shows a global heat map representation of the expression profile in the four groups.
257 genes (Group 1) were expressed at a greater extent in differentiated DN adipocytes as compared to SC ones and this differential expression was not observed in DN versus SC preadipocytes ( Figures 1B and C). to maintain network integrity. Other genes, such as ESR1, MT2A, LEPR2, IRF7, and AGT may be also pivotal for the network (Table 1A).
272 genes (Group 2) were expressed at a greater extent in DN as compared to SC preadipocytes and their expression remained elevated in DN adipocytes as compared to SC ones after white or brown differentiation, presumably contributing to higher propensity to browning ( Figures 3A   and B). The most affected pathways were complement and coagulation cascades, signaling by retinoid acid and interaction between L1 and ankyrins ( Figures 3E and F, Supplementary   Figure 4B by Gephi). As expected, pathways in this group were already revealed in the analysis of higher expressed genes in DN versus SC preadipocytes (Figures 1B-D (Table 1B). TBX1, a well-established beige marker gene 17 , also appeared in this group as an important network member.

Downregulated pathways in DN adipocytes
Next, we investigated the lower expressed genes in the DN samples. 270 genes (Group 3) were downregulated in differentiated DN adipocytes but not in SCs compared to preadipocyte differences and this differential expression was not seen in the preadipocyte state; expression of some of these genes may be non-permissive for browning. STRING network and pathway analysis showed that these genes were predominantly involved in remodeling of the ECM, cell adhesion and synthesis of GPI-anchored proteins pathways (Figures 4A and B, Supplementary   Figure 4C by Gephi). Cytokine signaling was also highlighted as several interleukin receptors (e.g. IL1RE, IL18RE, IL20RE), TGFB3 and KIT were expressed higher in the matured SC adipocytes (Figures 4C and D). COL1A1 gene played a central role in the network, but the GPC3, KRT7, CNR1, and SEMA3A were also important in sustaining a stable network structure (Table 2A).
250 genes (Group 4) were expressed at a lower extent in DN preadipocytes as compared to SC and their expression remained low in the DN adipocytes as compared to SC ones after white or 18 brown differentiation; these genes may be part of pathways which are not permissive for browning either. ECM organization pathways had outstanding dominance among the related pathways and the repeated appearance of the GPI-anchored protein pathway (SC preadipocytes and Group 3) emphasizes the importance of glycosylphosphatidylinositol anchored membrane proteins in the formation of the SC phenotype. In addition, the transcription regulation by homeobox and homeodomain and the neuroactive ligand-receptor interaction pathways also appeared in this group (Figures 4C and D, Supplementary Figure 4D by Gephi). When analyzing relationships, the POSTN gene was the most significant, but RUNX2, NCAM1, GATA2, EDIL3, and PAX3 could be also important for maintaining network connectivity. Based on their network position, the IRX1, IRX3, and GRIA3 genes expressed higher in SC adipocytes may have roles in restraining browning potential (Table 2B).

Shared PPARγ mediated gene expression patterns in DN and SC adipocytes
Both ProFAT and BATLAS estimations clearly showed that SC precursor cells also reacted to browning as compared to the white differentiation protocol by upregualting brown genes. The degree of upregulaton of these genes could be similar in DN and SC adipocytes and missed when DN and SC data were compared. Therefore, we analysed genes which responded similarly to brown differentiation protocol, that is long-term rosiglitazone exposure, in SC and DN adipocytes. Indeed, out of the 217 genes responding by upregulation to brown differentiation in either DN or SC adipocytes only 40 (Supplementary Table 3) were also present in Groups 1 and 2 (529 genes) described above. Eighty genes were expressed higher in brown adipocytes irrespective to their anatomical origin and they displayed enrichment of the PPAR signaling pathway (Figures 5A-C). Regarding the latter glycerol kinase, PCK1, CPT1B were present among BATLAS genes 22 , PLIN5 was linked to browning 42  reprogramming of a PPARγ super-enhancer 43 . We found that KLF11 was among the rosiglitazone upregulated genes ( Figure 5D) in either SC or DN adipocytes suggesting that it may be involved in the regulation of browning of neck area adipocytes.
Among the upregulated browning genes outside of the PPAR pathway, CIDEA was listed by ProFAT and also found in the Group 1 genes above, similarly to PM20D1 which was strongly influenced by natural genetic variations in humans 44 and CPT1B which was essential for betaoxydation of fatty acids. Gamma-butyrobetaine dioxygenase (BBOX1) catalyzes the formation of L-carnitine from gamma-butyrobetaine and has been reported to play a role in adipocyte browning 45 Figure   5D) also demonstrated that the majority of the mitochondrial genes were expressed at a lower extent in samples that carried the obesity-risk genotype. The genes CKMT1A and CKMT1B, which encode key mitochondrial creatine kinases and reported to be involved in UCP1independent thermogenesis 4 , were among those being significantly less expressed in FTO obesity-risk samples. All these data underline fundamental differences in the metabolic processes and mitochondrial activity between T/T and C/C adipocytes with significant deficiency in those obtained from the FTO obesity-risk donors.

STRING analysis of the DEGs influenced by the presence of FTO obesity-risk allele revealed
interaction network of gene products and significantly enriched pathways with defined biological function which, by being dysregulated or deficient, may significantly contribute to the manifestation of the obese phenotype. Among the genes that were expressed at lower levels in samples carrying the FTO obesity-risk allele, metabolic pathways were most affected including lipid metabolism, thermogenesis, carbon metabolism, oxidative phosphorylation, and degradation of certain amino acids such as valine, leucine and isoleucine (Figures 6D and E). There was a striking resemblance between the set of network pathways negatively affected by the C/C genotype and those which were demonstrable analyzing the ProFAT and BATLAS gene products ( Figures 6H-K). Of the 23 KEGG significantly enriched pathways defined by ProFAT genes, 21 are also found in BATLAS, and 19 of these are defined by genes that showed suppressed expression in FTO C/C samples. While the 29 KEGG pathways defined by BATLAS genes 25 were also found in pathways determined by the lower expressed genes in C/C samples. Analyzing the Reactome pathways, we found similarly high agreement (Supplementary Table 5A-B).

Interactome network analysis of the lower expressed genes in FTO obesity-risk genotype
suggests that under normal conditions, DECR1 (dienoyl-CoA reductase), the chemokine CCL2, LIPE, LDHB, SOD2, and ATP5B and J are central elements in maintaining connectivity (Table   3A). Genes that were linked to adipo -and thermogenesis such as SLC2A4, EHHADH, PPARGC1A, CPT1B, and FABP4 appeared with high betweenness scores. In addition to the metabolic genes, a cytokine ligand, CCL2, and a receptor, CXCR4, also appeared as important critical components of this interactome network.
Comparing the FTO status related gene expression patterns to upregulated DEGs in DN versus SC adipocytes, we found that out of the 529 genes (Groups 1 and 2) that were significantly higher in the DN samples, 33 genes were poorly expressed in the FTO obesity-risk genotype samples Since the presence of the FTO risk allele may lead to loss of restrain on the expression of genes related to thermogenesis and its regulation, we also looked for and found DEGs which showed higher expression in obesity-risk genotype samples (Figures 6F and G). Clearly, the expression of genes related to the organization of the ECM was the most prominent upregulated pathway from this respect. This is consistent with the finding that the ECM organization pathways has was also among the genes kept lowly expressed in DN and FTO T/T adipocytes but became higher expressed in the adipocytes with the obesity-risk genotype.

Discussion
Human ASCs from distinct fat depots have the potential to undergo a browning program 14 . One of the major sites that contain active BAT in humans is located in the neck, particularly in its deep regions 14,48 . When we compared SC and DN adipocyte progenitors from human neck, the latter had greater browning potential, in accordance with previous findings 18,19,39,49 . On the other hand, SC hASCs were also able to build up a significant thermogenic competency. This was supported by gene expression, morphological and functional properties of the differentiated adipocytes. Primarily, we found 1420 DEGs when the hASCs of SC and DN origins were compared. This relatively large number of DEGs and their designated pathways suggest that the two anatomically localized tissues have cells of different character at precursor level and presumably have a specific ECM structure that is maintained ex vivo and after differentiation.
The organization of the ECM is differently regulated: degradation processes are prominent in the SC samples, whereas in DN samples ECM reorganization and regulation of vascular development appear to be significant. The retinoid acid signaling pathway is also essential and determinant in DN cells both at progenitor and at differentiated states, whose potential role in inducing thermogenesis via increasing angiogenesis by activating VEGFA/VEGFR2 signaling have been suggested but not explored in details 50,51 . It has been also recognized that the development of the neuronal system is different at the two tissues 52 and the precursor cells may be partly responsible for this.
Many genes were already proposed to predispose cells to a higher or lower thermogenic potential at the precursor level 49,[53][54][55]  Majority of the gene expression studies aiming to determine the molecular, functional and developmental differences between heat-producing brown and fat-storing white adipocytes target mature adipocytes, primarily in mice but also in humans [15][16][17][18][19] . Interestingly, very few overlapping confirmatory results were obtained from these types of studies. It is increasingly recognized that human and mouse-based studies cannot be directly compared due to the differences in environmental exposure and the apparent difference in body mass to body surface area, which fundamentally influences body temperature homeostasis 56,57 . Therefore, humanized mouse conditions were recently optimized to generate more conclusive results 58 . It is also difficult to directly compare tissue studies with in vitro cellular assays and many factors could also contribute to differences observed between cellular models (e.g.: cell origin, passage number, immortality process, media composition, etc.) 59          Venn diagram shows significantly higher expressed genes number in samples expected higher thermogenic activity, and among of these genes several were lower expressed in FTO obesityrisk genotype samples (B) shows significantly lower expressed genes number in samples expected higher thermogenic activity, and among of these genes many expressed higher in FTO obesity-risk genotype samples. SC: subcutaneous; DN: deep-neck.