We collected blood and stool samples from all patients who were admitted to Youjiang Medical University for Nationalities between August 2019 and November 2020. Patients with hemiplegia, speech disorder, hemianesthesia, distortion of commissure and stroke history, with CT/MRI evidence (while the other neuroinflammation diseases such as encephalitis, tumor etc. Were excluded) were diagnosed and recruited into stroke group; stroke patients with depression, frustration and/or interest loss persistent over 1 week were diagnosed and recruited into PSD group. The HDRS score was employed to evaluate the severity of depression, to minimize the system error, only those patients scored between 7 and 17 (mild depression) was recruited. Patients ranged in age from 44 to 73 years (no significant difference was found in the HDRS score between PSD patents <59 years old and ≥59 years old, data not shown). This study was approved by the ethics committee of Youjiang Medical University for Nationalities. Written informed consent was obtained from the patients to allow for their anonymized data to be collected and analyzed for research purposes. Prior to conducting further genetic analyses, patients were stratified into three groups that included healthy controls (n = 30, group A), stroke patients (n = 34, group B), and PSD patients (n = 26, group C).

The UPARSE-OTU and UPARSE-OTUref algorithms of the UPARSE package were used for the sequence analysis. The alpha (within-sample) and beta (between-sample) diversities were analyzed using internal Perl scripts in which sequences with ≥97% similarity were assigned to the same OTU. We selected a representative sequence for each OTU and annotated the taxonomic information for each representative sequence using the RDP classifier. To calculate alpha diversity, we sparsed the OTU table and computed three metrics that included Chao1 that estimates species abundance, Observed Species that estimates the number of unique OTUs for each sample, and the Shannon index. The three metrics provided the basis for the rarity curves. For quality control, we observed the sequencing results with cluster density of 76 ± 8.5 K/mm2, and only those >85% of the clusters were used for down-streaming analysis, besides, a quality score of Q > 31 was assigned to 94.5% of all bases from both reads. Only the reads passed filtering were used for down-streaming analysis.

The relative abundance of bacterial diversity was displayed using a Krona plot. Prior to cluster analysis, principal component analysis (PCA) was utilized to reduce the dimensionality of the original variables using the QIIME package. We used unweighted UniFrac distances for principal coordinate analysis (PCoA) and unweighted arithmetic mean pair group method (UPGMA) clustering. The principal coordinates from complex multidimensional data were visualized using PCoA. UPGMA clustering is a hierarchical clustering method that can be used to interpret the distance matrix.

The 16S rDNA analysis was performed to identify and quantify microbes in the blood samples collected from the patients (Supplementary Table S1). We used the Statistical Analysis of Metagenomic Profiles software to examine differences across the three groups in regard to the abundances of individual microbial taxa found in the samples. A quantitative analysis of biomarkers within each group was conducted based on linear discriminant analysis effect size (LEfSe). This method was designed to analyze datasets where the number of species is much higher than is the number of samples, and this method can be used to determine the features that are most likely to explain the differences between biological classes by testing for statistical significance and biological consistency and also by estimating the effect sizes of predicted biomarkers.

The ‘limma’ package (Yu et al., 2019) in R was used to analyze DEGs between the case (atheroma plaque) samples and control samples. Statistical significance was set at p < 0.05. The ‘clusterProfiler’ package (Yu et al., 2012) in R was used to functionally analyze the enriched pathways in which the DEGs were involved based on Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Statistical significance was set at p < 0.05.

Gene set enrichment analysis (GSEA) was performed using GSEA software (Subramanian et al., 2005) with reference to the gene sets c2. cp.kegg.v6.2., symbols. gmt, and c5. bp.v7.0. entrez.gmt, all of which were downloaded from the Molecular Signature Database (http://software.broadinstitute.org/gsea/index.jsp) (Liberzon et al., 2015). Statistical significance was defined as a nominal p-value < 0.05 and a false discovery ratio <0.25.

For the control group, we identified 1,492 DEGs in stroke patients and 4,376 DEGs in patients with depression (Figure 1A). Messenger RNA (mRNA) expression pattern analysis revealed that 29 DEGs were upregulated and 4 DEGs were downregulated in both stroke and PSD patients (Figure 1B). These genes were used to define the molecular signature for PSD. In contrast, the DEGs that were upregulated in stroke patients but downregulated in PSD patients (or vice versa) were defined as the molecular signature of recovery from PSD. We observed that the expression pattern of these two types of DEGs exhibited the ability to distinguish stroke or PSD patients from control samples (Figures 1C,D). DEGs that were expressed only in PSD patients (not in stroke patients) were defined as depression-specific genes, and their expression pattern could be used to distinguish PSD patients from healthy controls (Figure 1E).

We identified various gene expression profiles that may contribute to depression. Genes exhibiting |log FC|_st> 0 and |log2FC|_De > 0 maintained high expression in stroke patients but low expression during recovery (Figure 1Fa). Genes with |logFC|_St < 0 and |logFC|_De < 0 maintained low expression during stroke and recovery (Figure 1Fb). Genes |logFC|_St > 0 and |logFC|_De < 0 “over-recovered” during stroke recovery based on the observation that their levels during recovery exceeded their levels in controls (Figure 1Fc). Genes exhibiting |logFC|_St < 0 and |logFC|_De > 0 also “over-recovered” during stroke recovery based on the observation that their levels during recovery fell below the levels in controls (Figure 1Fd).

To explore the driving factors that may induce depression after stroke, enrichment analysis was performed on three types of DEGs, including those common to stroke and depression, those whose expression differed substantially between stroke and recovery, and those whose expression was altered only in depression. The GO enrichment results revealed that these three types of DEGs were involved in biological processes (BPs) associated with metabolism, immunity, microbes, and nerves (Figure 2A). Depression-specific genes play a significant role in the activation of T cells and xenogenesis.

GSEA indicated that the BPs associated with the genes common to stroke and depression and also depression-specific genes were significantly associated with immunity and neuronal differentiation. Additionally, BPs associated with depression-specific genes were associated with neuroinflammation, immune inflammation, and apoptosis (Figure 2B). KEGG enrichment analysis revealed that these three types of DEGs were involved in pathways related to intestinal flora, metabolism, brain disease, and immunity (Figure 2C). GSEA also demonstrated that the KEGG pathways related to metabolism, brain disease, immunity, and intestinal flora were enriched in patients with stroke and depression (Figure 2D).

We observed overlaps in the intensity and retention time of each chromatographic peak in the quality control (QC) samples (Supplementary Figures S1A,B). The correlation coefficients between the QC samples were >0.9, thus indicating good repeatability (Supplementary Figure S1C). More than 80% of the peaks in the QC sample possessed a relative standard deviation of ≤30%, thus indicating the high stability of the instrument analysis system (Supplementary Figure S1D). Based on Hotelling’s T2 test, all QC samples were within the 99% confidence interval, thus indicating that the experiment was reproducible (Supplementary Figure S1E). Fluctuations in the QC samples were within ±3 standard deviations and were reflective of the normal fluctuation of the instrument (Supplementary Figure S1F).

Among the metabolites with chemical classifications, organic acids and derivatives were found in large concentrations, while organosulfur compounds were found in small concentrations (Figure 3A). Compared to the controls, we identified 3,443 dysregulated metabolites in stroke patients (p < 0.05) (Figure 3B). We also determined that these dysregulated metabolites could be used to distinguish between stroke patients and controls (Figure 3C).

Compared to the controls, there were several metabolic pathways that were significantly enriched by differential metabolites in stroke patients, including ascorbate and aldarate metabolism and also valine, leucine, and isoleucine biosynthesis (Figure 3D). Metabolic pathway enrichment analysis was performed based on DEGs and differential metabolites (Figure 3E). Furthermore, we observed extensive interactions between the differential metabolites and DEGs (Figure 3F). The KEGG global metabolic network analysis also revealed that the purine metabolism pathway possessed the most significant enrichment score (Figure 3G) and that DEGs and differential metabolites play an important role in purine metabolism (Figure 3H).

To identify the metabolic abnormalities in patients with PSD, differential expression analysis was performed using non-targeted metabolomics sequencing data. Patients with PSD exhibited 5,408 differential metabolites when compared to control samples and 1,749 differential metabolites when compared to stroke patients (Supplementary Figure S2A). There were 53 abnormal metabolites common to PSD and stroke patients, and these metabolites could represent the metabolic basis of stroke as a precondition for depression (Figure 4A).

Among the common abnormal metabolites, we observed a decrease in 21 differential metabolites and an increase in seven differential metabolites during the development of stroke into PSD (Figure 4B). We also examined the abundance of these metabolites in all three patient groups (Supplementary Figure S2B and Supplementary Figure S4C). For example, MI80T35_Acamprosate decreased in stroke patients but increased significantly during the process of excessive recovery during the development of PSD (Figure 4D). Furthermore, we determined that these different metabolites in PSD patients can help distinguish between PSD and stroke patients (Figure 4E). As these metabolites represent the occurrence of metabolic changes complicated by depression during stroke recovery, further analysis of specific metabolites observed in PSD patients may help to distinguish PSD patients from controls and stroke patients (Supplementary Figure S2C).

Metabolic pathway enrichment analysis was performed using common genes, genes linked to excessive recovery, and specific metabolites. Valine, leucine, and isoleucine biosynthesis pathways were significantly enriched (Figure 4F). A detailed examination of the pathways associated with DEGs, the KEGG/metabolic-related pathways, and the three types of metabolites are listed above (Figure 3B, Supplementary Figure S2D). Additionally, we mapped significantly enriched disease-related metabolic pathways (Supplementary Figures S2E–H, Supplementary Figure S3 and Supplementary Figure S4).

We observed extensive interactions between the DEGs and metabolites (Figure 4G). Our analysis of the influence of differential metabolites and DEGs on the KEGG global metabolite network revealed that four metabolic pathways were significantly enriched, including the “metabolic pathway,” inositol phosphate metabolism, purine metabolism, and glycolysis/gluconeogenesis (Figure 4H). We were able to identify specific metabolites and DEGs that were extensively involved in long-term depression, the MAPK signaling pathway, the PI3K-Akt signaling pathway, and pathogenicity/E. coli infection in patients with PSD (Supplementary Figure S5). Our observations of differential metabolites in patients with PSD may reflect metabolic abnormalities which could be influenced by the intestinal flora.

We performed α-diversity analysis to measure the diversity of the microbial communities within the samples. Reasonable amounts of sequencing data were acquired from the three groups, and species abundance was relatively high (Supplementary Figure S6). Therefore, we examined the biodiversity index of the sequencing volume of each sample at different sequencing depths, and we observed that the sequencing data volume of the three groups was large enough to reflect most microbial information in the samples (Supplementary Figure S6B). Our analysis also revealed high species abundance and species uniformity across the three groups (Supplementary Figure S6C). Additionally, we determined that this study possessed an appropriate sample size and high species richness (Supplementary Figure S6D).

In regard to the species composition of gut microbes in the samples, we observed that the abundances of each group were different at each classification level (phylum, class, order, family, genus, and species; Figure 5A). The abundance of each sample was different at each taxonomic level (phylum, class, order, family, genus, and species; Figure 5B). Additionally, we observed that the abundance and diversity of intestinal flora were highest in controls, and this was followed by stroke patients and then PSD patients (Figure 5C).

To test for significant differences in the microbial communities among samples (groups), a β-diversity analysis was performed. Species diversity was highest in the samples collected from stroke patients, lower in controls, and lowest in PSD patients (Figure 5D). Based on the linear discriminant analysis effect size (LDASe) analysis, we determined that communities differing significantly between controls and stroke patients could be classified into six categories, while those differing between controls and PSD patients could be classified into five categories (Figure 5E). Based on the linear discriminant analysis, we identified the microbial population that played an important role in each group (Figure 5F). Species with an LDA value >2 were considered to be statistically significant biomarkers among the groups (Figure 5G).

Additionally, STAMP difference analysis was performed to identify species that differed significantly in terms of abundance among the groups. We observed that the abundance of g_Erysipelotrichaceae UCG-003 in the control samples was significantly different from that in the stroke patients, while g_pseudomonas abundance was significantly different between controls and patients with PSD (Figure 5H).

When we compared the intestinal microbes between stroke patients and PSD patients using the LEfSe analysis, we observed that the communities exhibiting significant differences could be classified into two categories (Figure 6A). Based on LDA, we identified microbial groups that played an important role in each group of patients (Supplementary Figure S7A). Additionally, we identified biomarkers that were able to differentiate between the control and stroke patients or between the stroke and PSD patients (Supplementary Figure S7B). Based on the STAMP difference analysis, we determined that g_Erysipelotrichaceae UCG-014 abundance was significantly different between stroke and PSD patients (Supplementary Figure S7C).

Compared to the control group, we observed 65 common dysregulated intestinal microbes in patients with stroke and PSD (Figure 6B). During the development of stroke into PSD, 22 common dysregulated intestinal microbes underwent an initial upregulation (Figure 6Ca) that was followed by downregulation (Figure 6Cb). Of these, 12 common dysregulated microbes recovered after being upregulated and then returned to normal levels (Figure 6Cc). In contrast, nine other common dysregulated microbes recovered after downregulation but could not return to normal levels (Figure 6Cd). Furthermore, we identified a total of 474 depression-specific disordered microbes in PSD patients, and these microbes were able to distinguish among controls, stroke, and PSD patients (Supplementary Figure S7D). Additionally, a correlation analysis was performed to test the relationship between differential metabolites and differential microbes. There were 3,038 significant metabolite-microbe interactions, including 59 differential metabolites and 518 differential operational taxonomic units (Figure 6D). We also constructed a global regulatory network to explore the changes in microbes that affect phenotypic function (Figure 6E).