Single cell full-length transcriptome of human subcutaneous adipose tissue reveals unique and heterogeneous cell populations.

White adipose tissue (WAT) is a complex mixture of adipocytes and non-adipogenic cells. Characterizing the cellular composition of WAT is critical for identifying where potential alterations occur that impact metabolism. Most single-cell (sc) RNA-Seq studies focused on the stromal vascular fraction (SVF) which does not contain adipocytes and have used technology that has a 3' or 5' bias. Using full-length sc/single-nuclei (sn) RNA-Seq technology, we interrogated the transcriptional composition of WAT using: snRNA-Seq of whole tissue, snRNA-Seq of isolated adipocytes, and scRNA-Seq of SVF. Whole WAT snRNA-Seq provided coverage of major cell types, identified three distinct adipocyte clusters, and was capable of tracking adipocyte differentiation with pseudotime. Compared to WAT, adipocyte snRNA-Seq was unable to match adipocyte heterogeneity. SVF scRNA-Seq provided greater resolution of non-adipogenic cells. These findings provide critical evidence for the utility of sc full-length transcriptomics in WAT and SVF in humans.

Introduction

Adipocytes in white adipose tissue (WAT) serve as robust lipid reservoirs and as powerful secretory cells—both of which are critical to maintaining whole-body glucose homeostasis and insulin sensitivity (Goodpaster and Sparks, 2017). Clone-based analyses of adipocytes have postulated that transcriptionally and metabolically distinct adipocytes exist beyond the standard classification of white, brown, and beige/brite (Lee et al., 2017, 2019; Min et al., 2019); however, little is known about adipocyte heterogeneity in human WAT in vivo. While adipocytes account for >90% of WAT volume, they only represent <50% of WAT cellular content (Corvera, 2021); thus, overall metabolic homeostasis within WAT also depends on the large portion of non-adipocyte cells. Characterizing the cellular composition of WAT is critical for identifying where potential alterations occur in these cellular processes that impact WAT metabolism.

Recent advancements in single-cell (sc) RNA-Seq technology have enabled high-throughput transcriptional profiling of thousands of non-adipocyte cells derived from WAT (Acosta et al., 2017; Burl et al., 2018; Cho et al., 2019; Gu et al., 2019; Hepler et al., 2018; Hildreth et al., 2021; Merrick et al., 2019; Schwalie et al., 2018; Vijay et al., 2020). The majority of the studies have been performed in mice (Burl et al., 2018; Cho et al., 2019; Hepler et al., 2018; Merrick et al., 2019; Schwalie et al., 2018) or are subject to bias from prior fluorescence-activated cell sorting (FACS) (Cho et al., 2019; Hepler et al., 2018; Hildreth et al., 2021; Merrick et al., 2019). Investigating only the SVF completely omits the transcriptional profiling of adipocytes because this lipid-laden fraction is removed in the process. Due to the technical limitations in applying scRNA-Seq to isolated adipocytes, the transcriptional profiling of adipocytes can be achieved by extracting nuclei from isolated adipocytes (Rajbhandari et al., 2019) or whole WAT (Emont et al., 2022; Sárvári et al., 2021; Sun et al., 2020) and subjecting them to single-nuclei (sn) RNA-Seq or by performing spatial transcriptomics on sectioned WAT (Bäckdahl et al., 2021).

Importantly, snRNA-Seq in whole WAT permits a comprehensive interrogation of all cell types in WAT in their natural milieu and abrogates the need for manual separation of the tissue which can lead to loss of critical cell types. To date the cellular compositions of SVF, isolated adipocytes and whole WAT have not been compared from the same samples. Recent literature suggests that scRNA-Seq and snRNA-Seq are comparable in detecting cell types despite loss of mRNA from the cytoplasm and other organelles (Bakken et al., 2018). To date, all sc/sn RNA-Seq research in adipose tissue has been performed with technology that only permits 3′ or 5′ amplification which results in lower gene coverage per cell. We leveraged a full-length SMART-Seq technology which amplifies from both 3′ and 5′ ends to interrogate the compositional differences in WAT when comparing: 1) whole adipose tissue snRNA-Seq, 2) isolated adipocytes snRNA-Seq, and 3) SVF scRNA-Seq. This approach and level of sc/snRNA-Seq benefits the capture of more mapped transcript reads (Mamanova et al., 2021), i.e. enhanced gene coverage, that has specifically been applied in this study and significantly extends previous sc- and snRNA-Seq findings in WAT.

Results

We performed subcutaneous abdominal adipose tissue biopsies on two female participants (Table 1, STAR Methods). A portion of each biopsy was immediately snap frozen for future nuclei isolation for snRNA-seq. Another portion was digested with collagenase buffer and separated by differential centrifugation into adipocytes and SVF (see STAR Methods). Nuclei were extracted from adipocytes and subjected to snRNA-seq, whereas cells from the SVF underwent scRNA-Seq. All samples were processed using the full-length SMART-Seq technology on the iCELL8 platform (Takara Bio USA, San Jose, CA) (see STAR Methods).

Table 1

Participant characteristics

	Participant A	Participant B
Age (yrs)	26	44
Weight (kg)	105.6	71.97
BMI kg/m2	39.1	27.2
Fasting Glucose (mg/dL)	85	81
Fasting insulin, μIU/mL	14.2	21
HbA1c (%)	5.3	4.8
Total cholesterol (mg/dL)	253	225
LDL (mg/dL)	160	166
HDL (mg/dL)	59	39
Triglycerides (mg/dL)	170	102
ALT (units/L)	12	15
AST (units/L)	22	18

ALT, Alanine Aminotransferase; AST, Aspartate Transaminase; BMI, Body Mass Index; HbA1c, Hemoglobin A1C; HDL, High-Density Lipoprotein; LDL, Low-Density Lipoprotein.

snRNA-Seq of whole white adipose tissue covers major cell types and highlights adipocyte heterogeneity

Integration and clustering of whole WAT samples yielded nine clusters from 2253 nuclei (Figure 1A). Stem cells (PTPRC-/PECAM1-/CD34+ /PDGFRA+ /PDGFRB+) accounted for 7% of total WAT (Figure 1B). Pre-adipocytes were defined by their increased ATXN1, ZNF423, and CD38 expression and decreased expression of stem cell markers (PDGFRA/PDGFRB/DCN) and comprised three distinct clusters accounting for 37% of total WAT (Figure 1B). Endothelial cells were enriched with PECAM1, CDH5, and VWF and accounted for 11% of total WAT (Figure 1B). The immune cluster (3%) was enriched specifically with macrophage transcripts (FCGR1A/HLA-DRB1/HLA-DPA1/ITGAX/TREM2; Figure 1B). Adipocytes were initially identified with known marker genes (DGAT2/PLIN1/LEP/ADIPOQ/PPARG/LIPE/FABP4/SAA1) and had downregulation of pre-adipocyte (ATXN1/CD38) and stem cell markers (PDGFRA/PDGRFB/DCN) and accounted for 43% of whole WAT and comprised 3 clusters. To determine how these clusters differed, we performed differential gene expression and over-representation analyses on upregulated and downregulated genes (Tables S1 and S2). Adipocyte 2 from the whole WAT had upregulated genes related to mitochondrial processes and protein targeting to the cell membrane and downregulated genes related to synaptic membrane and metal ion transmembrane transporter activity (Figures 1C and Table S2). This cluster also had the highest expressions of mature adipocyte markers LEP, ADIPOQ, PPARG, and SAA1 (Figure 1B) and likely represents fully developed adipocytes. Adipocyte 3 had upregulation of genes related to RNA splicing, glucose homeostasis, and lipid metabolic processes and downregulation of genes related to extracellular matrix and lipoprotein regulation (Figures 1C and Table S2). In addition to upregulation of glucose homeostasis genes (MLXIPL/PIK3CA/HK2), Adipocyte 3 also had upregulation of insulin signaling genes (IRS1/IRS2/PIK3CA; Figure 1B; Table S1), indicative of a glycolytic phenotype. The concept of glycolytic adipocytes was established by Digirolamo et al. (1992) and later by Lee et al. (2017) who identified TBX15-expressing pre-adipocytes that had a more glycolytic phenotype and gave rise to more glycolytic-like adipocytes in vitro (Lee et al., 2017). In agreement, we observe upregulated expression of TBX15 in Adipocyte 3 (Figure 1C), thus supporting the existence of glycolytic adipocytes in human WAT. As expected with adipocytes, some features of lipid metabolism (PLIN1/ABHD5) were retained in Adipocyte 3. Spatial transcriptomics of human subcutaneous WAT recently identified an insulin-responsive adipocyte with high PLIN1 expression in agreement with our findings (Bäckdahl et al., 2021). Adipocyte 1 was significantly enriched for GO terms; lipid droplet, synaptic vesicle, and mitochondrial inner membrane and did not have any significant GO terms that were downregulated (Figure 1C; Table S2). The mitochondrial inner membrane GO term enrichment was less significant and had less contributing genes than Adipocyte 2 indicating that it was not as oxidative as the Adipocyte 2 cluster. Adipocyte 1 had an upregulation of genes related to early stages of adipogenesis (CEBPA/STAT5A; Figure 1B; Table S1) and lipid droplet formation (SYNGR2/PLIN4/CES1; Table S1) but did not have the highest expression of the quintessential adipocyte lipid droplet marker PLIN1. The transcriptional profiles of human mesenchymal stem cells undergoing adipogenesis have recently been determined (Yi et al., 2019). To confirm that Adipocyte 1 cluster was undergoing the early stages of adipogenesis, we compared the markers defined by this recent publication (2019) (Figure S1). In this list of genes, all are increased with adipogenesis except BDNF F2R, RAC2, and RAPGEF3, whose expressions decreased with adipogenesis. In agreement, we find that Adipocyte 1 has reduced expressions of all genes that should increase during adipogenesis—except LIPE—and increased expressions of genes that should decrease during adipogenesis (RAPGEF3, RAC2, and BDNF). In contrast, Adipocyte 2 and 3 had increased expression of all “increasing” adipogenesis genes, suggesting that Adipocyte 1 included adipocytes that were less mature than adipocytes in Adipocyte clusters 2 and 3.

Figure 1

Single nuclei RNA-Seq of frozen subcutaneous white adipose tissue (WAT)

(A) UMAP showing 9 clusters from 2253 nuclei.

(B) Dotplot showing average standardized expression of differentially expressed genes that distinguish cell population determined by Wilcoxon rank-sum test.

(C) Selected gene ontology (GO) terms over-represented in the three adipocyte clusters.

(D) Pseudotime trajectory of pre-adipocytes and adipocytes mapped to the UMAP.

(E) RNA velocity analysis of pre-adipocytes and adipocytes mapped to the UMAP.

(F) Boxplots showing median and minimum and maximum quartiles of adipocyte differentiation score (ADS) in each of the pre-adipocyte and adipocyte clusters.

(G) Boxplots showing median and minimum and maximum quartiles of ADS according to octiles along the pseudotime trajectory.

Adipocyte differentiation can be tracked in vivo with pseudotime and RNA velocity analyses using snRNA-Seq of whole WAT

Differentiation of pre-adipocytes into adipocytes is a critical step in adipogenesis and adipose tissue expansion. The majority of studies investigating adipogenesis have been restricted to in vitro studies or in vivo mouse models (Farmer, 2006; Kang et al., 2012; Sarjeant and Stephens, 2012; Tang and Lane, 2012). Notably, evidence is emerging in humans evaluating the rates of adipocyte formation in vivo in obesity and in response to lifestyle interventions (Ludzki et al., 2020; White et al., 2016, 2021); however, the transcriptional regulation of adipogenesis in humans in vivo has not been thoroughly investigated. We leveraged the use of previously validated pseudotime trajectory analyses (Merrick et al., 2019) and RNA velocity (la Manno et al., 2018) to track adipocyte differentiation of the pre-adipocytes and adipocytes from whole WAT snRNA-Seq (Figures 1D and 1E). The pseudotime trajectories initiated in the pre-adipocyte clusters and continued into the adipocyte clusters, ending at Adipocyte 3 (Figure 1D). To validate the pseudotime trajectory of adipocyte differentiation, we generated an adipocyte differentiation score (ADS) by calculating the average standardized expressions of genes implicated in adipocyte differentiation (Sarjeant and Stephens, 2012) (see STAR Methods). As expected, the ADS was lower in pre-adipocytes compared to differentiated adipocytes (Figure 1E). The ADS was lower at the start of the pseudotime and increased during the trajectory, thus representing enhanced adipocyte differentiation with the pseudotime (Figure 1F). An alternative way to tracking differentiation is with RNA velocity which takes into consideration the relative abundance of unspliced nascent mRNA and mature spliced mRNA (la Manno et al., 2018). This analysis is underlined by a simple model of transcriptional dynamics that posits during transcription there is first an upregulation of unspliced mRNA followed by a subsequent increase in spliced mRNA (la Manno et al., 2018; Zeisel et al., 2011). In agreement with this model, we show that pre-adipocytes have a higher ratio of unspliced mRNA to spliced mRNA in comparison to adipocytes, indicative of them going through transcriptional dynamics to transition to adipocytes (Figure S2A). Our RNA velocity map also revealed a trajectory ending in Adipocyte cluster 3 (Figure 1E) in agreement with the pseudotime trajectory.

The average number of genes per nuclei after filtering and integration for this analysis was 2996 ± 1230 (Figure S3A), which exceeds previous analyses on whole WAT in mice using 3′ DE technology (Sárvári et al., 2021). Therefore, using full-length SMART-Seq technology, we were able to annotate the major cell types observed in human WAT with relatively few cells and profile the transcriptional regulation of adipogenesis. We next determined how the identified cell types in whole human WAT from snRNA-Seq compared to the cellular compositions of isolated adipocytes and to the SVF using snRNA-Seq and scRNA-Seq, respectively.

Adipocyte snRNA-Seq shows lack of adipocyte heterogeneity

Integration and clustering of the isolated adipocyte fraction samples revealed five clusters from 2025 nuclei (Figure 2A), which averaged 3888 ± 1118 genes per nuclei (Figure S3B). We resolved two distinct clusters of adipocytes; however, contamination from non-adipocyte cells such as stem cells (21%), pre-adipocytes (25%), and immune cells (2%) were also observed (Figures 2A and 2B). The slow centrifugation speed (200g) to separate adipocytes from the SVF is routinely used to prevent rupture of adipocytes during the isolation procedure. While keeping adipocytes intact, however, it appears that complete separation from non-adipocyte cells that seem to bind to adipocytes following collagenase digestion is not achieved. The issue can potentially be exacerbated if the WAT is fibrotic, as it will contain vast quantities of fibrotic and mast cells (Divoux et al., 2010). To confirm the presence of cell types in our different analyses, we integrated and normalized all datasets and conducted pairwise gene expression correlation analyses on each major cell/nuclei types (Figure S4). Stem cell and pre-adipocytes from the isolated adipocyte fractions were highly correlated with the stem cell and pre-adipocytes from the SVF, thus confirming their identities (Figure S4).

Figure 2

Single nuclei RNA-Seq of isolated adipocytes derived from subcutaneous abdominal white adipose tissue (WAT)

(A) UMAP showing 5 clusters from 2025 nuclei.

(B) Dotplot showing average standardized expression of differentially expressed genes that distinguish cell population determined by Wilcoxon rank-sum test.

(C) Selected gene ontology (GO) terms over-represented in the two adipocyte clusters.

The immune cell cluster had high expressions of markers of pro-inflammatory macrophages (Figure 2B). Previous reports found that obesity is associated with increased proportions of pro-inflammatory macrophages forming crown-like structures around dying adipocytes in rodents (Cinti et al., 2005; Coats et al., 2017; Murano et al., 2008) and humans (Cinti et al., 2005). The small proportion (2%) of pro-inflammatory macrophages observed in the isolated adipocyte fraction may represent macrophages actively engulfing dying adipocytes that are also extracted during the slow centrifugation step of the isolation procedure. Consistent with this finding, macrophages have been observed in a 3D analysis of floating adipocytes separated from murine WAT by collagenase digestion and fractionated by density centrifugation (Ebke et al., 2014).

To test for adipocyte heterogeneity, we performed DEG and over-representation analyses (Tables S3 and S4). Adipocyte 1 had an upregulation of genes associated with lipid metabolism (LIPE/DGAT2) and mitochondrial capacity (NDUFV1/ATP5MC3) and downregulation of genes involved with synaptic membrane and metal ion transmembrane transporter activity (Figure 2C; Tables S3 and S4). Adipocyte 2 had an upregulation of genes related to co-translational protein targeting to membrane (RPL27/RPL29) which encode for ribosomal subunits (Figure 2C; Tables S3 and S4) and may be indicative of protein synthesis (Wolins et al., 2006) and growth. Adipocyte 2 had a downregulation of genes related to mitochondria and fatty acid metabolic process, in opposition to Adipocyte 1 (Figures 2C and Table S4). To date, no one has transcriptionally profiled isolated adipocytes from human WAT. Previous studies using adipocyte nuclei derived from inguinal WAT in mice have identified oxidative adipocytes but did not observe adipocytes with high amounts of ribosomal-encoding genes (Rajbhandari et al., 2019). Adipocyte 1 and 2 in the isolated adipocyte snRNA-Seq analysis both had similar enrichment to Adipocyte 2 from the whole WAT snRNA-Seq analysis, suggesting that both adipocyte clusters present in the adipocyte fraction are grouped together as Adipocyte 2 in the whole WAT. The adipocyte fraction was void of adipocytes that were more glycolytic or undergoing the early stages of adipogenesis, suggesting that more adipocyte heterogeneity is obtained with whole WAT analyses compared to isolated adipocyte fractions.

scRNA-Seq of SVF provides greater resolution of non-adipogenic cells

Integration and clustering of the two unsorted SVF scRNA-Seq samples revealed 16 clusters (Figure 3A) from 1776 cells, which were annotated using known gene expression markers (Figure 3B) (Briot et al., 2018; Ehrlund et al., 2017; Smyth et al., 2018; Tran et al., 2012; Trim et al., 2022; Vijay et al., 2020). The average number of genes per cell was 3201 ± 1450 (Figure S3C). Stem cells (PTPRC-/PECAM1-/CD34+ /PDGFRA+ /PDGFRB+) accounted for 28% of the total SVF. Stem cell 2 was subcategorized as hematopoietic-derived stem cells owing to higher expressions of CD34, CD59, and THY1 and a lower expression of PDGFRA (Gao et al., 2017) compared to stem cell 1. Stem cell 3 had upregulation of transcripts related to very early signs of adipogenesis (WISP2/SFRP2; Table S5) and a lower expression of CD34 (Figure 3B). Committed pre-adipocytes (9%) were identified by ATXN1, ZNF423, and CD38 and downregulation of stem cell markers (PDGFRA/PDGFRB/DCN) (Figure 3B). Differential expressions of DLK1, DGAT2, and CIDEC between the subclusters Pre-Ad 1 and Pre-Ad 2 suggested differing degrees of commitment along the adipogenic lineage. The SVF scRNA-Seq offered greater resolution and heterogeneity of stem cells and pre-adipocytes in comparison to whole WAT snRNA-Seq.

Figure 3

Single cell RNA-Seq of the stromal vascular fraction derived from subcutaneous abdominal white adipose tissue (WAT)

(A) UMAP plot showing 16 clusters from 1776 cells.

(B) Dotplot showing average standardized expression of differentially expressed genes that distinguish the cell populations determined by Wilcoxon rank-sum test.

(C) Pseudotime trajectory of stem cells and pre-adipocyte mapped onto SVF UMAP.

(D) RNA velocity analysis of the stem cells and pre-adipocytes mapped to the UMAP.

(E) Heatmap of aggregated expression of co-regulated genes for each module mapped to each cell-type from the pseudotime trajectory.

(F) Single-cell expression trajectories of highlighted genes along the pseudotime trajectory.

Endothelial cells were identified with expressions of VWF, PECAM1, and CDH5 and accounted for 16% of total cells (Figure 3A). In agreement with previous work (Vijay et al., 2020), we subclassified endothelial cluster 2 as fatty acid (FA) handling endothelial cells due to high expressions of FA handling transcripts (FAPB4/CD36/RGS5; Figure 3B), frequently associated with the microvasculature of WAT (Briot et al., 2018). Endothelial cluster 3 also showed high expressions of stem cell markers (PDGFRA/PDGRB) commonly observed with pericytes (Smyth et al., 2018). The endothelial cell proportions reported herein are consistent with a previous report (Trim et al., 2022) but slightly higher in comparison to other scRNA-Seq analyses in unsorted SVF (Vijay et al., 2020) and FACS-sorted SVF proportions (Ehrlund et al., 2017).

Immune cells (PTPRC+) comprised the majority of the remaining cell populations (38%) and were further subclassified as mast cells (Immune 1: KIT/TPSAB1/CMA1/CTSG), NK/T cells (Immune 2: NKG7/GNLY/CD8A/CD8B/CD3E), resident macrophages (Immune 4 and 5: CD14/LGMN/HLA-DP1/CD8/, CD163) and pro-inflammatory macrophages (Immune 6: ITGAX/HLA-DP1/CD86; Immune 7: FCGR1A/ITGAX/TREM2). Immune cell 3 had upregulated expressions of macrophages markers (LYVE1/CCL8/CCL4L4/CCL2), in addition to fibroblast markers (PDGFRA/PDGFRB) suggesting that these cells may be transitioning from hematopoietic-derived stem cells.

Cluster 4 had upregulated expressions of adipocyte-specific markers such as PLIN1 and vasculature markers such as VEGFB. A recent report demonstrated that endothelial cells may be capable of forming new pre-adipocytes that undergo differentiation (Tran et al., 2012). It is also plausible that these cells are adipocytes that have not generated sufficient lipid to be float during the isolation procedure and remain in the SVF. We consider that the cells identified in cluster 4 are therefore a mixture of early differentiating adipocytes or endothelial cells that are transcriptionally similar to each other in this dataset.

Commitment of mesenchymal stem cells to the adipogenic lineage can be tracked in vivo with pseudotime analysis using scRNA-Seq of SVF

Commitment of mesenchymal stem cells to the adipogenic lineage is an initial step in adipogenesis (Cawthorn et al., 2012); however, the transcriptional regulation of this process in humans in vivo has not been thoroughly investigated. Due to Pre-Adipocyte 2 having higher expressions of committed pre-adipocyte markers (DLK1/DGAT2/CIDEC), we hypothesized that the pre-adipocyte clusters differ by their degree of commitment to the adipogenic lineage. We leveraged the use of previously validated pseudotime trajectory analyses (Merrick et al., 2019) and RNA velocity (la Manno et al., 2018) to interrogate the transcriptional pathways of stem cells and committed pre-adipocytes within the SVF. A unidirectional pseudotime trajectory was produced, initiating from the stem 1 and terminating at Pre-Adipocyte 2 (Figure 3C). Mapping of RNA velocity revealed a similar directional pattern from stem cell to pre-adipocytes with RNA velocity being greatest throughout the pre-adipocyte populations (Figure 3D). Pre-adipocytes had a greater proportion of unspliced to spliced RNA in comparison to the stem cell populations indicating pre-adipocytes are going through an active transcriptional transition phase (Figure S1B). Module-based analyses were performed to summarize co-regulated mRNA transcripts throughout the pseudotime trajectory revealing 14 different modules (Figure 3E; Table S6). Module 9 (Table S6) was enriched with extracellular matrix genes (ECM-DCN/LUM/ANXA1/ANXA2) whose expressions decreased throughout the trajectory (Figure 3F), in line with previous findings (Zhang et al., 2007). We then identified a module that increased toward pre-adipocyte commitment (Module 7; Table S6) that contains the well-established pre-adipocyte marker DLK1 (Figure 3B). Transcription factors FOXA3 and ISL1—which have been implicated in adipocyte differentiation (Ma et al., 2014; Xu et al., 2013)—were among the top 25 top ranked genes in Module 7. FOXA3 cooperates with CEBPB and CEBPD to transcriptionally induce PPARG expression (Xu et al., 2013). In contrast, ISL1 reduces adipocyte differentiation by inhibiting PPARG expression through downregulation of BMP4 (Ma et al., 2014). PPARG was not included in any of the modules indicating that its change over the trajectory was not in concert with other genes. Mapping of PPARG revealed a steady expression throughout the trajectory (Figure 3F), confirming these pre-adipocytes have not yet initiated differentiation. Module 7, therefore, encapsulates a cohort of genes indicative of pre-adipocyte commitment culminating in Pre-Adipocyte 2 that is primed, but not yet triggered, toward adipocyte differentiation.

Discussion

We performed a comprehensive interrogation of sc/snRNA-Seq in whole WAT, isolated adipocytes, and SVF using in humans using full-length SMART-Seq technology. The benefit of using the whole WAT rather than isolated adipocytes is that it abrogates the need for manual separation of the tissue which can lead to loss of critical cell types, while also allowing for the potential to utilize archived samples from previous clinical studies. We captured all of the major cell populations including endothelial cells, immune cells, stem cells, pre-adipocytes, and adipocytes. Using the full-length SMART-Seq technology, we achieved robust gene coverage from the 5′ to 3′ amplification in each type of analysis. Each of our analyses exceeded the number of genes per cell/nuclei that have been seen in single cell adipose tissue research to date due to the depth of coverage for each cell (Emont et al., 2022; Sárvári et al., 2021; Vijay et al., 2020). Enhanced gene coverage permits more accurate profiles of the transcriptome and identification of cell types within the adipose niche. We posit that these detailed information are beneficial for characterizing the cellular compositions of WAT and for future identification of where potential alterations occur that impact cellular process of WAT metabolism. Interestingly, the number of genes per nuclei was similar to that observed in the SVF scRNA-Seq, thus indicating that robust cell identification can be also achieved at the nuclei level.

We were able to identify three distinct adipocyte clusters from the whole WAT analyses that were defined as 1) newly formed, 2) lipolytic/oxidative, and 3) glycolytic and in line with previous research (Bäckdahl et al., 2021; Lee et al., 2017, 2019). These three distinct adipocyte clusters were not identified in the snRNA-seq of the isolated adipocyte fractions. In the isolated adipocyte fractions, we resolved two adipocyte clusters that were a mixture of Adipocyte 2 cluster in the whole WAT RNA-Seq analyses. In the isolated adipocyte snRNA-Seq, there was also contamination from stem cells, pre-adipocytes, and immune cells. These contaminating cell fractions highlight a potential limitation of this routine isolation procedure. Due to the contamination and selection bias of adipocytes, we suggest the transcriptional regulation of adipocytes is more amply measured using whole WAT rather than prior isolation. Furthermore, we were able to perform trajectory analyses and RNA velocity of pre-adipocytes and adipocytes in the WAT to assess the transcriptional profile of adipocyte differentiation, highlighting a critical advantage of profiling all cell types together within the same tissue. Samples with varying BMI were used to encapsulate greater heterogeneity of cells. Increasing the nuclei number and diversifying the metabolic phenotype of the donors will likely enable sub-categorization of non-adipocyte cells. Increasing the number of samples with a lean and obese phenotype would permit analysis on how gene expression is potentially altered in cell subtypes with obesity.

Full-length scRNA-Seq on SVF provided robust data to deeply interrogate heterogeneity among immune cells, endothelial, pre-adipocytes, and stem cells, while profiling relatively few cells in comparison to previous literature (Vijay et al., 2020). The same resolution was not achieved in isolated adipocyte snRNA-Seq and whole WAT snRNA-Seq. The deep transcriptional profiling of SVF scRNA-Seq permitted analyses into the transcriptional regulation of stem cell to pre-adipocyte commitment in human SVF in vivo which is useful for researchers investigating the early stages of adipogenesis and the origins of pre-adipocytes (i.e. mesenchymal stem cell, endothelial cells, and/or pericytes).

RNA velocity analyses in both whole WAT snRNA-Seq and SVF scRNA-Seq showed pre-adipocytes have a greater abundance of unspliced mRNA indicating they are in transition. In the SVF, this transition period marks the commitment to the adipogenic lineage as committed pre-adipocytes. Whereas in the whole WAT, this transitional period is marking differentiation to adipocytes. It was notable that while pre-adipocytes consistently had greater abundance of unspliced mRNA in both analyses, there was a discordance in their abundance between analyses with the SVF pre-adipocytes showing greater unspliced mRNA in comparison to the whole WAT pre-adipocytes. The RNA velocity model was designed for single cell data under the assumptions there is constant RNA degradation and nuclear transport (Bergen et al., 2021). While RNA velocity has shown promising results in single nuclei data (Bergen et al., 2021; Marsh and Blelloch, 2020), the assumption of these models has not been conclusively verified (Bergen et al., 2021). Thus, RNA velocity shows similar patterns between scRNA-Seq and snRNA-Seq, but their data should not be compared between analyses. Nevertheless, the RNA velocity confirmed our pseudotime trajectory showing commitment to pre-adipocytes and differentiation of pre-adipocytes to adipocytes. Furthermore, we highlight that pre-adipocytes differ from other adipose progenitor cells as they are transitioning through the adipogenic trajectory and are marked by higher ratios of unspliced RNA.

In conclusion, we highlight that the cellular compositions and transcriptional profiles of adipogenesis can be evaluated with the use of full-length SMART-Seq sc technology when using nuclei extracted from frozen whole human WAT. Furthermore, non-adipogenic cell diversity can be thoroughly explored in human SVF with full-length SMART-Seq single cell technology while using relatively few cells.

Limitations of the study

The manuscript assesses the cellular composition of WAT when analyzed with different experimental techniques and highlight advantages used with snRNA-Seq on WAT and scRNA-Seq on SVF. We used participants with varying BMI to encapsulate varying cell types for this analysis. However, we were unable to perform analyses comparing lean and obese individuals. There may be more diverse cell populations present with varying metabolic diseases that were not profiled in the current study. Our focus was on subcutaneous abdominal adipose tissue as one of the most researched depots in human metabolic health. Other adipose tissue depots in humans such as gluteal-femoral subcutaneous adipose tissue and visceral adipose tissue may have different cellular compositions compared to what is identified in this manuscript. Lastly, we only used females in this current analyses and work should be done to extrapolate these findings to males.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Dr. Lauren Sparks (Lauren.Sparks@adventhealth.com).

Materials availability

The study did not generate new unique reagents.

Experimental model and subject details

Participants and tissue collection

Two healthy eumenorrheic, pre-menopausal females with no adverse cardio-metabolic disorders (Table 1) were recruited to the Translational Research Institute at AdventHealth to participate in the study. This study was approved by AdventHealth Institutional Review Board and carried out in accordance with the Declaration of Helsinki. Participants provided written informed consent for their study participation.

Subcutaneous abdominal adipose tissue biopsies were performed following an overnight fast (Divoux et al., 2018), using the tumescent lidocaine approach with a Mercedes aspiration cannula. Following removal of excess blood and connective tissue, the sample was cleaned with PBS. A portion (∼100mg) was immediately snap frozen for subsequent nuclei isolation for snRNA-seq. The remaining tissue was digested for isolation of SVF and adipocytes (see below).

Method details

SVF isolation and adipocyte collection

Tissue was digested with collagenase buffer (Type I collagenase in αMEM supplemented with 1%BSA) at 37°C for 45 min. Adipocytes were collected by centrifugation at 200g for 5 minutes and washed twice with PBS prior to nuclei isolation (see below). After the adipocytes were collected, the SVF was isolated by additional centrifugation at 500g for 5 minutes and washed with Red Blood Cell lysis buffer. The SVF solution was subsequently filtered through 100μm and 40μm strainer (BD Falcon) and counted with a countess II automated cell counter (Thermofisher Scientific, Waltham, MA). Cells were stained with propidium iodine and Hoechst 33342 (ReadyProbes Cell Viability Imaging Kit, Thermofisher Scientific). Following a wash, cells were resuspended in 1% BSA- nuclease free-water.

Nuclei isolation from adipocytes

Washed adipocytes (∼150μL) were homogenized in a glass dounce with 1ml of homogenization buffer (5mM MgCl2, 25mM Tris Buffer pH 8.0, 25mM KCL, 250mM sucrose, 1μm DTT, 1 x protease inhibitor, 0.2 U/μL SUPERase · In RNase Inhibitor (Thermofisher Scientific) in nuclease-free water). Triton-X100 (0.1% v/v) was added to the homogenate prior to incubation on ice for 30 min with regular vortexing. Samples were filtered through a 100μm and 40μm strainer (BD Falcon), centrifuged at 1,000g for 10 min at 4°C, resuspended in 1mL nuclei isolation medium (5mM MgCl2, 25mM Tris Buffer pH 8.0, 25mM KCL, 1 mM EDTA, 0.2 U/μL Ribolock RNAase inhibitor, 1% BSA in nuclease-free water) and centrifuged again at 1,000g for 10 min at 4°C. Nuclei were resuspended in 500μL nuclei isolation medium and filtered 10 x with a 25g syringe prior to being stained with Hoechst 33342 (ReadyProbes Cell Viability Imaging Kit, Thermofisher Scientific). Nuclei were subsequently counted with a countess II automated cell counter (Thermofisher Scientific).

Nuclei isolation from whole adipose tissue

Frozen whole adipose tissue (∼100mg) was pulverized under liquid nitrogen and homogenized in 2mL of homogenization buffer (5mM MgCl2, 25mM Tris Buffer pH 8.0, 25mM KCL, 250mM sucrose, 1μm DDT, 1 x protease inhibitor, 0.2 U/μL SUPERase · In RNase Inhibitor (Thermofisher Scientific) in nuclease-free water) with a glasscol homogenizer. The homogenate was incubated for 30 minutes on ice with the addition of Triton-X100 (0.1% v/v) with regular vortexing. Samples were filtered through a 100μm and 40μm strainer (BD Falcon), centrifuged at 2,700g for 10 min at 4°C, resuspended in homogenization buffer and recentrifuged again at 2,700g for 10 min at 4°C. The pellet was then re-suspended in 1mL nuclei isolation medium (5mM MgCl2, 25mM Tris Buffer pH 8.0, 25mM KCL, 1 mM EDTA, 0.2 U/μL Ribolock RNAase inhibitor, 1% BSA in nuclease-free water) and centrifuged at 2,700g for 10 min at 4°C. Following re-suspension in 500μL nuclei isolation medium, sample was filtered 10 x with a 25g syringe. Nuclei was stained with Hoechst 33342 (ReadyProbes Cell Viability Imaging Kit, Thermofisher Scientific) prior to counting with a countess II automated cell counter (Thermofisher Scientific).

Single-cell and single-nuclei RNA-Seq

Single-cell suspension (28L/mL) or single-nuclei suspension (40K/mL) was distributed into eight wells of a 384-well source plate (Takara Bio USA, San Jose, CA) and dispensed onto a iCELL8® 350v Chip (Takara Bio USA) using an iCELL8 MultiSample NanoDispenser (Takara Bio USA). Following dispense, the chip nanowells were imaged using the iCELL8 Imaging Station to identify nanowells containing a single nucleus/ a live single cell. Only these nanowells were subjected for downstream dispenses. After imaging, the chip was subjected to freeze-thaw to lyse the cells/nuclei, followed by a 3 min incubation at 72°C to denature the RNA. Selected cells were subjected to first-strand cDNA synthesis initiated by oligo dT primer (SMART-Seq iCELL8 CDS), followed by template switching with template switching oligo (SMART-Seq iCELL8 oligonucleotide) for 2nd strand cDNA synthesis, prior to unbiased amplification of full-length cDNA. Synthesized full-length cDNA was tagmented by tagment DNA enzyme 1 (TDE1, Illumina, San Diego, CA) and amplified using forward (i5) and reverse (i7) indexing primers. Each single nucleus/cell was indexed by a unique combination of 1 of 72 forward and 1 of 72 reverse indexing primers allowing for downstream identification. Collected cDNA was purified twice using a 1:1 proportion of AMPure XP beads (Beckman Coulter, Brea, CA), further amplified according to manufacturer’s instructions and purified again at a 1:1 proportion of AMPure XP beads. The resultant sequencing-ready cDNA library was assessed for concentration by fluorometer (Qubit, Thermofisher Scientific) and quality by electrophoresis (Agilent Bioanalyzer high sensitivity DNA chips). Libraries were sequenced with Illumina HiSeq 4000 to produce between 226 M and 241 M barcoded reads per library.

Quantification and statistical analysis

Initial analyses of the single cell and single nuclei libraries were performed using Mappa™ Analysis Pipeline (Takara Bio, USA). GRCh38 was used as the genome reference. Cell and gene filtering was performed in R package scran (Lun et al., 2016). Briefly, cells were filtered if they had; < 200 genes, < 10,000 read counts, > 20% mitochondrial reads or if cell complexity was < 0.8. A minimum threshold of 0.1 was applied to filter out low expressed genes. Imputation was performed with scImpute to account for dropout values (Li and Li, 2018). Sample integration and clustering was performed in Seurat (Hao et al., 2021). Samples were integrated based on analysis. Specifically, two SVF scRNA-Seq samples were integrated together, two adipocyte snRNA-Seq samples were integrated together and two adipose tissue snRNA-Seq samples were integrated together. Imputed counts were normalized with SCTransform (Hafemeister and Satija, 2019), while regressing out variation due to, cell cycle and mitochondrial gene expression. Once integrated, significant principal components were used to perform unsupervised K-nearest neighbor (KNN) graph-based clustering. Visualization was achieved with uniform manifold approximation and projection (UMAP). Differential expression analysis was performed using a Wilcoxon rank sum test with Seurat’s “FindMarkers” function with a FDR cut off of < 0.05, log2 FC > 0.25 or < -0.25 and expressed in >25% of cells/nuclei in that cluster. A hypergeometric test was used to assess over-representation of upregulated genes (log2 FC > 0.5) and downregulated genes (log2 FC < -0.5) from clusters using the clusterProfiler R package and the Gene -Ontology (GO) and KEGG database (Yu et al., 2012). Adipocyte differentiation score was calculated using the AddModuleScore() function in Seurat for genes KLF5, KLF4, KLF6, KLF15, CEBPA, CEBPB, CEBPD, PPARG, STAT5A, STAT5B & SREBF1. This function calculates the average expression of the program of genes subtracted by aggregated expression of control feature sets. Statistical details of experiments can be found in figure legends.

Pseudotime trajectory

Monocle 3 was used for trajectory analysis (Cao et al., 2019; Levine et al., 2015; Qiu et al., 2017; Traag et al., 2019; Trapnell et al., 2014). After analyzing the dynamic biological changes of each cell, an individual position of every single cell is plotted in a learned trajectory. Based on the clustered annotation and marker genes, we identified the root of the given trajectory and ordered cells along the pseudotime according to their developmental progress. For the SVF analysis we identified the root as Stem 1 and for the whole WAT we identified the root as Pre-Ad 1. Lastly, for SVF analyses we collected genes that change over the pseudotime and grouped them into gene modules.

RNA velocity

FastQ files were ran through cogentAP™ analysis pipeline (Takara Bio, USA) to generate an individual BAM file for each single cell/nuclei. Velocyto run_smart-Seq2 command was performed to annotate spliced and unspliced reads for each cell and to generates a sample level loom file (la Manno et al., 2018). RNA velocity was performed on the selected cell populations with a PCA reduction using scVelo (Bergen et al., 2020). 25 of the nearest neighbors were selected for slope calculation smoothing and gamma fit was performed on the top/bottom 2% quantiles of expression magnitudes. UMAP embeddings from the previous Seurat objects were used for visualization.

Pairwise correlation analysis

All 6 samples (2 SVF scRNA-seq, 2 adipocyte snRNA-Seq and 2 AT snRNA-Seq) were also integrated to allow comparison between analyses. Barcodes of cells/nuclei comprising the main cell types from each analysis were used to define the main cell populations. Correlation analyses were performed among these different cell populations.

Additional resources

Clinical trial registration: NCT04034706.

Reagent or Resource	Source	Identifier
Chemicals, peptides and recombinant proteins
Type I collagenase	Worthington	M2C13334
αMEM	Gibco	32561-037
BSA	Sigma-Aldrich	820452
Red blood cell lysis buffer	BioLegend	420301
Ready proves cell viability imaging kit	Thermo Fisher Scientific	R37610
MgCl2	Ambion	AM9530G
Tris Buffer pH 8.0	Thermo Fisher Scientific	AM9855G
KCL	ThermoFisher Scientific	AM9640G
Sucrose	Sigma-Aldrich	50389
DTT	Thermo Fisher Scientific	R0861
100x Protease inhibitor	Thermo Fisher	78437
SUPERaseIn RNase Inhibitor	Thermo Fisher Scientific	AM2695
Triton-X100	Fisher Scientific	AC327372500
Ribolock RNAse inhibitor	Thermo Fisher Scientific	EO0382
UltraPure™ 0.5M EDTA, pH 8.0	Gibco	15575020
Tagment DNA enzyme 1	Illumina	20034198
Beckman Coulter AMPURE XP KIT	Fisher Healthcare	NC9959336
Critical commercial assays
SMART-Seq® ICELL8® Application Kit – 5 Chip	TakaraBio USA	640221
Qubit™ 1X dsDNA Assay Kits, high sensitivity (HS) and broad range (BR)	Invitrogen	Q33230
High sensitivity DNA kit	Agilent	5067-4626
Deposited data
Raw and mapped single cell/ single nuclei RNAseq data	This paper	GSE189346
Software and algorithms
Mappa™ Analysis Pipeline	TakaraBio USA	https://www.takarabio.com/documents/User%20Manual/mappa%20Analysis%20Pipeline%20v1.0%20User%20Guide/mappa%20Analysis%20Pipeline%20v1.0%20User%20Guide_120319.pdf
Scran	(Lun et al., 2016)	https://bioconductor.org/packages/release/bioc/html/scran.html
SCImpute	(Li and Li, 2018)	https://github.com/Vivianstats/scImpute
Seurat	(Hao et al., 2021)	https://github.com/satijalab/seurat/
SCTransform	(Hafemeister and Satija, 2019)	https://github.com/satijalab/sctransform
Over-representation analysis with Gene-Ontology terms: clusterProfiler	(Yu et al., 2012)	https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
Monocle 3	(Cao et al., 2019; Levine et al., 2015; Qiu et al., 2017; Traag et al., 2019; Trapnell et al., 2014)	https://cole-trapnell-lab.github.io/monocle3/
Correlation analysis: Hmisc	https://www.rdocumentation.org/packages/Hmisc/versions/4.6-0	https://www.rdocumentation.org/packages/Hmisc/versions/4.6-0
Velocyto	Velocyto.org	Velocyto.org
scVelo	http://scvelo.readthedocs.io/#	http://scvelo.readthedocs.io/