Systems view of adipogenesis via novel omics-driven and tissue-specific activity scoring of network functional modules

The investigation of the complex processes involved in cellular differentiation must be based on unbiased, high throughput data processing methods to identify relevant biological pathways. A number of bioinformatics tools are available that can generate lists of pathways ranked by statistical significance (i.e. by p-value), while ideally it would be desirable to functionally score the pathways relative to each other or to other interacting parts of the system or process. We describe a new computational method (Network Activity Score Finder - NASFinder) to identify tissue-specific, omics-determined sub-networks and the connections with their upstream regulator receptors to obtain a systems view of the differentiation of human adipocytes. Adipogenesis of human SBGS pre-adipocyte cells in vitro was monitored with a transcriptomic data set comprising six time points (0, 6, 48, 96, 192, 384 hours). To elucidate the mechanisms of adipogenesis, NASFinder was used to perform time-point analysis by comparing each time point against the control (0 h) and time-lapse analysis by comparing each time point with the previous one. NASFinder identified the coordinated activity of seemingly unrelated processes between each comparison, providing the first systems view of adipogenesis in culture. NASFinder has been implemented into a web-based, freely available resource associated with novel, easy to read visualization of omics data sets and network modules.


SUPPLEMENTARY FILE 1
Although the goal or our analysis was to provide a systems view of adipogenesis, much of the previous research in this field focused on either individual genes or subsets of genes in a pathway. To compare results obtained in our study with the published literature, we also provide below a more detailed analysis of data by pathway, gene, or sets of genes within a pathway. The lists of all pathways at a time point (TP)

Oil Red O staining method
SGBS pre-adipocytes were cultured and induced to differentiate into adipocytes as described above.
For Oil Red-O staining, at specific time points, the cells were gently washed twice with PBS and were fixed by exposing the cells to 10% formalin for 45 minutes, following which the cells were washed three times with deionized water. Fixed cells were stained with a commercially available

NASFinder identified known pathways during adipogenesis: the Leptin Pathway receptor
The leptin transcriptional profile was analyzed across the differentiation process: a major function of mature adipocytes is the production of leptin 1, 2 . NASFinder identified changes in expression of genes in the leptin pathway at 48, 96, 192 and 384 hours (Supplement Fig. 2) after induction.
Leptin receptor (LEPR) was the source node for this pathway and LEPR itself was differentially expressed over controls at 96 (1.4x), 192 (1.5x), and 384 (1.02x) hours but not at 0 or 6 hours. The common differentially expressed genes at 48, 96, 192, or 384 were PRKAG1, PRKAG2, ACACA, LEPR, PRKAA1, and CPT1A (Supplement Table 1). The AMP kinase regulatory gene PRKAG2 is up-regulated at each time point compared to progenitor cells reflecting the increased energy demands of the differentiation process and the role of this subunit in glucose metabolism 3,4 .
ACACA, the rate-limiting step in long chain fatty acid synthesis found by NASFinder in the network, is up-regulated during adipogenesis compared to its expression pre-differentiation (Supplement Table 3). Down-regulation of CPT1A is consistent with a change in energy metabolism from fatty acid oxidation of the progenitor cell to glucose. Each illustrated pathway is obtained as the intersection between two networks, the first being the subnetwork of the adipocyte interactome identified by NASFinder as active on the basis of the experimental evidence, and the second being canonical leptin pathway. Contextual experimental information is provided by the 1-neighborhood nodes (see Fig. 2 of main text for more details). The leptin receptor (LEPR) here is represented as a blue triangle.

Time Point Analysis
The specific transcriptomic changes occurring during the adipocyte differentiation process were  Table 2 Fold Change of Genes Described for 6, 48, and 96 hr   Table 2) at this time point and specifically, the key regulatory genes IL1A, IL1B, IL6, and IL12A. Down-regulation of these genes would decrease synthesis of extracellular cytokines that could block adipogenesis 8 Table 2). Experimentally manipulated changes in the expression of the functionally related Arp2/3 have been shown to alter adipogenesis of 3T3-L1 cells due to disruption of normal cytoskeletal remodeling 16 .
The circadian rhythm pathway (NAS = 0.3) was strongly induced at 48 hours while 50% of genes involved in the WNT signaling pathway were down-regulated 20 . ARNTL (+1.7 foldarylhydrocarbon nuclear translocator protein, involved circadian rhythmicity) has been shown to suppress WNT signaling 21 . Other

hours versus Control
Based on work from our laboratory and others, 96 hours is about the half-way point to a fully differentiated adipocyte 22 Table 3). However, the 96-hour "snap-shot" of 31 (of a total of 59 genes) up-regulated in the Focal Adhesion network identified subsets of gene families that were oppositely regulated during this transition. For example, different members of the integrin, CD (cluster of differentiation), and laminin gene set were either up-regulated or down-regulated (Supplement Table 3) which may suggest that different combinations and interactions between these proteins alter complexes necessary for maintaining the shape of pre-adipocytes or the rounded morphology required for the adipogenesis process. Changes in cytoskeleton structure and its interactions with integral membrane proteins can be disrupted by mechanical stress or changes in the expression of key regulator genes such as calreticulin (rev in 24 ). Supplement Table 3.
Focal Adhesion genes at the 96h timepoint (0 neighborhood).  Table 4). Major components of this pathway include several anaphase promoting complex/cyclosome (APC/C) proteins as well as members of the proteasome (PSM) protein family (see Supplement Table 4). Interpreting how individual proteins contribute to the structure and regulatory properties of APC/C is challenging because of the structural complexity and interactions with components of diverse cellular processes 27 . Comparing the transcriptional regulation at 192 and 384 hours post induction (Supplement Table 4 Table 5). Transcription of PI3KR1 is done by PPARγ 18 and, in turn, this regulatory subunit helps orchestrate not only adipogenesis but also the mature adipocyte state through many related and unrelated pathways (see TP crosstab file).  The source node receptor for integrin 4 pathway (100% up-regulated), ITGA4 (integrin alpha 4), has not previously been described in adipocytes but is known to be associated with the ITGA6 in cells such as epithelial cells, Schwann cells etc 31 . ITGA6 was demonstrated to support the clustering of growth arrested preadipocytes on the basement membrane, a critical step in establishing contact inhibition and preventing the preadipocytes from re-entering the cell cycle 32 . The integrin 3 pathway was mediated via the TGFBR2 and TGFBR3 receptors which were 66.7 and 100 % down regulated, respectively, at 48 hours compared to the 6 hour time point. These pathways were marginally up-regulated in mouse 3T3L1 adipocytes 33 . While many of the pathways identified by NASFinder were independent of each other at the 0-neighborhood level, a number of sub-networks shared genes indicating coordinate regulation of these complexes ( Fig. 3 and 4

in main paper).
Transcriptomic results could provide the molecular detail of how various protein complexes participate in the forming the extracellular and intracellular processes for differentiation to occur with the caveat that transcription is not equivalent to protein levels 34 .
Additional networks that are altered related to cell structure include the Notch signaling-2 (NAS = 0.18) and Notch3 (NAS = 0.18) pathways, which were both down regulated by 50% relative to the 6 hour time point. The role of Notch signaling is still to be fully defined across adipogenesis (e.g., 35  While anti-inflammatory cytokines including IL1RA and IL13 (which were not DEG) have been shown to be augmented at earlier stages of differentiation (day 7 compared to days 14 and 21) 38 , the specific early (48h) changes in IL-4 and IL-5 anti-inflammatory pathways during adipogenesis have not been previously described. The TGF-β receptor signaling pathway (NAS = 0.15) was upregulated consistent with results that show TGF-β inhibits adipogenesis 39 .
A notable difference observed when comparing 48 hours versus 6 hours (time lapse) and 48 hours and 0 (time point) is the absences of sub-networks involved in translational machinery, indicating that these pathways were continuously expressed across this time period.

Time Lapse: 96 versus 48 hours
Differential gene expression decreased between 96 hours and 48 hours with 238 of 416 genes being up-regulated (57%) across 70 pathways (Fig. 3 in main paper). The up-regulated leptin pathway had the highest NAS (0.28) during this period (Fig. 3 and 4 Table 6) overlap and connect to cell signaling or are regulated by genes involved in transcription. These transcriptional changes are consistent with membrane restructuring that occurs during differentiation 34 . Table 6 Differentially Expressed Genes involved in Cell Membrane Function  SOS1 plays key roles in growth regulation through the RAS pathway and in regulating the structure of the cytoskeleton through interactions with actin 45 but has not been well studied in adipogenesis or the mature adipocyte.

Additional Figures and Tables
Supplement Figure 5. The plots of the empirical standard deviation versus the rank of the mean computed for the list of 19 variance stabilization / normalization methods tested. The plots were obtained using the diagnostic functions in the VSN package 26 .  Table 8. The combinations of transformation and normalization methods tested to select the best one for our dataset. Methods Implementation log2 + scale normalization limma log2 + quantile normalization limma log2 + cyclic loess normalization limma log2 + quantile normalization beadarray log2 + cubic splines normalization beadarray log2 + scale normalization beadarray log2 + rank invariant normalization beadarray vsn normalization beadarray vsn normalization vsn log2 + quantile normalization lumi log2 + robust spline normalization lumi log2 + simple scaling normalization lumi log2 + loess normalization lumi log2 + rank invariant normalization lumi cubic root + quantile normalization lumi cubic root + robust spline normalization lumi cubic root + simple scaling normalization lumi cubic root + loess normalization lumi

Summary
NASFinder was used to provide a systems-wide yet more detailed analysis of gene expression at 6,48,96,192, and 384 hours after induction of differentiation in pre-adipocytes. Developing methods and tools for systems analysis should provide more comprehensive understanding of the behavior of cellular processes of complex systems over time.