Introduction

Single cell RNA sequencing (scRNA-seq) has revolutionized various fields of biology by enabling high-resolution profiling of single-cell heterogeneity within complex tissues¹. It has not only facilitated the identification of new cell types and subtypes², but also uncovered alterations in cellular composition and gene expression across diverse biological contexts during development, disease progression and responses to environmental stimuli^3–5. However, experimental methods for scRNA-seq still present limitations and challenges. First, these methods require the use of fresh tissues, which can be difficult to obtain and thus lead to logistic complexities^6–9. Even when fresh tissues are available, the process of obtaining single cell suspension through enzymatic digestion and mechanical dissociation can introduce unwanted artifacts. Different types of cells can be digested to varying degrees and rates, resulting in biased cell recovery. Also, the tissue dissociation process can cause cellular stress, leading to reduced cell viability and altered transcriptional profiles. Lastly, certain cell types are unsuitable for scRNA-seq due to their unique cellular properties^10,11. Adipocytes are one example; their large size and excess lipid content render them incompatible with the microfluidic platforms used in the majority of scRNA-seq technologies.

To overcome these issues, single nucleus RNA sequencing (snRNA-seq) has been developed as an alternative strategy. While generating transcriptional profiles with high similarity to those by scRNA-seq^10,12, snRNA-seq offers advantages. First, this method utilizes nuclei that are promptly released from tissues during homogenization, obviating the need for an extended digestion process and thereby minimizing cell type biases and artifactual transcriptional alteration. Second, its compatibility with cryopreserved samples allows for greater flexibility in experimental designs and enables the utilization of samples stored in biobanks. Furthermore, snRNA-seq can be employed to study cell types that are less amenable to scRNA-seq, such as adipocytes and cardiomyocytes, as isolated nuclei no longer retain limiting cellular characteristics and can be loaded on to various platforms. Nevertheless, snRNA-seq also poses challenges to be considered. Once released from cells, nuclei are highly fragile and susceptible to mechanical stress, resulting in leakage or aggregation. Moreover, their enclosed RNA becomes increasingly vulnerable to RNases, leading to rapid RNA degradation. Nuclei isolated from in vivo tissues are typically more delicate than those from cells cultured in vitro, with potential variation depending on cell type. Hence, meticulous preparation of high-quality nuclei is crucial for successful snRNA-seq experiments.

Adipose tissue plays a pivotal role in energy metabolism, serving as a dynamic energy reservoir within the body. It exhibits remarkable plasticity in response to nutritional fluctuations¹³. In obesity, adipose tissue undergoes a significant remodeling process, expanding substantially to enhance its storage capacity. This process is orchestrated through coordinated actions among diverse cell types, including adipocytes, preadipocytes, endothelial cells, and various immune cells¹⁴. However, prolonged obesity leads to pathological changes and dysfunction in adipose tissue, characterized by disrupted lipid metabolism, impaired insulin signaling, inflammation and cell death. These changes are critical precursors to the onset of type 2 diabetes and other metabolic diseases¹⁵. The distribution of body fat in obesity is closely correlated with metabolic health. The expansion of visceral fat depots increases risks for metabolic disorders, whereas subcutaneous fat depots pose lesser risks or can even offer protection¹⁶. This is likely due to varying responses of different adipose tissue depots to obesity¹⁷. In obesity, visceral adipose tissues display enlarged hypertrophic adipocytes, along with hypoxia, elevated fibrosis and inflammation from increased immune cell infiltration. In contrast, subcutaneous adipose tissues have smaller adipocytes with lower levels of inflammation and hypoxia¹⁵. scRNA-seq has been used to characterize the dynamics of cell populations underlying adipose tissue remodeling, with a focus on the cell types within stromal vascular fractions, particularly preadipocytes^18–21. Yet adipocytes have been excluded from these studies due to technical limitations. More recently, snRNA-seq has been employed, revealing evidence of the cellular heterogeneity within adipocyte populations^22–24. Despite this progress, acquiring high-quality snRNA-seq data from adipose tissues remains highly challenging. Technical difficulties have posed a major hurdle to achieving accuracy, sensitivity, and efficacy in understanding adipose tissue remodeling, especially regarding mature adipocyte populations during obesity.

Here, we present a robust method for isolating nuclei for snRNA-seq, which enhances the quality of nuclear RNA and preserves nucleus integrity across various tissue types. Using this method, we have generated snRNA-seq data of superior quality compared to existing datasets, offering novel and precise insights into the specific changes within cell populations across different fat depots in response to obesity. By integrating bulk nuclear RNA-seq with adipocyte nuclei of different sizes, we have identified multiple subpopulations of adipocytes primarily categorized by size and functionality, which are influenced by obesity. Specifically, we have characterized dysfunctional hypertrophic adipocytes prevalent in visceral adipose tissues during obesity, exhibiting hallmark features of cellular stress and inflammation. All these adipocyte subpopulations contributed differentially to a range of biological pathways involved in the pathophysiology of adipose tissue during obesity.

Results

Robust nucleus isolation protocols protect RNA quality and nucleus integrity

To systemically address the issue of nuclear RNA degradation, we first assessed RNA quality in isolated nuclei from a range of mouse tissues. During nucleus isolation, no RNase inhibitors were added to evaluate intrinsic RNase activity in each tissue type. We observed a wide spectrum of RNA quality among the tissues examined (Figure 1A). The brain exhibited highly intact RNA, followed by the heart showing relatively good RNA quality. The liver displayed moderate RNA degradation. Among the adipose tissues, including brown adipose tissue (BAT), subcutaneous inguinal white adipose tissue (iWAT), and visceral epididymal white adipose tissue (eWAT), significant RNA degradation was detected, with eWAT showing the highest level of degradation. RNA was completely degraded in the pancreas (Figure 1A). These results indicate significant variation in the quality of RNA extracted from nuclei depending on tissue types. To limit RNA degradation, we incorporated recombinant RNase inhibitors during nucleus isolation and examined additional tissue types, including the spleen, the lung, and the kidney. Surprisingly, the addition of recombinant RNase inhibitors at a concentration of 0.4 U/µl had no discernible effect on RNA quality of any of those tissues (Figure 1B). Even when the concentration of RNase inhibitors was increased to 4 U/µl, it failed to yield any substantial enhancement in RNA quality (data not shown). Therefore, we screened various conditions and molecules known for their ability to inhibit RNase activity. Among these, the vanadyl ribonucleoside complex (VRC)²⁵ was identified as having a significant impact on RNA quality. High-quality RNA with minimal degradation was obtained in all the tested tissues except for the pancreas and its connecting tissue, the spleen (Figure 1C). Finally, we optimized the protocol by combining VRC with recombinant RNase inhibitors for further improvements, albeit to a modest extent (Figure 1D). In addition to RNA quality, achieving a single nucleus suspension is an important requirement for snRNA-seq sample preparation. Using VRC, the optimized protocol successfully retrieved a well-dispersed nucleus population without clumps in all tissue types (Figure 1E). Further analysis of nucleus morphology at a higher resolution revealed a smoother membrane surface and larger nucleus size when using the optimized protocol, as compared to the standard protocol (Figure 1F). Consistent with the preservative function of VRC proven in maintaining tissue morphology²⁶, these findings show that our nucleus isolation protocol effectively maintains nucleus structure without shrinkage. Given that current snRNAseq methods require prompt sample processing due to the rapid degradation of RNA in isolated nuclei, we further tested the temporal stability of RNA in nuclei isolated using our optimized protocol. After isolating nuclei from eWAT, we stored them at 4°C for different durations (2, 6 and 24 hours) before RNA extraction and quality analysis. In the standard protocol, nuclear RNA exhibited severe degradation within 2 hours post-homogenization, worsening over 24 hours. In contrast, nuclear RNA isolated using the optimized protocol remained intact at 2 and 6 hours, maintaining good quality even after 24 hours (Figure 1G). To determine whether RNA could remain intact with additional handling procedures, we conducted the same time-course RNA extraction with nuclei sorted through flow cytometry. The RNA remained stable, comparable to that of unsorted nuclei (Figure 1G). Taken together, our optimized protocol, employing both VRC and recombinant RNase inhibitors, yields high-quality nuclei while preserving RNA integrity and nuclear structure for extended procedures up to 24 hours.

Figure 1.

The optimized protocol improves the data quality of snRNA-seq

Considering the previous report on the inhibitory activity of VRC on reverse transcriptases²⁶, we investigated whether the inclusion of VRC could confound the accurate representation of RNA expression profiles within nuclei. By utilizing the Nuclear tagging and Translating Ribosome Affinity Purification (NuTRAP) mice designed to label nuclei from specific cell types²⁷, we sorted nuclei from adipocytes and non-adipocytes to conduct a comparative analysis of their gene expression profiles through cDNA synthesis followed by quantitative real time PCR (qRT-PCR). We observed strong PCR amplification and highly specific expression of adipocyte markers in adipocyte nuclei and immune/endothelial cell markers in non-adipocyte nuclei (Figure S1A). These results indicate that the use of VRC during nucleus isolation does not pose an issue in gene expression profiling.

Subsequently, we applied the optimized nucleus isolation protocol for snRNA-seq of mouse adipose tissues using the 10X Genomics Chromium Single Cell 3’ Gene Expression system. We employed a two-step tissue homogenization process involving pulverization and Dounce homogenizers to ensure complete liberation of nuclei (Figure 2A). Additionally, we included a fluorescence-activated nucleus sorting (FANS) step in the workflow to eliminate tissue debris and ambient RNA (Figure S1B). From eWAT and iWAT samples collected from mice fed either a chow or high fat diet (HFD), we generated high-quality cDNA with long fragment sizes and constructed robust sequencing libraries (Figure S1C, S1D). To evaluate our data quality, we benchmarked them against two published snRNA-seq datasets from mouse adipose tissues^22,23, based on the quality metrics determined by the standard Cell Ranger pipeline. With comparable sequencing depths (Figure 2B), our datasets displayed significantly higher counts of Unique Molecular Identifiers (UMIs) and genes per nucleus, exceeding counts from the other datasets by almost two-fold (Figure 2C, 2D). Also, the proportions of mitochondrial reads in our datasets represented less than 1% (mean 0.37%), which was significantly lower than in the other datasets (mean 38.7 and 2.31%; Figure 2E). Importantly, our method produced substantially higher fractions of reads in nuclei (Figure 2F), suggesting low ambient RNA contamination. This was consistent with the Barcode Rank Plot for our data, which exhibited a sudden drop of UMI counts (Figure 2G), indicating a clear distinction between intact nuclei and background barcodes during cell calling. Analysis using CellBender²⁸ also contrasted the minimal number of cells in the narrow transition zone in our datasets with the greater number of cells within the broader transition zone observed in the published samples (Figure S1E). These results indicate that our method reduces the contribution of ambient RNA to the expression profiles of single cells. Collectively, these findings suggest that our optimized protocol produces high-quality snRNA-seq data that outperforms other existing methods.

Figure 2.

A single nucleus atlas of different adipose tissues from lean and obese mice

To characterize the cellular dynamics within mouse adipose tissues during obesity, we integrated data from both eWAT and iWAT of chow- or HFD-fed wild-type male mice. After filtering nuclei based on the detected numbers of UMIs and genes, and subsequently removing doublets and ambient RNA, a total of 53,395 cells/nuclei were identified, distributed across 10 distinct cell types (Figure 2H). The major cell types present in adipose tissues included adipocytes, adipocyte progenitor cells (APC), macrophages/monocytes, and endothelial cells. Additionally, we identified less abundant cell types, such as pericytes (PC)/smooth muscle cells (SMC), mesothelial cells, dendritic cells, mast cells, T cells, and B cells. We determined distinct marker gene sets for each cell type, including both established and novel markers (Figure 2I, S2A). The proportions of these cell types in our data largely matched those reported in previous studies, except for the notable improvement in recovering vascular cells (endothelial cells and PC/SMC) (Figure S2B). This is likely due to our robust two-step tissue homogenization, which effectively releases nuclei from vasculature that would otherwise be mechanically resistant.

The cellular composition altered substantially in response to HFD feeding, with variations depending on the fat depot (Figure 2J). HFD resulted in a reduction in adipocytes and an increase in macrophage/monocyte fractions in both eWAT and iWAT, with more prominent changes in eWAT (Figure 2J, S2C). APC fractions were smaller in eWAT compared to iWAT under chow diet, with marginal effects under HFD. The population size of endothelial cells did not show notable differences between the conditions (Figure S2C). Most of the minor cell populations exhibited minimal differences across the conditions, except for dendritic cells which increased in eWAT in response to HFD, and mast cells which were more abundant in iWAT than eWAT (Figure S2C).

Obesity induces the transition of APC from an early state to a committed state in specifically eWAT

We next performed subclustering of each cell type found in mouse adipose tissues to gain insight into the cellular states and dynamics of cell subtypes during obesity, focusing on their unique gene expression profiles and proportional changes during obesity in different depots. Analysis of APCs revealed four different subpopulations (Figure 3A, 3B), which generally aligned with prior studies. APC1 and APC2 corresponded to ‘adipocyte progenitors,’ characterized by high expression of both Dpp4 and Pi16²¹, while displaying distinct gene expression profiles from each other (Figure S3A). Notably, APC2 exhibited higher expression of fibrotic genes such as Fn1 and Loxl1 (Figure S3B), resembling ‘fibro-inflammatory progenitors²⁰’. APC3 was identified as ‘adipose regulatory cells,’ expressing F3 (encoding CD142) and showed strong expression of secretory signaling proteins (Figure 3A, S3C), supporting its paracrine regulatory role¹⁸. APC4, the largest subpopulation, represented ‘committed preadipocytes,’ expressing lipid metabolism-related genes such as Lpl, Hsd11b1, and Fabp4 (Figure 3A, S3D). The distribution of APC subpopulations differed between adipose depots, with eWAT displaying a higher fraction of APC3 and APC4, while iWAT showed a larger fraction of APC1 (Figure 3C). Specifically in eWAT, there was an expansion of APC4 alongside a reduction in APC1 and APC2 populations (Figure 3C) during obesity, suggesting distinct tissue expansion mechanisms in each adipose tissue depot. eWAT appears to promote de novo adipogenesis through conversion of adipocyte progenitors to committed preadipocytes, as demonstrated in a prior study²⁹.

Figure 3.

Adipose tissue immune cells acquire pro-inflammatory phenotypes predominantly in eWAT during obesity

We identified seven subpopulations of macrophages and monocytes based on distinct marker expression profiles (Figure 3D, 3E, S4A). The two largest populations, MC1 and MC2, expressed Clec10a (encoding CD301) and Mrc1 (encoding CD206), indicative of alternatively-activated M2 macrophages (Figure S4A). However, MC2 stood out for its enrichment of genes associated with antigen presentation and integrin binding, including major histocompatibility complex (MHC) markers such as H2-Eb1 and H2-Aa (Figure S4B, S4C). These MC1 and MC2 cells corresponded perivascular-like (PVM) and non-perivascular-like macrophages (NPVM), respectively, as described by Sávári et al.²² MC3 represented ‘lipid-associated macrophages (LAMs),’ characterized by high expression of lipid metabolism genes (Figure S4D). We found two minor macrophage subpopulations, MC6 and MC7, which resembled previously described regulatory macrophages (RegM) and collagen-expressing macrophages (CEM), respectively (Figure 3D, S4A). Their prevalence showed depot-dependent differences, with eWAT demonstrating a higher abundance of RegM and iWAT containing more CEM (Figure 3F). MC4 and MC5, accounting for approximately 20% of the myelocyte populations, were identified as classical monocytes (CMO) and non-classical monocytes (NCMO), respectively (Figure 3D, S4A). In response to HFD feeding, only eWAT exhibited a marked reduction in MC1 and MC2 populations with a substantial increase in MC3 (LAMs) (Figure 3F).

Lymphocytes were clustered into four T cell subpopulations – regulatory (Treg), CD4+, CD8+, and natural killer T cells (NKT) – along with two B cell subpopulations (Figure 3G, 3H, S4E). The distribution of T cell subpopulations differed between adipose depots: CD4+ and CD8+ T cells were more abundant in iWAT, whereas Treg were more prevalent in eWAT (Figure 3I). Conversely, B cells exhibited relatively high homogeneity in proportions between fat depots under chow conditions (Figure 3G, 3I). The most remarkable change in lymphocytes during HFD feeding was the emergence of a small B cell subpopulation, BC2, observed in eWAT (Figure 3G), which highly expressed genes involved in B cell proliferation, such as Prkar2b and Mcc (Figure S4F). Collectively, our high-quality data enabled accurate and detailed characterization of diverse immune cell populations in adipose tissue, revealing their dynamic changes during obesity in a depot-dependent manner.

Obesity-induced changes in vascular cell subpopulations represent depot-dependent dynamics of adipose tissue remodeling

The higher recovery of vascular cells in our snRNA-seq data provided comprehensive insights into the vasculature of adipose tissues. We identified nine subpopulations of vascular cells (Figure 3J, 3K), grouped into three major types: blood endothelial cells (blood EC), mural cells, and lymphatic endothelial cells (LEC) (Figure S5A). Blood EC included venous (VenEC), arterial (ArtEC), and capillary endothelial cells (CapEC). CapEC was further divided into arterial capillary endothelial cells (AcapEC) and venous capillary endothelial cells (VcapEC) (Figure S5B), with the latter being the largest vascular cell type in adipose tissues (Figure 3L). We additionally recognized angiogenic endothelial cells (AngEC), comprising migratory tip cells expressing extracellular matrix (ECM) genes, and proliferative stalk cells with elevated expression of cell cycle genes (Figure S5B). Furthermore, we discovered an endothelial cell type characterized by high expression of lipid metabolism genes (LipEC) (Figure S5B). Following HFD feeding, an expansion of AngEC and LipEC was observed in both fat depots (Figure 3L), suggesting increased vascularization and lipid transport during obesity-induced adipose tissue expansion. Intriguingly, we noted that the VcapEC population drastically decreased specifically within eWAT during HFD feeding (Figure 3L), indicating impaired metabolite transport into the venous circulation, which may potentially contribute to the pathology of visceral adipose tissue.

In mural cells, pericytes (PC) and smooth muscle cells (SMC) were identified (Figure S5C), with the PC population notably larger than that of SMC (Figure 3L), consistent with their respective association with capillaries and larger vessels^30,31. Lastly, we identified LEC (Figure S5D), comprising 1-3% of the vascular cell population, which interestingly diminished following HFD feeding in both eWAT and iWAT (Figure 3L). The HFD-induced reduction of LEC, along with the aforementioned decline in VcapEC, implies that metabolic drainage into lymphatic vessels and the venous circulation could be a key pathological feature of visceral adipose tissues during obesity.

Adipocyte subpopulations embody diverse functional and cellular states influenced by obesity

snRNA-seq has recently begun to reveal the heterogeneity of adipocytes within white adipose tissue^22–24. In our dataset, we identified six subpopulations of adipocytes, characterized by a gradient of unique gene expression signatures among them rather than distinct cell types (Figure 4A, 4B). Ad1-Ad3 represented a profile of functional mature adipocytes with high expression of insulin signaling genes and lipid metabolism genes, with Ad2 demonstrating particularly robust expression of lipid metabolism genes (Figure S6A). Ad3, along with Ad4, also showed elevated expression of genes associated with cell junction, extracellular matrix (ECM) organization and calcium homeostasis (Figure S6B). Ad6 was distinguished by its pronounced expression of genes involved in unfolded protein binding, oxidative stress, cytoskeleton, antigen presentation, and adipokines such as Lep, Retn and Rbp4 (Figure S6C). In response to HFD feeding, Ad1 and Ad2 were significantly reduced in both eWAT and iWAT, while Ad3 showed an elevation in iWAT. Notably, Ad6 exhibited a substantial increase in eWAT and a moderate rise in iWAT (Figure 4C). These findings indicate that the identified adipocyte subpopulations may represent diverse functional and cellular states influenced by nutritional conditions.

Figure 4.

Adipocytes undergo hypertrophy, an increase in cell size, during obesity, accompanied by significant cytoskeletal rearrangement³². Thus, we questioned whether the classification of adipocyte subpopulations aligns with the spectrum of adipocytes during HFD-induced expansion, particularly regarding their size. To address this question, we sought to compare gene expression profiles between large and small adipocytes. Interestingly, we observed enlarged nuclei in hypertrophic adipocytes isolated from HFD-fed mice (Figure 4D), which correlated positively with adipocyte cell size in both fat depots (Figure 4E). Leveraging this finding, we profiled the size distribution of adipocyte nuclei, labeled with mCherry and GFP, isolated from NuTRAP mice under different nutritional states using flow cytometry. As indicated by the intensity of mCherry and GFP, eWAT and iWAT from chow-fed mice primarily displayed smaller adipocyte nuclei. In HFD-fed mice, eWAT predominantly contained larger adipocyte nuclei, while iWAT showed a heterogeneous size range from small to large nuclei (Figure 4F, 4G). Next, we sorted and separately collected adipocyte nuclei by size from the same iWAT of HFD-fed mice for RNA-seq analysis (Figure S6D) to define the genuine gene expression signature associated with cell/nucleus size, minimizing potential confounding factors such as diet or depot dependency. We identified adipocyte hypertrophy signature genes based on their enrichment in large adipocyte nuclei compared to smaller adipocyte nuclei (Figure 4H). They also exhibited higher overall expression levels in adipocyte nuclei in iWAT from HFD-fed mice compared to those from chow-fed mice (Figure 4H). When overlaying the expression of these hypertrophy signature genes across adipocyte subpopulations, we observed high expression in Ad3-Ad4 and Ad5-Ad6 (Figure 4I), suggesting that they may represent hypertrophic adipocytes. However, certain genes differentiated between Ad3-Ad4 and Ad5-Ad6. For example, Hdac9 was more enriched in Ad3, while Pde2a was in more enriched in Ad5-Ad6 (Figure 4I). Importantly, there were marked distinctions in the expression of insulin signaling- and lipid metabolism-related genes between Ad3 and Ad6, with Ad6 having substantially lower gene expression levels (Figure S6A). This implies that Ad3 represents ‘metabolically healthy hypertrophic (MHH)’ adipocytes, while Ad6 represents ‘metabolically unhealthy hypertrophic (MUH)’ adipocytes.

To further discern the obesity-driven dynamics of adipocyte subpopulations, we performed pseudotime trajectory analysis, starting with the Ad1-Ad2 subpopulations that were prevalent under chow conditions (Figure 4J). Two divergent trajectories were revealed: One trajectory moves toward Ad3 and Ad4, which may depict an ‘adaptive healthy transition.’ Conversely, the alternate trajectory progresses towards Ad5 and ultimately to Ad6, indicating a potential ‘maladaptive dysfunctional transition.’ In addition, the farthest position of Ad6 along the pseudotime scale (Figure 4J) indicates prolonged pathological conditions. Importantly, eWAT and iWAT displayed contrasting preferences in adipocyte trajectories during HFD feeding: the dysfunctional trajectory leading to Ad6 was dominant in eWAT while the adaptive trajectory toward Ad3 remained more prevalent in iWAT (Figure 4C). This supports the depot-dependent dynamics of adipocytes that contribute to distinct metabolic outcomes during obesity.

Adipocyte subpopulations differentially contribute to tissue-level changes during obesity

We next aimed to evaluate the functional role of individual adipocyte subpopulations on tissue-level alterations during obesity. First, we identified genes that were differentially expressed in response to HFD within each subpopulation in eWAT or iWAT. Subsequently, we classified these genes into four clusters based on their pattern of HFD-induced changes across adipocyte subpopulations (Figure 5A, Table S1). Gene cluster 1 represented a group of genes generally upregulated by HFD across all subpopulations, with the most significant induction observed in Ad3, followed by Ad4 and Ad2 to a lesser extent. The genes associated with cell junctions, ECM and the cytoskeleton included in this cluster (Figure 5B) may describe biological responses during HFD-induced adipose tissue remodeling, which were featured in the aforementioned adaptive transition, among Ad2-Ad4 subpopulations. On the other hand, Gene clusters 2 and 3 comprised genes up- and down-regulated, respectively, with changes most pronounced in Ad6 and to a lesser extent in Ad5, during HFD feeding. The induction of genes in cluster 2, as well as the repression of genes in cluster 3, reflect stress response, inflammation, and cell death in adipose tissues, along with impaired machinery regulating gene expression (e.g., chromatin remodeling and mRNA metabolic process) during HFD-induced adipose pathology (Figure 5B). Genes in cluster 4, commonly down-regulated by HFD across all adipocyte subpopulations, had enriched metabolic processes such as lipid metabolism, thermogenesis, and insulin phosphatidylinositol 3-kinase (PI3K) signaling (Figure 5B), indicating a general decline in metabolic activity and insulin sensitivity across the entire adipocyte population rather than specific subpopulations. Collectively, these results suggest that tissue-level biological changes in adipose tissues during obesity result from a combination of alterations within specific subpopulations and throughout the entire adipocyte population.

Figure 5.

Ad6 may represent stressed and dying adipocytes with transcriptional shutdown

Our snRNA-seq data unveiled a previously unrecognized adipocyte subpopulation, Ad6, predominantly present in eWAT during HFD feeding. This subpopulation represented a dysfunctional state of adipocytes, characterized by elevated expression of stress response and inflammation genes, coupled with the loss of metabolic genes. Using the identified profile as a basis, we sought to confirm the presence of the Ad6 adipocyte subpopulation at the tissue level through adipocyte or whole-mount immunostaining. We first visualized reactive oxygen species (ROS) within adipocytes isolated from adipose tissues of chow- or HFD-fed mice using the chemical ROS reporter H2DCFDA. Under HFD conditions, a significantly larger fraction of adipocytes exhibiting ROS was observed, with effects more prominent in eWAT than in iWAT, together with significantly higher intensity (Figure 6A, 6B). This provides evidence of Ad6 adipocytes experiencing oxidative stress. Next, we inspected Ad6 adipocytes by employing whole-mount immunostaining of adipose tissues for beta-2 microglobulin (B2M), a marker highly enriched in this subpopulation (Figure S6C). B2M was primarily detected in eWAT under HFD feeding, particularly within crown-like structures (CLS), as well as in the cytoplasm of adipocytes neighboring these CLS (Figure 6C, 6D). As dead adipocytes are detected within CLS³³, this result suggests Ad6 as stressed and dying cells due to inflammation, both within and outside the CLS. Another hallmark of Ad6 adipocytes was the depletion of genes crucial for lipid metabolism, including PPARγ (Figure S6A). To mitigate the potential for artifacts derived from ambient RNA, we directly visualized PPARγ proteins within adipocytes in adipose tissues through whole-mount immunostaining utilizing NuTRAP mice. While most GFP-labeled adipocyte nuclei robustly expressed PPARγ proteins under typical conditions, a significant proportion of adipocyte nuclei from eWAT in HFD-fed mice lacked PPARγ proteins (Figure 6E, 6F), confirming the presence of PPARγ-deficient Ad6 adipocytes in eWAT during obesity. Furthermore, the substantially fewer transcripts, represented by UMIs or genes, detected in the Ad6 subpopulation (Figure 6G), suggests a global loss of gene expression. Taken together, Ad6 represents a subpopulation of stressed and dying adipocytes with transcriptional shutdown that arise during obesity.

Figure 6.

Discussion

snRNA-seq has emerged as an alternative approach to overcome the limitations of scRNA-seq. Nevertheless, technical challenges in obtaining high-quality nuclei and RNA have remained as significant hurdles in the standard protocol. In this study, we present a robust protocol for nucleus isolation that effectively addresses these issues across various tissue types. To demonstrate the efficacy of our method, we applied it to characterize one of the most challenging tissue types to investigate, adipose tissue. Our snRNA-seq data demonstrates its superior quality compared to previously published datasets, as well as provides novel insights into the underlying mechanisms of mouse adipose tissue remodeling during obesity.

We observed that significant variations in RNase abundance among various tissue types affect the quality of RNA obtained from isolated nuclei using standard protocols. Notably, brain tissue, which has been primarily used in prior snRNA-seq studies, exhibited the lowest RNase activity. This implies that the standard nucleus isolation protocols employed in those studies may not readily translate to other tissues, particularly those with high RNase activities. Recombinant RNase inhibitor proteins are commonly used to control RNA quality during nucleus isolation. However, these inhibitors marginally improved RNA quality, likely due to their specific targeting of certain RNase types, such as RNase A, 1, and T1. Instead, through experiments, we discovered that VRC, a potent inhibitor capable of targeting a wide range of RNases, significantly enhanced RNA quality in isolated nuclei. VRC functions as a transition state analog, binding to the active sites of RNases and inhibiting their enzymatic activity^34,35. Despite its demonstrated versatility in RNA experiments^34–36, the application of VRC in gene expression studies has been limited due to its inhibitory effects on other RNA enzymes, including reverse transcriptase³⁵. Our method addressed this issue by reducing VRC from nucleus isolation buffer prior to loading onto microfluidic devices, as confirmed by the robust qRT-PCR and snRNA-seq results we obtained. Recently, we encountered a few protocols that have been introduced to enhance RNA quality by employing VRC or acidic buffers during nucleus isolation for snRNA-seq, particularly for kidney⁴² and pancreas³⁷ samples. However, their applicability to other tissue types remains untested. In contrast, our optimized protocol successfully preserved not only nuclear RNA quality in various tissue types, but also the integrity of the nuclear structure. We noted that VRC has preservative effects on isolated nuclei, resulting in smoother and more intact membranes in larger nuclei when compared to standard protocols not using VRC. This observation aligns with a recent study that demonstrated the preservative properties of VRC, protecting tissue morphology, as well as RNA, protein, and genomic DNA quality²⁶. More importantly, we demonstrated the capability of our protocol to retain nuclear RNA quality for up to 24 hours post-tissue processing, even after undergoing additional steps such as FANS. This advancement holds the potential to greatly improve experimental logistics. Considering that snRNA-seq experiments typically require completion of nucleus isolation and loading onto microfluidic devices (e.g., 10X Genomics Chromium controller) within a narrow time window on the same day, our protocol offers greater flexibility for these critical steps, allowing them to be performed on different days if necessary. This adaptability has the potential to unlock a wide array of applications for snRNA-seq techniques. In the clinical setting, for example, this will enable patient biopsies to be collected for diagnostic tests and then be analyzed via snRNA-seq at a different facility.

Our single nucleus atlas of mouse white adipose tissues generally aligns with previous studies while offering deeper insights into depot specificity and cellular composition dynamics during obesity. An important advancement in our dataset was the improved recovery of vascular cell populations, which were underrepresented in prior studies, likely due to our protocol’s effective tissue homogenization and the high quality of RNA preserved in isolated nuclei. This enabled precise characterization of diverse endothelial and related vascular cell types. We observed a marked decrease in capillary endothelial cells in eWAT specifically during HFD feeding, consistent with the previously documented capillary rarefaction and elevated hypoxia associated with adipose tissue expansion^13,38. We further noticed that the VcapEC population, which is responsible for transporting metabolites into the venous circulation, was primarily affected. Additionally, we observed a HFD-induced decline in the LEC population, which plays a role in removing fluid and macromolecules from the interstitial space. These findings suggest that defective drainage of metabolites and molecules into the venous circulation may be another key feature of dysfunctional adipose tissue in obesity, warranting further investigation.

Another noteworthy discovery from our snRNA-seq data was gaining deep insights into the spectrum of adipocyte subpopulations and their dynamics during obesity. Based on our observation of a strong correlation between adipocyte cell size and nucleus size, we developed a novel approach to profile the size distribution of adipocyte subpopulations using nuclei as surrogates, which yielded intriguing findings. Adipocyte cell size appeared small and relatively uniform under chow conditions in both eWAT and iWAT. During HFD feeding, adipocytes notably enlarged, particularly in eWAT, whereas iWAT exhibited a broader distribution with a mixture of small and large adipocytes. Therefore, by isolating adipocyte nuclei of various sizes from iWAT, we identified gene signatures that determine adipocyte hypertrophy. Overlaying this data with snRNA-seq data revealed that cell size is a key determinant of adipocyte clustering in mice. This finding contrasts with the heterogeneity of human adipocytes, where insulin sensitivity is shown to be a primary factor in adipocyte classification rather than cell size, as demonstrated by spatial transcriptomics mapping³⁹. Interestingly, ex vivo differentiated human adipocytes derived from various donors displayed distinct cellular phenotypes, consistent with the heterogeneity in adipocyte subpopulations identified by snRNA-seq²³. This suggests that the significant genetic variation among humans, potentially serving as a major driver of adipocyte clustering, overrides the size-dependent differences found in inbred mouse strains with homogeneous genetic backgrounds. Nevertheless, our findings integrate various characteristics of adipocytes reported from diverse biological contexts, broadening our understanding of distinct adipocyte subpopulations. For instance, despite both Ad1 and Ad2 subpopulations having small cell size and abundance under chow conditions, they exhibited distinct expression levels of lipid metabolism genes. This finding may be consistent with previous studies demonstrating mosaic lipid transport activity among small white adipocytes from Rhesus Macaques monkeys⁴⁰. Also, the unique enrichment of genes involved in calcium transport and ossification in hypertrophic Ad4 adipocytes may account for the frequent occurrence of calcium deposits observed in obese adipocytes⁴¹.

We further noted two distinct types of hypertrophic adipocytes developed by HFD feeding: MHH adipocytes (Ad3) abundant in iWAT, and MUH adipocytes (Ad6) prevalent in eWAT. We not only validated the presence of MUH adipocytes, primarily around CLS in obese eWAT with oxidative stress and inflammation markers, but discovered their potential association with prolonged obesity through pseudotime analysis. Interestingly, MUH adipocytes showed deficient expression of the master transcriptional regulator PPARγ as well as a substantial decline in metabolic gene expression. The significant decrease in transcript counts in Ad6 suggests that these cells may represent stressed or dying cells, akin to the ‘stressed lipid-scavenging’ adipocytes described by Sávári et al²². This also aligns with our earlier findings, which demonstrated the failure of adipocytes to maintain their identity during obesity³². Note that these cells may pose challenges for reliable detection due to severely compromised RNA contents. This highlights the effectiveness and precision of our robust snRNA-seq methodology, which could be valuable for profiling other disease states with impaired cell viability and gene expression.

Lastly, we elucidated the contribution of distinct adipocyte subpopulations to various established signatures of obese adipose tissue. For instance, Ad2-Ad4 adipocytes appeared to undergo tissue remodeling for enhanced adipose functionality, such as storage capacity. In contrast, the traditional features of dysfunctional adipose tissue, including inflammation, cell death, and stress response, were predominantly derived from Ad5-Ad6 adipocytes. They were also responsible for chromatin alterations and defective gene expression. On the other hand, the core functions of adipocytes, including lipid metabolism and insulin signaling, were generally impaired in all adipocyte subpopulations. The distinctions between obesity-induced features occurring in the entire adipocyte population versus specific adipocyte subpopulations may explain the outstanding effectiveness of PPARγ agonists, such as thiazolidinediones (TZDs), in improving adipose function and glucose homeostasis⁴², surpassing the effects of other anti-inflammatory or anti-stress reagents^43–45. Our discovery emphasizes the need to tailor therapeutic approaches based on the specific biological processes to be targeted to improve adipose tissue function during obesity.

In conclusion, this study introduces a robust nucleus isolation method for snRNA-seq, leading to enhanced data quality. This advancement has enabled the discovery of novel insights into the dynamics of adipose tissue remodeling during obesity. In addition, this method holds great promise for widespread applications across various research fields and clinical diagnoses, involving numerous types of biological samples.

Materials and methods

Animal studies

All animal experiments were performed according to a protocol approved by the Indiana University School of Medicine (IUSM) Institutional Animal Care and Use Committee (IACUC). The mice were housed under a 12h light-dark cycle at 22 °C with free access to food and water. To label adipocyte nuclei, NuTRAP mice (Jackson Laboratory, 029899) were crossed with Adipoq-Cre mice (Jackson Laboratory, 010803). In experiments assessing nuclear RNA quality, male C57BL/6J mice at 8-10 weeks of age were used for tissue collection. To generate dietinduced obesity models, male mice were subjected to either a HFD (60% kcal from fat, Research Diets D12492i) or a standard chow diet (control) from eight weeks of age for a duration of 10 weeks. Upon dissection, tissues were immediately snap-frozen in liquid nitrogen. During the dissection of eWAT, we carefully removed epididymal ducts from the adipose tissues, minimizing the risk of contamination. For iWAT, the lymph nodes were removed to avoid over-representation of lymphocytes.

Nucleus isolation

For each experiment, the nucleus preparation buffer (NPB) was freshly prepared at room temperature with the following components: 10 mM HEPES (pH 7.5), 1.5 mM MgCl2, 10 mM KCl, 250 mM sucrose, 0.1% NP-40, and 1 mM DTT. Next, vanadyl ribonucleoside complex (VRC; New England Biolabs) was added to the NPB at a concentration of 10 mM. The mixture was vigorously vortexed until the VRC was completely dissolved in the NPB. The NPB with VRC was then chilled on ice, and 0.4 U/µl of NxGen RNase inhibitor (Biosearch Technologies) and 1X Halt protease inhibitor (Thermo Fisher) were added. Once the NPB was ready, frozen tissue samples were pulverized using mortars and pestles in liquid nitrogen. The resulting tissue powders were immediately transferred to and homogenized in chilled Douncers placed on ice, using 5-7 ml of the NPB. The tissue homogenates were filtered through 100 µm cell strainers and centrifuged at 200 g for 10 minutes at 4°C. After removing the supernatant, the nucleus pellet was resuspended in 5 ml of PBS-N-VRC, which is calcium and magnesium-free phosphate-buffered saline (PBS) with 0.1% NP-40 and 5 mM VRC. Hoechst 33342 (Thermo Fisher) was added to the tissue homogenates at a concentration of 1 µg/ml to label nuclei. The mixture was then centrifuged at 200 g for 5 minutes at 4°C. Subsequently, the nucleus pellets were resuspended in PBS-N-VRC buffer and filtered through 40 µm cell strainers. The protocol for nucleus isolation is currently under the provisional patent application process.

Nuclear RNA extraction for quality analysis

The nucleus resuspension in PBS-N-VRC was centrifuged at 200 g for 10 minutes at 4°C, and the supernatant was discarded. To reduce VRC amount, the nucleus pellet was resuspended in PBS-N (PBS with 0.1% NP-40), and then centrifuged again at 200 g for 10 minutes. After removing the supernatant, the nucleus pellet was subjected to RNA extraction using TRIzol as per the manufacturer’s instructions. Briefly, the nucleus pellet was thoroughly resuspended in TRIzol (Invitrogen) by vigorous vortexing and subsequently mixed with chloroform. Following centrifugation at 13,000 g for 10 minutes at 4°C, the aqueous phase was collected and combined with isopropanol and 1.5 µg of GlycoBlue coprecipitant (Invitrogen). After a 10-minute incubation at room temperature, the mixture was centrifuged at 13,000 g for 10 minutes at 4°C. The RNA pellet was washed with 75% ethanol and then dried. Subsequently, the pellet was resuspended in nuclease-free water, and the RNA concentration was measured using a Qubit Fluorometer with the Qubit RNA High Sensitivity Assay Kit (Invitrogen). The quality of the RNA was assessed using a high sensitivity electrophoresis system, the Agilent Bioanalyzer. Since the automatically generated RNA integrity number (RIN) provided by the system did not accurately reflect the quality of nuclear RNA due to the presence of a substantial amount of small RNA molecules, we opted for a qualitative analysis of the RNA quality. For time-course experiments to assess RNA quality over time, the nuclear suspension in PBS-N-VRC was kept at 4 °C for the indicated duration prior to RNA extraction.

For visual inspection of nucleus quality, the nucleus resuspension was mounted on hemacytometer and visualized using a Zeiss Observer.Z1 fluorescence microscope using bright field and DAPI filters.

snRNA-seq with nucleus sorting

The nucleus resuspension in PBS-N-VRC was subjected to nucleus sorting using a BD FACSAria Fusion system. We set gating based on forward scatter (FSC), side scatter (SSC), and Hoechst fluorescence to eliminate debris and nucleus doublets/multiplets. Sorted nuclei were collected in PBS-N-VRC, then centrifuged at 200 g for 10 minutes at 4°C using a swing bucket rotator. After removing the supernatant, the nucleus pellet was resuspended in PBS-N, filtered through 20 µm cell strainers, and centrifuged at 200 g for 5 minutes at 4°C in a swing bucket rotator. Subsequent to removing the supernatant, the nucleus pellet was resuspended in PBS-N in volumes aimed at achieving a nucleus density of approximately 1000 nuclei/µl. The nuclei in the resuspension were counted using a hematocytometer or automated cell counters. Following the counting, 16,000 nuclei (with a target of 10,000) were loaded onto the 10x Genomics Chromium controller. Libraries were prepared using the Chromium Single Cell 3’ Reagent Kit v3.3 (10x Genomics), indexed using sample barcodes, and subsequently sequenced on an Illumina NovaSeq 6000 sequencer with paired-end dual indexing.

Bulk nuclear RNA-seq with nucleus sorting

For bulk nuclear RNA-seq, nuclei were sorted based on the aforementioned gating criteria following isolation. Sorted nuclei were washed in PBS-N as for snRNA-seq but skipping additional filtration. Nuclear RNA was extraction using TRIzol as stated above. Additionally, the extracted RNA was purified by on-column DNA digestion using TURBO DNase kit (Invitrogen). For library construction, 20-50 ng of purified RNA was processed by NEBNext rRNA Depletion Kit (New England BioLabs, Inc.) for removing ribosomal RNA, Next the RNA was converted to cDNA using Maxima Reverse Transcriptase (Thermo Fisher Scientific) and NEBNext mRNA Second Strand Synthesis kit. Following purification based on size selection using AMPure XP beads (Beckman Coulter), sequencing libraries were generated through tagmentation using Nextera XT DNA Library Preparation kit and subsequent PCR amplification. The quantity and quality of the resulting libraries were assessed using Qubit and Agilent Bioanalyzer, respectively. Finally, the libraries were sequenced on an Illumina NextSeq500.

H2DCFDA staining with isolated adipocytes

Fresh adipose tissues were minced and digested in PBS containing collagenase D (1.5 U/ml) and dispase II (2.4 U/ml) with shaking for 30-40 min at 37°C. After digestion, the tissue homogenates were filtered through 250 µm cell strainers and left to settle for 10 min. Floating adipocytes were collected, washed in PBS containing 1% fat-free bovine serum albumin (BSA), and incubated with H2DCFDA (Invitrogen, 20 uM) and Hoechst 33342 (Thermo Fisher) for 30 min in the dark at room temperature. The cells were washed in PBS containing 0.05% Tween20 and 0.1% fat-free BSA twice, placed on a chamber made on slides using 90% glycerol, and covered with coverslips. The slides were visualized using a Leica SP8 confocal microscope.

Whole-mount immunostaining

Fresh adipose tissues were fixed in 10% buffered formalin phosphate for 1 day and then washed in PBS. Subsequently, the fixed tissues were permeabilized in PBS containing 1% Triton X-100 for 1 hour, following by blocking in PBS containing 5% BSA and 0.05% Tween20, for another hour. The tissues were incubated with (1) anti-B2M (Proteintech) or (2) anti-GFP (Novus Biologicals) and anti-PAPRg (Cell Signaling Technology) at 1/200 dilution for 1 day. After washes in PBS/0.05% Tween20, the tissues were further incubated with Alexa Fluor conjugated antibodies at 1/100 dilution for 1 hour. For B2M staining, the tissues were additionally incubated with BODIPY 493/503 (Invitrogen) for 30 min. After washed twice in PBS/0.05% Tween20 including once with Hoechst, the tissues were mounted on 35 mm glass bottom dish using glycerol mounting buffer and covered with coverslips. The mounted tissues were visualized using a Leica SP8 confocal microscope.

Quantitative real time PCR

Total RNA was extracted from sorted nuclei using TRIzol, and the extracted RNA was then converted into cDNA using the High-Capacity cDNA Reverse Transcription Kit from Applied Biosystems. qRT-PCR was conducted using the SYBR Green PCR Master Mix on a QuantStudio5 system. The fold change was determined using the ΔΔCT method by comparing target gene expression with a reference gene 36B4 (Rplp0). The primers used for qRT-PCR are provided in Table S2.

Bioinformatics analysis

Sequencing, alignment, and quality control

Data was processed using Cell Ranger (10x Genomics, v7.1.0) to demultiplex raw sequencing reads to FASTQ format files, align them to the mm10 reference genome, and perform filtering, barcode counting, and UMI counting. CellBender v0.3.0²⁸ was applied to remove counts from ambient/background RNA. The CellBender output files (.h5 format) were read into R (version 4.3.1) using Seurat 4.3.0, and each sample was further filtered prior to downstream analyses based on the number of UMIs (ζ 500) and genes (ζ 250), mitochondrial ratio (< 15%), and the number of genes detected per UMI (> 0.8). At the gene level, genes detected less than 10 nuclei and mitochondrial genes were filtered. To identify and remove doublets, we employed DoubletFinder v2.0⁴⁶ using an estimated multiplet rate derived from the total number of recovered nuclei.

Integration, clustering, annotation, and subclustering

Gene counts were first normalized using SCTransform, and regressed on mitochondrial read ratio and cell cycle score by Seurat. For data integration among samples, we performed canonical correlation analysis (CCA) to identify shared sources of variation among samples, focusing on 3,000 genes. Next, we employed the FindIntegrationAnchors function to determine mutual nearest neighbors (MNNs), commonly referred to as anchors. Then, utilizing the identified anchors, the datasets were integrated using the IntegrateData function. For graphbased clustering, we visualized the integrated data using dimensionality reduction techniques, such as principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP), and clustered the nuclei with a resolution of 0.6. As a result, clusters were annotated as adipocytes, APC, mesothelial cells, endothelial cells, PC/SMC, macrophages, dendritic cells, mast cells, T cells, and B cells using their marker genes. For further investigation, the integrated data was subdivided into separate objects by cell type and reclustered to achieve higher resolutions, with minimal additional filtering applied to potential doublets. To understand distinct functionality of subclusters, if needed pathway analysis was performed on the marker genes using clusterProfiler v4.8.2⁴⁷. For visualization of gene expression levels within the specific cells, we used the R package called scCustomize⁴⁸ which generates plots where cells expressing the gene of interest can be superimposed on top of the negative cells with the argument of ‘order=TRU’.

Pseudobulk differential expression analysis and downstream analyses

To identify genes differentially expressed between chow and HFD conditions within particular cell types, we aggregated the counts and metadata to the sample level and compared the counts between conditions using edgeR v.3.42.4⁴⁹. Significance was determined based on false discovery rate (FDR) <0.05 and absolute log2 fold change >1. Differentially expressed genes in either eWAT or iWAT were clustered using K-means clustering through Morpheus (https://software.broadinstitute.org/morpheus), and each group of genes were passed to Metascape⁵⁰ for pathway analysis.

Single cell trajectory analysis

We used Monocle 3⁵¹ to infer trajectories within cell types, potentially representing the sequences of gene expression changes that each cell/nucleus undergoes as part of a dynamic biological process. This analysis was conducted using the Seurat object containing subclustering information.

Analysis of published data

For comparison of data quality, scRNA-seq or snRNA-seq datasets from previous studies were obtained from Emont et al.²³ (GSE176171) and Sávári et al.²² (GSE160729). The data were processed using Cell Ranger and read into R using Seurat. The median number of UMIs and genes detected per cell, the fraction of reads in cell, and the median number of mitochondrial read ratios were evaluated.

Statistical analysis

Statistical details are available in the figure legends. Unless otherwise noted, all bar graphs and scatterplots show both mean ± SEM and individual values. P values were calculated using GraphPad Prism v10.0.0. and represent the results of Student’s t test, one-way ANOVA with multiple comparisons, or linear regression test.

Data availability

The raw and processed data of the snRNA-seq (GSE241987) and bulk nuclear RNA-seq (GSE261417) in this paper are available in the Gene Expression Omnibus (GEO) repository.

Acknowledgements

We thank Dr. Hongyu Gao for providing guidance in setting up the computational analysis pipeline, and Dr. Charlie Dong for sharing equipment. We are grateful to the Flow Cytometry Core and the Center for Medical Genomics at Indiana University School of Medicine. This study was supported by National Institute of Diabetes and Digestive and Kidney Diseases (R01DK129289), and American Diabetes Association Junior Faculty Award (7-21-JDF-056) to H.C.R. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute.

Author contributions

H.C.R. conceived, designed, and supervised the project. H.C.R., J.S., O.S., J.W., K.K., and A.P. performed experiments and interpreted the data. J.S. performed computational data analysis with guidance from D.J.L, L.C.D, and J.K under supervision of H.C.R. J.S. and H.C.R. wrote the manuscript with inputs from all the other authors.

Declaration of interests

H.C.R. is the inventor on a provisional patent application that covers the nucleus isolation protocol for snRNA-seq. The remaining authors declare no competing interests.