Introduction

Therapeutic sleep deprivation (SD) rapidly induces robust antidepressant effects in a large proportion of major mood disorder patients suffering from a depressive episode1,2,3,4,5. The effects of the treatment are transient as relapse is usually observed after recovery sleep. Although SD is considered useful and recommended, it is rarely applied in the clinical routine6. The mechanisms through which SD exerts its antidepressant effects nevertheless offer important insights into the biological factors involved in depression and antidepressant response, and have been the focus of recent research7,8,9,10,11. Work in humans12,13,14,15 as well as animals16,17,18 consistently documents the effects of mistimed or insufficient sleep and sleep deprivation on circadian gene expression (such as CLOCK, ARNTL [BMAL1], PER1, PER2, PER3, etc.) as well as on genes involved in related biological processes such as inflammatory, immune and stress response19,20. A prominent hypothesis about the antidepressant mechanism underlying SD is that it restores circadian rhythmicity, which is often dysregulated in depression, via resetting clock gene transcription21,22. The well-controlled nature and rapidity of response to SD treatment23 renders it a promising context to investigate associated biological measures such as gene expression.

While no systematic investigation of gene expression changes in depressed patients undergoing SD treatment has been conducted to date, the study of gene expression in major depressive disorder (MDD) has raised the idea that genes associated with MDD are enriched for inflammation and immune response pathways, which may be linked to the sleep disturbances observed in depression24,25,26. Circadian rhythms are found in the majority of physiological processes and the immune system is no different, with alterations of these rhythms leading to disturbed immune responses27. The immune system and circadian clock circuitry crosstalk, with immune challenges and mediators, such as cytokines, also feeding back to affect circadian rhythms28. Cytokines, including chemokines, interferons and interleukins, are integral to sleep homeostat regulation and can modulate behavioural and physiological functions29.

We recently conducted a naturalistic study, which aimed to examine clinical and genetic factors predicting response to SD30. This was conducted in a sample of major mood disorder inpatients experiencing a depressive episode (n = 78) and healthy controls (n = 15). Briefly, 72% of patients responded to SD. Responders and non-responders did not differ in self/expert assessed symptom ratings or chronotype, but mood differed. Response was associated with lower age and later age at lifetime disease onset. Higher genetic burden of depression was observed in non-responders than healthy controls, with responders having intermediate risk scores.

The present study now aimed to examine gene expression changes accompanying SD and recovery sleep during a depressive episode using a longitudinal design, looking at changes in peripheral blood gene expression in the same sample. Gene expression changes occurring after SD and recovery sleep, as well as associated with response and case-control status, were explored. In a hypothesis-driven approach, expression patterns of circadian genes were investigated. For a systematic search, transcriptome-wide gene-by-gene and gene set enrichment analyses were performed, while longitudinal gene set analysis31 explored dynamics in expression trajectories at the functional gene set level.

Methods and materials

Participants

The present sample has been described elsewhere30. Seventy-eight inpatients (34 females; age mean ± standard deviation = 43.54 ± 14.80 years) presenting with a depressive episode (unipolar, n = 71; bipolar I, n = 6; and bipolar II, n = 1) and on stable medication for ≥5 days participated in this study. Depression was diagnosed according to ICD-10 criteria. Patients were recruited from consecutive admissions to the depression unit of the Department of Psychiatry and Psychotherapy of the Central Institute of Mental Health (CIMH), Mannheim, Germany. Prescribed medication included typical and atypical antidepressants, lithium, and adjunct therapies (anticonvulsants, antipsychotics and sleeping agents). Fifteen healthy controls (eight females; 40.53 ± 15.90 years) with no history of psychiatric/somatic disorders were recruited through an online advertisement on the CIMH website. The criteria for inclusion and exclusion were the same as the criteria for patients, except controls needed to lack psychiatric diseases, which was evaluated via Structured Clinical Interview for DSM-IV Axis II disorders (SCID-II) prior to SD. The investigation was carried out in accordance with the Declaration of Helsinki and approved by the local ethics committee. All participants provided written informed consent following a detailed explanation of the study.

Sleep deprivation

On Day 1, participants gave informed consent and entered the study (see Fig. 1). SD was conducted from Day 2–Day 3, in small groups of 1–5 participants under staff supervision. Participants were free to move around during the night of the SD protocol and were supervised by staff and occupied with activities such as games and walks to ensure wakefulness. Resting in bed was not allowed and SD was carried out on the ward under regular ward lighting conditions (no bright or dim lighting). Food intake was not restricted and some patients consumed snacks (i.e. bread) during SD. Participants underwent the same protocol irrespective of habitual sleep/wake timing; patients followed ward routines (i.e. lights on, lights out times) prior to inclusion to the study. Participants remained awake from ~0600 h on Day 2 to 1800 h on Day 3 (36 h). On Day 3, participants underwent recovery sleep from 1800–0100 h. Sleep phase advance was then carried out, shifting sleep 1 h forward each day until the patient’s regular sleep pattern was reached. Controls participated alongside patients. Response was assessed by the senior clinical researcher using the Clinical Global Impression Scale for Global Improvement/Change in the afternoon on Day 3.

Fig. 1: Experimental procedure.
figure 1

Patients entered the study on Day 1 and underwent sleep deprivation for ~36 h from Day 2 to Day 3 before undergoing recovery sleep. Response was assessed with the CGI-C in the afternoon of Day 3 before recovery sleep. Blood for gene expression was taken at the same time (0600–0730 h) on Days 2, 3 and 4 (T1, T2 and T3, respectively).

Full size image

Data collection

Blood samples were collected in RNA-stabilizing PAXgene tubes (Qiagen, Hilden, Germany), processed according to standard procedures, and stored at −80 °C until analysis. Blood was collected at the same time (between 0600–0730 h) on Day 2 before SD (T1), on Day 3 after SD (T2) and on Day 4 (T3) (see Fig. 1). The number of samples decreased over time points (Table 1) due to non-participation.

Table 1 Samples included for analysis at different time points.
Full size table

Sample preparation and analysis

RNA and microarray analyses

Laboratory analyses were performed using standard methods (see Supplementary Information, SI).

Gene expression data pre-processing

A Custom CDF Version 20 with ENTREZ-based gene definitions was used to annotate the Affymetrix GeneChip™ Human Gene 2.0 ST arrays used for gene expression profiling32. The Raw fluorescence intensity values were normalized applying quantile normalization and RMA background correction using SAS JMP11 Genomics, version 7 (SAS Institute, Cary, NC, USA). The final dataset comprised mRNA expression targeting 24,733 unique genes for each time point per participant.

Data analysis

Analyses were conducted in R (Microsoft R Open 3.4.2). Significance was set at FDR q < 0.05.

Gene-based analysis

Analyses were conducted using lme4 (Version 1.1–17). Linear mixed effects models with random intercepts were fitted to examine gene expression differences between T1 and T2 (‘effect of time point’). Three main models were fitted: effect of time point in ‘all patients’ (M1), in ‘responders vs. non-responders’ (M2) and in ‘patients vs. controls’ (M3). In all analyses, age and sex were included as covariates.

Models were specified as follows: for each transcript, likelihood ratio tests were calculated between two models (h1 vs h0). In both models, gene expression was specified as the dependent variable, covariates were specified as fixed effects and the individual was specified as a random effect. h1 contained the comparison of interest specified as a fixed effect while h0 was a reduced model without it. That is, in M1, the ‘effect of time point’ was the only difference between h1 and h0, while in M2 and M3, the effect of interest was the interaction of the comparison group status (i.e. responder/non-responder, and case/control) and ‘effect of time point’.

Additionally, whether differences in expression levels at T1 ‘baseline’ were informative about response and disease status was examined using linear models, controlling for sex and age.

Targeted examination of differential expression of circadian genes

To determine differential expression of circadian genes, a list of genes comprising both traditional clock genes (in order of highest signal-to-noise ratio predicting circadian rhythmicity): PER1, NR1D2, PER3, NR1D1, PER2, ARNTL, NPAS2, CLOCK, CRY2 and CRY1 as well as other top genes shown to predict circadian rhythmicity in human blood (DDIT4, CLEC4E, FKBP5, DAAM2, TPST1, IL13RA1, SMAP2, HNRNPDL, FOSL2, PER1, FLT3, CDC42EP2, TMEM88, NR1D2, RBM3) in ref. 33 was used to take a focused look at results of the gene-based analysis. Using a Monte-Carlo approach, the probability of obtaining at least the observed number of significant associations (p < 0.05 uncorrected) in a random gene set of the same size was calculated (see SI, Fig. S1a-c).

Gene set enrichment analysis

Gene Set Enrichment Analysis (GSEA)34 was used to determine whether differentially expressed genes offered biologically meaningful insights about SD. Ranked lists were created based on results of models M1-M3, using a signed log10 transformed p-value with sign denoting direction of change, as described elsewhere35. To allow a concise interpretation of the potentially widespread effects of SD interventions, the Hallmark gene set collection36 (MSigDB Version 6.2), comprising 50 gene sets representing specific well-defined biological states/processes displaying coherent expression, was used. The heme metabolism gene set was excluded due to the globin interference artefact (communication with ThermoFisher).

Longitudinal Gene Set Trajectory

Time-course Gene Set Analysis (TcGSA)31 (Version 0.12.1) was employed to examine gene expression dynamics over all three time points. TcGSA, employing mixed models, detects gene sets in which expression changes over time, taking between-gene and individual variability into account, with higher sensitivity and better interpretability than univariate individual gene analysis (for details, see ref. 31). As above, the Hallmark gene set collection was used. Significance was set at FDR q < 0.05. TcGSA was employed separately for all patients, controls, responders and non-responders. Data at T3 from two patients not following the recovery sleep protocol were excluded, while 8 participants had dropped out.

Results

Gene-based analysis

In all patients [M1], 4071 (2083 up, 1988 down) genes were significantly differentially expressed after SD. Significant differential changes in gene expression after SD between responders and non-responders [M2] were observed in 360 genes [150 up, 210 down] and in patients vs. controls [M3] in 495 genes [248 up, 247 down].

Table 2 shows top differentially expressed genes for these models. Tables S1.1-S1.3 show the number of genes differentially expressed, upregulated and downregulated for all models, and detailed lists of differentially expressed genes.

Table 2 Top 10 differentially expressed genes for each model (T2 vs. T1).
Full size table

At baseline, no genes were significantly differentially expressed between patients/controls and responders/non-responders (Table S1.4).

Target examination of differential expression of circadian genes

We observed significant differential expression of circadian genes in models M1, M2 and M3. (Table 3, for more details, see Tables S2.1–3). Baseline differences in circadian genes between responders/non-responders, and patients/controls did not reach FDR q < 0.05, but achieved nominal significance (see Table S2.4).

Table 3 Differential expression in circadian genes after SD.
Full size table

Monte Carlo simulations confirmed that significantly more circadian genes were differentially expressed than random gene sets of the same size (see SI).

Gene set enrichment analysis

GSEA notably found enrichment in immune response related pathways (see Tables S3.1–3).

For M1, M2 and M3, 12, 23 and 11 gene sets were significantly positively enriched, respectively (FDR q < 0.05). The TNFα signalling via NFKβ gene set had the strongest positive enrichment in all models, while Inflammatory Response was also consistently among the significantly positively enriched gene sets.

Given the enrichment observed in Tumor Necrosis Factor Alpha (TNFα) and immune pathways, selected genes prominent in immune processes26 were further examined (Table S4).

Time-course gene set analysis

TcGSA results mirrored and extended gene-based analysis results:

In responders, 48 gene sets varied significantly (the model for one gene set, G2M Checkpoint, did not converge) (see Fig. 2a). Descriptively, in comparison to T1, responders showed a spectrum of differential gene expression at T2, with strong upregulation observed (TNFα Signalling via NFKβ, IL6-JAK-STAT3-Signaling, Inflammatory Response, Angiogenesis) and maintained until T3. The Interferon Gamma Response and Interferon Alpha Response gene sets showed the strongest upregulation at T3.

Fig. 2
figure 2

Heatmaps of estimated dynamics from significant gene sets in a responders, b non-responders, c patients, and d controls. The median gene expression over subjects is used for each trend. Each trend is zeroed at T1 to represent baseline expression. Each row is a group of genes having the same trend inside a gene set, while the columns are time points. Trends are hierarchically clustered. Trends become red as median expression is upregulated or blue as it is downregulated compared to the baseline value at T1. The colour key represents the median of the standardized estimation of gene expression over the group of participants for a given trend in a significant gene set. Non-converging gene sets were excluded. For non-responders, the gene set for peroxisome was excluded from visualisation due to non-homogeneous expression within the set.

Full size image

In non-responders, 44 gene sets varied significantly (the model for one gene set, Allograft Rejection, did not converge) (see Fig. 2b), the majority of gene sets were downregulated at T2, with upregulation of the abovementioned gene sets only observed at T3.

In patients, expression was observed to vary significantly in 49 gene sets (see Fig. 2c). Descriptively, upregulation was observed in immune system related gene sets in two main patterns; (1) strong upregulation at T2, sustained but weakening at T3 (i.e. Inflammatory Response, IL6-Jak-Stat3-Signalling and TNFα Signalling via NFKβ) and (2) light upregulation at T2 followed by stronger upregulation at T3 (e.g., Interferon Gamma Response, Interferon Alpha Response).

In controls, 21 gene sets varied significantly (see Fig. 2d). Descriptively, at T3, immune/inflammation related gene sets were upregulated (e.g. Inflammatory Response, IL6-Jak-Stat3-Signalling, TNFα Signalling via NFKβ). The majority of other gene sets were downregulated at T2, followed by upregulation at T3.

Tables S5.1–4 list significant TcGSA gene sets.

Discussion

Here, we report the first longitudinal transcriptome-wide study of SD treatment in major mood disorder patients suffering from a major depressive episode. Widespread differential gene expression was observed after SD; circadian genes were differentially expressed, and enrichment in pathways related to immune function, inflammatory response and sleep regulation was observed.

It has been hypothesized that SD exerts its antidepressant effect by resetting disturbed clock gene functioning thought to be a feature of depression10. Consequently, we examined changes in expression of circadian genes. All individuals, i.e. patients and controls, showed significant differential expression of circadian genes after SD. In responders, many significant gene expression changes were observed while in non-responders only a few nominally significant changes were found. The following observation may be of interest: three circadian genes, PER1, CLEC4E and SMAP2 (FDR q = 0.083 approaching significance, 0.048, and 0.039, respectively) showed strong differential gene expression, i.e. increased expression, in responders versus non-responders after SD. While the roles of CLEC4E and SMAP2 are yet unexplored in this context (however, see below for additional discussion of CLEC4E), similar PER1-related findings have been previously reported; consistent with the present findings, PER1 expression is observed to be increased in SD responders and decreased in non-responders10. Animal studies have shown that sleep deprivation or prolonged wakefulness enhanced Per1 expression in several brain regions21,37, and that quetiapine increased Per1 expression in the amygdala38.

The present results are interesting in light of a study in human post-mortem brain tissue of ~12,000 ranked genes according to robustness of circadian rhythmicity across six brain regions; the top ranked circadian genes in that study were ARNTL, PER2, PER3 and NR1D1 and dysregulation in MDD vs controls were observed in these genes22. The present study observed that these same genes were the most affected (of the traditional clock genes) as a result of SD, supporting the idea that SD may act on the dysregulation of these genes in depression.

The present findings cannot yet demonstrate that clock gene dysregulation/normalization is at the core of the SD mechanism. However, they suggest that (genes of) the molecular clocks of responders/non-responders, as well as patients/controls respond differently to chronobiologic stimuli, which may be associated with treatment outcomes.

Among the top 10 genes most significantly differentially expressed after SD in patients [M1] were: (i) KLF6, observed to be upregulated in individuals after experimental restriction in a genome-wide association study of sleep duration39; (ii) SIPA1L1, expression of which was found to be increased in a study examining changes in military personnel at baseline and after improvement of sleep40; (iii) NHSL2, the function of which is unknown, but which is located in a genomic region found to be differentially hydroxymethylated in a study of sleep deprivation in mice41; (iv) TREM1, which has immune function12 and (v) TSPAN2, of which a study looking for gene expression changes under fluoxetine in rats found hippocampal upregulation42. Among the top 10 genes observed in responders vs non-responders [M2] is BECN1 (Beclin1)—which (together with FKBP5, see below) is shown to be involved in priming autophagic pathways to set the stage for antidepressant action43,44. Circadian rhythms of autophagic proteins, including beclin1 have been linked to sleep disturbances45. Also in the top 10 in M2 were DNER, and GSR which have functions related to immune response46,47. Between patients and controls [M3], top genes included EZH2, which is reported to have a close relationship to IL-648. ZBTB16 is implicated in human sleep duration39 and sleep deprivation in animal studies49, and CISH is shown to be involved in immune-related processes50,51. Among other differentially expressed genes are ones evidently associated with antidepressive intervention, such as FKBP552 where we observed significant upregulation in patients vs. controls after SD (see Tables S1.1, S2.1). FKBP5 plays an important regulatory role in stress response52,53,54,55; circadian rhythm abnormalities have been linked to the stress response system, and circadian clock genes can both regulate and be regulated by rhythms of glucocorticoid release56.

Pathway analyses showed a global effect of SD on gene expression; immunological, inflammatory response and sleep regulation involved pathways were most strongly affected. These findings are of interest in light of previous gene expression studies linking immune function to MDD57,58; SD may affect pathways involved in the MDD etiology. In patients and especially responders, in the Inflammatory Response (genes related to cytokines, growth factors, cell differentiation markers and transcription factors)36, IL6-JAK-STAT3-Signalling (aberrantly hyperactivated in patients with cancer and chronic inflammatory conditions)59, and TNFα Signalling via NFKβ (cell proliferation, differentiation, apoptosis, neuroinflammation mediated cell death) gene sets, strong upregulation was observed at T2. Prior findings have shown associations between the immune system and depression, suggesting that causal pathways exist from immune dysregulation and inflammation to MDD24,25,60,61,62,63. This raises the question of how SD and the immune system interact and whether SD counteracts and/or enhances depression-related immunological processes.

Interestingly, and in contrast to responders, mainly weaker downregulation was observed at T2 in non-responders, with upregulation of immune-related gene sets only observed at T3. This differential function of immune-related processes (i.e. blunted system responsiveness) might be associated with non-response to treatment. This preliminary finding must be further explored in larger sample sizes.

Of note, downregulation of MYC targets pathways after SD was consistently observed in both responders and non-responders but not controls. MYC, well known as an oncogene64, acts as a mediator and coordinator of cell behaviour which also inversely modulates the impact of the cell cycle and circadian clock on gene expression65. It has been shown that dysregulation of MYC disrupts the molecular clock by inducing dampened expression and oscillation of CLOCK-BMAL1 master circadian transcription factor66. Removal of clock repression by MYC may play a role in the interaction of SD with depression and further investigation of these potential mechanisms is warranted.

Sleep and immunity are connected by anatomical and physiological bases67. The role of cytokines in the brain is complex and remains to be fully understood; while a full discussion of their roles is beyond the scope of this work, Table S4 shows cytokines and substances which have been implicated in both immune and sleep regulatory processes26. The present findings, observing upregulation after SD, support reports linking levels of pro-inflammatory cytokines to depression and suggesting involvement in disease pathogenesis26,68,69.

One important clue to SD response may lie in the observed TNFα expression patterns. TNFα is a pro-inflammatory cytokine controlling expression of inflammatory genetic networks; in addition to many immune-related functions in the brain, it influences whole organism function, including sleep regulation26,61. After SD, patients had significantly decreased expression of TNF compared to controls (FDR q = 0.01); upregulation was observed in controls and downregulation in patients. It should be noted here that the sample sizes were unbalanced, potentially introducing bias in the result; however, upregulated TNF in controls is consistent with reports of sleep deprivation-induced increases in TNFα levels in healthy people70. The decrease observed in patients, and especially responders, may inform the mechanism of SD response-increased TNFα concentrations are reported to be a marker of depression and TNFα administration is reported to induce depressive symptoms71.

Sleep-wake cycles are accompanied by changes in circulating immune cell numbers and disturbing the circadian cycle affects immune response72,73. Acute SD reportedly affects diurnal rhythmicity of cells, mirroring the immediate immune stress response74,75. Consistent with a stress-like response, upregulation of expression of IL6 and IL1B, thought to be involved in depression76 and regulation of circadian rhythm/sleep homeostasis13,69,77, was observed after SD. Given the tight coupling between circadian and immune systems, it may be that like the ‘resetting’ of circadian genes, the altered rhythmicity of the stress immune response system in depression is transiently normalized, leading to antidepressant effects; chronotherapeutic approaches that extend SD effects may be prolonging this normalization21.

In support of the idea of involvement of stress response associated with SD is the fact that several of the differentially expressed genes are known to be driven by glucocorticoid signalling. One possibility is that differential expression observed has been induced by stress responses associated with SD; it is not possible in the present protocol to disentangle stress and SD. For example, while as indicated above CLEC4E is not explored in the SD context, it is known to be involved in immune function and also to be driven by glucocorticoids and is thus related to the stress response78. Animal studies have shown that by removing glucocorticoid signalling via adrenalectomy, PER1, for example, was no-longer affected by SD, and responds to the stress induced by SD rather than SD itself79.

There are several limitations to our study. The sample sizes of non-responders and controls were limited, with less power potentially contributing to the lower number of significant genes. The differences in sample sizes might furthermore create bias in the analyses of differential gene expression. While the statistical approach (linear mixed models) used in the single-gene analysis and underlying TcGSA is robust against unbalancedness and missing data, caution should still be exercised when considering the comparative results, as mentioned above. Next, although blood was sampled once per day on three consecutive days, a significant but manageable load for patients, the amount of data is sparse for statistical estimation. To better leverage longitudinal data and methods, a denser sampling scheme (e.g. every few hours) will be required to attain a more refined understanding of underlying SD mechanisms. It should be noted that the recovery sleep episode was a total of 7 h in length (i.e. 1800–0100 h) after a 36 h homeostatic buildup of sleep pressure; this may not be long enough for homeostasis to recover and thus SD could have still had an effect on gene expression at T3. However, it is still unclear how circadian and homeostatic factors interact to regulate sleep processes80 and to address this in more detail, a different experimental design would have been required. Also, while food intake was not given special emphasis in this study, timing of food intake can be an entrainment signal for the circadian clock81. In the present study, some participants had small or simple meals (i.e. bread) and food intake was not restricted; possible effects cannot be excluded in the present design. Another direction, which would be informative about circadian phase, would be to include a marker of the central clock phase, e.g. melatonin secretion, also sampled at high temporal resolution. The sample studied was a naturalistic sample recruited from consecutive admissions and treated following standardized clinical guidelines. Medication regimens were tailored to the individual based on specific need, resulting in a variety of therapies used. While the effects of particular medications were not examined, to control for effects on results of the study, as an inclusion criterion, it was stipulated that the patient had to be on stable medication for at least 5 days prior to SD. Given the variety of medications used it was not possible to test for associations with drug response with sufficient statistical power. However, there was no apparent difference in substance class across response status. Nevertheless, robust effects of SD were observed, perhaps in part attributable to the consistent methodology for applying SD treatment in a relatively large cohort for SD. Considering the fact that depression is a heterogeneous phenotype accompanied by a heterogeneous immunophenotype61, means that even larger sample sizes will be needed to substantiate the present findings. Finally, we investigated peripheral tissue transcriptome-wide expression changes; although an easily obtained, valuable proxy, the correlation of expression in the blood with expression in the brain, where depression is thought to act, is imperfect and requires further study82. On the other hand, it is precisely these cells that best represent the current status of the immune system and the inflammatory response, the gene expression of which appears at the centre of SD. In addition, the crosstalk of the organs and biorhythm tuning is also mediated via the blood system.

Our findings affirm and emphasize the close relationship between circadian, immune and sleep systems in depression at the transcriptomic level, but the directionality of cause-effect remains unclear. Circadian, immune and sleep dysregulation may precede, accompany or come as a result of depression; they now represent targets for treatment, which have the ability to influence clinical outcomes. Closer investigation of these systems with larger sample sizes and denser, longer-term sampling schemes will be key to disentangling and understanding the multi-level interactions occurring.