Microstructural deficits of the thalamus in major depressive disorder

Abstract Macroscopic structural abnormalities in the thalamus and thalamic circuits have been implicated in the neuropathology of major depressive disorder. However, cytoarchitectonic properties underlying these macroscopic abnormalities remain unknown. Here, we examined systematic deficits of brain architecture in depression, from structural brain network organization to microstructural properties. A multi-modal neuroimaging approach including diffusion, anatomical and quantitative MRI was used to examine structural-related alternations in 56 patients with depression compared with 35 age- and sex-matched controls. The seed-based probabilistic tractography showed multiple alterations of structural connectivity within a set of subcortical areas and their connections to cortical regions in patients with depression. These subcortical regions included the putamen, thalamus and caudate, which are predominantly involved in the limbic-cortical-striatal-pallidal-thalamic network. Structural connectivity was disrupted within and between large-scale networks, including the subcortical network, default-mode network and salience network. Consistently, morphometric measurements, including cortical thickness and voxel-based morphometry, showed widespread volume reductions of these key regions in patients with depression. A conjunction analysis identified common structural alternations of the left orbitofrontal cortex, left putamen, bilateral thalamus and right amygdala across macro-modalities. Importantly, the microstructural properties, longitudinal relaxation time of the left thalamus was increased and inversely correlated with its grey matter volume in patients with depression. Together, this work to date provides the first macro–micro neuroimaging evidence for the structural abnormalities of the thalamus in patients with depression, shedding light on the neuropathological disruptions of the limbic-cortical-striatal-pallidal-thalamic circuit in major depressive disorder. These findings have implications in understanding the abnormal changes of brain structures across the development of depression.

criteria, and met the following criteria: 1) age range of 18 to 65 years; 2) for male/female who is capable of fertility, s/he agreed to use effective contraceptive methods during the study period and within one month after the end of the study to ensure the effective contraception for himself/herself or sexual partner. The exclusion criteria were as follows: 1) participants met the DSM-V diagnostic criteria for another psychiatric illness, including schizophrenia, bipolar disorder, anxiety disorders, obsessive-compulsive disorder, physical symptoms, organic mental disorders or depression caused by hypothyroidism, substance abuse or dependence; 2) participants who have received electroconvulsive therapy (ECT) within three months before screening; 3) participants with any other unsuited conditions considered by investigators that could not participate the study. The severity of depressive symptoms was assessed using the Hamilton Depression Rating Scale (HAMD). Because of the high comorbidity of depression and anxiety 1 , we assessed the severity of anxiety for MDD patients using the Hamilton Anxiety Rating Scale (HAMA) and examined the influence of anxiety in subsequent analyses. All HCs were cognitively normal, free of neurological disease, and had no cognitive complaints. For MRI data quality checks, data of one HC were excluded from diffusion MRI (dMRI) data analysis due to the pool image quality, resulting in 56 MDD patients and 34 HCs included in the dMRI data analysis.

MRI data acquisition
High-resolution structural images were acquired with a three-dimensional (3D) T1-weighted (T1w) fast spoiled gradient-echo (FSPGR) sequence with TE = 2.9 ms, TR = 6.7 ms, FA = 12°, field of view (FoV) = 256 × 256 mm, 192 slices, and voxel size = 1 mm 3 . The dMRI data was collected using a single-shot spin-echo echo-planar imaging (EPI) sequence (TR = 8724 ms, TE = 81.4 ms, FA = 90°, voxel size = 2 × 2 × 2 mm 3 ). The diffusion weighting was isotropically distributed along 64 directions (b = 1000 s/mm 2 ). Ten non-diffusion weighted (b = 0) images were acquired at the beginning of the dMRI session. The quantitative MRI (qMRI) data was collected using the protocols described in a previous study 2 . In brief, the quantitative macromolecular tissue volume (MTV) and longitudinal relaxation time (T1) were measured from four spoiled gradient echo (SPGE) images with flip angles of 4°, 10°, 20°, 30° (TR = 14 ms, TE = 2 ms) at a 1 × 1 mm 2 in-plane resolution with a slice thickness of 2 mm. We collected four additional spin-echo inversion recovery (SEIR) scans with an EPI read-out, a slab inversion pulse, and spectral spatial fat suppression to remove field inhomogeneities. The SEIRs were acquired with a TR of 3 s, echo time set to minimum full, and 2x acceleration. The inversion times were 50, 400, 1200, and 2400 ms with a 2 × 2 mm 2 in-plane resolution and a slice thickness of 4 mm.

Data preprocessing
The dMRI data were preprocessed using the FMRIB's Diffusion Toolbox (FDT) from FSL package 3,4 (version 5.0.9; https://fsl.fmrib.ox.ac.uk/fsl). Briefly, diffusion-weighted images were first corrected for eddy-current-induced distortions and head movements by using affine transformation of each diffusionweighted image to non-diffusion-weighted (b = 0) image. The diffusion tensor models were then fitted at each voxel of the brain on distortion-corrected images to yield three eigenvectors and three eigenvalues.

Network node definition
To determine the nodes of structural network, an automated anatomical labeling (AAL) atlas was used to parcellate the whole cerebral cortex into 90 regions/nodes (Supplementary Table 1). The mask for each brain region was extracted and transformed into each participant's diffusion space. In brief, ac-pc aligned T1w image in individual participant' s space was first normalized to the ICBM 152 template in Montreal Neurological Institute (MNI) space by applying nonlinear transformation between T1w structural space and standard space. Second, the T1w image was coregistered to participant' s nondiffusion-weighted (b = 0) image by using linear transformation between diffusion space and the T1w structural space. Here, relevant transformation matrices and their inverses were derived and concatenated to produce transformation matrices between diffusion and standard space. The concatenated transformation matrices were then used to warp AAL labels from MNI space to native diffusion space. This allowed us to perform inter-subject comparison in diffusion space and produce an individual-specific connectivity matrix. To this end, 90 parcellated regions for each participant in native diffusion space were produced, representing 90 nodes in the structural network, and used for structural network construction.

Network construction
To reconstruct the structural network of white matter tracts between 90 AAL regions, probabilistic tractography was performed on the preprocessed dMRI dataset using PANDA pipeline 5 (http://www.nitrc.org/projects/panda) in MATLAB2018b. We used a within-voxel probabilistic diffusion model to build up distributions on diffusion parameters at each voxel of the dataset 6 . The probability distribution estimates the multi-fiber tract orientation as well as its uncertainty in each voxel and guides multiple fiber samples starting from a seed voxel to a specified target region.
Probabilistic tractography was then applied by drawing 5000 individual streamlines (maximum step: 2000, step length: 0.5mm, curvature threshold: 0.2) in each voxel in native diffusion space. Structural connectivity between two regions was measured by calculating connectivity probability from each seed region to each of the other target regions. In the present study, each AAL region was selected as the seed and the remaining 89 regions as targets. The tracts were terminated as soon as leaving the brain or reaching a particular target region. This resulted in a 90 × 90 connectivity matrix with inter-regional connection probability, representing a structural network for each participant. For this matrix, the probability from the i-th region to the j-th region is not necessarily equivalent to the one from j to i. However, these two probabilities are highly correlated across the cerebral cortex. Thus, a unidirectional 3 symmetric matrix for each participant was acquired by averaging these two probabilities. In native diffusion space, the resulted structural network is composed of nodes, representing a pair of regions, and edges, representing structural connectivity among the nodes, i.e., the region (i) and another region (j) were connected through an edge (eij = [i, j]), in case of at least one fiber existed between them 7,8 .
Consequently, given the definition of nodes and edges for the structural connectivity network (matrix), the connectivity probability was calculated as an edge between pair of regions 6,[9][10][11] .

NBS analysis
The NBS approach was applied for group comparisons by using NBS Connectome (http://www.nitrc.org/projects/nbs/). The NBS is a non-parametric statistical method to deal with the multiple comparisons problems by controlling for the FWE rate in a weak sense while performing mass univariate hypothesis testing on all connections. With the NBS, the null hypothesis was evaluated at a level of interconnected subnetworks rather than individual connections. Briefly, a two-sample t-test at each edge was conducted independently to test for significant differences in the value of connectivity between the two groups. A primary threshold (t ≥ 2.5, which corresponds to uncorrected p ≤ 0.01) was then applied to define a set of supra-threshold connections. This step identified all the possible connected components in the matrix at the uncorrected level, followed by computing the size of the connected components which link supra-threshold connections with a set of nodes in the network. To estimate the statistical significance of size for each component, the null distribution of maximal connected component size was empirically derived using non-parametric permutation testing with 5000 times. Supra-threshold connections that were significant at the level of p < 0.01 were reported. The BrainNet viewer (http://www.nitrc.org/projects/bnv/) was used to visualize and display the significant connections 12 .

Structural connectivity analysis of the thalamus with the whole brain.
A seed-based probabilistic tractography was conducted to evaluate white matter structural connectivity of the thalamus with the whole brain in patients with MDD. The steps followed similar procedures with the seed-target probabilistic tractography as described above (see Network construction) but only each unilateral thalamus was defined as a seed and voxels of the rest of the whole brain were set as targets. The unilateral thalamus seed mask adapted for the analysis was derived from the individual T1 parcellation based on AAL atlas. A two-sample t-test was conducted to compare the thalamus-centered structural connectivity between MDD patients and HCs. The multiple comparison correction was performed using threshold-free cluster enhancement (TFCE), a non-parametric method to enhance detectability of neuroimaging signal on cluster-based thresholding by integrating information of cluster extend and height into voxel-wise statistical inference 13 . The TFCE-based analysis was performed with 5000 permutation tests and the resulting TFCE-statistic map was considered significant at cluster-wise p < 0.05, FWE-corrected.

Voxel-based morphometry analysis
To detect macroscopic abnormalities of brain volumes in MDD, voxel-based morphometry (VBM) analysis was performed using VBM8 toolbox (http://dbm.neuro.uni-jena.de/vbm8/) implemented in SPM 12 (Statistical Parametric Mapping, http://www.fil.ion.ucl.ac.uk/spm/software/spm12) software running on MATLAB R2018 for Linux. All ac-pc aligned T1w images were segmented into the gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) in native space. The segmented images were then normalized to MNI space with a dimension of 121 × 145 × 121 and a spatial resolution of 1.5 × 1.5 × 1.5 mm using high-dimensional DARTEL normalization. The normalized segmented images for GM and WM were subsequently modulated by Jacobian determinants derived from the spatial normalization, which allows to estimate the GM volume (GMV) and WM volume (WMV) while controlling for overall brain size. After preprocessing, sample homogeneity was checked to identify outliers which were two standard deviations outside of the mean. Finally, these normalized, segmented, modulated volumes were spatially smoothed with a 6-mm full-width half maximum (FWHM) Gaussian kernel. To avoid possible edge effects around the margin between different tissue types, all voxels with a probability value < 0.2 (absolute threshold; range 0-1) for both GM and WM were excluded. For statistical analysis, a two-sample t-test was conducted to identify significant differences in GMV and WMV between MDD patients and HCs, respectively, and the TFCE method with 5000 permutation tests was used for correcting the multiple comparisons. The resulted TFCEstatistic map was thresholded at p < 0.05, FWE-corrected.

Surface-based morphometry analysis.
The synthetic T1w images generated from mrQ analysis pipeline were then processed for surface-based morphometry (SBM) analysis by using Freesurfer v5.3.0 (http://surfer.nmr.mgh.harvard.edu) with standard recon-all procedures. The Freesurfer pipeline performs the volumetric GM and WM segmentation, providing several automated cortical parcellations and assigning a neuroanatomical label to each location based on prior atlas. The processing steps included skull stripping, atlas registration, spherical surface registration and parcellation. In this study, cortical parcellation and subcortical segmentation were based on the Freesurfer built-in atlas 14,15 , resulting in a total of 82 distinct brain regions, including 74 cortical, 7 subcortical regions (the thalamus, caudate, putamen, pallidum, hippocampus, amygdala and nucleus accumbens) per hemisphere and the brainstem. The resulted segmented brain tissues were used for subsequent region of interest (ROI) analysis of qMRI modality.
After the recon-all procedures of all participants, cortical thickness maps were smoothed using Gaussian Kernel with FWHM of 10 mm. The vertex-wise group comparison for cortical thickness was performed in Freesurfer using general linear model (GLM) with age, gender and years of education as covariates. Moreover, the vertex-wise results were Bonferroni corrected for multiple comparisons with a threshold of 2, corresponding to p < 0.01. Cluster-wise statistical threshold of p < 0.05 was used after Bonferroni correction. A statistical thickness difference map was constructed using -log10 (p value).
Cluster-wise p values (CWP) were reported. The volumes of subcortical regions were automatically derived from Freesurfer's output report for group comparisons.

Data preprocessing
Both the SPGE and SEIR scans were processed by using the mrQ software package (https://github.com/mezera/mrQ) in MATLAB to produce the evaluation of MTV and quantitative T1 maps for each participant. The RF coil bias was corrected by combining SPGE scans with a set of lowresolution unbiased SEIR-EPI scans, producing accurate proton density (PD) and unbiased T1 fits across the brain. The MTV maps were estimated by calculating the fraction of non-water volume in each voxel while CSF voxels were approximated as entirely filled with water 1 . The synthetic T1w images which were spatially matched with MTV maps and optimized for segmentation of brain tissues, were also computed by using mrQ package and then processed using Freesurfer for ROI analysis of qMRI modality (see Surface-based morphometry analysis).

ROI analysis for qMRI data.
ROIs were segmented and extracted from Freesurfer's standard pipeline for each participant and saved as volumetric binary masks. All of these masks were generated in native MTV space and matched well with participants' synthetic T1w images through manual inspection. The individual ROI masks were then applied on corresponding MTV and T1 maps for each participant, and averaged MTV and T1 values across voxels within each ROI were computed.

Correlation analysis of qMRI measurements.
Regions with group differences in MTV/T1 values were included for this analysis. The MTV/T1 values within each region were respectively correlated with its GMV and FA values across the two groups.
Briefly, each region was defined as a ROI using a volumetric 1-mm-diameter sphere centered on its peak MNI coordinates with the strongest statistical significance based on the results of VBM analysis.
To extract FA values within each ROI, FA maps of all participants were aligned onto the FMRIB58_FA_1mm template in standard space by using nonlinear registration 16 . The sphere ROI masks were then applied to the normalized FA maps for all participants. The FA values across voxels within each sphere ROI were averaged and extracted. The GMV values within each ROI were averaged and extracted from the preprocessed GM volumes generated from the VBM analysis for all participants. The Pearson's correlation coefficient (r) and p value between the MTV/T1 and GMV and FA values for each ROI were calculated. The statistical significance was set at p < 0.05, FWE corrected.
We also correlated qMRI measurements with HAMD and HAMA scores to test the influence of anxiety on microstructural alternations in MDD. Because that data of some MDD patients were missing, the correlation analysis was only conducted based on the existing scale data (n = 36). HAMA and HAMD scores of each MDD patient nearest to the time of MRI scanning were selected. The demographic data (age, gender and education years) of all MDD patients were included as covariates. 6

Supplementary structural network analysis
A supplementary structural network analysis was performed to evaluate the effect of participants' demographic data on our dMRI results. The procedures performed in the supplementary analysis were similar to the structural network analysis mentioned above (see Network construction) but the effect of participants' age, gender and education years were removed before two-sample t-test.

Supplementary VBM analysis
Since there were significant differences in education years between the two groups, we additionally conducted a supplemental VBM analysis to remove the effect of participants' demographic data on the GMV/WMV measures such that further validate the outcome of our study. The procedures of the supplemental VBM analysis were similar to those mentioned above (see Voxel-based morphometry analysis) but participants' demographic data (age, gender and education years) were classified as covariates and subsequently regressed out.

Alternations of the thalamus-centered structural connectivity in MDD patients
To identify specific pattern of thalamic structural connectivity in MDD patients, we extracted the bilateral thalamus-to-whole-brain connections identified by the NBS analysis (Supplementary Figure   1A). These connections formed a thalamus-centered network, primarily involving the bilateral  Figure 1B). We found that the target regions structurally connecting to the thalamus were located in the cortical regions e.g., the bilateral middle frontal gyrus, middle temporal gyrus, dorsolateral superior frontal gyrus, precentral and postcentral gyrus, precuneus, inferior temporal gyrus, superior temporal gyrus, superior medial frontal gyrus, insula, inferior orbitofrontal gyrus and anterior cingulate and paracingulate gyri as well as subcortical regions including the bilateral putamen, hippocampus, parahippocampal gyrus, caudate and pallidum. No significant decreased structural connectivity was survived after multiple comparison correction in MDD patients relative to HCs.

The correlation of qMRI measurements across 2 groups
The correlation analysis showed that T1s were negatively correlated with GMV values in the THA.L (r = -0.349, p = 0.002; Supplementary Figure 2) across two groups. There was no significant correlation observed between T1s of the THA.L with HAMA (r = -0.259, p = 0.127) or HAMD scores (r = -0.069, p = 0.689) in MDD patients.

Results of supplementary structural network analysis
By using the NBS method, we only identified a single network after correcting for all the covariates suggest altered structural connectivity in MDD patients is unlikely to be effect by participants' age, gender and education years.

Results of supplementary VBM analysis
The results of the supplementary VBM analysis were comparable to those in the main text (Supplementary Figure 4 and Table 3). A total of 23 clusters showing significant reduced GMV were 8 identified in MDD patients relative to HCs, while WMV reductions were found in eight clusters (Supplementary Figure 4A). Additionally, the reduced GMV values of the key regions e.g., the bilateral thalamus, bilateral insula, left parahippocampal gyrus (PHG.L), right amygdala (AMYG.R), left inferior orbitofrontal gyrus (ORBinf.L) and ORBsup.L were still observed even participants' demographic data were regressed out. The GMV values of these key regions were also extracted from the identified clusters for visualization (Supplementary Figure 4B)