Single-cell transcriptional networks in differentiating preadipocytes suggest drivers associated with tissue heterogeneity

Abstract

White adipose tissue plays an important role in physiological homeostasis and metabolic disease. Different fat depots have distinct metabolic and inflammatory profiles and are differentially associated with disease risk. It is unclear whether these differences are intrinsic to the pre-differentiated stage. Using single-cell RNA sequencing, a unique network methodology and a data integration technique, we predict metabolic phenotypes in differentiating cells. Single-cell RNA-seq profiles of human preadipocytes during adipogenesis in vitro identifies at least two distinct classes of subcutaneous white adipocytes. These differences in gene expression are separate from the process of browning and beiging. Using a systems biology approach, we identify a new network of zinc-finger proteins that are expressed in one class of preadipocytes and is potentially involved in regulating adipogenesis. Our findings gain a deeper understanding of both the heterogeneity of white adipocytes and their link to normal metabolism and disease.

Fig. 1. PAGODA and t-SNE reveals at least two clusters in differentiating preadipocytes.

a Single-cell RNA seq was performed on differentiating preadipocytes beginning at 80% confluency and differentiated for 7 days. PAGODA was used to determine the optimal cell clustering based on the genes driving the heterogeneity. The result was plotted by t-SNE. A total of 2092 are shown. b Single-cell expression profiles in were analyzed to determine the genes and pathways driving the heterogeneity. PAGODA was used to perform weighted principle component analysis on pre-defined (e.g., Gene Ontology terms) and de novo gene sets. The gene sets were scored on their significance. Correlated gene sets were coalesced in order to reduce redundancy. The heatmap of the significance gene sets shows a few de novo gene sets captured major aspects (i.e., non-redundant principal components) of heterogeneity. c Differential gene expression was performed between the left and right clusters for each day of differentiation in differentiating preadipocytes beginning at 100% confluency. Gene set enrichment analysis was performed on the differentially expressed genes and the top 10 up- and down-regulated pathways sorted by z-score are shown for day 7. d Differential gene expression was performed between the left and right clusters of day 7 adipocytes. The genes shown are previously-studied adipocyte markers for lineage or differentiation stage. For each gene, the maximum likelihood estimate and 95% confidence interval of the log2 expression ratio (right cluster over left cluster) is shown. The brown adipocyte markers MYF5 and UCP1 were not detected in any cells.

Fig. 2. Clonal preadipocyte cell lines cluster into distinct metabolic signatures.

RNA seq was performed on 35 (9 clones each from 4 subjects; 1 library failed QC) clonally expanded immortalized human neck white subcutaneous preadipocytes and analyzed with PAGODA. a Principal Component Analysis of subject-corrected gene expression profiles. b Top gene sets capturing significant aspects of transcriptional heterogeneity among the 36 clones. Each of these aspects had several genes marking potential subpopulations of clones within the clusters. Basal oxygen consumption (OCR) and basal extracellular acidification rate (ECAR) were measured in preadipocytes via the Seahorse XF Analyzer. Glucose uptake was measured with radiolabeled 2-deoxy-glucose in preadipocytes. The levels of PPARG were measured after 18–21 days of differentiation. c Venn diagram of the differentially expressed genes between clusters A and B in the clonal cell line data versus the clusters in day 0 cell in the single-cell data (see Fig. 1A). d The seven most extreme clones of cluster A compared against the seven most extreme clones of cluster B show differences in glucose uptake (FWER ~ 0.059) and differentiation capacity (FWER ~ 0.095). Asterisks indicate FWER < 0.1 as assessed by a one-way ANOVA followed by a Bonferroni correction (N = 7 clonal cells lines). Bars indicate mean ± s.e.m. e The four phenotypes were correlated (Spearman) to every gene in the clones to generate phenotype correlation vectors. The top ten genes for each phenotype are shown in a heatmap. f The four phenotype correlation vectors (oxygen consumption, extracellular acidification, glucose uptake, and PPARG expression after differentiation) were then correlated (Spearman) to the single-cell gene expression profiles for each cluster. The distribution of correlation coefficients was compared between the left and right cluster for each day. Boxplots are centered on the median, the interquartile range (IQR) spans the 25–75% percentile, and the whiskers extend to 1.5 times the IQR above the 75% percentile (maximum) and below the 25% percentile (minimum). Points indicate correlation coefficients of individual cells. Asterisks indicate FWER < 0.05 as assessed by an ANOVA followed by a Bonferroni correction.

Fig. 3. Subnetwork detection reveals a cluster of ZNFs in adipogenesis.

a An overview of the subnetwork detection algorithm where the expression of each gene can be described as the sum of underlying biological networks. b The subnetwork detection algorithm was used to reveal a set of cells primarily found in the left cluster. The detected subnetwork had 30 connected genes. A total of 4319 cells are shown. c Pathway enrichment revealed these 80 genes were enriched in KRAB domains and C2H2 zinc-fingers. d The mean of the max expression of the ZNFs (ZNF264, ZNF490, ZNF587, and ZNF714) in the detected cells showed a negative correlation to the adipogenic marker FABP4. Line indicates mean expression, gray shaded area indicates 95% confidence interval. e RNA-seq of whole preadipocyte cultures during differentiation show that the ZNFs do not decrease over time. Bars indicate mean ± s.e.m. N = 2 independent experiments.