figure a

Introduction

All mammalian species possess a circadian timing system consisting of a central brain clock in the suprachiasmatic nucleus (SCN) and clocks in peripheral tissues such as muscle, liver and adipose tissue [1, 2]. The fundamental clock mechanism is the autonomous transcriptional–translational feedback loop consisting of core clock genes that oscillate with a period duration of ~24 h [3]. The central clock in the SCN is synchronised to the 24 h rhythm of the environment via retinal input, and synchronises the peripheral clock rhythms via hormones, the autonomic nervous system and the regulation of food intake and body temperature [1]. The circadian timing system and energy metabolism are closely intertwined. At the molecular level, the transcriptional–translational feedback loop has reciprocal connections with oscillations in NAD+, ATP and the cellular redox state [3]. At the cellular level, cultured adipocytes show circadian rhythms in glucose uptake [4], insulin sensitivity [5] and lipolysis [6], and the necessity of the clock gene Arntl for normal adipogenesis [7]. In healthy humans, diurnal rhythms in glucose tolerance with increased glucose tolerance at the onset of the active period were described as long ago as the 1970s [8]. The observation that this diurnal rhythm of glucose tolerance is diminished in humans with impaired glucose tolerance or diabetes mellitus was one of the first clues linking disturbed circadian rhythms to metabolic disease [9]. By definition, circadian rhythms are different from diurnal rhythms in that circadian rhythms are rhythms with a period of approximately 24 h that are regulated by the endogenous circadian timing system and that persist under constant conditions, while diurnal rhythms are rhythms observed under non-constant conditions, for example, during a 24 h light/dark cycle.

According to the circadian desynchrony hypothesis, desynchrony between the various clocks in the circadian timing system and behavioural sleep/wake and fasting/feeding rhythms may be involved in the pathophysiology of obesity and type 2 diabetes [2, 10]. In support of the circadian desynchrony hypothesis, the amplitude of rhythms of adipose clock gene expression is reduced in mouse models of type 2 diabetes [11, 12]. Moreover, specific ablation of the clock gene Arntl in skeletal muscle causes insulin resistance [13], and specific Arntl ablation in adipose tissue causes obesity [14].

In humans, there is a correlation between polymorphisms in the clock genes ARNTL [15] and CRY2 [16] and the risk of type 2 diabetes. It has been observed that the amplitude of clock gene expression rhythms is reduced in peripheral blood leucocytes in humans with type 2 diabetes [17], but these cells are probably not directly involved in the pathophysiology of insulin resistance.

Despite this circumstantial support of the desynchrony hypothesis, the question of whether diurnal gene expression rhythms are altered in metabolic tissues of humans with type 2 diabetes remains to be answered. One previous study detected no differences in diurnal rhythms of gluteal subcutaneous adipose tissue core clock gene expression between individuals with type 2 diabetes and healthy individuals as assessed by quantitative PCR (qPCR) of a limited number of genes [18]. Since another study reported that 25% of the genes in human adipose tissue show a diurnal rhythm [19], we adopted a genome-wide approach to perform a comprehensive rhythmic gene expression analysis. Obesity and type 2 diabetes are strongly correlated, and the majority of patients with type 2 diabetes are obese [20, 21]. In the current study, we therefore compared the in vivo diurnal gene expression rhythms in subcutaneous adipose tissue in obese individuals with type 2 diabetes vs healthy lean control individuals.

Methods

Participants

We included six obese individuals with type 2 diabetes according to the 2010 American Diabetes Association (ADA) criteria [22], BMI 25–40 kg/m2, male sex, aged 30–75 years old, and using no glucose-lowering drugs other than metformin. Inclusion criteria for the six healthy control participants were: BMI ≤25 kg/m2, male sex and age 30–75 years. Exclusion criteria for both groups were: any acute or chronic disease (other than type 2 diabetes) that impairs food digestion, absorption or metabolism, obstructive sleep apnoea, and shift work or crossing more than one time zone in the month prior to study participation. Participants were recruited by announcements at public locations, in local newspapers and patient magazines. We provided study information to 46 individuals with type 2 diabetes, nine individuals were screened, three were excluded (one had obstructive sleep apnoea and two declined to participate). Six individuals with type 2 diabetes provided written informed consent and completed the study. We provided study information to 37 healthy control individuals, six of whom were screened, and all six provided written informed consent and completed the study.

Study protocol

The study was approved by the Institutional Review Board of the Amsterdam UMC (location AMC), and performed according to the Declaration of Helsinki of October 2004. The study was registered at the Netherlands Trial Registry (number NTR3234) and was performed between February 2012 and March 2013 at the Department of Endocrinology and Metabolism of the Amsterdam UMC.

For a detailed description of the study protocol, refer to the electronic supplementary material (ESM). In brief, at baseline, participants recorded food intake and sleep times for 3 days. Subsequently, participants visited the clinical research unit on the morning of day 1. For each participant, an individual zeitgeber time (ZT) 0:00 was determined by the average wake-up time from the baseline sleep log. The individuals with type 2 diabetes were instructed to stop metformin use from 2 days prior to the measurements until study end (a total of 5 days). We provided participants with three identical liquid meals (Ensure Plus [6.3 kJ/ml, 54% energy from carbohydrates, 29% from fat and 17% from protein; Abbott Nutrition, Columbus, OH, USA]) per day at evenly spaced fixed time points (ZT 0:30, ZT 6:00 and ZT 11:30), starting from lunch on day 1. Daily energy intake was set at 105 kJ/kg body weight. Participants returned home and then came back to the clinical research unit on the evening of day 2. Participants slept undisturbed in darkness (0 lux) during their habitual sleep times. On day 3 at ZT 0:00, room lights were turned on at ~150 lux at eye level; participants remained in a semi-recumbent position. Superficial subcutaneous adipose tissue samples were obtained by vacuum suction with a 15 gauge needle, from the four peri-umbilical quadrants, in a random order, on day 2 at ZT 15:30, and on day 3 at ZT 0:15, ZT 5:45 and ZT 11:15. Starting from ZT 0:00 on day 3, blood samples were obtained at 30 min intervals (15 min intervals in the first postprandial hour) until ZT 15:00.

Plasma measurements

Plasma glucose concentrations were assessed with the glucose oxidation method with a Biosen glucose analyser (EKF Diagnostics, Barleben, Germany). Plasma cholesterol and lipids were measured with a Cobas 8000 modular analyser (Roche Diagnostics, Rotkreuz, Switzerland), plasma insulin and cortisol using a chemiluminescent immunometric assay on an Immulite 2000 system (Siemens, Breda, the Netherlands), and plasma NEFA using an enzymatic calorimetric method [NEFA-HR(2) Assay, Wako Chemicals, Neuss, Germany].

Adipose tissue RNA sequencing

For a detailed description of the adipose tissue RNA sequencing and qPCR techniques, see the ESM. In summary, RNA isolation was performed according to the RNeasy-Mini protocol for animal tissues and cells on the Qiacube (Qiagen, Hilden, Germany). cDNA was constructed using the Ovation RNA-Seq System V2 (NuGen, Manchester, UK) and the product was purified with the MinElute PCR Purification kit (Qiagen). cDNA was sheared by sonication using a Covaris sonicator (Woburn, MA, USA) and subsequently size selected by double AMPure XP bead purification (Agencourt, Beckman Coulter, Brea, CA, USA). Barcoded adaptor ligated library construction was performed using the 5500 SOLiD Fragment Library Core Kit (Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA). Sequencing was performed on the SOLiD 5500 wildfire sequencer (Life Technologies). The relative expression of PER1, PER2, PER3 and ARNTL (BMAL1) was validated by qPCR (for primer sequences, see ESM Table 1).

Statistical analysis

Baseline measurements, ambulant measurements and plasma values

Normally distributed data are expressed as means ± SD (for baseline variables) or standard errors of the mean (SEM) (for outcome variables); non-normally distributed data are expressed as medians (interquartile range).

Incremental AUCs (iAUCs) for postprandial plasma glucose and insulin were calculated using the trapezoid rule with GraphPad Prism for Windows (version 5.01, GraphPad Software, San Diego, CA, USA), where the postprandial period was defined as the time up to 210 min after meal onset. Actigraphy data were analysed according to the method described by te Lindert and Van Someren [23].

All further statistical analyses were performed with IBM SPSS Statistics (version 21, IBM, Armonk, NY, USA). Baseline between-group differences were assessed with the Mann–Whitney U test for non-normally distributed variables, and the independent samples two-sided t test for normally distributed variables.

AUC and fasting values were log10-transformed or loge-transformed to achieve a normal distribution if necessary and subsequently analysed with a linear mixed-effects model for repeated measurements, with time, group, and time×group interaction as fixed effects using heterogeneous compound symmetry (CSH) as covariance structure. Subsequent differences between time points were assessed with the Wilcoxon signed rank test for non-normally distributed variables and the paired samples two-sided t test for normally distributed values. Graphs were made with GraphPad Prism.

RNA sequencing data analysis

Alignment of sequence reads was performed in colorspace for the individual samples against the GRCh37/hg19 reference genome using the Lifescope aligner (v2.5.1) from Applied Biosystems (Thermo Fisher Scientific) applying the ‘whole transcriptome mapping’ workflow with a mapping quality cut-off (bamgen.mqv.threshold) set to 20. Counts were obtained using the corresponding genome annotation (as provided by Applied Biosystems) and applying the default values of the Lifescope workflow. The obtained counts were used as input for the consecutive analysis using the statistical computing environment R (version 3.1.0, the R foundation, Vienna, Austria). Only features with counts per million mapped reads (CPM) of more than one in two or more samples were retained in the analysis (16818 genes). Counts were normalised by applying the trimmed mean of M-values method [24]. Dispersion estimation and differential expression analyses were performed using edgeR (v3.6.7 [25]).

Differences in expression over time in either healthy controls or patients were determined using a generalised linear model likelihood ratio test with a nested factorial design with a fixed main effect for group (control/type 2 diabetes), and nested interactions of group with participant and time point. Rhythmic expression was analysed by determining differential expression at either ZT 15:30, ZT 5:45 or ZT 11:15 versus the ZT 0:15 baseline using a likelihood ratio test for all three relevant coefficients simultaneously (analogous to an ANOVA F-test for a normal linear model). For rhythmic genes, amplitude was defined as the difference between highest and lowest log2 CPM value. Between-group differences in amplitude were assessed for two gene groups of interest (genes that are rhythmic in both groups, and the core clock genes) with a related samples Wilcoxon signed rank test. A four-time-point analysis is not suitable for sinusoidal curve fitting [26, 27]. Therefore, acrophases were determined manually by categorising genes as ‘increase over the day’, ‘decrease over the day’, ‘peak at mid-day’ or ‘trough at mid-day’.

Differences in expression between healthy controls and patients at individual time points were determined using a generalised linear model likelihood ratio test on the four corresponding contrasts for a design matrix with a fixed main effect for each combination of group and time point. Resulting p values were corrected for multiple testing using the Benjamini–Hochberg FDR adjustment. From the genes with an FDR<0.1, we identified genes with a different mesor by selecting those genes that were consistently up or downregulated at all time points. Original gene symbols were reannotated using biomaRt (v2.20.0, Ensembl v77 [28]). Heatmaps of the mean expression values per time point of the selected genes were created using the heatmap.2 function of gplots (v2.17.0, https://cran.r-project.org/web/packages/gplots/index.html) with row-wise scaling and gene-wise clustering with complete linkage and Pearson correlation distance.

The raw sequencing data and supplementary files for rhythmic expression analysis have been deposited in NCBI Gene Expression Omnibus (GEO series accession number GSE104674); a description of the GEO datasets can be found in the ESM Methods.

Enrichment and upstream regulator analysis

Enrichment and upstream regulator analysis was performed with Ingenuity Pathway Analysis (Qiagen, June 2016 release). For this analysis, gene names were changed into their alternative gene names when necessary in order to be recognised by Ingenuity, using the HGNCHelper R package, followed by manual correction using the NCBI Gene database (see tab 1 of GEO dataset 2, accession number GSE104674). Genes were included if they reached the significance level FDR<0.1. Genes detected in our RNA-Seq experiment were used as the background set. For rhythmic genes, enriched canonical pathways were identified with Fisher’s exact test with the significance level –log10[P]>3. For up- or downregulated genes, average fold change values were used as differential expression input values, and enriched canonical pathways were identified with Fisher’s exact test with the significance level –log10[P]>3 and │z│≥2. For upstream analysis, we identified regulators with the significance level –log10[P]>3 and │z│≥2. Supplementary files with Ingenuity Pathway Analysis results have been deposited under GEO accession number GSE104674; a description of the GEO datasets can be found in the ESM Methods.

Results

The characteristics of the individuals with type 2 diabetes and the healthy lean control individuals are shown in Table 1. Sleep/wake logs, food diaries and actigraphy measurements showed no significant between-group differences in compliance or sleep quality, except for a 10 min difference between prescribed breakfast time and actual breakfast time on day 2 in the individuals with type 2 diabetes (ESM Table 2).

Table 1 Participant characteristics
Full size table

RNA sequencing of subcutaneous adipose tissue samples detected a total number of 16,818 expressed genes. Diurnal rhythms of biological processes are characterised by their amplitude, acrophase, period and mesor (Fig. 1a). Our four-time-point approach enabled us to assess amplitude, acrophase and mesor of all genes in the transcriptome. We first analysed the core clock genes. The clock genes ARNTL and the CLOCK orthologue NPAS2, which are part of the positive limb of the transcriptional–translational feedback loop, showed the expected increase in expression over the day in both groups. CLOCK showed a significant rhythm in healthy control individuals, but not in individuals with type 2 diabetes (Fig. 1b–d). The clock genes PER1, PER2, PER3 and CRY2, which are part of the negative limb of the transcriptional–translational feedback loop, showed the expected decrease over the day in both groups. CRY1 showed no significant rhythm in healthy control individuals, but did show a significant rhythm in individuals with type 2 diabetes (Fig. 1e–i). The core clock genes overall showed a decreased amplitude in individuals with type 2 diabetes compared with healthy control individuals (Wilcoxon signed rank test, p = 0.025, Fig. 1j). Data for four clock genes (PER1, PER2, PER3 and ARNTL [BMAL1]) were validated by qPCR (ESM Fig. 1a–h).

Fig. 1
figure 1

 Core clock genes showed reduced diurnal amplitudes in individuals with type 2 diabetes. (a) Circadian rhythms and diurnal rhythms can be described in terms of amplitude, acrophase, period and mesor. (b–d) Positive clock genes. (b) CLOCK showed a significant rhythm in healthy control individuals (FDR 0.004), but not in individuals with type 2 diabetes (FDR 0.83). (c) NPAS2 and (d) ARNTL (BMAL1) showed a diurnal rhythm in both groups (FDR <0.001). (e–i) Negative clock genes. (e) PER1, (f) PER2, (g) PER3 and (i) CRY2 showed a diurnal rhythm in both groups (FDR<0.1). (h) CRY1 showed no significant rhythm in healthy control individuals (FDR 0.120), but did show a significant rhythm in individuals with type 2 diabetes (FDR 0.004). (j) Overall, the amplitude of expression was reduced in individuals with type 2 diabetes vs healthy individuals (Wilcoxon signed rank test, p=0.025). The units for RNA sequencing data are log2 CPM normalised expression values. Note that y-axes of RNA sequencing data have different scales. Data are means ± SEM. n=6 per time point

Full size image

Next, we analysed the number of genes with a significant diurnal rhythm with a generalised linear model likelihood ratio test, using a false discovery rate (FDR) <0.1 as a cut-off. Out of the 16,818 genes, we identified three gene groups of interest: 1237 genes with a significant rhythm in healthy control individuals only, 119 genes with a significant rhythm in individuals with type 2 diabetes only, and 184 genes with a significant rhythm in both the control individuals and individuals with type 2 diabetes (Figs 2a and 3; see tab 1 and 3 of GEO dataset 1, accession no. GSE104674). Rhythmic genes included genes encoding enzymes, transcription regulators, kinases, transporters and membrane receptors (ESM Fig. 1i, j).

Fig. 2
figure 2

 Diurnal rhythms of adipose tissue gene expression were reduced in individuals with type 2 diabetes compared with healthy control individuals. (a) Venn diagram showing the number of rhythmic genes in each group (significance level FDR<0.1). (b) Scatterplot of differences in average amplitude (control – T2D) of the 1237 genes with a significant rhythm in healthy control individuals, the 184 genes with a significant rhythm in both control individuals and individuals with type 2 diabetes and the 119 genes with a significant rhythm only in the individuals with type 2 diabetes. (c, d) Enrichment analysis of (c) genes that were rhythmic in healthy control individuals only and (d) genes that were rhythmic in both groups, performed with IPA. The ratios indicate the number of genes from a canonical pathway with a significant diurnal rhythm, as a proportion of the total number of genes in the canonical pathway. (e–h) Examples of genes with a significant rhythm in healthy control individuals but not in individuals with type 2 diabetes: the genes included (e) LPL (encoding lipoprotein lipase), (f) PPARGC1A (encoding peroxisome proliferator-activated receptor γ, coactivator 1 α), (g) ACACA (encoding acetyl-CoA carboxylase α) and (h) ACSL1 (encoding acyl-CoA synthetase long-chain family member 1). (i–k, f) Examples of genes with an overall up- or downregulation (i.e. a different mesor) in individuals with type 2 diabetes compared with healthy control individuals. Upregulated genes included (i) FGF1 (encoding fibroblast growth factor 1), (j) LEP (encoding leptin) and (k) CYP19A1 (encoding aromatase). Downregulated genes included (f) PPARGC1A, (l) ARNTL2 and (m) SIRT1 (encoding sirtuin 1). (n–p) Among the genes with a phase difference are (n) PDK4 (encoding pyruvate dehydrogenase kinase 4), (o) PFKFB3 (encoding 6-phosphofructo-2-kinase) and (p) IRS2 (encoding insulin receptor substrate 2). The units for RNA sequencing data are log2 CPM normalised expression values. Data are means ± SEM. n=6 per time point. T2D, type 2 diabetes

Full size image
Fig. 3
figure 3

 Heatmaps showing reduced diurnal rhythms in adipose tissue gene expression in individuals with type 2 diabetes compared with those in healthy control individuals. Heatmaps of the standardised average log2 CPM normalised expression values per time point and group (columns) for (a) 1237 genes that were rhythmic in healthy control individuals only, (b) 184 genes that were rhythmic in both groups, and (c) 119 genes that were rhythmic in type 2 diabetic individuals only. n=6 per time point. T2D, type 2 diabetes

Full size image

We performed an enrichment analysis on these gene groups of interest. The 184 genes with significant rhythmicity in both healthy controls and patients were enriched in canonical pathways, including ‘circadian rhythm signaling’, ‘TR/RXR activation’ and ‘adipogenesis pathway’ (Fig. 2d; see tab 2 in GEO dataset 2 accession no. GSE104674). These 184 genes that showed rhythmicity in both groups showed a reduced diurnal amplitude in individuals with type 2 diabetes compared with healthy controls (Wilcoxon signed rank test, p = 0.002, Fig. 2b). For the 184 genes with a rhythm in both groups, diurnal acrophase was similar between groups for 171 genes (93%) with 50 genes showing an increase over the day, 68 genes showing a decrease over the day, 35 genes showing a peak at mid-day and 18 genes showing a trough at mid-day. Thirteen genes showed a phase difference, including the insulin-sensitive genes PDK4, PFKFB3 and IRS2 (Fig. 2n–p).

The 1237 genes with a significant diurnal rhythm in healthy control individuals only, were enriched in canonical pathways, including ‘PPARα/RXRα activation’ (including the gene LPL), ‘AMPK signaling’ (including the genes PPARGC1A and ACACA) and ‘cAMP-mediated signaling’ (Fig. 2c, e–g; see tab 2 in GEO dataset 2, accession no. GSE104674). The 119 genes with a significant rhythm in individuals with type 2 diabetes only were not enriched for any canonical pathways (ESM Fig. 1k, l; see tab 2 in GEO dataset 2, accession no. GSE104674).

To identify additional pathways potentially involved in the pathophysiology of type 2 diabetes, we analysed genes that were consistently either up- or downregulated at all time points (i.e. with a different mesor, Fig. 1a). We identified 1392 genes (out of 16,818 genes) that were either upregulated (708 genes) or downregulated (684 genes) at all time points in individuals with type 2 diabetes compared with healthy controls (FDR<0.1) (see tab 2 in GEO dataset 1, accession no. GSE104674). Enrichment analysis with ingenuity pathway analysis (IPA) revealed the upregulated canonical pathway ‘LPS/IL-1 mediated inhibition of RXR function’ (–log10[p value] = 3.6, z score 3.2, ratio 30/172) and the downregulated canonical pathway ‘EIF2 signaling’ (–log10[p value] = 9.1, z score –3.4, ratio 41/165). Specifically, consistently upregulated genes included FGF1, the leptin gene LEP and the aromatase gene CYP19A1. Downregulated genes included ARNTL, ARNTL2 and SIRT1 (Figs 1d and 2i–m; see tab 3 in GEO dataset 2, accession no. GSE104674).

We performed IPA upstream regulator analysis to identify potential regulators involved in the pathophysiology of type 2 diabetes. The upstream analysis predicted increased activity of the transcriptional regulators ETS homologous factor (EHF) and mouse double-minute 2 homologue (MDM2), the cytokine TNF, the transporter synoviolin 1 (SYVN1) and the growth factor TGFβ1 in individuals with type 2 diabetes compared with healthy control individuals. Regulators with predicted decreased activity were the transcriptional regulator forkhead box A1 (FOXA1) and the protein complex 5′ AMP-activated protein kinase (AMPK) (see tab 4 in GEO dataset 2, accession no. GSE104674).

On day 2 of the fixed meal schedule, the healthy control individuals showed a diurnal rhythm in postprandial glucose excursions, with an increase in the glucose iAUC from breakfast to lunch, and no change in the iAUC from lunch to dinner. In contrast, the individuals with type 2 diabetes showed a decrease in the glucose iAUC from breakfast to lunch, and no change in the iAUC from lunch to dinner (linear mixed-effects model, time×group interaction p = 0.004) (Table 2).

Table 2 Glucose and insulin iAUC values
Full size table

On the third day, the diurnal rhythm of plasma glucose iAUC at the clinical research unit showed a difference in rhythm between healthy control individuals and individuals with type 2 diabetes (linear mixed-effects model, time×group interaction p = 0.006) that was in line with the results in the ambulant setting (Fig. 4a and Table 2). Plasma insulin iAUC values showed no significant difference in rhythm (linear mixed-effects model, time×group interaction p = 0.194), nor an overall effect of time (linear mixed-effects model, time p = 0.522) (Fig. 4b and Table 2).

Fig. 4
figure 4

 Plasma values at the clinical research unit after ingestion of three identical meals. (a) The diurnal rhythm of postprandial glucose excursions was reduced in individuals with type 2 diabetes compared with that in healthy control individuals, despite higher plasma glucose values in individuals with type 2 diabetes. n=5 for 6 time points, n=6 for 70 time points. (b) Insulin levels were increased in individuals with type 2 diabetes compared with healthy control individuals, but neither group showed a diurnal rhythm in postprandial insulin excursions. n=5 for 7 time points, n=6 for 69 time points. (c) Plasma NEFA levels were higher in individuals with type 2 diabetes with a slower recovery compared with healthy control individuals, resulting in a decrease of preprandial NEFA levels over the day. n=4 for 1 time point (T2D ZT 3:30), n=5 for 5 time points, n=6 for 70 time points. (d) Plasma cortisol rhythm was not different between groups. n=6 per time point. Data are means ± SEM

Full size image

Preprandial plasma NEFA showed a different pattern over the day between the groups. Whereas in the healthy control individuals there was a slight increase over the day, in the individuals with type 2 diabetes there was a decrease in preprandial NEFA levels over the day (linear mixed-effects model, time×group interaction p = 0.033). Trough postprandial NEFA levels were higher in the individuals with type 2 diabetes vs healthy control individuals (linear mixed-effects model, group p = 0.015) (Fig. 4c and ESM Table 3).

Plasma cortisol concentrations showed the expected diurnal pattern, with high morning levels and a subsequent decrease over the day (linear mixed-effects model, time p < 0.001). There was no difference between groups in diurnal rhythm (linear mixed-effects model, time×group interaction p = 0.185) or baseline levels (independent samples two-sided t test, p = 0.112) (Fig. 4d).

Discussion

In the present paper, we provide the first transcriptome atlas of diurnal gene expression rhythms in white adipose tissue of obese individuals with type 2 diabetes and healthy lean control individuals. In individuals with type 2 diabetes, the amplitude of core clock gene oscillations was reduced, and the total number of genes with a significant diurnal oscillation was more than fourfold lower compared with that in healthy control individuals. Canonical pathways with reduced rhythmicity in individuals with type 2 diabetes included ‘AMPK signaling’ and ‘cAMP-mediated signaling’. Canonical pathways with an elevated mesor in individuals with type 2 diabetes compared with control individuals were enriched in ‘LPS/IL-1 mediated inhibition of RXR function’, whereas genes with a decreased mesor in individuals with type 2 diabetes were enriched in ‘EIF2 signaling’. Upstream analysis of genes with a different mesor identified potential pathophysiological transcriptional regulators, such as SYVN1 and FOXA1. The reduction in the rhythmicity of adipose tissue gene expression was paralleled by a loss of the diurnal rhythm in postprandial glucose excursions.

Previous rodent studies [11,12,13,14, 29] and human polymorphism studies [15, 16] suggested that reduced diurnal rhythms of gene expression may be a predisposing factor for development of insulin resistance. Here we show for the first time that, in line with this hypothesis, the diurnal rhythm of subcutaneous adipose tissue gene expression is reduced in obese individuals with type 2 diabetes. It is known that increased adipose tissue lipolysis causes elevated plasma NEFA levels, which induces liver and muscle insulin resistance [30]. Patients showed a loss of rhythm for many genes in the canonical pathways ‘AMPK signaling’ and ‘cAMP-mediated signaling’ which both regulate lipolysis via protein phosphorylation [31, 32]. The reduced diurnal rhythms in these genes may reflect reduced synchrony between metabolic gene expression and the diurnal rhythm of food intake, which may contribute to increased daytime lipolysis in individuals with type 2 diabetes. However, since we did not assess protein phosphorylation status, this hypothesis remains to be proven in follow-up experiments.

The lower amplitudes of clock gene oscillations are in line with rodent models of type 2 diabetes showing reduced diurnal rhythms in adipose tissue clock gene expression [11, 12]. However, our finding of reduced clock gene amplitude contrasts with a previous study that found no difference in clock gene oscillations between obese individuals with type 2 diabetes, healthy obese individuals and healthy lean individuals [18]. This difference may be related to the chosen subcutaneous adipose tissue depot (abdominal in our study vs gluteal [18]), since we know from rodent studies that clock gene rhythms may differ between white adipose tissue depots [33]. In addition, the accumulation of subcutaneous adipose tissue in the abdominal subcutaneous compartment probably has a greater contribution to the development of insulin resistance than the gluteal subcutaneous compartment [34,35,36].

Since we used an observational case–control design, a limitation of our study is that we cannot draw conclusions on a causal relationship between the observed alterations in diurnal gene expression rhythms and metabolic disease. The observed loss of rhythmic gene expression in patients with type 2 diabetes may be due to reduced intracellular clock functioning, metabolic disturbances, altered rhythms in other neuro-hormonal signals, or a combination of these factors. For example, we considered the hormone cortisol as a potential candidate, because cortisol affects circadian clock gene expression [37], but we found no difference in the cortisol rhythm between healthy control individuals and individuals with type 2 diabetes.

We compared obese individuals with type 2 diabetes with healthy lean control individuals, since obesity and type 2 diabetes show a strong correlation in the insulin-resistant subtype of the metabolic syndrome [21]. Consequently, a second limitation of our study design is that we cannot deduct whether the differences in gene expression are mainly due to obesity or due to type 2 diabetes, and additional studies are needed to assess which factor is most relevant.

We standardised the food intake schedule for 3 days. Rodent studies have shown that the timing of food intake is an important zeitgeber for adipose tissue clock gene and metabolic gene expression [38,39,40], and in humans clock and metabolic gene expression in adipose tissue responds to feeding status [19]. Therefore, it is possible that the restriction of the food intake period and the reduction of the number of eating episodes in our study design, strengthened clock and metabolic gene expression rhythms in both individuals with type 2 diabetes and healthy control individuals compared with free-living conditions. On the other hand, with the current experimental set-up we can exclude the possibility that the differences observed are due to (acute) effects of differences in sleep/wake or feeding/fasting cycles.

Several physiological features of diurnal gene expression were conserved in the individuals with type 2 diabetes. Genes with preserved rhythmic expression showed virtually no phase differences between individuals with type 2 diabetes and healthy control individuals. Exceptions included the insulin-dependent genes (PDK4, IRS2 [41], and PFKFB3). The phase changes in these genes are probably due to prolonged postprandial insulin signalling, since we observed that insulin levels of individuals with type 2 diabetes do not return to baseline, in contrast to healthy control individuals.

There were 119 genes that showed a rhythm in the individuals with type 2 diabetes but not in healthy control individuals. A major proportion of these 119 genes encode for proteins with cellular functions, including enzymes, transcriptional regulators and kinases. However, as a group, these 119 genes were not enriched in any canonical pathways, which prohibits general statements about the function of this gene group. In individuals with type 2 diabetes, the increased amplitude of the expression of genes responding to lipid metabolism (such as the sterol transporter ABCG1) or insulin (such as the cyclic nucleotide transporter ABCC5) with a diurnal peak or trough expression at ZT 0:15 may be explained by the larger difference between ZT 0:15 and the other time points in terms of plasma insulin, preprandial NEFA levels and, probably, intracellular metabolite levels in individuals with type 2 diabetes vs healthy control individuals.

Parallel to the reduced amplitude in gene expression rhythms, we found a reduced diurnal rhythm in postprandial glucose excursions in individuals with type 2 diabetes. The reduced rhythm in postprandial glucose excursions was expected based on earlier observations [9] and may be related to alterations in the diurnal rhythm of insulin sensitivity [42], although the underlying mechanism remains to be elucidated. If the reduction in diurnal adipose tissue gene expression rhythms can be extrapolated to muscle and liver tissue, this may offer a potential explanation, given that the molecular clock affects muscle [13] and liver [43] insulin sensitivity [2].

In rodent models of obesity, the natural compound nobiletin exerted a beneficial metabolic effect through RAR related orphan receptor (ROR)α- and RORγ-mediated enhancement of the Clock gene rhythm [44]. Since we observed a reduction of the CLOCK rhythm in individuals with type 2 diabetes, future studies should investigate whether nobiletin has beneficial effects in humans with obesity and type 2 diabetes.

To identify additional potential targets for the treatment of type 2 diabetes, we analysed the 1392 genes with a different mesor, i.e. the genes that show either upregulation or downregulation at all time points in individuals with type 2 diabetes compared with healthy control individuals. For example, the downregulated canonical pathway ‘EIF2 signaling’ is involved in translation initiation, and a previous rodent study reported that induction of EIF2 signaling protects against endoplasmic reticulum stress in type 2 diabetes [45]. Potential upstream regulators with increased activity included TNF and TGFβ1, consistent with literature on the role of these inflammatory mediators in the pathophysiology of insulin resistance [46, 47]. Among the potential upstream regulators with increased activity that may represent adipocyte treatment targets was SYVN1, which has a role in the regulation of endoplasmic reticulum protein degradation, and was shown to increase body weight by interacting with peroxisome proliferator-activated receptor γ, coactivator 1β (PGC1β) in mice [48]. Potential upstream regulators with decreased activity were FOXA1, a transcription factor that has been implicated in the development of obesity [49], and the energy sensor AMPK. These potential treatment targets warrant further investigation, since our study design cannot distinguish whether altered expression levels are due to obesity or type 2 diabetes. If future studies show that altered expression is specific for type 2 diabetes, preclinical therapeutic studies are justified.

In conclusion, we provide evidence of reduced diurnal gene expression in adipose tissue of obese individuals with type 2 diabetes, compared with healthy lean control individuals. Core clock genes showed reduced expression amplitudes, the total number of rhythmic genes was strongly reduced and the rhythm in metabolic pathways was lost in individuals with type 2 diabetes. We hypothesise that the reduced diurnal expression amplitudes may contribute to the development of insulin resistance via increased adipose tissue lipolysis. Further investigation is needed to explore potential targets for the treatment of obesity and/or type 2 diabetes, including clock enhancement, and induction of EIF2 and FOXA1.