Gut microbiota in Parkinson's disease: Temporal stability and relations to disease progression.

BACKGROUND: Several publications have described differences in cross-sectional comparisons of gut microbiota between patients with Parkinson's disease and control subjects, with considerable variability of the reported differentially abundant taxa. The temporal stability of such microbiota alterations and their relationship to disease progression have not been previously studied with a high-throughput sequencing based approach.
METHODS: We collected clinical data and stool samples from 64 Parkinson's patients and 64 control subjects twice, on average 2·25 years apart. Disease progression was evaluated based on changes in Unified Parkinson's Disease Rating Scale and Levodopa Equivalent Dose, and microbiota were characterized with 16S rRNA gene amplicon sequencing.
FINDINGS: We compared patients to controls, and patients with stable disease to those with faster progression. There were significant differences between microbial communities of patients and controls when corrected for confounders, but not between timepoints. Specific bacterial taxa that differed between patients and controls at both timepoints included several previously reported ones, such as Roseburia, Prevotella and Bifidobacterium. In progression comparisons, differentially abundant taxa were inconsistent across methods and timepoints, but there was some support for a different distribution of enterotypes and a decreased abundance of Prevotella in faster-progressing patients.
INTERPRETATION: The previously detected gut microbiota differences between Parkinson's patients and controls persisted after 2 years. While we found some evidence for a connection between microbiota and disease progression, a longer follow-up period is required to confirm these findings.

Evidence before this study

The cause of idiopathic Parkinson's disease has been hypothesized to involve an external agent, for example a pathogen, and one potential route of entry for such an agent could be via the gastrointestinal system. The interactions of gut microbes and the central nervous system, also known as the microbiota–gut–brain axis, have recently become a topic of intense research. Within the past 5 years, a total of twelve studies from three continents have reported differences in the composition of gut microbiota of patients with Parkinson's disease when compared to non-parkinsonian control subjects, showing promise for this novel field of research.

Added value of this study

Our study is the first to use a high-throughput sequencing based approach to explore gut microbiota of Parkinson's patients at two different timepoints. We show that the differences detected at baseline can be replicated at a follow-up timepoint 2 years later, and that there might be changes in gut microbiota composition in patients with faster disease progression.

Implications of all the available evidence

The consistent differences in gut microbiota between Parkinson's patients and control subjects could lead to new diagnostic or therapeutic modalities.

Alt-text: Unlabelled Box

Introduction

The early non-motor symptoms of Parkinson's disease (PD), such as hyposmia and gastrointestinal (GI) disorders [1,2], have led to the hypothesis that the disease could originate outside of the central nervous system, for example in the olfactory bulb or the enteric nervous system [3]. Research comparing nasal microbiota of PD patients and control subjects has not revealed notable differences [4,5]. In contrast, several studies have suggested that patients' gut microbiota differ from controls' [[5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16]], although the differentially abundant taxa reported in them vary considerably. This could be due to differences in subject populations or methodology, such as PCR primers, sequencing platforms, and statistical tools. Nevertheless, some microbial community alterations, including a decreased abundance of the family Prevotellaceae, the genus Prevotella, and the species Prevotella copri [6,9,12,13], and an increase in Akkermansia and Verrucomicrobiaceae [[5], [6], [7],9,11,12,15], Bifidobacterium and Bifidobacteriaceae [9,11,13], as well as Lactobacillus and Lactobacillaceae [6,11,13], have been detected multiple times.

Aside from a recent disease progression study using a qRT-PCR-based assay [17], all PD gut microbiota publications have been case-control studies with one timepoint. They have used varying approaches to control for the effects of potential confounders, such as diet and medications. Diet influences the gut microbiome [18]. Since it has been hypothesized to be an important determinant of the abundance of Prevotellaceae [18], it could be associated to the decrease of that family seen in PD [6,7]. Additionally, PD medications can affect gut microbiota [6,11]. They are a particularly important confounder when studying disease progression, since progression is measured based on symptoms, which respond to medications.

In the present study, we explore the gut microbiota of a previously recruited group of PD patients and control subjects at baseline and 2 years later, while also considering their diet, medications, and other clinical variables.

Materials and methods

Study subjects and clinical data

152 age and sex matched subjects (76 PD patients, 76 control subjects) originally recruited for a pilot study in Parkinson's disease and gut microbiota [6] were invited to a follow-up appointment on average 2·25 (SD: ± 0·20) years later. Out of the original subjects who returned, nine were excluded, while five subjects whose samples were not used in the pilot study for various reasons were included at follow-up, bringing the total number of subjects to 128 (64 PD patients, 64 control subjects; Table 1, Fig. 1). The study was approved by the ethics committee of the Hospital District of Helsinki and Uusimaa. All participants gave informed consent.

Table 1

Changes in subject inclusion/exclusion after pilot study.

n (recruited for pilot study)	152
Excluded from pilot study, not included in current study:
Excluded from pilot due to insufficient read count, drop-out from follow-up	2
Excluded from pilot due to lack of matching subject	1
Excluded from pilot study, included in current study:
Not in pilot study due to insufficient read count	2
Not in pilot study because sample received too late	2
Not in pilot study: nasal polyps; deemed eligible for follow-up study	1
n (pilot study)	144
Excluded after pilot study:
Exitus after pilot study	1
Drop-outs	11
Clinical exclusion: diagnosis changed to Lewy body dementia	1
Clinical exclusion: recent surgery	3
Technical exclusion: insufficient reads from baseline sample	1
Technical exclusion: missing sample	1
Outlier exclusion: microbial community outlier	3
n (follow-up study)	128

Fig. 1

NMDS ordination of samples excluded as microbiome outliers.

The subjects filled several questionnaires concerning non-motor symptoms, including the Wexner constipation score [19], Rome III IBS questionnaire [20], Non-Motor Symptoms Scale (NMSS) [21], Swallowing Disturbance Questionnaire (SDQ) [22], Sialorrhea Clinical Scale for PD (SCS-PD) [23], 15-item Geriatric Depression Scale (GDS-15) [24], REM sleep behavior disorder screening questionnaire (RBDSQ) [25] and Sniffin’ Sticks 16-item smell identification score [26]. The subjects' dietary habits were evaluated at follow-up with a 163 item Food Frequency Questionnaire (FFQ) with 9 frequency response options (based on [27]). The severity of the patients' parkinsonian symptoms was assessed with the Unified Parkinson's Disease Rating Scale (UPDRS) [28], and their total medication load was calculated using the Levodopa Equivalent Dose (LED) [29]. We also determined tremor and postural instability and gait difficulty (PIGD) symptom scores and derived motor phenotypes (postural instability and gait difficulty (PIGD), tremor dominant (TD), or mixed (MX)) [30]. At follow-up, we also collected stool consistency information using the Victoria Bowel Performance Scale (BPS) [31]; the patients kept a diary of their scores for one week leading up to sampling, and we used the weekly average stool consistency values for comparisons.

For comparisons of microbiota and disease progression, we excluded patients who were on Deep Brain Stimulation at either timepoint, those with missing UPDRS or LED values, and one patient whose LED score had decreased considerably between the two timepoints due to adjustment of overmedication symptoms. This left a subset of 56 PD patients for progression analyses. To classify patients into stable or progressed, we calculated changes in UPDRS I-III score (in the ON state) and LED between baseline and follow-up, divided each value by the number of days between appointments, z-transformed these two variables, and added them up. Based on the distribution of samples on this progression scale, we chose the 3rd quartile as a cut-off, resulting in 41 stable and 15 progressed patients (Fig. 2). We also used this subset of patients for additional analyses contrasting PD phenotypes, additionally excluding patients with a mixed (MX) phenotype, which resulted in a data set with 21 TD patients and 28 PIGD patients at baseline, and 20 TD patients and 35 PIGD patients at follow-up.

Fig. 2

Histogram of the progression variable based on change in UPDRS I-III and LED.

Legend: Solid vertical line represents the median, dashed vertical line represents the 3rd quartile, which was used to categorize subjects into “stable” and “progressed”.

Sequencing and sequence analysis

Stool samples for both timepoints were collected at home by the study subjects into collection tubes with pre-filled DNA stabilizer (PSP Spin Stool DNA Plus Kit, STRATEC Molecular), and stored in the refrigerator until transport (for up to 3 days). Once received at the clinic, they were transferred to −80 °C. To minimize potential technical differences between baseline and follow-up samples, we reanalysed the baseline samples together with the follow-up samples, including new DNA extractions, PCR, and sequencing. Thus, the baseline samples, which had been frozen after collection, then thawed for sequencing at the time of the pilot study, re-frozen, and stored at −80 °C since, were thawed for a second time for this follow-up study.

We extracted bulk DNA from stool samples with the PSP Spin Stool DNA Plus Kit (STRATEC Molecular). Each extraction batch included one blank sample to assess potential contamination. The V3-V4 regions of the 16S rRNA gene were amplified following a previously published protocol [4], with the following changes: we used two technical replicates (25 μL reactions) per patient sample, and a mixture of the universal bacterial primers 341F1–4 (5′ CCTACGGGNGGCWGCAG 3′) and 785R1–4 (5′ GACTACHVGGGTATCTAATCC 3′) [4] with partial Illumina TruSeq adapter sequences added to the 5′ ends (F1; ATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT, F2; ATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTgt, F3; ATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTagag, F4; ATCTACACTCTTTCCCTACACGACGCTCTTCCGATCTtagtgt and R1; GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT, R2; GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTa, R3; GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTtct, R4; GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTctgagtg). The additional nucleotides (small letters) are introduced for mixing in sequencing. The two-step PCR and subsequent quantification, pooling, and purification were done as described previously [4]. The obtained PCR amplicon pool was checked using Fragment Analyzer (Advanced Analytical Technologies Inc., Ankeny, IA, USA). Every PCR batch included a blank sample (no added DNA template) to assess potential contamination. Finally, the PCR products were sequenced with Illumina MiSeq (v3 600 cycle kit), with 325 bases for the forward and 285 bases for the reverse read.

The raw sequence data contained a total of 34 701 899 sequence reads. These sequences are available at the European Nucleotide Archive with the accession number PRJEB27564. Primers were removed from the reads using cutadapt (version 1.8.3) [32]. The reads were paired, quality trimmed, taxonomically classified, and clustered into OTUs following mothur's Standard Operating Procedure (SOP) for MiSeq (make.contigs: mothur 1.38.1, the rest of the workflow: mothur 1.39.5; SOP version last updated on 4 April 2018) [33,34]. The following changes were made to the SOP parameters: insert = 40 and deltaq = 10 in make.contigs; maxlength = 450 in the first screen.seqs step; start = 6428 and end = 23440 in the second screen.seqs step; diffs = 4 in pre.cluster. Additionally, singleton sequences were removed with split.abunds (cutoff = 1) before running classify.seqs. The references used were the full-length SILVA alignment release 128 for align.seqs and the RDP 16S rRNA reference (PDS) version 16 for classify.seqs. We inspected the extraction and PCR blank samples, which overall had low amounts of sequence reads (Supplementary material: R Markdown), and since these did not suggest any overall problems with contamination, we deleted these samples before downstream analyses. The final sequence data set (without the blanks) consisted of 18 867 278 good quality reads (median [IQR]: 73 078 [50 032–98 685]).

Statistics

All statistical comparisons and data visualization were performed with R (v. 3.5.1) [35]. The full analysis workflow is available in the supplementary materials (Supplementary material: R Markdown). In all comparisons, p-values <0·05, or adjusted p-values <0·05 in the case of multiple comparison corrected tests, were considered significant. False discovery rate [36] was used for multiple comparison correction. For basic comparisons of potentially confounding clinical variables, we used either Student's t-test, Wilcoxon signed rank test or Fisher's exact test depending on the type and distribution of each variable, and these comparisons were not corrected for multiple comparisons. Figures were plotted with ggplot2 (v. 3.1.0), except for the diet principal component analysis (PCA) biplot (which additionally utilized ggfortify v. 0.4.5) and Euler plots (drawn with eulerr v. 5.0.0).

Food and nutrient variables from the FFQ were adjusted for energy intake (divided by total energy intake in kilocalories and multiplied by 1000). For variable-specific comparisons, these continuous variables were split into categories by quintiles. To look for dietary patterns, we used PCA of a hand-picked set of 31 non-overlapping, continuous, energy-adjusted and z-transformed food items.

The R-package phyloseq (v. 1.26.0) [37] was used for microbiota data handling and calculating alpha diversity indices (observed richness, Shannon index and inverse Simpson index). Wilcoxon rank-sum test and linear regression were used for statistical testing of alpha diversity differences. Differences in Firmicutes/Bacteroidetes and Prevotella/Bacteroides ratios were also compared with the Wilcoxon rank-sum test. Enterotyping was run with the reference-based online tool [38], and distributions of enterotypes were compared with the chi-square test, for which simulate.p.value = TRUE was used when comparing PD patients subsetted according to progression, as the group sizes were small.

Beta-diversity comparisons were done with vegan (v. 2.5–3) [39] using Bray-Curtis dissimilarity of data subsampled to the lowest amount of sequences in a sample (2201), on three different taxonomic levels (OTU, genus and family). PERMANOVA was run with the command adonis2, with the parameters by = “margins” and perm = 9999 for all comparisons, except for the Timepoint + Parkinson + single confounder as well as the diet variable tests, which were run with 999 permutations. To narrow down the lists of potential confounders for testing, we focused on those clinical variables that differed significantly between groups (PD/control or progressed/stable), and a few common confounders that did not (age, sex, body mass index). Variables that were correlated with PD status (|Pearson's r| ≥ 0·5; Table 2) were only used for the within PD group comparisons; additionally, GDS 15 (Pearson's r = 0·491, close to the cutoff and with high values mainly present in PD patients at follow-up) and catechol-O-methyl transferase (COMT) inhibitor use (Pearson's r = 0·307, but with no controls using the medication) were also only used in the within PD group comparisons. To explore potential confounders for the PD/control comparisons, we first tested each variable for microbial community effects in a model that included timepoint (baseline or follow-up), PD status, and the variable in question. The variables that were significant (adonis2 p < 0·05 on at least one taxonomic level) in these single-confounder comparisons were then tested together in a combined adonis2 model. As an alternative test, we used the envfit function (run with permutations = 9999), fitting variables that were significant on at least one level with adonis2 onto a Nonmetric Multidimensional Scaling (NMDS) ordination performed with the metaMDS function (run with try = 500). Beta-diversity analysis for the progression comparisons was performed in a similar manner, with the categorical disease progression variable (progressed/stable) in place of the PD status variable.

Table 2

Correlations with PD status for potential microbiome confounders.

Variable	Pearson correlation coefficient with PD	Use for PD status comparisons	Use for progression comparisons
Rome III 9–15 sum score	0.432	Yes	No
RLS	0.419	Yes	No
Wexner total	0.418	Yes	No
Rome III IBS criteria fulfilled	0.280	Yes	No
Medication anticholinergic	0.161	Yes	No
Diet PC1	0.145	Yes	No
BMI	0.094	Yes	No
Age at stool collection	0.061	Yes	No
Sex	0.016	Yes	No
MMSE total	−0.159	Yes	No
Medication ACE-I/ARB	−0.160	Yes	No
Medication Ca channel blockers	−0.165	Yes	No
Medication warfarin	−0.205	Yes	No
Medication statin	−0.332	Yes	No
History TIA / ischemic stroke	−0.364	Yes	No
Medication dopamine agonist	0.807	No	Yes
LED	0.753	No	Yes
Medication MAO inhibitor	0.736	No	Yes
NMSQuest total	0.714	No	Yes
Medication dopa	0.662	No	Yes
NMSS total	0.607	No	Yes
RBDSQ	0.580	No	Yes
SDQ total	0.529	No	Yes
SCS PD total	0.520	No	Yes
GDS 15	0.491	No	Yes
Medication COMT inhibitora	0.307	No	Yes
Sniffin’ Sticks	−0.769	No	Yes

Table legend: RLS: Restless Legs Syndrome, IBS: Irritable Bowel Syndrome, PC: Principal Component, BMI: Body Mass Index, MMSE: Mini-Mental State Examination, ACE-I/ARB: angiotensin-converting-enzyme inhibitor / angiotensin II receptor blocker, TIA: Transient Ischemic Attack, LED: Levodopa Equivalent Dose, MAO: Monoamine Oxidase, NMSQuest: Non-Motor Symptoms Questionnaire, NMSS: Non-Motor Symptoms Scale, RBDSQ: REM Sleep Behavior Disorder Screening Questionnaire, SDQ: Swallowing Disturbance Questionnaire, SCS-PD: Sialorrhea Clinical Scale for PD, GDS 15: 15-item Geriatric Depression Scale, COMT: catechol-O-methyl transferase.

COMT inhibitor variable used only for progression although |Pearson's r| < 0.5, since only PD patients take this medication.

For differential abundance comparisons, we used ANCOM (v. 2.0; unpublished version shared online) [40,41], DESeq2 (v. 1.22.1) [42], and random forests ([43], packages randomForest (v. 4.6–14) [44] for the classification, rfUtilities (v. 2.1–3) [45] for estimating the significance of the classification, and rfPermute (v. 2.1.6) [46] for assessing the significances of specific taxa). All comparisons were performed on OTU, genus, and family levels, with data trimmed to include taxa that had more than one read in at least 1/10 of samples (26 samples for PD status comparisons, 11 for progression comparisons, and 10 for PD phenotype comparisons); additionally, OTUs were required to have at least 1000 sequence reads altogether. ANCOM and random forests were run separately for baseline and follow-up, classifying data by PD status or progression category. For ANCOM, we adjusted the PD/control comparisons for Rome III score and BMI, and the progression comparisons for COMT inhibitor use, and chose the less stringent multiple comparison correction option (multcorr = 2) and the 0.90 cutoff (prev.cut = 0.90) for proportion of zeroes. The results for the 0.6 detection level column were considered significant. DESeq2 comparisons were corrected for the same confounders and were run with the parameters fitType = “parametric” and sfType = “poscounts”. For the PD status comparisons, we used a model that was additionally corrected for subject (model: Rome III score + BMI + PD : subject + timepoint * PD) and extracted timepoint-specific contrasts from the results. Subjects that lacked baseline BMI information (2 PD, 3 control) were excluded from this analysis, which resulted in an unbalanced data set. Because of this and to assess the robustness of the results, the PD status DESeq2 comparison was performed with a leave-one-out approach, excluding each one of the 62 PD patients in turn, with the average FDR-adjusted p-value of the 62 rounds as the final result. DESeq2 comparisons for progression were run separately for each timepoint, correcting for COMT inhibitor medication (model: COMT + Progression). PD phenotype comparisons (TD vs PIGD, excluding patients with MX type) were run with ANCOM and DESeq2, separately for each timepoint, and not corrected for confounders due to the small number of subjects included in this subgroup analysis.

Results

Clinical and diet data

Control subjects and PD patients were matched for age and sex in the pilot study [6], and those who returned for the follow-up still had similar distributions for these variables, as well as body mass index (BMI) (Table 3). As expected, PD patients had higher scores on the non-motor symptoms scale for Parkinson's disease (NMSS, p < 0·001 at both timepoints), and the Rome III IBS questionnaire (p < 0·001 at both timepoints). Controls more commonly had a history of transient ischemic attack (TIA) or ischemic stroke (p < 0·001 at both timepoints) and were on several medications more often than patients: statins (p < 0·001 at both timepoints), warfarin at baseline (p (baseline) = 0·017, p (follow-up) = 0·076), and calcium channel blockers (CCB) at follow-up (p (baseline) = 0·271, p (follow-up) = 0·035) (Table 3).

Table 3

Clinical variables in Parkinson's patients and control subjects at each timepoint.

Variable	Timepoint	Control subjects (% (n) / mean ± SD / median [IQR])	Parkinson's patients (% (n) / mean ± SD / median [IQR])	p-value	Test
n		64	64
Age (at first stool collection)		64.45 ± 6.9	65.2 ± 5.52	0.499	t
Sex (% males)		50 (32)	51.56 (33)	1.000	Fisher
BMI	Baseline	26.23 [24.1–28.05]	26.51 [24.25–29.36]	0.319	Wilcox
	Follow-up	26.94 [24.32–28.64]	27.24 [23.95–30.08]	0.572	Wilcox
History of TIA or ischemic stroke	Baseline	37.5 (24)	6.25 (4)	<0.001	Fisher
	Follow-up	37.5 (24)	7.81 (5)	<0.001	Fisher
Medication: calcium channel blockers	Baseline	15.62 (10)	7.81 (5)	0.271	Fisher
	Follow-up	20.31 (13)	6.25 (4)	0.035	Fisher
Medication: statins	Baseline	56.25 (36)	21.88 (14)	<0.001	Fisher
	Follow-up	50 (32)	20.31 (13)	<0.001	Fisher
Medication: warfarin	Baseline	14.06 (9)	1.56 (1)	0.017	Fisher
	Follow-up	15.62 (10)	4.69 (3)	0.076	Fisher
NMSS total	Baseline	8 [4–12]	40 [23.75–55]	<0.001	Wilcox
	Follow-up	6 [2–10.25]	40 [19.75–58.75]	<0.001	Wilcox
Rome III constipation-defecation score (sum of items 9–15)	Baseline	2 [1–4]	6 [2.75–11]	<0.001	Wilcox
	Follow-up	2 [0–3]	8 [2.75–11]	<0.001	Wilcox
Rome III IBS criteria fulfilled	Baseline	7.81 (5)	23.44 (15)	0.027	Fisher
	Follow-up	7.81 (5)	35.94 (23)	<0.001	Fisher

Table legend: Statistically significant p-values are marked in bold italic font. SD: Standard Deviation, IQR: Interquartile Range, BMI: Body Mass Index, TIA: Transient Ischemic Attack, NMSS: Non-Motor Symptoms Scale, IBS: Irritable Bowel Syndrome.

Contrasting stable and progressed patients (Table 4), defined based on between-timepoint change in Unified Parkinson's Disease Rating Scale (UPDRS) I-III sum and total medication load calculated using the Levodopa Equivalent Dose (LED), there were no differences in basic demographics. The UPDRS I-III sum was stable or even decreased slightly between timepoints in stable patients (p = 0·004) and increased in progressors (p = 0·004). While there was no statistically significant difference in UPDRS I-III sum between groups at baseline (p = 0·201), there was at follow-up (p = 0·020). LED was similar for both groups at both timepoints and increased in both between timepoints (p (stable) < 0·001; p (progressors) = 0·001). Medications that differed between progression groups at baseline were acetylsalicylic acid (p = 0·037), statins (p = 0·010) and ropinirole (p (dose (mg)) = 0·048). At follow-up, progressors used more COMT inhibitors (p (entacapone (mg)) = 0·027; p (yes/no) = 0·051) and less pramipexole (p (dose (mg)) = 0·001).

Table 4

Clinical variables in stable and progressed Parkinson's patients at each timepoint.

Variable	Timepoint	Stable patients (n = 41)		Progressed patients (n = 15)
		% (n) / median [IQR]	p-value (baseline vs follow-up)	% (n) / median [IQR]	p-value (baseline vs follow-up)	p-value (progressed vs stable)	Test
Sex (% males)		53.66 (22)		46.67 (7)		0.765	Fisher
Age at stool collection		65 [61–69]		65 [62.5–66.5]		0.802	Wilcox
Age PD diagnosed		60 [58–64]		60 [57–65]		0.963	Wilcox
BMI	Baseline	26.25 [24.15–29.8]	0.225	26.53 [25.28–28.24]	0.890	0.923	Wilcox
	Follow-up	27.28 [24.27–30.21]		26.9 [23.89–29.89]		0.608
UPDRS I to III score total	Baseline	45 [35–52]	0.004	36 [33–46.5]	0.004	0.201	Wilcox
	Follow-up	38 [32–49]		48 [41–49]		0.020
UPDRS IV total	Baseline	2 [0–3]	0.062	1 [0.5–2]	0.026	0.685	Wilcox
	Follow-up	1 [0–5]		1 [1–6]		0.277
Hoehn & Yahr stage	Baseline	2.5 [2–2.5]	0.735	2.5 [2–2.5]	0.021	0.453	Wilcox
	Follow-up	2.5 [2–3]		2.5 [2.5–3]		0.038
Rome III IBS criteria fulfilled	Baseline	24.39 (10)	0.467	6.67 (1)	0.330	0.255	Fisher
	Follow-up	34.15 (14)		26.67 (4)		0.751
Jankovic subtypes
	Baseline	MX: 9.76 (4), PIGD: 53.66 (22), TD: 36.59 (15)	0.412	MX: 20.00 (3), PIGD: 40.00 (6), TD: 40.00 (6)	0.029	0.481	Fisher
	Follow-up	MX: 2.44 (1), PIGD: 53.66 (22), TD: 43.9 (18)		MX: 0 (0), PIGD: 86.67 (13), TD: 13.33 (2)		0.070
Medications
LED (mg)	Baseline	420 [205–505]	<0.001	340 [210–585]	0.001	0.774	Wilcox
	Follow-up	505 [362–662]		604 [440–811.75]		0.177
Levodopa (%)	Baseline	53.66 (22)	0.171	46.67 (7)	0.128	0.765	Fisher
	Follow-up	70.73 (29)		80 (12)		0.735
Levodopa (mg)	Baseline	100 [0–400]	<0.001	0 [0−300]	0.034	0.715	Wilcox
	Follow-up	300 [0–500]		350 [100–475]		0.708
COMT inhibitors (%)	Baseline	12.2 (5)	1.000	6.67 (1)	0.080	1.000	Fisher
	Follow-up	12.2 (5)		40 (6)		0.051
Entacapone (mg)	Baseline	0 [0–0]	0.850	0 [0–0]	0.204	0.617	Wilcox
	Follow-up	0 [0–0]		0 [0–750]		0.027
Acetylsalicylic acid (%)	Baseline	17.07 (7)	1.000	46.67 (7)	0.450	0.037	Fisher
	Follow-up	17.07 (7)		26.67 (4)		0.461
Statins (%)	Baseline	12.2 (5)	1.000	46.67 (7)	0.710	0.010	Fisher
	Follow-up	14.63 (6)		33.33 (5)		0.142
Pramipexole (mg)	Baseline	0.26 [0–1.05]	0.169	0 [0–0.52]	0.098	0.160	Wilcox
	Follow-up	0.26 [0–1.57]		0 [0–0]		0.001
Ropinirole (mg)	Baseline	0 [0–0]	0.887	8 [0–8]	0.904	0.048	Wilcox
	Follow-up	0 [0–0]		0 [0–8]		0.062

Table legend: Statistically significant p-values are marked in bold italic font. IQR: Interquartile Range, BMI: Body Mass Index, UPDRS: Unified Parkinson's Disease Rating Scale, IBS: Irritable Bowel Syndrome, TD: tremor dominant, PIGD: postural instability and gait difficulty, MX: mixed disease phenotype, LED: Levodopa Equivalent Dose, COMT: catechol-O-methyl transferase.

Regarding dietary data, there were no significant differences in intakes of any dietary items between patients and controls, or stable and progressed patients (Supplementary results). However, the first principal component (PC1) of a Principal Component Analysis (PCA) seemed to reflect diet healthiness (Table S1), with PD patients more often on the unhealthy side (Fig. 3). We kept PC1 as a potential confounder to be assessed in further analyses.

Fig. 3

Biplot of principal component analysis of diet data, showing principal components 1 and 3.

Legend: A. Component loadings; B. PD patients and controls, with 90% confidence ellipses.

Microbiota data

The 16S rRNA gene amplicon data contained 2836 OTUs, 198 genera and 77 families. The most common taxa were similar for the PD and control groups at both timepoints, with Ruminococcaceae and Lachnospiraceae dominating both groups' microbiota (Fig. 4A). Subsetting PD patients by progression status did not suggest major differences in microbial communities between these groups (Fig. 4B).

Fig. 4

10 most common bacterial families at each timepoint.

Legend: A. All study subjects by PD status; B. PD patients by progression status.

Considering commonly used ratios of specific bacteria, the Firmicutes/Bacteroidetes ratio did not differ between patients and controls, but Prevotella/Bacteroides was higher in controls (p (baseline) = 0·052, p (follow-up) = 0·011). Comparing progressed and stable patients, Firmicutes/Bacteroidetes differed significantly at baseline (p = 0·012) but not at follow-up, while Prevotella/Bacteroides did not differ between groups.

Enterotype analysis suggested that PD patients in general, and progressed patients in particular (Fig. 5), were overrepresented in the Firmicutes-dominated enterotype. This difference was statistically significant at both timepoints when contrasting controls and all patients (p (baseline) = 0·044; p (follow-up) = 0·025) or controls and progressed patients (p (baseline) < 0·001; p (follow-up) = 0·043), and at baseline when contrasting stable and progressed patients (p (baseline) = 0·016; p (follow-up) = 0·291). Enterotype distribution did not differ significantly between stable patients and controls (p (both timepoints) > 0·3).

Fig. 5

Distributions of enterotypes in the subject groups.

Legend: Showing control subjects and PD patients classified by progression (excluding PD patients without a progression classification).

Alpha diversity

We explored alpha diversity (microbial community richness and evenness) using three different indices (observed richness, Shannon, inverse Simpson). There was no difference between timepoints (p > 0·1, all indices), controls and patients (p > 0·6, all indices, both timepoints), stable and progressed patients (p > 0·2, all indices, both timepoints), or PD phenotypes (postural instability and gait difficulty (PIGD) vs tremor dominant (TD); p > 0·3, all indices, both timepoints). To look for effects of potential confounders, we calculated their correlations with each index (Table 5). The only variable with significant adjusted p-values in these comparisons was CCB use (adjusted p (Shannon) = 0·042; adjusted p (inverse Simpson) = 0·039). BMI and history of ear, nose and throat (ENT) surgery had a significant unadjusted p for all indices. We also found interactions between PD status and BMI in linear regression for the Shannon and inverse Simpson indices, and between PD status and average Victoria Bowel Performance Scale (BPS) score (follow-up only) when modelling observed richness (Fig. 6, Supplementary results, Table S2). Dietary variables were not associated with differences in alpha diversity (Supplementary results).

Table 5

Correlations for alpha diversity measures with variables of interest and confounders.

	Observed richness			Shannon			Inverse Simpson
	Pearson correlation	p-value	Adjusted p-value	Pearson correlation	p-value	Adjusted p-value	Pearson correlation	p-value	Adjusted p-value
Timepoint	0.104	0.098	0.496	0.042	0.506	0.832	0.055	0.382	0.873
PD status	0.005	0.942	0.990	−0.023	0.718	0.908	−0.026	0.678	0.915
Meds CCB	0.136	0.030	0.266	0.216	0.001	0.042	0.217	<0.001	0.039
History ENT surgery	0.169	0.007	0.184	0.189	0.002	0.096	0.163	0.009	0.181
BMI	−0.200	0.001	0.057	−0.159	0.012	0.190	−0.129	0.041	0.417
History CAD	0.105	0.096	0.496	0.168	0.007	0.153	0.169	0.007	0.181
Tobacco 100 in life	0.092	0.140	0.630	0.167	0.008	0.153	0.140	0.025	0.403
Selegeline mg	−0.043	0.494	0.970	−0.124	0.048	0.487	−0.163	0.009	0.181
History IBS	0.149	0.017	0.232	0.116	0.064	0.517	0.131	0.036	0.417
History hyper-thyreosis	−0.034	0.591	0.970	−0.143	0.022	0.299	−0.110	0.078	0.576
History appendectomy	−0.054	0.392	0.970	−0.133	0.034	0.388	−0.110	0.078	0.576
History lactose intolerance	0.129	0.039	0.280	0.116	0.063	0.517	0.108	0.086	0.579
History hernia repair	0.013	0.835	0.990	0.101	0.108	0.729	0.135	0.031	0.417
History hypothyreosis	0.161	0.010	0.199	0.078	0.212	0.832	0.031	0.618	0.894
Wexner total	0.151	0.015	0.232	0.036	0.569	0.832	0.020	0.744	0.972
Rome III 9–15 sum score	0.212	0.001	0.052	0.031	0.616	0.832	0.010	0.872	0.985
Rome III IBS criteria fulfilled	0.128	0.040	0.280	0.033	0.597	0.832	−0.007	0.913	0.985
Meds thyroxine	0.143	0.022	0.258	0.051	0.412	0.832	0.010	0.874	0.985
Food PC1	−0.127	0.042	0.280	−0.087	0.164	0.832	−0.055	0.379	0.873
Pramipexole mg	−0.139	0.026	0.265	−0.028	0.657	0.845	−0.001	0.985	0.985

Table legend: Showing only variables with a significant p-value in at least one comparison, and the PD status and Timepoint variables. Statistically significant p-values are marked in bold italic font. CCB: Calcium Channel Blockers, ENT: Ear, Nose and Throat, BMI: Body Mass Index, CAD: Coronary Artery Disease, IBS: Irritable Bowel Syndrome, PC: Principal Component.

Fig. 6

Interactions in linear regression of alpha diversity, PD status and confounders.

Legend: A. Shannon index and BMI; B. Inverse Simpson index and BMI, C. Observed richness and Victoria Bowel Performance Scale.

Beta diversity

Regarding beta diversity (between-sample community dissimilarity), subjects' microbial communities did not differ between timepoints in any comparison, but they did between patients and controls (p (PD status) < 0·017, all comparisons; Table 6A–C, Fig. 7). The confounding variables with the most consistent effects were constipation scores (Rome III: p ≤ 0·001, all comparisons; Wexner: p ≤ 0·004 in all single-confounder comparisons, dropped from further comparisons due to collinearity) and BMI (p < 0·025 in all comparisons except envfit test at genus and family level); two medication variables, CCB and ACE-I/ARB (angiotensin-converting-enzyme inhibitor / angiotensin II receptor blocker), and PC1 from the diet analysis were also significant in more than one comparison (Table 6B-C, Fig. 8). Separate beta diversity comparisons for dietary variables revealed no strong dietary confounders (Table S3, Supplementary results).

Table 6

Beta diversity comparisons of PD patients, control subjects and confounding variables.

A. Adonis: Timepoint + PD status (no confounders)
Level	p-value (Timepoint)	p-value (PD status)
OTU	0.898	<0.001
Genus	0.544	<0.001
Family	0.647	<0.001

Table legend: Statistically significant p-values are marked in bold italic font. OTU: Operational Taxonomic Unit, BMI: Body Mass Index, ACE-I/ARB: angiotensin-converting-enzyme inhibitor / angiotensin II receptor blocker, PC: Principal Component, CCB: Calcium Channel Blocker, RLS: Restless Legs Syndrome, TIA: Transient Ischemic Attack, MMSE: Mini-Mental State Examination, IBS: Irritable Bowel Syndrome.

Fig. 7

NMDS ordination illustrating microbial community differences between PD patients and control subjects.

Legend: Also showing centroid locations for microbial genera reported as differentially abundant between the two groups in previous studies.

Fig. 8

Genus-level NMDS ordination plots of the main confounding variables.

Legend: A. BMI, diet and Rome III score; B. Calcium channel blockers (CCB) medication; C. ACE-I/ARB medication.

Comparing progressed and stable patients, there were no beta diversity differences between timepoints nor for progression (Table 7A–C), while COMT inhibitor use had a very significant community effect (p ≤ 0·001, all comparisons except the envfit test on genus and family levels; Table 7B–C). Other confounders which were significant in more than one comparison were statins, SCS-PD total, and NMSQuest. There were no differences for beta diversity between PD phenotypes (Supplementary results).

Table 7

Beta diversity comparisons of stable and progressed PD patients and confounding variables.

A. Adonis: Timepoint + progression category (no confounders)
Level	p-value (Timepoint)	p-value (progression)
OTU	1.000	0.245
Genus	0.945	0.344
Family	0.850	0.142

Table legend: Statistically significant p-values are marked in bold italic font. OTU: Operational Taxonomic Unit, SDQ: Swallowing Disturbance Questionnaire, COMT: catechol-O-methyl transferase, NMSQuest: Non-Motor Symptoms Questionnaire, SCS-PD: Sialorrhea Clinical Scale for PD, GDS 15: 15-item Geriatric Depression Scale, NMSS: Non-Motor Symptoms Scale, RBDSQ: REM Sleep Behavior Disorder Screening Questionnaire, MAO: Monoamine Oxidase, PIGD: postural instability and gait difficulty.

Differential abundance of microbial taxa

We used three differential abundance comparison methods, all of which suggested several taxa as differing significantly between PD patients and controls (ANCOM: 2 families, 1 genus and 6 OTUs for baseline and 3 families, 3 genera and 8 OTU for follow-up; random forests: 4 families, 5 genera and 33 OTUs for baseline, 3 families, 9 genera and 29 OTUs for follow-up; DESeq2: 3 families, 5 genera and 2 OTUs at baseline, 3 families and 5 genera at follow-up; Fig. 9, Supplementary material: Table S4). All random forest classifiers were significantly better than chance (p ≤ 0·046 for all comparisons; Table S5). Combining the results of the three methods, a handful of taxa overlapped at one or both timepoints, mainly the families Bifidobacteriaceae, Prevotellaceae, and Puniceicoccaceae, and the genera Bifidobacterium, Roseburia, Prevotella, and Clostridium XIVa (Table 8, Fig. 9, Fig. 10). A few other taxa reported as differentially abundant in previous studies were significant only according to random forests: the family Lachnospiraceae and the genus Blautia at baseline, and Lactobacillaceae and Lactobacillus at follow-up. Considering the confounders included in DESeq2, 2 OTUs, the genera Clostridium and Coprococcus, and the family Clostridiaceae 1 were associated with BMI, and Bifidobacterium and Bifidobacteriaceae with the Rome III score (Supplementary material: Table S4C). Bifidobacteriaceae were not associated with PD status according to DESeq2 at either timepoint (Table 8). Finally, the DESeq2 model enabled the detection of microbes with a difference in between-timepoint trends when comparing controls and PD patients, but we found no such taxa.

Fig. 9

Differentially abundant taxa between control subjects and PD patients according to three different tools.

Legend: A. Baseline; B. Follow-up.

Table 8

Summary of differential abundance results contrasting PD patients and control subjects.

			Baseline					Follow-up
	Level	Previously published	Relative abundance (%, mean ± SD)		p-value			Relative abundance (%, mean ± SD)		p-value
			Control	PD	ANCOM	Random Forests	DESeq2	Control	PD	AN-COM	Random Forests	DESeq2
Bifidobacteriaceae	Family	Yes	2.629 ± 2.719	6.528 ± 8.573	s.s.	0.129	0.151	2.189 ± 3.531	5.919 ± 8.256	s.s.	0.040	0.150
Prevotellaceae	Family	Yes	4.345 ± 9.22	0.725 ± 2.12	s.s.	0.010	0.002	2.972 ± 4.882	1.395 ± 3.474	s.s.	0.079	0.002
Rikenellaceae	Family	Yes	2.201 ± 2.035	3.494 ± 2.899	n.s.	0.030	0.075	2.479 ± 2.559	2.836 ± 2.131	n.s.	0.386	0.172
Lachnospiraceae	Family	Yes	22.478 ± 10.241	16.48 ± 9.121	n.s.	0.020	0.990	21.787 ± 8.964	17.753 ± 8.198	n.s.	0.416	0.995
Pasteurellaceae	Family	Yes	0.036 ± 0.097	0.009 ± 0.027	n.s.	0.089	0.625	0.059 ± 0.28	0.014 ± 0.052	n.s.	0.158	0.585
Lactobacillaceae	Family	Yes	0.032 ± 0.124	0.28 ± 1.199	n.s.	0.129	0.993	0.04 ± 0.159	0.226 ± 0.867	n.s.	0.040	0.995
Puniceicoccaceae	Family	No	0.045 ± 0.103	0.01 ± 0.029	n.s.	0.010	0.993	0.032 ± 0.061	0.01 ± 0.03	s.s.	0.010	0.995
Bifidobacterium	Genus	Yes	2.628 ± 2.717	6.524 ± 8.57	s.s.	0.089	0.244	2.188 ± 3.53	5.917 ± 8.256	s.s.	0.040	0.410
Roseburia	Genus	Yes	7.014 ± 6.944	3.588 ± 4.096	n.s.	0.020	0.147	6.683 ± 6.162	4.395 ± 5.298	s.s.	0.010	0.121
Prevotella	Genus	Yes	4.187 ± 9.003	0.659 ± 2.102	n.s.	0.050	0.042	2.85 ± 4.854	1.306 ± 3.355	s.s.	0.030	0.043
Anaerotruncus	Genus	Yes	0.056 ± 0.079	0.071 ± 0.084	n.s.	0.079	0.986	0.119 ± 0.501	0.098 ± 0.163	n.s.	0.693	0.990
Blautia	Genus	Yes	2.419 ± 1.516	1.829 ± 1.598	n.s.	0.040	0.742	2.63 ± 1.847	2.429 ± 2.334	n.s.	0.644	0.927
Lactobacillus	Genus	Yes	0.031 ± 0.125	0.278 ± 1.197	n.s.	0.356	0.989	0.04 ± 0.159	0.221 ± 0.862	n.s.	0.030	0.996
Clostridium XlVa	Genus	No	1.971 ± 2.631	1.83 ± 3.265	n.s.	0.168	<0.001	2.136 ± 2.589	1.501 ± 2.076	n.s.	0.020	<0.001
Otu0003 Roseburia	OTU	Yes	6.506 ± 6.71	3.115 ± 3.742	n.s.	0.020	0.342	6.144 ± 5.933	4.109 ± 5.027	n.s.	0.020	0.275
Otu0007 Bifidobacterium	OTU	Yes	1.563 ± 1.808	4.579 ± 6.456	s.s.	0.079	0.613	1.518 ± 3.155	4.029 ± 6.243	s.s.	0.089	0.692
Otu0024 Blautia	OTU	Yes	0.904 ± 0.875	0.649 ± 0.903	n.s.	0.020	1.000	1.057 ± 1.14	0.712 ± 0.884	n.s.	0.020	1.000
Otu0027 Ruminococcus	OTU	Yes	0.679 ± 0.959	0.437 ± 0.895	n.s.	0.089	1.000	0.832 ± 1.515	0.604 ± 1.23	n.s.	0.960	1.000
Otu0030 Alistipes	OTU	Yes	0.404 ± 0.861	0.955 ± 1.392	s.s.	0.010	1.000	0.711 ± 1.627	0.623 ± 1.134	n.s.	0.535	1.000
Otu0036 Roseburia	OTU	Yes	0.483 ± 0.757	0.468 ± 1.279	n.s.	0.089	1.000	0.534 ± 0.771	0.282 ± 0.689	n.s.	0.030	1.000
Otu0041 Alistipes	OTU	Yes	0.278 ± 0.247	0.448 ± 0.499	n.s.	0.099	1.000	0.33 ± 0.342	0.412 ± 0.37	n.s.	0.713	1.000
Otu0055 Prevotella	OTU	Yes	0.344 ± 1.243	0.197 ± 1.048	n.s.	0.129	1.000	0.493 ± 1.758	0.342 ± 1.413	n.s.	0.079	1.000
Otu0062 Blautia	OTU	Yes	0.446 ± 0.576	0.266 ± 0.358	n.s.	0.030	0.049	0.328 ± 0.46	0.217 ± 0.299	n.s.	0.198	0.111
Otu0098 Bacteroides	OTU	Yes	0.105 ± 0.239	0.29 ± 0.599	n.s.	0.030	0.499	0.111 ± 0.22	0.231 ± 0.407	n.s.	0.079	0.371
Otu0109 Ruminococcus	OTU	Yes	0.343 ± 1.142	0.122 ± 0.844	n.s.	0.465	1.000	0.369 ± 1.175	0.001 ± 0.002	s.s.	0.010	1.000
Otu0110 Ruminococcus	OTU	Yes	0.189 ± 0.332	0.136 ± 0.279	n.s.	0.020	1.000	0.205 ± 0.519	0.14 ± 0.322	n.s.	0.139	1.000
Otu0131 Bacteroides	OTU	Yes	0.183 ± 0.654	0.079 ± 0.256	n.s.	0.406	0.998	0.209 ± 0.613	0.049 ± 0.205	s.s.	0.020	0.995
Otu0363 Lactobacillus	OTU	Yes	0.019 ± 0.107	0.056 ± 0.225	n.s.	0.792	1.000	0.016 ± 0.122	0.006 ± 0.035	n.s.	0.050	1.000
Otu0379 Alistipes	OTU	Yes	0.001 ± 0.003	0.041 ± 0.2	n.s.	0.317	1.000	0.004 ± 0.033	0.047 ± 0.236	n.s.	0.020	1.000
Otu0464 Lactobacillus	OTU	Yes	0.002 ± 0.009	0.004 ± 0.016	n.s.	0.446	1.000	0.001 ± 0.006	0.044 ± 0.236	n.s.	0.010	1.000
Otu0468 Faecalibacterium	OTU	Yes	0.016 ± 0.026	0.009 ± 0.019	n.s.	0.723	1.000	0.021 ± 0.033	0.007 ± 0.016	n.s.	0.050	1.000
Otu0513 Anaerotruncus	OTU	Yes	0.006 ± 0.006	0.011 ± 0.011	n.s.	0.030	0.998	0.006 ± 0.005	0.012 ± 0.019	n.s.	0.208	0.995

Table legend: This table includes (1) taxa that are differentially abundant according to more than one method at either timepoint, and (2) taxa that have been reported as differentially abundant in previous studies and have p < 0.1 with at least one method at either timepoint in this study. Statistically significant p-values are marked in bold italic font. OTU: Operational Taxonomic Unit, SD: Standard Deviation, s.s.: statistically significant, n.s.: not significant.

Fig. 10

Genera and families that differ significantly between control subjects and PD patients.

Legend: Showing taxa that differ at either timepoint according to more than one method, or have been reported as differing in a previous publication and have p < 0.1 at one timepoint according to at least one method.

Analyses contrasting progressed patients to stable ones also resulted in lists of several differentially abundant taxa (Supplementary material: Table S6A—C). However, random forest models failed to differentiate progressors and stable patients better than chance; the only comparison with a significant p-value was OTU level at follow-up (p = 0·044), and even there, the difference between actual and randomly permuted error values was small (Table S7). The differential abundance results of the three methods did not overlap substantially. One low-abundance, unclassified Lachnospiraceae OTU was significant according to all methods at follow-up; DESeq2 and random forests agreed on one Streptococcus OTU and the family Streptococcaceae at baseline, and three OTUs, the genus Asteroleplasma, and the family Anaeroplasmataceae at follow-up; and a single Bifidobacterium OTU, which was more common in progressors at follow-up, was supported by both ANCOM and random forests (Table 9). The genus Prevotella was more abundant in stable patients at both timepoints according to DESeq2, while the family Prevotellaceae only differed at follow-up (Fig. 11; Table 9; Supplementary results: Table S6C). There were as many or more taxa associated with COMT inhibitors in the DESeq2 results as there were with progression, although these were also inconsistent across timepoints (Supplementary results: Table S6C). Among them, Bifidobacterium was more common in COMT users than non-users at follow-up.

Table 9

Summary of differential abundance results contrasting progressed and stable PD patients.

		Baseline					Follow-up
		Relative abundance (%, mean ± SD)		p-value			Relative abundance (%, mean ± SD)		p-value
	Level	Stable	Progressed	AN-COM	Random Forests	DESeq2	Stable	Progressed	AN-COM	Random Forests	DESeq2
Streptococcaceae	Family	0.207 ± 0.584	0.614 ± 1.316	n.s.	0.010	0.001	0.299 ± 0.904	0.398 ± 0.74	n.s.	0.663	0.176
Anaeroplasmataceae	Family	0.426 ± 1.792	0.371 ± 1.125	n.s.	0.198	0.966	0.133 ± 0.596	0.311 ± 1.05	n.s.	0.010	<0.001
Prevotella	Genus	0.964 ± 2.571	0.147 ± 0.493	n.s.	0.554	<0.001	1.966 ± 4.04	0.187 ± 0.678	n.s.	0.386	<0.001
Asteroleplasma	Genus	0.191 ± 1.022	0.302 ± 1.113	n.s.	0.455	0.767	0.133 ± 0.596	0.31 ± 1.051	n.s.	0.050	<0.001
Otu0148 Bifidobacterium	OTU	0.06 ± 0.252	0.143 ± 0.498	n.s.	0.465	0.862	0.051 ± 0.125	0.744 ± 1.324	s.s.	0.010	0.330
Otu0327 Lachnospiraceae unclassified	OTU	0.013 ± 0.032	0.104 ± 0.282	n.s.	0.050	0.965	0.021 ± 0.091	0.178 ± 0.357	s.s.	0.010	0.016
Otu0118 Ruminococcaceae unclassified	OTU	0.084 ± 0.102	0.246 ± 0.221	n.s.	0.010	0.094	0.135 ± 0.261	0.165 ± 0.14	n.s.	0.040	0.838
Otu0166 Ruminococcaceae unclassified	OTU	0.069 ± 0.102	0.166 ± 0.203	n.s.	0.020	0.531	0.082 ± 0.097	0.194 ± 0.227	n.s.	0.040	0.603
Otu0222 Phascolarctobacterium	OTU	0.091 ± 0.281	0 ± 0	n.s.	0.436	0.002	0.043 ± 0.137	0 ± 0	n.s.	0.911	0.000
Otu0111 Streptococcus	OTU	0.127 ± 0.369	0.46 ± 1.187	n.s.	0.050	0.007	0.217 ± 0.643	0.254 ± 0.501	n.s.	0.743	0.438
Otu0042 Coprococcus	OTU	0.344 ± 1.26	0.002 ± 0.002	n.s.	0.733	<0.001	0.511 ± 1.44	0.004 ± 0.011	n.s.	0.525	0.000
Otu0268 Desulfovibrio	OTU	0.059 ± 0.168	0.034 ± 0.131	n.s.	0.475	<0.001	0.042 ± 0.102	0.001 ± 0.002	n.s.	0.792	0.009
Otu0115 Lachnospiraceae unclassified	OTU	0.167 ± 0.404	0.607 ± 2.296	n.s.	0.921	0.018	0.04 ± 0.086	0.776 ± 2.121	n.s.	0.010	0.014
Otu0049 Ruminococcaceae unclassified	OTU	0.26 ± 0.707	0.595 ± 0.919	n.s.	0.040	0.864	0.286 ± 0.663	0.56 ± 0.896	n.s.	0.030	0.297
Otu0241 Clostridiales unclassified	OTU	0.029 ± 0.054	0.044 ± 0.064	n.s.	0.446	0.698	0.015 ± 0.022	0.084 ± 0.159	n.s.	0.040	0.020
Otu0084 Clostridium IV	OTU	0.313 ± 0.978	0.335 ± 0.726	n.s.	0.485	0.864	0.071 ± 0.091	0.275 ± 0.345	n.s.	0.010	0.001

Table legend: This table includes (1) bacterial taxa that are differentially abundant at one timepoint according to more than one method, or (2) at both timepoints according to one method. Statistically significant p-values are marked in bold italic font. OTU: Operational Taxonomic Unit, SD: Standard Deviation, s.s.: statistically significant, n.s.: not significant.

Fig. 11

Genera and families that differ significantly between stable and progressed PD patients.

Legend: Showing taxa that differ at one timepoint according to more than one method, or at both timepoints according to a single method.

When contrasting PD phenotypes (TD and PIGD), ANCOM detected no differentially abundant families or genera, and only 2 OTUs at baseline, while DESeq2 detected longer lists of taxa (Supplementary results; Table S8).

Discussion

Several studies have demonstrated that in healthy subjects, intersubject variability is greater than temporal variability in individuals [47,48]. Based on our results, this seems to also hold true for PD patients. In contrast to the lack of differences between timepoints, the microbial communities clearly differed between PD patients and controls, which is in line with previous studies of PD and gut microbiota: they have consistently reported differences in beta diversity between PD patients and controls [ [5], [6], [7],[10], [11], [12], [13], [14], [15], [16]].

The results of differential abundance comparison methods based on different statistical approaches can vary widely [49], offering one explanation for the varying results of the specific taxa reported in previous publications on PD and gut microbiota [[5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16]]. Other potential sources of variation include differences in the subjects' demographics and clinical details, technical differences in the sampling, storage and sequencing protocols, and different sample sizes. Our study ranks among the larger publications when it comes to the number of subjects, but it is still possible that we have missed some subtle effects due to not having enough samples.

In the present study, we used three different tools to look for taxa that differ between PD patients and controls, and our results underlined the variation between methods: the lists of detected taxa overlapping between methods were short. The only taxa significant at both timepoints according to multiple methods were Prevotella and Prevotellaceae [6,9,12,13]. Some taxa reported to differ between PD patients and control subjects by earlier studies, such as Akkermansia [5,7,11,12] and Verrucomicrobiaceae [[5], [6], [7],11,13,15], did not differ significantly in our comparisons, while many others, including Bifidobacterium and Bifidobacteriaceae [9,11,13], Lactobacillus and Lactobacillaceae [6,11,13], Lachnospiraceae [7,11,15] and Roseburia [7,11], were inconsistent across tools and timepoints. An additional, new taxon that was detected as significant at both timepoints in our study was Puniceicoccaeae, but as this family mainly contains environmental species [50] and was present at a very low abundance (mean ± SD (%): 0·025 ± 0·065), it seems unlikely to be of interest.

The baseline samples compared in this analysis were a subset of our earlier study [6], but the new results differed from our previous ones. For example, the new analyses did not support a difference in Lactobacillaceae at baseline. As the earlier study used different PCR primers, sequencing platform, and statistical tools, it is not unexpected that the results also differ. Additionally, the baseline stool samples were stored longer and thawed twice, compared to the follow-up samples which were stored for a briefer period and only thawed once. Although we could detect no overall differences between timepoints in alpha or beta diversity comparisons, some of the variation in lists of differentially abundant taxa could relate to this.

Regarding disease progression, the differentially abundant taxa supported by more than one method did not overlap between timepoints aside from two unclassified Lachnospiraceae OTUs. According to DESeq2, Prevotella was more abundant in stable subjects at both timepoints, and Prevotellaceae at follow-up. We also contrasted patients with different disease phenotypes (TD and PIGD). Although the difference in Enterobacteriaceae between phenotypes detected in our pilot study with the metastats method [6] was significant in the resequenced baseline data when tested with DESeq2, there was no significant difference for this family at follow-up. The family Anaeroplasmataceae and the genus Asteroleplasma were significant in both progression and PD phenotype comparisons, but both taxa had values close to zero in nearly all samples, so these are likely to be false positives caused by the outlier subjects.

One previous publication on gut microbiota and PD progression, also with a 2-year follow-up period, found that stable and deteriorated patients had similar microbiota at each timepoint, although Bifidobacterium was less abundant in deteriorated patients at baseline [17]. Several taxa, including Prevotella, decreased between timepoints, and the abundances of Bifidobacterium and B. fragilis were negatively correlated to UPDRS I [17]. In our study, we saw no overall decreasing trend for Prevotella, although the genus was more abundant in stable subjects. Bifidobacterium was more abundant in COMT inhibitor users at follow-up. The previous study assessed only selected microbial taxa and did not consider the changes in medication over time in their progression classification [17]. To get a more universal measure of PD symptom severity than pure motor function, and to reduce dependency on the rater, we chose to use the sum of UPDRS I-III as a symptom severity measure for progression calculations [51]. For practical reasons, UPDRS-III assessments in medication OFF-state were only available for the follow-up timepoint of our cohort. Therefore, we chose to use UPDRS-III in medication ON-state for our calculations [51]. Thus, UPDRS scores are to some degree influenced by PD medication load. A dichotomous classification into progressed and stable, only based on this score, could falsely classify patients with progressed disease as stable e.g. if they had been undertreated at baseline, and this was compensated during the follow-up interval by overproportional medication increase. This may have influenced the results of the previous longitudinal study on microbiota in PD [17]. If disease severity is assessed in ON-state, the progression classifier must at the same time account for the PD medication load that may make motor symptom severity appear less than it actually is. The severity measure that we used to classify our patients by progression is driven by symptom severity and medication load and thus increase in either one will increase the probability of a patient being classified as a progressor. Even though not commonly used previously, our measure is based on a widely used indicator of disease severity (UPDRS I-III in ON-state) and is adjusted for medication load, which we deem essential in this context.

Bearing in mind the methodological differences, the discordant results of the two follow-up studies are unsurprising. Taken together, they seem to imply that there is no strong microbial signal associated with PD progression. However, measuring progression is challenging. A 2-year follow-up period is short, and the number of patients classified as progressed in our comparisons was low (15 patients). A longer follow-up period and a larger patient cohort might be required to observe progression-related changes in microbiota.

The lowered abundance of the family Prevotellaceae in PD was a key finding of our pilot study [6]. Considering the results of our follow-up analysis, it remains a taxon of interest. Based on our FFQ comparisons, the abundance of this taxon was not strongly affected by dietary factors. The difference in Prevotella abundance also remained significant when correcting for Rome III constipation score, although this does not entirely exclude the possibility that the results are confounded by the constipation commonly seen in PD [52]. Prevotella is a common colonizer of the human gut, and the dominant taxon in one of the three suggested enterotypes [38]. It has been linked to numerous medical conditions, its role varying from beneficial to detrimental, possibly depending on the species or strain [53]. Two PD studies have reported a lowered abundance particularly for the species P. copri [[10], [12]]. Additionally, one study found the genus less abundant in PD patients with IBS-like symptoms [1]. We have no strain level information in the present study, as 16S rRNA gene amplicon sequencing lacks the resolution for reliably defining species or strains. Hopefully, future research into the features of different Prevotella taxa will shed more light on their potential relationship to PD.

Numerous factors, such as BMI, sex, diet, and medications, can influence gut microbiota [11,18,54]. Our subjects represent a small and fairly uniform population, and patients and controls were matched for age and sex, but some differences between groups are inevitable: for example, PD patients have more GI symptoms than the general population [1,2,52], and take PD-specific medications. Additionally, stroke or TIA were more common among our controls, who also used CCBs, statins and warfarin more often than patients. Diversity comparisons suggested that CCBs might influence gut microbiota diversity, but with so few subjects using them (baseline: 15, follow-up: 17), evaluating the significance of this effect is difficult. The same was true for COMT inhibitors: they were associated with differences in beta diversity and many differentially abundant taxa but were only used by a small subset of patients (baseline: 6, follow-up: 11). Nevertheless, as they have been previously linked to microbiota changes [6,11], future studies of PD and microbiota should consider COMT inhibitors as an important confounder.

Since diet could underlie some microbiota differences detected in earlier studies, we also collected FFQ data. PC1 from a dietary PCA implied a less healthy diet for patients and was associated with beta diversity differences in single-confounder comparisons, but in multiple-confounder models it was only significant on OTU level with one test, while other variables (PD status, Rome III score, BMI) remained significant on multiple taxonomic levels with two different tests. Overall, diet did not appear to be a major confounder. Dietary fibre, which is thought to influence gut microbiota (and specifically the abundance of Prevotella), was included in all diet comparisons, with separate variables for total fibre, insoluble fibre, and soluble fibre. There were no differences in fibre consumption between PD patients and controls, and fibre was not associated with significant differences in alpha or beta diversity. Another potentially important confounder, the use of probiotics, was equally insignificant in all comparisons. This is likely to reflect the fact that our data set is small, particularly considering the semiquantitative nature of FFQs, and may not have enough statistical power to catch the microbial influences of these diet variables.

Intriguing incidental findings were the interactions between PD and BMI in relation to alpha diversity indices that include richness and evenness, and PD and stool consistency for observed richness. An inverse correlation between BMI and alpha diversity has been reported in several previous publications [[55], [56], [57]]; in our analyses, this effect was noticeably weaker in PD patients. Additionally, the inverse correlation between the BPS stool consistency scale and richness, analogous to a previously reported inverse correlation between richness and the Bristol Stool Scale [58], was also stronger in control subjects than PD patients. These results offer support for PD patients having abnormal gut microbial communities that do not follow trends seen in analyses with healthy subjects. Statistical modelling that goes beyond linear regression between alpha diversity and a few key confounders could offer further information regarding the complex relationships between PD and the other variables associated with changes in alpha diversity.

Conclusions

In this study, we have shown that the differences in gut microbiota of PD patients and controls persist at follow-up sampling 2 years later. As all previous studies have only included data from a single timepoint, this is an important first confirmation of replicability of such findings. Our lists of differentially abundant taxa between patients and controls included many previously reported bacteria (such as Bifidobacterium, Prevotella, Lactobacillus, and Roseburia), although the results varied considerably between statistical tools. Progressed PD patients had a Firmicutes-dominated enterotype more often than stable patients or control subjects. Additionally, Prevotella, a genus already shown to be less abundant in PD patients compared to controls, also appears less abundant in patients with faster disease progression. This further underlines the potential importance of this genus in PD. A longer follow-up period might be warranted for capturing the trends in microbial community changes during disease progression.