Stressful events induce long-term gut microbiota dysbiosis and associated post-traumatic stress symptoms in healthcare workers fighting against COVID-19.

OBJECTIVE: The microbiota-gut-brain axis is a key pathway perturbed by prolonged stressors to produce brain and behavioral disorders. Frontline healthcare workers (FHWs) fighting against COVID-19 typically experience stressful event sequences and manifest some mental symptoms; however, the role of gut microbiota in such stress-induced mental problems remains unclear. We investigated the association between the psychological stress of FHW and gut microbiota.
METHODS: We used full-length 16S rRNA gene sequencing to characterize the longitudinal changes in gut microbiota and investigated the impact of microbial changes on FHWs' mental status.
RESULTS: Stressful events induced significant depression, anxiety, and stress in FHWs and disrupted the gut microbiome; gut dysbiosis persisted for at least half a year. Different microbes followed discrete trajectories during the half-year of follow-up. Microbes associated with mental health were mainly Faecalibacterium spp. and [Eubacterium] eligens group spp. with anti-inflammatory effects. Of note, the prediction model indicated that low abundance of [Eubacterium] hallii group uncultured bacterium and high abundance of Bacteroides eggerthii at Day 0 (immediately after the two-month frontline work) were significant determinants of the reappearance of post-traumatic stress symptoms in FHWs. LIMITATIONS: The lack of metabolomic evidence and animal experiments result in the unclear mechanism of gut dysbiosis-related stress symptoms.
CONCLUSION: The stressful event sequences of fighting against COVID-19 induce characteristic longitudinal changes in gut microbiota, which underlies dynamic mental state changes.

Introduction

The ongoing pandemic of coronavirus disease 19 (COVID-19) is still a global pandemic that has affected more than 200 countries. Healthcare workers, who relieve suffering and save lives on the frontline, are the most important force for solving this severe and urgent health crisis. Therefore, the World Health Organization highlights the fact that the protection of healthcare workers is a key step towards keeping patients safe. Besides the significantly higher proportion of COVID-19 infection in healthcare workers than in the general population, the mental health problems suffered by healthcare workers due to the pandemic have received increasing attention(Kang et al., 2020; Pappa et al., 2020; Wang et al., 2020, 2021; Yao and Xing, 2020).

The COVID-19 pandemic has exposed frontline healthcare workers (FHWs) to extraordinary levels of psychological stress. The mental burden placed on them includes constant fear of infection, facing social stigmatization, as well as feelings of depression, helplessness, exhaustion, frustration, guilt, and self-accusation, all of which are caused by overwhelming workloads and unceasing patient deaths(Chen et al., 2020; Kang et al., 2020). According to Elliot and Eisdorfer's stressors classification system (Segerstrom and Miller, 2004), stressful event sequences are the main type of stressor that FHWs encounter when treating COVID-19. This stressor is a focal event giving rise to a series of related challenges, but the affected individuals have a clear sense that, at some point in the future, these challenges will subside (Segerstrom and Miller, 2004). COVID-19 has markedly increased the prevalence of several types of psychological distress and mental health symptoms (Pappa et al., 2020; Wang et al., 2020, 2021; Yao and Xing, 2020). Compared with the general population, FHWs manifest more severe degrees of depression, anxiety, insomnia, and distress(Huang et al., 2020; Lai et al., 2020; Liu et al., 2020; Pappa et al., 2020; Song et al., 2020). Although several cross-sectional studies have characterized the mental health status of healthcare workers since the outbreak of COVID-19, a longitudinal study investigating the immediate and lasting effects of pandemic-related stressful event sequences on the mental health of FHWs has not been conducted.

The microbiome-gut-brain axis is one of the key pathways mediating adverse environmental stimuli-induced plastic changes in neural structure, function, and behavior (Burokas et al., 2017; Foster et al., 2017). Preclinical studies have demonstrated the extensive impact of stressful events on the composition and functional potential of host gut microbiota (De Palma et al., 2014; Montiel-Castro et al., 2013), such as reduced lactobacilli, beneficial lactic acid bacteria, and increased Clostridia in mice exposed to different types of stress (Bailey et al., 2011; Knowles et al., 2008). In contrast to numerous rodent studies, human studies regarding the reshaping of stressors on gut microbiota are limited and have only focused on brief naturalistic stressors such as academic examinations or negative events (Knowles et al., 2008; Michels et al., 2019). The alterations of gut microbiota after exposure to stressful event sequences in the fight against COVID-19, and their relationships with altered mental status in healthcare workers, are currently unclear. Gut microbiota dysbiosis is a common pathogenic basis underpinning various neuropsychiatric disorders, including autism, anxiety, depressive disorder, and schizophrenia (Zheng et al., 2016; Zhu et al., 2020a, 2020b). Manipulation of gut microbiota via probiotics, antibiotics, or germ-free feeding conditions significantly modulates stressful event-induced behavioral outcomes in rodents (Burokas et al., 2017; Lyte et al., 2020). In addition, several types of probiotics significantly improve stress-induced anxiety- and depressive-like behaviors in mice (Li et al., 2018; Stenman et al., 2020). In humans, probiotics also display some beneficial effects on mental health, for instance, altering emotional bias in healthy volunteers (Schmidt et al., 2015), and alleviating stress and anxiety in stressed adults (Chong et al., 2019). Therapeutic strategies through the "gut-brain" axis have been proven to be effective for future treatment practices in psychiatric medicine. Therefore, depicting the structure and function of the gut microbiota of healthcare workers who suffer from pandemic-related stressful event sequences is crucial for understanding the effects of stressful events on the human gut-brain axis and identifying therapeutic targets of gut microbes.

Here, we firstly aim to explore the alterations of gut microbiota in FHWs with significant stress-related symptoms. Secondly, we want to identify the disturbing microbes playing a crucial role in long-term post-traumatic stress in FHWs through longitudinal association analysis between microbiota and mental status.

Materials and methods

Study design and participants

This is a publicly-registered clinical study described in the ClinicalTrials (No.: NCT04443075) and approved by the Ethics Committee of the First Affiliated Hospital of Xi'an Jiaotong University (KYLLSL-2020–043). This study aimed to investigate the impact of stressful event sequences in the treatment of COVID-19 on the mental health and gut microbiota of FHWs. Gut microbiota and psychotic status were compared between 71 FHWs who had fought COVID-19 for two months in isolation wards in Wuhan, China, and 104 second-line healthcare workers (SHWs) who treated non-COVID-19-infected patients in routine hospitals. Moreover, dynamic changes in gut microbiota and psychotic status after exposure to the frontline medical work were delineated via a longitudinal investigation in FHWs at four time-points: immediately after they finished treatment and left the isolation wards (Day 0), after a two-week quarantine in a hotel (Day 14), four weeks after their return to normal life (Day 45), and a half year after the frontline work (Day 180). The process is shown schematically in Supplementary Fig. 1.

All participants were required (1) to be between 18 and 50 years old, (2) not to have taken antibiotics within 3 months before sample collection, (3) to have a body mass index (BMI) between 17.5 and 30. Participants were excluded if they fulfilled the following criteria: (1) had serious cardiovascular, hematologic, and/or endocrine disease; (2) had a history of cancer; (3) had active gastrointestinal diseases and/or serious systemic diseases; (4) had a history of brain organic diseases and/or developmental delay; (5) had psychiatric disorders such as mood disorder and anxiety disorder, (6) pregnant or lactating, (7) drank alcohol in the past week (liquor > 250 mL or beer > 1 bottle) or the previous day (liquor > 50 mL or beer > 50 mL). All participants signed the informed consent after knowing the study details.

Psychological stress symptoms assessment

Depression, anxiety, sleep, psychopathology symptoms, post-traumatic stress symptoms, and somatic symptoms of participants were evaluated using the 9-item Patient Health Questionnaire (PHQ-9), the 7-item Generalised Anxiety Disorder Scale (GAD-7), the Pittsburgh Sleep Quality Index (PSQI), the Symptom Checklist 90 (SCL-90), the Impact of Event Scale-Revised (IES-R), and the 15-item Patient Health Questionnaire (PHQ-15), respectively. Details of those psychological scales, lifestyle, and dietary habits assessment are described in the Supplementary Information.

Fecal samples

Fresh fecal samples were collected and immediately stored in commercial tubes with fecal DNA preservation agents (MGIEasy Stool Sample Collection Kit, BGI, Wuhan, China) and then stored at −80 °C until DNA was extracted. Bacterial genomic DNA was extracted using the E.Z.N.A.® Stool DNA Kit (Omega Bio-Tek, Norcross, GA, U.S.). The V1-V9 region of the bacteria 16S rRNA gene was amplified by PCR (95 °C for 2 min, followed by 27 cycles at 95 °C for the 30 s, 55 °C for 30 s, and 72 °C for 60 s and a final extension at 72 °C for 5 min) using primers 27F 5′-AGRGTTYGATYMTGGCTCAG-3′ and 1492R 5′-RGYTACCTTGTTACGACTT-3′, where a barcode is an eight-base sequence unique to each sample. Amplicons were extracted from 2% agarose gels and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, U.S.) according to the manufacturer's instructions.

Full-length 16S rRNA gene sequencing

Circular consensus sequence (CCS) libraries were prepared using SMRTbell™ Express Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, CA) and followed by immediate treatment with the Enzyme Clean Up Kit (PN: 101–843–100). The appropriate fractions for sequencing runs were identified on the Femto Pulse System (Agilent). After pooling the desired size fractions, the final libraries were further cleaned up and concentrated using AMPure PB beads (Pacific Biosciences PN:100–265–900). Finally, all libraries were checked for concentration using Qubit™ 1X dsDNA HS Assay Kit (Thermo Fisher PN: Q33231), and final size distribution was confirmed on the Femto Pulse. Purified SMRTbell libraries from the Zymo and HMP mock communities were sequenced on dedicated PacBio Sequel cells using the Sequencing Kit 2.0 chemistry. Purified SMRTbell libraries from the pooled and barcoded samples were sequenced on a single PacBio Sequel cell.

Acquisition of clear reads

PacBio raw reads were processed using the SMRT Link Analysis software version 9.0 to obtain demultiplexed CCS reads with the following settings: minimum number of passes = 3, minimum predicted accuracy = 0.99. Raw reads were processed through SMRT Portal to filter sequences for length (<800 or >2500 bp) and quality. Sequences were further filtered by removing barcode, primer sequences, chirmas, and sequences if they contained 10 consecutive identical bases.

Clear reads were clustered into amplicon sequence variants (ASVs) at 100% similarity (identical) using the Deblur denoising algorithm, which removes noise due to sequencing error (Amir et al., 2017). Then single-read clusters were considered as poor-quality ASVs and discarded in further analysis. Clusters of identical sequences allowed us to detect microbial changes at fine scale resolution. The phylogenetic affiliation of each 16S rRNA gene sequence was analyzed by the uclust algorithm within the usearch v11 software package (https://www.drive5.com/usearch/) against the silva (SSU132) 16S rRNA database (Edgar, 2018). The optimal identity threshold of ASVs for approximated species was > 98.6% in this study (Edgar, 2018; Johnson et al., 2019).

Statistical analyses

Statistical analyses were performed using R (v4.0.2). Demographic analyses were performed using the R package CBCgrps (v2.8.1). Alpha and beta diversity analyses were performed using the vegan package (v2.5–7). Bray-Curtis distance matrix based on the composition at the ASV level was calculated in vegan and subsequently used to perform principal coordinates analysis (PCoA). Permutational multivariate analysis of variance (PERMANOVA) was performed by adonis() function in the vegan package to compare Bray-Curtis distance between two groups. Fuzzy c-mean clustering was performed to depict different trajectories of gut microbiota using R package Mfuzz (v2.48.0). We calculated the within-cluster sum of squares for a range of cluster numbers, and the optimal number was chosen using the ‘elbow’ method. A random forest model with 10-fold cross-validation was performed (Feng et al., 2015) (R package randomForest, v4.6–14) to identify psychological scale-associated bacteria by the relative abundance of microbial species against four main psychological scales: IES-R, GAD-7, PHQ-9, PSQI.

Health planes were constructed to further analyze the dynamics of different microbial clusters, as previously described (Pan et al., 2019). In brief, principal coordinates analysis of the Bray-Curtis distance was performed at the genus level, and then samples of second-line healthcare workers were fitted to a two-dimensional plane using the least-squares method on the first three principal coordinates. The Euclidean distances from samples of frontline healthcare workers at each time point to these planes were calculated.

Linear mixed models including age, sex, and BMI as covariates were computed in R using the lm base function after log10-transformation and Z-score scaling of the data. Compositional Lotka-Volterra (cLV) analyses were performed in Python (v3.8.8) to infer the network interactions of stress-associated bacteria using the online scripts published before (Joseph et al., 2020). Bacteria associated with stress at Day 0 were identified as perturbations. P-values of estimated effects were computed by the bootstrap method. Directed networks were visualized using Cytoscape software (v3.8.0). P-values and false discovery rates (FDR) by Benjamini-Hochberg multiple testing corrections of < 0.05 were considered statistically significant.

Results

Stressful event sequences caused serious mental symptoms and disrupted the gut microbiome

First, psychological status and gut microbiota were compared between 71 FHWs and 104 SHWs. There were no significant differences in age, sex, and BMI between FHWs and SHWs (all P > 0.05, Supplementary Tables 1 and 2). At Day 0, FHWs presented significant post-traumatic stress symptoms, sleep disorder, anxiety, and depression (all P < 0.001, Fig. 1 A). In addition, psychopathology symptoms and somatic symptoms in FHWs were more severe than those in SHWs (all P < 0.001, Supplementary Fig. 2).

Fig. 1

Stressful event sequences cause serious mental symptoms and disrupt the gut microbiome in FHWs. (A) Total IES-R, PSQI, GAD-7, and PHQ-9 score in frontline healthcare workers (FHWs, n = 71) at Day 0 and second-line healthcare workers (SHWs, n = 104). (B) Chao1 index and the number of observed amplicon sequence variates (ASVs) in FHWs at Day 0 and SHWs. (C) PCoA based on the Bray-Curtis matrix. Boxes represent the 25th–75th percentile of the distribution; the thick line in the middle of the box indicates the median; whiskers extend to values with 1.5 times the difference between the 25th and 75th percentiles, and outliers are represented as dots. P values across multiple boxplots are calculated using the Mann-Whitney U test(A, B). IES-R, Impact of Event Scale-Revised; PSQI, Pittsburgh Sleep Quality Index; GAD-7, the 7-item Generalized Anxiety Disorder scale; PHQ-9, the 9-item Patient Health Questionnaire. PCoA, principal coordinates analysis.

We analyzed the gut microbiome from healthcare workers using full-length 16S rRNA sequencing. On average, 22,350 (standard deviation: 10,407) reads per sample were generated; 68.66% of these sequences are in the length of 1451–1500, suggesting that most of the reads are valuable (Supplementary Fig. 3A). Totally these reads were clustered into 40,809 ASVs, and we identified 6174 ASVs (with relative abundance higher than 0.001% and present in more than 10% of samples, Supplementary Table 3) for further analysis. The distributions of discarded and remaining ASVs were shown in Supplementary Fig. 3B. Alpha diversities measured using the Chao 1 index and observed ASVs were significantly lower in FHWs than in SHWs (all P < 0.05, Fig. 1B). Frontline work caused gut dysbiosis, characterized by compositional differences between FHWs and SHWs (P= 0.0001, Fig. 1C). Because SHWs were free to have meals they liked and wanted at home, while FHWs ate what was supplied by the local government at irregular times, there were differences in the regularity, balance, and structure of the diets between the two groups (all P < 0.01, Supplementary Table 2). Testing the influence of these variables on the gut microbiome, we found that BMI, diet regularity, and staple food structure affected the compositional difference, but stress exposure was indeed the strongest factor (R 2 = 0.024, P= 0.0002, Supplementary Fig. 4, Supplementary Table 4).

Gut dysbiosis persisted and different microbes followed discrete trajectories during the half-year follow-up

To confirm whether disturbed gut microbiome recovered during the follow-up, we compared the alpha and beta diversities in FHWs with those in SHWs and found that Chao 1 index and number of observed ASVs temporarily recovered at Day 14 but continued to decrease after Day 14 and reached the lowest level at Day 180 (Supplementary Fig. 5A). In addition, Simpson and Shannon index decreased significantly at Day 180 (all P < 0.0001, Supplementary Fig. 5B). Compositional differences in gut microbiota between FHWs and SHWs persisted during the follow-up (all P= 0.0001, Supplementary Fig. 5C).

To characterize the dynamics of the gut microbiome in FHWs, we performed fuzzy c-mean clustering using samples from 52 FHWs who completed 180 days of follow-up based on the mean relative abundance of microbial genera. This analysis identified six clusters of longitudinal trajectories of microbial fluctuation (Supplementary Fig. 6): (1) genera continuing to increase from Day 14 to Day 180; (2) genera with highest relative abundance at Day 14; (3) genera that increased from Day 45 to Day 180; (4) genera with the lowest relative abundance at Day 45; (5) genera with the highest relative abundance at Day 45; and (6) genera with the continued decline (Fig. 2 , Supplementary Table 5).

Fig. 2

Different microbes follow discrete trajectories during the half-year follow-up. The cluster of longitudinal trajectories used the mean relative abundance of bacteria at the genus level. For detailed data of microbes in each cluster, see Supplemental Table 5.

To further analyze the longitudinal changes of each cluster in FHWs, we compared the distance from each sample to its Day 0 status with that from each sample to the “healthy plane” (HP) based on each sample's PCoA coordinates of the Bray-Curtis distances within each cluster. All clusters, except for cluster 3, were closer to HP at each time point, indicating that their microbial profile was similar to the HP. Of note, cluster 3 at Day 14 exhibited a similar distance to Day 0 and the HP (Supplementary Fig. 7), suggesting that distinctive changes happened at Day 45. Moreover, both distances to Day 0 and the HP increasing with time indicated that each cluster of FHWs did not change toward Day 0 status or the HP but toward another disease status that could be referred to as post-traumatic stress because FHWs had post-traumatic stress at Day 180.

Longitudinal changes of psychological impact were correlated with gut microbiome changes in healthcare workers

To determine if the psychological impact in FHWs returned to baseline, we compared psychological scales in FHWs at each time point with those in SHWs. This analysis revealed that most stress-associated symptoms such as anxiety and depression recovered at Day 14 (all P > 0.05, Fig. 3 A,B), and did not reappear after Day 14. IES-R recovered at Day 14 (P > 0.05, Fig. 3C) and PSQI recovered at Day 45 (P > 0.05, Fig. 3D). It should be noted that both post-traumatic stress and sleep disorders in FHWs reappeared at Day 180 (both P < 0.05, Fig. 3C,D).

Fig. 3

Longitudinal changes of psychological impact are correlated with gut microbiome in healthcare workers. (A–D) Total GAD-7, PHQ-9, IES-R, and PSQI scores in frontline healthcare workers (FHWs) and second-line healthcare workers (SHWs) at Day 0, Day 14, Day 45, and Day 180, respectively. (E) Stress-associated species selected by the cross-validated random forest models. Heatmap shows each bacterium's normalised mean decrease accuracy (equivalent to importance) against psychological scales in the random forest models. C1-C6 means the cluster number of each genus to which the shown species belong. Boxplot illustration is provided in Fig. 1. P values across multiple boxplots are calculated using a one-tailed Mann-Whitney U test, with SHWs as the reference group (A). IES-R, Impact of Event Scale-Revised; PSQI, Pittsburgh Sleep Quality Index; GAD-7, the 7-item Generalized Anxiety Disorder scale; PHQ-9, the 9-item Patient Health Questionnaire.

To identify the species associated with stress symptoms, we used random forest regression of species against each psychological scale in the complete series at four time points. Five repeats of 10-fold cross-validation led to the optimal selection of 24 species for GAD-7, 33 for IES-R, 27 for PHQ-9, and 12 for PSQI, respectively (Supplementary Fig. 8). In total, 59 species were identified (Fig. 3E, Supplementary Table 6). Most species were from the genus [Eubacterium] eligens group, Bacteroides, Faecalibacterium, Lachnospiraceae NK4A136 group, and Streptococcus in Cluster 6 (genera with continuing decline). Other common species were from Bifidobacterium and Sutterella in Cluster 4 (genera with lowest relative abundance at Day 45), and the Lachnospiraceae ND3007 group in Cluster 3 (genera that increased from Day 45 to Day 180).

Stress-associated species induced the reappearance of post-traumatic stress in healthcare workers through bacteria interaction networks

Sleep disorder in FHWs at Day 180 was possibly caused by many factors, but post-traumatic stress symptoms were specifically induced by the fight against COVID-19 and were a key mental disturbance in FHWs. Therefore, to examine the role of the gut microbiome in the reappearance of post-traumatic stress symptoms, we explored the associations between species changes and psychological scale changes using two steps. First, we performed random forest regression by the relative abundance of species at Day 0 against changes in IES-R scores from Day 45 to Day 180. Four species were selected by the model; also, significant differences between FHWs at Day 0 and SHWs were noted (all corrected P < 0.01, Fig. 4 A, Supplementary Table 7). These four species were possible inducers of post-traumatic stress. Second, we used linear mixed models, with age, sex, and BMI as covariates, to analyze associations between microbial changes and changes in IES-R scores from Day 45 to Day 180 (delta association). Ten species with P < 0.05 were identified using the mixed models; these species were present in more than 20% of samples (Fig. 4B, Supplementary Table 8).

Fig. 4

Stress-associated species induced reappearance of post-traumatic stress in healthcare workers through bacteria interaction network. (A) Relative abundance (log10) of four species selected using the random forest model; significant differences between frontline healthcare workers (FHWs) at Day 0 and second-line healthcare workers (SHWs). (B) Delta associations between IES-R scores and ten species that were significant at P < 0.05 and also present in more than 20% of samples. *P < 0.05, ***P < 0.001. (C) Network of microbe-microbe interactions from Day 0 to Day 14, from Day 14 to Day 45, from Day 45 to Day 180, respectively. The three microbial interaction networks between four inducers (orange nodes) and ten species that caused the reappearance of post-traumatic stress symptoms (purple nodes) were reduced for visualization from the complete network in Supplementary Fig. 9. Red edges indicate facilitation (standardized effect > 0.5), and blue ones indicate inhibition (standardized effect < –0.5). The width of the line reflects the strength of the relationship. Solid lines indicate P values of the associations less than 0.05, and dash lines indicate more than 0.05, respectively. P values across boxplots are calculated using the Mann-Whitney U test with Benjamini-Hochberg correction (A). See also Supplemental Tables 7–9. IES-R, Impact of Event Scale-Revised.

We observed no significant associations between the four inducers selected using random forest regression and changes in IES-R scores (all P > 0.05, Supplementary Table 8). Therefore, we inferred that the four inducers might affect the reappearance of post-traumatic stress through microbe-microbe network interaction. To model this, we performed cLV analysis using four inducers as perturbation and ten species associated significantly with IES-R scores as disturbed microbes at time points from Day 0 to Day 14 (isolation period), Day 14 to Day 45 (recovery period), and Day 45 to Day 180 (reappearance period), respectively (Supplementary Table 9). The inferred network shows that Bacteroides eggerthii inhibited the growth of stress-associated species during the isolation period, although the P-value of these correlations computed by the bootstrap method was not significant (Fig. 4C, Supplementary Fig. 9). During the recovery period and the reappearance period, the [Eubacterium] hallii group uncultured bacterium had a positive effect, and Bacteroides eggerthii inhibited the growth of stress-associated species ([Eubacterium] xylanophilum group uncultured bacterium, a beneficial bacterium, was inhibited by Bacteroides eggerthii and Lactobacillus salivarius was facilitated by [Eubacterium] hallii group uncultured bacterium), suggesting that the [Eubacterium] hallii group uncultured bacterium had a beneficial effect on the bacterial interaction network during the reappearance of post-traumatic stress symptoms but Bacteroides eggerthii had the opposite effect.

Additionally, we identified two species, Dialister invisus DSM 15,470 and Acinetobacter sp. CIP 101,966, which had predictive value in the reappearance of sleep disorder in FHWs at Day 180 (Supplementary Table 10) and also correlated negatively with the changes of PSQI score (Supplementary Table 11) using the two-step method. Dialister invisus DSM 15,470 was depleted in FHWs at Day 0 (Supplementary Table 7), which suggests its beneficial effect in sleep disorders. Acinetobacter sp. CIP 101,966 were enriched in FHWs at Day 0, but it had no significant antagonistic effect with Dialister invisus DSM 15,470 in the interaction networks (Supplementary Fig. 10, Supplementary Table 12); it had few connections, suggesting that its role was not as important as Dialister invisus DSM 15,470.

Discussion

Our investigations in FHWs provide a unique opportunity to understand the impact of a special type of chronic stress exposure on the human gut-brain axis. As chronic stress is harmful to human health, randomized controlled human studies using stress factors as exposure factors are not allowed by medical ethics principles. Investigations on survivors and victims who have encountered disasters, such as earthquakes and massacres, are mainly used to base studies of the subsequent development of mental health problems and related human neurobiology of stress (Dyb et al., 2014; Lai et al., 2017). As all patients infected with COVID-19 in China were transferred to isolation wards, FHWs who were involved in their treatment over two months experienced typical stressful event sequences. The increased prevalence of mental symptoms in the healthcare workers in our study and those in several other studies (Huang et al., 2020; Pappa et al., 2020; Song et al., 2020; Yao and Xing, 2020), suggest that healthcare workers who treated COVID-19-infected patients faced vast amounts of stress. As the city where we recruited our participants was locked down before our study started, SHWs who experienced similar stress levels at work and in life before and after the outbreak of COVID-19 were included as controls. Although the two cohorts of healthcare workers were not randomly assigned, they had the same occupations and similar demographic features and socioeconomic status; thus, differences in gut microbiota can be attributed to exposure to stressful event sequences.

Based on reliable comparisons with SHWs, the gut microbiome of FHWs was characterized by lasting decreased alpha diversity. Low bacterial diversity occurs in various psychiatric diseases, including major depressive disorder and generalized anxiety disorder (Huang et al., 2018; Jiang et al., 2018). A previous study consistently detected a decreased alpha diversity in stress reflected by negative events (Michels et al., 2019). To note, our data revealed that low alpha diversity in FHWs persisted for half a year. Moreover, microbial community structure was long-term disturbed by stressful sequence events in FHWs based on PCoA analysis. Nevertheless, the symptoms of anxiety and depression in FHWs could be recovered at Day 14, suggesting that the symptoms of anxiety and depression may due to the absence of rest without interruption as reported before (Chen et al., 2020).

We observed six longitudinal trajectories of disturbed bacteria. Some beneficial bacteria such as [Eubacterium] hallii group in cluster 1 and Lachnospiraceae ND3007 group in cluster 3 gradually recovered over time; both of them belong to Lachnospiraceae, a kind of butyric-producing bacteria. In addition, Enterococcus faecalis in cluster 3 were reported to promote social activity and reduce corticosterone levels in mice following social stress in a recent study (Wu et al., 2021). Cluster 2, 4 and 5 only change briefly and would return to the baseline level over time. However, several bacteria in these clusters were also related to mental health. For example, [Ruminococcus] gnavus in cluster 2 was related to PSQI. Members of [Ruminococcus] gnavus group could degrade gastrointestinal mucins, which may lead to the breakdown of mucus and thereby increased permeability of the gut (Lyte et al., 1998). A previous study reported that stress decreases Muc2 synthesis in the goblet cells and subsequently reduces the goblet cell number via Notch signaling suppression, which impairs the intestinal barrier (Shigeshiro et al., 2012). Three species of Bifidobacterium in cluster 4 were associated with depression. Bifidobacterium spp. are the most commonly used probiotics (Gismondo et al., 1999), and probiotics containing Bifidobacterium strain may contribute to ameliorate chronic stress-induced abnormal brain plasticity (Ait-Belgnaoui et al., 2014). Previous study has shown that inoculation of Bifidobacterium infantis in germ-free mice could rapidly induce c-Fos activation in the paraventricular nucleus and reversed the exaggerated HPA stress response, suggesting that B.infantis probably inference brain via a neural-mediated pathway (Sudo et al., 2004). In addition, Klebsiella pneumoniae in cluster 5 was related to anxiety. Klebsiella pneumoniae is one of the oral pathobionts and could also promote colitis (Kitamoto et al., 2020).

Importantly, most of the stress-related bacteria were in cluster 6 with relative abundance continuously decreasing. This cluster mainly included Faecalibacterium, [Eubacterium] eligens group, and Bacteroides. [Eubacterium] eligens group and Faecalibacterium prausnitzii were related to anti-inflammatory effects in previous studies (Chung et al., 2017; Quevrain et al., 2016). Bacteroides spp. has been reported as important gamma-aminobutyric acid producers (Strandwitz et al., 2019). Host energy metabolism and immune functions critically depend on butyrate as a potent regulator. Experimental studies document butyrate-producing bacteria involved in stress-induced behavioral changes in rodents (Bangsgaard et al., 2012) enhancing the integrity of the gut barrier (Sherry et al., 2010), and oral supplementation with butyrate can suppress inflammatory responses and alleviate stress-induced brain-gut axis alterations in mice (van de Wouw et al., 2018). In our study, cluster 6 also contained several species from Lachnospiraceae and Roseburia that belonged to butyrate-producing bacteria and associated with depression and post-traumatic stress. Lachnospiraceae members could also produce short chain fatty acid that promote host serotonin biosynthesis in the gastrointestinal tract and impact gastrointestinal motility (Yano et al., 2015). These results suggest that gut microbiota may contribute to the intestinal permeability-increasing caused by stress, thus allowing bacteria to translocate across the intestinal mucosa and directly access both immune cells and neuronal cells of the enteric nervous system (Malan-Muller et al., 2018; Simeonova et al., 2018). The function of these health-associated bacteria indicates that gut microbiota may influence stress symptoms in FHWs through neurotransmitter, immune pathways, and microbial translocation, which is also supported by previous studies (Bailey et al., 2011; Jarbrink-Sehgal and Andreasson, 2020).

Intriguingly, we found the reappearance of post-traumatic stress symptoms in FHWs, which is consistent with a previous study that reported significantly higher post-traumatic stress from 13 to 26 months after the severe acute respiratory syndrome (SARS) outbreak in healthcare workers who treated SARS patients (Maunder et al., 2006). We identified four potential inducers involved in the flashback (the key symptom of post-traumatic stress). High abundance of Bacteroides eggerthii and low abundance of [Eubacterium] hallii group uncultured bacterium at Day 0 were essential determinants of flashback from Day 45 to Day 180; these were antagonistic in the inferred network interaction. We inferred that [Eubacterium] xylanophilum group uncultured bacterium were inhibited by Bacteroides eggerthii, which induced the appearance of flashback. By contrast, beneficial bacteria Lactobacillus salivarius facilitated by [Eubacterium] hallii group uncultured bacterium may help the recovery of flashback. The [Eubacterium] hallii group has been identified as the dominant producer of butyrate in the colon (Flint et al., 2015), and the observation of its decrease post-stress invites further investigation on its potential value during stress exposure. These microbial biomarkers for psychological flashbacks are potential therapeutic targets for post-traumatic stress symptoms.

Several limitations need to be noted. First, to disturb the medical staff as little as possible, we did not investigate the changes of bacterial metabolites in feces and blood such as short-chain fatty acids and the effects of the crucial bacteria in the stressed mice model, which are needed to help elucidate the mechanisms in the gut-brain communication as reported recently (Jarbrink-Sehgal and Andreasson, 2020). Second, not using negative controls or positive controls in our study may lead to neglecting the influence of the factors of DNA extraction kits, sampling method, contaminations, and sequencing methods (Qian et al., 2020). Despite these limitations, long-term alterations in gut microbiota following the stressful event sequences suggest their sustainable negative impact on the health of healthcare workers, and stress-associated bacteria identified in this study help to design prevention and treatment strategies for post-traumatic stress symptoms in frontline healthcare workers.

In conclusion, this study illustrates the stressful event sequences-induced characteristic profile of gut microbiota and dynamic changes in bacterial composition in humans, which provides a better understanding of the gut bacterial architecture of stress-induced mental disorders and highlights potential targets for future interventions.

Data availability statement

The 16S sequencing data are publicly available on Genome Sequence Archive (GSA) under the study accession numbers PRJCA006021.

Code availability statement

R scripts demonstrating how to reproduce all findings shown in the main figures are available at https://github.com/gaofengjie-2019/Gao2021

Funding source

This work was funded by the (No. 81771471); National Outstanding Youth Science Fund Project of National Natural Science Foundation of China (82022023); Basic Research Project of Natural Science Fund of Shaanxi Province (No. 2016JQ8026).

CRediT authorship contribution statement

Fengjie Gao: Data curation, Investigation, Formal analysis, Writing – review & editing. Ruijin Guo: Visualization, Formal analysis, Writing – review & editing. Qingyan Ma: Data curation, Investigation. Yening Li: Data curation, Investigation. Wei Wang: Visualization, Investigation. Yajuan Fan: Data curation, Writing – review & editing. Yanmei Ju: Formal analysis. Binbin Zhao: Data curation. Yuan Gao: Investigation, Writing – review & editing. Li Qian: Investigation. Zai Yang: Writing – original draft. Xiaoyan He: Writing – original draft. Xiaoying Jin: Writing – original draft. Yixin Liu: Writing – original draft. Yuan Peng: Writing – original draft. Ce Chen: Data curation, Investigation. Yunchun Chen: Investigation, Writing – review & editing. Chengge Gao: Investigation, Writing – review & editing. Feng Zhu: Visualization, Formal analysis, Writing – review & editing. Xiancang Ma: Visualization, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no competing interests.