Transcriptional Regulatory Mechanisms in Adipose and Muscle Tissue Associated with Composite Glucometabolic Phenotypes.

OBJECTIVE: Tissue-specific gene expression is associated with individual metabolic measures. However, these measures may not reflect the true but latent underlying biological phenotype. This study reports gene expression associations with multidimensional glucometabolic characterizations of obesity, glucose homeostasis, and lipid traits.
METHODS: Factor analysis was computed by using orthogonal rotation to construct composite phenotypes (CPs) from 23 traits in 256 African Americans without diabetes. Genome-wide transcript expression data from adipose and muscle were tested for association with CPs, and expression quantitative trait loci (eQTLs) were identified by associations between cis-acting single-nucleotide polymorphisms (SNPs) and gene expression.
RESULTS: The factor analysis identified six CPs. CPs 1 through 6 individually explained 34%, 12%, 9%, 8%, 6%, and 5% of the variation in 23 glucometabolic traits studied. There were 3,994 and 929 CP-associated transcripts identified in adipose and muscle tissue, respectively; CP2 had the largest number of associated transcripts. Pathway analysis identified multiple canonical pathways from the CP-associated transcripts. In adipose and muscle, significant cis-eQTLs were identified for 558 and 164 CP-associated transcripts (q-value < 0.01), respectively.
CONCLUSIONS: Adipose and muscle transcripts comprehensively define pathways involved in regulating glucometabolic disorders. Cis-eQTLs for CP-associated genes may act as primary causal determinants of glucometabolic phenotypes by regulating transcription of key genes.

Introduction

Dysregulation of transcript expression in tissues is linked to the pathophysiology of obesity, insulin resistance, and type 2 diabetes (T2D) (1). Genetic and genomic components of these complex processes are often studied via quantitative endophenotypes. These traits are correlated, and pleotropy among gluco-metabolic endophenotypes has been reported. Several transcriptome-wide analyses identified association of gluco-metabolic traits (e.g. insulin sensitivity, body mass index, % fat mass, HDL-cholesterol) with expression levels of transcripts in human muscle, adipose, liver, pancreatic islet, and blood cells (2-5). These studies were successful in defining biological pathways and mechanisms involved in the pathophysiology of T2D and obesity (1;6). However, individual measures are often modestly correlated and may reflect only a portion of true underlying biological process. Thus, studies focused on single traits may not adequately capture differences in gluco-metabolic phenotypes between individuals similar in one trait but different in others. For example, two individuals may have the same BMI, but their % fat mass, waist-to-hip ratio (WHR) and insulin sensitivity can differ substantially so that each has a different overall metabolic status, translating into differences in disease predisposition (7;8). Approaches that test each endophenotype separately are also liable to reductions in statistical power due to multiple testing penalties. Thus, applying methods that combine correlated endophenotypes into composite phenotypes (CPs) capturing underlying gluco-metabolic constructs are likely to provide novel insight into pathophysiological and molecular processes involved in these disorders.

This study applied a factor analysis (FA)-based dimension reduction approach to combine 23 glucose homeostasis, anthropometric and lipid quantitative traits into a set of uncorrelated gluco-metabolic dimensions or “CPs” in African Americans without diabetes. Using the CPs as outcomes, genome-wide transcript expression data from adipose and muscle tissue were analyzed to identify associated transcripts and biological processes that may molecularly define the gluco-metabolic CPs. The expression quantitative trait (eQTL) analysis integrated genome-wide transcript expression and genotype data to identify CP-associated transcripts whose expression levels are influenced by genetic variants.

Materials and Methods

Human subjects

This study was completed at the Wake Forest School of Medicine (WFSM) Clinical Research Unit, and approved by the WFSM Institutional Review Board. All participants provided written informed consent. The study utilized gluco-metabolic phenotype and multiomic data from 256 unrelated and non-diabetic individuals from the “African American Genetics of Metabolism and Expression” (AAGMEx) cohort (9). Participants were healthy, self-reported African Americans residing in North Carolina aged 18-60 years with a body mass index (BMI) between 18 and 42 kg/m2. A standard 75-g oral glucose tolerance test (OGTT) was used to exclude individuals with diabetes. Height, weight, waist and hip circumference were measured, body fat determined by bioelectrical impedance analyzer, and fasting blood drawn for DNA isolation and biochemical analyses at the screening visit. At a second visit, subcutaneous adipose tissue and skeletal muscle biopsies were collected under overnight fasting condition (see supplementary methods). Frequently sampled intravenous glucose tolerance test (FSIGT) was performed to evaluate insulin sensitivity and secretion by minimal model analysis (10). Clinical, anthropometric, and physiological characteristics of the AAGMEx cohort have been described (9) and are shown in Table 1. Participants had a broad range of gluco-metabolic characteristics suitable for capturing the composite multi-dimensional structure of these phenotypes.

Table 1

Summary of 23 anthropometric, obesity, serum lipid, and glucometabolic phenotypes and orthogonal factor loadings in AAGMEx cohort

Trait [Unit]	N	Mean	StdDev	Rotated Factor Pattern***
				Factor1	Factor2	Factor3	Factor4	Factor5	Factor6
Weight [kg]	256	85.0	18.6	0.885	0.204	0.181	0.085	-0.025	0.305
Height [cm]	256	171.1	10.0	0.149	-0.053	0.175	-0.065	-0.015	0.850
BMI [kg/m2]	256	29.0	5.5	0.898	0.249	0.089	0.130	-0.015	-0.175
Waist [cm]	254	97.0	15.0	0.880	0.255	0.246	0.151	-0.043	-0.026
Hip [cm]	252	105.7	12.0	0.885	0.178	-0.033	0.095	-0.087	-0.174
WHR	252	0.92	0.08	0.385	0.205	0.444	0.147	0.043	0.160
Fat [%]§	232	33.2	9.6	0.501	0.160	-0.051	0.140	-0.097	-0.699
Cholesterol-total [mg/dL]	255	176.8	36.4	0.058	0.044	0.185	0.928	-0.035	-0.073
Triglycerides [mg/dL]	255	84.1	43.0	0.068	0.219	0.829	0.334	-0.135	0.052
HDL-cholesterol [mg/dL]	255	56.1	15.7	-0.406	-0.133	-0.475	0.241	0.203	-0.123
VLDL-cholesterol [mg/dL]	255	16.9	8.7	0.070	0.220	0.828	0.336	-0.137	0.049
LDL-cholesterol [mg/dL]	255	103.8	32.3	0.237	0.055	0.215	0.816	-0.099	-0.034
HbA1c [%]	256	5.6	0.3	0.070	0.417	0.052	0.308	-0.277	0.216
Fasting Insulin(I0) [mu/L]	256	10.1	8.4	0.279	0.842	0.114	-0.030	-0.030	-0.140
Insulin-120min (I120) [mu/L]*	256	62.6	78.0	0.130	0.481	0.281	0.262	-0.213	-0.355
Fasting Glucose(G0) [mg/dl]	256	91.2	9.4	0.076	0.672	0.146	0.084	-0.160	0.406
Glucose-120min(G120) [mg/dl]*	255	100.8	29.8	0.058	0.484	0.144	0.342	-0.356	-0.113
MATSUDA index*	249	6.2	6.7	-0.296	-0.602	-0.123	0.045	0.184	0.208
HOMA-IR	256	2.3	2.0	0.282	0.858	0.160	0.000	-0.039	-0.058
SI [×10-4. (mu/l)ˆ-1.minˆ-1]**	233	4.0	3.3	-0.429	-0.163	-0.212	0.189	0.540	0.324
AIRG [mu.lˆ-1.min]**	233	773.0	641.6	0.403	-0.037	0.425	-0.273	0.316	-0.467
DI**	233	2276.6	1511.5	-0.079	-0.209	0.021	-0.097	0.868	-0.130
SG [minˆ-1]**	233	0.019	0.010	0.042	-0.092	-0.149	-0.094	0.799	0.049

Percent fat mass determined by bioelectrical impedance analyzer;

from 75-gm Oral glucose tolerance test;

from insulin modified (0.03 U/kg) frequently sampled intravenous glucose tolerance test (FSIGT), units are taken from MINMOD Millennium program;

principal component extraction of 23 phenotypes followed by varimax rotation generated orthogonal factor loadings. Factor loadings with an absolute value >0.4 are bolded. Factors are denoted as composite phenotypes (CP) for biological interpretation.

Laboratory measures and physiological phenotypes

see Supplementary Methods

Gene expression analysis and genotyping

see Supplementary methods

Statistical analysess

Quality control

Quality control of phenotype, gene expression and genotype data has been reported (9;11) and briefly described in Supplementary Methods.

Composite Phenotype (Factor) analysis

Values of the 23 gluco-metabolic traits likely reflect latent constructs with shared variation. To capture the various dimensions of the gluco-metabolic traits reflecting these latent constructs, a factor analysis based on the covariance matrix was computed using principal component extraction and varimax rotation via ‘PROC FACTOR’ in SAS; varimax rotation generates orthogonal (independent) factors, denoted here as CPs. All factors with eigenvalues >1.0 were retained and the proportion of variance explained was recorded. From these factors, each representing a latent construct, the factor loadings with an absolute value >0.4 were retained (i.e., <-0.4 or >0.4) and a linear combination (i.e., weighted mean) of these loadings computed (Table 1). The resulting scores for factors 1–5 were natural logarithm transformed, standardized (i.e., subtracting the mean and dividing by the standard deviation) and remaining outliers were winsorized (see supplementary methods). Factor 6 did not require natural logarithm transformation and was standardized and winsorized. Thus, the resulting CPs approximately follow a standard normal distribution and were used in subsequent analyses. These six CPs represent six unique dimensions of the gluco-metabolic domain.

To test for an association between the CPs and expression levels, a linear regression model was computed where the standardized CPs were modeled as the outcome and the log2 of the transcript expression was the predictor of interest. Models included age, gender, and African ancestry proportion as covariates. Admixture estimates were computed using the program ADMIXTURE (12). Expression of a transcript associated with a CP at uncorrected-p<0.001 was considered for subsequent analyses (e.g., pathway analysis).

cis-eQTLs

For transcripts associated with one of the six CPs, a cis-eQTL analysis (i.e., within ±500kb around the respective transcript expressed in ≥90% of participants) was computed. For each transcript associated with a CP, a linear regression was computed with the log2 transformed expression value as the outcome and an additive genetic model for the SNP as implemented in the R-package MatrixEQTL (13), with age, gender, and African ancestry proportion as covariates. Cis-eQTLs with a false discovery rate (FDR)-corrected p-value (Q-value) <0.01 (or 1.0%) were considered significant.

Bioinformatic analysis

see Supplementary methods

Results

Composite phenotypes in AAGMEx cohort

The factor analysis identified six orthogonal CPs (eigenvalues>1.0) that cumulatively explained 74% of the variation in these 23 gluco-metabolic traits (Figure 1A, B). Factors 1 through 6 individually explained 34%, 12%, 9%, 8%, 6% and 5% of the variation. The factor loadings are reported in Figure 1C, and loadings with absolute values greater than 0.4 (i.e., <-0.4 and >0.4) are highlighted in Table 1.

Figure 1

Loadings and explained variance from factor analysis of the gluco-metabolic phenotypes in AAGMEx cohort

Line graphs shows eigenvalues (A) and variance explained (B) by each factor and radar plots (C) show the corresponding factor loadings for 23 phenotypes.

Factor 1 (CP1) exhibited positive loadings for weight, BMI, waist measurement, hip circumference, % fat mass, and acute insulin response (AIRG), and negative loadings for insulin sensitivity (SI) and HDL-cholesterol. Factor 2 (CP2) exhibited positive loadings for fasting and 2hr insulin, fasting and 2hr glucose, HbA1c, and negative loading for OGTT derived insulin sensitivity (Matsuda ISI). Factor 3 (CP3) was defined by positive loadings for waist-to-hip ratio, serum triglyceride (TG), TG-rich very low density lipoprotein (VLDL)-cholesterol and AIRG, and negative loading for HDL-cholesterol. Factor 4 (CP4) captured a cholesterol dimension independent of the TG-based Factor 3, with positive loadings for total and LDL cholesterol. Factor 5 (CP5) exhibited positive loadings for SI, disposition index (DI), and glucose effectiveness (SG). Factor 6 (CP6) was defined by positive loadings for height and fasting glucose and negative loadings for % fat mass and AIRG. Thus, the six CPs partitioned traits into complex constructs.

Transcripts associated with gluco-metabolic CPs

Expression levels of 3994 transcripts in subcutaneous adipose tissue were significantly associated (uncorrected-p<0.001) with at least one of the six CPs (Table S1). CP1 was associated with expression level of 1925 transcripts in adipose (Table 2). Transcripts most strongly associated included ORM1-like protein 3 (ORMDL3/ORMDL sphingolipid biosynthesis regulator 3), transmembrane 7 superfamily member 2 (TM7SF2), and thymocyte nuclear protein 1 (THYN1). In humans, the ORMDL3 gene shows highest transcript level expression in liver and adipose tissue. The ORMDL3 expression in adipose was positively associated with insulin sensitivity (SI, β=0.77, p=7.01×10-8) and negatively associated with BMI (β= -0.86, p=5.50×10-24) in this cohort, and was replicated in an independent study in Caucasians (METSIM cohort, BMI β= -0.429, p=6.79×10-36) (14). In vitro studies in human and mouse cells suggest that downregulation of ORMDL3 increase ceramide, a sphingolipid metabolite involved in inflammatory processes, and potentially involved in pathophysiology of obesity, insulin resistance and asthma (15). Among the 1925 transcripts associated with CP1, 161 were uniquely associated (at p<0.001 threshold), while 1764 were also associated with some of the other CPs (Figure 2). Compared to CP1, CP2 explained a much smaller fraction of total variation (34% vs 12%) in the 23 gluco-metabolic phenotypes. However, of the six CPs, expression levels of the highest number of transcripts (3337 transcripts) were associated with CP2. Fasting insulin levels contributed a high loading (0.842) to CP2 and may have influenced the expression level of adipose transcripts studied after overnight fasting. Transcripts most strongly associated with CP2 include alpha-2-glycoprotein 1 zinc-binding (AZGP1), ubiquitin carboxyl-terminal esterase L1 (UCHL1), and galactosidase, beta 1 (GLB1). The Venn diagram (Figure 3) enumerates the shared transcripts across the CPs. Figure 3 demonstrates the notably larger number of transcripts (1384 Entrez gene, 41.2%) uniquely associated with CP2. The smallest number of adipose transcripts was associated with CP4 (47, p<0.001). At a more stringent threshold of FDR corrected p-value <0.01 no adipose transcript remained significantly associated with CP4 or CP5.

Table 2

Summary of adipose and muscle tissue transcripts associated with the six orthogonal gluco-metabolic phenotype dimensions derived from 23 gluco-metabolic traits

Tissue	Significance threshold	CP1	CP2	CP3	CP4	CP5	CP6
N for Adipose		202	246	225	250	230	207
Adipose*	P <0.001	1925	3337	1579	47	245	499
Adipose*	FDR-P <0.01	1350	3135	984	0	0	43
N for Muscle		198	241	220	245	225	203
Muscle*	P <0.001	299	606	248	62	41	133
Muscle*	FDR-P <0.01	33	197	15	0	2	0
Both tissue[#]		148 (134)	210 (177)	112 (97)	0	5 (5)	41 (37)
Adipose Cis-eGene[§]	Factor association p<0.001 & eQTL q-value<0.01	265	463	218	7	33	55
Muscle Cis-eGene[§]	Factor association p<0.001 & eQTL q-value<0.01	61	106	57	6	9	24

N, number of samples used in analysis (N is variable due to either missing phenotype or expression data),

values in these rows indicate number of transcripts (probes) significantly associated with each composite phenotype (CP) at given threshold.

expression of number of transcripts in both adipose and muscle tissue is associated (p<0.001 in one tissue and P<0.01 in other tissue) with a composite phenotype. Number in parenthesis shows the number of transcripts showing same effect direction (β) in both tissues.

values in these rows indicate number of transcripts (probes) significantly associated with composite phenotype (p<0.001) and are cis-eGenes (FDR <1% or q<0.01).

Figure 2

Distribution of adipose and muscle tissue transcripts associated with top six gluco-metabolic composite phenotypes

Stacked Bar graph shows number of transcripts (represented by probes in Illumina expression arrays) associated (uncorrected-p<0.001) uniquely with each (red) or shared (blue) composite gluco-metabolic phenotype (CP).

Figure 3

Venn diagram showing overlap among adipose and muscle tissue transcripts associated with gluco-metabolic composite phenotypes

Number indicates the count of unique Entrez id genes uniquely associated (uncorrected-p<0.001) with a composite phenotype or overlapping with genes associated with other composite phenotypes.

Compared to adipose tissue, expression levels of fewer skeletal muscle transcripts were associated with CPs. A total of 929 transcripts in muscle were associated (p<0.001) with at least one of the six CPs (Table S2). Among the 299 CP1-associated muscle transcripts, growth factor receptor-bound protein 14 (GRB14), pleckstrin homology-like domain family A member 3 (PHLDA3), and transmembrane protein 192 (TMEM192/FLJ38482) were most strongly associated. The GRB14 transcript level in muscle was positively associated with CP1 (β= 1.27, p= 6.6×10-9). GRB14 protein interacts with insulin receptors and insulin-like growth-factor receptors, and likely has an inhibitory effect on receptor tyrosine kinase signaling and, in particular, on insulin receptor signaling, and may play a role in signaling pathways that regulate growth and glucose metabolism (16). GRB14 knockout mice show improved insulin sensitivity and several genome-wide association studies identified association of SNPs near GRB14 with obesity (waist-to-hip ratio, % fat mass), fasting insulin and T2D (16-18). Similar to adipose, the highest number of muscle tissue transcripts (606 transcripts) was associated with CP2 (Table 2). The smallest numbers of muscle transcripts were associated with CP5 (41 transcript, p<0.001), and only two genes, solute carrier family 25 member 20 (SLC25A20; mitochondrial carnitine/acylcarnitine translocase) and angiopoietin-like 4 (ANGPTL4), remained significantly associated with CP5 at FDR-p<0.01. No muscle transcript remained significantly associated with CP4 or CP6 at FDR-p <0.01.

Expression of a subset of transcripts in both adipose and muscle was associated with CPs. CP1 was associated with 148 transcripts in adipose and muscle (Table 2), with 134 showing directional concordance (increased expression with greater obesity-insulin resistance). For example, expression of GRB14 in both adipose and muscle was positively associated with CP1 (β=1.31, p=1.46×10-5 in adipose and β= 1.27, p=6.6×10-9 in muscle). Similarly, CP2 was associated with 210 transcripts in adipose and muscle, with 177 showing directional concordance and 33 showing directional discordance. Expression level of genes involved in ribosome function (e.g. RPS17, RPL10A, RPL17, RPL22) and translation initiation (EIF2A, EIF3F) in adipose and muscle were negatively associated with CP2. Expression level of 12 transcripts involved in mitochondrial function in adipose (e.g. ECH1, ETFA, ACOT2, CPT2) was negatively (inversely) associated with CP2, while their expression in muscle was positively associated. Expression of five transcripts, CYP1A1, BCKDHB, PER3, SREBF1, and ANGPTL4 in both tissues was significantly associated with CP5 and exhibited the same effect direction. The angiopoietin-like 4 (ANGPTL4) expression level was negatively associated with CP5 (β = -0.77, p =3.27×10-5 in adipose and β =-0.89, p =2.37×10-7 in muscle) that capture efficient insulin sensitivity upon glucose loading. Among the three FSIVGT-derived glucose-homeostasis traits (SI, DI and SG) that define CP5, ANGPTL4 expression was most strongly associated with DI (β= -12.67, p= 5.58×10-7 in muscle and β= -12.5, p=8.61×10-6 in adipose). Studies in mouse models have shown that ANGPTL4 is a glucocorticoid receptor primary target gene that promotes lipolysis in adipocytes, inhibits extracellular lipoprotein lipase, and triggers interorgan communication (19). Increased glucocorticoid level during fasting induces ANGPTL4 expression. ANGPTL4 mediated lipolysis in adipocyte activates ceramide synthesis in the liver and induces whole-body insulin resistance by stimulating the activities of the downstream effectors of ceramide, protein phosphatase 2A and protein kinase Cζ (20).

Pathway enrichment analysis identifies salient biological process linked to gluco-metabolic CPs

Ingenuity pathway analysis (IPA) identified significant enrichment of biological pathways among genes linked to the transcripts associated with the six CPs. Genes annotated in oxidative phosphorylation and mitochondrial dysfunction pathways were enriched among the first three CPs, CP1, CP2 and CP3 -associated adipose transcripts (Figure 4, Table S3). The oxidative phosphorylation pathway was most strongly enriched among CP2-associated adipose transcripts (50 genes, B-H p-value = 1.0×10-15), but was not significantly enriched among CP2-associated muscle transcripts. In adipose tissue, expression of nearly all transcripts in these two pathways were negatively (inversely) associated with CP1, CP2 and CP3. In muscle, genes in oxidative phosphorylation and mitochondrial dysfunction pathways were also enriched among CP1 and CP3-associated transcripts. In contrast to adipose, expression levels of oxidative phosphorylation pathway transcripts in muscle were positively associated with CP1 and CP3 (Table S4). Genes annotated in the EIF2 (Eukaryotic Initiation Factor-2) signaling, a pathway involved in protein synthesis, were most strongly enriched among adipose tissue transcripts associated with CP1 and CP3 (B-H p =2.51×10-7 and 1.15×10-7). The EIF2 signaling pathway was strongly enriched among muscle transcripts associated with CP1, CP2 and CP3 (B-H p =2.0×10-8-3.98×10-22) and transcript profile indicate repression of this pathway (activation z-score<-2; Figure 5). Expression level of most transcripts in this pathway in adipose and muscle was inversely associated with CP1, CP2 and CP3. Genes annotated in pathways regulating translation and cellular metabolic state based on nutrient availability (e.g. Regulation of eIF4 and p70S6K Signaling pathway, and mTOR signaling pathway) were also enriched among adipose and muscle transcripts associated with CP1, CP2 and CP3.

Figure 4

Transcripts in subcutaneous adipose tissue are associated with composite gluco-metabolic phenotypes and enriched for salient biological pathways

A) Heat map shows hierarchical clustering of -log10 p-values for 3994 adipose tissue transcripts (each row indicate probe for a transcript) associated (p<0.001) with composite phenotypes. B) Enrichment and C) activation of genes in biological pathways among six composite phenotype-associated adipose transcripts based on ingenuity (IPA) pathway comparison analysis are shown as heat maps.

Figure 5

Transcripts in skeletal muscle tissue are associated with composite gluco-metabolic phenotypes and enriched for salient biological pathways

A) Heat map shows hierarchical clustering of -log10 p-values for 929 muscle tissue transcripts (each row indicate probe for a transcript) associated (p<0.001) with composite gluco-metabolic phenotypes. B) Enrichment and C) activation of genes in biological pathways among six composite phenotype-associated muscle transcripts based on ingenuity (IPA) pathway comparison analysis are shown as heat maps.

CP-associated transcripts in adipose were also enriched for genes determining fatty acid, amino acid (including branched chain amino acids valine leucine and isoleucine), and bioactive amine concentrations (adrenaline, noradrenaline, serotonin, and dopamine). CP1 had the strongest positive loading for BMI. Corroborating our previous findings on obesity (21), CP1-associated transcripts were enriched for endoplasmic reticulum (ER) stress induced unfolded protein response pathway (11 genes, B-H p = 0.031). CP4-associated adipose and muscle transcripts were not enriched for any biological pathways on IPA analysis or Gene Ontology categories by DAVID analysis. The CP5-associated transcripts in adipose were only marginally enriched for triacylglycerol biosynthesis pathway (B-H p =0.038), while CP5-associated transcripts in muscle were enriched for superoxide radicals degradation (B-H p =0.003) and triacylglycerol biosynthesis (B-H p =0.049). In adipose, five triacylglycerol biosynthesis pathway genes (GPAM, LPIN1, DGAT2, DGAT1, and ELOVL6) were positively associated with CP5, while in muscle, two genes in this pathway (ABHD5 and PLPP1) were negatively associated. DGAT1 is an ER-localized diacylglycerol O-acyltransferase (DGAT) enzyme and during adipocyte lipolysis it mediates triglyceride synthesis by fatty acid re-esterification; this protects the ER from lipotoxic stress and related adipose tissue inflammation (22).

Adipose tissue transcript profiles for CP2-associated genes displayed repression of Rho GDP-dissociation inhibitor (RhoGDI, activation z score = -3.18) signaling, LXR/RXR-activation (z = -2.83), and PPAR signaling (z = -2.23) pathway. The LXR/RXR-activation and PPAR signaling pathway were also repressed among CP1 and CP3-asssociated genes in adipose, while transcript profiles for CP5-associated genes indicate significant activation of the PPAR signaling pathway (z =2). The CP2-associated transcripts indicated strong activation of inflammation-related pathways in adipose, including Fcγ receptor-mediated phagocytosis in macrophages and monocytes (z =4.7), Tec kinase signaling (z =4.24), integrin signaling (z =4.14), TREM1 signaling (z =4.02), leukocyte extravasation signaling (z =4.01), IL-8 signaling (z =3.53), dendritic cell maturation (z = 4.33), and inflammasome pathway (z =2.64). Many of these inflammatory related pathways were also activated among CP1 and CP3-associated adipose transcripts, while CP6-associated transcripts show repression of these pathways (Figure 4).

Expression of a subset of gluco-metabolic CP-associated transcripts is dependent on regulatory genetic polymorphisms

To develop putative causal models, we integrated genotype information (SNPs with MAF≥0.01) and gene expression data through expression quantitative trait (eQTL) analysis to identify cis-regulatory variants (cis-eSNPs) in modulating the expression of CP-associated transcripts in adipose and muscle. In adipose, significant cis-eQTLs were identified for 558 CP-associated transcripts (q-value <0.01, Table S5). In muscle, significant cis-eQTLs were identified for 164 CP-associated transcripts (q-value <0.01; Table S6). Twenty-four CP-associated transcripts were cis-eGenes in both adipose and muscle. Among the CP-associated transcripts, TGF beta-inducible nuclear protein 1 (TINP1/NSA2) had the strongest cis-eSNP in adipose (rs6873912, β=0.391, p=33.47×10-67), while Abelson helper integration site 1 (AHI1) has the strongest cis-eSNP in muscle (rs7772705, β= 0.518, p= 1.63×10-60). Among the adipose cis-eGenes (FDR 1%), dicarbonyl/L-xylulose reductase (DCXR) was most significantly associated with CP2 (β= -2.22, p= 1.11×10-15), while among muscle cis-eGenes prostaglandin D2 synthase (PTGDS) was most significantly associated with CP2 (β= 1.14, p= 1.1×10-11). The top 10 CP-associated cis-eGenes or genetically regulated transcripts in each tissue based on average ranking for CP association p-value and eQTL p-value are shown in Table 3. The top average ranking transcripts included membrane-spanning 4-domains subfamily-A member-6A (MS4A6A) and galectin-related protein (LGALSL/GRP/HSPC159) in adipose and muscle, respectively. The expression of an isoform of MS4A6A (NM_022349.2) in adipose was positively associated with CP2 (β=0.93, p=1.66×10-11) and common minor allele (MAF= 0.33) of SNP rs597982_C was associated with reduced transcript expression (β= -0.45, p= 3.52×10-17).

Table 3

Top gluco-metabolic composite phenotype-associated cis-eGenes in adipose and muscle.

Probe id	Symbol	eQTL							Factor association
		cis-eSNP*	Chr	A1	A2	MAF	beta	P-value	CP§	Beta	P-value
in Adipose
ILMN_1797731	MS4A6A	rs597982	11	C	A	0.333	-0.452	3.52×10-17	2	0.93	1.66×10-11
ILMN_1674069	TOMM7	rs2240726	7	G	A	0.294	0.120	2.16×10-16	2	-3.39	1.74×10-9
ILMN_1665132	CD36	rs3211938	7	G	T	0.100	-0.710	2.11×10-14	2	-0.74	4.33×10-9
ILMN_2072178	ECHDC3	rs200943982	10	T	C	0.401	-0.281	1.94×10-9	2	-1.12	4.85×10-13
ILMN_2364384	PPARG	rs3856806	3	T	C	0.075	-0.754	2.81×10-28	2	-0.95	1.68×10-7
ILMN_1774949	PIGP	rs2298682	21	G	A	0.317	-0.254	1.75×10-24	2	-1.47	1.52×10-7
ILMN_1786105	PCBD1	rs16928023	10	G	A	0.146	-0.220	7.56×10-10	2	-1.75	3.59×10-10
ILMN_1663538	CLYBL	rs2281756	13	G	A	0.266	-0.140	5.14×10-9	2	-2.11	3.75×10-11
ILMN_1720303	OSTM1	rs9372177	6	A	G	0.383	-0.134	8.74×10-10	2	2.09	1.33×10-9
ILMN_1690982	DDT	rs79966373	22	G	C	0.178	-0.302	2.81×10-14	2	-1.25	1.32×10-7
in Muscle
ILMN_1673548	LGALSL	rs10496115	2	C	T	0.253	-0.182	1.87×10-12	2	-1.25	4.01×10-7
ILMN_1666471	UQCRQ	rs36093416	5	G	C	0.016	-0.388	1.38×10-12	2	1.82	1.53×10-6
ILMN_1690125	PDLIM7	rs199995933	5	C	T	0.344	-0.268	1.66×10-22	2	0.80	2.73×10-5
ILMN_1807455	DHRS7	rs384139	14	A	G	0.483	-0.144	3.74×10-9	2	1.11	3.48×10-7
ILMN_1667494	SPTB	rs55845128	14	A	G	0.022	-0.418	1.18×10-8	2	-1.33	4.85×10-8
ILMN_1673788	CDV3	rs115222502	3	T	C	0.020	-0.780	1.01×10-26	2	-1.03	4.09×10-5
ILMN_1675797	EPDR1	rs115983529	7	C	A	0.027	-0.354	1.74×10-7	2	1.42	7.98×10-10
ILMN_1757631	DBNDD1	rs116113385	16	A	T	0.020	0.472	3.58×10-8	2	1.15	4.51×10-7
ILMN_1762747	RPL15	rs9855481	3	C	T	0.229	-0.113	2.08×10-8	3	-1.88	8.17×10-7
ILMN_1807106	LDHA	rs4596	11	C	G	0.151	-0.222	1.12×10-10	3	1.02	2.62×10-5

Top 10 composite phenotype-associated genetically regulated transcripts in each tissue based on average ranking for phenotype association p-value and eQTL p-value are shown.

result for the most significantly associated SNP with the transcript level are shown;

result for most significant association of the transcript with a composite phenotype (CP) are shown.

Expression of genes predicted by a genome-wide association study (GWAS) for BMI are associated with CPs

GWAS in large well powered cohorts identified many - loci associated with increased risk of obesity and other gluco-metabolic traits. For example, Locke et al. (23) identified 97 genome-wide significant (p<5×10-8) loci for BMI. However, most of these trait-associated SNPs are in the non-coding region of the genome, and cannot directly implicate the “culprit- gene”. Thus, Locke et al. (23) used Data-driven Expression prioritization Integration for Complex Traits (DEPICT) (24) to predict and prioritize genes in an expanded set of 511 BMI-associated (p<5×10-4) genomic regions. DEPICT predicted 989 potential causal genes in BMI-associated genomic regions. Among the DEPICT-predicted BMI-genes, expression of 127 and 19 genes in adipose and muscle, respectively, were associated with CPs in our AAGMEx cohort (Table S7).

Discussion

Existing genome-wide transcriptomic studies have tested association of transcript levels in tissues with single anthropometric, glucose homeostasis, and lipid traits (2;4;9;14). Some employed covariate adjustment strategies to account for confounding effects of correlated traits (e.g. SI adjusted for BMI) (3;9). These strategies cannot fully capture variation across multiple traits simultaneously. Specifically, many of these individual measures are partial manifestations of underlying latent gluco-metabolic phenotypes and applying a method that combines correlated endophenotypes into CPs capturing the underlying gluco-metabolic construct is more likely to provide novel insight into the pathophysiological and molecular processes involved in T2D and obesity. This study used factor analysis to identify and partition 23 measures of obesity and glucose metabolism into six orthogonal dimensions of gluco-metabolic CPs. For example, CP1 explained 34% of the variation in the 23 gluco-metabolic measures in this African American cohort (AAGMEx) and comprehensively captured the obesity and FSIVGT-derived glucose-homeostasis phenotypes. Availability of detailed phenotype data for AAGMEx participants enabled us to capture the composite multi-dimensional structure of gluco-metabolic phenotypes. We believe transcripts associated with the CPs most comprehensively define the repertoire of, including novel, biological pathways involved in genetic regulation of gluco-metabolic traits.

Focusing on the top six CPs, a total of 3994 associated transcripts were identified in subcutaneous adipose; only 929 transcripts in muscle were similarly associated. Thus, transcriptional dysregulation involved in determining gluco-metabolic phenotypes appears more pervasive in adipose. Although CP1 (reflecting a composite obesity-insulin resistance phenotype) explained the largest proportion of variation in the 23 measures, expression levels of the largest number of transcripts in adipose and muscle were most strongly associated with CP2 (reflecting a composite hyperinsulinemic-insulin resistance phenotype). Fasting insulin and HOMA-IR index had the largest loadings for CP2. Fasting insulin is higher in insulin-resistant subjects and alters adipose and muscle gene expression and mediates cross-talk between tissues involved in glucose-homeostasis (25). Short term experimental hyperinsulinemia (2hr-hyperinsulinemic euglycemic clamp) induced transcriptional response of 230 genes in adipose of subjects with obesity, and the difference in response was distinct in insulin-sensitive compared to insulin-resistant individuals with obesity (26). Among insulin responsive genes, transcript levels of 45 genes in adipose were associated with CP2 in this study, including genes with known roles in adipose development (RORC, AACS, PPARGC1A, LRP5), AMPK signaling pathway (IRS2, PFKFB3, PPARGC1A, PIK3R1), and T2D (IRS2, DBP, HMOX1, PPARGC1A, PIK3R1, DDIT3, LRP5). In this study, adipose and muscle tissue samples were collected for gene expression analysis following an overnight fast, and participants in this cohort had a broad range of fasting insulin concentrations. Our data suggests a role for plasma insulin level in determining the fasting transcriptome in adipose and muscle. Additional studies are required to resolve temporal and mechanistic connections between hyperinsulinemia, obesity and insulin resistance.

Expression of a subset of transcripts in both adipose and muscle were associated with the six gluco-metabolic CPs. Inverse correlation of genes in pathways involved in protein synthesis (EIF2 signaling), regulation of translation and cellular metabolic state based on nutrient availability (eIF4 and p70S6K Signaling pathway, and mTOR signaling) with CP1 and CP3suggest concordant downregulation of these pathways in adipose and muscle. However, CP1 and CP3-associated genes show discordant regulation of oxidative phosphorylation pathway genes, and CP5-associated genes suggest discordant regulation of triacylglycerol biosynthesis pathway in adipose and muscle. Enrichment of CP2-associated adipose genes in various inflammation-related pathways supports the tissue specific activation of these pathways. Some of the gluco-metabolic CP-associated genes identified here (e.g., ORMDL3, MS4A6A) are involved in asthma and Alzheimer’s disease, suggesting common transcriptional mechanisms across diseases. Precise triggers of adipose tissue inflammation are poorly understood (27); our data supports involvement of multiple potential mechanisms. In adipose, PPAR signaling was repressed among CP1 and CP3-asssociated genes, while transcript profiles for CP5-associated genes indicate significant activation of this pathway. CP5 captures a dimension measuring efficient insulin sensitivity upon intravenous glucose loading, and CP1 reflects a combined obesity-insulin resistance phenotype. Together, these genome-wide transcriptomic and biological pathway analyses define the repertoire of biological pathways involved in regulating distinct dimensions of obesity and glucose homeostasis. Our study used only adipose and muscle tissue to define transcriptional mechanisms determining CPs. Other metabolic tissues are of interest but are not readily accessible in the clinical setting.

Recent GWAS approaches successfully identified genetic loci associated with gluco-metabolic phenotypes; however, identification of precise causal genes in those loci typically remains elusive. Most studies considered the genes closest to the sentinel SNP as the effector gene. For example, FTO was considered the causal gene in the most significant and highly replicated BMI-associated locus on chromosome 16 (28). Recent studies refute this conclusion. Functional genetic analyses, including eQTL and chromatin interaction analysis, suggest that BMI-associated SNPs in the FTO-locus contribute to obesity by regulating expression of the IRX3 and IRX5 genes in pre-adipocytes or brain (29;30). IRX3 is located ~513Kb from the BMI-associated SNPs. In a similar fashion, adipose and muscle transcript levels are key molecular phenotypes associated with composite gluco-metabolic traits, and act proximal to actions of DNA sequence variants. Therefore, the present study focused on identifying transcriptional mechanisms associated with composite gluco-metabolic phenotypes. Our previous studies showed that a subset of gluco-metabolic trait GWAS-identified SNPs are cis-eSNPs (11;31). Herein, we demonstrate that expression of a subset of CP-associated transcripts is determined by cis-eSNPs, and CP-associated transcripts are among the genes predicted by bioinformatics analysis of GWAS-implicated BMI loci. As an alternative to GWAS, this approach provides more direct evidence for putative causal genes and novel genetically-regulated mechanisms determining gluco-metabolic phenotypes.

Our data implicates thousands of genes in biological processes determining gluco-metabolic phenotypes. However, it is likely that a subset of these processes is due to reactive changes in response to primary causal mechanisms. This study cannot conclusively differentiate causal effects from reactive effects based solely on transcriptomic data. Naturally occurring genetic variants, including SNPs in our genome, determine gene expression levels in tissues by controlling transcriptional regulation. Thus, regulatory SNPs may act as primary initiators determining gluco-metabolic phenotypes via roles in modulation of transcript level (SNP →Transcript→ Phenotype) in tissues important for glucose homeostasis. The eQTL analysis in this cohort identified cis-eQTLs for a subset of CP-associated transcripts in adipose and muscle. These CP-associated cis-eGenes may act as key derivers in transcriptional regulatory mechanisms involved in determining gluco-metabolic phenotypes.

Conclusions

Adipose and muscle transcripts associated with composite phenotypes comprehensively define the repertoire of biological pathways involved in regulating distinct dimensions of obesity and glucose homeostasis. The cis-eSNPs may act as primary initiators influencing obesity and glucose homeostasis by regulating transcript levels of a subset of genes in adipose and muscle. Further computational analysis and in vitro functional studies will be required to prioritize these genes and validate the causal regulatory role of the key derivers in remodeling transcriptional regulatory networks relevant to glucose homeostasis.