# Bariatric Surgery Induces Disruption in Inflammatory Signaling Pathways Mediated by Immune Cells in Adipose Tissue: A RNA-Seq Study

Christine Poitou, Claire Perret, François Mathieu, Vinh Truong, Yuna Blum, Hervé Durand, Rohia Alili, Nadjim Chelghoum, Véronique Pelloux, Judith Aron-Wisnewsky, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Christine Poitou, Claire Perret, François Mathieu, Vinh Truong, Yuna Blum, et al., Bariatric Surgery Induces Disruption in Inflammatory Signaling Pathways Mediated by Immune Cells in Adipose Tissue: A RNA-Seq Study. PLoS ONE, 2015, 10 (5), pp.e0125718. 10.1371/journal.pone.0125718. hal-01227880

HAL Id: hal-01227880

https://hal.sorbonne-universite.fr/hal-01227880

Submitted on 12 Nov 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# Bariatric Surgery Induces Disruption in Inflammatory Signaling Pathways Mediated by Immune Cells in Adipose Tissue: A RNA-Seq Study

Christine Poitou¹,²*, Claire Perret¹,³*, François Mathieu¹,³, Vinh Truong¹,³, Yuna Blum⁴, Hervé Durand¹,³, Rohia Alili¹,², Nadjim Chelghoum⁵, Véronique Pelloux¹,², Judith Aron-Wisnewsky¹,², Adriana Torcivia⁶, Jean-Luc Bouillot⁷, Brian W. Parks⁴, Ewa Ninio¹,³, Karine Clément¹,², Laurence Tiret¹,³*

¹ Institute of Cardiometabolism And Nutrition (ICAN), Assistance Publique-Hôpitaux de Paris, Pitié-Salpêtrière Hospital, Nutrition Department, F-75013, Paris, France, ² Sorbonne Universités, University Pierre et Marie Curie (UPMC), Institut National de la Santé et de la Recherche Médicale (INSERM) UMR_S 1166, Nutriomics team, F-75005, Paris, France, ³ Sorbonne Universités, University Pierre et Marie Curie (UPMC), Institut National de la Santé et de la Recherche Médicale (INSERM) UMR_S 1166, Genomics and Pathophysiology of Cardiovascular Diseases team, F-75013, Paris, France, ⁴ Department of Medicine/Division of Cardiology, University of California Los Angeles, Los Angeles, California 90095, United States of America, ⁵ Sorbonne Universités, University Pierre et Marie Curie (UPMC), Post-Genomic Platform of Pitié-Salpêtrière (P3S), F-75013, Paris, France, ⁶ Assistance Publique-Hôpitaux de Paris, Department of Visceral Surgery, Pitié-Salpêtrière Hospital, F-75013, Paris, France, ⁷ Assistance Publique-Hôpitaux de Paris, Department of General, Digestive and Metabolic Surgery, Ambroise-Paré Hospital, F-92100, Boulogne-Billancourt, France

* These authors contributed equally to this work.

* laurence.tiret@upmc.fr

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Poitou C, Perret C, Mathieu F, Truong V, Blum Y, Durand H, et al. (2015) Bariatric Surgery Induces Disruption in Inflammatory Signaling Pathways Mediated by Immune Cells in Adipose Tissue: A RNA-Seq Study. PLoS ONE 10(5): e0125718. doi:10.1371/journal.pone.0125718

**Academic Editor:** Susanne Kaser, Medical University Innsbruck, AUSTRIA

**Received:** November 14, 2014

**Accepted:** March 17, 2015

**Published:** May 4, 2015

**Copyright:** © 2015 Poitou et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All expression files are available from the GEO database (accession number GSE65540).

**Funding:** The research was supported by the Fondation pour la Recherche Médicale (FRM ING20121226269) (F.M.), the Agence Nationale de la Recherche, France (ANR 09 GENO 106 01) (V.T.), the Assistance Publique-Hôpitaux de Paris (APHP) and the Direction of Clinical Research (PHRC MICROBARIA). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Abstract

### Background

Bariatric surgery is associated to improvements in obesity-associated comorbidities thought to be mediated by a decrease of adipose inflammation. However, the molecular mechanisms behind these beneficial effects are poorly understood.

### Methodology/Principal Findings

We analyzed RNA-seq expression profiles in adipose tissue from 22 obese women before and 3 months after surgery. Of 15,972 detected genes, 1214 were differentially expressed after surgery at a 5% false discovery rate. Upregulated genes were mostly involved in the basal cellular machinery. Downregulated genes were enriched in metabolic functions of adipose tissue. At baseline, 26 modules of coexpressed genes were identified. The four most stable modules reflected the innate and adaptive immune responses of adipose tissue. A first module reflecting a non-specific signature of innate immune cells, mainly macrophages, was highly conserved after surgery with the exception of *DUSP2* and *CD300C*. A second module reflected the adaptive immune response elicited by T lymphocytes; after surgery, a disconnection was observed between genes involved in T-cell signaling and mediators of
**Competing Interests:** The authors have declared that no competing interests exist.

the signal transduction such as *CXCR1*, *CXCR2*, *GPR97*, *CCR7* and *IL7R*. A third module reflected neutrophil-mediated inflammation; after surgery, several genes were dissociated from the module, including *S100A8*, *S100A12*, *CD300E*, *VNN2*, *TUBB1* and *FAM65B*. We also identified a dense network of 19 genes involved in the interferon-signaling pathway which was strongly preserved after surgery, with the exception of *DDX60*, an antiviral factor involved in RIG-I-mediated interferon signaling. A similar loss of connection was observed in lean mice compared to their obese counterparts.

## Conclusions/Significance

These results suggest that improvements of the inflammatory state following surgery might be explained by a disruption of immuno-inflammatory cascades involving a few crucial molecules which could serve as potential therapeutic targets.

## Introduction

The widespread increase of obesity in the past few decades has become a major global health challenge [1]. A main burden of obesity lies in its interconnection with a number of metabolic and cardiovascular diseases including type 2 diabetes mellitus (T2D), dyslipidemia, hypertension and atherosclerosis.

Adipose tissue is increasingly recognized as a highly active organ involved in numerous metabolic, hormonal and immune processes, which plays a pivotal role in the development of obesity-related complications. A central mechanism linking obesity to metabolic dysfunctions is adipose inflammation which is caused by infiltration of immune cells into fat depots and increased production and secretion of pro- and anti-inflammatory factors into the circulation [2]. Besides macrophages which are the most abundant immune cells in fat depots, nearly all immune cell populations are modulated by the metabolic state of adipose tissue [3].

Bariatric surgery is the most effective therapy of severe forms of human obesity. It is associated with improvements in co-morbidities, including T2D, insulin resistance and cardiovascular risks [4, 5]. Several studies also reported changes in the inflammatory state both in adipose tissue and blood, particularly at early stage during the first year post-surgery [6]. However, the pathophysiological mechanisms behind these beneficial effects are poorly understood. Identifying the underlying molecular mechanisms that are critical to metabolic dysregulation and may contribute to surgery-induced improvements would be of considerable interest for targeting affected pathways.

Gene expression profiling has been instrumental in uncovering several biological pathways affected by weight loss. Transcriptomic studies of adipose tissue have highlighted major modifications induced by weight loss in lipid metabolism, extracellular matrix (ECM) organization and inflammatory processes [7–11]. However, all previous studies were based on the microarray technology which has been shown to be affected by several drawbacks. Recently the RNA-sequencing (RNA-seq) technology has emerged as a powerful and reliable tool for profiling gene expression which overcomes most limitations of microarrays [12].

In this study, we analyzed RNA-seq expression profiles in subcutaneous adipose tissue (SAT) from 22 obese women undergoing bariatric surgery. Expression levels were assessed before and 3 months after surgery. We investigated modifications in SAT transcriptome by applying both single-gene and network approaches. This combination of approaches provided
insights not only into major changes in transcriptome profiles, but also into underlying molecular mechanisms causing disruption of biological pathways. Our study suggested that improvement of the inflammatory state after bariatric surgery might be explained by a break in the coexpression between inflammatory signaling pathways and a few crucial molecules involved in chemotaxis and activation of immune cells.

## Results

The mean age of the 22 women participating in the study was 37.5 ± 8.3 years. Bariatric surgery resulted in a decrease of the mean body mass index (BMI) from 46.5 to 39.9 kg/m² and a significant improvement of all biological and metabolic parameters as expected (Table 1).

The mean number of read pairs produced by sequencing was 40.8 million pairs per sample (min: 35.5; max: 51.3). After mapping to the human reference transcriptome and filtering out genes not or weakly expressed (<10 counts in average at both times), 15,972 of the 23,329 annotated genes were included in analysis. In average, these 15,972 genes totalized 28.5 million read counts per sample (min: 23.8; max: 39.3).

## Effect of bariatric surgery on single gene expression levels

At a 5% false discovery rate (FDR), 1214 genes (7.6%) were differentially expressed (DE) three months after surgery (S1 Table). In the top list of upregulated expressions, there were several genes involved in mRNA processing and translation initiation such as heterogeneous nuclear ribonucleoproteins (HNRNP symbol), ribosomal proteins (RP symbol) and eukaryotic translation initiation factors (EIF and EEF symbols). The most significantly downregulated gene was *STMN2* which encodes stathmin-like 2, a protein primarily known for its role in dynamics of

<table><caption>Table 1. Descriptive characteristics of the study population (n = 22 women).</caption><thead><tr><th></th><th>Before surgery</th><th>3 months after surgery</th><th>P-value</th></tr></thead><tbody><tr><td>Age (years)</td><td>37.5 ± 8.3</td><td></td><td></td></tr><tr><td>BMI (kg/m²)</td><td>46.5 ± 5.6</td><td>39.9 ± 4.9</td><td>&lt;10<sup>-4</sup></td></tr><tr><td>Type 2 Diabetes</td><td>3 (13.6%)</td><td>2 (9.1%)</td><td>NS</td></tr><tr><td>Metabolic syndrome</td><td>11 (50.0%)</td><td>8 (36.4%)</td><td>NS</td></tr><tr><td>Total cholesterol (mmol/L)</td><td>4.83 ± 0.70</td><td>4.44 ± 0.82</td><td>0.007</td></tr><tr><td>LDL-cholesterol (mmol/L)</td><td>3.14 ± 0.58</td><td>2.90 ± 0.70</td><td>0.07</td></tr><tr><td>HDL-cholesterol (mmol/L)</td><td>1.11 ± 0.32</td><td>1.03 ± 0.30</td><td>0.07</td></tr><tr><td>Triglycerides (mmol/L)</td><td>1.30 ± 0.59</td><td>1.13 ± 0.43</td><td>0.08</td></tr><tr><td>Glucose (mmol/L)</td><td>5.12 ± 0.65</td><td>4.81 ± 1.15</td><td>0.25</td></tr><tr><td>Insulin (µUI/ml)</td><td>21.6 ± 10.9</td><td>12.2 ± 6.73</td><td>2x10<sup>-4</sup></td></tr><tr><td>CRP (mg/L)</td><td>9.36 ± 4.34</td><td>4.81 ± 2.43</td><td>0.005</td></tr><tr><td>IL-6 (pg/mL)</td><td>5.40 ± 3.11</td><td>4.06 ± 1.69</td><td>0.09</td></tr><tr><td>Leptin (ng/mL)</td><td>85.4 ± 35.5</td><td>39.6 ± 20.8</td><td>&lt;10<sup>-4</sup></td></tr><tr><td>Adiponectin (µg/mL)</td><td>4.59 ± 1.94</td><td>5.56 ± 3.04</td><td>0.051</td></tr><tr><td>Fat mass (kg)</td><td>61.9 ± 10.5</td><td>50.9 ± 8.3</td><td>&lt;10<sup>-4</sup></td></tr><tr><td>Lean mass (kg)</td><td>59.4 ± 6.4</td><td>53.5 ± 5.2</td><td>&lt;10<sup>-4</sup></td></tr><tr><td>Adipocyte volume (pL)</td><td>936.7 ± 188.4</td><td>823.1 ± 168.9</td><td>0.007</td></tr><tr><td>White blood cells (10<sup>3</sup>/mm<sup>3</sup>)</td><td>8.34 ± 2.26</td><td>6.38 ± 1.75</td><td>&lt;10<sup>-4</sup></td></tr></tbody><tfoot><tr><td colspan="4">Data are mean ± SD or n (%)</td></tr></tfoot></table>

doi:10.1371/journal.pone.0125718.t001
microtubules during neuronal growth [13] but also thought to be involved in the regulation of
the fate of adipocyte differentiation [14].

Expression of the *LEP* gene encoding leptin, a key adipose-derived hormone whose serum level positively correlates to fat mass, was reduced by 31% (FDR = 0.001). *LOX*, a gene encoding lysyl-oxidase involved in collagen cross-linking was also downregulated (FDR = 5.7 10⁻⁴). By contrast, expression of the *ADIPOQ* gene encoding adiponectin, another adipokine inversely related to obesity-associated insulin resistance, did not significantly vary (FDR = 0.37). Expression of *FTO*, the most widely replicated gene in genome-wide association study of obesity, was not affected by surgery (FDR = 0.99). Several gene expressions that were reported in the literature to be modified by weight loss [7, 9–11] were found differentially expressed in the present dataset, including *ABCC6*, *ALDOC*, *B4GALT6*, *CCND2*, *CIDEA*, *EGFL6*, *ESR1*, *FADS1*, *FADS2*, *FASN*, *NQO1*, *SCD*, *SFRP4*, *SRPX* and *UCHL1*.

As the surgical procedure included 15 Roux-en-Y gastric bypass (RYGB) and 7 adjustable gastric banding (AGB), we tested whether differential expression was influenced by the type of procedure. No significant interaction between time and procedure was detected (all FDR > 0.24), even though there was a general trend towards larger differences in RYGB than in AGB. Results from the full negative binomial regression model including all covariates and interaction are shown in S2A Table and results for RYBG patients only are shown in S2B Table. Despite a smaller sample size, there were more genes that were differentially expressed in the RYGB group than in the whole population (1375 vs 1214 at 5% FDR).

Gene Set Enrichment Analysis (GSEA) of genes differentially expressed by time

We applied GSEA [15] to explore functional pathways enriched among the most differentially expressed genes by time. For this analysis, genes were ranked according to their signed t-statistics obtained from the negative binomial model. Among genes upregulated at T3, 28 gene sets had a significant enrichment score at 5% FDR, 7 of them remaining after excluding overlapping gene sets (Fig 1A, S3 Table). The 2 most significantly enriched categories ("Ribosomes" and "mRNA splicing") reflected activation of the transcriptional and protein synthesis machineries. The category "Antigen processing/presentation" was driven by *HSPA1A* and *HSPA1B*, 2 members of the HSP70 protein family which mediate the folding of newly translated proteins, and *PSME1*, an interferon (IFN) γ-induced gene implicated in immune proteasome assembly and required for efficient antigen processing. The category "IFN-γ signaling" was driven by *PTPN2* and *IFNGR2* which both modulate IFN-γ signal transduction. These 2 categories were also enriched in genes of the class II/III major histocompatibility complex. The 3 remaining categories reflected downregulation of SMAD and MDA5 signaling pathways, and were all driven by *RPS27A*, a ribosomal gene encoding an ubiquitin gene [16] possibly involved in protein degradation.

When focusing on genes downregulated at T3, 11 gene sets had a significant enrichment score, 7 of them remaining after excluding overlapping gene sets (Fig 1B, S4 Table). A first category was composed of genes involved in ECM organization, including several collagen and metalloproteinase genes. Two categories included genes involved in the biosynthesis of cholesterol and unsaturated fatty acids, respectively, and one category corresponded to the ATP-binding cassette family of proteins involved in cholesterol transport. A category "Packaging telomere ends" included several histone genes (HIST symbol) as well as *TINF2* encoding shelterin, a protein involved in maintaining telomere integrity [17]. The two last categories concerned metabolism of porphyrin and chlorophyll, and metabolism of starch and sucrose, respectively.
A

**Genes upregulated at T3**

<figure><img src="image_5.png" /><figcaption>**Fig 1. Pathways enriched in genes upregulated (A) and downregulated (B) after surgery-induced weight loss.** The x axis indicates the significance of the GSEA score test for enrichment in the given pathway.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g001

B

**Genes downregulated at T3**

<figure><img src="image_6.png" /></figure>

A last analysis considering genes ranked according to their *absolute* t-statistics indicated 6 enriched gene sets, 2 of them being not detected by analyses separating up- and downregulated genes: one category included genes involved in bile acid metabolism (9 genes); the other category included genes involved in activation of the NF-κB pathway by TAK1 (8 genes, the top one being IRAK2) (S5 Table).

**Weighted gene coexpression network analysis (WGCNA) at baseline**

To get deeper insight into the biology of SAT, we sought to identify networks of coexpressed genes by applying WGCNA [18] to the 10,000 most variable gene expressions at baseline. Genes not selected for this analysis, i.e. the 5972 least variable genes, were highly enriched in nuclear genes (P = 7.1 10⁻⁵⁰).

Twenty-seven gene coexpression modules were detected by WGCNA. The grey module comprising genes not being assigned to any module was not further considered. The 26 remaining modules, whose size varied from 56 to 1901 genes, could be regrouped into 4 large clusters (Fig 2A). The proportion of DE genes within each module varied from 0.5% to 28.4%. Interestingly, the 5 modules showing the highest proportions of DE genes (>20%) tended to cluster together (Fig 2A, red square).

When examining the relationships between modules and phenotypic traits measured at baseline, a structure into 4 blocks reflecting the 4 aforementioned clusters was apparent (Fig 3). Clusters 2 and 3 had a mirror-like structure, cluster 2 being positively associated to obesity- and inflammation-related traits (BMI, fat mass, lean mass, adipocyte volume, C-reactive protein [CRP], interleukin[IL]-6, leptin and adiponectin) and negatively related to lipids and diabetes-related traits (total, LDL and HDL cholesterol, triglycerides, metabolic syndrome [MetS], glucose, insulin and Homeostasis Model Assessment of Insulin Resistance [HOMA-IR] index), while cluster 3 showed associations in the opposite direction. Clusters 1 and 4 showed a less
A: Clustering of modules eigengenes

<figure><img src="image_7.png" /></figure>

**Fig 2.** Gene coexpression networks at baseline identified by WGCNA. (A) Hierarchical clustering of module eigengenes. The open red square indicates the 5 modules with the highest proportions of differentially expressed genes (>20%). (B) Stability of modules by Jackknife analysis. The x axis indicates the mean Jaccard index and the y axis the proportion of genes correctly classified (i.e. attributed to the proper module in ≥ 18 out of 22 sub-analyses). The ellipse surrounds the most stable modules.

doi:10.1371/journal.pone.0125718.g002

B: Stability of modules

<figure><img src="image_8.png" /></figure>

clearcut structure. The complete lists of genes belonging to the different modules are provided in S6 Table.

The stability of modules was assessed by the Jackknife technique. Stability criteria that were used were the Jaccard index, the proportion of correctly classified genes and the mean eigen-gene correlation (see Methods). Worthy of note, modules including the highest proportion of DE genes showed a relatively low stability (S7 Table). For further biological interpretation, we focused on the most stable modules as shown in Fig 2B.

Biological significance of RNA-seq modules at baseline

We examined modules according to several criteria: 1) the Gene Ontology (GO) functional categories significantly enriched; 2) the detailed exploration of the network formed by the top 30 hubgenes displayed as circle plot; and 3) the Spearman's rank correlation of the module eigen-gene with baseline phenotypic traits. For each module, gene module membership, stability score and genes belonging to enriched categories are given in S8 Table.

The *turquoise* module was the largest module. It was enriched for the GO process "immune response" (FDR = 7.9 10<sup>-16</sup>) and was associated to a very dense network of hubgenes (Fig 4). The top hubgenes were *MAPK13*, a member of the MAP kinase signal transduction pathway, and *FERMT3*, which encodes a protein of the kindlin family involved in integrin activation and cell adhesion. This module included many inflammatory and immune proteins such as CD molecules (including *CD14* and *CD33*, 2 monocyte markers), chemokines and proteins of the complement system. This module may reflect the general obesity-associated immune response characterized by monocyte recruitment and macrophage infiltration. The eigengene positively correlated to lipid variables including total cholesterol (r = 0.53) and triglycerides (r = 0.68).

The *magenta* and *tan* modules were close to each other in the hierarchical clustering (Fig 2A). The *magenta* module was a highly connected module enriched for the GO category "cell
<figure><img src="image_9.png" /><figcaption>**Fig 3.** Correlations between module eigengenes and patient's characteristics at baseline. Clinical and biological traits (rows) were regrouped in 2 categories (obesity and inflammation-related traits, glucose and lipid homeostasis-related traits). Module eigengenes (columns) were clustered as in Fig 2. In the heatmap, red color represents positive correlation and green color negative correlation.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g003

surface receptor linked signal transduction" (FDR = 2 10⁻⁹) which comprised many genes coding for cytokine receptors as well as CD molecules. The top 30 hubgenes (Fig 4) were present in all sub-sample analyses, indicating a high stability of the network. Among them several genes were involved in the lymphocyte-activation pathway: *LCK* encoding the lymphocyte-specific protein tyrosine kinase, *CD79A* and *MS4A1* (alias *CD20*)—two B-lymphocyte antigens strongly linked one with each other, and *CARD11*—a molecule with a caspase domain essential for T-cell activation. The module also contained *CD8A* and *CD8B* characterizing the specific subpopulation of CD8⁺ T cells that play an essential role in the initiation and propagation of adipose inflammation in rodents [19]. The *magenta* module was additionally enriched for the GO class "hemoglobin complex" (P = 1.8 10⁻⁶). This enrichment was due to the inclusion of several hemoglobin genes (*HBA1*, *HBA2*, *HBB*, *HBD*, *HBM*, *HBQ1*, *HBG1*, *HBG2*) and erythrocyte-specific genes (*AHSP*, *ALAS2*, *SLC4A1*). The *magenta* module eigengene tended to be positively related to all obesity-related traits, although only the correlation with adipocyte volume reached significance (r = 0.49).
<figure><img src="image_10.png" /><figcaption>**Fig 4. Network circle plots of immune-response modules at baseline.** The figure depicts the connectivity patterns between the top 30 hubgenes of each module. Genes are ranked in the anticlockwise direction according to decreasing module membership. The size of each black circle indicates the gene connectivity (number of connected genes). An orange line indicates a positive correlation and a blue line a negative correlation between any 2 genes. The *turquoise* module reflects a general immune response, the *magenta* module reflects lymphocyte activation, the *tan* module reflects neutrophil activation and the *plum* module reflects the interferon signaling pathway.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g004

The *tan* module was enriched for the GO category "defense response" (FDR = 6.5 10⁻¹¹) which included several receptors of the innate immune response such as chemokine receptors or Toll-like receptors. Several of these proteins act as monocyte chemoattractants. The dense network of the 30 hubgenes (Fig 4) was driven by 2 genes, *S100A8* and *S100A9* (calgranulin A and B, respectively), which form a heterodimer *S100A8/A9* (or calprotectin) originally discovered as an immunogenic protein expressed and secreted by neutrophils [20]. The module also contained the gene encoding *S100A12* (calgranulin C). Calgranulins bind to pattern recognition receptors and serve as trigger factors of inflammation in concert with *IL1B*, another member of the module [21]. The 3rd hubgene of the module was *CSF3R* which encodes the receptor
<figure><img src="image_11.png" /><figcaption>**Fig 5.** Network circle plots of other modules at baseline. See legend of Fig 4. The *darkmagenta* reflects ectoderm, the *skyblue* and *steelblue* modules reflect the extracellular matrix and the *violet* module reflect lipid metabolism.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g005

for colony stimulating factor 3, the principal cytokine controlling neutrophil development and function [22]. The *tan* module might then reflect coexpression of genes characterizing the inflammation associated to neutrophil activation and monocyte infiltration of adipose tissue. As the *magenta* module, the *tan* module tended to be positively related to obesity-related traits, but in addition it significantly correlated to MetS-related traits, in particular triglycerides levels (r = 0.46).

The *darkmagenta* module was a small module strongly enriched in the GO biological process "ectoderm development" (FDR = 1.8 10⁻¹⁹). Several genes of this module coded for proteins involved in the differentiation of keratinocytes, such as keratins, intermediate filaments and proteins of the cornified envelope (Fig 5). The top hubgene *KRT15* has been implicated in breast and lung cancers [23, 24]. Among the major contributors of this module were also *FLG* coding for fillagrin, a protein that aggregates keratin intermediate filaments, and *TFAP2A*, a
transcription factor expressed in the ectoderm. A small sub-network was formed by *KRT2*, *SBSN*, *LGALS7B*, *LY6D*, *KRTDAP* and *CALML5*, all genes coding for proteins involved in differentiation of keratinocytes. The *darkmagenta* module was negatively associated to CRP levels (r = -0.46).

The *paleturquoise* module was a large module moderately enriched for the GO biological process "transcription" (FDR = 2.5 10⁻²). Thirty-eight genes in the module were zinc finger proteins (ZNF symbol) which are major transcriptional regulators whose role in adipogenesis is increasingly recognized [25]. The top hubgene (data not shown) was *NIPBL* encoding delangin, a member of the cohesin complex which plays a key role in transcription [26]. Other hubgenes were *ARFGEF2*, an exchange factor involved in the trans-Golgi trafficking of proteins [27] and *NCOA2* (alias *TIF2*) which modulates energy metabolism [28]. The *paleturquoise* module negatively correlated to obesity traits.

The *skyblue* and *steelblue* modules were close to each other in the hierarchical clustering (Fig 2A) and they were both enriched for the GO component "extracellular region" (FDR = 6.4 10⁻³ and 1.8 10⁻⁵, respectively). These modules contained several genes coding for ECM constituents, in particular collagens, as well as proteins involved in cell adhesion, cell-cell interaction, peptidase activity and cytokine activity. The *skyblue* module was associated to a network of weakly interconnected hubgenes (Fig 5). The strongest correlation was between *PAMR1—a* gene reported to modulate phagocytosis, a process contributing to normal tissue remodeling [29], and *SCARA5—a* scavenger receptor with a collagen-like domain suspected to be involved in adipogenesis [30]. The *skyblue* module positively correlated to glucose (r = 0.45) and negatively to obesity-related traits.

By contrast, the *steelblue* module was associated to a network of strongly interconnected hubgenes (Fig 5) which exhibited a great stability. Among the most highly connected genes were *PPL* (periplakin)—a component of the cornified envelope, *PDGFRA—the* platelet-derived growth factor receptor, *DCN* (decorin)—a component of connective tissue, *C1R* and *C1S* — two subcomponents of the complement molecule C1, and *CTSK* (cathepsin K) which is involved in ECM remodeling and was strongly associated to *ECM1* (extracellular matrix protein 1). As the *skyblue* module, the *steelblue* module negatively correlated to obesity-related traits, especially BMI, adipocyte volume and leptin (r = -0.46, -0.48 and -0.51, respectively).

The *violet* module was enriched in the GO processes "lipid biosynthesis" (FDR = 6.3 10⁻⁹), "cholesterol metabolism" (FDR = 8.2 10⁻⁶) and "fatty acid metabolism" (FDR = 1.4 10⁻⁵). The top hubgenes were highly interconnected (Fig 5). In particular, strong correlations were observed between *INSIG1* that plays a key role in the control of cholesterol synthesis and *ACLY*, *FAD5* and *IDI1*, all involved in cholesterol and fatty acid synthesis. Another small highly interconnected sub-network was constituted of *ELOVL5*, *ME1* and *SPTSSA*, the 2 former being involved in fatty acid synthesis and the latter in sphingolipid synthesis. *SREBF1*, a transcription factor that binds to the sterol regulatory element-1 (SRE-1) and regulates lipid metabolism was also a top hubgene. Among the top 30 hubgenes, 16 were significantly downregulated at T3, a result consistent with the aforementioned enrichment of lipid pathways in differentially expressed genes. The *violet* module positively correlated to CRP (r = 0.49) and negatively, although not significantly, to lipid traits.

Finally, the *plum* module, despite a moderate stability (Fig 2B), deserved a special interest because of its strong enrichment for the GO process "response to virus" (FDR = 1.4 10⁻⁴). Almost all positive hubgenes of the network encoded proteins involved in the IFN-induced pathway and this sub-network was highly stable. *CRY1* encoding a circadian protein was negatively correlated to IFN-induced genes (Fig 4). The *plum* module positively correlated to adiponectin (r = 0.49).
<figure><img src="image_12.png" /><figcaption>**Fig 6. Modifications of immune-response networks after surgery-induced weight loss.** Heatmaps representing pairwise correlation matrices of gene expressions are displayed for each module. At T3, genes are ranked in the same order as in T0. Only genes having the highest membership to each module are displayed, except for the *plum* module where all genes are shown.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g006

Remarkably, the 3 most stable modules (*turquoise*, *magenta* and *tan*), to which could be added the IFN-related *plum* sub-module, were all reflecting innate and adaptive immune pathways. Because of the key role of immune cells in the low-grade inflammation associated to obesity [2], we focused further analyses on these modules.

## Modifications of immune response-related modules after surgery

We investigated whether the improvement of the inflammatory profile following surgery might be reflected by modifications of immune cell pathways. For this purpose, correlation heatmaps of immune response-related modules were compared between T0 and T3. For each module, we focused on genes having the highest contribution (module membership) to the module. A general view of the network changes is displayed in Fig 6 while detailed heatmaps at T3 are given in S1, S2, S3 and S4 Figs.

The *turquoise* module, reflecting a general immune response, conserved a highly dense connected network at T3, except a few genes losing association with the network (S1 Fig). The
most striking one was *DUSP2* which encodes a member of the dual specificity protein phosphatase subfamily that regulates MAPK signaling and has a critical role in inflammatory response of macrophages [31]. Other genes disconnected from the module at T3 were *KCNN4* coding for a component of a calcium-activated potassium channel in T-lymphocytes, and *CD300C* coding for a member of the CD300 family of receptors that regulate immune cell processes [32].

The *magenta* module associated to lymphocyte expression showed profound perturbations after surgery (S2 Fig). The dense matrix of expressions at T0 was split into several disconnected subsets at T3. A first highly preserved sub-module comprised several hubgenes which characterized lymphocyte activation. The only gene disconnected from this block at T3 was *CMTM5*—a gene belonging to the chemokine-like factor gene superfamily with putative tumor suppressor activity. A second subset contained genes involved in signal transduction such as chemokine receptors (*CXCR1*, *CXCR2*, *CCR7*), G-protein-coupled receptor (*GPR97*) and interleukin receptors (*IL7R*, *IL18RAP*). This group exhibited looser correlations at T3, especially *CXCR1* and *CXCR2* were no longer associated with other module members. Worthy of note, the membership of hemoglobin and erythrocyte-specific genes to the module at T0 disappeared at T3. This loss of correlation was accompanied by a marked decrease of hemoglobin gene expressions at T3. These observations may be explained by a higher blood cell contamination at baseline.

The *tan* module reflecting neutrophil expression also showed important modifications after surgery (S3 Fig). Although *S100A8* and *S100A9* were still coexpressed at T3, *S100A8* exhibited a lower pattern of correlations with other genes than *S100A9*, suggesting a dissociation of the dimer *S100A8/A9* involved in neutrophil activation. Similarly *S00A12*—the 3rd calgranulin gene—lost association with other module members. Several other hubgenes which were densely connected at T0 were no longer coexpressed at T3, such as *VNN2* encoding vanin 2—a protein involved in migration of neutrophils, *TUBB1* (tubulin β1 involved in cell migration), *FAM65B* which dampens chemokine-induced T cell migration [33], *APOBEC3A* which mediates antiviral immunity, and *FPR2* which encodes a receptor mediating activation of neutrophils.

A small *tan* sub-module highly conserved at T3 included several genes involved in the immune response mediated by dendritic cells (*CLEC7A*, *CLECL1*, *TLR1*, *TLR6*, *LY75*, *PIK3CG*, *SIGLEC10* and *LILRA5*). Interestingly, *CD300E*, another member of the CD300 receptors which, contrary to the ubiquitously leukocyte-expressed *CD300C*, is specifically expressed in monocytes and dentritic cells, completely lost its association with the sub-module (S3 Fig).

Within the *plum* module, the highly interconnected cluster of genes reflecting the IFN-signaling pathway was strongly preserved at T3. There were, however, two notable exceptions: *DDX60* and *LURAP1* (S4 Fig). *DDX60* positively regulates *DDX58* and *IFIH1*-dependent type I interferon response to viral infection [34]. *LURAP1* (previously *C1orf190*) [35] is an activator of the NF-κB signaling pathway which plays a crucial role in obesity-induced inflammation [36].

We next investigated to what extent the main genes cited above were also differentially expressed after surgery (Table 2). Almost all these genes showed a tendency towards downregulation at T3, and 9 of them exhibited a nominally significant decrease of expression after surgery (P<0.05).

Lastly, we checked whether modifications were dependent of the type of surgical procedure. For this purpose, we compared the module membership of genes (i.e. the correlation with their module eigengene) at T0 and T3, separately in RYBG and AGB. This analysis showed that the decrease of correlation of *CXCR1*, *CXCR2*, *GPR97*, *S100A8* and *S100A12* to their module eigen-gene was very similar after RYBG and AGB (Fig 7).
<table><caption>Table 2. Differential expression by time of the main genes disconnected from immune-response modules after surgery.</caption><thead><tr><th>Gene ID</th><th>Module</th><th>Time coefficient</th><th>Nominal P-value</th><th>FDR</th><th>Mean counts T0</th><th>Mean counts T3</th></tr></thead><tbody><tr><td>S100A8</td><td>Neutrophils</td><td>-0.708</td><td>0.003</td><td>0.047</td><td>348</td><td>165</td></tr><tr><td>CXCR1</td><td>Lymphocytes</td><td>-0.678</td><td>0.005</td><td>0.058</td><td>33</td><td>18</td></tr><tr><td>S100A12</td><td>Neutrophils</td><td>-0.815</td><td>0.005</td><td>0.062</td><td>28</td><td>12</td></tr><tr><td>FPR2</td><td>Neutrophils</td><td>-0.849</td><td>0.007</td><td>0.069</td><td>17</td><td>7</td></tr><tr><td>CD300E</td><td>Neutrophils</td><td>-0.467</td><td>0.009</td><td>0.084</td><td>92</td><td>58</td></tr><tr><td>FAM65B</td><td>Neutrophils</td><td>-0.329</td><td>0.017</td><td>0.119</td><td>162</td><td>126</td></tr><tr><td>CXCR2</td><td>Lymphocytes</td><td>-0.385</td><td>0.025</td><td>0.149</td><td>76</td><td>54</td></tr><tr><td>CCR7</td><td>Lymphocytes</td><td>-0.264</td><td>0.034</td><td>0.179</td><td>82</td><td>69</td></tr><tr><td>IL7R</td><td>Lymphocytes</td><td>-0.442</td><td>0.044</td><td>0.209</td><td>44</td><td>29</td></tr><tr><td>S100A9</td><td>Neutrophils</td><td>-0.324</td><td>0.092</td><td>0.314</td><td>721</td><td>457</td></tr><tr><td>APOBEC3A</td><td>Neutrophils</td><td>-0.212</td><td>0.104</td><td>0.337</td><td>78</td><td>66</td></tr><tr><td>DUSP2</td><td>Macro/mono</td><td>-0.240</td><td>0.158</td><td>0.420</td><td>65</td><td>54</td></tr><tr><td>VNN2</td><td>Neutrophils</td><td>-0.198</td><td>0.167</td><td>0.431</td><td>58</td><td>46</td></tr><tr><td>CMTM5</td><td>Macro/mono</td><td>-0.173</td><td>0.232</td><td>0.509</td><td>11</td><td>10</td></tr><tr><td>GPR97</td><td>Lymphocytes</td><td>-0.186</td><td>0.432</td><td>0.690</td><td>16</td><td>13</td></tr><tr><td>LURAP1</td><td>Interferon</td><td>0.070</td><td>0.582</td><td>0.794</td><td>29</td><td>33</td></tr><tr><td>DDX60</td><td>Interferon</td><td>-0.024</td><td>0.621</td><td>0.817</td><td>549</td><td>536</td></tr><tr><td>KCNN4</td><td>Macro/mono</td><td>-0.083</td><td>0.627</td><td>0.821</td><td>142</td><td>143</td></tr><tr><td>CD300C</td><td>Macro/mono</td><td>0.071</td><td>0.688</td><td>0.857</td><td>67</td><td>72</td></tr><tr><td>IL18RAP</td><td>Lymphocytes</td><td>-0.049</td><td>0.735</td><td>0.880</td><td>41</td><td>43</td></tr><tr><td>TUBB1</td><td>Neutrophils</td><td>0.069</td><td>0.741</td><td>0.883</td><td>100</td><td>106</td></tr></tbody></table>

doi:10.1371/journal.pone.0125718.t002

## The IFN-signaling module in obese versus lean mice

For sake of validation, we applied WGCNA to transcriptomic data of adipose tissue from obese mice fed a high fat/high sucrose (HF/HS) diet. Among 9 detected modules, one module

<figure><img src="image_13.png" /><figcaption>Fig 7. Module membership of some genes of interest before (T0) and after surgery (T3), separately in RYGB and AGB. The module membership of a gene represents the correlation to its module eigengene. For each module, the eigengene at T3 was calculated on the same gene set as at T0.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g007
<figure><img src="image_14.png" /><figcaption>**Fig 8.** Sub-networks of expression around *Ddx60* and *Dhx58* genes in obese versus lean mice. Comparison of obese mice fed a high-fat/high-sucrose (HF/HS) diet and lean mice fed a chow diet (microarray data). Only edges corresponding to significant correlations are shown (P<0.05 and cor>0.4). The thicker and shorter edges correspond to stronger correlations. Networks were visualized using the gplot function of the sna R package.</figcaption></figure>

doi:10.1371/journal.pone.0125718.g008

including 63 genes was highly enriched for genes of the IFN-signaling pathway (P = 1.6 10⁻⁶). In particular, this module included 9 genes common to the human *plum* module: *Cmpk2*, *Ddx60*, *Ifi44*, *Ifih1*, *Ifit1*, *Irf7*, *Mx1*, *Oas2* and *Rsad2*. This module showed a high preservation in lean mouse fed a chow-diet (S5 Fig). However, similar to what was observed in humans after bariatric surgery, most of the connections of *Ddx60* in the obese state were lost in the lean state (Fig 8A). A similar pattern was observed for *Dhx58*, another regulator of the *DDX58* and *IFIH1*-mediated signaling cascades (Fig 8B).

Discussion

In the present study, the impact of bariatric surgery on SAT transcriptome was studied by two approaches: first by measuring the modifications in single gene expressions after surgery, and second by investigating the perturbations in networks of coexpressed genes. The first approach informs about the difference in the abundance of a given molecule, while the second approach provides insight into mechanisms of rewiring of transcriptional networks in response to surgery-induced weight loss. While the two approaches are complementary, the latter one has the advantage of providing a more informative picture of the dynamic changes in the gene regulatory networks [37].

Using a stringent 5% FDR, we identified 1214 genes whose expression significantly varied between T0 and T3. Upregulated genes were mostly involved in the basal cellular machinery and reflected increased protein synthesis and decreased degradation as already reported [8, 11]. This enhanced biogenesis is probably a requirement of the organism for counteracting the caloric restriction caused by the surgical procedure. The upregulation of the IFN-γ induced pathway may be explained by the fact that immunoproteins are thought to be highly demanding in terms of energy and may be particularly affected by energy restriction [38]. As already known from microarray studies [7, 8, 11], downregulated genes were enriched in major metabolic functions of adipose tissue such as ECM organization and metabolism of lipids, fatty
acids and carbohydrates. The enrichment of downregulated genes in the KEGG category "Packaging telomere ends" is of particular interest given the inverse link reported between adiposity and telomere length [39].

By WGCNA, we identified 26 modules of coexpressed genes at baseline, 8 of them showing a great stability by Jackknife analysis. Worthy of note, modules having the largest proportion of DE genes were among those having the lowest stability. This finding may be explained by the fact that DE genes are genes that are less constrained in terms of variability than those that are tightly regulated within systems sustaining homeostasis of the organism. In keeping with this interpretation, we have previously shown that gene expressions that are under the influence of genetic variants are more likely to be also modulated by non-genetic factors [40].

The 3 most stable modules at T0 (*turquoise*, *magenta* and *tan*) were related to the immune response. Immune modules are thought to reflect the low-grade chronic inflammation associated to obesity. This immune signature shared by many common chronic diseases and observed across several tissues has been shown to be driven by a few key genes [41]. Worthy of note, all the predicted key drivers of the "inflammatome" regulatory network described in human adipose tissue [41] were among the major contributors of the *turquoise* module.

The 3 immune response-related modules reflected the diversity of immune cells involved in adipose tissue inflammation. The *turquoise* module was a large module reflecting a non-specific signature of innate immune cells, mainly macrophages which are the most abundant immune cells in adipose tissue, but also monocytes. The *magenta* module reflected the adaptive immune response elicited by T lymphocytes which have an essential role in the propagation of adipose inflammation [19]. The *tan* module characterized the activity of neutrophils which play a primary role in the initiation of inflammation by recruiting macrophages, dendritic cells and lymphocytes to fat depots [3].

Three months after surgery, the 3 immune response-related modules exhibited distinct patterns of modifications. The *turquoise* module was strongly preserved with very few exceptions of disconnected genes at T3. Among them, *DUPS2* is of particular interest for its potential as a therapeutic target to modulate immune effectors functions [31]. The dissociation of *CD300C* from the *turquoise* module, parallel to that of *CD300E* from the *tan* module, is also worthy of note knowing that CD300 molecules are important modulators of the immune system [32].

In the *magenta* module, the most striking feature was the disconnection between genes involved in T-cell signaling and receptors mediating the signal transduction such as *CXCR1*, *CXCR2*, *CCR7*, *IL7R* and *GPR97*.

In the *tan* module, the major modification induced by surgery was a downregulation of *S100A8* and *S100A12* and their dissociation from other genes of the module. These proteins are secreted by neutrophils upon activation of the immune system and participate in neutrophil-mediated inflammation [21]. Several other genes losing their module correlation were involved in neutrophil chemotaxis and migration (*FPR2*, *VNN2*, *TUBB1*, *FAM65B*). These results suggest a critical role of these factors in dampening the migration and the activity of neutrophils in adipose tissue after surgery.

Interestingly, several of these genes (*S100A8*, *S100A12*, *CXCR1*, *CXCR2*, *CCR7*, *IL7R*, *FPR2*, *CD300E* and *FAM65B*), besides being disconnected from their module after surgery, were also significantly downregulated, suggesting that the encoded molecules—mostly receptors and ligands—play a key role in the rewiring of inflammatory transcriptional networks after surgery-induced weight loss.

The present study also highlighted a small network of strongly interconnected genes which comprised 19 IFN-related genes. A similar module was identified in circulating monocytes of the same obese subjects (data not shown) and was also present in mice fed an obesogenic diet. One member of the network was *IRF7*, a master regulator of type I IFN-dependent immune
response [42]. IRF7 has been identified as the key driver of an inflammatory network associated to type 1 diabetes [43], which greatly overlaps with the IFN network found in the present study. IRF7 has also been shown to contribute to diet-induced obesity and insulin resistance [44]. Interestingly, the strongest negative contributor of the plum module—reproducibly found in almost all sub-analyses—was CRY1, a circadian gene. The negative association of CRY1 with IFN-induced genes is consistent with the protective role of CRY1 in the modulation of inflammation by the circadian clock [45]. A striking finding was the dissociation of DDX60 from several IFN-induced genes after surgery. A similar result was observed when lean mice were compared to obese mice. DDX60 has been recently described as a novel antiviral factor involved in RIG-I-mediated IFN signaling [34]. The present finding suggests that DDX60 may represent a potential therapeutic target for attenuating obesity-associated inflammation.

Several limitations of the present study should be acknowledged. First, it was based on a limited sample size, even though this is the largest study, to our knowledge, using the RNA-seq technology for assessing the transcriptome of human adipose tissue. Second, the visceral fat compartment is known to become more pro-inflammatory during obesity that the subcutaneous fat compartment. However, visceral fat sampling needs an invasive surgical biopsy and therefore could not be obtained in the post-surgical follow-up. Third, adipose tissue represents a heterogeneous mixture of cells composed of adipocytes and stromal vascular cells, and the cell distribution may differ across individuals and before and after surgery. The lack of tissue samples for immunohistochemical or cytometry analyses did not allow us to quantify the abundance of the different types of immune cells. Fourth, adipose tissue remodeling after bariatric surgery is a very dynamic process, and analyzing the expression profile at a single postsurgical time point only represents a snapshot that may not illustrate the full picture of adipose remodeling. Fifth, to keep gender homogeneity, the study was restricted to women who represent the vast majority of patients undergoing bariatric surgery, but this may have introduced a source of variability induced by menstrual cycle variations and/or altered sex hormones on fat biology. Finally, functional experiments targeting the genes in the pathways highlighted by the present study are needed to firmly establish their role in the modulation of inflammation after bariatric surgery.

# Methods

## Subjects

The study population initially included 24 obese women undergoing bariatric surgery recruited at Pitié-Salpêtrière Hospital (Nutrition Department, Institute of Cardiometabolism and Nutrition, ICAN). Following usual guidelines, all patients met the criteria for needing obesity surgery: these included a BMI ≥40 kg/m² or ≥35 kg/m² with at least one comorbidity (hypertension, T2D, dyslipidemia or obstructive sleep-apnea syndrome). Preoperative evaluation included a detailed medical history, a physical examination, and assessment of nutritional, metabolic, cardiopulmonary, vascular and psychological factors. The weight of the subjects had been stable (variation < 2 kg) for at least 3 months prior to surgery. Subjects did not demonstrate evidence of acute or chronic inflammatory disease, infectious diseases, cancer and/or known alcohol consumption (> 20g per day) at any time point. Clinical and biological parameters were assessed prior to surgery (T0) and at 3 months after surgery (T3).

## Ethics statement

All subjects gave their written informed consent to participate in the study. The Ethics Committee of the "Comité de Protection des Personnes Ile de France 1" (Hôtel-Dieu Hospital)
approved the clinical investigation. The study was conducted in accordance with the Helsinki Declaration and was registered in a public trials registry: ID Number NCT01454232.

Fat mass was measured by whole-body fan-beam dual energy x-ray absorptiometry (DXA, Hologic Discovery W, software v12.6; Hologic, Bedford, MA). Subjects were considered to have diabetes if they were taking an antidiabetic treatment or if fasting glycemia was ≥7 mmol/l on at least two different occasions. To estimate the patients cardiometabolic risk, we used the thresholds established by the International Diabetes Federation (IDF, 2005) for the definition of MetS. Subjects were included in the MetS group if they had at least two of the four factors: 1) raised triglyceride levels (≥ 1.7 mmol/L) or specific treatment; 2) reduced HDL-cholesterol < 1.29 mmol/L) or specific treatment; 3) raised blood pressure (systolic ≥ 130 or diastolic ≥ 85 mmHg) or treatment of previously diagnosed hypertension; 4) raised fasting plasma glucose (≥ 5.6 mmol/L) or previously diagnosed T2D.

SAT samples were obtained by needle biopsy from the peri-umbilical area under local anesthesia (1% xylocaïne) and venous blood samples were collected in the morning after a 12-hour fast for routine determination of blood cell counts, glucose, insulin, triglycerides, total cholesterol, LDL-cholesterol, HDL-cholesterol, leptin, adiponectin and IL-6 as previously described [46]. Insulin resistance was assessed using the HOMA-IR index: insulinemia (mUI/L) × fasting blood glucose (mmol/L)/22.5. High sensitive CRP was measured using an immunoturbidimetric assay on a Cobas 501 analyser (Roche Diagnostics, Meylan, France).

## mRNA preparation and sequencing

Total RNA was extracted from SAT using the *mirVana* miRNA Isolation Kit (Ambion). 500 ng of total RNA were used to prepare the library with TruSeq RNA sample Prep Kit V2 (*Illumina*). Briefly, mRNA was selected from total RNA by poly-T oligo-attached magnetic beads and fragmented. Reverse transcription was done with random hexamers by SuperScript II reverse transcriptase. From resulting cDNA, classical library preparation steps included end repair, adenylate 3 ends, adaptors ligation with specific index and 15 cycle-enrichment amplification performed according to manufacturing protocol. The generated libraries were submitted to quality control on High Sensitivity DNA chip from the *Agilent* Bioanalyzer 2100 (Agilent Technologies, Wokingham, UK) and quantified by fluorimetry using the Quant-iT dsDNA Assay Kit (*Life Technology*). A dilution at 10nM was performed for each library and pools of 12 libraries were constituted. Each pool included the libraries from 3 individuals for the 2 time points, each library being identifiable by a specific index. The resulting pools were loaded on a High Sensitivity DNA chip (*Agilent*) and quantified by integration. 9 pM of this solution was used for cluster generation using the TruSeq PE Cluster Kit v3 (*Illumina*) on a CBot instrument. A 75 bp paired-end run was performed on the *Illumina* HiSeq2000 platform with the TruSeq SBS kit V3 (*Illumina*). Sequencing was performed on two 8-laned flow cells. Each flow cell included 4 pools sequenced twice in 2 different lanes. The two flow cells were processed in the same conditions at 2 months of interval (runs 1 and 2, respectively).

## RNA-seq data pre-processing

Reads of bad quality were filtered out using `fastq_illumina_filter`. Trimming was performed using sickle with default parameters (quality threshold = 20, length threshold = 20). Orphan reads were excluded. Quality of reads was checked using FastQC. Paired-end reads were mapped to the human transcriptome using Bowtie2 v2.1.0 (http://bowtie-bio.sourceforge.net/bowtie2/). The reference human transcriptome was generated from the UCSC reference genome (hg19) and the *Illumina* iGenome annotation file of transcripts using RSEM [47]. Quantification of transcript abundance was performed using eXpress [48]. The data used for
analysis were the *eff_counts* which give the estimated number of fragments mapping to a given target sequence adjusted for fragment and length biases [49]. Analyses were performed at the gene level by summing isoform counts for each gene.

## Statistical analysis

After exclusion of 2 individuals failing to be sequenced because of bad quality or insufficient amount of RNA, analysis was performed in 22 subjects in the R environment (v3.0.2). Gene coverage was assessed using the RSeQC package (http://rseqc.sourceforge.net/). The plot showed a well-documented 3' bias due to polyA selection, but there was no clear relationship between the skewness of the coverage and the RNA quality assessed by RIN values (S9 Table, S6 Fig). We further checked by principal component analysis whether samples with low RIN may behave as outliers, but no outlier was detected. Genes with low expression levels, defined by a mean number of counts < 10 at both times, were filtered out. We also excluded a few genes annotated as micro RNAs (MIR symbol), leaving out 15,972 genes for analysis.

**Statistical model for differential analysis by time.** For each gene, the difference of expression level between T0 and T3 was estimated using a negative binomial regression model fitted on read counts. We used the R hglm package [50] which allows fitting fixed and correlated random effects resulting from repeated time measures. Counts were normalized using the size factor estimated by the R DESeq package [51]. This factor accounts for variations of library size among samples and corrects for the bias induced by highly expressed genes. It was introduced as an offset in the negative binomial model. Fixed effects were time (T3 vs T0), age, ethnicity (European vs non European), run, RIN and RNA concentration. Additionally we tested whether time effects differed according to the type of surgical procedure by introducing the surgical procedure (RYGB vs AGB) and an interaction term time\*procedure in the model. Parameters and their standard errors were estimated by the h-likelihood and tested by means of the t-test. To account for multiple testing, we used the FDR [52] and considered as significant an FDR < 0.05.

**Gene Set Enrichment Analysis (GSEA).** The GSEA software (v2.0.12) was used to investigate whether the most differentially expressed (DE) genes between T0 and T3 shared common biological functions [15]. The analysis was performed on the list of all genes ranked either by their signed or by their absolute t-statistics obtained from the negative binomial model. Parameters were set to their default value. Gene sets tested for enrichment were taken from the MSigDB database (BioCarta, KEGG and REACTOME gene sets). The statistical significance of the enrichment score was assessed by performing 1000 permutations. Overlap among gene sets was examined by a leading edge analysis (www.broadinstitute.org/gsea/). We considered as significant a FDR < 0.05.

**Weighted Gene Coexpression Network Analysis (WGCNA).** The R WGCNA package (v1.34) was used to identify networks of coexpressed genes at T0 [18]. Prior to analysis, data were transformed using the VST transformation proposed by the DESeq package. The resulting variable was further adjusted for run, RIN and RNA concentration and the 10,000 most variable gene expressions at T0 were used for analysis. We used the adjacency matrix to specify the connection strength $a_{ij}$ between genes *i* and *j* ($a_{ij} = |\text{cor}(x_i, x_j)|^\beta$). The power $\beta$ was fixed to 6 after applying the function pickSoftThreshold. The topological overlap measure was used to define modules of interconnected genes.

Modules were merged using the Dynamic Tree Cut function [53]. Parameters were fixed to their default values. Spearman's rank correlation coefficients between module eigengenes and phenotypic traits were represented by a correlation heatmap. The module membership (MM)
measured the correlation of a given gene expression profile with its module eigengene, genes with the highest rank of MM having the strongest contribution to the module.

Enrichment analysis of modules in functional GO categories was performed using the DAVID software (v6.7) [54] with the list of the 15,972 genes expressed in SAT taken as reference. A FDR < 0.05 was considered significant. For each module, the network of the 30 genes having the highest absolute MM in the module, thereafter named hubgenes [18], was displayed using the R function circlePlot.

Stability of modules was assessed by the Jackknife technique. The WGCNA algorithm was repeatedly run on the 22 subsets obtained by systematically leaving out 1 individual from the dataset and in each sub-analysis, we selected the module having the maximum overlap with each original module by means of the Jaccard index [55]. The mean and standard deviation of the maximum Jaccard indices were calculated over the 22 subsets of n-1 individuals. In order to select the most stable modules, for each gene *i* we computed the number $k_i$ of times it was correctly classified in the best matching module of the different sub-analyses. We considered that a gene was steadily clustered if it was attributed to the right original module in ≥ 18 of the 22 sub-analyses. We repeated the same procedure by using the correlation between eigengenes instead of the Jaccard index, i.e. we selected the module having the maximum eigengene correlation with the original module.

## Microarray data in obese versus lean mice

These data have been previously described [56]. Briefly, mice from the Hybrid Mouse Diversity Panel including more than 100 inbred mouse strains were purchased from the Jackson Laboratory and bred in control (chow diet) and “obese” conditions at University of California, Los Angeles. Male mice of the obese group were maintained on a chow diet (Ralston Purina Company) until 8 weeks of age when they were given a high-fat/high-sucrose (HF/HS) diet (Research Diets-D12266B) with the following composition: 16.8% kcal protein, 51.4% kcal carbohydrate, 31.8% kcal fat. Age-matched control male mice were fed a chow diet (18% kcal fat) during 16 weeks. Total RNA was extracted from flash frozen epididymal adipose samples using Qiazol and the Adipose Tissue RNAeasy kit (Qiagen) according to the manufacturer's instructions. Gene expression levels were assessed using *Affymetrix Mouse Genome HT_MG-430A* arrays and data were processed as previously described [57]. To investigate the preservation of gene expression modules between HF/HS and chow diet-fed mice, we applied the same pipeline of analysis as in humans. Modules were first detected in the HF/HS group by WGCNA and the within-module correlations were then compared to those in the chow diet group.

## Supporting Information

**S1 Fig. Modifications of the *turquoise* network after surgery-induced weight loss.** The heat-map representing the pairwise correlation matrix of gene expressions is shown. Only genes having the highest membership to each module are displayed. For purpose of comparison, genes are ranked in the same order as at T0.
(PDF)

**S2 Fig. Modifications of the *magenta* network after surgery-induced weight loss.** See legend
of S1 Fig.
(PDF)
**S3 Fig. Modifications of the *tan* network after surgery-induced weight loss.** See legend of S1 Fig.

(PDF)

**S4 Fig. Modifications of the *plum* network after surgery-induced weight loss.** See legend of S1 Fig.

(PDF)

**S5 Fig. Heatmaps of the interferon-related module in obese versus lean mice.** Obese mice were fed a high-fat/high-sucrose (HF/HS) diet and lean mice were fed a chow diet. Heatmaps represent the pairwise correlation matrix of gene expressions assessed by microarray.

(PDF)

**S6 Fig. Plot of gene coverage of the 44 samples (22 individuals at 2 times).** (A) Reads coverage according to position in the gene from 5' to 3'; (B) RIN values of the samples ranked by increasing 3' skewness of coverage.

(TIF)

**S1 Table. Comparison of expression levels between T0 and T3 for the 15,972 genes expressed in the subcutaneous adipose tissue.**

(XLSX)

**S2 Table. Parameter estimates from the full negative binomial regression including age, ethnic origin, run, RIN, RNA concentration, type of surgical procedure, time and time-*procedure*.**

(XLSX)

**S3 Table. Lists of upregulated genes contributing to enrichment in KEGG or REACTOME pathways by GSEA.** The first sheet shows the GSEA report and each following sheet corresponds to an enriched pathway.

(XLSX)

**S4 Table. Lists of downregulated genes contributing to enrichment in KEGG or REACTOME pathways by GSEA.** The first sheet shows the GSEA report and each following sheet corresponds to an enriched pathway.

(XLSX)

**S5 Table. Lists of genes contributing to enrichment in KEGG or REACTOME pathways in the GSEA analysis combining up- and downregulated genes.** The first sheet shows the GSEA report and each following sheet corresponds to an enriched pathway.

(XLSX)

**S6 Table. Description of modules identified by WGCNA at baseline.**

(XLSX)

**S7 Table. Stability of modules identified at baseline assessed by the Jackknife technique leaving out 1 individual at each time.**

(XLSX)

**S8 Table. List of genes in the most stable modules.** Each sheet of the Excel file corresponds to a different module.

(XLSX)

**S9 Table. RIN values of SAT samples at T0 and T3.**

(XLSX)
**Acknowledgments**

We thank Drs. Florence Marchelli and Valentine Lemoine who were involved in patient recruitment at the Center for Research on Human Nutrition, Pitié-Salpêtrière Hospital, Paris, France, and Michèle Guerre-Millo, Joan Tordjman, Sébastien André, Fabrizio Andrelli and Isabelle Dugail for their critical interpretation of the results.

**Author Contributions**

Conceived and designed the experiments: C. Poitou C. Perret EN KC LT. Performed the experiments: C. Perret HD RA. Analyzed the data: FM VT YB NC VP LT. Contributed reagents/materials/analysis tools: JAW AT JLB BWP. Wrote the paper: C. Poitou EN KC LT.

**References**

1. Ng M, Fleming T, Robinson M, Thomson B, Graetz N, Margono C, et al. Global, regional, and national prevalence of overweight and obesity in children and adults during 1980–2013: a systematic analysis for the Global Burden of Disease Study 2013. Lancet 2014; 384:766–781. doi: 10.1016/S0140-6736 (14)60460-8 PMID: 24880830

2. Mraz M, Haluzik M. The role of adipose tissue immune cells in obesity and low-grade inflammation. J Endocrinol. 2014; 222:R113–127. doi: 10.1530/JOE-14-0283 PMID: 25006217

3. Ferrante AW Jr. The immune cells in adipose tissue. Diabetes Obes Metab. 2013; 15 Suppl 3:34–38. doi: 10.1111/dom.12154 PMID: 24003919

4. Schauer PR, Bhatt DL, Kirwan JP, Wolski K, Brethauer SA, Navaneethan SD, et al. Bariatric surgery versus intensive medical therapy for diabetes—3-year outcomes. N Engl J Med. 2014; 370:2002–2013. doi: 10.1056/NEJMoa1401329 PMID: 24679060

5. Zhang H, Wang Y, Zhang J, Potter BJ, Sowers JR, Zhang C. Bariatric surgery reduces visceral adipose inflammation and improves endothelial function in type 2 diabetic mice. Arterioscler Thromb Vasc Biol. 2011; 31:2063–2069. doi: 10.1161/ATVBAHA.111.225870 PMID: 21680898

6. Poitou C, Dalmas E, Renovato M, Benhamo V, Hajduch F, Abdennour M, et al. CD14dimCD16+ and CD14+CD16+ monocytes in obesity and during weight loss: relationships with fat mass and subclinical atherosclerosis. Arterioscler Thromb Vasc Biol. 2011; 31:2322–2330. doi: 10.1161/ATVBAHA.111.230979 PMID: 21799175

7. Dahlman I, Linder K, Arvidsson Nordstrom E, Andersson I, Liden J, et al. Changes in adipose tissue gene expression with energy-restricted diets in obese women. Am J Clin Nutr. 2005; 81:1275–1285. PMID: 15941876

8. Henegar C, Tordjman J, Achard V, Lacasa D, Cremer I, Guerre-Millo M, et al. Adipose tissue transcriptional signature highlights the pathological relevance of extracellular matrix in human obesity. Genome Biol. 2008; 9:R14. doi: 10.1186/gb-2008-9-1-r14 PMID: 18208606

9. Kolehmainen M, Salopuro T, Schwab US, Kekalainen J, Kallio P, Laaksonen DE, et al. Weight reduction modulates expression of genes involved in extracellular matrix and cell death: the GENOBIN study. Int J Obes (Lond). 2008; 32:292–303. PMID: 17848939

10. Bouchard L, Rabasa-Lhoret R, Faraj M, Lavoie ME, Mill J, Perusse L, et al. Differential epigenomic and transcriptomic responses in subcutaneous adipose tissue between low and high responders to caloric restriction. Am J Clin Nutr. 2010; 91:309–320. doi: 10.3945/ajcn.2009.28085 PMID: 19939982

11. Johansson LE, Danielsson AP, Parikh H, Klintenberg M, Norstrom F, Groop L, et al. Differential gene expression in adipose tissue from obese human subjects during weight loss and weight maintenance. Am J Clin Nutr. 2012; 96:196–207. doi: 10.3945/ajcn.111.020578 PMID: 22648723

12. Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet. 2011; 12:87–98. doi: 10.1038/nrg2934 PMID: 21191423

13. Sobczak A, Debowska K, Blazejczyk M, Kreutz MR, Kuznicki J, Wojda U. Calmyrin1 binds to SCG10 protein (stathmin2) to modulate neurite outgrowth. Biochim Biophys Acta. 2011; 1813:1025–1037. doi: 10.1016/j.bbamcr.2010.12.023 PMID: 21215777

14. Chiellini C, Grenningloh G, Cochet O, Scheideler M, Trajanoski Z, Ailhaud G, et al. Stathmin-like 2, a developmentally-associated neuronal marker, is expressed and modulated during osteogenesis of human mesenchymal stem cells. Biochem Biophys Res Commun. 2008; 374:64–68. doi: 10.1016/j.bbrc.2008.06.121 PMID: 18611392
15. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102:15545-15550. PMID: 16199517

16. Kirschner LS, Stratakis CA. Structure of the human ubiquitin fusion gene Uba80 (RPS27a) and one of its pseudogenes. Biochem Biophys Res Commun. 2000; 270:1106-1110. PMID: 10772958

17. Longhese MP, Anbalagan S, Martina M, Bonetti D. The role of shelterin in maintaining telomere integrity. Front Biosci (Landmark Ed). 2012; 17:1715-1728. PMID: 22201831

18. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008; 9:559. doi: 10.1186/1471-2105-9-559 PMID: 19114008

19. Nishimura S, Manabe I, Nagasaki M, Eto K, Yamashita H, Ohsugi M, et al. CD8+ effector T cells contribute to macrophage recruitment and adipose tissue inflammation in obesity. Nat Med. 2009; 15:914-920. doi: 10.1038/nm.1964 PMID: 19633658

20. Gebhardt C, Nemeth J, Angel P, Hess J. S100A8 and S100A9 in inflammation and cancer. Biochem Pharmacol. 2006; 72:1622-1631. PMID: 16846592

21. Kessel C, Holzinger D, Foell D. Phagocyte-derived S100 proteins in autoinflammation: putative role in pathogenesis and usefulness as biomarkers. Clin Immunol. 2013; 147:229-241. doi: 10.1016/j.clim.2012.11.008 PMID: 23269200

22. Panopoulos AD, Watowich SS. Granulocyte colony-stimulating factor: molecular mechanisms of action during steady state and 'emergency' hematopoiesis. Cytokine. 2008; 42:277-288. doi: 10.1016/j.cytok.2008.03.002 PMID: 18400509

23. Cimino D, Fuso L, Sfiligoi C, Biglia N, Ponzone R, Maggiorotto F, et al. Identification of new genes associated with breast cancer progression by gene expression analysis of predefined sets of neoplastic tissues. Int J Cancer. 2008; 123:1327-1338. doi: 10.1002/ijc.23660 PMID: 18561318

24. Sanchez-Palencia A, Gomez-Morales M, Gomez-Capilla JA, Pedraza V, Boyero L, Rosell R, et al. Gene expression profiling reveals novel biomarkers in nonsmall cell lung cancer. Int J Cancer. 2011; 129:355-364. doi: 10.1002/ijc.25704 PMID: 20878980

25. Wei S, Zhang L, Zhou X, Du M, Jiang Z, Hausman GJ, et al. Emerging roles of zinc finger proteins in regulating adipogenesis. Cell Mol Life Sci. 2013; 70:4569-4584. doi: 10.1007/s00018-013-1395-0 PMID: 23760207

26. Dorsett D. Cohesin, gene expression and development: lessons from Drosophila. Chromosome Res. 2009; 17:185-200. doi: 10.1007/s10577-009-9022-5 PMID: 19308700

27. Christis C, Munro S. The small G protein Arl1 directs the trans-Golgi-specific targeting of the Arf1 exchange factors BIG1 and BIG2. J Cell Biol. 2012; 196:327-335. doi: 10.1083/jcb.201107115 PMID: 22291037

28. Picard F, Gehin M, Annicotte J, Rocchi S, Champy MF, O'Malley BW, et al. SRC-1 and TIF2 control energy balance between white and brown adipose tissues. Cell. 2002; 111:931-941. PMID: 12507421

29. Jeon H, Go Y, Seo M, Lee WH, Suk K. Functional selection of phagocytosis-promoting genes: cell sorting-based selection. J Biomol Screen. 2010; 15:949-955. doi: 10.1177/1087057110376090 PMID: 20660795

30. Menssen A, Haupl T, Sittinger M, Delorme B, Charbord P, Ringe J. Differential gene expression profiling of human bone marrow-derived mesenchymal stem cells during adipogenic development. BMC Genetics. 2011; 12:461. doi: 10.1186/1471-2164-12-461 PMID: 21943323

31. Wei W, Jiao Y, Postlethwaite A, Stuart JM, Wang Y, Sun D, et al. Dual-specificity phosphatases 2: surprising positive effect at the molecular level and a potential biomarker of diseases. Genes Immun. 2013; 14:1-6. doi: 10.1038/gene.2012.54 PMID: 23190643

32. Borrego F. The CD300 molecules: an emerging family of regulators of the immune system. Blood. 2013; 121:1951-1960. doi: 10.1182/blood-2012-09-435057 PMID: 23293083

33. Rougerie P, Largeteau Q, Megrelis L, Carrette F, Lejeune T, Toffali L, et al. Fam65b is a new transcriptional target of FOXO1 that regulates RhoA signaling for T lymphocyte migration. J Immunol. 2013; 190:748-755. doi: 10.4049/jimmunol.1201174 PMID: 23241886

34. Miyashita M, Oshiumi H, Matsumoto M, Seya T. DDX60, a DEXD/H box helicase, is a novel antiviral factor promoting RIG-I-like receptor-mediated signaling. Mol Cell Biol. 2011; 31:3802-3819. doi: 10.1128/MCB.01368-10 PMID: 21791617

35. Jing Z, Yuan X, Zhang J, Huang X, Zhang Z, Liu J, et al. Chromosome 1 open reading frame 190 promotes activation of NF-kappaB canonical pathway and resistance of dendritic cells to tumor-associated inhibition in vitro. J Immunol. 2010; 185:6719-6727. doi: 10.4049/jimmunol.0903869 PMID: 21048106

36. Cildir G, Akincilar SC, Tergaonkar V. Chronic adipose tissue inflammation: all immune cells on the stage. Trends Mol Med. 2013; 19:487-500. doi: 10.1016/j.molmed.2013.05.001 PMID: 23746697
37. Hudson NJ, Dalrymple BP, Reverter A. Beyond differential expression: the quest for causal mutations and effector molecules. BMC Genomics. 2012; 13:356. PMID: 22849396

38. Speakman JR, Mitchell SE. Caloric restriction. Mol Aspects Med. 2011; 32:159-221. doi: 10.1016/j.mam.2011.07.001 PMID: 21840335

39. el Bouazzaoui F, Henneman P, Thijssen P, Visser A, Koning F, Lips MA, et al. Adipocyte telomere length associates negatively with adipocyte size, whereas adipose tissue telomere length associates negatively with the extent of fibrosis in severely obese women. Int J Obes (Lond). 2014; 38:746-749. doi: 10.1038/ijo.2013.175 PMID: 24100424

40. Zeller T, Wild P, Szymczak S, Rotival M, Schillert A, Castagne R, et al. Genetics and beyond—the transcriptome of human monocytes and disease susceptibility. PLoS One. 2010; 5:e10693. doi: 10.1371/journal.pone.0010693 PMID: 20502693

41. Wang IM, Zhang B, Yang X, Zhu J, Stepaniants S, Zhang C, et al. Systems analysis of eleven rodent disease models reveals an inflammatome signature and key drivers. Mol Syst Biol. 2012; 8:594. doi: 10.1038/msb.2012.24 PMID: 22806142

42. Honda K, Yanai H, Negishi H, Asagiri M, Sato M, Mizutani T, et al. IRF-7 is the master regulator of type-I interferon-dependent immune responses. Nature. 2005; 434:772-777. PMID: 15800576

43. Heinig M, Petretto E, Wallace C, Bottolo L, Rotival M, Lu H, et al. A trans-acting locus regulates an anti-viral expression network and type 1 diabetes risk. Nature. 2010; 467:460-464. doi: 10.1038/nature09386 PMID: 20827270

44. Wang XA, Zhang R, Zhang S, Deng S, Jiang D, Zhong J, et al. Interferon regulatory factor 7 deficiency prevents diet-induced obesity and insulin resistance. Am J Physiol Endocrinol Metab. 2013; 305: E485-495. doi: 10.1152/ajpendo.00505.2012 PMID: 23695216

45. Curtis AM, Bellet MM, Sassone-Corsi P, O'Neill LA. Circadian clock proteins and immunity. Immunity. 2014; 40:178-186. doi: 10.1016/j.immuni.2014.02.002 PMID: 24560196

46. Poitou C, Coussieu C, Rouault C, Coupaye M, Cancello R, Bedel JF, et al. Serum amyloid A: a marker of adiposity-induced low-grade inflammation but not of metabolic status. Obesity (Silver Spring). 2006; 14:309-318. PMID: 16571858

47. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011; 12:323. doi: 10.1186/1471-2105-12-323 PMID: 21816040

48. Roberts A, Pachter L. Streaming fragment assignment for real-time analysis of sequencing experiments. Nat Methods. 2013; 10:71-73. doi: 10.1038/nmeth.2251 PMID: 23160280

49. Roberts A, Trapnell C, Donaghey J, Rinn JL, Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. 2011; 12:R22. doi: 10.1186/gb-2011-12-3-r22 PMID: 21410973

50. Ronnegard L, Shen X, Alam M. hglm: A Package for Fitting Hierarchical Generalized Linear Models. The R Journal. 2010; 2:20-28.

51. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010; 11: R106. doi: 10.1186/gb-2010-11-10-r106 PMID: 20979621

52. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. 1995; Series B 57:289-300.

53. Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics. 2008; 24:719-720. PMID: 18024473

54. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009; 4:44-57. doi: 10.1038/nprot.2008.211 PMID: 19131956

55. Fuxman Bass JI, Diallo A, Nelson J, Soto JM, Myers CL, Walhout AJ. Using networks to measure similarity between genes: association index selection. Nat Methods. 2013; 10:1169-1176. doi: 10.1038/nmeth.2728 PMID: 24296474

56. Parks BW, Nam E, Org E, Kostem E, Norheim F, Hui ST, et al. Genetic control of obesity and gut microbiota composition in response to high-fat, high-sucrose diet in mice. Cell Metab. 2013; 17:141-152. doi: 10.1016/j.cmet.2012.12.007 PMID: 23312289

57. Bennett BJ, Farber CR, Orozco L, Kang HM, Ghazalpour A, Siemens N, et al. A high-resolution association mapping panel for the dissection of complex traits in mice. Genome Res. 2010; 20:281-290. doi: 10.1101/gr.099234.109 PMID: 20054062