Association of Rev-erbα in adipose tissues with Type 2 diabetes mellitus amelioration after gastric bypass surgery in Goto-Kakizaki rats

Rui Zhang, Caifeng Yan, Xinrong Zhou, Bangguo Qian, Fuqiang Li, Yidan Sun, Chen Shi, Bing Li, Shigeru Saito, Katsuhisa Horimoto, Huarong Zhou

Abstract

We estimated the key molecules related to Type 2 diabetes mellitus (T2DM) in adipose, liver, and muscle tissues, from nonobese diabetic Goto-Kakizaki (GK) rats and their Wistar controls, by computationally analyzing the expression profiles in open source data. With the aid of information from previous reports, Rev-erbα in adipose tissue emerged as one of the most plausible candidates. Here, in animal models, including GK rats surgically treated to ameliorate T2DM, we examined the association of Rev-erbα in adipose tissue with T2DM progression. After analyses of the Rev-erbα mRNA expression in the adipose tissue of our animal models, we compared the Rev-erbα protein expression levels in the adipose, liver, and muscle tissues of GK and Wistar controls at the ages of 1 mo (M), 3M, and 6M. The Rev-erbα protein levels in adipose tissue showed a distinctive pattern, with the negative correlation of an increasing trend in GK rats, and a decreasing trend in Wistar rats during aging, from those in liver and muscle tissues. Moreover, dysregulation of the circadian Rev-erbα expression in the adipose tissue of 6-mo-old GK rats was also observed. In particular, we ameliorated T2DM in GK rats by gastric bypass surgery, and revealed that T2DM amelioration in diabetic GK rats was associated with improved circadian Rev-erbα expression, in a comparison between the surgically treated and untreated GK rats. The roles of Rev-erbα in adipose tissue were further investigated by observations of Rev-erbα-related molecules, with reference to previous reports.

  • Type 2 diabetes mellitus
  • Rev-erbα
  • adipose tissue
  • network analysis
  • gastric bypass surgery

in the past two decades, the worldwide incidence of diabetes has dramatically increased (20). As one of the most complex human diseases, diabetes has its roots in genetics, but it is also heavily influenced by the environment (32). It is generally accepted that Type 2 diabetes mellitus (T2DM), the majority of all cases, has two basic metabolic defects: insulin resistance and insulin-secreting β-cell defects. Some controversy still exists as to whether insulin resistance or inadequate insulin secretion occurs first in the pathogenesis of T2DM (25). In the transition stage, pancreatic β-cells are able to compensate for secretion defects and/or insulin resistance by increasing insulin production and secretion, to maintain normal or near-normal glucose levels. However, hyperinsulinemia per se can increase insulin resistance, which requires more insulin secretion from β-cells. As insulin resistance worsens or more defects in “exhausting” β-cell insulin secretion occur, together, these abnormalities lead to impaired glucose tolerance and further impaired fasting hyperglycemia, which are prediabetic phenomena. The metabolic sequences that eventually lead to the diagnosis of T2DM, which is called the time of “no return,” occur over years (3, 14, 35).

Significant progress related to diabetes at the molecular level has been made, due to intensive research during the last few decades. Indeed, the participation of some important molecules has been confirmed by in vivo gene modifications in animal models and clinical studies via integrated physiology. Furthermore, the translation of the research at the molecular, cellular, and physiological levels into practice has achieved better control of blood glucose, and thus reduced mortality in people with diabetes (16). However, despite the wealth of knowledge about diabetes, the incidence of T2DM has been continuously increasing for the last two decades (20).

Because of the diverse abnormalities in multiple genes from different tissues, it is generally difficult to identify the genes that are vital to the development and progression of T2DM. Here, to overcome this problem, we integrated the previous information with high throughput data, by a systematic combination of computational analyses and experiments. We applied our approach to an animal model, Goto-Kakizaki (GK) rats, for evaluating phenomena at the molecular and cellular levels in the whole body. The GK rat is a nonobese spontaneous diabetic rodent model with congenital β-cell defects, and insulin resistance has been established in GK rats at 8 wk, under both basal state and clamp conditions (7). Since GK rats display the spontaneous onset of T2DM and exhibit the two basic defects of human T2DM (19), we used three groups of GK rats, ages 1 mo (1M), 3M, and 6M, representing the early, middle, and stable stages of T2DM, respectively, in this study (18).

Here, we report that the Rev-erbα protein levels in adipose tissue are associated with T2DM in the GK animal model. First, we applied our computational procedure (33) to the expression profiles from three tissues—adipose, liver, and muscle—of GK and Wistar rats at 4 wk in open source data (2, 28, and 51), and we found one candidate for a regulatory relationship, retinoic acid receptor-related orphan receptor alpha (RORA)-Rev-erbα, in adipose tissue. We then proceeded to experimentally examine the plausibility of the Rev-erbα association with T2DM in our animal models. An examination of the Rev-erbα mRNA levels in our animal models revealed a distinct Rev-erbα protein expression pattern in adipose tissue from those in liver and muscle, in comparisons between GK rats and Wistar controls with aging. Moreover, dysregulated circadian Rev-erbα expression in adipose tissues from 6M GK rats was also observed. In particular, we found a negative correlation between the Rev-erbα protein levels in adipose tissue and the degree of T2DM in GK rats after gastric bypass surgery (42). Interestingly, the dysregulated circadian Rev-erbα protein expression in GK rats was partially reversed after bariatric surgery. In addition, the roles of Rev-erbα in the molecular and cellular mechanisms underlying T2DM progression were further surveyed by observations of Rev-erbα-related molecules, with reference to previous knowledge about T2DM.

MATERIALS AND METHODS

Data Analyzed for Selecting the Candidate Genes Related to T2DM in GK Rats

We analyzed three sets of gene expression data, measured in adipose (51), liver (2), and muscle (28) tissues from GK and Wistar rats of a Japanese colony. In terms of the circadian rhythm and the data collection, the animals were killed between 1.5 and 3.5 h after the lights were turned on (2, 28, and 51). The data were cited from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/projects/geo/) database (GSE13269, GSE13271, and GSE13270). The three data sets comprised 31,099 probes, measured using the Affymetrix Microarray Suite 5.0 (Affymetrix), and they were further reduced into 14,506 genes, for five samples of male spontaneously diabetic GK rats and Wistar controls at each of five time points (4, 8, 12, 16, and 20 wk of age). GK rats younger than 3 wk of age are considered prediabetic. Hyperglycemia is detected after 4 wk, with some variations between 5–14 wk of age (18). T2DM is stably established in animals older than 18 wk of age. In this analysis, the five analysis periods were reclassified into three periods: 4-wk, 8–12-wk, and 16–20-wk periods. In this reclassification, the periods are composed of 5, 10, and 10 samples, representing the early, middle, and stable stages of diabetes, respectively. Indeed, this enhanced breakdown was employed because diabetes progression is generally discussed from simply the dichotomy of nondiabetes and diabetes, while only considering the dynamic ranges of gene expression measurements and small numbers of animal samples per period, as reported in a previous diabetes study (56).

Overview of Computational Procedures to Narrow the Key Molecules

The computational procedure to identify the candidates of key molecules responsible for phenotype changes was described previously (33) (also see appendix). The overview is schematically shown in Fig. 1.

Fig. 1.

Overview of the computational procedure. Schematic description of the computational procedure to narrow down the key molecules from the gene expression data. First, we estimate the network structures by two methods: network screening (NS) and modified path consistency algorithm (mPCA). The binary structures from each of the network structures are then selected on the basis of whether the genes regulated by each transcriptional factors (TF) are included in the gene expression signature and whether the set of genes regulated by each TF is functionally enriched for the known gene set database (MSigDB) (45). The binary relationships that are commonly identified by the two methods are finally selected.

First, we performed network analyses by two methods. One is a knowledge-based approach (named “network screening”; NS), which constructs functional networks and then estimates the consistency of the network structures with the measured data (39). This is a powerful method to detect significantly activated functional networks in a particular phenotype, but it depends on previous knowledge. The other is a data-driven inference approach, which estimates the conditional independency between the genes by calculating the partial correlation coefficients. The calculation algorithm was previously developed (44), and it was modified to analyze redundant biological data (named “modified path consistency algorithm”; mPCA) (33, 40). Second, we used the two methods to extract each set of binary relationships of transcriptional factors (TFs) and their regulated genes among the detected network structures. Third, we further selected the binary relationships from the two sets, under the conditions that the regulated genes were included in the gene expression signature obtained by the t-test with 5% significant probability, and that the set of genes in the expression signature, which were regulated by one TF, was enriched for each of the 871 gene sets in the C2 collection of MSigDB (45) with 5% probability significance. Thus, the binary relationships were selected from the original networks in terms of two aspects, the difference in gene expression and the functional enrichment in known gene sets. Finally, we defined the TFs as the key molecules that were commonly included in the two sets of binary relationships. Note that once the key TFs were selected, the genes regulated by the key TFs were automatically selected by the above procedure.

Experiments

Animals.

Wistar and GK rats from a Japanese colony were obtained from the National Rodent Laboratory Animal Center, Shanghai. Animals were housed in a temperature- and humidity-controlled facility with strict light control. The lights were turned on at 7 AM [zeitgeber time (ZT) 0] and were turned off at 7 PM (ZT12). To match the time of tissue harvest with the microarray analysis, unless otherwise specified, the animals in this report were killed between 9 and 10 AM, which corresponds to ZT2. However, in the experiments measuring the circadian Rev-erbα expression, the groups of Wistar and GK animals were killed at ZT2, ZT7, ZT10, and ZT12. The 1M-old GK rats received an intraperitoneal injection of 6-bromoindirubin-3′-oxime (BIO) (100 mg/kg ip; Sigma, St. Louis, MO) twice per week for 2M and were killed at 3M of age. All experiments were approved by the Institutional Animal Care and Use Committee at Shanghai Institutes of Biological Science, China Academy of Science, in compliance with the National Institutes of Health guidelines and the Guide for the Care and Use of Laboratory Animals.

Surgical techniques.

Operations were performed on 3M-old GK rats, which were killed at 6M of age. After overnight fasting and under continuous isoflurane anesthesia, the abdomen was opened by midline incision. The stomach was mobilized by opening the lesser sac through the greater omentum. The vagus nerve and the accompanying right and left gastric vasculature were gently dissected from the lesser curve, and the body of the stomach was divided between two hemostats. A biliopancreatic limb extending 16 cm from the ligament of Trietz was transected. The distal segment was anastomosed to the gastric remnant, and the proximal segment was drained into 30 cm of the “Roux” limb, using side-to-side anastomosis.

Postoperatively, the animals recovered in a warm box before returning to the animal facility. All animals were provided buprenorphine and ampicillin for 3 days. The rats in the surgical groups were only allowed to drink purified water for 12 h after surgery. A liquid diet containing 5% glucose and 0.2% KCl was provided for the next 48 h. The rats were slowly accustomed to ingesting normal laboratory chow thereafter. Normal laboratory chow was provided to all experimental groups.

Insulin tolerance tests.

Insulin tolerance tests were performed in 16-h fasted rats at 11 wk postoperation. The rats received 1 UI/kg NovoLog regular insulin, injected intraperitoneally. Blood glucose levels were measured at 0, 15, 30, 60, and 120 min via the tail vein, using a Roche Accu-Chek blood glucose meter.

Oral glucose tolerance tests.

Oral glucose tolerance tests were measured at 12 wk postoperation. Tests were performed in 16-h fasted rats. The rats were given 2 g/kg d-glucose orally, and glucose was measured at 0, 15, 30, 60, 90, and 120 min from tail vein blood. The area under the curve (AUC) was calculated for estimating the glucose tolerance.

Histochemistry.

Adipose tissues were fixed in 4% paraformaldehyde and paraffin-embedded. Sections were subjected to standard hematoxylin-and-eosin staining.

Tissue extraction.

The epididymal fat pad, liver, and the lateral surface of the hind leg skeletal muscle were removed from each rat at 1M, 3M, and 6M, and frozen at −80°C. On the day of the analysis, the same amounts of frozen tissues were cut and immediately homogenized in ice-cold extraction buffer with complete proteinase inhibitor mixture (Sigma). Triton X-100 (1% final concentration) was added to the homogenates, which were centrifuged at 1,000 g at 4°C, resulting in the separation of a fat layer. The supernatants, thus, generated were filtered to remove any residual fat, and the filtrates were centrifuged at 10,000 g. The resulting supernatant solutions were used for Western blot analyses.

Western blot analyses.

The protein concentration was determined with a DC protein assay kit (Bio-Rad Laboratories, Hercules, CA). Aliquots containing equal amounts of protein (40 μg) were subjected to SDS-PAGE and transferred to PVDF membranes. To compare the 1M, 3M, and 6M samples, we excised the bands corresponding to the same protein from different SDS-PAGE gels, and transferred them onto one PVDF membrane. After transfer, the membrane was blocked by 5% BSA for 1 h at room temperature and incubated with Rev-erbα (1:1,000; Sigma), β-actin (Cell Signaling Technology, Danvers, MA) and tubulin (Cell Signaling) antibodies overnight at 4°C, with gentle shaking. The membrane was then washed and incubated with the secondary antibody for 1 h at room temperature. The target protein band was detected with ECL by a FUJI LAS4000 bioimager.

RNA isolation and quantitative RT-PCR.

Total RNA was extracted from the frozen adipose tissue by the TRIzol reagent (Invitrogen, Carlsbad, CA), according to the manufacturer's protocol. cDNA was prepared by reverse transcription of 1 μg RNA, using an iScript cDNA Synthesis Kit (Bio-Rad, Tokyo, Japan), and real-time PCR was performed using the SYBR Green real-time PCR master mix (Toyobo, Osaka, Japan) and a StepOne PCR machine (Applied Biosystems, Foster City, CA). The following primers were synthesized by Biosune (Shanghai, China): For rat Rev-erbα (forward, 5′-AGAATGTTCTGCTGGCATGTCC-3′ and reverse, 5′-CTCGGAATGCATGTTG TTCAGG-3′); GAPDH (forward, 5′-ACAGCAACAGGGTGGTGGAC3-′ and reverse, 5′-TTT GAGGGTGCAGCGAACTT-3′); EF1A (forward, 5′-TGGTGGTTACCTTTGCTC3′ and reverse, 5′-ACTGATCTGGCCTGGATG-3′). All real-time PCRs were performed under the following conditions: 95°C for 5 min and then 40 cycles of 15 s at 95°C, 30 s at 60°C, and 30 s at 72°C. To accurately normalize real-time quantitative RT-PCR data, we chose both GAPDH and EF1A as internal controls and estimated the relative gene expression level for that sample by the geometric mean of the two control genes (47). Separate analyses in which relative expression values were normalized with the relative expression values of each control gene yielded similar results.

Statistical Analysis of Measured Data

Data are expressed as means ± SE. Significant differences between groups were determined using ANOVA followed by the Scheffé post hoc test. A P value of less than 0.05 was considered significant.

RESULTS

Candidates of Key TF Regulatory Relationships for T2DM in Three Tissues

According to the procedure outlined in Fig. 1, we detected the candidates of TFs and their regulated genes that were activated at an early stage of diabetes in adipose, liver, and muscle tissues. The distributions of the candidates of TFs and their regulated genes in the three tissues, determined by the two methods, are depicted in Fig. 2. Overall, the number of candidates detected by NS was smaller than that by mPCA. This was expected, since NS extracted the TF candidates among the restricted numbers of known regulatory relationships, while mPCA extracted the candidates from the mathematically inferred network structures, without any restrictions.

Fig. 2.

Venn diagrams of TFs detected by NS and mPCA by using the adipose, liver, and muscle tissue expression data in the 4-wk period. The numbers of detected TFs are indicated in the figure, and the total numbers of each tissue are indicated in parentheses.

Among the candidates of TFs and their regulated gene pairs, we searched for the TFs that were specifically detected in each tissue and that were commonly detected by the two methods. As shown in Fig. 2, 1 TF by NS and 23 TFs by mPCA were specifically detected in adipose tissue, and 1 TF by NS, RORA (RAR-related orphan receptor-α), was also found among the 23 TFs detected by mPCA. In liver, 22 TFs by NS and 30 TFs by mPCA were specifically detected, and then 1 TF, ELK3 (ETS-like transcription factor 3), was found commonly by the two methods. Since no TFs were specifically detected by NS in muscle, no candidates were listed. As a result, two TFs—RORA in adipose tissue and ELK3 in liver—were found. We analyzed the top 2 TF candidates to identify the candidate with the highest potential. For this purpose, we referred to the significant appearance of RORA and ELK3 by the above two methods at the middle and stable stages. According to this TF selection rule, RORA appeared in adipose tissue at 4 wk, in liver at 8 wk and 12 wk, and in muscle at 16 wk and 20 wk, and ELK3 appeared in liver at 4 wk and in muscle at 16 wk and 20 wk, after disappearing from the three tissues at 8 wk and 12 wk (data not shown). Although it is difficult to judge the influence of genes on diabetes progression by only the pattern of the network analyses, the expression of RORA seemed to be clearer than that of ELK3 during diabetes progression. Therefore, we decided to focus on RORA to further investigate its role in diabetes progression.

Next, we proceeded to narrow down the candidates of the genes regulated by RORA. Our procedure for the consideration of the expression signature and the enrichment of gene sets automatically detected only two candidates of RORA-regulated genes in adipose tissue, Rev-erbα and PRKCQ [calcium-independent, phospholipid- and diacylglycerol-dependent serine/threonine-protein kinase] (Table 1). As reported previously (4, 12, 48), both genes are closely related to diabetes, suggesting the promising performance of our computational procedure. Among the two RORA-regulated genes, Rev-erbα has been validated as an antiobesity drug target (43), but more information is needed to clarify the association of Rev-erbα in adipose tissue with diabetes progression in vivo. In contrast, the association of PRKCQ with insulin resistance is well known and has been further confirmed in a human study, as an SNP of PRKCQ on 10p15 was reportedly independently associated with Type 1 diabetes (10). Thus, we selected Rev-erbα as the candidate for further investigations of its association with T2DM, by molecular and physiological experimental approaches in GK animal models.

View this table:
Table 1.

Candidate genes regulated by RORA in adipose at 1M, identified by the two methods

Examination of the Rev-erbα mRNA Levels in Adipose Tissue in Animal Models

To test whether the inference of Rev-erbα from the open source is valid in our animal models, we first compared the Rev-erbα mRNA expression levels in adipose tissues from 1M-, 3M-, and 6M-old GK and Wistar rats in our models (Fig. 3). The only difference in Rev-erbα expression was observed at 1M between GK and Wistar rats. Interestingly, in GK rats, the difference in the Rev-erbα expression levels steadily increased along with diabetes progression at 1M, 3M, and 6M, while in Wistar rats, the difference was observed between 1M and 3M, but not between 3M and 6M. In other words, the Rev-erbα mRNA expression level in adipose tissue from GK rats was different throughout the early, middle, and stable stages of diabetes, and the difference between GK and Wistar rats was observed at 1M. In summary, the Rev-erbα mRNA expression levels at 3M and 6M were not different between GK and Wistar rats, but the association of Rev-erbα with diabetes progression was suggested from the expression changes accompanying aging in our animal model. In the following, we will focus on the role of Rev-erbα in adipose tissue during diabetes progression in our model, from the viewpoint of its physiological levels, compared with those in liver and muscle.

Fig. 3.

Rev-erbα mRNA levels in adipose tissues of Wistar and Goto-Kakizaki (GK) rats. The amount of Rev-erbα mRNA was significantly increased in the adipose tissue of 1-mo-old (1M) GK rats, compared with 1-mo-old Wistar rats. The amounts of Rev-erbα mRNA were significantly increased in the adipose tissues of 3- and 6-mo-old (3M, 6M) GK rats, compared with the 1-mo-old GK rats. Each bar represents the means ± SE (n = 6 or 7). Significance is indicated as follows: **P < 0.01; and *P < 0.05.

Rev-erbα Protein Levels with Animal Aging

We measured the Rev-erbα protein levels in adipose, liver, and muscle tissues in GK and Wistar rats during aging along with diabetes progression (early diabetes stage of 1M, middle stage of 3M, and stable stage of 6M) (Fig. 4, A–C). The protein levels in GK rats were distinctive between the three tissues along the diabetic states, compared with those of the Wistar rats (see also Table 2).

Fig. 4.

Rev-erbα protein expression levels in adipose, liver, and muscle tissues of Wistar and GK rats during animal aging. The ZT2 Rev-erbα protein expression levels at 1M, 3M, and 6M of age in the three tissues of Wistar and GK rats were measured by Western blots. Samples of the same tissue from different aged animal groups were transferred onto one PVDF membrane and stained with antibodies specific to the target proteins. The protein levels of Rev-erbα were normalized by those of β-actin for adipose and liver tissues and by tubulin for muscle tissue. A: protein levels in adipose tissue with Western blot results. B: protein levels in liver. C: protein levels in muscle. Each bar represents the means ± SE (n = 6 or 7). Significance is presented in Table 2.

View this table:
Table 2.

Relative Rev-erbα protein levels with statistical significance, identified by comparisons between and within Wistar and GK rats

As shown in Fig. 4A, the Rev-erbα protein levels in adipose tissue from the diabetes-prone GK rats showed a trend to increase as the animals aged, while in contrast, those from the Wistar rats showed a declining trend. In Table 2, the differences in the protein levels of Rev-erbα normalized by those of β-actin or tubulin, were statistically tested between the three diabetic states. Indeed, the Rev-erbα protein level in GK rats significantly increased between 3M and 6M, while that in Wistar rats significantly decreased between 1M and 3M. In addition, a significant difference between GK and Wistar rats was detected at 6M, as a result of the increment in the Rev-erbα protein level in GK rats and its decrement in Wistar rats during aging. Interestingly, in the 1M-old GK and Wistar rats, the protein levels in Fig. 4A did not reflect the significant difference in mRNA levels in Fig. 3. This reminded us of the report that during adipogenesis, Rev-erbα gene expression initially declines and subsequently increases, while the Rev-erbα protein levels are nearly the opposite, increasing early in adipogenesis and then markedly decreasing during the 3T3-L1 preadipocyte differentiation process (48). At any rate, the Rev-erbα protein levels in adipose tissue showed a distinctive pattern during aging, in the comparison of GK and Wistar rats, suggesting the possible association of Rev-erbα with diabetes in GK rats.

The Rev-erbα protein levels in liver did not significantly change in GK rats but showed gradual increments in Wistar rats with aging (Fig. 4B). A significant difference in the Rev-erbα protein levels between GK and Wistar rats was found at 6M, and that within the Wistar rat group, Rev-erbα was between 3M and 6M, due to the increment in Wistar rats (Fig. 4B and Table 2). The importance of the increment of Rev-erbα protein levels in Wistar rats and its absence in GK rats in liver at 6M, in terms of diabetes progression, seems to agree with the fact that Rev-erbα is a strong negative regulator of Bmal1 in the hypoglycemia phenomenon in liver-specific Bmal1 knockout mice (22).

In contrast to the distinctive behaviors of the Rev-erbα protein levels in adipose and liver tissues between GK and Wistar rats, relatively static behaviors were observed in muscle tissue (Fig. 4C). Indeed, in muscle, the Rev-erbα protein levels gradually increased with no significant differences between GK and Wistar rats during aging (Table 2). The behaviors implied that the Rev-erbα protein levels in muscle have no positive association with diabetes development and progression in GK rats before 6M of age.

Rev-erbα Protein Levels Along Circadian Time Points

We further investigated the distinctive Rev-erbα protein levels in GK rats by focusing on the circadian pattern of Rev-erbα expression in adipose tissues at 6M, which is the age when significant differences in the Rev-erbα protein levels were detected between GK and Wistar rats (Fig. 4A). For this purpose, we compared the mRNA and protein levels of Rev-erbα at four circadian time points between 6M-old GK and Wistar rats (Fig. 5).

Fig. 5.

Dysregulation of the circadian rhythm of Rev-erbα in adipose tissues from diabetic GK rats. 6M-old GK and Wistar rats were killed at ZT2, ZT7, ZT10, and ZT12 (n = 3 per group at each time point). RT-PCR and Western blot analysis were performed to investigate the expression of circadian Rev-erbα mRNA and the protein levels in adipose tissues. A: relative Rev-erbα mRNA levels were normalized to the geometric averaging of two internal control genes EF1A and GAPDH. B: Rev-erbα protein levels were normalized to those of β-actin. Significance within groups is indicated as follows: **P < 0.01; *P < 0.05. Significance between groups is not shown. The Rev-erbα protein levels in Wistar rats at ZT2 were significantly lower than those at any time points in GK rats, while the Rev-erbα protein levels in GK rats at ZT2 are significantly higher than those at any time points in Wistar rats.

The Rev-erbα mRNA expression levels at all of the time points were not different between the GK and Wistar rats (Fig. 5A). When ZT2 was compared with ZT12, however, the Rev-erbα mRNA levels were significantly increased in the adipose tissues of Wistar rats, and this finding was quite consistent with previous studies of liver tissues in normal mice (46). In GK rats, there was no significant difference in the Rev-erbα mRNA levels at the different time points, but the highest values appeared at ZT7. In contrast, the Rev-erbα protein levels showed different patterns between GK and Wistar rats (Fig. 5B). In Wistar rats, the Rev-erbα protein levels significantly increased between ZT2 and ZT10 and started to decrease at ZT12 after gradually reaching the peak at ZT10. These results are consistent with the Rev-erbα protein circadian oscillation in mouse liver, which was barely detectable at ZT2 and peaked around ZT10 (36). In GK rats, the Rev-erbα protein level peaked at ZT2, and after the decrease from ZT2 to ZT10, it slightly increased at ZT12. A comparison between the two groups revealed that the Rev-erbα protein level in Wistar rats at ZT2 was significantly lower than those at any time points in GK rats, while Rev-erbα protein level in GK rats at ZT2 was significantly higher than those at any time points in Wistar rats. In conclusion, the differences observed in the Rev-erbα protein levels reflected a phase-shift effect as well as the absolute expression difference at ZT2.

Distinctive Associations of Rev-erbα in Adipose Tissues with Diabetes Parameters Between GK Rats after Gastric Bypass Surgery and Untreated GK Rats

We further examined the association of Rev-erbα in adipose tissue with diabetes in vivo. We ameliorated diabetes in GK rats by gastric bypass surgery and compared the Rev-erbα mRNA and protein levels between GK rats with or without gastric bypass surgery (Fig. 6).

Fig. 6.

Association of Rev-erbα with diabetes improvement in GK rats after gastric bypass surgery. A: we examined the improvement of diabetes in surgically treated GK rats (GK-S; n = 5), compared with the untreated GK group (GK; n = 6), by measuring the physiological parameters of diabetes (see details in the text). B: representative pictures of adipose tissue morphology from surgically treated and untreated GK rats. Open arrow denotes lymphocyte infiltration, while solid arrows denote nucleus. C: Rev-erbα mRNA and protein levels were measured in adipose tissue from surgically treated and untreated GK rats at ZT2 (n = 5 or 6 per group) and ZT10 (n = 4 per group). The abbreviations are as follows: GK, untreated GK rats; GK-S, surgically treated GK rats; AUC, area under curve. Significance within groups is indicated as follows: **P < 0.01; and *P < 0.05.

First, we examined the improvement of diabetes in GK rats at 3M after gastric bypass surgery, in terms of body weight and several diabetes parameters, including fasting blood glucose, insulin levels, and blood glucose levels after oral glucose and insulin challenges. As seen in Fig. 6A, the body weights were significantly decreased after bypass surgery. Although the fasting blood glucose levels were not significantly different between the GK rats after gastric bypass surgery (GK-S group) and the untreated group (GK group), the fasting insulin levels were significantly lower in the surgical group, indicating improved insulin resistance in the GK-S group. After the oral glucose challenge, the blood glucose levels in the GK-S group at 60 min and 120 min, as well as the AUCs of glucose, were significantly lower than those in the GK group. The decrements of blood glucose levels after insulin injection in the GK-S group were also significantly faster than the GK group. These diabetes parameters indicated that the degree of the diabetes was improved in GK rats at 3M after gastric bypass surgery.

In addition to the diabetes parameters, the adipose tissue was compared between the GK-S and untreated GK groups. As shown in Fig. 6B, the adipocytes in the GK-S group are smaller than those of the GK group. As indicated by the black arrows, lymphocyte infiltration was frequently seen in the GK group, but it was much less apparent in the adipose tissue from the GK-S group. It is well known that inflammation in adipose tissue is related to diabetes (5, 34) and that Rev-erbα is correlated with inflammation (27). These data suggested that the inflammation in adipose tissue was decreased in the GK-S group, and this may also be associated with improved diabetes after gastric bypass surgery.

The Rev-erbα mRNA and protein levels along the circadian time points in adipose tissues from the GK and GK-S groups were further compared at ZT2 and ZT10 (Fig. 6C). Interestingly, both levels indicated that the diabetes improvement in GK-S group. For the mRNA levels, there was no significant difference at ZT2 and ZT10, when comparing GK with GK-S, and this result is similar to that in the comparison between Wistar and GK rats in Fig. 5A. In terms of the Rev-erbα protein level, the Rev-erbα protein levels at ZT2 in the adipose tissue from the GK-S group were significantly lower, compared with those from the untreated GK rats at ZT2. However, the Rev-erbα protein levels at ZT2 from the GK-S group were not different, compared with those at ZT10 from the GK group. Furthermore, the trends of the Rev-erbα protein levels in the GK-S and GK groups between ZT2 and ZT10 were similar to those of the Wistar and GK rats in Fig. 5B. Indeed, the Rev-erbα protein levels significantly increased in the GK-S group and significantly decreased in the GK group, consistent with the finding that the Rev-erbα protein levels increased in Wistar rats and decreased in GK rats. These results indicated that the Rev-erbα protein level in the improved diabetic GK rats was shifted toward that in the Wistar rats. In summary, the parallel relationship between the Rev-erbα protein levels and the diabetes parameters indicates that not only the lower protein levels of Rev-erbα at ZT2, but also the circadian Rev-erbα protein levels are associated with improved diabetes in the gastric bypass-treated GK rats.

Associations of Rev-erbα-Related Molecules with Diabetes Improvements in GK Rats

Heme was recently identified as the physiological ligand for Rev-erbα (37, 53), and weight loss and potentially insulin-sensitizing effects were observed by treating diet-induced obese mice with synthetic Rev-erb agonists derived from heme (43). To investigate the relationship between diabetes and heme in surgically treated GK rats, we measured the mRNA expression of aminolevulinic acid synthase 1 (ALAS1), the rate-limiting enzyme in heme biosynthesis (21, 50). As shown in Fig. 7A, the ALAS1 mRNA expression levels in Wistar rats were lower than those in GK rats at all of the circadian time points, especially at ZT2 and ZT10. This trend was also observed in the comparison between untreated and surgically treated GK rats. Indeed, the expression of ALAS1 mRNA in the surgically treated GK rats was lower than that of the untreated rats. These findings suggested that the diabetes improvement may be coordinated with heme biosynthesis.

Fig. 7.

Relationships between Rev-erbα and its related molecules in diabetes, especially in terms of diabetes improvement of GK rats after gastric bypass surgery. We investigated the Rev-erbα-related molecules, ALAS-1, PPARγ, IL-6, and GSK-3β, and evaluated their relationships to Rev-erbα in terms of diabetes progression. A: circadian ALAS-1 mRNA levels were measured in adipose tissues from Wistar and GK rats (n = 3 per group at each time point), as well as surgically treated and control GK rats (n = 4–6 per group at each time point), normalized to the geometric averaging of EF1A and GAPDH. Significance within groups is indicated as follows: **P < 0.01; *P < 0.05. B: relationship between the mRNA levels of Rev-erbα and PPARγ in adipose from 1M, 3M, and 6M diabetic GK (n = 17) and Wistar control (n = 17) rats. Each line was estimated by linear regression: R = 0.78 for GK rats and R = 0.05 for Wistar rats. C: comparison between IL-6 mRNA levels in adipose tissue at ZT2 from surgically treated (n = 5) and untreated GK rats (n = 6). D and E: effect of Rev-erbα inhibition in adipose tissue on insulin tolerance in the whole body. BIO (100 mg/kg) was injected into 4 wk-old GK rats (n = 3), and the injections were continued for 10 wk, until the animals were 3 mo old.

Recently, the role of Rev-erbα through other molecules, such as PPARγ, has been highlighted in the link between inflammation and insulin resistance (15, 23). Furthermore, Rev-erbα also functions in the regulation of the inflammatory response, involving molecules such as IL-6 and NF-κB, in A7r5 vascular smooth muscle cells (8, 27). Therefore, we analyzed the mRNA levels of these molecules in our animal model. The average mRNA levels of Rev-erbα and PPARγ in adipose tissues from 1M, 3M, and 6M diabetic GK and Wistar rats were not different. However, when we plotted the Rev-erbα and PPARγ mRNA levels of individual rats, we found a positive correlation in diabetic GK rats, as shown in Fig. 7B. However, less correlation was found between the mRNA levels of Rev-erbα and PPARγ in adipose tissues from Wistar rats. The difference between the GK and Wistar rats suggested that Rev-erbα might be linked to primary insulin resistance via PPARγ in adipose tissue. Previously, we found a correlation between the decreased levels of Rev-erbα protein and inflammation improvement in the surgically treated rats, as shown in Fig. 6C. Interestingly, we also found decreased IL-6 mRNA levels in the surgically treated rats, compared with those in the untreated rats (Fig. 7C).

Increased glycogen synthase kinase-3 (GSK-3) in fat tissue is implicated as a potential factor contributing to the development of insulin resistance and T2DM in high-fat diet-fed C57BL/6J mice (13), and GSK-3 phosphorylates and stabilizes Rev-erbα, thus protecting Rev-erbα from degradation (52). We treated 1M GK rats for 2 mo with a specific inhibitor of GSK-3, BIO (26) and found a dramatic decrease in the Rev-erbα protein levels in adipose tissue and improved insulin tolerance in the BIO-treated animals, as shown in Fig. 7, D and E.

DISCUSSION

In this study, we selected one TF and its regulated gene pair in adipose tissue, RORA-Rev-erbα, as a plausible candidate associated with diabetes in GK rats, by applying our computational procedure to the open-source gene expression data of diabetic GK and Wistar rats. After checking the Rev-erbα mRNA levels by our original animal models, the changes in the Rev-erbα protein levels in adipose, liver, and muscle tissues were measured during aging at 1M, 3M, and 6M. By comparing the Rev-erbα protein levels in the three tissues, a distinctive pattern emerged in adipose tissue: an unusual increment in GK rats, and a drastic decrement in Wistar rats with aging. The circadian Rev-erbα protein levels were also compared between GK and Wistar rats at 6M, the time when the Rev-erbα protein levels showed the biggest difference. The circadian Rev-erbα expression was nearly opposite during the day time between the GK rats and Wistar controls, and a significantly low value was detected at ZT2 in Wistar rats, compared with the other measured time points in GK rats. To investigate the specific association of Rev-erbα protein levels in adipose tissue to diabetes, we ameliorated diabetes in 3M GK rats by gastric bypass surgery. After detailed monitoring of the ameliorated condition, we investigated the relationship between the Rev-erbα protein levels and diabetes, by comparing the surgically treated and untreated GK rats at two circadian time points. Interestingly, the trend of the Rev-erbα protein levels in the surgically treated GK rats was similar to that of the Wistar rats, suggesting that the circadian expression of Rev-erbα is related to the degree of diabetes improvement. These analyses revealed that Rev-erbα is associated with diabetes in the adipose tissue of GK rats, as a potential factor in the molecular mechanisms underlying T2DM for diabetes progression through circadian dysregulation.

A relationship between the circadian clock and metabolism, and especially a causal link between the disruption of circadian rhythms and metabolic complications, has been suggested from numerous association studies demonstrating that shiftworkers display increased risks of developing certain metabolic diseases, such as obesity and T2DM (38). Indeed, the molecular basis responsible for the maintenance of circadian rhythms has been identified in adipose tissue (17) and pancreatic islets (24), and furthermore, the circadian regulation of glucose and lipid metabolism has been well established (29–30, 54). As a result of multiple mechanisms, including protein ubiquitination, the relative levels of the BMAL1-CLOCK heterodimer and the CRY-PER heterodimer mRNAs and proteins are expressed out of phase with each other, resulting in a self-maintaining feedback loop that oscillates within a period of ∼24 h. In the negative-feedback loop, Rev-erbα (43), in addition to the nuclear corepressor (NCoR)/histone deacetylase 3 (HDAC3) complexes (1), represses the transcription of BMAL1. The whole body BMAL1 knockout displays altered circadian plasma glucose oscillations and glucose intolerance; however, the specific knockdown in liver shows hypoglycemia and increased insulin sensitivity. Interestingly, NCoR mutant mice with increased BMAL1 activity have a lean phenotype and ameliorated metabolism in vivo (1). With these situations in mind, the following results in our study may provide a clue to the relationship between the circadian clock and metabolism. In Fig. 5B, we observed high circadian expression of the Rev-erbα protein levels in adipose tissue from 6M-old Wistar rats, with a similar pattern to that in mouse liver (36), and in contrast, we found a clear phase-shift of the Rev-erbα protein expression with nearly the opposite pattern in GK rats. The circadian protein expression after bypass surgery is more similar to the normal pattern in Wistar rats, with significantly decreased expression at ZT2 and an increment from ZT2 to ZT10, as shown in Fig. 6C. Interestingly, the Rev-erbα mRNA levels in GK and Wistar rats are not significantly different between serially harvested adipose tissues, as shown in Fig. 5A, consistent with the fact that obesity and T2DM are not related to the mRNA levels of core circadian genes from human white adipose tissue (31). As a result, our observations, including the changes in the Rev-erbα mRNA and protein levels in adipose tissue between GK and Wistar rats in the circadian time points of 6M-old rats (Fig. 5, A and B) and the shifting of the Rev-erbα protein levels in the surgically treated GK rats to those in Wistar rats (Fig. 6C), may reflect the link between the circadian clock and diabetes in the adipose tissue of diabetes-prone animals.

We observed a significant decrease in the ALAS1 mRNA levels when we compared the Wistar and GK-S groups with the GK group, as shown in Fig. 7A. Since ALAS1 displays negative feedback regulation by heme (50), the decrease in ALAS1 in the Wistar and GK-S groups, which is correlated with an increase in the heme concentration, also suggests the decrease of Rev-erbα, with reference to the fact that heme is the physiological ligand for Rev-erbα (37, 53). Thus, these results are quite consistent with the increase of Rev-erbα in GK rats (Fig. 4A) and the diabetes improvement in the GK-S group (Fig. 6A). Thus, the synthetic Rev-erbα agonists are promising antiobesity drugs, with weight loss and potentially insulin-sensitizing effects in diet-induced obese mice (43).

Recent work has shed light on the cellular and molecular mechanisms through which adipose tissue metabolism causes primary insulin resistance in skeletal muscle (15, 23) and liver (8). Indeed, chronic inflammation in adipose tissue triggers insulin resistance in skeletal muscle, by impairing triglyceride deposition. At the molecular level, PPARγ downregulation by TNF-α in adipose cells (15), which decreases triglyceride deposition and increases lipolysis in adipose cells, impairs triglyceride storage. Furthermore, the IL-6 released from adipose tissue is also involved in the development of insulin resistance, in LIKK mouse liver (8). Rev-erbα appears to be a regulator of PPARγ in 3T3-L1 cells, and the constitutive expression of the Rev-erbα protein prevents adipogenesis by inhibiting PPARγ induction (48). Rev-erbα also regulates the inflammatory response, and Rev-erbα reportedly upregulated proinflammatory cytokines, such as IL-6 and NF-κB, in A7r5 vascular smooth muscle cells (27). In our animal model, we found a positive correlation between the Rev-erbα and PPARγ mRNA levels in adipose tissue from diabetic GK rats, but no correlation in Wistar rats in Fig. 7B. Although the regulation of Rev-erbα and PPARγ is very complicated, with the mutual regulation of Bmal1, RORA, and PGC-1, in addition to the protein degradation process (15, 23), Rev-erbα might be linked to primary insulin resistance via PPARγ in adipose tissue. We also found the coordination of the decreased Rev-erbα protein levels with inflammation improvement (Fig. 6C) and decreased IL-6 mRNA levels (Fig. 7C) in surgically treated rats, compared with those levels in untreated rats. These results are consistent with the finding that IL-6 levels were decreased by 72% in mice treated with synthetic Rev-erbα agonists (43).

As shown in Fig. 7, D and E, we observed improved insulin resistance after GSK-3 treatment by using BIO. Although the BIO experiments only provide another association of decreased Rev-erbα protein levels with antedated diabetes, the elevation of Rev-erbα protein levels by GSK-3β may be one of the potential mechanisms contributing to the development of diabetes in diabetic-prone mice, by the increased GSK-3β activity in fat tissue (8).

Our computational approach revealed RORA and ELK3 as key TFs contributing toward T2DM development in GK rats, and we focused on a gene regulated by RORA, Rev-erbα to investigate its association with T2DM progression. Although an investigation of the genes regulated by ELK3 is beyond this work, it may be worthwhile to study the contribution of ELK3 to T2DM development. ELK3 belongs to the ternary complex factor subfamily of the ETS-domain transcription factor family, which includes the cofactors for serum response factor activated by MAPK signal-induced phosphorylation. Indeed, MAPK/ERK inhibition acts via the ETS-1 transcription factor to induce insulin resistance (55).

Perspectives and Significance

In the present work, by a systematic combination of computational analyses and experiments, we revealed the upregulation of Rev-erbα protein expression in adipose tissue with T2DM progression, as well as dysregulation of the circadian rhythm of Rev-erbα in GK rats, especially its association with T2DM amelioration after gastric bypass surgery. In addition, the link between heme and Rev-erbα also revealed by the observation of upregulation of ALAS1, which indicates that the decrease in heme is physiologically significant in diabetic GK rats. Since clock-driven lipoprotein lipase expression was reported to be important in the balance of energy metabolism in Rev-erbα−/− mice (11), which might be directly relevant to the results obtained, our work suggests that the role of Rev-erbα in adipose tissue is one of the key molecules potentially causing metabolic syndromes by circadian disruption. To confirm this hypothesis, adipose tissue-specific Rev-erbα knockout mice in a diabetic background should be established in the future.

Appendix

Detailed Description of our Method to Narrow the Candidates of Key Molecules Responsible for Phenotype Changes

Recently, we developed a computational procedure by combining two methods, to find the key molecules that are activated under particular conditions (33). One is a knowledge-based approach, which estimates the consistency of the network structures with the measured data among the known networks (named “network screening”) (39). Thus, any binary relationships between a transcription factor and its regulated genes, which are detected by network screening, are guaranteed by experiments in previous reports. However, this approach has a pitfall. Unfortunately, the network structures are restricted because of our incomplete knowledge of the networks. To compensate for this restriction, we use another approach to find activated regulatory networks, by inferring the network structures from the measured data based on a mathematical model. We adopted the path consistency algorithm (44) as a mathematical model, in which the conditional independency is considered by calculating the partial correlation coefficients. The original algorithm often stopped in applications to gene expression data that frequently show redundancy, such as similar patterns of expression data among many genes. We partially modified the algorithm to make it applicable to redundant data (40). Indeed, the performance of the combination of the two methods was successfully confirmed to computationally find the master regulators of a brain tumor (41), which were previously identified by combined computational and experimental approaches (9). Furthermore, the analyses of the expression profiles in the livers of diabetic rats revealed several candidates of key molecules along the diabetes progression pathway (33). The following sections provide details about the two computational methods and the enrichment analysis that was used in the computational procedure.

Network Screening

The candidates of the active regulatory networks were detected by network screening (NS) (39). Here, we will briefly summarize the network screening performed in the present study.

First, the regulatory network sets were generated in the same manner as in the previous study (33), as follows. The mouse binary relationships compiled in the TRANSFAC database (49) were used. According to the correspondence between the mouse and rat gene ids, 3,015 binary relationships of 1,507 genes between 503 TFs and 1,123 regulated genes were identified. On the basis of these binary relationships, transcriptional networks were constructed according to the functional gene sets previously defined in the Molecular Signatures Database (MSigDB) (45). In each gene set, the regulated genes in the binary relationships were searched, and if at least one gene was found in the gene set, then the corresponding binary relationships were regarded as a regulatory network characterized by the gene set. In the present study, the reference network comprised 1,760 regulatory networks characterized by biological functions that are composed of 1,195 genes. The numbers of TFs and regulated genes were 335 and 860, respectively.

Then, we calculated the graph consistency probability (GCP) (39), which expressed the consistency of a given network structure with the monitored expression data of the constituent genes in this study. The consistency of a directed acyclic graph, G(Vi, Ej), where Vi is a vertex (I = 1, 2, . . . , nv) and Ej is an edge (j = 1, 2, . . . , ne) in the graph, and the joint density function f (Xi), corresponding to Vi for the graph G with the measured data, is quantitatively expressed by the logarithm of the likelihood, based on the Gaussian graphical model (GN, Gaussian Network); i.e., l(G0)=lni=1nvf(Xi|pa{Xi})=12i=1nvj=1ni{1σi2k=1m(xikj=1niβijxkj)2+ln(2πσi2)},A1 where pa{Xi} is the set of variables corresponding to the parents of Vi in the graph, xik is the measured value of Xi, at the kth point, and ni is the number of variables corresponding to the parents of Vi. Since the likelihood depends on the graph size, we designed a simple procedure to transform the likelihood to the probability for the expression of the graph consistency with the data (32). First, we generated Nr networks under the condition that the networks shared the same numbers of nodes and edges as those of the given networks. Then we defined GCP, as follows, GCP=NsNr,A2 where Ns is the number of networks with larger log-likelihoods than the log-likelihood of the tested network. In the present study, Nr was set to 2,000, and the GCP significance of the given network was set to 0.05.

Modified Path Consistency Algorithm

The modified path consistency (mPCA) algorithm (44) is an algorithm to infer a causal graph composed of two parts: the undirected graph inference by a partial correlation coefficient and the following directed graph construction by the orientation rule. The present method partially exploits the first part of the PC algorithm for the inference of the network structures. A simple example of the PC algorithm is illustrated, as follows:

We assume that five variables, X1, X2, X3, X4, X5, have the following five relationships:

  • 1) X1X2,

  • 2) X2 ∐ (X1, X4),

  • 3) X3X4|(X1, X2),

  • 4) X4∐(X2, X3)|X1, and

  • 5) X5∐(X1, X2)|(X3, X4),

where the symbol, ∐, in the above relationships, means the independence between variables. The PC algorithm reconstructs the above relationships, as follows.

  1. Prepare a complete graph, C, between the five variables.

  2. Test the correlation between two variables by calculating the zeroth-order of the partial correlation coefficient (Pearson's correlation coefficient). From the test, two variable pairs, (X1, X2) and (X2, X4), are excluded, due to the relationships i and ii.

  3. Test the correlation between three variables by calculating the first-order of the partial correlation coefficient of the variable pairs, given one variable. Then, one variable pair, (X3, X4), is further excluded from the updated graph by 2), due to iii and iv.

  4. Test the correlation between four variables by calculating the second-order of the partial correlation coefficient of the variable pairs, given two variables. Then, two variable pairs, (X1, X5) and (X2, X5), are excluded, due to iv.

  5. We could not find any edges adjacent to the three edges in the updated C. Thus, the algorithm naturally stops. Finally, the five relationships emerged completely.

In general, the (m-2)th order of the partial correlation coefficient is calculated between two variables, given (m-2) variables; i.e., rij,rest, between Xi and Xj, given the “rest” of the variables, {Xk} for k=1, 2, . . . , m, and ki, j, and after calculating the (m-2)th order of the partial correlation coefficient, the algorithm naturally stops. However, the algorithm does not usually request the (m-2)th order of the correlation coefficient for the natural stop. This happens because after the variables are excluded, the adjacent variables are often not found, even in the calculation of the lower orders of partial correlation coefficients.

In the actual expression profile data, many genes frequently show profiles with similar patterns. This makes the numerical calculation of correlation coefficients difficult, due to the multicolinearity between the variables. The original PC algorithm accidentally stops, if only one correlation between a pair of variables shows a violation of the numerical calculation. However, in a biological sense, the gene pairs that cause the accidental stop can be interpreted as a case of their high association with each other, in terms of gene expression. Thus, we modified the original PC algorithm to prevent it from accidentally stopping with the highly associated gene pairs, as follows (40). If the calculation of any order of the partial correlation coefficient between the variables is violated, then the corresponding pair of variables is regarded as being dependent. For example, if the first-order correlation coefficient, rij,k, cannot be calculated numerically, due to the multicolinearity between Xi and Xj, then the edge Xi-Xj is kept without the statistical test. The other parts remain unchanged in the modified algorithm. Note that the above modification ensures that the algorithm will naturally stop for the data, including a high correlation.

As seen in the original algorithm, the output is not unique, depending on the calculation order of pairs (44). A permutation test for the calculation order is a convenient way to partly resolve this issue. In this study, the estimation without permutation was empirically adopted as the first approximation, based on the successful estimations of the relationships in our previous studies (40). In addition, one of the most remarkable features of the PC algorithm is that the algorithm removes the pseudocorrelations between the variables (genes) by considering the higher-order partial correlations. If we have the measurement data for a complex network, then we frequently face the more serious issue of the pseudocorrelation, rather than the correlation level. The merit of the PC algorithm may be its ability to identify the real relationships between TFs and their regulated genes.

Definition of Key Molecule Candidates by the Combination of Network Screening and Network Inference

We first referred to the two sets of networks obtained by the network screening (39) and the network inference (40). From each network set, the binary relationships between the TFs and their regulated genes were extracted, only if the regulated genes were included in the expression signature, which is the ensemble of genes with significant differences in gene expression, as statistically estimated by the false discovery rate (FDR) test for multiple comparisons (FDR< 0.05) (6). In the extraction of TFs and their regulated genes, the TF was also cited from the TRANSFAC database (49), but the expression degree of the TF was not considered, due to the small expression changes, even under different conditions. Only the regulated genes that were estimated to directly bind TFs were extracted. The number of genes in the gene expression signatures of the initial stage (period of 4 wk) was 1,582.

Then, we defined the key molecule (KM) candidates from the binary relationships by one criterion. Here, the specificity simply means that the TF emerged only in one tissue of the GK networks, but not in the other tissue networks. Note that in the selection of the TFs, we only selected those that were estimated to regulate the genes, including the expression signature, to consider the enrichment of the regulated genes in the signature. Finally, we selected the KM candidates that were commonly found in the two sets.

Enrichment Analyses

The period-specific signatures are characterized by known gene sets with biological functions. Here, we adopted the gene sets in CP (canonical pathways from KEGG, BICARTA, and REACTOME) of C2 (curated gene sets) in MSigDB (45), which contained 871 gene sets that were manually curated from previous reports. We then estimated the enrichment probability of the genes in the expression signature for each gene set, based on the hypergeometric probability. When the gene set is composed of k genes, and l genes were detected in the period-specific signature, the probability is obtained by P(Xl)=1i=0l(Mi)(NMki)(Nk) where M and N are the total number of genes in the period-specific signature, and the total number of genes in the gene sets, respectively. In this study, we set the significance probability to 0.01.

GRANTS

This work was partly supported by Chinese Academy of Sciences (CAS) Visiting Professorship for Senior International Scientists (Grant No. 2012T1S0012, KH), Major State Basic Research Development Program of China 973 Program (No. 2011CB504003, HZ), National Natural Science Foundation of China (Grant Nos. 81070657 and 61134013, HZ) and NN-CAS Research Foundation (No. NNCAS-2009-1, HZ).

DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

AUTHOR CONTRIBUTIONS

Author contributions: R.Z., C.Y., X.Z., B.Q., and S.S. performed experiments; R.Z., C.S., B.L., and K.H. analyzed data; R.Z. drafted manuscript; C.Y., X.Z., F.L., Y.S., and K.H. interpreted results of experiments; F.L., Y.S., and S.S. prepared figures; K.H. and H.Z. conception and design of research; K.H. and H.Z. edited and revised manuscript; K.H. and H.Z. approved final version of manuscript.

ACKNOWLEDGMENTS

We thank Dr. Go Nagamatsu of Keio University for helpful discussions and comments on the manuscript. The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding agency.

Footnotes

  • * R. Zhang, C. Yan, and X. Zhou are co-first authors.

REFERENCES

  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 26.
  27. 27.
  28. 28.
  29. 29.
  30. 30.
  31. 31.
  32. 32.
  33. 33.
  34. 34.
  35. 35.
  36. 36.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
  51. 51.
  52. 52.
  53. 53.
  54. 54.
  55. 55.
  56. 56.
View Abstract