Differential DNA methylation in recovery from shift work disorder

The human DNA methylome is responsive to our environment, but its dynamics remain underexplored. We investigated the temporal changes to DNA methylation (DNAme) in relation to recovery from a shift work disorder (SWD) by performing a paired epigenome-wide analysis in an occupational cohort of 32 shift workers (25 men, age = 43.8 ± 8.8 years, 21 SWD cases). We found that the effect of vacation on DNAme was more prominent in the SWD-group as compared to controls, with respect to the amount of significantly differentially methylated positions (DMPs; Punadj < 0.05) 6.5 vs 3.7%, respectively. The vast majority (78%) of these DMPs were hypomethylated in SWD but not in controls (27%) during the work period. The Gene Ontology Cellular component “NMDA glutamate receptor” (PFDR < 0.05) was identified in a pathway analysis of the top 30 genes in SWD. In-depth pathway analyses revealed that the Reactome pathway “CREB phosphorylation through the activation of CaMKII” might underlie the recovery. Furthermore, three DMPs from this pathway, corresponding to GRIN2C, CREB1, and CAMK2B, correlated with the degree of recovery (Punadj < 0.05). Our findings provide evidence for the dynamic nature of DNAme in relation to the recovery process from a circadian disorder, with biological relevance of the emerging pathways.


Scientific Reports
| (2021) 11:2895 | https://doi.org/10.1038/s41598-021-82627-0 www.nature.com/scientificreports/ ameliorate, and normal sleep-wake function should be restored. As many as one-third of shift workers develop SWD [25][26][27] , thus making it a common condition that negatively affects work performance and disrupts family and social life [28][29][30] . Importantly, the vulnerability to experience SWD symptoms varies remarkably between individuals: some shift workers tolerate and recover from shift work better than others [31][32][33] , which is at least partially likely to result from genetic influences 34 . Despite a growing interest to study the biology of tolerance and adaptation to shift work, little is known about the molecular processes that underlie recovery from shift work among those with SWD. In order to explore the dynamic nature of DNAme and to investigate biological processes underlying recovery from SWD, we conducted enrichment analyses of the top findings of the methylation changes for the individuals suffering from SWD. We performed a paired epigenome-wide analysis of DNAme and, for the first time, evaluated the significance of a vacation-specific effect for each CpG in the Illumina 450 K array, separately for the SWD group and the controls. In post-hoc analyses, we assessed the association between the changes in DNAme and the degree of recovery and reported SWD-specific DNA methylation changes that may regulate the restoration of a healthy state.

Results
Methylome-wide analyses of the paired Airline data. Our study focused on 32 shift workers (21 cases with SWD, 11 controls, 7 women, age = 43.8 ± 8.8 years) with DNA samples collected during a working period and after a vacation. Based on the estimated change in recovery symptoms for each shift worker, we divided SWD individuals into three groups: well-recovered, recovered, and poorly recovered (See Table 2 in the Supplementary Methods for the characteristics of the SWD participants in the groups).
To estimate the effects of being on vacation, as opposed to working on DNAme, we performed EWAS employing a linear regression model (see "Methods" and Supplementary Methods). The effect of vacation was stronger in the SWD group as compared to controls. In SWD, 6.5% (28,419/433,479) of all CpGs were DMPs using an uncorrected P < 0.05 (469 CpGs with an FDR-corrected P < 0.3). In controls, 3.7% (16,056/433,479) of the CpGs were DMPs using an uncorrected P < 0.05 and none of the CpGs remained significant after the multiple testing correction procedure. The resulting quantile-quantile (QQ) plots for all models are shown in the Fig. 1. For  We next used the same linear regression model to estimate the effect of vacation in the three recovery groups of SWD individuals separately. Of the DMPs, only 397, 336, and 6,090 (less than 1.5%) were significant in the poorly recovered, recovered, and well-recovered SWD individuals, respectively. However, at 83% (331/397), 88% (299/336), and 79% (4,836/6,090) for the three groups, respectively, the DMPs were hypomethylated during work, as compared with vacation, similar to the overall trend of the SWD group.
Considering the prominent effect of vacation on methylation in the SWD group, we focused on these findings in analysis of biological relevance: gene ontology and exploration of the vacation-sensitive pathways. For the subsequent analysis, we opted to pool the SWD individuals, due to a limited power to analyze sub-groups separately and due to increased hyper-vs hypomethylation as a result of vacation. We selected the top 30 genes (see Supplementary Table S1) ranked by significance (unadjusted P value < 1E-5, FDR-corrected P value ≤ 0.16) and subjected this set of genes to gene set enrichment analyses.
Gene ontology enrichment of the vacation-sensitive gene set. Gene ontology analyses. To identify the affected biological pathways, we explored three gene set ontology (GO) libraries in Enrichr 35,36 : GO Biological Process, GO Molecular Function, and GO Cellular Component 2018. The top finding from the GO Cellular Component 2018, "NMDA selective glutamate receptor complex" (GO:0017146), was found to be statistically significant at P < 0.05 even after correction for multiple testing (P = 8 × 10 -5 , Benjamini-Hochberg (BH) -corrected P = 0.03). The top terms for GO Molecular Function 2018 and GO Biological processes 2018 were "NMDA glutamate receptor activity" (GO:0004972, P = 4 × 10 -5 , BH-corrected P = 0.05) and "Excitatory chemical synaptic transmission" (GO:0098976, P = 4 × 10 -5 , BH-corrected P = 0.23), respectively. The findings from GO enrichment analyses are summarized in Supplementary Table S2. Thus, the analyses revealed enrichment of biological processes related to the NMDA glutamate receptor complex in all three ontology libraries.
In-depth analysis of the vacation-sensitive pathways with reactome. We further explored the pathways involving the vacation-sensitive DMPs by using Reactome. The hierarchical structure of Reactome allows the dissection of the signals into detailed molecular processes 37 . Analysis of the top 30 genes using Reactome 2016 resulted in similar findings as those with GO, adding further evidence for the involvement of changes in glutamate transmission in relation to recovery from SWD (Supplementary Table S3 and Table 1). Furthermore, the results indicated that six out of the ten top findings (P value < 0.01) are hierarchically organized under the same parent term "Activation of NMDA receptor upon glutamate binding and postsynaptic events". This www.nature.com/scientificreports/ parent pathway is also included in the top-ranked pathway, "CREB phosphorylation through the activation of CaMKII" (Supplementary Fig. S1).
For the in-depth pathway analyses, we included CpGs corresponding to all genes listed in each of the Reactome top 10 pathways (genes: n = 5-39 per pathway and CpGs: n = 73-1136 per gene, see Table 1). To investigate the DNAme changes from the paired EWAS study (work vs. vacation), we calculated the DMP enrichment and hypomethylation scores in the CpGs from each pathway both in the SWD and the control groups (Table 1). In the SWD group, we found the highest enrichment score (8.3%, 38/458) for the pathway "CREB phosphorylation through the activation of CaMKII" (CREB-CaMKII pathway). This pathway also showed hypomethylation in 74% of the DMPs (28/38) (for details, see Supplementary Table S4). In controls, the enrichment and hypomethylation scores for this pathway numbered 4.1% (19/458) and 26.3% (5/19), respectively. The QQ plots for CpGs from this pathway are displayed in Fig. 2a (SWD group) and 2b (controls), which show that the significant DMPs are frequent in cases but not in controls. Thus, our results provide evidence for a relatively robust involvement of CpGs from genes of the CREB-CaMKII pathway specific to the SWD group. www.nature.com/scientificreports/ In order to visualize the magnitude, direction, and significance of the DNAme changes induced by vacation for this pathway, we utilized a scatter plot ('volcano plot') visualization (Fig. 2c). The largest changes in hypermethylation, along with the highest statistical significance in SWD group were observed for the DMPs corresponding to the genes encoding NMDA receptor subunits GRIN2C and GRIN2D. We also discovered DMPs corresponding to 13/15 genes (GRIN2A, GRIN2B, GRIN2C, GRIN2D, CAMK2A, CAMK2B, CAMK2D, CAMK2G, CREB1, CALM1, DLG4, ACTN2, AKAP9) in this pathway in the SWD group, suggesting that virtually all genes in this pathway had one or more CpGs affected by the vacation.

In-depth analysis of DMPs from the CREB-CamKII pathway. Correlations of vacation M-values with
the recovery in SWD group. We further investigated whether the changes in the methylation values (M-values) of the 38 DMPs from the CREB-CaMKII pathway correlate with the degree of recovery in the SWD group based on data from on-site questionnaires. We found 3/38 DMPs to correlate at P < 0.05 with the degree of recovery, none of which remained significant after multiple testing correction (Table 2 and Supplementary Table S5). The change in M-values correlated inversely with the degree of recovery for cg13823003 (GRIN2C), cg05019488 (CREB1), and cg18848222 (CAMK2B), as all three sites showed hypomethylation during working period as compared with vacation. Figure 3 illustrates the biggest absolute change in the M-values of individuals with the highest degree of recovery.
SWD vs. controls during work. We then compared the methylation levels of the 38 DMPs in the SWD group and the controls during the working period. In the SWD group, 28/38 of the sites were hypomethylated and, respectively, 10/38 were hypermethylated, during the working period, as compared with vacation. The majority of the hypomethylated (20/28) sites showed lower M-values in the SWD group as compared with the controls; accordingly, the majority (9/10) of the hypermethylated showed higher M-values in SWD (one site, cg13823003 from GRIN2C, with a nominally significant difference at P < 0.05, analysis of variance (ANOVA) (Supplementary Table S6).
The results of the in-depth analysis of the CREB-CaMKII pathway are summarized in the Table 3. Of all the DMPs identified, the findings with cg13823003, corresponding to gene GRIN2C, were the most robust and exhibited all the following properties: a nominally significant P value < 0.05 in the secondary analyses of recovery; ANOVA analysis indicated an altered behavior during the working period; and the site was among the 30 most significant sites for the pathway analysis.
Global DNAme patterns across the recovery groups of SWD group. In order to explore global DNAme patterns across the recovery groups of SWD individuals and to compare them with DNAme patterns of Table 2. Correlation between the methylation changes with the changes in recovery symptoms: the three identified DMPs and the corresponding genes. a P value, uncorrected P value, 2-tailed; q value, FDR-corrected P value. www.nature.com/scientificreports/ controls, we performed a t-SNE mapping of the global DNA methylome profiles (Fig. 4). Globally the methylation profiles of the well-recovered individuals resembled the profiles of the controls, while the poorly recovered and recovered individuals cluster separately at specific regions of the DNA methylome space.

Discussion
In order to gain insight into the dynamic nature of DNAme, we explored the effect of vacation on DNAme compared to a working period in shift workers. Our hypothesis was that among the shift workers with SWD, recovery from insufficient sleep would be reflected in the changes of DNAme patterns in white blood cells. This study used paired data from an occupational cohort of shift workers, with and without SWD, for genome-wide analysis of DNAme from whole blood samples, in combination with information from sleep diaries and questionnaires collected during work and after vacation. To the best of our knowledge, this is the first study to longitudinally examine dynamic DNAme changes in relation to a manifesting sleep disorder and recovery from it.
The main finding of this study is that the effect of vacation on DNAme was more prominent in the SWD group, as compared to the controls, with dynamic changes observed in 6.5% of the DMPs. The observed changes included a gain in methylation induced by vacation, which was more remarkable for the individuals suffering www.nature.com/scientificreports/ from SWD, as compared with the controls. An overall trend towards global DNAme hypomethylation was earlier reported in several epigenome-wide association studies conducted in the cohorts of shift workers [38][39][40] . Moreover, this finding is in accordance with our previous study 12 , where we observed a loss of methylation in individuals suffering from insufficient sleep. In the current study, we further confirmed that insufficient sleep during work is associated with hypomethylation, while the recovery from SWD symptoms during vacation corresponded to the restoration of methylation. Gene enrichment analyses of the DMPs which showed the largest changes after vacation in the SWD group revealed pathways related to the activity of a glutamatergic NMDA receptor. These pathways included GRINsubunits of NMDA receptor, several Ca 2+ /calmodulin-dependent protein kinases (CaM kinases), and the transcription factor CREB1. This study's focus on longitudinal observations of the shift workers enabled us to narrow down a wide group of pathways to the specific molecular process involving only a few genes. In depth analysis revealed the Reactome pathway "CREB phosphorylation through the activation of CaMKII" to have the highest enrichment of methylome changes. This pathway includes the genes encoding for cAMP responsive element binding protein 1 (CREB1), calcium/calmodulin dependent protein kinase II beta (CAMK2B), and glutamate ionotropic receptor NMDA type, subunit 2C (GRIN2C).
Recovery dynamics in shift workers in the SWD group were found to be highly individual, with some shift workers reporting full recovery, while others still mentioned sleepiness and insomnia even after the vacation. We found that the largest changes in recovery symptoms were accompanied by the largest DNAme changes of three CpGs from the CREB-CaMKII pathway. One of them, a CpG site cg13823003 from GRIN2C, was remarkably robust through our study, showing statistical significance both in ANOVA analysis during the working period and in the secondary analyses of recovery. This site was also among the top DMPs obtained from the paired EWAS. According to the Ensembl Regulatory Build in Ensembl (http://www.ensem bl.org/info/genom e/funcg en/index .html) and Illumina 450 K annotation, cg13823003 locates at the CpG Island in the GRIN2C promoter and CTCF-binding site. In the same CpG Island in the TSS200 region, we found two adjacent CpG sites, which were differentially hypermethylated during work in SWD cases (See Supplementary Table S7 for details). Interestingly, one of these sites, cg09722397, was identified earlier in our previous study in a sample of adolescents suffering from depression and comorbid insomnia 41 . In coherence with the present study, where this site was hypermethylated in SWD cases during work, it was hypermethylated also among the adolescent cases, as compared to age and gender-matched healthy controls. However, the opposite direction of methylation observed in these three CpGs located at the promoter of GRIN2C (two hyper-and one hypomethylated) presents a challenge for an attempt to propose the possible consequences for gene expression.
The members of the CREB-CAMK pathway that we identified in the present study belong to the Ca 2+ -dependent hyperpolarization pathway, which has been implicated in the regulation of sleep duration 42 . Several other studies exploring knockouts of the genes activating the Ca 2+ -dependent hyperpolarization pathway resulted in the same conclusion: the activity of these genes is associated with the cortical capacity to evoke slow-wave oscillations, which, in turn, is proportional to animal sleep duration [43][44][45] .
As our study was carried out in the DNA samples from peripheral blood, we explored what is known about possible functions of such key members of CREB-CAMK pathway, as GRIN-subunits of NMDA receptor, CaM kinases, and CREB1 in the white blood cells. www.nature.com/scientificreports/ Human lymphocytes express various receptors, including NMDA glutamate receptors. The activation of these receptors may be involved in lymphocyte survival and differentiation of T-cells, and modulations of immune functions, as suggested by studies using NMDA receptor-binding, structural glutamate analogues 46 .
CRE elements have been identified in the promoters of many T-cell genes, indicating that transcription factor CREB might be involved in processes, such as activation and proliferation of T and B cells, as well as monocyte survival 47 . Less is known about calcium-calmodulin pathway, though some rodent studies evidenced involvement of CaM kinases 48 and calmodulin in the cytokine production, T cell activation and proliferation 49 . Based on our findings from the white blood cells, it is possible to suggest that the vacation-induced changes in the activity of some members of CREB-CAMK pathway might be regulated via DNA methylation. The observed alterations in DNAme of GRIN-subunits, CREB, and CaM kinases may indicate immune function-related differences between individuals with SWD and controls. Several previous studies have shown an association between continuous sleep insufficiency, common among the shift workers, and alterations in host defence and immunity [50][51][52][53][54][55] .
Our analysis of the global DNAme profiles in the recovery groups of individuals with SWD showed no clear differences in the degree of changes in poor versus well-recovered groups. This supports our findings from epigenome-wide association study and gene enrichment analyses, suggesting that the vacation did not affect DNAme in a genome-wide fashion, but affected the CpGs locally in specific pathways.
There are several limitations in this study. Firstly, the sample size was limited, and larger samples are required to detect statistically robust findings on a larger scale. However, a carefully a priori selected study design with two data collection points allowed us to detect statistically significant differences by applying a statistical model and assessing differences separately for the groups. Furthermore, the SWD and non-SWD status was assessed by both objective and subjective measurements, as in addition to the questionnaires, we performed both actigraphy monitoring and evaluated sleep diaries, as required by the International Classification of Sleep disorders-Third Edition (ICSD-3) criteria 56 . The second limitation is that in the EWAS study of the whole blood samples we did not perform adjustment for the blood cell type composition. However, prior to analyses, we estimated that our groups did not differ in their blood cell type composition, nor did this composition change from working period to vacation in the same individual. Thus, most of the significant changes in DNAme observed after the vacation are unlikely to be explained by the cell type composition. Thirdly, we lacked objective measurements of sleep insufficiency and relied on self-reported data in the post-hoc analyses. Nevertheless, it is known that there are large individual differences in shift work tolerance [31][32][33] and self-reported data allowed us to capture such variability among the shift workers and identify specific sites that are sensitive for the degree of recovery. We also lacked the information regarding the chronotype and sleep/wake habits of our subjects. The fourth limitation relates to the statistically significant differences in age observed between SWD cases and controls. Since age has a known impact on the DNAme, we included age as a covariate in our model. The fifth limitation of our study was that we could not assess the functional significance of the identified DNAme changes due to absence of corresponding gene expression data. The last potential limitation is related to the possible order effect of this study: we investigated the changes occurring from the SWD-symptoms towards the recovery, but the changes from recovered to the diseased state could be asymmetric.
Despite these weaknesses, our study is one of the few existing paired EWASes investigating short-term changes in DNAme patterns in humans. We observed prominent methylation changes within several months in the SWD group and depicted a molecular brain-specific pathway associated with these changes. Furthermore, we distinguished specific molecules that might indicate the degree of recovery of a shift worker with SWD. Further studies in larger cohorts are required to fully uncover the mechanisms and validate their health-related potential.

Conclusion
Our findings suggest that the DNAme of shift workers shows a specific response to vacation, particularly in individuals suffering from SWD symptoms. The observed differences in methylation patterns take place in processes related to NMDA receptor activity. In particular, our analysis identified the CpG sites in GRIN2C, CREB1, and CAMK2B as putatively important indicators of recovery in a shift worker with SWD symptoms. Further, our results demonstrate that even a two-week period of vacation exerts changes in the human DNA methylome -a highly dynamic and complex regulatory mechanism that we are only beginning to understand.

Methods
Study design and participants. Paired measurements. Whole blood samples, sleep diaries, and questionnaires on SWD symptoms were collected twice in the lab for each participant. The 'working period' measurement was performed during the shift worker's working period: blood samples were collected between 7 a.m. and 11 a.m. from a healthy participant, with no infection episode in the previous 7 days. On the same morning, the participant completed the questionnaire in the lab. 'Work' sample collection dates varied between May 2012 and June of 2013. The 'vacation' measurement was performed after at least two weeks of vacation (at the last day of vacation). The blood sample was collected between 7 a.m. and 11 a.m., and each participant completed the questionnaire. 'Vacation' sample collection dates varied between June 2012 and August 2013.
Participants. 32 shift workers aged 27-60 years (22% women) from a Finnish airline company were involved 31 . The work schedule during the last 12 months included regular night shifts (a minimum of 3 h of work between 11 p.m. and 6 a.m.) and/or early morning shifts (starting by 6 a.m.). www.nature.com/scientificreports/ required by the International Classification of Sleep disorders-Third Edition (ICSD-3) criteria 56 . 21 shift workers (24% women, age = 41.0 ± 8.9 years) constituted the SWD group: a participant reported symptoms of insomnia and/or sleepiness "often/continuously" during working period only and had a reduced total sleep time. 11 shift workers (18% women, age = 49.2 ± 9.1 years) lacking the symptoms were deemed as controls. The SWD status assessment is described in detail in our previous studies 12,31 .
Recovery groups in the SWD group. To assess the degree of recovery in the SWD group, we used two questions from the on-site questionnaire that measured important aspects of recovery for shift workers: (A) "How often do you not feel fresh after sleep?" and (B) "How often do you feel daytime sleepiness?". The change in recovery symptoms for each shift worker in the SWD group was estimated from the difference as follows: Tables 1 and 2 in Supplementary Methods summarize the participants' characteristics.
Infinium human methylation 450K bead chip methylation measurements. The DNA extraction, CpG methylation, and quality control procedures are described in detail in Lahtinen et al. 12 . DNA was extracted from the whole blood samples using standard methods (See Supplementary Methods for details) and methylation measurements were performed using Infinium HumanMethylation450k BeadChip (Illumina, Inc., San Diego, CA, USA). Methylation data preprocessing and quality control were performed using R package "minfi" v1.18.4 57 in R software v3.6.1 (https ://www.r-proje ct.org/). 433,479 probes passed the quality control.

Data analyses. Methylome-wide paired analysis.
To estimate an effect of being on vacation as opposed to working on the methylation of each of the n probes in the m individuals, we employed the following linear regression model: where y i, 0, j and y i, 1, j are the methylation of the i:th individual at the j:th probe when on work and vacation, respectively; h i is a flag indicating if the individual belongs to the SWD group or control (h i = 1 for a control and h i = 0 for the SWD group); v i is a flag indicating if the measurement is from a vacation period or not (v i = 0 for the working period, v i = 1 for vacation); D j are the nuisances such as age, sex, plate, alcohol consumption, and smoking status; and ε i,1,j are the errors, assumed to be independent and identically distributed zero-mean normal random variables (i.e. the effects were found using the ordinary least squares procedure). In this model, β h=0, j and β h=1, j represent the group average methylation for the SWD and control groups (at the j:th probe) when working, respectively (after removing the nuisance effects); β v-w, h=0, j and β v-w, h=1, j represent the effect of vacation (versus work) on the methylation in the two groups; and β D, j the nuisance effects. The significance of the effect of vacation (versus work) on the methylation of the j:th probe, β v-w, h, j , was tested using a variance ratio (F) test. The relevant null hypotheses considered no effect on each, either, or both groups, that is, the null hypothesis β v-w, h=0, j = 0 and β v-w, h=1, j = 0 tested for an effect on either group; β v-w, h=0, j = 0 tested for an effect in the SWD group, while β v-w, h=1, j = 0 tested for an effect in the control group; and the hypothesis β v-w, h=0, j = 0 or β v-w, h=1, j = 0 (which can be evaluated using a composite null model) tested for an effect in both groups. Finally, the acquired P values were adjusted to control false discovery rate (FDR) in multiple hypothesis testing using the Benjamini-Hochberg procedure. We defined a CpG site to be hyper-or hypomethylated based on the value of the beta coefficient in our tests, with hypermethylation characterized by a positive value. To estimate the relative proportions of cell type components in the whole blood we used a statistical method described by Houseman et al. 58 . T-test was used to investigate the differences in the cell type compositions between and within the groups.
In order to obtain the gene names associated with the CpGs, we used Infinium HumanMethylation450k Bead-Chip annotation data (Illumina, Inc., San Diego, CA, USA) and R software v3.6.1 (https ://www.r-proje ct.org/). Pathway analyses using Enrichr. Enriched terms were explored using Enrichr (https ://amp.pharm .mssm.edu/ Enric hr/) 36 . We focused on three GO libraries: GO Biological Process, GO Molecular Function, and GO Cellular Component 2018. A list of genes ranked by P values served as the input and P < 0.05 (after Benjamini-Hochberg correction) was used to determine significance.
In-depth pathway analysis using Reactome 2016. For the top 10 Reactome 2016 pathways ranked by unadjusted P values, we extracted the full gene set from Reactome 2016 library and calculated the (1) enrichment score (%) and (2) hypomethylation score (%), for both the SWD and control groups. License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.