We conducted a prospective cohort study in Spain as part of an ongoing longitudinal multi-national study [13] (https://mentalnet.cl/en/home/). We collected data through an online survey at two timepoints. Baseline assessments were performed from 24th April to 22nd June 2020 (during the initial pandemic outbreak in Spain). Follow-up assessments took place between 26th January and 8th March 2021, which in Spain coincided with the third pandemic wave and with administration of COVID-19 vaccines for the majority of HCWs.

The study participants were HCWs aged 18 years and older, recruited from different outpatient and inpatient healthcare facilities, with clinical and non-clinical duties, and not necessarily involved in the direct care of COVID-19 patients. Participation was voluntary. The snowball sampling strategy was as follows. During both the baseline and follow-up assessments, key stakeholders (e.g., hospital managers, heads of worker unions) from healthcare facilities located in the study regions (Andalusia, Madrid, and Murcia) forwarded the survey link to all HCWs. Participants were also asked to forward the survey to peers in order to enhance response rates. In addition, during the follow-up period, we sent email or telephone survey reminders to baseline participants. Participants’ personal and clinical data remained stored in separate datasets (e.g., pseudonymized), following European General Data Protection Regulation Standards. Baseline assessments are described elsewhere [14]. In this manuscript, we focus on participants who were assessed at follow-up only and on participants assessed at follow-up who had been assessed also during baseline procedures. We hereafter refer to these two subgroups as partial and full respondents, respectively.

All procedures contributing to this work comply with the Helsinki Declaration of 1975, as revised in 2013. It received approval from the Hospital La Paz Ethics Committee in Madrid (Madrid, Spain), study ID 4099, and was ratified by the local ethics committees from the participating sites.

Baseline assessments included the following COVID-19-related exposures: direct involvement in the care of COVID-19 patients (yes, no), adequate access to personal protective equipment, fear of getting infected, and fear of infecting loved ones (all rated from 0 to 3).

Both baseline and follow-up assessments included the following sociodemographic and mental health outcome variables.

Sociodemographic variables: Age in years, gender (male, female), educational level (primary, secondary, or university studies), and type of job. We collapsed job types into the following categories: physicians, nurses, health technicians (e.g., nurse, X-ray, and laboratory technicians), ancillary workers (e.g., security staff, drivers, administrative staff, and cleaning staff), other HCWs (e.g., clinical psychologists, physiotherapists, and biologists), and residential support workers (e.g., from mental health assisted living facilities, nursing homes).

Mental health outcomes: Psychological stress, as measured by the validated Spanish version of the 12-item General Health Questionnaire (GHQ-12) [15]; and probable depression symptoms, as measured by the validated Spanish version of the 9-item Patient Health Questionnaire (PHQ-9) [16]. We used widely accepted thresholds for detecting people screening positive for psychological distress (GHQ-12 higher than 2 points) [17, 18] and for depression (PHQ-9 score higher than 9 points) [19]. In addition, follow-up assessments also included posttraumatic stress disorder (PTSD) symptoms, as measured by the Spanish 5-item version of the Primary Care PTSD Screen for DSM-5 (PC-PTSD-5), where a total score higher than 2 points suggests probable PTSD [20]. The scale has been translated, but not validated in Spain.

Cronbach’s alphas were 0.87 (95 percent CI: 0.86, 0.88) for the GHQ-12 total score; 0.89 (95 percent CI: 0.88, 0.89) for the PHQ-9 total score; and 0.70 (95 percent CI: 0.68, 0.72) for the PC-PTSD-5. To control for region-level cumulative COVID-19 incidence, we calculated region-specific 14-day cumulative incidence rates 2, 4, 6, and 8 weeks after the start of the follow-up period and, as rates were stable over time, classified regions as “high” or “low” incidence depending on whether average cumulative incidence over time points fell under or over 750 cases per 100,000 based on visual examination of region-specific cumulative incidence rates (Supplementary Figures S1, S2).

First, we removed baseline respondents who provided informed consent but did not go on to initiate the survey (n = 95). We reported categorical variables as frequencies and valid percentages, and continuous and interval variables as either mean (standard deviation [SD]), or median (interquartile range [IQR]). Descriptive statistics were calculated separately for full and partial respondents.

Then, we used multivariable mixed-effects linear and binary logistic regression models to explore the associations between baseline variables, including sociodemographic characteristics (i.e., age, gender, and educational level) and COVID-19-related exposures (i.e., direct involvement in the care of COVID-19 patients, adequate access to personal protective equipment, fear of getting infected, and fear of infecting loved ones), and follow-up mental health outcomes (i.e., psychological distress, depression symptoms, and PTSD symptoms), defined both as continuous questionnaire scores and dichotomous variables. We conducted sensitivity analyses adjusted by baseline assessments of the follow-up outcome under consideration. We used baseline GHQ-12 score for the model where follow-up PTSD was the outcome, as we did not have estimates for the latter in the baseline assessment. The GHQ-12 score, an instrument that has good convergent validity with the PC-PTSD-5 and accurately detects PTSD in primary care settings [21]. To ease the interpretability of effect estimates from binary logistic regression models (odds ratios), we dichotomized the following baseline exposure variables: adequate access to personal protective equipment (adequate vs. inadequate), fear of getting infected (not or slightly afraid vs. considerably or extremely afraid), and fear of infecting loved ones (not or slightly afraid vs. considerably or extremely afraid).

Next, we used baseline and follow-up mental health outcomes to categorize respondents into three mental health trajectories, separately for psychological distress and for depression, according to whether they screened negative at baseline and follow-up (asymptomatic stable), positive at baseline and negative at follow-up (recovering), or positive or negative at baseline and positive at follow-up (persistently symptomatic/worsening). For instance, if a respondent screened negative in the GHQ-12 and positive in the PHQ-9 at baseline, and subsequently screened positive in the GHQ-12 and negative in the PHQ-9 at follow-up, they would belong to the persistently symptomatic/worsening trajectory for psychological distress and to the recovering trajectory for depression. We selected these trajectory categories because of their potential implications for clinical practice.

Finally, we explored the association between baseline exposures and longitudinal psychological distress and depression trajectory membership, using multinomial regression models where asymptomatic stable was considered the reference category. Baseline exposures were age group, gender, educational level, direct involvement in the care of COVID-19 patients, adequate access to personal protective equipment (adequate vs. inadequate), fear of getting infected (not or slightly afraid vs. considerably or extremely afraid), and fear of infecting loved ones (not or slightly afraid vs. considerably or extremely afraid).

All models were adjusted for confounding based on prior causal knowledge, using direct acyclic graphs and backdoor criteria [22]. Region-level cumulative incidence of COVID-19 was entered in all models as a fixed factor. Likert-type variables were included as continuous for main analyses—following an accepted analytical approach [23]. (In a set of sensitivity analyses, we included them as dummy variables to assess whether main results were robust to potential model misspecification.) We did not impute missing data. All analyses were performed using packages dplyr, gtsummary, flextable, ggplot2, psych, multinom of R Studio for Mac (Version 1.2.5042).

Of 1,807 respondents who answered the survey at follow-up (between 26th January and 25th March 2021), 1,471 (81.4%) completed the entire survey, with a median response time of 21 min. Most missing data pertained to the last section of the questionnaire, suggesting that data missingness was driven by survey extension and hence largely random. Respondents who did and did not complete the survey were comparable in terms of mean age (42 vs. 40 years, respectively) and gender distribution (78% vs. 74% female, respectively). Response rates were estimated across facilities and job types and ranged from 2.7% to 100% (see Supplementary Tables S1, S2).

There were 1,058 (59%) partial respondents (i.e., assessed only at follow-up) and 749 (41%) full respondents (i.e., assessed at both baseline and follow-up). Of note, this indicates that we retained 32% of the 2,370 original baseline respondents for follow-up assessments (see Figure 1). Sociodemographic characteristics of follow-up respondents, overall and divided into full and partial respondents, are shown in Table 1. In short, full respondents were more frequently female and more likely to have completed university studies than partial respondents. Also, while most full respondents were physicians or nurses, partial respondents included a larger proportion of residential support workers.

Overall, 56% of follow-up respondents screened positive for psychological distress, 21% for probable depression, and 51% for PTSD. Psychological distress, probable depression, and PTSD were more frequent among younger and female respondents, and respondents with higher educational levels—with substantial heterogeneity across specific job types (see Supplementary Table S3). Notably, follow-up mental health outcomes were comparable between full and partial respondents, with similar mean (SD) GHQ-12 score (3.9 [3.5] vs. 3.8 [3.4], respectively), mean (SD) PHQ-9 score (6.3 [5.1] vs. 6.4 [3.4], respectively), proportion of respondents screening positive for psychological distress (56% vs. 55%, respectively), and proportion of respondents screening positive for probable depression (21% vs. 21%, respectively).

Comparisons between baseline and follow-up mental health outcomes among full respondents are shown in Table 2. The proportion of respondents screening positive for psychological distress and probable depression decreased, respectively, from 74% to 56% and from 28% to 21%. Figure 2 shows the distribution of trajectories of depression symptoms and psychological distress over time, overall and across baseline covariates. Trajectories show that, in terms of depression symptoms, 66% respondents remained asymptomatic/stable, 15% recovered, and 19% remained symptomatic or worsened over time. In terms of psychological distress, 18% respondents remained asymptomatic/stable, 26% recovered, and 56% remained symptomatic or worsened over time. The distribution of trajectories was heterogeneous across baseline covariates.

Table 3 shows crude and adjusted estimates of the association between baseline exposures and follow-up mental health outcome scores among full respondents. Overall, women had higher scores (i.e., worse mental health) than men, and job-related factors such as direct involvement in the care of COVID-19, inadequate access to protective equipment, or fear of infecting oneself or loved ones were associated with higher negative mental health outcome scores–especially for PTSD, and with higher odds of testing positive for psychological distress, probable depression, and PTSD (see Supplementary Table S4). We repeated all models, first including further adjustment by baseline mental health outcome scores (see Supplementary Tables S5, S6) and then including Likert-type variables as dummy variables: results did not change.

Table 4 shows crude and adjusted associations between baseline exposures and trajectory membership for psychological distress and probable depression. Women showed higher symptom variability over time than men, as indicated by women’s higher odds of belonging to both the recovering and the persistently symptomatic/worsening categories for both psychological distress and probable depression. In terms of psychological distress, reporting inadequate access to personal protective equipment was associated with persistently symptomatic/worsening category membership. In terms of probable depression, fear of infecting loved ones was associated with recovering category membership. Results did not change in models including Likert-type variables as dummy variables.