Induction of Olfaction and Cancer-Related Genes in Mice Fed a High-Fat Diet as Assessed through the Mode-of-Action by Network Identification Analysis

The pathophysiological mechanisms underlying the development of obesity and metabolic diseases are not well understood. To gain more insight into the genetic mediators associated with the onset and progression of diet-induced obesity and metabolic diseases, we studied the molecular changes in response to a high-fat diet (HFD) by using a mode-of-action by network identification (MNI) analysis. Oligo DNA microarray analysis was performed on visceral and subcutaneous adipose tissues and muscles of male C57BL/6N mice fed a normal diet or HFD for 2, 4, 8, and 12 weeks. Each of these data was queried against the MNI algorithm, and the lists of top 5 highly ranked genes and gene ontology (GO)-annotated pathways that were significantly overrepresented among the 100 highest ranked genes at each time point in the 3 different tissues of mice fed the HFD were considered in the present study. The 40 highest ranked genes identified by MNI analysis at each time point in the different tissues of mice with diet-induced obesity were subjected to clustering based on their temporal patterns. On the basis of the above-mentioned results, we investigated the sequential induction of distinct olfactory receptors and the stimulation of cancer-related genes during the development of obesity in both adipose tissues and muscles. The top 5 genes recognized using the MNI analysis at each time point and gene cluster identified based on their temporal patterns in the peripheral tissues of mice provided novel and often surprising insights into the potential genetic mediators for obesity progression.


Introduction
Microarray analysis has enabled the use of whole-genome expression profiling to understand the mechanisms underlying obesity and metabolic complications and to identify key genetic mediators. Statistical approaches used to analyze microarray data can be classified into 2 major categories: methods that identify differentially expressed genes [1,2] and those that classify genes according to the functional dependency (e.g., hierarchical clustering) [3]. Although microarray analysis has yielded some promising results, it is not a very practical method considering the fact that identification of genes directly affected by a condition is difficult from the hundreds to thousands of genes that exhibit changes in expression. To overcome this problem, Berneardo et al. developed a model-based approach that accurately distinguishes a compound's targets from the indirect responders [4]. This approach, namely, the mode-of-action by network identification (MNI), involves the reverse engineering of a network model of regulatory interactions in an organism of interest by using a training dataset of whole-genome expression profiles. The MNI algorithm has been applied successfully to identify disease mediators as well as drug targets by studying gene-expression data from yeast [4], humans (A. Ergun and J.J. Collins, unpublished data), bacteria, and other organisms (X.H., unpublished data). Differential expression can be studied from a static or temporal viewpoint. In a static experiment, the arrays are obtained irrespective of time, essentially taking a snapshot of gene expression. On the other hand, in a temporal experiment, the arrays are collected over a time course, facilitating the study of the dynamic behavior of gene expression. Most previously obtained microarray datasets were static, that is, the results obtained on the basis of the measurement of gene expression at a single time point [5]. Since the regulation of gene expression is a dynamic process, it is important to identify and characterize the changes in gene expression over time. Therefore, numerous time-series microarray experiments have been performed to study such biological processes such as abiotic stress, disease progression, and drug responses [6][7][8].
Microarray analysis for studying the mechanisms underlying obesity was first reported by Soukas et al. in 2000 [9]. They used approximately 6,500 murine genes in pairs of adipose tissues in ob/ ob mice and wild-type lean mice. Subsequently, many such studies were conducted: more than 30 microarray approaches have been exploited in assessing the changes in gene expression in the adipose tissues, liver, hypothalamus, skeletal muscles, small intestines, and kidneys of lean and obese animals or human subjects. A frequent limitation of these studies is that they are not time-resolved and do not necessarily provide information of an end-point or disease stage. Considerably less is known about the key genetic mediators of HFD-induced obesity and the dynamics of changes in metabolic processes related to this condition. To gain more insight into the genetic mediators associated with the onset and progression of diet-induced obesity and metabolic diseases, we studied the molecular changes in response to the HFD by using an integrative time-resolved approach.

Ethics statement
All animal experiments were performed in accordance with the Korean Food and Drug Administration (KFDA) guidelines. Protocols were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of the Yonsei Laboratory Animal Research Center (YLARC) (Permit #: 2011-0061). All mice were maintained in the specific pathogen-free facility of the YLARC.

Animals and diets
Five-week-old male C57BL/6N mice were obtained from Orient Bio (Gyeonggi-do, South Korea). All animals were housed in specific pathogen-free conditions, with 2162.0uC temperature, 5065% relative humidity, and a 12 h-light/12 h-dark cycle. From a week before the diet intervention was started, all animals were fed standard chow. At the beginning of the study, mice were divided into 2 groups: (1) control group fed the normal diet (ND, n = 40) and (2) a group fed the high-fat diet (HFD, n = 40). Mice were provided food and water ad libitum. The body weight and food intake were monitored throughout the study. At 2, 4, 8, and 12 weeks after the initiation of the study, 10 animals from each group were killed. Tissues were snap-frozen immediately in liquid nitrogen and stored at 280uC until further processing.

RNA extraction for microarray analysis
Total RNA was extracted from the epididymal and subcutaneous fat tissues and gastrocnemius muscle of each mouse by using Trizol (Invitrogen, CA, USA), according to the manufacturer's recommendations. Concentrations and purity of RNA samples were determined using a Nano Drop ND-1000 spectrophotometer (Nano Drop Technologies, Inc., Wilmington, DE, USA). RNA preparations were considered suitable for array hybridization only if samples showed intact 18S and 28S rRNA bands and displayed no chromosomal peaks or RNA degradation products. The integrity of the RNA samples was determined using a Bioanalyzer 2100 System (Agilent Technologies, Palo Alto, CA, USA).

Microarray hybridization and data analysis
Equal amounts of total RNA were pooled from 10 mice in each experimental group and subjected to microarray experiments in triplicate. For analysis, 2 mg total RNA was labeled and amplified using the Universal Linkage System antisense RNA (aRNA) labeling kit (Kreatech Diagnostics, Amsterdam, The Netherlands). The Cy5-labeled aRNAs were resuspended in 10 mL of hybridization solution (GenoCheck, Korea). The labeled aRNAs were hybridized to the NimbleGen mouse whole genome 12-plex array (Roche NimbleGen, Inc., WI, USA) that contained 60-mer probes representing 42,576 genes (average 3 probes per target). The arrays were scanned using a GenePix 4000B microarray scanner. The data were extracted from the scanned images using NimbleScan software version 2.4 (Roche NimbleGen), and the Robust Multichip Average algorithm was used to generate gene expression values. The normalized and log-transformed intensity values were then analyzed using GeneSpring GX 10 (Agilent Technologies, Santa Clara, CA, USA) and GenePlex (Istech, Inc., Seoul, South Korea). The details of labeling, hybridization, scanning, and normalization of the data are provided on the NimbleGen website (http://www.nimblegen.com). Gene expres-sion levels between the ND and HFD samples were assessed by comparing the average expression ratios of each group. Hierarchical clustering was performed in GeneSpring GX 7.3.1 software (Agilent Technologies, Santa Clara, CA), using average gene expression values under HFD condition divided by the median of ND gene expression, per time-point.

MNI algorithm
We constructed a compendium dataset consisting of hundreds of expression profiles in the organism of interest; that the expression profiles were downloaded from the Gene Expression Omnibus, a public repository of microarray studies. The MNI algorithm was applied, using the method developed by Xing et al. [10], and was configured to output the top 200 mediators for each sample and generate the associated Z-scores for those probe sets. The Z-score for probe sets that were not within the list of the top 100 probe sets identified as mediators for a given sample were set to zero. To identify a characteristic list of genes within each group, the Z-scores across samples and probe sets for corresponding genes were averaged and ranked. The top 100 genes within that list were selected to be reported as significant genetic mediators. A higher average Z-score is an indication of higher number of occurrences of a gene on the lists generated by the MNI algorithm in each group. The 100 highest ranked genes were classified according to the biological process in which they are involved as per the criteria established by the GO.

Effect of HFD feeding on visceral adiposity
The body weight gains of mice fed the 2 diets over the 12-week period are shown in Figure 1A. The difference in body weight between the 2 groups continued to increase over the course of experimental feeding: the difference was about 45% by 12 weeks. The increase in body weight associated with the HFD was partially attributed to the expansion of visceral adipose tissues. The masses of the epididymal, perirenal, mesenteric, and retroperitoneal fat pads of the mice fed HFD for 12 weeks were 42%, 40%, 54%, and 42%, respectively; the difference in the masses was larger in the HFD-fed mice than in the ND-fed group ( Figure 1B-F). Moreover, HFD-fed mice exhibited significant reductions in the wet weights of the gastrocnemius (213%) and soleus (216%) muscles at 12 weeks compared with those in the ND-fed mice ( Figure 1G and H).
Transcription response of WAT and muscle to HFD during the 12-week time-course Gene expression profiling in the WAT and muscle of mice was assessed through the oligonucleotide microarray analysis. Among 25,291 genes on the NimbleGen Mouse Whole Oligo 12-plex chip used in this study, 21,890 genes (86%) were identified as known genes. After determination of the temporal effects of the HFD across 12-week time-course, we focused on dissecting the HFD specific effects on the transcriptome of epididymal and subcutaneous fats and muscle. Microarray data were analyzed by hierarchical clustering of enriched functional groups of genes (based on Gene Ontology) and the major results are graphically illustrated in a heat map ( Figure 2). The HFD elicited distinct changes in gene expression in epididymal and subcutaneous fats and muscle of mice over time, and most significant changes were shown in epididymal fat tissue. Specifically, prominent expression changes were observed at the early phase (week 2 to week 4) and the enrichment of lipid metabolism and inflammatory processes were significant among the up-regulated HFD-responsive genes, whereas G-protein coupled receptor protein signaling pathway and electron transport were most significant among the downregulated HFD-responsive genes in the epididymal fat tissue.

MNI analysis of the time course treatment with the HFD
To elucidate the time course and metabolic processes underlying obesity progression induced by the HFD, we determined the gene expression profiles of the epididymal and subcutaneous fat tissues and gastrocnemius muscle of mice by using oligonucleotide microarray analysis. Each of these data was queried against the reconstructed network (MNI algorithm), and the resulting       potential genetic mediators in each case were ranked according to the Z-score statistic. The lists of top 5 potential genetic mediators for obesity progression in the epididymal and subcutaneous fat tissues and gastrocnemius muscle of mice fed the HFD for 2, 4, 8, and 12 weeks are shown in Tables 1, 2, 3. The most characteristic genes across all tissues in the list were associated with cancer; the genes in this category included Nek11, Gli2, Tmem46, Mep1b, Ccdc109b, Rab23, Patz1, and Hdac9. The second representative functional theme was related to olfactory transduction, and these genes included Olfr1181, Olfr1173, Olfr855, Olfr1056, Olfr716, and Tmem16b. To validate the microarray results quantitatively, we analyzed the mRNA expression levels of top-ranked genes by realtime PCR. In all cases, a strong correspondence between the microarray data and the real-time PCR results was observed ( Figure 3). We also measured the basal expression levels of selected genes including several olfactory receptors in the epididymal fat tissues of ND-or HFD-fed mice, using real-time PCR. The results indicated that the basal expression levels of highly ranked olfactory genes (Olfr1181, Olfr513, Olfr960, and Olfr1245) were comparable to those of top five genes (Gli2, Gucy2c, Atp8b3, and Tmem46) identified by MNI analysis at week 4 in the epididymal adipose tissue ( Figure S1).

Functional analysis of the highly ranked genetic mediators
We next focused on the GO-annotated pathways that were significantly overrepresented among the highly ranked genetic mediators. For our analysis, we subjected the 100 highest ranked genes identified by MNI analysis in the epididymal and subcutaneous fat tissue and gastrocnemius muscle of mice with diet-induced obesity to pathway analysis based on the GO biological process annotations (Tables 4, 5, 6). We found that the olfactory transduction was highly enriched in the epididymal and subcutaneous fat tissue and gastrocnemius muscle of the HFD-fed mice compared to the ND-fed mice at all time points. Even the second representative functional theme of epididymal fat was related to cancer at all time points. The pathways thought to be associated with obesity progression in the epididymal fat as per the MNI analysis included Wnt signaling pathway, melanogenesis, chemokine signaling pathway, focal adhesion, MAPK signaling pathway, purine metabolism, regulation of actin cytoskeleton, Table 7. List of genes that exhibited decreasing ranking across time point as revealed by the MNI analysis, with a peak at 2 week in the epididymal adipose tissue of mice. neuroactive ligand-receptor interaction, and extracellular matrix (ECM)-receptor interaction. In the subcutaneous fat, other pathways identified by the MNI analysis for obesity progression included calcium signaling pathway, gonadotropin-releasing hormone (GnRH) signaling pathway, axon guidance, cell cycle, and tyrosine metabolism. In the gastrocnemius muscle, besides the olfactory transduction mentioned above, the over-represented groups identified according to GO biological processes for obesity progression were those involved in the various cellular processes such as neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, pathways associated with cancer, insulin signaling pathway, pathways associated with colorectal cancer, adipocytokine signaling pathway, type II diabetes mellitus, and cell adhesion molecules.

Representative time-course profile clusters
We subjected the 40 highest ranked genes identified by MNI analysis at each time point in the epididymal and subcutaneous fat tissues and gastrocnemius muscle of mice with diet-induced obesity to clustering based on their temporal pattern. Figure 4 shows genes that were observed to have decreasing ranking across time point, with a peak at 2 week. Biological processes controlled by genes in this cluster included regulation of insulin resistance ( (Table 19).

Discussion
Animals use their olfactory system to monitor the chemical environment for molecules that reveal food sources or toxic substances and signal the presence of predators [11]. There are numerous olfactory receptors of different types, with as many as 1,000 in the mammalian genome that represent approximately 3% of the human genomeentire genetic information [12]. The information related to odor gathered by these olfactory receptors is funneled through a common signaling pathway. When an olfactory receptor binds to its odorant, it activates a single species of G protein, the olfactory trimeric G protein (Golf), which in turn activates the olfactory isoform of adenylate cyclase (AC3) [12]. Converging evidence has demonstrated that the olfactory system is a target for hormones related to metabolism and food-intake regulation; it adapts its function to nutritional needs by promoting or inhibiting food foraging [13]. Recent studies have found that obese patients display decreased olfactory acuity [14] and are significantly more likely to have absolute olfactory dysfunction or anosmia [15]. Furthermore, Simchen et al. showed that the abilities to detect and identify odors have been found to decrease as body mass index (BMI) increases in subjects less than 65 years old, independent of any linkage to food odor or gender [16].
Recently, the elements of olfactory-like chemosensory signaling have been found to also present in nonolfactory tissues such as testis [17], brain [18], and heart [19]. To our knowledge, this is the first study that shows a differential mRNA expression and high MNI ranking of olfactory receptors in the epididymal and subcutaneous fat tissues and muscles between ND and HFD-fed mice (Tables 1, 2, 3). These results imply that the olfactory receptors and the molecules involved in olfactory transduction might be the mediators of HFD-induced obesity progression in the peripheral tissues. This hypothesis is supported by the fact that the increased cAMP production by AC3 activates cAMP responsive element binding (CREB) protein, leading to increased adipogenesis in an obese mouse model. Furthermore, mice lacking AC3, which is a downstream regulator of olfactory receptors, exhibit obesity that is apparently caused by low locomotor activity, hyperphagia, and leptin insensitivity [20]. In future studies, it will be intriguing to further investigate the role of individual olfactory receptors in peripheral tissues, such as the pancreas, liver, muscle, and fat, to better understand the activation process of these signaling pathways and their physiological roles. Cancer-related genes such as Nek11, A4gnt, Srp9, Gli2, Gucy2c, Lsm1, Duoxa1, Lasp1, Ret, Bex2, Vav3, Kcnrg, Tle6, Rab23, Dcc, Rassf2, Perp, Pdgfr1, Lin28, Gstm1, Safb2, Tmem46, and Hdac9 were remarkably overrepresented in time-course clusters identified by the MNI analysis in the epididymal and subcutaneous fat tissues and gastrocnemius muscle of mice with diet-induced obesity (Figs. 4,5,6,7,8). Of these 23 genes, 6 are known breast cancerrelated genes (Lsm1, Duoxa1, Ret, Bex2, Rassf2, and Safb2). Lsm1 is a transforming oncogene that is amplified and overexpressed in breast cancer [21] and might affect either cell cycle progression or apoptosis [22]. Duoxa1, which was originally identified as a numbinteracting protein, was recently shown to function as a maturation factor in breast cancer [23]. Ret exhibits both estrogenand retinoic acid-dependent transcriptional modulation in breast cancer [24]. Bex2 has a significant role in promoting cell survival and growth in breast cancer cells [25,26], and Rassf2 might function as a tumor suppressor gene in in vitro cell migration and cell cycle progression [27]. The expression of Safb2 protein, which functions as estrogen receptor co-repressor and growth inhibitor, was lost in approximately 20% of breast cancers [28]. Many studies have attempted to determine the relationship between diet and breast cancer. Dietary fat is a source of endogenous estrogen and has been suggested as a possible risk factor for breast cancer Table 9. List of genes that exhibited decreasing ranking across time point as revealed by the MNI analysis, with a peak at 2 week in the muscle of mice.  [29]. To our knowledge, this is the first study showing an association between these 6 genes involved in breast cancer development and HFD-induced obesity in a rodent model. Colon cancer-related genes such as Nek11, Gucy2c, Srp9, Tle6, and Pdgfrl were also overrepresented in the time-course clusters identified by the MNI analysis in the epididymal and subcutaneous fat tissues and gastrocnemius muscle of mice with diet-induced obesity (Figs. 4, 5, 6, 7, 8). Nek11, a member of the NIMA-related kinase family, phosphorylates Cdc25a and controls its degradation; Cdc25a phosphorylation is required for cell cycle progression in colorectal cancer cells [30]. Gucy2c and Srp9 have been shown to be overexpressed in colorectal cancer cells and were recently shown to function as a candidate biomarker for colon cancer [31,32]. Tle6 is recurrently overexpressed in human colon cancer and enhances cell proliferation, colony formation, migration, and xenograft tumorigenicity [33]. Pdgfrl acts as a tumor suppressor and inhibits the growth of colorectal cancer cells [34]. Epidemiological studies indicate that both high body weight and high body mass index (BMI) were significantly associated with an increased colon cancer risk. Intra-abdominal visceral obesity, high plasma glucose levels, HbA1C, and C-peptide were also found to be associated with increased risk of colorectal cancer [35][36][37][38]. The current study showed that the above-mentioned genes that are involved in the regulation of colon cancer might play a genetic role in the development of obesity. No mechanistic insights have been reported to explain the relationship between the regulation of cancer-related genes in the adipose tissue or muscle and cancer susceptibility. It could be probable that the changes in the expression of cancer-related genes in the adipose tissue may accompany the regulation of same genes in epithelial tissues such as breast or colon.
Genes that were found to have the highest rank at the early phase and return to baseline after several weeks might be considered genetic mediators of acute-phase response in metabolic processes related to HFD-induced obesity. Dusp12 was one of the 58 genes that were observed to have decreasing ranking during the development of obesity, with a peak at 2 week. Previous studies identified several single nucleotide polymorphisms in this gene associated with type 2 diabetes in different populations, including Caucasians and Chinese [39]. Dusp12 is a glucokinase-associated protein that participates in glycolysis in the liver and dephosphorylation of cytoplasmic glucokinase in the pancreatic beta cells [40]. Therefore, Dusp12 might play a role in the regulation of glycolysis during the early stages of obesity. When glycolysis was decreased, whole-body glucose disposal was also reduced, indicating a decrease in glucose utilization in the peripheral tissues in response to the HFD. The latter likely results from an impaired glucose transport that precedes impaired insulin signaling.
Def6 and Mapk9 were one of the 145 genes that were found to have the highest rank at the intermediate time points of 4 or 8 weeks during the development of obesity. Def6, a novel type of activator for Rho GTPase, is expressed in myeloid cells, and disruption of Def6 expression leads to defects in toll-like receptor 4 (TLR4) signaling and innate immune responses [41]. Rho GTPases have been shown to be recruited to the cytosolic domain of TLR and the closely related interleukin 1 receptor (IL-1R) and to regulate the production of proinflammatory cytokines [42,43]. In the present study, the high MNI ranking of Def6 in the subcutaneous adipose tissue of HFD-fed mice suggested that it might participate in the regulation of obesity-induced inflammation through TLR4 signaling. Mapk9, which is ubiquitously expressed, can invoke transcription factors such as c-Jun and many other apoptosis-related proteins [44]. Interestingly, recent studies have shown that the knockdown of Mapk9 leads to reduced serum levels of glucose, insulin, and homeostatic model assessment and therefore reverses insulin resistance in HFD-fed mice [45]. These findings provide supporting evidences to the high MNI ranking of Mapk9 associated with HFD-induced obesity observed in the present study. However, further studies are required to elucidate the precise function of Mapk9 in the development of HFD-induced type 2 diabetes.
Smad7, Adhfe1, and Pyy are one of the 65 genes that showed increasing ranking during the development of obesity, with a peak at 12 weeks. Smad7 was initially characterized as a factor induced by shear stress in vascular endothelial cells [46]. Only recently, new functions of Smad7 were elucidated: it inhibits transforming growth factor-b (TGF-b)-activated responses [46]. TGF-b is known to inhibit adipose differentiation of preadipocyte cell lines and primary cultures [47] and to block adipogenesis in vivo [48]. This suggests that Smad7 enhances adipogenesis through the inhibition of TGF-b signaling. Adhfe1 was characterized as a hydroxyacid-oxoacid transhydrogenase that catalyzes the conversion of c-hydroxybutyrate to succinic semialdehyde [49]. Recently, Adhfe1 was suggested to play a role in adipocyte differentiation. The expression of Adhfe1 transcript is tightly linked to the phenotype of mature adipocytes both in vivo and in vitro, although the mechanisms underlying Adhfe1-mediated regulation of adipogenesis remain poorly understood [50]. Pyy, which is expressed and secreted in endocrine intestinal cells, plays a role in reducing appetite and caloric intake [51]. Recently, plasma Pyy concentra-tions were found to be decreased in both obese humans [52] and diet-induced obese mice [53]. These studies might suggest that Smad7 and Adhfe1 play a role in obesity by amplifying the aggressive effect of adipogenesis.
Camk2g and Tmem67 are one of the 8 genes that exhibited a constant high MNI ranking from 2 to 12 weeks. The increase of cytosolic Ca 2+ in the beta cells is central to the initiation of insulin secretion under physiological conditions [54]. Recent findings suggest that Camk2g involved in the regulation of calcium in the islet beta cells is a candidate gene for type 2 diabetes [55]. The Tmem67 gene mediates a fundamental developmental stage of ciliary formation and epithelial morphogenesis [56]. In addition, defects in the Tmem67 gene resulted in Meckel syndrome type 3, Joubert syndrome type 6, and nephronophthisis 11, which show many clinical phenotypic similarities, including hepatic fibrosis [56,57]. Consumption of fat-rich diets seems to play an important role in the pathogenesis of hepatic steatosis and its progression to fibrosis [58]. The constantly high MNI ranking of Tmem67 from 2 to 12 weeks associated with a HFD suggests that Tmem67 might participate in the development of hepatic fibrosis.
In summary, this study is the most comprehensive investigation of the gene expression patterns conducted using a time-resolved approach to gain insight into the development of HFD-induced obesity in a mouse model. A reverse-engineered gene network was used for the first time for the identification of key genetic mediators and pathways that have been implicated in the initiation and advancement of obesity. We highlighted the sequential induction of distinct olfactory receptors and stimulation of cancer-related genes during the development of obesity. To our knowledge, the proposed changes in the olfactory transduction machinery as per the MNI ranking have not been previously reported. These putative mechanisms clearly need further investigation. The top 5 genes recognized through the MNI analysis at each time points (2,4,8, and 12 weeks) and gene clusters identified based on their temporal patterns in the 3 different tissues (visceral and subcutaneous adipose tissues and muscle) of mice need special attention as potential genetic mediators for obesity progression. Figure S1 The basal expression levels of some target genes identified by MNI analysis. Quantitative real-time PCR analysis of the basal expression on highly ranked olfactory genes and top 5 genes at week 4 in the epididymal adipose tissues