Integrative network analysis reveals different pathophysiological mechanisms of insulin resistance among Caucasians and African Americans

Background African Americans (AA) have more pronounced insulin resistance and higher insulin secretion than European Americans (Caucasians or CA) when matched for age, gender, and body mass index (BMI). We hypothesize that physiological differences (including insulin sensitivity [SI]) between CAs and AAs can be explained by co-regulated gene networks in tissues involved in glucose homeostasis. Methods We performed integrative gene network analyses of transcriptomic data in subcutaneous adipose tissue of 99 CA and 37 AA subjects metabolically characterized as non-diabetic, with a range of SI and BMI values. Results Transcripts negatively correlated with SI in only the CA or AA subjects were enriched for inflammatory response genes and integrin-signaling genes, respectively. A sub-network (module) with TYROBP as a hub enriched for genes involved in inflammatory response (corrected p = 1.7E-26) was negatively correlated with SI (r = −0.426, p = 4.95E-04) in CA subjects. SI was positively correlated with transcript modules enriched for mitochondrial metabolism in both groups. Several SI-associated co-expressed modules were enriched for genes differentially expressed between groups. Two modules involved in immune response to viral infections and function of adherens junction, are significantly correlated with SI only in CAs. Five modules involved in drug/intracellular transport and oxidoreductase activity, among other activities, are correlated with SI only in AAs. Furthermore, we identified driver genes of these race-specific SI-associated modules. Conclusions SI-associated transcriptional networks that were deranged predominantly in one ethnic group may explain the distinctive physiological features of glucose homeostasis among AA subjects. Electronic supplementary material The online version of this article (doi:10.1186/s12920-015-0078-0) contains supplementary material, which is available to authorized users.


Background
Ethnic differences exist in prevalence of common metabolic disorders, including type 2 diabetes (T2D) and metabolic syndrome [1]. In the United States, the recently estimated age-adjusted prevalence of T2D in adults is 7.6% among non-Hispanic Caucasians (European Americans or CA), and 14.9% among non-Hispanic African Americans (AA) [2]. Clinical (visceral adiposity and waist circumference) and biochemical (markers of inflammation) factors only partly explain these ethnic disparities [3][4][5]. The physiology of glucose homeostasis in African Americans is distinctive, with more pronounced insulin resistance, higher insulin secretion, and reduced hepatic insulin clearance compared to CAs [6]. These ethnic differences are significant even after accounting for age, gender, and obesity. Differences in evolutionary history have distinctly shaped the genetic (and epigenetic) architecture of AA and CA populations and may explain the ethnic differences in the stabilization points (canalization) of insulin sensitivity and other glucose homeostasis phenotypes. A recent study suggested that this ethnic stratification could be implicated in the different natural courses of T2D onset and may lead to differences in T2D prevalence [7]. Identifying molecular mechanisms that can explain the observed ethnic differences in insulin sensitivity will be helpful in developing suitable prevention, surveillance, and treatment strategies for the pathogenesis of T2D. However, the genetic cause or the physiological mechanisms for greater insulin resistance among AAs are unknown [8].
Transcript expression in human tissues is correlated with insulin sensitivity and linked to dysregulation of genes involved in the pathophysiology of obesity, insulin resistance, and T2D [9]. Lymphoblastoid cell line transcript profiles of different ethnic groups from the Hap-Map project also identified significant ethnic differences in expression of transcripts [10][11][12][13]. Comparisons of gene expression in subcutaneous adipose and skeletal muscle from a small cohort of CA and AA subjects by our laboratory and others indicated ethnic differences in transcript expression; some of these transcripts are associated with insulin sensitivity [14,15]. However, no analyses have identified transcriptional mechanisms associated with insulin sensitivity and related metabolic traits that predominate in one ethnic group compared to another.
Expression of transcripts involved in the same biological function tend to be co-regulated by similar factors (genetic or environmental) and can be identified as distinct network modules, where genes within a module are more highly interconnected (correlated) with each other than genes in other modules [16]. Thus, we hypothesize that distinctive physiological differences, including differences in insulin sensitivity (S I ) between CAs and AAs, can be explained at least partially by gene expression network modules in tissues involved in glucose homeostasis. These modules may show characteristic patterns of insulin resistance-associated derangements predominantly in one ethnic group compared to others.
In this study, we systematically analyzed genome-wide expression data from subcutaneous adipose tissue of 99 CA and 37 AA non-diabetic subjects. A pressing task in analyzing omics data like those generated by this study is to identify and visualize a global landscape of interactomes that contribute to clinical endpoints such as disease onset and progression. Weighted gene coexpression network analysis (WGCNA) has emerged as a way to solve this problem through identification of gene modules comprised of highly interconnected genes over a gene-gene interaction heatmap [17]. WGCNA has been widely used to identify pathways and gene targets for a variety of common human diseases such as cancer [18,19], atherosclerosis [20,21], Alzheimer's disease [16], obesity and diabetes [22,23]. A key contribution of this work is the integration of differential expression analysis, gene-trait correlation analysis, and gene network analysis to identify both the coexpressed gene modules that contribute to insulin sensitivity and the key causal regulators. Utilizing this integrative network approach, we identified transcript co-expression modules correlated with S I and other related metabolic phenotypes in each ethnic group. Several modules enriched for important biological pathway genes showed differential connectivity between ethnic groups. Our study successfully delineated transcriptional networks associated with insulin sensitivity and other glucose homeostasis traits that are deranged predominantly in one ethnic group, and may explain differences in insulin sensitivity between CAs and AAs.

Ethics statement
All participants provided written informed consent under protocols originally approved (IRB #53028) by the University of Arkansas for Medical Sciences (UAMS). Our current study was approved by the Institutional Review Board at Wake Forest School of Medicine (IRB#00011083).

Study subjects
We used subcutaneous adipose tissue samples from 136 non-diabetic individuals. These subjects were recruited previously; detailed methods for recruitment, physical examination, physiological experiments, and obtaining biopsies have been published [24]. In brief, Caucasian (European-American) or African-American men and women who were generally in good health, between 19 and 60 years of age, and had a BMI between 19 and 45 kg/m 2 were recruited. All participants had a screening visit, at which time height, weight, and waist and hip circumference were measured; body fat was determined by dual x-ray absorptiometry (DXA) scan (Hologic QDR-4500); fasting blood samples for lipid measurements were taken; and a standard 75-g oral glucose tolerance test (OGTT) was done, with measurement of glucose and insulin at baseline and 30-minute intervals for 2 hours. OGTT data were analyzed by homeostatic model assessment [25]. Based on their OGTT results, individuals who did not have diabetes were selected for this study.
At the second visit, adipose biopsies were obtained from subjects under fasting conditions. Insulin modified (0.03 U/kg) frequently sampled intravenous glucose tolerance tests (FSIGT) were performed as described previously [14,24]. We used the MINMOD Millennium program to analyze FSIGT data to determine insulin sensitivity (S I ) and acute insulin response (AIR G ) [26]. Biopsies were obtained from abdominal subcutaneous fat near the umbilicus under local anesthesia (Lidocaine) using a Bergstrom needle. Each biopsy is made by opening the Bergstrom needle and applying suction using a 120 ml syringe. Samples were rinsed immediately in sterile normal saline, quick-frozen in liquid nitrogen, and stored at −80°C for further use. Demographic characteristics of our cohort are shown in Table 1.

Laboratory measurements
Insulin was measured by the UAMS Clinical Research Center core laboratory using an immuno-chemiluminometric assay (Invitron Limited, Monmouth, UK). Plasma glucose was measured by glucose oxidase methods at LabCorp, Inc. (Burlington, NC). Plasma triglyceride, total cholesterol, and HDL cholesterol concentrations were measured directly (at LabCorp) by enzymatic colorimetric methods using an automated Hitachi Cobas C system and assay kits from Roche. LDL cholesterol and VLDL cholesterol concentrations were calculated indirectly by the Friedewald equation.

RNA extraction
Total RNA was isolated from whole adipose tissue using the RNAeasy Lipid Tissue Mini kit (QIAGEN Inc-USA, Valencia, CA). The quantity and quality of the isolated total RNA samples were determined by ultraviolet spectrophotometry (Nanodrop, Thermo Scientific, Pittsburgh, PA) and electrophoresis (Experion nucleic acid analyzer, BioRad Laboratories, Inc., Hercules, CA), respectively. High-quality RNA with RIN (RNA integrity number) >8 was used for genome-wide transcriptome analysis.

Microarray studies
Genome-wide transcriptome analysis and initial array processing were done at the Center for Human Genomics Core Laboratory (Wake Forest School of Medicine) using HumanHT-12 v4 Expression BeadChip (Illumina, San Diego, CA) whole genome gene expression arrays according to the vendor-recommended standard protocol [24]. In brief, biotin-labeled cRNA was prepared from 500 ng of total RNA using the Illumina TotalPrep RNA amplification kit, and 750 ng of biotinylated cRNA was hybridized to the chip in a hybridization oven at 58°C for 18 hours. Chips were processed through a series of posthybridization washes, blocked in blocking buffer, and moved to streptavidin-Cy3-containing detection buffer. Chips were then washed, dried, and scanned in the Illumina BeadArray Reader. All RNA samples were processed using the same lot of chips and reagents, and all experiments were performed within a month to avoid technical variability. Raw expression intensity was background subtracted and normalized by the average normalization algorithm as implemented in GenomeStudio Gene Expression Module v1.0 application software (Illumina). Normalized data were used for further analysis. Data for 17,434 probes annotated as an Entrez gene and with a detection p-value of ≤0.01 in 90% of the samples were used for final analysis. All expression data has been deposited in Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih. gov/projects/geo/) under accession number GSE65221.

Real-time qPCR validation
Total RNA (1 μg) was reverse transcribed using a Qiagen QuantiTect reverse transcription kit (QIAGEN Inc., Valencia, CA) according to the manufacturer's protocol. Transcript-specific oligonucleotide primers for quantitative real-time PCR (qRT-PCR) were designed to capture splice variants captured by the probe sequences in the array [14,24]. Selected transcripts were measured by qRT-PCR using Power SYBR green chemistry (Applied Biosystems, Inc., Foster City, CA) and normalized to the expression of human ribosomal protein, large P0 (RPLP0 or 36B4) gene [24]. The standard curves were generated for absolute quantification using pooled cDNA from the samples assayed.

Statistical and bioinformatic data analysis
To identify differential expression of individual transcripts between CA and AA subjects, normalized geneexpression data from adipose tissue were first adjusted using linear regressions for age and gender (expressioñ age + gender). The residuals following regression were analyzed by Student's T-test to identify transcripts with differential expression between CA and AA samples. Differentially expressed transcripts were defined as those with a ≥1.2 fold change and a p value ≤ 0.05. The residual expression values after covariate adjustments (for age and gender) were used for further testing, including network construction and gene-expression physiological phenotypic correlations. We used Spearman's correlation coefficient to identify correlation of each transcript with S I . False discovery rates (FDR) were also calculated; a transcript was considered associated with S I if FDR ≤0.05.
The CA and AA data sets from all 17,434 expressed probes were independently processed through the weighted gene co-expression network analysis (WGCNA) [17,27]. The weighted network analysis began with a matrix of the Pearson correlations between all gene pairs, and then converted the correlation matrix into an unsigned adjacency matrix using a power function, so that the resulting adjacency matrix, i.e., the weighted coexpression network, is approximately scale-free. To explore the modular structures of the co-expression network, the adjacency matrix was further transformed into a topological overlap matrix (TOM) [17]. Because topological overlap between two genes reflects both their direct interaction and their indirect interactions through all other genes in the network, this approach helps create more cohesive and biologically more meaningful modules. To identify modules of highly coregulated genes, we used average linkage hierarchical clustering to group genes based on the topological overlap of their connectivity, followed by a dynamic cut-tree algorithm to dynamically cut clustering dendrogram branches into gene modules. To distinguish between modules, each module was assigned a unique color identifier. The gene with the highest intramodular connectivity was considered the "hub gene" (most connected gene).
To examine how each gene module was related to gluco-metabolic traits (e.g., S I ) we first performed principal component analysis for each module, and then computed module-trait correlation between the first principal component (module eigengene) and each trait [16,28]. The significance (p value) and false discovery rate (FDR) of each correlation were also calculated. The FDR was estimated through random permutation of sample names of the trait data. A module was associated with a trait if the FDR was ≤0.05. Overrepresentation of canonical pathways and biological processes in modules was measured through Fisher's exact test and corrected for numbers of modules and functional categories tested. Ingenuity Pathway Analysis (IPA) also was used to reveal enrichment of canonical pathways within selected modules.
To quantify differences in adipose tissue transcript network organization between CA and AA subjects, we employed a modular differential connectivity (MDC) metric [16]. In brief, MDC is the ratios of the connectivity of all gene pairs in a module from CA subjects, to that of the same gene pairs from AA subjects. MDC is a continuous measure ranging from 0 to infinity. MDC > 1 indicates gain of connectivity or enhanced co-regulation between genes, whereas MDC < 1 indicates loss of connectivity or reduced co-regulation between genes. The statistical significance of the MDC metrics was computed through the FDR procedure. The significance or FDR of the MDC statistic can be accessed by permuting the data underlying the two networks. We estimated FDR based both on shuffled samples (i.e., networks with nonrandom nodes but random connections) and shuffled gene labels (i.e., networks with random nodes but nonrandom connections), and then selected the larger value as the final FDR estimate.
To further identify key regulators (drivers) genes of the modules identified by WGCNA, we applied the key driver analysis [16] to the module-based unweighted coexpression networks derived from ARACNE (Algorithm for the Reconstruction of Accurate Cellular Networks) [29]. ARACNE was used first to identify significant interactions between genes in each module based on their mutual information, and then remove indirect interactions through data processing inequality (DPI). For each ARACNE-derived unweighted network, we further identified the key regulators by examining the number of N-hob neighborhood nodes (NHNN) for each gene. For a given network, let μ be the numbers of NHNNs and d be the out degrees for all the genes. Genes with a number of NHNNs greater than μ þ σ μ ð Þ are nominated as key regulators [16]. This criterion identified genes with number of NHNNs significantly above the corresponding average value.
We also tested conservation of S I -associated modules between ethnic groups by analyzing overlap of module members and enrichment of modules for common or ethnic-specific S I -associated genes. Fisher's exact test was used to quantitate significant overlap and enrichment. Modules in CAs that show significant overlap (FET_P ≤ 3.1E-06, i.e., Bonferroni corrected p value of 0.05) with a module in AAs and are enriched for S I -associated genes in both ethnicities were considered conserved, while modules enriched for genes associated with S I in only one ethnic group were considered ethnic-specific.
We used two different methods to control for multiple testing issues: (1) the highly conservative Bonferroni correction method, to constrain the study-wise significance level; and (2) an empirical false discovery rate (FDR) method, which constrains the overall rate of false positive events. For instance, for the enrichment of functional categories in modules, we corrected for the number of modules and Gene Ontology processes and other categories tested. For the enrichment of expression clinical trait correlation signatures in modules, we used a combination of FDR <0.05 and nominal P < 0.05.

Adipose tissue transcripts show ethnic differences in expression
We identified significant differential expression (fold change 1.2 and p < 0.05) of 655 transcript probes between CA and AA subjects. These differentially expressed probes included 225 up-regulated and 193 down-regulated mRNA transcripts (Additional file 1: Table S1) of RefSeq genes with known function. Genes with >2-fold higher expression in AA subjects included crystallin-beta 2 (CRYBB2), splicing factor 1 (SF1), and C-type lectin domain family 4, member A (CLEC4A). Genes with >2 fold lower expression included diabetes-associated protein in insulin-sensitive tissue (DAPIT/USMG5), macrophage receptor with collagenous structure (MARCO), and dehydrogenase/reductase member 4 like 1 (DHRS4L1) ( Figure 1A). We also validated differential expression of MARCO by qRT-PCR ( Figure 1B) and detected a positive correlation of this gene with BMI in CA subjects ( Figure 1C). Annotation of these differentially expressed genes using DAVID [30] indicated overrepresentation of these genes in important functional categories (Additional file 1: Table S2). Highly represented categories included genes involved in signal transduction and cell communication, endoplasmic reticulum, regulation of blood pressure, and the carboxylic acid biosynthetic process. In IPA analyses, enrichment of genes involved in glutathione-mediated detoxification and NRF2-mediated oxidative stress response pathways were among those differentially expressed in AA subjects.
We identified significant correlations of 1,625 transcripts with insulin sensitivity (S I ) in both CA and AA subjects (812 and 813 transcripts showed positive and negative correlations, respectively with S I ). Some transcripts showed significant correlation (FDR ≤ 0.05) with S I only in CA (n = 2,797 transcripts) or only in AA subjects (n = 1,983 transcripts) ( Figure 2). These transcripts were enriched for important gene ontology categories and pathways (Additional file 1: Table S3 and Figure 2). Transcripts negatively correlated with S I only in CA subjects (n = 1,367 transcripts) were enriched for inflammatory response pathway genes (75 genes, corrected p-value = 6.90E-09). Transcripts negatively correlated with S I only in AA subjects (n = 932 transcripts) were enriched for integrin signaling pathway genes (24 genes, corrected p-value = 8.90E-06).

Co-expression modules in adipose are associated with S I and related metabolic phenotypes in Caucasians
Network distance coupled with hierarchical clustering and dynamic tree cutting methods implemented in WGCNA identified 93 co-expression networks (modules) in adipose tissue of CA subjects (module sizes from 1,957 to 10 transcripts). Modules are denoted by colors ( Figure 3A). Four modules in CA subjects included >1,000 probes, the largest being the "turquoise" module, with 1,957 transcript probes. Transcripts that could not be included in any module (3,777 of 17,434 probes) were classified as "gray". Module membership of each transcript, simple adjacency matrix-based connectivity, and TOM-based connectivity measures of those transcripts within each module are shown in Additional file 1: Table S4.
Correlation analysis of module eigengenes (principal components of the module) with gluco-and cardiometabolic traits of CA subjects identified some strong module-trait correlations. As shown in Figure 3B, some modules were correlated with multiple phenotypes (also see Additional file 1: Table S5). The clustering of S I and triglyceride (TG) indicates a significant overlap of modules correlated with these two traits in CA subjects.
Modules showing significant positive correlation with S I showed significant negative correlation with TG in CA subjects. Module eigengene of 27 modules showed significant correlation with S I (p ≤ 0.05 and FDR ≤ 0.05) ( Figure 3C). Eight large modules (>25 gene members) showed strong correlation (positive or negative correlation ≥0.4) with either FSIGT-derived (S I and AIRg) [26] or OGTT-derived (HOMA-IR and Matsuda index) [25] glucose homeostasis traits ( Table 2). The "yellow" module showed strong negative correlation (r = −0.426, P = 4.95E-04) with S I and was strongly enriched for genes involved in inflammatory response (corrected p = 1.7E-26). The TYR-OBP gene (TYRO protein tyrosine kinase binding protein) was the most connected hub gene of this co-expression network module based on TOM-based intramodular connectivity, consistent with the finding of the prominent role of TYROBP in other complex diseases [16,31].
The "pink" and "light yellow" modules, which are both enriched for mitochondrion genes are positively correlated with S I with r = 0.411 (P = 3.53E-04) and 0.413 (P = 7.3E-06), respectively. Two nuclear-encoded mitochondrial genes, SDHB (Succinate dehydrogenase complex, subunit B) and DCI (Dodecenoyl-Coenzyme A delta isomerase or ECI1), were the hub genes for these two modules. The "light yellow" module included genes involved in carboxylic acid metabolism (corrected p = 5.7E-06), while the "pink" module was enriched for branched-chain amino acid (valine, leucine, and isoleucine) degradation genes (corrected p = 2.70E-14, see Additional file 1: Table S6).
Co-expression modules in adipose are associated with S I and related metabolic phenotypes in African Americans WGCNA of adipose tissue transcripts of 37 African American subjects using the same parameters as described above for Caucasians (beta = 5, scale R^2 = 0.81, trunc. R^2 = 0.99, hiercutoff =0.99 and minModuleSize =10, Additional file 2: Figure S1) identified 173 co-expression modules (module size ranges from 1980 to 10 probes, Figure 4A, Additional file 1: Table S7). Three modules in AA subjects included >1,000 probes, the largest being the "turquoise" module, with 1980 transcript probes. In AA subjects, clustering of module-trait correlations indicated the strongest overlap between S I -and AIRgassociated modules ( Figure 4B). S I was correlated with the eigengene of 34 modules (p ≤ 0.05 and FDR ≤ 0.05) ( Figure 4C and Additional file 1: Table S8). Sixteen large modules showed strong correlations (positive or negative correlation ≥0.4) with either FSIGT-or OGTTderived glucose homeostasis traits in AA subjects  Table S3 for further details.
To validate our findings, we compared our co-expression modules with the coexpression network modules from an independent human adipose gene expression data (called deCODE data) from 640 Caucasian subjects [23]. Twelve modules in the AA network and nineteen modules in the CA network significantly (Bonferroni corrected FET

Some metabolic trait-associated modules show ethnic differences
We tested the enrichment of modules for genes differentially expressed in AA compared to CA subjects. Nineteen co-expression modules in adipose of CA subjects were enriched for differentially expressed genes (enrichment p ≤0.05). Among these modules, five showed significant correlations with S I (Additional file 1: Table S10A).   The "red" module in CA subjects included the largest number of genes (50 genes, p =3.17E-26) upregulated in African Americans, correlated with AIR G (r = 0.276, p = 0.0476), and was enriched for genes involved in cytoskeletal protein binding and adherens junction (enrichment p-value <3.0E-10). Similarly, 18 transcript modules in adipose of AA subjects were enriched for genes differentially expressed in African Americans compared to Caucasians, and five were significantly correlated with S I (Additional file 1: Table S10B). The "magenta" network module in AA adipose tissue was enriched for genes differentially up-regulated in AA compared to CA subjects (p = 0.0084, includes 11 up-regulated transcripts), showed a significant positive correlation with AIR G (r = 0.51; p = 0.006), and was enriched for genes involved in unsaturated fatty acid biosynthesis (corrected p = 0.0037) and glucose metabolism (corrected p = 0.007). We used the modular differential connectivity (MDC) metric [16] to quantify differences in adipose tissue transcript network organization between CA and AA subjects. This analysis suggested significant differential connectivity (FDR < 0.2) in 59 CA adipose tissue-derived modules (Additional file 1: Table S11), almost all showing gain of connectivity (MDC > 1). Among these modules, 13 showed significant correlation with S I in CA subjects and were enriched for genes in functional categories important for insulin resistance (Table 4). Similar analyses of AA adipose tissue-derived networks identified MDC in 11 modules; 10 had high gain of connectivity (MDC> > 1) and 1 had loss of connectivity (MDC < 1). Only one of these modules was correlated with insulin sensitivity (Additional file 1: Table S12).
We identified 16 and 15 conserved S I -associated network modules in CAs and AAs, respectively (FDR ≤ 0.05) (Table 5A and B, Figures 3C and 4C). These modules showed significant overlap (FET_p value ≤3.1E-06) of module gene members in both ethnicities, and were significantly enriched for genes correlated with S I in both ethnicities. We also identified five and nine ethnicspecific S I -associated network modules, respectively, in CAs and AAs (FDR ≤ 0.05) ( Table 6). These modules were significantly enriched for genes correlated with S I in only one ethnicity (Figures 3C and 4C)  modules (three in CA and five in AA) also showed no significant overlap (FET_p value > 3.1E-06) of gene members between ethnicities, which indicates the presence of highly ethnic-specific mechanisms associated with insulin sensitivity. Among these ethnicity-specific S I -associated modules, "bisque" is most significantly enriched for transcripts negatively correlated with S I in CAs (p =1.85E-16). This module, which has a strong differential connectivity (MDC = 11.16), is enriched for the interferon signaling pathway (in IPA analysis B-H P-value = 5.01E-20, Figure 6A, B). Several interferon-inducible genes involved in immune response to viral infection (including OAS2, ISG15, STAT1, and ADAR) were identified as "drivers" in this sub-network ( Figure 6C) based on the ARACNE [29] and key driver analysis [16,31]. The pro-inflammatory T-cell-derived cytokine interferon gamma (IFNγ) promotes the macrophage M1 phenotype. IFNγ also modulates adipocyte metabolic functions. McGillicuddy et al. (2009) demonstrated that IFNγ may play an important role in T-cell modulation of dietinduced obesity, insulin resistance, and type 2 diabetes via activation of the adipocyte JAK-STAT pathway [32]. IFNγ attenuates insulin signaling, lipid storage, and differentiation in human adipocytes, an effect likely mediated via sustained STAT1 activation [33]. Indeed, 19 genes of this module were associated with metabolic and cardiovascular disorders based on the Genetic Association Database (http://geneticassociationdb.nih.gov/). Specifically, IRF7 and NLRC5 are associated with type 2 diabetes, and SP110 and ZBP1 are associated with waist circumference and body fat distribution, respectively.
Among the 9 AA-specific and S I -associated modules, the "light yellow" module is most enriched for the transcripts positively correlated with S I in AAs (p = 6.95E-09).    The top driver gene of this sub-network ( Figure 7A), YWHAE (14-3-3-epsilon protein), may play a role in the regulation of S I by modulating the interaction between the insulin receptor and IRS1. This module was enriched for oxidoreductase activity (p = 2.90E-04), and was marginally enriched for genes in ERK5 (p = 0.007) and vascular endothelial growth factor (VEGF) signaling pathways (p = 0.018) ( Figure 7B) on IPA analysis. Both pathways include the YWHAE gene. The top interaction network based on the IPA database also includes YWHAE (14-3-3) ( Figure 7C). We identified 35 genes from this module in the genetic association database, including genes associated with type 2 diabetes, obesity, and insulin resistance (i.e. PNPLA2, SOD2, ZNF716, VAMP3, SHC1, IL18, DLC1, CYFIP1, and ASPSCR1).

Discussion
The primary goal of our study was to measure transcriptional modulation related to S I and related metabolic traits in adipose tissue, and to examine how these might relate to molecular mechanisms of insulin resistance in African Americans compared to European Americans. Our analyses revealed distinct ethnic-specific, S I -associated configurations and biological functions of these co-expression modules. Our study provides critical insight into the transcriptional mechanisms and molecular processes underlying Figure 6 The "bisque" module is a Caucasian-specific insulin sensitivity (S I )-associated network module. A) Enrichment of the "bisque" module genes for canonical pathways based on IPA analysis, B) Genes in this module are involved in the interferon signaling pathway, and are negatively correlated with S I (shown in green color), and C) Network plots of the "bisque" module in adipose of CA subjects. The red nodes are "drivers"; the size of a node is proportional to the number of nodes in the node's 2-layer neighborhood.
insulin sensitivity that could be implicated in T2D onset and may explain differences in T2D prevalence between CA and AA individuals.
We identified significant differential expression of 655 transcript probes between CA and AA subjects, which represent 418 genes with known function. A significant enrichment of these genes in important functional categories, including the glutathione-mediated detoxification and oxidative stress response pathways, was observed. Several isoforms of glutathione S-transferases (GSTT2, GSTT2B, GSTM3, and GSTM5) showed higher expression, while microsomal glutathione S-transferase 1 (MGST1) and superoxide dismutase 2 (SOD2) showed lower expression in AA subjects compared to CA subjects. The increased GST transcript expression we observed may indicate a stronger oxidative stress response capacity (in response to environmental factors) [34] in African Americans. However, further molecular genetic studies will be required to define its role in relation to insulin resistance and obesity.
Intriguingly, AA subjects show lower expression of a transcript splice form of USMG5 (DAPIT). DAPIT is a phylogenetically conserved ATP synthase (mitochondrial and lysosomal) that was down-regulated at the mRNA level in skeletal muscle of streptozotocin (STZ)-treated diabetic rats [35]. Its protein expression was up-regulated in epididymal adipose tissue of STZ-treated mice [36]. However, this gene was not associated with insulin sensitivity or BMI in our study. Compared to CA subjects, AA subjects also show significantly lower expression of the MARCO gene, which was positively correlated with BMI in CA subjects. The MARCO gene is a scavenger receptor expressed on macrophages and dendritic cells. Macrophages from MARCO-deficient mice exhibited lower interleukin-12 (IL-12) production in responses to stimulation with lipopolysaccharide and interferon-γ. MARCO −/− mice have severely impaired ability to clear bacteria from the lungs and increased mortality during bacterial infection, possibly due to impaired Th1 polarization [37,38]. Further studies will be required Figure 7 The "light yellow" module is an African American-specific insulin sensitivity (S I )-associated network module. A) Network plots of the "light yellow" module in adipose of AA subjects. The red nodes are "drivers"; the size of a node is proportional to the number of nodes in the node's 2-layer neighborhood, B) Enrichment of module genes for canonical pathways in IPA analysis, C) The most significant network among the module members based on IPA interaction network analysis involves the YWHAE (14-3-3) gene. Genes positively or negatively correlated with S I in this interaction network are shown in red and green, respectively.
to reveal the species-specific and tissue-specific role of the DAPIT and MARCO genes.
We detected associations between several transcript modules, S I, and related metabolic phenotypes. Analysis to identify enrichment of these trait-associated models for functional categories of genes revealed molecular mechanisms that may cause insulin resistance in both ethnic groups, and indicated some pathophysiological mechanisms that are restricted to one ethnic group. The "black" module in adipose of AA subjects overlapped significantly with the "light yellow" module identified in CA subjects. The DCI (or ECI1) was the hub gene for both modules. This gene encodes a key mitochondrial enzyme involved in beta-oxidation of unsaturated fatty acids [39]. These modules also showed strong enrichment for genes involved in mitochondrial carboxylic acid and branched chain amino acid metabolism. Positive correlations of these modules enriched for mitochondrial function with insulin sensitivity, and negative correlations with BMI in both AA and CA subjects, indicate that a common molecular mechanism may be involved in obesity-associated insulin resistance.
Previous studies had indicated a role of inflammatory and immune response pathway genes in S I and obesity. However, these studies were restricted to Caucasian subjects. In line with these findings, we identified negative correlations between S I and the "yellow" module in adipose of CA subjects (r = −0.426, p = 4.95E-04) and a strong positive correlation of this module with BMI (r = 0.607, p = 4.03E-09). This module was strongly enriched for genes involved in inflammatory response (corrected p = 1.7E-26). IPA indicated enrichment of genes involved in Fcγ receptor-mediated phagocytosis in macrophages and monocytes (p = 2.51E-12). The TYROBP is the hub gene for this transcript network module. This gene encodes a transmembrane signaling polypeptide containing an immunoreceptor tyrosinebased activation motif in its cytoplasmic domain. TYROBP is involved in signaling events that mediate crosstalk between dendritic cells and natural killer cells. Recently, it was identified as a central regulator of a gene expression network in brain tissue involved in late-onset Alzheimer's disease [16,40]. The "yellow" module in CA subjects showed the strongest overlap with the "turquoise" module of AA subjects (452 genes, p = 5.01E-169). This "turquoise" module showed enrichment for immune cell activation and the module eigengene showed significant correlation with S I (r = −0.315, FDR =2.84E-02) in AA subjects. We found enrichment of inflammatory response pathway genes (p-value = 6.90E-09) among genes that negatively correlated with S I only in CA subjects. Interestingly, a recent study in an animal model suggested that local pro-inflammatory response in adipose tissue is adaptive, required for proper adipose tissue remodeling and expansion to enable safe storage of excess nutrient and metabolic homeostasis [41].
Many adipose co-expression modules were enriched for genes differentially expressed between the two ethnicities or had differential connectivity (MDC) among members of the network module. Notably, the "bisque" module showed strong differential connectivity (MDC = 11.16) and was highly enriched for transcripts negatively correlated with S I only in CAs. Driver genes of this subnetwork (OAS2, ISG15, STAT1, and ADAR) are involved in immune response to viral infections, and expression of these genes is inducible by interferon. A recent personal-omics analysis showed transcriptional modulation of genes in pathways involved in glucose regulation of insulin secretion occurred shortly after viral infection, and was predicted as involved in elevation of glucose and HbA1c in a Caucasian subject [42]. A recent in vitro analysis showed significant ADAR genotype (and haplotype)-dependent modulation of IFNγ-producing cells (or IFNγ-Elispot responses upon viral infection) in CAs but not in AAs [43]. We showed that ADAR is one of the 7 CA-specific driver genes in this module. IRF7, another CA-specific driver of this module, was modulated by lipopolysaccharide, interferon, and virus-responsive cisregulatory genetic variants in monocytes and dendritic cells [44,45]. IRF7 modulates expression of several downstream genes in this pathway as a master trans-regulator. The "light yellow" module was enriched for transcripts correlated with S I only in AAs. The top driver gene of this sub-network, YWHAE (14-3-3-epsilon protein), may play a role in the regulation of S I by modulating the interaction between the insulin receptor and IRS1. However, the functions of many other genes in this sub-network are not clear.
Our study cohort included far more Caucasian subjects (N = 99) than African American subjects (N = 37). Unequal racial distribution in our cohort was influenced by the population demographics of Arkansas (recruitment site) and volunteer participation bias. However, to our knowledge, this is the largest single cohort of adipose tissue gene expression related to glucose homeostasis that includes phenotypically well-characterized subjects of both races. While sample size is an important factor for all analyses (e.g. differential expression, gene-trait correlation and network construction), the integrative network analysis was developed in this study to minimize the impact of a relatively small sample size, to obtain a more comprehensive understanding of differentially expressed genes and the network context under which individual genes operate. The highly significant overlap between the results from the two ethnic groups studied here and an independent but much larger study (the deCODE study) in Caucasian subjects suggests that the integrative network analysis approach was appropriate in balancing for the impact of the small sample size.

Conclusions
By comparing co-expression networks in adipose tissue of metabolically well-characterized CA and AA subjects, we found derangements in transcriptional networks of genes involved in mitochondrial energy metabolism pathways that lead to insulin resistance in both ethnicities. We also found evidence that other molecular mechanisms involved in modulating glucose homeostasis are restricted predominantly to subjects from one ethnic origin. Distinct genetic (and epigenetic) regulatory architecture, including differences in the frequency of expression regulatory alleles among populations, may determine the structure of co-expression networks observed in adipose tissues of these subjects [46,47]. Further studies will be required to identify how genetic and epigenetic factors determine the structure of co-expression networks in adipose tissue that modulate glucose homeostasis and related physiological traits.

Additional files
Additional file 1: Table S1. Genes showing differential expression in adipose tissue of African Americans compared to Caucasians. Table S2. Enrichment of functional category and biological pathway among genes showing differential expression in adipose tissue of African Americans compared to Caucasians. Table S3. Enrichment of functional category and biological pathway among genes correlated with insulin sensitivity (S I ) in adipose tissue. Table S4. Weighted gene co-expression network analysis (WGCNA) identified 93 co-expression networks (modules) in adipose tissue of Caucasian subjects. Table S5. Correlation of module eigengenes with gluco-and cardio-metabolic phenotypes in Caucasian subjects. Table S6. Co-expression modules in adipose tissue of Caucasians are enriched for genes in functional categories and biological pathways. Table S7. Weighted gene co-expression network analysis (WGCNA) identified 173 co-expression networks (modules) in adipose tissue of African American subjects. Table S8. Correlation of module eigengenes with gluco-and cardio-metabolic phenotypes in African American subjects. Table S9. Co-expression modules in adipose tissue of African Americans are enriched for genes in functional categories and biological pathways. Table S10. A) Adipose tissue derived transcript modules in CA are enriched for genes differentially expressed in African Americans. B) Adipose tissue derived transcript modules in AA are enriched for genes differentially expressed in African Americans. Table S11. Co-expression modules in adipose tissue of Caucasians show differential connectivity. Table S12. Co-expression modules in adipose tissue of African Americans show differential connectivity.