The study consisted of a screening visit during which the diagnostic interviews and questionnaires were administered, and a test session held between 12.00 and 18.30 during which participants performed two runs of the 5 min repeated mental arithmetic stress task separated by 15 min of break. During the entire task, participants were wearing an ECG. See Figure 1 for an overview of the study design and below for detailed procedures.

Cardiac electrical activity (ECG) and respiratory movement were recorded using a 10-channel isolated bioamplifier (NeXus-10 MKII, Mind Media BV, Netherlands). Three disposable electrodes (Kendall™, H66LG electrodes) were placed in a lead II configuration on the chest to record ECG at sampling frequency of 1,024 Hz. Simultaneously, respiratory (chest) movement was recorded by placing a strain gauge belt over the abdomen, with the central part just above the navel and a sampling rate of 32 Hz was used, as instructed by the manufacturer. Correct functioning of the recordings was checked using the BioTrace+ Software (Mind Media BV, Netherlands).

The ECG signal was processed and R-wave peaks were detected using Kubios software version 2.0 (49). A visual inspection of raw data for possible movement artifacts and ectopic beats was performed per segment (i.e., each 5-min interval). Where applicable, an inbuilt artifact correction based on cubic spline interpolation was applied. Segments with more than 5% of beats in need for correction were excluded from analysis (one segment). Data was excluded for one patient showing a high number of ectopic beats (arrhythmia), and two segments of data were lost for two participants due to equipment failure. The respiratory movement signal was processed using custom-written algorithms in MATLAB software (R2018b, Mathworks, Inc., Natick, MA, United States), see Gholamrezaei et al. (50), and breathing rate was measured for each segment.

Heart rate variability analysis was performed using the Kubios software on segments with duration of 5 min according to the guidelines (51). HR, RMSSD, and high-frequency HRV (HF-HRV) were extracted for each segment. RMSSD is a time-domain HRV index reflecting beat-to-beat variation in HR, and HF-HRV is a frequency-domain HRV index reflecting variability of HR at spontaneous breathing frequency (i.e., respiratory sinus arrhythmia or RSA). Both indexes are mainly derived by cardiac parasympathetic (vagal) activity as is shown by studies using vagolytic (e.g., high dose atropine) (52–55) or vagotonic (e.g., low dose atropine) agents (56–58) to inhibit or increase cardiac parasympathetic stimulation, respectively (59). We decided to use RMSSD as the primary HRV index in our analysis for specific reasons: Evidence has shown that RMSSD is less influenced by breathing behavior, particularly the breathing frequency, compared to HF-HRV (54, 55, 60). Controlled breathing is recommended for a reliable assessment of cardiac vagal activity when using frequency-domain measures (54, 55). Therefore, RMSSD is a more suitable measure of cardiac vagal activity when controlled breathing is not possible/feasible, as was the case in the current study. Notwithstanding, to provide more comparable data with previous studies, we decided to also analyze HF-HRV while taking breathing frequency into account. Results of HF-HRV and for completeness, low frequency (LF-HRV) are provided in the Supplementary Information D.

Statistical analyses were performed in R (61) and matching of participants and controls was performed using SPSS version 25 (62). The power analysis on simulated data can be found in Supplementary Information B. We conducted a repeated-measures mixed model with random intercept, a phase (baseline, stress, recovery) × group (MDD/control) interaction and PSS × group interaction as independent variables, and either subjective stress rating, or HR/HRV as dependent variables. Assumptions of mixed models were checked by visual inspection of q–q-plots and histograms. Due to slight heterogeneity of residuals, HR and HRV indexes were transformed using the logarithmic (log) transformation for all statistical analyses. The initially planned interaction of PSS × phase to test the moderating effect was collinear with the group × phase interaction and was therefore not explored in the overall model. Nevertheless, to investigate this effect we conducted a separate model within the groups (either MDD or HC) with a PSS × phase. Significant time × group interactions were followed-up by planned contrasts adjusted for multiple testing using Benjamini–Hochberg correction between baseline and stress, between stress and recovery and between baseline and recovery of each run. Since no difference occurred between the baseline and the control condition (see Supplementary Information C), and the baseline phase was more comparable to the recovery phases, the control condition was not considered in the model.

To investigate classification accuracy of physiological parameters, we used binary logistic regression with leave one out cross validation (LOOCV) from the caret package (50).

Patients and controls were female, between 18 and 60 years old, matched for age (using fuzzy matching, tolerance up to ±5 years) with a similar BMI. Menstrual phases (contraception, post-menopausal, follicular, or luteal phase) did not differ significantly between groups (X² = 3.54, df = 3, p = 0.315). The mean HRSD score for patients was 24.63 indicating a severe form of depression (51). Controls were more likely to be non-smokers (X² = 13.55, df = 2, p = 0.001) and were significantly more active than patients (W = 634, p < 0.001). Additionally, patients had significantly higher depression severity, rumination- and perceived stress levels in the last month (all p < 0.001, see Table 1). None of the control population and around half (15 of 27) of the patients were taking SSRIS/SNRIs/SARIS (n = 11) or atypical antidepressants (n = 4). Four patients were exclusively treated with benzodiazepines.

While antidepressant medication intake was limited to SSRIs/SNRIs/SARIs and atypical antidepressants and smoking immediately prior to the experiment (at least for 1 h) was prohibited, these parameters may constitute an important confounder. We therefore tested whether the demographic variables that were different between controls and patients (i.e., medication, smoking, and physical activity levels) contributed to HR/HRV group differences at baseline. No main or moderating effect occurred for smoking (all p > 0.231) or physical activity (all p > 0.765). Antidepressant medication had no main effect on HR at baseline (t = −0.33, p = 0.748), but significantly decreased RMSSD (t = −2.46, p = 0.021).

As a manipulation check, we tested whether subjective stress induction was achieved. A significant group by phase interaction occurred (F = 3.11, df = 4, p = 0.016). Post hoc contrasts of phase by group showed the expected increase of SS in healthy controls, during both runs (compared to baseline t = 7.57, p.adj < 0.001 and t = 6.64, p.adj < 0.001), returning to approximate baseline levels during recovery (t = −1.48, p.adj = 0.154 and t = −0.37, p.adj = 0.712). Congruent with the control’s stress response, in patients, SS increased in response to the first stressor (t = 6.21, p < 0.001). However, while there was a subtle recovery after stress (t = −0.89, p.adj = 0.030), SS levels did not return to baseline during the recovery period (t = −3.95, p < 0.001) and stayed significantly elevated throughout the remainder of the experiment (all p < 0.004 compared to baseline), indicating incomplete subjective stress recovery (Figure 2 and Supplementary Table 1).

Mean levels of HR and RMSSD reactivity can be found in Supplementary Table 2. The final model included a group × phase interaction (and respective main effects), a group × PSS interaction (and respective main effects) and antidepressant medication. The PSS × group interaction (F = 10.41, p = 0.002) and the group × phase interaction (F = 2.99, p = 0.019), and antidepressant medication (F = 5.78, p = 0.020) were significant. Patients with MDD had overall higher HR throughout the experiment (Figure 3A). Post hoc contrasts revealed that in controls, HR increased significantly during the first (t = −5.11, p.adj < 0.001) and second stressor compared to baseline (t = −4.01, p.adj = 0.001) and returned to approximate baseline levels (or slightly lower levels) during first (t = 0.89, p.adj = 0.421) and second recovery phase (t = 2.02, p.adj = 0.080) compared to baseline (see Table 2 for exact statistics).

In patients with MDD, HR significantly increased during the first stressor (t = −2.90, p.adj = 0.011 compared to baseline) and returned to baseline levels during recovery (t = 0.43, p.adj = 0.667). However, no significant reaction to the second stressor was observed, neither compared to baseline, nor the first recovery period (t = −0.67, p.adj = 0.534 and t = −1.09, p.adj = 0.346, respectively). During the second recovery phase, a non-significant reduction of HR occurred compared to the second stressor (t = 2.16, p.adj = 0.062) (Figure 3A). Exploration of the significant group × PSS interaction showed an inverse U shape at baseline, where high stress levels were associated with higher HR in healthy controls, but with lower HR in patients with depression (see Supplementary Figure 1). For an overview of planned contrasts per phase for patients and controls, see Table 2. Of note, in an additional exploratory model within the MDD group, the main effect of PSS (F = 19.16, p = 0.002), but not the phase × PSS interaction (F = 1.22, p = 0.306) was significant, indicating that chronic stress influenced overall HR levels but did not (additionally) alter the reactivity in patients with MDD. No effect of PSS (F = 1.26, p = 0.271) or the phase × PSS interaction (F = 0.370, p = 0.830) occurred for healthy controls. It should however be noted that healthy controls had in general much lower stress levels than patients with MDD and an exact disentanglement of group and stress effects was therefore not possible.

For RMSSD reactivity, group (F = 13.98, p < 0.001), phase (F = 6.36, p < 0.001) and a group × PSS interaction (F = 4.26, p = 0.043) were significant. Antidepressant medication use was not a significant predictor (F = 0.35, p = 0.557), but since RMSSD is known to be influenced by this measure, was kept in the model. Post hoc contrasts by group revealed that controls showed significantly lower RMSSD levels during the second stress task when compared to the first recovery phase (t = 3.18, p.adj = 0.010) and significant recovery after the second stress exposure (t = −4.81, p < 0.001). Reactivity of the first stressor compared to baseline was not significant (t = −1.92, p.adj = 0.168). Patients with depression showed no significant reactivity at any time point (all p.adj > 0.170). The data is displayed in Figure 3B and Supplementary Table 2.

To assess if, and which phase of HR reactivity could best be used to distinguish depressed patients from controls, we performed binary logistic regression with LOOCV. The best distinction of patients and controls was achieved using HR with RMSSD during recovery. Supplementary Table 3 shows that, as a unique predictor, HR outperformed HRV indexes for classification. During recovery, combining these indexes led to a slightly better classification than for HR alone. Predictive accuracy was lowest during stress (65.45 and 63.63%) and highest during recovery from stress (first recovery 75.47% and second recovery 72.22%, respectively). Receiver operating characteristic (ROC) curves comparing various conditions can be found in Figure 4.

As expected from the literature (64–66), our data clearly show overall higher subjective stress ratings in patients with MDD compared with healthy controls, both over the past month (as measured by the PSS) and during all phases of the experimental manipulation (as measured with the VAS). We here also replicate the finding that patients with depression show higher basal levels of HR and lower levels of HRV than controls (29), thereby confirming the representability of our sample. Most intriguing, however, is the finding that patients could not subjectively recover from stress and showed blunted reactivity to the second stressor.

Physiological habituation, conceived as significantly attenuated responses to a second stress exposure, can be influenced by many factors. In their review, Hughes et al. (67) conclude that healthy controls often show habituation of HR to stressors comparable to the here used MAT, but that this process can be disrupted by the presence and observation of the participants by the experimenter (60, 61, 68, 69). Since we designed the task to contain psychosocial stress elements including observation by the experimenter, it is not surprising that healthy controls do not show (significant) adaptation to the stressor. In our opinion, the psychosocial stress elements are crucial to test relevant stressors since psychosocial stress is highly important for depression pathogenesis. Opposed to subjective stress reactivity, HR levels in patients with depression did fully recover from the first stress exposure but did not show any reactivity to the second stressor. This may suggest a rapid habituation in MDD patients, despite the above-mentioned psychosocial stress elements, or alternatively may point to a dysfunctional or exhausted stress response. The first scenario would suggest that the observed physiological non-responsivity is an adaptive strategy to reduce energy expenditure. From this perspective, while essential to mount sufficient stress responses in the face of acute danger, a fast adaptation to a similar stressor would be very beneficial for energy resources. Yet, adaptive habituation would fail to account for the chronic high stress levels, increased HR and consistently blunted HR variability that have been observed in MDD (18, 70) and is therefore an unlikely hypothesis.

The physiological response could thus perhaps find a more integrative explanation when conceived as the consequence of allostatic load, as introduced by McEwen and Stellar (19, 21, 22, 71, 72). The authors propose that an initially adaptive mechanism to mobilize energy resources (allostasis) can lead to considerable wear and tear if it persists long term, as is the case with chronic stress (allostatic load). In line with these observations is a recent study: In a large population of healthy young adults, high levels of stress also led to blunted HR responses and additionally led to lower habituation over time (24), although the authors did not account for depressive symptomatology. In our sample, the group of healthy people was too small to investigate this differential effect in detail. However, we did observe a qualitatively different effect of chronic stress on baseline HR in healthy controls compared to patients with depression: While chronic stress increased HR in the controls, it decreased HR in patients with depression. It is therefore essential to account for this relationship in studies involving patients with MDD and would be highly informative to assess the effect of depressive symptoms in a large healthy working population.

The here observed discrepancy of subjective and physiological responses to stress in depression has also been found by others. In their study, Söder et al. found a higher discrepancy of recorded HR and reported subjective stress in patients with depression compared to controls during a noise stressor (73). This discrepancy of higher self-reported stress and relatively lower HR reactivity may contribute to the maintenance of chronic stress levels and deserves further investigation.

The repeated stress paradigm as tested in this study, albeit an artificial exposure, still seems highly relevant for the daily experience of patients with clinical depression, who often report a continuous feeling of high subjective stress and suffer from repeated minor stressors with which they seem unable to cope adequately. This inability to cope with repeated stressors might be somehow related to the decoupling of the biological stress response, caused by exhaustion. A next step to test this hypothesis would be to test patients with depression in daily life, carefully recording stressors, subjective responses, and physiological parameters using wearable devices.

In response to the first stressor, patients with MDD showed increased subjective stress levels of approximately equal magnitude as controls, and a corresponding increase in HR. However, RMSSD showed nearly no variation throughout all phases of the task within the MDD population. While healthy controls also did not show vagal withdrawal to the first stressor, this observation could possibly be attributed to anticipatory stress levels at baseline. This is supported by the apparent recovery after the first stressor and significant reactivity to and recovery from the second stressor.

Ample research has documented blunted HR and HRV responses to stress in both participants with depressive symptoms (70) and patients with clinical depression (18). While we did not observe altered HR reactivity in patients with depression in response to the first stressor, the response to the second stress exposure was blunted, although the blunting could perhaps better be described as exhausted stress response (see section “Blunted Heart Rate Responses in Major Depressive Disorder: Adaptive Habituation or System Exhaustion?”). On the other hand, the association between depression and blunted cardiac vagal reactivity was confirmed and may point to consistently exhausted cardiac vagal modulation.

Blunted reactivity to stress in depression is not a novel concept in stress research: several reviews have investigated the matter with a focus on HPA axis reactivity. The most recent meta-analysis reveals that active, self-relevant stressors provoke similar, though sex-specific changes in cortisol release: whereas female participants with MDD show blunted reactivity to stress, male participants with MDD show increased responses (74). This corroborates findings with regard to the here described blunted cardiac vagal modulation in response to stress in our sample of female participants with MDD and highlights the importance of a careful study design accounting for relevant variables such as sex. The biological background behind this blunted reactivity alongside elevated HR to the first stressor may point toward primarily sympathetic activation rather than parasympathetic withdrawal during stress. This could be highly relevant for current therapies such as vagus nerve stimulation for depression, since restoring the balance between parasympathetic and sympathetic nervous systems should evaluate both sympathetic and parasympathetic output to quantify efficacy. However, as nerve activity was not assessed directly in this study, it is impossible to draw conclusions on this hypothesis.

Despite the lack of (complete) subjective recovery after the first stress exposure, the patient group did show significant physiological recovery. This is interesting, as these results may inform about a possible process contributing to the maintenance of subjective stress levels: a failure to subjectively recover from minor stressors, as observed here, may lead to continuous maintenance of perceived stress and ultimately, increase the impact estimation of a stressor (66). This altered stress perception may also be the cause for the change in physiological parameters at baseline: e.g., chronic stress may increase basal HR levels and reduce parasympathetic activity.

The “why’s” and “how’s” of these physiological changes and their link to subjective stress can find a possible explanation in the literature of “perseverative cognition”: Brosschot et al. proposed that particularly rumination (both retrospective and prospective), as frequently encountered in depression, can negatively impact the consequences of psychological stress (75). As such, repetitive or “perseverative” thoughts about stress, and not necessarily the stressor itself, would lead to persistent physiological changes such as a chronically increased HR, blood pressure and altered cortisol levels (69). While the present study was not designed to directly test hypotheses related to perseverative cognition, it is conceivable that it impacts these physiological parameters. Particularly patients with depression did have high rumination scores. Yet, in this study, perceived stress rather than rumination scores increased model performance for the group with depression considerably, while rumination and perceived stress did not improve model fit for controls (data not shown). Therefore, preservative cognition may not be the most relevant factor for physiological changes in depression when accounting for chronic stress levels.

In the present study, the most accurate classification of patients could be achieved based on HR, with slightly improved classification when combined with RMSSD. Differences were largest during recovery from stress and lower accuracy was achieved during stress. Others have recently published results on classifying individuals with depression based on cardiovascular parameters, with similar accuracy as here reported. To illustrate, Kuang et al. found high classification accuracy (86.4%) with RMSSD throughout a series of physiological ANS stimuli (36). In a stress task, Byun et al. also found the clearest distinction with a HRV feature (SDNN) during stress recovery (34). Of note, although HR or HF-HRV were not the most important features reported by Byun et al., we here specifically focused on these parameters, given their biological relevance and relative robustness regarding interpretability.

The usefulness of an objective biomarker is obvious in a time where remote monitoring for the clinic becomes tangible. While this additional data does not replace diagnoses or contact with a clinician it may help making informed decisions for possible treatment augmentation, in patients with partial response to therapy, or for secondary prevention in recovered patients. These findings could also be particularly interesting in the light of cardiovascular health. This is certainly important in light of the increased risk for cardiovascular disease observed in depression (76). While we here did not have longitudinal data to ascertain whether HR could be used as a state marker, others have shown that in remitted patients, cardiovascular parameters normalized toward the control group (77). More recently, Hartmann et al. showed that in their brief longitudinal study, indeed, HF-HRV (significantly) and RMSSD (non-significantly) increased in patients with depression in the course of a 2-week antidepressant treatment regimen, emphasizing the need for further longitudinal studies to assess the feasibility of cardiovascular measurements as state markers for depression (78).

Firstly, the generalizability of this study is limited to female participants. Due to the well-documented sex differences in autonomic activity and sample sizes necessary to control for this effect, this was a necessary choice. However, menstrual status is known to have important effects on autonomic (re)activity. While recorded in this study, and an approximately equal proportion of patients and controls in each menstrual phase, analyses were underpowered to include this factor as covariate. Larger studies should include menstrual cycle as covariate and explore whether our findings are also applicable to men. It should also be mentioned that subjective stress experience and physiological reactivity were not always congruent. This well-known apparent incongruence (79) can arise from confounding factors such as age (80), sex (81) or the nature of the stressor. Other, external factors in larger samples may show different results. For instance, it has been shown that early life adversity can have profound effects on (physiological) stress reactivity (82), and future studies should account for this variable. In addition, habituation might better be tested in a stress paradigm without social stressors. However, this was not yet clear when this study was conceived in 2015, and social stress may actually be crucial to test ecologically relevant habituation (i.e., social compared to non-social stressors). A frequently reported measure for HRV analyses is LF-HRV. LF-HRV has previously been interpreted as reflecting cardiac sympathetic innervation. Yet, in resting conditions (as is the case here for parts of our experiments, the LF band reflects baroreflex activity and not cardiac sympathetic innervation (83). Interpretation of LF-HRV thus is complex and is not suitable to investigate our study aims. However, for the sake of completeness, it is reported in Supplementary Section. Finally, this is a cross-sectional study and therefore by definition, our results are correlational. Large-scale, longitudinal studies are needed to further investigate the potential of HR and HRV during stress as state markers for depression and its usefulness for the clinic.

We here show that stress reactivity is functional in patients with depression but appears exhausted in face of recurrent stress. The specific design of this study enabled us, for the first time, to experimentally demonstrate physiological (HR/HRV) exhaustion of the stress response in depression, lending support to Mc Ewen’s conception of “allostatic load” as destructive consequence of chronic stress, and in particular its importance for depression. Importantly, our results also show that re-exposure to a stressor of the same nature can be used to study exhaustion. Future studies should therefore replicate these findings and extend them to male subjects.