Study population

The study sample was recruited from an ongoing project enrolling consecutive patients with early PD diagnosed according to the modified diagnostic criteria of the UK Parkinson’s Disease Society Brain Bank at the Movement Disorders Unit of the First Division of Neurology at the University of Campania “Luigi Vanvitelli” (Naples, Italy). In this study, drug-naïve PD patients with a modified Hoehn and Yahr (mH&Y) stage ≤ 2.5 were included. Exclusion criteria were: (1) PD onset before age 40 years; (2) any previous treatment with dopaminergic, anticholinergic, antidepressant, or other centrally acting drugs, to rule out a potential effect on functional connectivity from these agents; (3) relevant cognitive impairment; (4) indwelling catheter, renal dialysis, cardiac failure requiring diuretics, urinary tract infection, prostate or bladder cancer, uncontrolled bladder outlet obstruction, pelvic organ prolapse, and previous urogynecologic surgery; and (5) any other clinically significant medical condition.

Presence/absence of urinary symptoms were assessed by means of the Nonmotor Symptom Scale (NMSS) - domain 7. (Chaudhuri and Schapira 2009) Patients were considered to present urinary symptoms (PD-urinary+) if they scored ≥ 1. Consequently, patients with NMSS - domain 7 = 0 were labeled as PD-urinary−.

Thirty-eight healthy age and sex-matched controls (HCs) with no family history of PD or parkinsonism, and no history of urologic disorders or urinary symptoms were also included in the study. History of any neurological symptoms as well as medical and surgical significant comorbidities were also considered as exclusion criteria for controls.

All the subjects signed their written informed consent. The study was approved by the ethics committee of the University of Campania “Luigi Vanvitelli”, Naples, Italy.

Clinical evaluation - motor and nonmotor symptoms

Disease severity and motor performance were assessed using the mH&Y stages(Hoehn and Yahr 1967) and the Unified Parkinson’s Disease Rating Scale (UPDRS) part III (UPDRS-III)(Goetz et al. 2008). The presence and severity of other nonmotor symptoms were assessed by means of the NMSS.

Clinical evaluation – neuropsychological and behavioral symptoms

All patients underwent an extensive neurological and neuropsychological assessments as previously reported (De Micco et al. 2021b; Micco et al. 2021a) (see Supplementary information for more details).

Statistical analysis of clinical data

Demographic data between PD subgroups and HCs were compared using ANOVA models. T-tests were used to compare motor, nonmotor, neuropsychological, and behavioral variables between PD-urinary+ and PD-urinary−. Chi-square was used to assess differences in the distribution of categorical variables. Cohen’s d was used to explore the effect size for T-test. Analyses were all Bonferroni-corrected for multiple comparisons. Analyses were performed with SPSS version 23 (SPSS Inc., Chicago, IL, USA).

Imaging parameters

Magnetic resonance images were acquired on a General Electric 3 Tesla MRI scanner equipped with an 8-channel parallel head coil. fMRI data consisted of 240 volumes of a repeated gradient-echo echo planar imaging T2*-weighted sequence (TR = 1508 ms, axial slices = 29, matrix = 64 × 64, field of view = 256 mm, thickness = 4 mm, interslice gap = 0 mm). During the functional scan, subjects were asked to simply stay motionless, awake, with their eyes closed. Three-dimensional high-resolution T1-weighted sagittal images (GE sequence IR-FSPGR, TR = 6988 ms, TI = 1100 ms, TE = 3.9 ms, flip angle = 10, voxel size = 1 × 1 × 1.2 mm3) were acquired for registration and normalization of the functional images as well as for voxel-based morphometry (VBM) analysis.

FMRI preprocessing

fMRI data preprocessing was performed in Matlab® (The MathWorks, Inc., www.mathworks.com) with the toolbox Data Processing Assistant for Resting-State fMRI (DPARSF, Yan and Zang 2010, http://rfmri.org/DPARSF), which is based on Statistical Parametric Mapping (SPM, http://www.fil.ion.ucl.ac.uk/spm) and Data Processing & Analysis of Brain Imaging (DPABI, Yan et al. 2016, http://rfmri.org/DPABI). Each individual rs-fMRI time series was first corrected for the different slice scan acquisition times by specifying the number of slices, the slice acquisition order and the reference slice. The alignment of the first volume of each subject time-series to the corresponding anatomical 3D-T1w image was obtained via affine transformation; then, all T1w images from all subjects were normalized to the MNI space with the non-linear diffeomorphic DARTEL approach; e14last, the coregistered functional data were normalized to the MNI space with the transformations obtained during the DARTEL procedure and functional scans were resampled to 3 × 3 × 3 mm voxel sizes.

To reduce the residual effects of head motion, as well as the effects of respiratory and cardiac signals, second-order motion and physiological nuisance correction were performed using a linear regression approach: the regression model included 24 motion-related predictors (Friston et al. 1996), with 6 head motion parameter time-series, their first-order derivatives and the 12 corresponding squared parameter time-series; the mean time-courses from a white matter mask and a cerebrospinal fluid mask (as obtained from 3D-T1w spatial segmentation) were added as two additional predictors. In order to account for residual motion-related spikes, an additional spike-related regressor was created from the framewise displacement time-series, i.e. a predictor with a value of 1 at the time points of each detected spike and a value of 0 elsewhere (Lemieux et al. 2007; Satterthwaite et al. 2013). Finally, the image time series were band-pass filtered between 0.01 Hz and 0.5 Hz and spatially smoothed with an isotropic 6-mm full width at half maximum (FWHM) Gaussian kernel.

To minimize the potential effects of head motion and possibly exclude subjects exhibiting excessive amounts of motion, we applied severe inclusion criteria: the six estimated head motion parameters (3 translation and 3 rotation) were considered and subjects exhibiting head translations > 3 mm and/or head rotations > 3 degrees were excluded from consecutive analyses. Then, the mean framewise displacement value (FD) was estimated as an additional measure of total instantaneous head motion (Power et al. 2012) and the percentage of spike-corrupted volumes in each time-series was calculated. Potential spike-corrupted volumes were identified where the FD value exceeded a threshold of 0.5 mm; at this stage, subjects for whom the percentage of corrupted volumes exceeded 50% in the scan were also excluded from the analyses.

Rs-fMRI time series were imported in BrainVoyager and transformed to the Talairach space for seed-based analysis.

Seed-based connectivity analysis

A seed-based analysis was performed to study FC from striatal loop areas, including bilateral limbic, cognitive, and sensorimotor regions, which were defined according to previous works (Biondetti et al. 2021) to the entire brain. Six ROIs were defined from the anatomical brain atlas “striatum-con-label-thr50-3sub-1 mm” which is freely available for download from the FSL website (see https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases/striatumconn for all details). Two versions of the atlas exist: the first version comprises of 3 subdivisions (limbic, executive, and sensorimotor), whereas the second comprises of 7 subdivisions (limbic, executive, rostral-motor, caudal-motor, parietal, occipital, and temporal). More specifically we used the probabilistic connectivity striatal atlas with subdivision into three subregions according to cortical-striatal anatomical connections from right and left putamen, right and left caudate and right and left ventral striatum, resulting in six ROIs in total that were used as seed for the FC analysis.

To calculate the functional connectivity maps corresponding to each selected ROI, the mean regional time course was extracted from all ROI voxels and this was correlated against the time-courses from all voxel of the brain. Separate correlation (r) maps were produced for each subject of each group and ROI. Fisher’s transform z = 0.5 Ln [(1 + r)/(1– r)] was applied to these correlation maps before entering a second-level random-effects statistical analysis where the main and differential effects for the two studied groups were summarized as t-statistic maps. For these computations, we used an in-house written Matlab script which was based on an open source Matlab tool called Neuroelf which is freely available at https://neuroelf.org. This tool provides an easy and convenient interface from NIFTI data processed in any software to the commercial software BrainVoyager (Brain Innovation B. V., Maastricht, The Netherlands, www.brainvoyager.com) that we used here to perform the second-level voxel-based analysis and ultimately overlay the resulting statistical contrast (t-statistics) maps onto an high resolution standard T1 template image for display at the voxel- and cluster-level statistical thresholds as resulting from the correction for multiple comparisons. This analysis was carried out by treating the individual subject map values as random observations at each voxel, thereby the classical analysis of variance (ANOVA) was performed at each voxel to map the whole-brain distribution of the seed-based functional connectivity for the difference among groups. From this model, for each seed, statistical contrasts were derived for the following comparisons: PD-all vs. HC, PD-urinary+ vs. HCs, PD-urinary− vs. HCs, PD-urinary+ vs. PD-urinary−. An inclusive mask was created from a standard T1 template image, to define the search volume for multiple comparisons correction: this mask included the entire brain after excluding the striatal loop ROIs that were used as seed for the analysis. To correct for multiple comparisons in the voxel-based analysis, regional effects resulting from the voxel-based comparative tests in the search volume, were only accepted for compact clusters surviving the joint application of a voxel- and cluster-level threshold chosen with a nonparametric randomization approach. Namely, an initial voxel-level threshold was set to p = 0.001 (uncorrected) and a minimum cluster size was estimated after 1000 Monte Carlo simulations that protected against false positive cluster up to 5% (p = 0.05, cluster-level corrected) (Forman et al. 1995; Eklund et al. 2016; Anderkova et al. 2017).

Individual FC z-scores for both the patients’ subgroups were extracted from regions identified in the above analyses. A univariate analysis of variance was performed between the individual FC z-scores in PD-urinary+ and PD-urinary−, running NMSS items scores and use of medications for urinary disturbances as covariates.

Partial correlation analysis

Partial correlation coefficients were computed between imaging (i.e., FC z-scores) and clinical data (i.e., NMSS - domain 7 scores) in PD-urinary+. Analyses were adjusted for age and sex. A p < 0.05 was considered statistically significant. Analyses were performed with SPSS version 23 (SPSS Inc. Chicago, IL).

VBM analysis (see Supplementary information).