Participants were recruited from the outpatient clinic for movement disorders and were diagnosed with adult-onset ICD following the criteria (3, 22) by two senior neurologists. Patients who matched the following criteria were included in the study: (1) isolated cervical dystonia; (2) no history of botulinum toxin treatment, related drug therapy, or operation in the previous 3 months; (3) no history of serious physical or neuropsychiatric disorders; and (4) right-handedness. The patients were excluded according the following criteria: (1) age <18 years; (2) genetic findings associated with the onset of dystonia, and (3) had medical implants contraindicated for undergoing MRI. Finally, thirty-five patients diagnosed with ICD and twenty-eight age and sex matched healthy controls were recruited. All participants completed safety questionnaires for MRI and provided written informed consent in accordance with the Declaration of Helsinki. This study was approved by the Ethics Committee of First Teaching Hospital of Tianjin University of Traditional Chinese Medicine.

Motor symptom severity in ICD patients was assessed using the motor subscore of the Toronto western spasmodic torticollis rating scale (TWSTRS) (23). Disease-specific quality of life was evaluated using the Craniocervical dystonia questionnaire-24 (CDQ-24) (24), which is specifically designed for patients with craniocervical dystonia and is sensitive to both motor and psychosocial burden. To assess mood symptoms, we used the Hamilton depression rating scale (HAMD) (25), and the self-rating anxiety scale (SAS) (26). Importantly, the HAMD and SAS assessments were administered to both ICD patients and healthy controls, enabling group-level comparisons of depressive and anxiety symptoms, and allowing us to explore the nonmotor symptom burden characteristic of dystonia.

All participants underwent high-resolution T1-weighted structural MRI and resting-state fMRI scans on a 3 T GE Discovery MR750 scanner. All subjects underwent scans with foam padding and earplugs to minimize head movement and reduce scanner noise. A high-resolution T1-weighted brain volume (BRAVO) 3D MRI sequence with 128 contiguous sagittal slices was performed with the scan parameters of repetition time (TR) = 8.2 ms; echo time (TE) = 3.2 ms; inversion time = 450 ms; field of view (FOV) = 256 × 256 mm²; slice thickness = 1.0 mm; no gap; flip angle (FA) = 12°; matrix = 256 × 256; and voxel size = 1 mm³. The resting-state fMRI scans were performed by an echo planar imaging (EPI) sequence with the scan parameters of TR = 2000 ms, TE = 30 ms, FA = 90°, matrix = 64 × 64, FOV = 240 × 240 mm², slice thickness = 3 mm, gap = 1 mm. Each brain volume comprised 36 axial slices, and each functional run contained 188 volumes. During the fMRI scans, all of the subjects were instructed to keep their eyes closed, to relax and to move as little as possible.

VBM preprocessing was performed using the CAT12 toolbox with the accompanying methodology: bias correction, segmentation, the creation of population-specific tissue templates, spatial normalization using the DARTEL technique, and smoothing with an 8 mm × 8 mm × 8 mm full-width. After these preprocessing steps, we acquired the normalized, modulated, and smoothed GMV images, and each voxel represented volume information. Total intracranial volume (TIV) was computed for each participant.

Resting-state fMRI data were preprocessed automatically using the DPABI software (27). The first 10 volumes were discarded to allow for signal equilibrium. Preprocessing steps included slice-timing correction, realignment for head motion correction, and exclusion of participants with head movements exceeding 2 mm in translation or 2° in rotation. Images were then spatially normalized to the MNI EPI template via SPM12 and resampled to 3 × 3 × 3 mm³ voxels, followed by smoothing with an 8 mm FWHM Gaussian kernel. The resulting time series underwent linear detrending and band-pass filtering (0.01–0.08 Hz). Nuisance covariates, including Friston-24 head motion parameters, white matter signal, cerebrospinal fluid signal, and global signal, were regressed out. Framewise displacement (FD) value was calculated for each participant (28), and scrubbing was performed by removing time points with FD > 0.2 mm to minimize motion-related artifacts.

We conducted the voxel-based comparisons to distinguish the brain regions that showed group differences in GMV by utilizing the two-sample t-test, with age, gender, FD, TIV, age-onset and disease duration as covariates. The assessment of whether this group difference is significant was performed using a non-parametric permutation-based two-sample t-test (with 5,000 permutations) and further employing threshold-free cluster enhancement (TFCE) combined with false discovery rate (FDR) correction at p < 0.05 to control for false positives due to multiple comparisons at the voxel level.

Seed-based functional connectivity analysis was conducted using the regions showing significant GMV differences as region of interests (ROIs). Thus, ROI-based FC analysis was performed to calculate the FC maps of each ROI using the DPABI software. For each subject, the correlation coefficients between the mean time series of each ROI and that of each voxel of the whole brain were computed and converted to z-values, using Fisher’s r-to-z transformation to improve the normality. Similar to the VBM analysis, the two-sample t-test analysis was performed to test the group difference between the rsFCs of each of the ROIs, age, gender, FD, TIV, age-onset and disease duration were entered as covariates. The correction for multiple comparisons was performed using FDR correction at p < 0.05, cluster size >30.

In addition to neuroimaging analyses, we conducted supplementary statistical analyses in R v4.4.3¹ to investigate clinical variable relationships and group-level differences. Group comparisons of age, HAMD, and SAS between ICD patients and healthy controls were performed using two-sample t-tests, while group differences in sex were evaluated using chi-square (χ²) tests. These analyses were conducted using the stats package. Statistical significance was defined as two-tailed p < 0.05.

We computed Pearson correlation coefficients among clinical measures, including TWSTRS motor score, CDQ-24, HAMD, and SAS, to explore the associations between motor symptoms, disease-related quality of life, and mood disturbances. Correlation matrices were generated using the psych package, and significance was set at two-tailed p < 0.05.

To examine brain behavior associations, we extracted mean GMV from the significant GMV differences identified in the VBM analysis. These regional GMV values of ICD patients were then correlated with clinical scores (TWSTRS, CDQ-24, HAMD, SAS) using partial Pearson correlation (controlling for age, sex, FD, TIV, age-onset and disease duration) via the ppcor package in R.

To further explore the potential mediating role of significant GMV differences regions in the relationship between disease status and mood symptoms, we performed mediation analyses using structural equation modeling. Specifically, we tested whether GMV in regions showing significant group differences mediated the association between group status (ICD patients vs. healthy controls; coded as a binary variable) and mood symptom severity, as assessed by HAMD and SAS scores.

The mediation models were constructed with group status as the independent variable (X), regional GMV as the mediator (M), and mood symptoms (HAMD or SAS) as the dependent variable (Y). Age and gender were included as covariates in all models. Standardized estimates were calculated using the lavaan package in R v4.4.3. The significance of indirect effects was assessed via bootstrapping procedures (5,000 resamples), and p-values <0.05 were considered statistically significant.

A total of 63 participants were included in the analysis, comprising 35 patients with ICD and 28 healthy controls. The demographic and clinical characteristics of the participants are summarized in Table 1. There was no significant difference in age, gender and HAMD scores between the ICD group and the healthy controls group. However, the ICD group showed significantly higher anxiety levels compared to healthy controls, as measured by the SAS (p < 0.001). As TWSTRS and CDQ-24 are disease-specific scales, they were only administered in the ICD group, yielding mean scores of 30.90 ± 13.69 and 53.66 ± 18.57, respectively.

Compared to healthy controls, ICD patients exhibited significantly decreased GMV in the left paracentral lobule (PCL, peak MNI coordinates: x = −9.3, y = −22.1, z = 54; cluster size = 196 voxels; t = −3.33; Cohen’s d = −0.84, Figure 1A), and right middle temporal gyrus (MTG, peak MNI coordinates: x = 49.1, y = −14.2, z = −18; cluster size = 227 voxels; t = −3.29; Cohen’s d = −0.83, Figure 1B). These differences remained significant after TFCE and FDR correction (p < 0.05). In addition, there were no regions demonstrating a significant increased GMV in ICD patients.

Seed-based rsFC analysis was conducted using the left PCL and the right MTG as ROIs with significantly reduced GMV in ICD patients. Group-level comparisons revealed both decreased and increased functional connectivity patterns in ICD patients compared to healthy controls. All results were significant after FDR correction at p < 0.05, cluster size > 30.

Compared to the healthy controls, decreased rsFC was observed between left PCL and the left rectus (REC, peak MNI coordinates: x = −7, y = 49, z = −20; cluster size = 38 voxels; t = −3.46; Cohen’s d = −0.88, Figure 2A), right temporal pole (TP, peak MNI coordinates: x = 31, y = 5, z = −23; cluster size = 67 voxels; t = −3.31; Cohen’s d = −0.84, Figure 2B), increased rsFC was identified between left PCL and the bilateral thalamus (THA, peak MNI coordinates: x = −4, y = −21, z = 3; cluster size = 36 voxels; t = 3.00; Cohen’s d = 0.76, Figure 2C) in the ICD patients.

Similarly, we found decreased rsFC between right MTG and the left TP (peak MNI coordinates: x = −30, y = 4, z = −24; cluster size = 47 voxels; t = −2.73; Cohen’s d = −0.69, Figure 3A), increased rsFC between left MTG and the right middle frontal gyrus (MFG, peak MNI coordinates: x = 35, y = 12, z = 45; cluster size = 145 voxels; t = 2.52; Cohen’s d = 0.64, Figure 3B) in ICD patients. These findings suggest disrupted connectivity between temporal/frontal regions and key affective and executive hubs, with both hypo- and hyper-connectivity potentially reflecting maladaptive reorganization in response to underlying structural deficits.

We conducted correlation analyses to investigate the relationships among clinical variables and their associations with regional GMV alterations in patients with ICD. As shown in Figure 4A, the TWSTRS motor score was significantly positively correlated with CDQ-24 (r = 0.60, p < 0.001), HAMD (r = 0.51, p < 0.01), and SAS (r = 0.35, p < 0.05), indicating that greater motor severity was associated with poorer quality of life and higher levels of depressive and anxiety symptoms. The CDQ-24 score showed robust positive correlations with both HAMD (r = 0.77, p < 0.001) and SAS (r = 0.57, p < 0.001), highlighting the strong contribution of mood disturbances to the perceived disease burden in ICD. Furthermore, HAMD and SAS scores were also correlated with each other (r = 0.55, p < 0.001). In terms of brain-behavior relationships, partial correlation analyses controlling for age and sex revealed that GMV in the left PCL was negatively correlated with TWSTRS (r = −0.35, p = 0.039, Figure 4B) and SAS (r = −0.42, p = 0.013) (Figure 4B). These findings underscore the close interplay between structural brain changes and clinical manifestations in ICD, and reinforce the role of affective symptoms in shaping disease burden.

We examined whether GMV in the left PCL or right MTG mediated the association between group status and mood symptoms (HAMD and SAS). Among the four tested models, only the pathway involving GMV in the left PCL and anxiety (SAS) demonstrated a significant indirect effect (Figure 4C). Specifically, patients with ICD showed significantly reduced GMV in the left PCL compared to controls (path a: standardized β = −0.570, p < 0.001, 95% CI [−0.709, −0.432]). In turn, lower PCL GMV was associated with higher anxiety scores (path b: standardized β = −0.464, p < 0.001, 95% CI [−0.705, −0.223]). The indirect effect (ab) was significant (standardized β = 0.265, p < 0.001, 95% CI [0.125, 0.405]), yielding a proportion mediated of approximately 53.8%. These findings indicate that GMV alterations in the left PCL statistically account for a substantial portion of the group difference in anxiety symptoms. No significant statistical mediation effects were observed in the models involving HAMD scores or MTG/PCL and SAS scores or MTG (all p > 0.05).