Role of NR3C1 and SLC6A4 methylation in the HPA axis regulation in burnout

Background: Work-related stress and burnout have become major occupational health concerns. Dysregulation of HPA axis is considered one of the central mechanisms and is potentially moderated through epigenetics. In the present study, we aim to investigate epigenetic regulation of the HPA axis in burnout, by focusing on salivary cortisol and cortisone and DNA methylation of the glucocorticoid receptor gene (NR3C1) and the serotonin transporter gene (SLC6A4). Methods: We conducted a cross-sectional study with 59 subjects with burnout and 70 healthy controls recruited from the general population. All participants underwent a clinical interview and psychological assessment. Saliva samples were collected at 0, 30 and 60 min after awakening and were used to quantify cortisol and cortisone. Pyrosequencing was performed on whole blood-derived DNA to assess DNA methylation. Results: There were no between-group differences in cortisol levels, whereas burnout participants had higher levels of cortisone. Job stress was associated with increased cortisol and cortisone. We observed both increased and decreased NR3C1 and SLC6A4 methylation in the burnout group compared to the control group. Some of these methylation changes correlated with burnout symptoms dimensionally. Increased methylation in a specific CpG in the SLC6A4 promoter region moderated the association between job stress and burnout. DNA methylation in this CpG was also associated with increased cortisol. In addition, average methylation of NR3C1 was negatively associated with cortisone levels. Limitations: This is a cross-sectional study and therefore no conclusions on causality could be made. Conclusions: We provide first evidence of changes in DNA methylation of NR3C1 and SLC6A4 in burnout, which were further associated with cortisol and cortisone. Further, increased cortisol and cortisone seemed to reflect job stress rather than burnout itself.


Introduction
Work-related stress has become a major occupational health concern in industrialized countries, not only affecting the well-being of employees but also posing an enormous socio-economic and public health burden (Hassard et al., 2017). Chronically persisting work-related stress is a risk factor for developing chronic conditions such as cardiovascular diseases (Kivimäki and Kawachi, 2015) but also mental health problems such as burnout (Xanthopoulou et al., 2007), depression (Ahola and Hakanen, 2007) and chronic fatigue syndrome (Coetzee et al., 2019).
Occupational burnout in particular has received a great amount of media attention as well as increased interest in the research community (Heinemann and Heinemann, 2017). Even though there is still a certain degree of vagueness around its definition and clinical assessment (Cox et al., 2005), three symptoms are commonly accepted as the core dimensions of burnout, namely exhaustion, depersonalization (cognitive distance toward job by a cynical attitude) and reduced professional efficacy (Maslach et al., 2001). The relationship between work-related stress and burnout has been explored though various theoretical models, one of them being the Job Demands-Resources (JD-R) model . The JD-R model focuses on describing how an imbalance between job demands and job resources the individual has at disposal to cope with those demands lead to increased job strain (Demerouti et al., 2001). In the long term, job strain in turn is a major risk factor for the development of burnout. Moreover, the causal relation between the three dimensions of burnout and the evolution of its symptoms has been described through different models. The most commonly used model proposes that burnout is initially presented with high levels of exhaustion, which then leads to elevated levels of depersonalization (cynicism), in turn leading to low levels of personal accomplishment .
However, knowledge on the psychobiological processes through which work-related stress contributes to burnout development still remains limited. Identification of biological pathways linking the exposure with adverse health effects is particularly relevant for the establishment of psychobiological plausibility of the observed association (Siegrist and Li, 2017) and could help develop strategies to protect and promote well-being of employees (Kerr et al., 2020).
Since burnout develops as a response to chronic work stress, it is hypothesised that the two main stress response systems: the hypothalamic-pituitary-adrenal (HPA) axis and the autonomic nervous system (ANS) could play a role in burnout (Kudielka and Wüst, 2010). However, studies investigating HPA axis and ANS biomarkers in work-related stress and burnout gave inconclusive findings (Siegrist and Li, 2017) (Danhof-Pont et al., 2011) (Jonsdottir and Dahlman, 2019), indicating that the biological changes related to these two systems might be more complex than initially hypothesized and there are other potential mediating and moderating intermediary factors involved. As cortisol is the key product of the HPA axis activation in response to stress, the biggest part of research on burnout biomarkers focused on assessing this hormone (Ehlert et al., 2001). Since cortisol secretion follows a circadian rhythm, a common measure of cortisol in stress-related research is the cortisol awakening response (CAR), which captures a peak in cortisol within the first hour after awakening (Pruessner et al., 1997). However, even though initial findings of blunted CAR in burnout seemed promising (Pruessner et al., 1999), they failed to be replicated as other studies also reported increased CAR in burnout or absence of any changes (Rothe et al., 2020).
In addition, since the salivary cortisol is immediately converted to its inactive metabolite cortisone, relevance of salivary cortisone has gained interest recently. Salivary cortisone was shown to better correlate with the unbound fraction of serum cortisol than salivary cortisol (Debono et al., 2018) (Mezzullo et al., 2016) and was associated with subjective and autonomic stress measures in healthy individuals (Bae et al., 2019). Still, there are no studies investigating the association between salivary cortisone and work-related stress and burnout, to the best of our knowledge.
Another possible reason for inconsistency in findings on HPA biomarkers is the potential involvement of other more stable molecular mechanisms in its dysregulation, which could be more reflective of chronic stress than acute biomarkers such as cortisol (Rothe et al., 2020). One of the main such molecular mechanisms suggested in the HPA-axis dysregulation are epigenetic markers. Studies indicate that stressful life events can lead to DNA methylation changes in stress-related genes and such changes were also observed in stress-related mental disorders, such as depression (Lee and Sawa, 2014) (Li et al., 2019). A large body of research focused on epigenetic regulation of the glucocorticoid receptor (GR) gene -NR3C1, particularly the region 1F, which overlaps with the nerve growth factor-inducible protein A (NGFI-A) binding site and therefore plays an important role in the regulation of gene expression (Palma-gudiel et al., 2015). Despite the large interest in the biological research on HPA-axis regulation in work-related stress and burnout in recent years, no studies were performed to investigate the epigenetic regulation of NR3C1 in this context. In addition, another gene that has gained attention in the epigenetic regulation of the HPA axis is the serotonin transporter gene -SLC6A4. This gene has traditionally been investigated in the context of mood (dys)regulation and depression due to its role in regulation of serotonin levels in the synaptic cleft (Booij et al., 2013). However, recent studies indicate that SLC6A4 can also play a role in dysregulation of HPA axis activity and cortisol secretion (Ouellet-Morin et al., 2012) . Interestingly, changes in DNA methylation of this gene were previously identified in nurses exposed to high stress working conditions (Alasaari et al., 2012). However, no studies on clinically diagnosed burnout population have been performed.
With the present study, we aim to bridge the existing literature gaps, by focusing on two main objectives. First, we aim to compare DNA methylation patterns of stress-related genes (NR3C1 and SLC6A4) between people with burnout and healthy controls and test whether DNA methylation of these genes is associated with job stress and self-reported burnout symptoms. Second, we aim to investigate the underlying epigenetic mechanisms of HPA axis dysregulation in burnout, by measuring salivary cortisol and cortisone and their link with DNA methylation.

Study population
For the present study, we recruited a well-characterised burnout group and a healthy control group, in a cross-sectional design. A detailed description of the recruitment procedure, inclusion and exclusion criteria has been previously described . Briefly, information about the study was disseminated using media channels and flyers, after which all potential candidates were asked to fill in an online screening tool with information about the burnout diagnosis and comorbidity. People who met inclusion/exclusion criteria in the screening were invited for a clinical interview with a psychologist, where burnout was assessed using the Dutch practice guidelines for managing adjustment disorders in occupational and primary health care (van der Klink and van Dijk, 2014). Moreover, all participants were screened for major depressive disorder and general anxiety disorder, by using MINI (Sheehan et al., 1998) and those fulfilling criteria were excluded. Finally, 129 subjects were included in the study (N=59 with burnout and N=70 healthy control subjects). All subjects were Caucasian.
This study was approved by the commission for Medical Ethics of the UZ Leuven (S59567) and all subjects signed an informed consent priory to inclusion in the study.

Psychological assessment
A clinical interview was performed to assess burnout, inclusion and exclusion criteria and collect data on variables related to burnout (e.g. medication use). To assess burnout dimensionally, all subjects were asked to fill in a validated Dutch version of the Maslach Burnout Inventory -General Survey (MBI-GS) (Maslach et al., 1996), which is the Utrecht Burnout Scale-A (UBOS-A) (W. Schaufeli and Dierendonck, 2001). MBI-GS is intended for burnout assessment in all occupations, and includes the following dimensions: exhaustion, cynicism and reduced professional efficacy. To assess job demands, job resources and job stress, we applied the Short Inventory to Monitor Psychosocial Hazards (SIMPH), which has been previously validated in the Belgian context (Notelaers et al., 2007). To measure job demands, we used the SIMPH scales workload, emotional load and role conflict. To measure job resources, we applied the scales social support, variety, skill utilization, autonomy, participation and role ambiguity (Vandenbroeck et al., 2017). A cumulative job demands and job resources score were calculated as the sum of individual scales. Balance between job demands and job resources was calculated as the ratio between these two variables and was termed "job stress". To test whether the mentioned scales are grouped in two main components (job demands and job resources), we conducted principal component analysis, which confirmed the existence of the two main components (Appendix A).
Finally, a questionnaire on socio-demographic data, smoking habits and alcohol consumption was applied to assess confounding effect of these variables.

Cortisol and cortisone awakening response
All subjects were instructed to collect saliva samples at home, on a regular day using Cotton Salivettes® (Cortisol, code blue) provided by Sarstedt (Nümbrecht, Germany). To assess cortisol and cortisone awakening response, subjects were asked to collect three samples on the same day: immediately upon awakening (while still lying in bed), 30 and 60 min thereafter and note down the exact time of awakening. Participants were asked to avoid drinking, eating and teeth brushing until collecting all the samples. The importance of protocol adherence for validity of the obtained results was clearly explained and subject were asked to report any deviations from the protocol. Samples were stored at 4 • C until delivery to the laboratory, where they were immediately centrifuged at 2000 x g for 5 min, transferred to 1.5-mL tubes and stored at − 80 • C until analysis. Quantification of salivary cortisol and cortisone was performed using highly sensitive and specific UPLC-MS/MS method, previously developed and validated by our research group (Bakusic et al., 2019). Saliva samples for cortisol and cortisone measurements were provided from 108 participants (45 in the burnout group and 63 in the control group).

DNA methylation analysis
On the day of clinical assessment, blood was withdrawn into two EDTA tubes. One EDTA tube was immediately transferred to a medical laboratory (Laboratoriumgeneeskunde, UZ Leuven) for white blood cell subset analysis. The other EDTA tube was immediately stored on -80 • C until analysis. DNA extraction was performed using QIAamp DNA Blood Mini Kit (Qiagen, Hilde, Germany) and quantity and purity of DNA were determined using a NanoDrop spectrophotometer. Next, DNA samples were bisulphite-converted using the EZ-96 DNA Methylation-Gold™ Kit (#D5008, Zymo Research) according to the manufacturer's protocol. Bisulphite-converted DNA sequences of interest were then amplified using polymerase chain reaction (PCR). DNA methylation analyses were carried out by pyrosequencing, using a PyroMark Q24 instrument (Qiagen) and the results were analyzed using the PyroMark analysis 2.0.7 software (Qiagen). The primer sequences were based on those previously published in the literature (Alexander et al., 2018) (Vangeel et al., 2015)   and were designed to target the whole CpG island of NR3C1 1F region and part of the CpG island overlapping with the SLC6A4 promoter region. Samples were randomized prior to DNA extraction and a positive control (highly methylated DNA) (Qiagen) was included in each analysis.
Blood samples were missing for eight subjects and therefore DNA methylation analysis was performed on 121 subjects. Two additional samples were excluded because they did not pass the quality control, which finally resulted in DNA methylation analysis of 52 subjects with burnout and 67 healthy controls.
A detailed protocol with all analyzed amplicons, PCR and sequencing primers is provided in the Appendix B.

Statistical analysis
Socio-demographic and clinical characteristics (job stress, burnout symptomsexhaustion, cynicism and professional efficacy) of the burnout and the control group were compared using an independent sample t-test for continuous variables, and Chi-Square test for categorical variables. Cortisol and cortisone data were analysed using a mixed model for repeated measures, with cortisol/cortisone as the outcome variable, measurement (0, 30, 60min after awakening) as a fixed factor, person as a random factor and age, gender and time of awakening as covariates. Other predictors of interest were added to the model, as appropriate. Cortisol and cortisone values were log-transformed prior to analysis to reduce skewness and approximate normality.
To compare DNA methylation levels between groups (burnout vs. control) a non-parametric Mann-Whitney U test was applied because NR3C1 1F and SLC6A4 methylation data were not normally distributed. To test associations between DNA methylation and clinical variables, we performed partial correlation analysis, controlling for appropriate covariates. In case both the correlation with job stress and burnout was observed with any of the biological parameters (cortisol/cortisone levels or DNA methylation of NR3C1 and SLC6A4), we tested the moderation effect, by using the PROCESS macro of SPSS version 21.0, developed by Hayes (Hayes and Scharkow, 2013). In the moderation analysis, we tested whether cortisol/cortisone or DNA methylation moderated the association between job stress as a predictor variable and burnout (dichotomous -yes/no) as the outcome variable. Bonferroni correction for multiple testing was used where appropriate.
All statistical analyses were performed using SPSS software package, version 26.0. All tests were two-sided, and the significance level was set at 0.05.

Demographic and clinical characteristics
Socio-demographic and clinical characteristics of the study sample (N=129) are presented in Table 1  . The two groups were comparable in terms of sex distribution and smoking habits, but differed in age (burnout group was older on average). As expected, participants in the burnout group scored higher on exhaustion and cynicism scales of MBI-GS and lower on professional efficacy. In addition, burnout subjects had higher levels of job demands, lower levels of job resources as well as higher job stress (calculated as job demands/job resources ratio).

Cortisol and cortisone awakening response
Group comparison revealed no significant differences in the morning cortisol levels (B=-0.020, F=0.27, p=0.601). However, burnout participants had significantly higher cortisone levels (B=-0.73, F=6.61, p=0.011), when controlled for age, gender and time of awakening. This difference was driven by the overall higher cortisone levels rather than difference in specific time points (non-significant group*time interaction). Between-group differences in cortisol and cortisone are presented in Fig. 1 and Table C1.
In the dimensional analysis on the overall sample (Table C2), cortisol and cortisone levels were not associated with burnout symptoms (exhaustion, cynicism, professional efficacy, all ps>0.05), except from a borderline association between exhaustion and cortisone (B=0.003, F=3.80, p=0.052). In addition, both cortisol and cortisone were positively associated with job stress (cortisol: B=0.362, F=5.64, p=0.018; cortisone: B=0.324, F=8.03, p=0.005). Since there was no significant interaction effect with the group (burnout vs. control, all ps>0.05), we did not perform further group-specific associations between cortisol/ cortisone and dimensional symptom measures.

NR3C1 and SLC6A4 methylation
When comparing the total average NR3C1 methylation levels between the groups, we observed no significant difference (p=0.146). However, when looking at average methylation of specific amplicons, we observed increased methylation of amplicon 1 (p=0.028) and decreased methylation of amplicon 3 (p=0.002) in the burnout group. Comparison of individual CpGs within these amplicons after Bonferroni correction showed increased methylation of CpG21 in amplicon 1 (p=0.002) and decreased methylation of CpG30 in amplicon 3 (p=0.004) in the burnout group. Both of these CpGs overlap with the NGFI-A-binding sites. Fig. 2 illustrates between group differences in NR3C1 methylation. All the identified differences remained significant when controlled for age, gender, smoking and white blood cell count (data not presented).
Regarding the SLC6A4 gene (Fig. 3), there were no group differences in the total average methylation level (p=0.477) nor between different amplicons (all ps>0.05, all data presented in Table C3). Looking at the individual CpGs within these amplicons, participants with burnout had increased methylation of CpG8 (p=0.017) and decreased methylation of CpG16 (p=0.026) and CpG22 (p=0.010), but none remained significant after Bonferroni correction for multiple testing. Between-group comparison of average methylation levels, individual amplicons and CpGs of both genes is provided in Table C3.
To reduce multiple testing, we only kept the total average methylation of both genes and the differentially methylated regions in the between-group comparison in further analysis. For NR3C1, we kept those CpGs that were differentially methylated and overlapped with the NGFI-A-binding sites (CpG21, CpG30 and CpG31) whereas for SLC6A4 we kept all the three differentially methylated CpGs (CpG8, CpG16 and CpG22), since these regions contain multiple binding sites for different transcription factors (Messeguer et al., 2002). Consequently, we did not perform Bonferroni correction for these analyses.
As confirmatory analysis, we tested the associations between methylation and clinical measures in the overall sample, by performing partial correlation analysis, controlled for age, gender, smoking and blood cell count. In these analyses, methylation of NR3C1 amplicon 1 and 3 were associated with burnout symptoms (exhaustion, cynicism) and job stress (only amplicon 1) whereas methylation of CpG8 in the SLC6A4 was significantly associated with all burnout dimensions (exhaustion, cynicism, professional efficacy) as well as job stress (Table C4). In addition, methylation of this CpG was a moderator of the effect of job stress on burnout as the dichotomous outcome (X 2 =5.577, p=0.0182), as illustrated in Fig. 4.

Associations between cortisol/cortisone and DNA methylation
Next, we tested whether DNA methylation of NR3C1 and SLC6A4 was associated with cortisol/cortisone levels. We found that the total average NR3C1 methylation was negatively associated with cortisone in the overall sample (B=-0.066, F=6.83, p=0.009), with significant group*methylation interaction (p=0.044). Group specific analysis showed that this association was significant in the burnout group (B=-0.148, F=14.681, p<0.001), but not in the control group (B=0. 005, F=0.03, p=0.867).
Regarding the SLC6A4 gene, total average methylation was not associated with cortisol/cortisone levels (p>0.05). However, methylation of CpG8 was positively associated with cortisol levels (B=0.049, F=4.88, p=0.028), without group*methylation interaction (p=0.080). Other regions of NR3C1 and SLC6A4 were not associated with cortisol or cortisone levels (all ps>0.05, all data presented in Table C5). Associations between NR3C1 and SLC6A4 methylation and cortisol/cortisone levels are given in Table C5.

Discussion
In the present study, we aimed to investigate NR3C1 and SLC6A4 methylation in burnout and their regulation of the HPA axis, by simultaneously assessing salivary cortisol and cortisone and DNA methylation in people with burnout and healthy controls.
First, we found no group differences in salivary cortisol but we observed increased cortisone levels in the burnout group. It is worth mentioning that even though not significant, between-group differences in the cortisol level were in the same direction as effects of cortisone (increased levels in the burnout group), descriptively. Absence of changes in cortisol has been previously reported in several studies (Menke et al., 2014(Menke et al., ) (Österberg et al., 2009 (Sjörs et al., 2012). However, other studies also reported increased (De Vente et al., 2003) (Grossi et al., 2005) or decreased (Oosterholt et al., 2015) (Pruessner et al., 1999) cortisol levels, although there was a large heterogeneity in the way burnout was assessed across these studies. The inconsistency of these findings as well as the fact that burnout is a process developing over a longer period of time and involving different phases lead to the proposition of different models to explain cortisol dysregulation in different stages of burnout development. A recent hypothetical phase model suggests 4 phases in burnout developmentengaged, strained, cynical and burned-out, which are based on different levels of exhaustion, cynicism, vigor and dedication (Morera et al., 2020). This model suggests a gradual increase in cortisol levels throughout the first three phases, finally resulting in the HPA axis exhaustion and a consequent hypocortisolemia in the last phase (burned-out). The reason why we did  Fig. 4. Moderation effect of SL6A4 (CpG8) methylation on the association between job stress and burnout as the outcome variable. Burnout is expressed as dichotomous outcome (burnout vs healthy), based on the clinical assessment. Methylation slopes are presented as mean-1SD (0.89), mean (1.58), and mean+1SD (2.28). The outcome variable is expressed as log-odds (logits) or the linearized probability (odds ratios) of group membership (burnout group) derived from logistic regression. not observe group differences in CAR could be the fact that our burnout participants differed in the exhaustion and cynicism levels, however our sample was relatively small for further subgroup analysis. Interestingly, we did observe a positive association between cortisol and job stress in the overall sample, which was driven by the control group. This indeed might indicate that HPA axis hyperreactivity reflect high levels of strain and a risk of developing a burnout. Future studies on bigger sample size could further investigate cortisol patterns in people in different burnout stages.
Interestingly, we did observe increased salivary cortisone in the burnout group compared with healthy controls. Similar to cortisol, cortisone levels were positively associated with job stress, but we also observed a borderline correlation with exhaustion. To the best of our knowledge, morning cortisone levels were not assessed in burnout so far, and therefore we cannot compare these findings. However, salivary cortisone was shown to highly correlate with self-reported anxiety state and autonomic parameters (heart rate) after laboratory stress induction in healthy adults (Bae et al., 2019). In the same study, salivary cortisone had higher discriminatory power than cortisol to differentiate people who underwent a laboratory stress task from those undergoing a placebo task. These finding together with our results suggest that salivary cortisone could be an interesting marker to be considered in future studies on burnout, at least in addition to cortisol.
Next, our study is the first one to demonstrate differential DNA methylation patterns in NR3C1 and SLC6A4 between subjects with burnout and healthy controls. We observed both increased and decreased DNA methylation patterns of NR3C1. More specifically, we identified increased methylation of amplicon 1, which was mainly driven by methylation of one specific CpG, which overlaps with the transcription factor (NGFI-A) binding site and therefore can be of relevance for gene expression. On the other hand, the region of amplicon 3 showed decreased methylation in burnout, both assessed as the average methylation as well as multiple CpGs overlapping with the transcription factor binding sites. Previously, increased methylation of NR3C1 was associated with other stress-related condition, such as major depression (Nantharat et al., 2015) (Roy et al., 2017)   and was suggested as one of the mechanisms contributing to glucocorticoid resistance after prolonged HPA axis hyperactivity (Pariante and Miller, 2001). On the other hand, decreased NR3C1 methylation of the 1F region was previously observed in patients with chronic fatigue syndrome (CFS) (Vangeel et al., 2018) and veterans with posttraumatic stress disorder (PTSD) (Yehuda et al., 2015) and was further associated with higher NR3C1 expression and increased GR sensitivity (Yehuda et al., 2015). Even though we observed changes in opposite directions (both increased and decreased), the larger part of the assessed CpGs was hypomethylated and therefore our findings might suggest that burnout is, similar to CFS and PTSD, followed by increased GR sensitivity. However, GR hypersensitivity would expected to be followed by decreased cortisol/ cortisone levels, due to increased negative feedback of the HPA axis, which is not what we found in our study. In contrast, we observed increased cortisone levels in participants with burnout and a negative association between the average NR3C1 methylation and the cortisone levels in the burnout group. Unfortunately, we did not perform measures of mRNA expression/ GR sensitivity to make further conclusions on this. Another possibility is that the changes in opposite directions we captured reflect different stages of HPA axis dysregulation in people with burnout, following the idea of different stages in the HPA dysregulation elaborated in one of the previous paragraphs. However, more research on stability and evolution of DNA methylation changes in this gene in the burnout population are necessary to further explore this hypothesis.
When assessing SLC6A4 methylation, we did not observe betweengroup differences in the average levels, but we did observe increased methylation of one specific CpG. It is important to mention that this significance did not survive correction for multiple testing, however it was associated with job stress, burnout symptoms and cortisol levels, which can unlikely all be attributed to false positive results. Previously, SLC6A4 methylation was assessed in a study with nurses exposed to low and high stress environment (Alasaari et al., 2012) and decreased methylation was associated with high stress, which is in contrast with our findings. However, the authors assessed different region of SLC6A4 promoter, which makes the findings difficult to compare. Interestingly, in our previous study, we observed hypermethylation of the same CpG in patients with major depression , and so did several other authors (Li et al., 2019), which suggests that potentially common mechanisms are involved in dysregulation of the serotonergic system in burnout and depression. Moreover, in our study, SLC6A4 methylation enhanced the effect of high job stress on burnout. Since our study was cross-sectional, we cannot draw conclusions on causality, however this moderation is an interesting starting point for future longitudinal studies. Finally, increased SLC6A4 methylation was also associated with increased CAR in our study population, which is in line with the increasing number of studies demonstrating the importance of DNA methylation changes in this gene for HPA axis (dys)regulation (Ouellet-Morin et al., 2012)   (Alexander et al., 2019). The exact intermediary mechanism linking SLC6A4 methylation and HPA axis still need to be elucidated in further studies. The region of the gene where we observed DNA methylation changes contains multiple transcription binding sites of importance for the glucocorticoid system (such as GR-α) (Messeguer et al., 2002), which might present one of the links between the two systems. In other words, the increased methylation in this region might lead to lower transcriptional activation of SLC6A4 by glucocorticoids. Whereas in vitro studies demonstrated the effect of glucocorticoids on SLC6A4 expression (Glatz et al., 2003), the exact role of DNA methylation in this process still needs to be clarified.
The present study has several strengths and limitations. The main strength of the study is a good characterization of the burnout group, by using a clinical interview and strict inclusion/exclusion criteria. In previous studies, burnout was often assessed only based on self-reported symptoms and this was pointed out as one of the main limitations in biological research (Bianchi et al., 2015). In addition, we performed white blood cell count and included this as a covariate in the DNA methylation analysis, which was rarely done in previous studies even though the number of leucocytes is known to affect DNA methylation (Adalsteinsson et al., 2012). Concerning limitations, our study was cross-sectional and therefore no conclusions on causal associations could be made. In addition, the study was powered for DNA methylation analysis and therefore it is possible that some of the negative findings involving cortisol and cortisone are a result of low power, especially when assessing between-group differences. Namely, as continuous variables allow for more power, it is possible that this is the reason why we observed associations between work stress, which was assessed as a continuous variable, and cortisol/cortisone whereas we failed to see any between-group differences (burnout vs. control, dichotomous variable) in these outcomes. Finally, we assessed peripheral DNA methylation, which might not completely reflect the central nervous system. However, usefulness of peripheral markers in stress-related research has been recognized, especially when looking at DNA methylation of CpG-enriched promoter regions, which was shown to be highly correlated across peripheral and brain tissues (Davies et al., 2012).
To conclude, we found the first evidence of DNA methylation changes of NR3C1 and SLC6A4 in burnout, some of which were also associated with job stress and burnout symptoms. In addition, lower NR3C1 methylation and higher SLC6A4 methylation were associated with increased salivary cortisol and cortisone concentrations in the same population. We recommend future studies to investigate HPA axis dysregulation in different burnout phases and recovery. In addition, simultaneous assessment of potential intermediary mechanisms, such as GR sensitivity and potentially other biological endpoints such as inflammatory markers could provide more insight into the complex molecular regulation of HPA axis in burnout.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest
None.