Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Inter-Tissue Gene Co-Expression Networks between Metabolically Healthy and Unhealthy Obese Individuals

  • Lisette J. A. Kogelman,

    Affiliation Department of Large Animal Sciences, Faculty of Health and Medical Sciences, University of Copenhagen, Frederiksberg, Denmark

  • Jingyuan Fu,

    Affiliations University of Groningen, University Medical Center Groningen, Department of Genetics, Groningen, The Netherlands, University of Groningen, University Medical Center Groningen, Department of Paediatrics, Groningen, The Netherlands

  • Lude Franke,

    Affiliation University of Groningen, University Medical Center Groningen, Department of Genetics, Groningen, The Netherlands

  • Jan Willem Greve,

    Affiliation Department of General Surgery, Atrium Medical Center Parkstad, Heerlen, the Netherlands

  • Marten Hofker,

    Affiliation Pediatrics Molecular Genetics, University of Groningen, University Medical Center Groningen, Department of Genetics, Groningen, The Netherlands

  • Sander S. Rensen,

    Affiliation Department of Surgery, NUTRIM School of Nutrition and Translational Research in Metabolism, Maastricht University, Maastricht, The Netherlands

  • Haja N. Kadarmideen

    hajak@sund.ku.dk

    Affiliation Department of Large Animal Sciences, Faculty of Health and Medical Sciences, University of Copenhagen, Frederiksberg, Denmark

Abstract

Background

Obesity is associated with severe co-morbidities such as type 2 diabetes and nonalcoholic steatohepatitis. However, studies have shown that 10–25 percent of the severely obese individuals are metabolically healthy. To date, the identification of genetic factors underlying the metabolically healthy obese (MHO) state is limited. Systems genetics approaches have led to the identification of genes and pathways in complex diseases. Here, we have used such approaches across tissues to detect genes and pathways involved in obesity-induced disease development.

Methods

Expression data of 60 severely obese individuals was accessible, of which 28 individuals were MHO and 32 were metabolically unhealthy obese (MUO). A whole genome expression profile of four tissues was available: liver, muscle, subcutaneous adipose tissue and visceral adipose tissue. Using insulin-related genes, we used the weighted gene co-expression network analysis (WGCNA) method to build within- and inter-tissue gene networks. We identified genes that were differentially connected between MHO and MUO individuals, which were further investigated by homing in on the modules they were active in. To identify potentially causal genes, we integrated genomic and transcriptomic data using an eQTL mapping approach.

Results

Both IL-6 and IL1B were identified as highly differentially co-expressed genes across tissues between MHO and MUO individuals, showing their potential role in obesity-induced disease development. WGCNA showed that those genes were clustering together within tissues, and further analysis showed different co-expression patterns between MHO and MUO subnetworks. A potential causal role for metabolic differences under similar obesity state was detected for PTPRE, IL-6R and SLC6A5.

Conclusions

We used a novel integrative approach by integration of co-expression networks across tissues to elucidate genetic factors related to obesity-induced metabolic disease development. The identified genes and their interactions give more insight into the genetic architecture of obesity and the association with co-morbidities.

Introduction

Obesity, characterized by an excessive accumulation of adipose tissue in the body, has major consequences for human health, like type 2 diabetes (T2D) and nonalcoholic steatohepatitis (NASH). However, it has now been acknowledged that (extremely) obese individuals may also be metabolically and cardiorespiratory fit, so called metabolic healthy obese (MHO) [1, 2]. It is estimated that 10 − 25 percent of the obese individuals are MHO [1]. The excessive accumulation of adipose tissue in case of obesity disturbs the endocrine balance. Already in 1968, it was indicated that the functioning of insulin metabolism is dependent upon adipose cell size and that adaptive functioning of adipose cells is linked to the metabolic condition of individuals [3]. The expandability of the adipose tissue to be able to store large amounts of fat may be an important factor determining obesity-induced metabolic disturbances [4]. However, expandability is not an unlimited process; in fact, adipose tissue storage capacity may become saturated, resulting in excess of fat “overspilled” to non-adipose tissues and subsequent lipotoxicity which can lead to metabolic syndrome [5]. In such cases, obesity results in elevated levels of free fatty acids (FFA) affecting the pancreatic beta cells, and in the secretion of a group of adipose tissue derived cytokines, the adipokines. The direct effect of FFA is thought to be the result of activation of multiple intracellular signals in the beta cell, eventually leading to apoptosis and reduced insulin secretion [6]. To date, involved genetic factors and pathways are largely unknown.

Here, we aim to find genetic and molecular mechanisms important for human obesity-induced disease development by comparing inter-tissue gene co-expression of MHO and metabolically unhealthy obese (MUO) individuals. A study in mice showed important dynamic inter-tissue crosstalk in obesity development, with a key role for inflammatory pathways [7]. Also, the inter-tissue crosstalk of adipose tissue and skeletal tissue has shown its importance in obesity [8] and insulin resistance [9]. Weighted Gene Co-expression Network Analysis (WGCNA) [10] clusters genes based on gene-gene interactions and is used to unravel genetic mechanisms of complex diseases, including obesity [1114]. We here use this approach to create an inter-tissue network, giving the potential of studying the inter-tissue gene co-expression and hereby gaining insight into the genetic architecture across tissues and pathogenesis of obesity and comorbidities. This approach has previously led to the successful identification of key drivers of coronary artery disease [15]. In this study, we investigate the co-expression patterns of insulin-related genes in severely obese individuals, some of whom are MHO. The expression profiles of four tissues (liver, muscle, subcutaneous adipose tissue, and visceral adipose tissue) are investigated using a novel integrative inter-tissue approach. We furthermore look for potential causal genetic variants for metabolic differences under similar obesity state, by integrating, modeling, and analyzing genomic and transcriptomic data by an eQTL approach.

Materials and Methods

Genomic and transcriptomic data from the study population

Subjects.

In this study, 60 severely obese individuals (42 males, 18 females) who underwent elective bariatric surgery were included. All individuals were free from acute or chronic inflammatory or degenerative diseases, and did not consume a high amount of alcohol (>10 g/day) or anti-inflammatory drugs. This study was approved by the Medical Ethical Board of Maastricht University Medical Centre, in line with the guidelines of the 1975 Declaration of Helsinki. Informed consent in writing was obtained from each subject personally. All subjects were deeply phenotyped for a number of metabolic measures, e.g. systemic glucose, cholesterol, and C-reactive protein concentrations.

Expression profiling.

During elective bariatric surgery samples were taken from the liver, muscle (musculus rectus abdominis), subcutaneous adipose tissue (SAT; abdominal), and visceral adipose tissue (VAT; omentum majus). RNA was hybridized to Illumina HumanHT12 BeadChips and scanned on the Illumina BeadArray Reader. Raw intensity data were extracted using Illumina’s BeadStudio Gene expression module v3.2. The samples were not pooled. All raw data are publicly available at Gene Expression Omnibus (GSE22070). Raw expression data were normalized by performing a quantile normalization and log2 transformation.

Genotyping.

After eight hours of fasting on the morning of elective bariatric surgery, venous blood samples were obtained. DNA was extracted and genotyped using HumanOmniBeadChips (Illumina) and imputed using the GIANT release from the 1000 Genomes project (5,763,069 unique SNPs).

More information about the study population, expression profiling, and genotyping can be found in previous publications [1619].

Filtering of insulin-related genes for inter-tissue network construction

Insulin-related genes were filtered from the expression dataset. First, GO-terms related to insulin were found using AmiGO 2 (http://geneontology.org), by using the search term “insulin”, resulting in the detection of 92 GO-terms. We then used Biomart [20] to detect genes that have been associated with those GO-terms, resulting in 648 unique genes corresponding to 885 transcripts.

Inter-tissue behavior of individual genes

Potential biologically important inter-tissue genes are expected to show altered inter-tissue co-expressions, and were therefore identified by the difference of the inter-tissue co-expression between the MHO and MUO subnetworks. To identify those important inter-tissue genes, we calculated the correlation between tissues for each insulin-related gene in two subnetworks (MHO and MUO). This resulted in a correlation between each tissue-pair (liver-muscle, liver-SAT, liver-VAT, muscle-SAT, muscle-VAT, SAT-VAT) for each insulin-related gene in both the MHO and MUO subnetwork. Following this, between each tissue-pair the absolute difference between MHO and MUO subnetworks was calculated for each gene. Then, the sum of tissue-pair differences (inter-tissue co-expression) in both networks was calculated.

Inter-tissue co-expression network analysis

To gain insight in gene-gene interactions we constructed weighted gene co-expression networks, using the Weighted Gene Co-expression Network Analysis (WGCNA) R-package [10]. The networks were constructed across tissues for both MHO and MUO individuals by calculating Pearson’s correlation coefficients between each gene-pair across all tissues. We created the adjacency matrix with diagonal blocks representing within-tissue correlations (tissue-specific) and the off-diagonal blocks representing inter-tissue co-expression (Fig 1). The inter-tissue co-expression was computed by taking the ith gene in the tth tissue and correlated with all other Gi genes in a Tt different tissue and repeated for all genes in G and all tissues in T. Here, T is a vector with the number of tissues and G a vector with the number of genes in a given tissue. The mean and standard deviation of each block showed no significant difference between the within-tissue and inter-tissue blocks, showing no biasing difference in connection strength between within-tissue blocks and inter-tissue blocks. Modules were detected by first calculating the dissimilarity Topological Overlap Measure (TOM), which was used for clustering. Modules were defined as branches of the cluster dendrogram, which were detected using the Dynamic Tree Cut library for R [21], with a minimum module size of 50 genes.

thumbnail
Fig 1. The adjacency matrix to construct the within- and inter-tissue gene network.

The matrix presents the adjacency for L = Liver, M = Muscle, S = Subcutaneous Adipose Tissue, V = Vat, for i till n genes. In red the within-tissue blocks, in black the inter-tissue blocks.

https://doi.org/10.1371/journal.pone.0167519.g001

eQTL mapping

The eQTL mapping was performed using the eQTL-mapping-pipeline (v.1.2.4) developed by the Department of Genetics, University Medical Centre Groningen, The Netherlands, which can be found at http://github.com/molgenis/systemsgenetics/tree/master/eqtl-mapping-pipeline [22]. We used previously normalized expression data of the selected insulin-genes, as described above. The SNPs were filtered based on call rate (> 0.95), Hardy-Weinberg equilibrium (P > 1E-4), and minor allele frequency (MAF > 0.05). We performed the cis-eQTL analysis, whereby the distance between SNP and probe was set on 1 Mb on either side of the probe. Detected p-values were corrected for multiple-testing by permutation testing (n = 10), and eQTLs were considered to be significant with a FDR < 0.05. The eQTL mapping was performed with the expression data of each tissue.

Functional annotation and visualization

Genes in detected modules were retained for gene set enrichment analysis (GSEA) based on their Module Membership (MM; correlation of gene with module eigengene). GSEA was performed using HumanMine (http://humanmine.org), detecting overrepresented GO-terms and KEGG pathways. Visualization of modules was performed in Cytoscape [23].

Results and Discussion

All individuals in this study were severely obese (BMI > 35) and underwent bariatric surgery. Deep metabolic phenotyping resulted in an overview of the metabolic state of the individuals, showing that nearly half were MHO (defined as having neither T2D nor NASH). Among the MUO individuals, 18 individuals suffered from T2D (7 males, 11 females), 27 suffered from NASH (8 males, 19 females) and of them, 13 individuals had both T2D and NASH (4 males, 9 females). Descriptive statistics of the metabolic phenotypes showed a significant difference (P < 0.05) between MHO and MUO individuals for glucose, glycated hemoglobin (HbA1c), FFA, and aspartate transaminase (ASAT) (Table 1). Those metabolic phenotypes did not show a significant difference between males and females (P>0.05), though we did find a significant difference for BMI, waist-hip ratio (WH-ratio), insulin, total cholesterol and low-density lipoprotein (LDL) levels (P<0.05). A significant difference in age was found: the MHO individuals were younger than the MUO individuals, which is in agreement with a study that found a decreasing prevalence of metabolic healthy obesity with age [24].

thumbnail
Table 1. Descriptive statistics with mean and standard deviation of the individuals.

https://doi.org/10.1371/journal.pone.0167519.t001

Inter-tissue co-expression focussing on individual genes

Dobrin et al. [25] showed that investigating inter-tissue co-expression can lead to the detection of genes that are related to tissue-specific changes in diseases and that this is representing cross-tissue communication. Therefore, we identified genes that are altered in inter-tissue co-expression between MHO and MUO individuals, as they might pinpoint to genes involved in the crosstalk of tissues with respect to obesity-induced comorbidities. We detected two altered genes between the MHO and MUO subnetwork: IL1B and IL-6. Both genes showed, most of them significant, stronger inter-tissue co-expression in the MUO subnetwork than in the MHO subnetwork.

Interleukin 1-β (IL1B) is an important cytokine mainly produced by activated macrophages. IL1B showed a strong correlation between liver and adipose tissues in the MUO subnetwork, with lower correlations among the other tissues and in the MHO subnetwork (Table 2). As obesity leads to an increase of macrophage activation in liver and adipose tissue [26], IL1B mRNA levels will consequently increase. Studies have shown the importance of IL1B in the development of obesity-induced insulin resistance [27, 28], and a central role for IL1B has been suggested in macrophage-adipocyte crosstalk-related blocking of insulin action in adipose tissue [29]. Moreover, the function of IL1B in the liver-adipose tissue crosstalk has been studied in mice [30]. Recent studies showed that there was no difference in systemic IL1B levels between MHO and MUO individuals [31, 32]. Here, we do not find differences in mRNA IL1B levels, but the stronger inter-tissue correlations in the MUO subnetwork indicates a potential role for IL1B in obesity-induced metabolic diseases.

thumbnail
Table 2. Correlation of IL1B mRNA expression between liver and adipose tissues in different groups of subjects, with p-values representing the significance of across-tissue correlation.

https://doi.org/10.1371/journal.pone.0167519.t002

A similar pattern across tissues was found for Interleukin-6 (IL-6) with strong correlations between liver and adipose tissues (Table 3). IL-6 is a cytokine and a myokine, meaning it is also secreted by skeletal muscle during contraction [33]. The many functions of IL-6 include roles in immunological responses and glucose metabolism [34, 35]. IL-6 has been proposed as independent predictor of T2D [36] and the same study also showed a significant interaction between IL1B and IL-6. In contrast to IL1B, significantly decreased IL-6 mRNA levels were shown in MHO vs. MUO individuals [32]. Surprisingly, our data show no significant difference in IL-6 mRNA levels between MHO and MUO individuals.

thumbnail
Table 3. Correlation of IL-6 mRNA expression between liver and adipose tissues in different groups of subjects, with p-values representing the significance of across-tissue correlation.

https://doi.org/10.1371/journal.pone.0167519.t003

Inter-tissue network analysis with a focus on IL1B and IL-6

We constructed two large gene co-expression networks using insulin-related genes in four tissues, resulting in a network for the MHO and for the MUO individuals. Clustering within those networks resultantly can lead to tissue-specific clusters or clusters with genes across tissues. Clustering of the genes in the MHO and MUO network is visualized using a gene dendrogram (S1 Fig). We focussed on the network characteristics of the identified genes (IL1B and IL-6) as these are potentially important in disease development. IL1B and IL-6, expressed in each of the four tissues, did not cluster together across tissues or within the MHO subnetwork. However, they did cluster together in the MUO network, within each tissue. We further investigated those modules and genes, focussing on altered correlations between MHO and MUO subnetworks, to get more insight into the mechanisms involved in disease development. We confirmed that those modules were not age and/or gender dependent, as the module eigengenes of the modules were not significantly correlated with age or gender. While we focused on overall co-expression patterns across tissues in MHO ad MUO individuals, it could be expected that co-expression patterns may differ among patients suffering from either NASH or T2D, however, due to the limited sample size we were not able to home in on them.

The Greenyellow module in the MUO network, containing both the IL1B and IL-6 gene, clustered the expression profiles of 29 genes (MM > 0.6) of which 27 are coming from the liver. Two genes in this module (PKM and SLC6A5) are coming from SAT and VAT, respectively. GSEA revealed the NOD-like receptor (NLR) signalling pathway as most significant KEGG pathway (P = 0.002) and GO-terms were associated with signal transduction (e.g. negative regulation of signal transduction, P = 3.55E-5) and cytokine receptor binding (P = 6.59E-4). Those findings are in concordance with the fact that excessive amounts of adipose tissue result in an increased release of cytokines, e.g IL-6, IL1B, and CCL2, causing inflammation through, for example, the NLR signalling pathway [37].

Pearson’s correlation coefficients of genes in the Greenyellow module were visualized in Cytoscape for both sub-networks (Fig 2A). Visualization shows that genes in the Greenyellow module are strongly correlated with each other in the MUO subnetwork, while the same genes are less strongly correlated in the MHO subnetwork. As expected, strong co-expression is detected between IL1B and IL-6, with a slight increase in strength between MHO and MUO individuals (0.50 vs. 0.71, respectively). In this module, IL1B does not show large alterations in co-expression with other genes between MHO and MUO individuals. However, IL-6 shows and altered correlation with the cytokine suppressor genes SOCS2 (r = 0.01 vs. r = 0.74) and SOCS3 (r = 0.47 vs. r = 0.75), and furthermore, the correlation between SOCS2 and SOCS3 is slightly altered (r = 0.44 vs. r = 0.69). Previous studies have shown that obesity impairs JAK-STAT3 signalling as a result of elevated IL-6 levels, leading to elevated expression of for example SOCS3 in white adipose tissue, liver, and muscle [38]. By binding to insulin receptor substrates, SOCS3 impairs insulin action and elevated SOCS3 levels are associated with insulin resistance [39]. The importance of the association between SOCS3 and IL-6 has been shown before, as the absence of SOCS3 leads to altered effects of IL-6 resulting in a strong inhibiting effect on macrophages and dendritic cells, resembling an inflammatory response [40]. Also the PTPN1 gene, encoding PTP1B, is associated with insulin resistance by negatively regulating insulin signalling [41]. Due to its effects on both insulin and leptin, it has been suggested as drug target for obesity and diabetes [42]. This gene showed a highly altered co-expression with IL-6 in the module between MHO and MUO subnetworks: no correlation was found between PTP1B and IL-6 in the MHO subnetwork, while they were moderately correlated in the MUO subnetwork (r = 0.53).

thumbnail
Fig 2. Network visualization of inter-tissue network modules.

Visualization of the genes in the (A) Greenyellow module (from MUO subnetwork) in the MHO and MUO individuals, and (B) Turquoise module (from MUO subnetwork) in the MHO and MUO individuals. Genes coming from liver are coloured yellow, from muscle orange, from VAT blue and from SAT green. IL1B and IL-6 are bordered red. Edges are coloured based on their correlation on a red-green scale representing a negative-positive correlation.

https://doi.org/10.1371/journal.pone.0167519.g002

Detection of potential causal genetic variants using an eQTL mapping approach led to the detection of a noteworthy eQTL in this module: the SLC6A5 gene, expressed in VAT. SLC6A5 shows generally low correlation coefficients with other genes in the MHO subnetwork, but several strongly correlated genes in the MUO subnetwork. Correlations with IL1B and IL-6 are only moderate in MUO (0.28 and 0.41, respectively) and low in MHO (-0.13 and -0.05, respectively) subnetworks. The strongest correlations of SLC6A5 within the MUO subnetwork are found with ANXA1 (r = 0.66), CCL2 (r = 0.63), MARCKS (r = 0.61), and SOCS3 (r = 0.61); all of which are not correlated in the MHO subnetwork. SLC6A5 is encoding the neuronal glycine transporter 2 (GlyT2) and its action involves protein kinase C (PKC) pathways [35]. Even though SLC6A5 itself has not been associated with obesity, its regulating properties via PKC pathways might affect the metabolic state of obese individuals. For example, ANXA1 (associated with adiposity [43]) is an important mediator in the neuroendocrine system and PKC-dependent mechanisms are essential for its activity [44]. Likewise, the previously discussed SOCS3 has been linked to involvement in PKC pathways [45]. We therefore suggest a potential causal role for SLC6A5 in regulatory processes related to obesity-induced metabolic diseases. This is supported by the fact that a mutation in this gene causes, among others, a decreased body weight in mice [46].

The Turquoise module in the MUO subnetwork (Fig 2B) clustered the expression profiles of 26 unique genes (MM > 0.9), all coming from the muscle. Both IL1B and IL-6 were filtered out due to the strict threshold on the MM, but were included in comparison of correlation strengths between the MHO and MUO subnetworks. The correlations among them were altered between the MHO and MUO subnetwork (r = 0.25 vs r = 0.75, resp.). One well-known obesity-genes is Leptin (LEP), encoding leptin, which is also called the “satiety hormone”. In skeletal muscle, it promotes energy dissipation and prevents fatty acid accumulation [47]. LEP shows strong positive and negative correlations with all other Turquoise module-genes (absolute correlation > 0.8) in the MUO subnetwork, while being low with many of the genes in the MHO subnetwork (absolute correlation < 0.4). The largest differential correlation for LEP was found between LEP and PTPRE (0.22 vs. 0.91 in MHO vs. MUO). PTPRE (Protein Tyrosine Phosphatase, Receptor Type, E) negatively regulates insulin signaling in skeletal muscle and has been suggested to play a role as negative feedback regulator of leptin signaling via JAK2 [8]. Furthermore, strong alterations were found for the correlation between LEP and IL1B (r = -0.15 vs. r = 0.81 in MHO vs. MUO). It has been shown that IL1B is necessary for induction of leptin during inflammation [12]. Our data suggest that this induction is not occurring in MHO individuals. Another interesting gene is Adiponectin (ADIPOQ), encoding the hormone adiponectin that enhances skeletal muscle insulin sensitivity and has been suggested as a drug target for obesity and T2D [13]. The correlation between ADIPOQ and IL1B differed to a large extent between MHO and MUO (r = -0.19 vs r = 0.68 in MHO vs. MUO). In adipocytes, but not in skeletal muscle, it has been shown that IL1B reduces adiponectin production, thereby negatively affecting insulin sensitivity [24, 29]. No eQTLs were detected in this module.

The Black module in the MUO subnetwork (Fig 3A) clustered the expression profiles of 30 genes (MM < 0.60), all coming from SAT, except for two that were coming from VAT: PRKAG2 (also co-expressed in SAT) and PIK3R1. The correlation between IL1B and IL-6, and their correlation with cytokine suppressor genes, was similar between MHO and MUO subnetworks in SAT. The two genes that are expressed in VAT, PRKAG2 and PIK3R1, are showing altered co-expression with IL-6 and IL1B. The co-expression of PIK3R1 with IL1B and IL-6 in the MHO subnetwork is very low (-0.16 and -0.20, respectively), but strongly negative in the MUO subnetwork (both -0.67). PIK3R1 encodes part of the phosphatidylinositol 3-kinase (PI3K) enzyme, and PI3K signaling plays an important role in insulin signaling. PIK3R1 is involved in the mediation of insulin sensitivity and the inflammatory response in adipose tissue [31]. The altered co-expression between the insulin genes and PIK3R1 indicates an important role in disease development that has been identified by others as well [48]. The correlation of PRKAG2 with IL1B and IL-6 is around zero in the MHO, but moderately positive in the MUO subnetwork (0.53 and 0.54, respectively). PRKAG2 is part of the AMP-activated protein kinase (AMPK) pathway, which is important in energy regulation, e.g. glucose homeostasis and insulin sensitivity [49]. Due to its properties, it is a drug target for metabolic syndrome.

thumbnail
Fig 3. Network visualization of the inter-tissue network modules.

Visualization of the genes in the (A) Black module (from MUO subnetwork) in the MHO and MUO individuals, and (B) Salmon module (from MUO subnetwork) in the MHO and MUO individuals. Genes coming from liver are coloured yellow, from muscle orange, from VAT blue and from SAT green. IL1B and IL-6 are bordered red. Edges are coloured based on their correlation on a red-green scale representing a negative-positive correlation.

https://doi.org/10.1371/journal.pone.0167519.g003

In this module we found one eQTL: the IL-6 receptor (IL-6R). Many of the correlations between IL-6R and other genes were not altered between MHO and MUO subnetworks. However, the correlation with SOCS3 was elevated (0.29 vs 0.67, respectively) and with PIK3R1 (in VAT) was decreased (0.01 vs -0.58, respectively). Surprisingly, the correlations of IL-6 with those genes were unaltered, implicating an important role for the receptor of IL-6, also based on previous discussion of results for the JAK-STAT pathway, in obesity-induced disease development. Several altered phenotypes have been detected in mutations in the IL-6 receptor in mice, e.g. insulin resistance and inflated inflammatory action [50]. However, in a large human survey, a significant correlation with BMI and weight loss within morbidly obese individuals could only be found with IL-6R expression levels in liver, and not in omental adipose, subcutaneous adipose and stomach tissue [51].

The Salmon Module in the MUO subnetwork (Fig 3B) clustered the expression profiles of 33 genes (MM < 0.6). These genes were all coming from VAT, except for one gene (Syntaxin 4, STX4) that was coming from both SAT and VAT. Several genes from previously discussed modules were also present in this module. The co-expression of IL1B and IL-6 with SOCS2 and SOCS3 is similar to the co-expression found in SAT (Black module), meaning that their correlation is only altered in the liver. Also the co-expression between IL1B and IL-6 is unaltered in SAT. The co-expression between STX4 in SAT and VAT is stronger in the MUO than in the MHO subnetwork (0.28 vs 0.58, respectively). The co-expression of STX4 in SAT is altered with IL1B, the IL1B receptor antagonist (IL1RN), IL6, and the IL-6 receptor (IL6R). In all four cases, the co-expression is moderate to strongly positive in the MUO subnetwork and low negative co-expressed in the MHO subnetwork. However, the co-expression of STX4 in VAT with the four genes that were expressed in VAT is similar between MHO and MUO subnetworks, suggesting a role in disease development for STX4 only in SAT.

In this module we found one eQTL: PTPRE, expressed in VAT. PTPRE is a negative regulator of insulin signalling in muscle [47]. In a study with obese individuals, it was shown that PTPRE was differentially expressed and methylated before and after bariatric surgery, and that NASH phenotype was negatively correlated with the bariatric reconstruction [52]. Besides the detection of PTPRE as eQTL, the PTPRE gene is present in all four modules that were further investigated due to presence of IL1B and IL-6. Homing in on the cross-tissue talk of this gene, the main change is found between the correlation in liver and VAT (0.05 vs 0.55 in MHO vs. MUO). In the Salmon module PTPRE shows altered co-expression with the previously discussed ISL1 gene (0.07 vs -0.60 in MHO vs MUO), encoding the insulin gene enhancer protein ISL-1. Many of the correlations of ISL1 with other genes in this module are altered, with strong negative correlations in the MUO subnetwork. Previously, it has been shown that ISL1 is expressed in VAT and negatively correlated with BMI and abdominal fat [34]. They also showed that expression was reduced in obese mice but reduced in lean insulin sensitive mice. In the Greenyellow module, with mostly liver genes, alterations between MHO and MUO individuals are found within the correlations of the expression of PTPRE in liver with e.g. SLC6A5 (VAT), SLC6A13 (liver) and PKM (SAT). It is remarkable that the module consists of mostly genes co-expressed in the liver, but that PTRPE shows altered co-expression with the two genes that are coming from SAT and VAT, suggesting an important function for PTPRE as a link between tissues in disease development resulting from obesity. To date, no studies have linked PTPRE with any of those genes.

In summary, we have shown the integration of gene co-expression networks across tissues in MHO and MUO individuals to identify genes and pathways related to obesity-induced disease development. The results provide important insights into genomics of adipose tissue expandability, lipotoxicity, and eventual comorbidities such as T2D and NASH. In a co-expression network setting, we detected IL-6 and IL1B as key genes for inter-tissue gene co-expression differences related to metabolic state. By investigating the modules in which IL-6 and IL1B were co-expressed, we detected many altered co-expressed genes and pathways that might be important in obesity-induced inflammation and comorbidity development. Even though IL-6 and IL1B mRNA levels were not altered between MHO and MUO individuals in our study, their co-expression with other genes indicates a potentially important role in obesity-induced development of metabolic disturbances. The chosen approach gave us insight into co-expression between genes in a network setting, but could not give information about directions in the network. However, an eQTL mapping approach was chosen to detect genes that affect mRNA levels and thereby affecting health status. This led to the identification of several genes (PTPRE, IL-6R, and SLC6A5) which might have an important role in insulin-related pathways of obese individuals. Functional validation of those genes is needed to identify their potential role in the development of obesity-induced comorbidities.

Supporting Information

S1 Fig. Clustering of the insulin-related genes.

Clustering is based on the dissimilarity Topological Overlap Measure (TOM) within the A) metabolically healthy obese (MHO) and B) metabolically unhealthy obese (MUO) network. Modules are detected using the Dynamic Tree Cut algorithm and presented by the color-coded bar under the dendrogram.

https://doi.org/10.1371/journal.pone.0167519.s001

(PNG)

Acknowledgments

Lisette J.A. Kogelman was funded by a post-doctoral fellowship within Biochild project (http://biochild.ku.dk/) with grant (0603-00457B) from the Programme Commission on Individuals, Disease and Society under the Danish Council for Strategic Research (Innovationsfonden). Haja N. Kadarmideen thanks EU-FP7 Marie Curie Actions—Career Integration Grant (CIG-293511) for partially funding his time spent on this research. The authors thank Dr. Tom Michoel for fruitful discussions and input about the methodology regarding across-tissue network construction. Prof. Wim Buurman is thanked for his role in setting up the cohort of severely obese subjects.

Author Contributions

  1. Conceptualization: LJAK HNK.
  2. Data curation: JF LF JWG MH SSR.
  3. Formal analysis: LJAK HNK.
  4. Funding acquisition: HNK.
  5. Methodology: LJAK HNK.
  6. Project administration: JF LF JWG MH SSR.
  7. Supervision: JF LF JWG MH SSR HNK.
  8. Visualization: LJAK.
  9. Writing – original draft: LJAK.
  10. Writing – review & editing: LJAK JF LF JWG MH SSR HNK.

References

  1. 1. Bluher M. The distinction of metabolically 'healthy' from 'unhealthy' obese individuals. Current opinion in lipidology. 2010;21:38–43. pmid:19915462
  2. 2. Lavie CJ, De Schutter A, Milani RV. Healthy obese versus unhealthy lean: the obesity paradox. Nat Rev Endocrinol. 2015;11(1):55–62. pmid:25265977
  3. 3. Salans LB, Knittle JL, Hirsch J. The role of adipose cell size and adipose tissue insulin sensitivity in the carbohydrate intolerance of human obesity. J Clin Invest. 1968;47(1):153–65. pmid:16695937
  4. 4. Medina-Gomez G, Virtue S, Lelliott C, Boiani R, Campbell M, Christodoulides C, et al. The link between nutritional status and insulin sensitivity is dependent on the adipocyte-specific peroxisome proliferator-activated receptor-gamma2 isoform. Diabetes. 2005;54(6):1706–16. pmid:15919792
  5. 5. Slawik M, Vidal-Puig AJ. Adipose tissue expandability and the metabolic syndrome. Genes & Nutrition. 2007;2(1):41–5.
  6. 6. Eldor R, Raz I. Lipotoxicity versus adipotoxicity—The deleterious effects of adipose tissue on beta cells in the pathogenesis of type 2 diabetes. Diabetes Research and Clinical Practice. 2006;74(2, Supplement):S3–S8.
  7. 7. Samdani P, Singhal M, Sinha N, Tripathi P, Sharma S, Tikoo K, et al. A Comprehensive Inter-Tissue Crosstalk Analysis Underlying Progression and Control of Obesity and Diabetes. Scientific Reports. 2015;5:12340. pmid:26202695
  8. 8. Rousso-Noori L, Knobler H, Levy-Apter E, Kuperman Y, Neufeld-Cohen A, Keshet Y, et al. Protein tyrosine phosphatase epsilon affects body weight by downregulating leptin signaling in a phosphorylation-dependent manner. Cell Metab. 2011;13(5):562–72. pmid:21531338
  9. 9. Argiles JM, Lopez-Soriano J, Almendro V, Busquets S, Lopez-Soriano FJ. Cross-talk between skeletal muscle and adipose tissue: a link with obesity? Medicinal research reviews. 2005;25(1):49–65. pmid:15389734
  10. 10. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.
  11. 11. Min JL, Nicholson G, Halgrimsdottir I, Almstrup K, Petri A, Barrett A, et al. Coexpression Network Analysis in Abdominal and Gluteal Adipose Tissue Reveals Regulatory Genetic Loci for Metabolic Syndrome and Related Phenotypes. PLoS Genetics. 2012;8(2):e1002505. pmid:22383892
  12. 12. Dewey FE, Perez MV, Wheeler MT, Watt C, Spin J, Langfelder P, et al. Gene coexpression network topology of cardiac development, hypertrophy, and failure. Circulation Cardiovascular genetics. 2011;4(1):26–35. pmid:21127201
  13. 13. de Jong S, Boks MP, Fuller TF, Strengman E, Janson E, de Kovel CG, et al. A gene co-expression network in whole blood of schizophrenia patients is independent of antipsychotic-use and enriched for brain-expressed genes. PLoS One. 2012;7(6):e39498. pmid:22761806
  14. 14. Kogelman LJA, Cirera S, Zhernakova D, Fredholm M, Franke L, Kadarmideen H. Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA Sequencing in a porcine model. BMC Med Genomics. 2014;7(1):57.
  15. 15. Talukdar Husain A, Foroughi Asl H, Jain Rajeev K, Ermel R, Ruusalepp A, Franzén O, et al. Cross-Tissue Regulatory Gene Networks in Coronary Artery Disease. Cell Systems. 2016;2:196–208. pmid:27135365
  16. 16. Fu J, Wolfs MGM, Deelen P, Westra H-J, Fehrmann RSN, te Meerman GJ, et al. Unraveling the Regulatory Mechanisms Underlying Tissue-Dependent Genetic Variation of Gene Expression. PLoS Genet. 2012;8(1):e1002431. pmid:22275870
  17. 17. Fehrmann RS, Jansen RC, Veldink JH, Westra HJ, Arends D, Bonder MJ, et al. Trans-eQTLs reveal that independent genetic variants associated with a complex phenotype converge on intermediate genes, with a major role for the HLA. PLoS Genet. 2011;7(8):e1002197. pmid:21829388
  18. 18. Wolfs M, Rensen S, Bruin-Van Dijk E, Verdam F, Greve J-W, Sanjabi B, et al. Co-expressed immune and metabolic genes in visceral and subcutaneous adipose tissue from severely obese individuals are associated with plasma HDL and glucose levels: a microarray study. BMC Med Genomics. 2010;3(1):34.
  19. 19. Bonder MJ, Kasela S, Kals M, Tamm R, Lokk K, Barragan I, et al. Genetic and epigenetic regulation of gene expression in fetal and adult human livers. BMC Genomics. 2014;15:860. pmid:25282492
  20. 20. Haider S, Ballester B, Smedley D, Zhang J, Rice P, Kasprzyk A. BioMart Central Portal—unified access to biological data. Nucleic Acids Research. 2009;37(suppl 2):W23–W7.
  21. 21. Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics. 2008;24(5):719–20. pmid:18024473
  22. 22. Westra H-J, Peters MJ, Esko T, Yaghootkar H, Schurmann C, Kettunen J, et al. Systematic identification of trans eQTLs as putative drivers of known disease associations. Nature Genetics. 2013;45(10):1238–43. pmid:24013639
  23. 23. Shannon P, Markiel A, Ozier O, Baliga N, Wang J, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504. pmid:14597658
  24. 24. van Vliet-Ostaptchouk JV, Nuotio ML, Slagter SN, Doiron D, Fischer K, Foco L, et al. The prevalence of metabolic syndrome and metabolically healthy obesity in Europe: a collaborative analysis of ten large cohort studies. BMC endocrine disorders. 2014;14:9. pmid:24484869
  25. 25. Dobrin R, Zhu J, Molony C, Argman C, Parrish ML, Carlson S, et al. Multi-tissue coexpression networks reveal unexpected subnetworks associated with disease. Genome Biology. 2009;10(R55).
  26. 26. Weisberg S, McCann D, Desai M, Rosenbaum M, Leibel R, Ferrante A. Obesity is associated with macrophage accumulation in adipose tissue. J Clin Invest. 2003;112:1796–808. pmid:14679176
  27. 27. Jager J, Gremeaux T, Cormont M, Le Marchand-Brustel Y, Tanti JF. Interleukin-1beta-induced insulin resistance in adipocytes through down-regulation of insulin receptor substrate-1 expression. Endocrinology. 2007;148(1):241–51. pmid:17038556
  28. 28. Mojtaba E, Mahdi K. Serum interleukin-1 beta plays an important role in insulin secretion in type II diabetic. International Journal of Biosciences. 2011;1(3):93–9.
  29. 29. Bing C. Is interleukin-1β a culprit in macrophage-adipocyte crosstalk in obesity? Adipocyte. 2015;4(2):149–52. pmid:26167419
  30. 30. Nov O, Shapiro H, Ovadia H, Tarnovscki T, Dvir I, Shemesh E, et al. Interleukin-1β Regulates Fat-Liver Crosstalk in Obesity by Auto-Paracrine Modulation of Adipose Tissue Inflammation and Expandability. PLoS ONE. 2013;8(1):e53626. pmid:23341960
  31. 31. Ahl S, Guenther M, Zhao S, James R, Marks J, Szabo A, et al. Adiponectin Levels Differentiate Metabolically Healthy vs Unhealthy Among Obese and Nonobese White Individuals. The Journal of clinical endocrinology and metabolism. 2015;100(11):4172–80. pmid:26401592
  32. 32. Marques-Vidal P, Velho S, Waterworth D, Waeber G, von Kanel R, Vollenweider P. The association between inflammatory biomarkers and metabolically healthy obesity depends of the definition used. European journal of clinical nutrition. 2012;66(4):426–35. pmid:21952696
  33. 33. Keller P, Keller C, Carey AL, Jauffred S, Fischer CP, Steensberg A, et al. Interleukin-6 production by contracting human skeletal muscle: autocrine regulation by IL-6. Biochem Biophys Res Commun. 2003;310(2):550–4. pmid:14521945
  34. 34. Schobitz B, Pezeshki G, Pohl T, Hemmann U, Heinrich PC, Holsboer F, et al. Soluble interleukin-6 (IL-6) receptor augments central effects of IL-6 in vivo. FASEB journal: official publication of the Federation of American Societies for Experimental Biology. 1995;9(8):659–64.
  35. 35. Fornes A, Nunez E, Alonso-Torres P, Aragon C, Lopez-Corcuera B. Trafficking properties and activity regulation of the neuronal glycine transporter GLYT2 by protein kinase C. Biochem J. 2008;412(3):495–506. pmid:18341477
  36. 36. Spranger J, Kroke A, Möhlig M, Hoffmann K, Bergmann MM, Ristow M, et al. Inflammatory Cytokines and the Risk to Develop Type 2 Diabetes: Results of the Prospective Population-Based European Prospective Investigation into Cancer and Nutrition (EPIC)-Potsdam Study. Diabetes. 2003;52(3):812–7. pmid:12606524
  37. 37. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57. pmid:19131956
  38. 38. Wunderlich CM, Hovelmeyer N, Wunderlich FT. Mechanisms of chronic JAK-STAT3-SOCS3 signaling in obesity. Jak-stat. 2013;2(2):e23878. pmid:24058813
  39. 39. Olson AL. Insulin resistance: cross-talk between adipose tissue and skeletal muscle, through free fatty acids, liver X receptor, and peroxisome proliferator-activated receptor-alpha signaling. Hormone molecular biology and clinical investigation. 2013;15(3):115–21. pmid:25436738
  40. 40. Hunter CA, Jones SA. IL-6 as a keystone cytokine in health and disease. Nat Immunol. 2015;16(5):448–57. pmid:25898198
  41. 41. Bluher M. The distinction of metabolically 'healthy' from 'unhealthy' obese individuals. Current opinion in lipidology. 2010;21(1):38–43. pmid:19915462
  42. 42. Min JL, Nicholson G, Halgrimsdottir I, Almstrup K, Petri A, Barrett A, et al. Coexpression Network Analysis in Abdominal and Gluteal Adipose Tissue Reveals Regulatory Genetic Loci for Metabolic Syndrome and Related Phenotypes. PLoS Genet. 2012;8(2):e1002505. pmid:22383892
  43. 43. Kosicka A, Cunliffe AD, Mackenzie R, Zariwala MG, Perretti M, Flower RJ, et al. Attenuation of plasma annexin A1 in human obesity. FASEB journal: official publication of the Federation of American Societies for Experimental Biology. 2013;27(1):368–78.
  44. 44. John C, Cover P, Solito E, Morris J, Christian H, Flower R, et al. Annexin 1-dependent actions of glucocorticoids in the anterior pituitary gland: roles of the N-terminal domain and protein kinase C. Endocrinology. 2002;143(8):3060–70. pmid:12130572
  45. 45. Wallerstedt E, Smith U, Andersson CX. Protein kinase C-delta is involved in the inflammatory effect of IL-6 in mouse adipose cells. Diabetologia. 2010;53(5):946–54. pmid:20151299
  46. 46. Bogdanik LP, Chapman HD, Miers KE, Serreze DV, Burgess RW. A MusD Retrotransposon Insertion in the Mouse Slc6a5 Gene Causes Alterations in Neuromuscular Junction Maturation and Behavioral Phenotypes. PLoS ONE. 2012;7(1):e30217. pmid:22272310
  47. 47. Aga-Mizrachi S, Brutman-Barazani T, Jacob AI, Bak A, Elson A, Sampson SR. Cytosolic protein tyrosine phosphatase-epsilon is a negative regulator of insulin signaling in skeletal muscle. Endocrinology. 2008;149(2):605–14. pmid:18006633
  48. 48. Beretta M, Bauer M, Hirsch E. PI3K signaling in the pathogenesis of obesity: The cause and the cure. Advances in biological regulation. 2015;58:1–15. pmid:25512233
  49. 49. Schultze SM, Hemmings BA, Niessen M, Tschopp O. PI3K/AKT, MAPK and AMPK signalling: protein kinases in glucose homeostasis. Expert reviews in molecular medicine. 2012;14:e1. pmid:22233681
  50. 50. Wunderlich FT, Ströhle P, Könner AC, Gruber S, Tovar S, Brönneke HS, et al. Interleukin-6 Signaling in Liver-Parenchymal Cells Suppresses Hepatic Inflammation and Improves Systemic Insulin Action. Cell Metabolism. 2010;12(3):237–49. pmid:20816090
  51. 51. Greenawalt DM, Dobrin R, Chudin E, Hatoum IJ, Suver C, Beaulaurier J, et al. A survey of the genetics of stomach, liver, and adipose gene expression from a morbidly obese cohort. Genome Research. 2011;21(7):1008–16. pmid:21602305
  52. 52. Ahrens M, Ammerpohl O, von Schönfels W, Kolarova J, Bens S, Itzel T, et al. DNA Methylation Analysis in Nonalcoholic Fatty Liver Disease Suggests Distinct Disease-Specific and Remodeling Signatures after Bariatric Surgery. Cell Metabolism. 2013;18(2):296–302. pmid:23931760