A comprehensive search was conducted on the Sequence Read Archive (SRA) database¹ and European Nucleotide Archive (ENA) database² using the following keywords: “Parkinson’s disease,” “blood,” “brain,” “gut,” “feces,” “16S,” “rRNA,” “RNA-Seq,” “transcriptome,” “microbiota,” and “microbiome.” Gene expression microarray data of the brain in PD patients were retrieved in the Gene Expression Omnibus (GEO) database. Only samples from Homo sapiens were included in this study. All metadata sets were downloaded from the SRA database or GEO database according to the unique Project ID of each study. The geographic distribution of all included studies (Figure 1A) and workflow (Figure 1B) for this study was presented in Figure 1.

The 16S rRNA and RNA-Seq data downloaded from the ENA were processed with QIIME2 2020.8.0 (Bolyen et al., 2019) and Kraken 2.1.1 (Wood et al., 2019) according to the standardized pipeline, respectively. Briefly, the quality of each sequence read was checked by trim-galore (V.0.6.6), and a base with Phred score less than 20 was removed from further analysis. For each 16S rRNA dataset, all reads were input into QIIME2 as demultiplexed sequences, which were then denoised to generate amplicon sequence variants (ASVs) using DADA2 (Callahan et al., 2016). The cluster-features-closed-reference method was used to pick up OTUs. The classify-consensus-vsearch method was applied to assign taxonomy using 99% Greengenes database 13.8 and 99% SILlVA 138 as reference. Recent Shotgun metagenomic data (PRJNA433459) (Qian et al., 2020), which were re-analyzed with Kraken2 in this study, were further used to verify the results obtained from 16S rRNA. For each RNA-Seq dataset, all reads were classified with Kraken2 using minikraken_8GB_20200312 database as reference, and Bracken was then applied to count the relative abundance of taxonomy based on the read length (Wood et al., 2019).

The same RNA-Seq data used above were also applied to detect gene expression levels in the brain. RNA-Seq and microarray data were processed with two web tools, GREIN³ (Al Mahi et al., 2019) and ImaGEO⁴ (Toro-Domínguez et al., 2019), respectively. Briefly, unique GSE ID obtained from the GEO was used as query in GREIN and ImaGEO individually, and data were processed with the default parameters. A meta-analysis was conducted automatically to identify DEGs in the ImaGEO after microarray data were processed. All statistical analyses were conducted with R (version 4.0.4) (Chambers, 2008) unless otherwise stated, and P value less than 0.05 was considered statistically significant.

To alleviate batch effects (hypervariable region, sequencing platform, and country, etc.), microbiota analysis was performed at the genus level. To exclude the impacts of bacterial contamination, those genera previously reported as potential contaminants or widely distributed in the tubes and agents (Supplementary Data 1) (Salter et al., 2014; Davis et al., 2018; Eisenhofer et al., 2019) were removed. To identify contaminants in blank tubes and agents, Salter et al. routinely process the a negative ‘blank’ controls without adding sample template alongside human-derived samples by 16S rRNA and shotgun metagenomics (Salter et al., 2014). A list of contaminant genera was detected in sequenced negative ‘blank’ controls (Supplementary Data 1), which were removed from microbiota analysis in our study. In addition, as gut microbiota is the main source of bacteria in human, we also excluded genera (“pan-contaminants”) exclusively found in the blood and brain. In this study, all data analyses were conducted using the microbiota genera data after removing pan-contaminants unless otherwise stated. The relative abundance of bacteria at the Phylum level was calculated and plotted using the tax_stackplot within amplicon (V.1.1.5) package (Liu Y.-X. et al., 2021). A heatmap of genera with a relative abundance equal or more than 0.5% was also plotted for each sample types (feces, blood, and brain) using the pheatmap (V.1.0.12) package (Kolde and Kolde, 2015). Since the sequencing depth and genera abundance among feces, blood, and brain samples was different, rarefied genera data were only used for diversity analysis. For the α-diversity analysis in the feces and blood or brain (RNA-Seq), 3,000, and 1,000 counts were extracted randomly from each sample without replacement to estimate Observed, Chao1, Shannon, and Simpson evenness indices using the phyloseq (V.1.34.0) package (McMurdie and Holmes, 2013). For β-diversity analysis, the relative abundance of genera profiles were used to calculate β-diversity with the following parameters (dist = “bray,” method = “PCoA,” Micromet = “Adonis”) in amplicon package (Liu Y.-X. et al., 2021).

To assess the difference in microbiota function between PD patients and controls, we used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) to determine the predicted function of gut microbial communities from 16S rRNA sequencing data (Douglas et al., 2020). The KEGG Orthology (KO) metagenome prediction data were obtained using the picrust2_pipeline.py script in PICRUSt2. After that, the predicted KO abundances data were categorized into KEGG pathway abundances of Level 2 and Level 3. Finally, a meta-analysis was performed to identify the different KEGG pathways of Level 2 and Level 3 between PD patients and controls.

The genus-level abundance data were processed with qiime2R (V.0.99.4) and phyloseq package (Bisanz, 2018). Since pooling summary statistics from each study was more robust in alleviating between-study heterogeneity than pooling raw data directly from all studies (Duvallet et al., 2017), the genus-level abundance data from each study were analyzed with DESeq2 (Love et al., 2014). The summary statistics (log2 fold-changes and 95%CIs) of a single genus from each study were then pooled with a fixed-effect meta-analysis using the meta (V.4.18-2) package (Schwarzer and Schwarzer, 2012). Due to only one study available for blood 16S rRNA, no meta-analysis was performed. The significantly different genera between PD patients and controls across brain, blood, and feces samples were summarized to obtain the common genera/genus with the same change direction, and a corresponding heatmap was plotted using the ggplot2 (V.3.3.3) package (Wickham, 2011).

For gene expression data with RNA-Seq (GREIN processed), raw counts data from each study were used to calculate log2 fold-changes of each gene using the DESeq2 (V.1.30.1) package (Love et al., 2014). After that, a meta-analysis was performed to identify the DEGs between PD patients and controls using the meta package (Schwarzer and Schwarzer, 2012).

To explore the potential interaction between Blautia genus and host genes, a correlation analysis was performed between the Blautia genus and each host gene (RNA-Seq) using the corrplot (V.0.89) package with Spearman method (Wei et al., 2017). Correlation coefficients with a magnitude above 0.3 were considered to have a correlation.

Host genes associated with Blautia genus (r > 0.3) were enrolled in enrichment analyses, including Gene Ontology (GO) (Gene Ontology Consortium, 2004), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), and Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) using the clusterProfiler (V.3.18.1) (Yu et al., 2012) and enrichplot (V.1.10.2) package (Yu, 2018). Briefly, the gene symbols with combined effect sizes were input as a query gene list and a P value cutoff was set to 0.05. Meanwhile, further enrichment analyses were also conducted based on DEGs (RNA-Seq or RNA-Seq and Microarray) correlated to Blautia genus. Moreover, to explore the relationship between Blautia-associated genes and diseases category (Piñero et al., 2016), Blautia-associated DEGs (overlapped between RNA-Seq and microarray) were also used to conduct diseases enrichment analysis with the disgenet2r (V.0.99.2) package in DisGeNET database (Gutiérrez-Sacristán et al., 2016).

To verify the importance of Blautia genus in discriminating PD patients from controls, we used the microbiota data to build a random forest (RF) model with stratified 5-fold cross-validation (CV) using the SIAMCAT (V.1.10.0) package (Wirbel et al., 2021). The features used for model building consist of patients’ metadata as well as genera with a relative abundance greater than 0.1% across all the samples. The metadata features consisted of age, gender, body mass index (BMI), or hypervariable region (16S rRNA). The RF models were built with 1,000 estimator trees and default parameters to select the proportion of features for each tree. The stratified 5-fold cross-validation was executed to configure training and testing data sets. The top 30 features from the top-performing model were selected as “important features.” Finally, the receiver operating characteristic curve (ROC) and precision-recall curve (PRC) with area under curve (AUC) was used to assess the performance of models.

A total of 16 studies, including ten for feces, one for blood, and five for the brain, with 2,883 participants were included in this study before data cleaning. All the Project IDs and details of metadata were provided in Table 1 (Supplementary Table 1 and Supplementary Data 2). Genera reported as potential contaminants and non-detectable in fecal microbiome were termed as “pan-contaminants” and were removed before microbiota analysis. For the feces, 311 genera were identified across 2,252 individual samples. After removing contaminants (Supplementary Data 3), 197 genera retained and 28 genera had an average relative abundance greater than 0.5% across all samples (Figure 2A). For the blood, 328 genera were identified across 86 individual samples with 100 retained after excluding pan-contaminants (Supplementary Data 4), and 14 genera had an average relative abundance greater than 0.5% across all samples (Supplementary Figure 1A). For the brain, 1,197 genera were identified across 376 individual samples with 116 retained after removing pan-contaminants (Supplementary Data 5), and 18 genera had an average relative abundance greater than 0.5% across all samples (Supplementary Figure 1B). For fecal studies, samples overlapping between PRJEB14674 and PRJNA601994 were de-duplicated, and 1,914 samples were finally included. There were four replicated runs of each sample in PRJNA557205, and the RNA-Seq data from the first run (run accession IDs) of each sample were used in both microbiota analysis and transcriptome analysis (Li et al., 2020).

Microbiota composition data without reported contaminants showed that the main flora in feces were Firmicutes and Bacteroidetes, while major flora in blood and brain were Proteobacteria, followed by a small part of Actinobacteria (Figure 2B). In fecal samples, patients with PD harbored a higher relative abundance of Firmicutes, Actinobacteria, and Synergistetes, but a lower relative abundance of Bacteroidetes (Figure 2B and Supplementary Figure 2A), contributing to an increased ratio of Firmicutes to Bacteroidetes (F/B) in the feces of patients with PD. The Firmicutes phylum in PD brain were almost depleted after excluding the pan-contaminants (Supplementary Figure 2A).

The results of α-diversity analysis revealed that the α-diversity of fecal microbiota was significantly higher in PD patients than controls (Figure 2C), but no significant difference was found in the blood and brain (Supplementary Figures 2B,D). For β-diversity, the results of principal coordinates analysis (PCoA) with Bray–Curtis distance showed that the bacterial communities in the feces differed significantly between PD patients and controls (Figure 2D) (p < 0.001, PERMANOVA by Adonis), but not in the blood and brain (Supplementary Figures 2C,F).

The results of PICRUSt2 predictions showed that top ten altered KEGG pathways of Level 2 were mainly associated with neurodegenerative diseases, infectious diseases, cancer, transport, and catabolism etc. (Supplementary Figure 3A). Data from the top ten altered KEGG pathways of Level 3 showed an increased abundance of betalain biosynthesis, amoebiasis, and hepatitis C, while decreased in beta-lactam resistance, Ubiquitin system, isoquinoline alkaloid biosynthesis, beta-alanine metabolism, etc. (Supplementary Figure 3B). A full list of all the changed KEGG pathways of both Level 2 and Level 3 were presented in Supplementary Tables 2, 3, respectively.

In the feces, the pooled results of taxonomical classification by Greengenes showed that 32 genera were significantly different between PD patients and controls (P < 0.05, Supplementary Table 4). Nine of the 32 different genera were depleted, while 23 were enriched in PD. Meanwhile, the pooled results of taxonomical classification by SILVA showed that there were 23 genera significantly altered in PD patients (P < 0.05, Supplementary Table 5 and Supplementary Data 6), and 15 of them overlapped with the different genera identified using Greengenes alignment (Supplementary Figure 4). Using shotgun metagenomics data (PRJNA433459), we found that 377 genera were significantly different between PD patients and controls (P < 0.05, Supplementary Data 7), with 251 genera depleted and 126 enriched in PD, respectively. In the blood, the results showed that all four different genera were significantly decreased in PD patients compared to controls (P < 0.05, Supplementary Table 6). In the brain, a total of 13 genera was found to be significantly different between PD patients and controls (P < 0.05, Supplementary Table 7), with 8 genera depleted and 5 genera enriched in PD (Figure 3).

The Halomonas and Finegoldia genus was depleted in the brain, but enriched in the feces of PD patients (P < 0.05) (Figure 3). Meanwhile, the Dorea genus was increased in the blood, but decreased in the feces (P < 0.05). The Bacteroides genus were depleted in the feces and brain of patients with PD (P < 0.05) (Figure 3). After comparing the different genera among feces, blood and brain samples, we found that only Blautia genus consistently altered (depleted) across all sample types in PD (P < 0.05) (Figure 3). Further analysis at species level using gut shotgun metagenomics data showed that Blautia argi, Blautiacoccoides, Blautia sp. SC05B48, and Blautiahansenii species were decreased in the feces of patients with PD as compared to controls (Supplementary Table 8).

To investigate the potential role of Blautia genus in PD, Spearman correlation analyses between Blautia genus and brain genes were conducted by integrating the microbiome and host transcriptome data obtained from the same raw RNA-Seq data. The results showed that 17,216 genes were associated with Blautia genus (P < 0.05), with an absolute r value ranging from 0.12 to 0.479. After removing those genes with an absolute r value less than 0.3, we found 3,583 genes remained associated with Blautia genus, with 1,775 and 1,808 genes having positive and negative relationship with Blautia genus (Supplementary Data 8), respectively. Biotype distribution analysis using the biomaRt (V.2.46.3) package showed that, besides protein coding class, there were 78 and 18 genes mapped to lncRNA and miRNA class, respectively (Supplementary Table 9). The top nine genes labeled with gene symbol in Figure 4A had an absolute r value greater than 0.45. Eight of the nine genes belong to protein encoding gene, while only MALAT1 was a lncRNA.

To assess whether host genes associated with Blautia genus were involved in the pathogenesis of PD, we then implemented three functional enrichment analyses, including KEGG, GO, and GSEA. The results of KEGG analysis showed that these genes were mainly related neurodegenerative diseases, including PD, AD, Huntington’s disease (HD), and amyotrophic lateral sclerosis (ALS), and metabolism related diseases, including diabetic cardiomyopathy (DCM) and non-alcoholic fatty liver disease (NAFLD) (Figure 4B). Data from GO enrichment analysis showed that these genes were mainly mapped to proteins transportation, immune response, energy metabolism, and mitochondrial function (Supplementary Figure 5). Meanwhile, the above results were further confirmed by GSEA-based KEGG and GO enrichment analysis (Supplementary Figures 6A,B).

To identify the DEGs in the brain of PD patients, a meta-analysis was performed using RNA-Seq and microarray data, respectively. Using brain RNA-Seq data, we found that 9,563 DEGs were significantly different between PD patients and controls (Supplementary Data 9), with 4,618 decreased genes and 4,945 increased genes, respectively. A total of 1,393 DEGs were significantly associated with Blautia genus (r > 0.3 and P < 0.05) (Figure 5A), with 569 downregulated genes and 824 upregulated genes (Figure 5B) individually. Among the 1,393 DEGs, 759 DEGs were positively correlated with Blautia genus, while 634 DEGs were negatively associated with Blautia genus (Supplementary Data 10).

Using brain microarray data, a total of 3,237 genes were detected in microarray and over 98% of the genes were also identified in RNA-Seq (Supplementary Figure 7). Only 67 genes were found exclusively in microarray, but 32,474 gene transcripts identified in RNA-Seq were not covered in microarray (Supplementary Figure 7A). There were 460 DEGs significantly changed in patients with PD using microarray data (Supplementary Data 11). Finally, only 36 DEGs correlated to Blautia genus were found in both RNA-Seq and microarray data (Figure 5A), and 26 out of the 36 genes were negatively associated with Blautia genus (Supplementary Figure 8), respectively.

Since 1,357 DEGs associated Blautia genus were exclusively identified in RNA-Seq data, but not covered in microarray data, thus, in order to avoid losing potential important information, we first performed an enrichment analysis based on RNA-Seq derived DEGs. The results of GO enrichment analysis showed that the 1,393 DEGs associated Blautia genus were mainly located in mitochondrial membrane, participating in regulating energy metabolism and mitochondrial function (Supplementary Figure 7B and Supplementary Data 12). Data from KEGG analysis demonstrated that the 1,393 DEGs were mainly mapped to PD, HD, DCM, and NAFLD (Figure 5C and Supplementary Data 13), which was further confirmed by GSEA-based KEGG analysis (Supplementary Figures 9A–D and Supplementary Data 14). Meanwhile, GSEA-based GO analysis showed that these DEGs were mainly targeted in inflammatory response (Figures 6A,B), mitochondrial function (Figure 6C), and miRNA binding (Figure 6D and Supplementary Data 15).

Moreover, we also conducted a GO analysis using the 36 DEGs overlapping between RNA-Seq data and microarray data. The result showed that the 36 DEGs were mainly mapped to pH regulation, protein degradation, and mitochondrial function (Supplementary Figure 10). Several of the 36 DEGs, including KCNJ6, ACOX1, USP9X, SMAD2, SLC8A1, PIK3C3, PDK1, have been previously reported to be implicated in the pathogenesis of PD (Rott et al., 2011; di Val Cervo et al., 2017; Muñoz et al., 2020). Disease enrichment analysis performed in Disgenet2 database showed that the 36 DEGs were mainly associated with nervous system diseases and digestive system diseases (Supplementary Figure 11).

To further manifest the role of Blautia genus in discriminating PD from control, we constructed a stratified 5-fold cross-validation RF model using the microbiota data of feces, blood, and brain, individually. The top 30 important genera contributing to the classifier (RF models, over 80% weights) in the feces were all decrease in PD. The Blautia genus was one of the important genera in this classifier (Figure 7A), which was further confirmed in the blood and brain (Supplementary Figures 12A, 13A). The classifier achieved an AUC of 0.704 and 0.787 for ROC and PRC in discriminating PD patients from controls (Figures 7B,C), respectively. The AUC value of ROC and PRC was 0.777 and 0.768 for the RF model using microbiota data from blood samples (Supplementary Figures 12B,C), and the AUC value of ROC and PRC was 0.689 and 0.701 for the RF model using microbiota data from brain samples (Supplementary Figures 13B,C).

Some limitations in our study should be addressed here. First, this is the first study using the Kraken2 to extract microbiota reads in the brain RNA-Seq data. Although Kraken2 has been applied to identify microbiota reads in tumor RNA-Seq data with high sensitivity and specificity (Poore et al., 2020), it is still hard to rule out the possibility of a false positive result owing to bacterial contamination. We have tried our best to remove the potential contamination genera to minimize the impact of bacterial contamination on the final results. Second, it is worthy to note there is obvious heterogeneity among the collected microbiome datasets across studies. As no optimal way is currently available to correct batch effects, only genus-level data was used in this study for microbiome analysis to alleviate heterogeneity (Wang et al., 2020). This leads to a loss of many OTUs, which may affect the results of diversity analysis. Third, microbiota identified in the blood and brain might originate from the bacterial debris, whose DNA/RNA can also transport in the blood and brain. However, there is no effective method now to distinguish them from live bacteria. Forth, the transcriptome and microbiota data obtained in this study were based on the same RNA-Seq data and the correlation threshold (r > 0.3) was relatively low, which may contribute to an over-estimation of the correlation and enrichment analysis. Fifth, since we use “pan-contaminants” to remove bacteria not found in the feces, these may result in the loss of a lot of information and introduce potential bias to the final results. Sixth, only one study was included for microbiota analysis in blood samples. Taken together, our results in this study need to be interpreted very cautiously.