Introduction
Major depressive disorder (MDD) is a highly heterogeneous syndrome characterized by significant intersubject variability in symptom manifestations, clinical trajectories, etiologies, and treatment responses (Bondar, Caye, Chekroud, & Kieling, Reference Bondar, Caye, Chekroud and Kieling2020; Drysdale, Grosenick, & Downar, Reference Drysdale, Grosenick and Downar2017; Krishnan & Nestler, Reference Krishnan and Nestler2008; Nguyen, Harder, Xiong, Kowalec, & Hägg, Reference Nguyen, Harder, Xiong, Kowalec and Hägg2022). This clinical heterogeneity complicates individualized diagnosis and treatment evaluations, and mechanism research in MDD. Neuroimaging studies have shown that MDD exhibits substantial intersubject variations in neuroimaging metrics, such as gray matter morphology (Chen et al., Reference Chen, Peng, Sun, Kuang, Li, Jia and Gong2016). This high variability hampers the ability to reach consistent findings in MDD research. A critical reason for this inconsistency is the reliance on case–control designs that focus on group-level effects and overlook intersubject heterogeneity (Wolfers et al., Reference Wolfers, Doan, Kaufmann, Alnæs, Moberget, Agartz, Buitelaar, Ueland, Melle, Franke, Andreassen, Beckmann, Westlye and Marquand2018). To address this heterogeneity, many attempts have been made, such as subtype analysis based on symptomatology or neuroimaging phenotype and even inferring subject-level neuroimaging differences (Han et al., Reference Han, Cui, Zheng, Li, Zhou, Fang, Sheng, Wen, Liu, Wei, Chen, Chen, Cheng and Zhang2023a; Wen et al., Reference Wen, Fu, Tosun, Veturi, Yang, Abdulkadir, Mamourian, Srinivasan, Skampardoni, Singh, Nawani, Bao, Erus, Shou, Habes, Doshi, Varol, Mackin and Sotiras2022). However, as a brain connectome disorder, MDD is characterized by abnormal structural and functional connections among large-scale brain systems (Fornito, Zalesky, & Breakspear, Reference Fornito, Zalesky and Breakspear2015; Kaiser, Andrews-Hanna, Wager, & Pizzagalli, Reference Kaiser, Andrews-Hanna, Wager and Pizzagalli2015). The distribution of intersubject variability in brain connectome remains largely understudied.
Brain function depends on effective structural and functional connectomes between distinct brain systems. Neuroimaging studies typically assess the brain connectome using techniques such as intrinsic functional connectivity and structural covariance (Alexander-Bloch, Giedd, & Bullmore, Reference Alexander-Bloch, Giedd and Bullmore2013; Fox & Raichle, Reference Fox and Raichle2007). In healthy individuals, brain connectome architecture exhibits considerable intersubject variability, distinguishing individuals from others and underlying differences in cognition and behavior (Finn et al., Reference Finn, Shen, Scheinost, Rosenberg, Huang, Chun, Papademetris and Constable2015; Smith et al., Reference Smith, Nichols, Vidaurre, Winkler, Behrens, Glasser, Ugurbil, Barch, Van Essen and Miller2015). This variability is modulated by genetic factors (Gao et al., Reference Gao, Elton, Zhu, Alcauter, Smith, Gilmore and Lin2014). In a healthy brain, connectome variability shows a specific distribution pattern, with higher variability in the heteromodal association cortex and lower variability in unimodal cortices (Finn et al., Reference Finn, Shen, Scheinost, Rosenberg, Huang, Chun, Papademetris and Constable2015, Smith et al., Reference Smith, Nichols, Vidaurre, Winkler, Behrens, Glasser, Ugurbil, Barch, Van Essen and Miller2015). This pattern is consistent across adults (Mueller et al., Reference Mueller, Wang, Fox, Yeo, Sepulcre, Sabuncu, Shafee, Lu and Liu2013) and neonates (Gao et al., Reference Gao, Elton, Zhu, Alcauter, Smith, Gilmore and Lin2014; Stoecklein et al., Reference Stoecklein, Hilgendorff, Li, Forster, Flemmer, Galie, Wunderlich, Wang, Stein, Ehrhardt, Dietrich, Zou, Zhou, Ertl-Wagner and Liu2020) and is influenced by neuropsychiatric disorders (Dumlu, Ademoglu, & Sun, Reference Dumlu, Ademoglu and Sun2020; Sun et al., Reference Sun, Liu, Ma, Duan, Wang, Xu, Xu, Xu, Wang, Tang, He and Xia2021). For example, patients with schizophrenia demonstrate higher intersubject variability of functional connectome (IVFC) than healthy controls in bilateral sensorimotor, visual, auditory, and subcortical regions. Additionally, altered IVFC is associated with clinical heterogeneity (Sun et al., Reference Sun, Liu, Ma, Duan, Wang, Xu, Xu, Xu, Wang, Tang, He and Xia2021; Zhang, Xu, Ma, Qian, & Zhu, Reference Zhang, Xu, Ma, Qian and Zhu2024). In temporal lobe epilepsy, altered IVFC is linked to impairment in attention, memory, and executive functioning (Dumlu et al., Reference Dumlu, Ademoglu and Sun2020). These findings have potential implications for the understanding of clinical heterogeneity and individualized clinical decision-making in these disorders. To our knowledge, only three studies have investigated altered intersubject variability of brain connectome in MDD (Gai et al., Reference Gai, Chu, Li, Guo, Ma, Shi, Che, Zhao, Dong, Li, Xie and Mao2024a; Hou et al., Reference Hou, Jiang, Li, Liu, Hou, Yin, Zhang, Zhang, Xie, Zhang, Kong and Yuan2023; Sun et al., Reference Sun, Liu, Ma, Duan, Wang, Xu, Xu, Xu, Wang, Tang, He and Xia2021). These studies have two main limitations: they use a single dataset with a small sample size, limiting statistical power and reproducibility, and they focus solely on functional connectome, neglecting the structural connectome.
In this study, we aimed to investigate intersubject variability of functional (IVFC) and structural covariance connectome (IVSC) using a large, multi-site dataset comprising 1,276 patients with MDD and 1,104 matched healthy controls (HCs). First, we identified abnormal patterns of IVSC and IVFC in patients with MDD. Next, we examined associations between these abnormal patterns and clinical features, such as age of onset, medicine status, and symptom severity. Additionally, we explored the potential molecular underpinnings of abnormal patterns of IVSC and IVFC. Finally, we conducted reproducibility analyses to evaluate the robustness and reproducibility of these abnormal patterns and performed disease specificity analysis to determine whether these patterns were specific to MDD compared with other psychiatric disorders.
Materials and methods
Imaging dataset and preprocessing
The MDD imaging datasets used in this study were sourced from the REST-meta-MDD consortium (http://rfmri.org/REST-meta-MDD) (Chen et al., Reference Chen, Lu, Li, Li, Wang, Castellanos, Cao, Chen, Chen, Cheng, Cui, Deng, Fang, Gong, Guo, Hu, Kuang, Li, Li and Yan2022; Yan et al., Reference Yan, Chen, Li, Castellanos, Bai, Bo, Cao, Chen, Chen, Chen, Cheng, Cheng, Cui, Duan, Fang, Gong, Guo, Hou, Hu and Wang2019), which includes data from 25 research sites and comprises 1,276 patients diagnosed with MDD (age: 36.23 ± 21.38 years, females: 63.71%) and 1,104 HCs (age: 36.15 ± 24.55 years, females: 58.06%). Data collection for all research sites was approved by local research ethics committees. Patients were diagnosed by experienced psychiatrists based on DSM-IV criteria for MDD. The severity of symptoms in patients was assessed using the Hamilton Depression Rating Scale (HAMD) (Hamilton, Reference Hamilton1960). Among the patients with MDD, 538 were first-episode individuals, 282 were recurrent cases, 447 people were medication-naive and 408 were treated. Further details about demographic and clinical characteristics, as well as image acquisition protocols are thoroughly described in previous studies (Chen et al., Reference Chen, Lu, Li, Li, Wang, Castellanos, Cao, Chen, Chen, Cheng, Cui, Deng, Fang, Gong, Guo, Hu, Kuang, Li, Li and Yan2022; Yan et al., Reference Yan, Chen, Li, Castellanos, Bai, Bo, Cao, Chen, Chen, Chen, Cheng, Cheng, Cui, Duan, Fang, Gong, Guo, Hou, Hu and Wang2019).
Structural and functional MRI images were preprocessed using standard analysis protocols across sites to minimize analytic variations. The preprocessing steps of structural and functional MRI images are detailed elsewhere (Chen et al., Reference Chen, Lu, Li, Li, Wang, Castellanos, Cao, Chen, Chen, Cheng, Cui, Deng, Fang, Gong, Guo, Hu, Kuang, Li, Li and Yan2022; Yan et al., Reference Yan, Chen, Li, Castellanos, Bai, Bo, Cao, Chen, Chen, Chen, Cheng, Cheng, Cui, Duan, Fang, Gong, Guo, Hou, Hu and Wang2019). In this study, we directly utilized the released gray matter volume (GMV) maps obtained via voxel-based morphometry analysis (Ashburner, Reference Ashburner2009) and extracted averaged time courses of each brain region defined in the Craddock 200 brain atlas from the REST-meta-MDD consortium (http://rfmri.org/REST-meta-MDD) (Craddock, James, Holtzheimer, Hu, & Mayberg, Reference Craddock, James, Holtzheimer, Hu and Mayberg2012). We used the time courses without global signal regression for the reason being its implication in the pathogenesis of mental disorders (Dong et al., Reference Dong, Wang, Zhou, Chang, Qiu, Feng, He, Lei and Chen2023; Han et al., Reference Han, Wang, He, Sheng, Zou, Li, Yang, Guo, Fan, Guo, Lu, Cui and Chen2019; Yang et al., Reference Yang, Murray, Repovs, Cole, Savic, Glasser, Pittenger, Krystal, Wang, Pearlson, Glahn and Anticevic2014).
Construction of structural covariance and functional network and harmonization
The construction of intra-individual structural covariance networks followed the methodology outlined in a previous study (Yun et al., Reference Yun, Boedhoe, Vriend, Jahanshad, Abe, Ameis, Anticevic, Arnold, Batistuzzo, Benedetti, Beucke, Bollettini, Bose, Brem, Calvo, Cheng, Cho, Ciullo, Dallaspezia and Kwon2020). The mean gray matter volumes (GMVs) of 200 brain regions, defined in the Craddock 200 brain atlas, were corrected for age, sex, and site effects (Craddock et al., Reference Craddock, James, Holtzheimer, Hu and Mayberg2012). The resulting residuals were transformed to Z-scores using mean and standard deviation (SD) values of each brain region calculated from healthy controls (Vuoksimaa et al., Reference Vuoksimaa, Panizzon, Chen, Fiecas, Eyler, Fennema-Notestine, Hagler, Franz, Jak, Lyons, Neale, Rinker, Thompson, Tsuang, Dale and Kremen2016). A measure of joint variation between 200 morphometric features, represented as edge weights (distributed between 0 and 1), was calculated to form the network (Yun, Jang, Kim, Jung, & Kwon, Reference Yun, Jang, Kim, Jung and Kwon2015). This procedure yielded a 200 × 200 structural covariance connectome matrix for each subject.
Functional networks were constructed by calculating pairwise Pearson’s correlation coefficients between brain regions, with the resulting values transformed into Fisher’Z scores. This process produced a 200 × 200 functional connectome matrix for each subject.
To minimize the site effect introduced by differences in scanners and acquisition protocols, the ComBat method was independently performed in structural and functional connections (Johnson, Li, & Rabinovic, Reference Johnson, Li and Rabinovic2007). ComBat is a technique widely used in neuroimaging studies to remove unwanted site variation while preserving biological variability. This approach helps mitigate the impact of site differences on imaging metrics (Fortin et al., Reference Fortin, Parker, Tunç, Watanabe, Elliott, Ruparel, Roalf, Satterthwaite, Gur, Gur, Schultz, Verma and Shinohara2017).
Group differences in intersubject variability of the structural and functional connectome
The intersubject variability of the structural and functional connectome (IVSC/IVFC) was calculated to measure the inter-individual deviation of regional connection profiles for each group separately (Mueller, Wang, Fox, Yeo, & Liu, Reference Mueller, Wang, Fox, Yeo and Liu2013). The IVSC/IVFC for region k in subject i was dented as
$ {V}_{i,k} $
and was computed are follows:
$$ {V}_{i,k}=1-\mathrm{mean}\left(\mathrm{corr}\left({C}_{i,k},{C}_{j,k}\right)\right)\; for\;j=1,2,\dots n $$
For each subject i, we first calculated the Pearson’s correlation coefficients between the connection profiles of region k (
$ {C}_{i,k} $
, the k-th column of the connection matrix of subject i) and those of every other subject j (
$ {C}_{j,k} $
). Then,
$ {V}_{i,k} $
was estimated by subtracting the averaged correlation coefficient from 1. A larger
$ {V}_{i,k} $
indicates a greater deviation of regional connection profiles from its own group.
A two-tailed two-sample t-test was performed to investigate the group differences in IVFC and IVSC, correcting for age, sex, and site. Statistical significance was set at p < 0.05 with Bonferroni correction for multiple comparisons. For the comparison of IVFC, mean framewise displacement (FD) was additionally corrected due to its significant effect on measures of functional connectivity (Ciric et al., Reference Ciric, Rosen, Erus, Cieslak, Adebimpe, Cook, Bassett, Davatzikos, Wolf and Satterthwaite2018).
Associations among distributions of IVSC, gray matter morphological atrophy, and abnormal patterns of regional structural covariance connectome
We investigated the relationships among distributions of IVSC, gray matter morphological atrophy, and abnormal patterns of regional structural covariance connectome. Specifically, we examined whether brain regions with more severe gray matter morphological atrophy exhibited higher intersubject variability or vice versa. Regional gray matter morphological atrophy (unthresholded t-static vector) was obtained by comparing harmonized regional GMVs between HCs and MDD with s two-sample t-test. Abnormal patterns of regional structural covariance connectome (unthresholded t-static vector) were obtained by comparing averaged intra-individual SCs between MDD and HCs using a two-sample t-test. We then calculated paired Spearman’s correlation coefficients among these patterns.
Given the significant correlations between abnormal IVSC and gray matter atrophy, and between abnormal IVSC and IVFC (see results), we further hypothesized that IVSC mediates the relationship between gray matter atrophy and abnormal IVFC. To test this, we employed a standard three-variable mediation analysis using the Mediation Toolbox (https://github.com/canlab/MediationToolbox) (Wager, Davidson, Hughes, Lindquist, & Ochsner, Reference Wager, Davidson, Hughes, Lindquist and Ochsner2008). More details about this model were described elsewhere (Wager et al., Reference Wager, Davidson, Hughes, Lindquist and Ochsner2008). Mediation analysis assesses whether the relationship between two variables is explained by a third variable (the mediator). In this model, gray matter atrophy was the independent variable, IVFC was the dependent variable, and IVSC was the mediator. The significance of the mediation model was estimated using the bias-corrected bootstrap approach with 10,000 random samplings.
Associations between IVSC and IVFC and clinical features
To examine the impact of episodicity (first episode versus recurrent) and medication status (treated versus untreated) on the distributions of IVSC and IVFC, we conducted a two-sample t-test to compare subgroup differences. We also analyzed IVSC and IVFC differences between male and female patients, using a two-tailed two-sample t-test. Statistical significance was set at p < 0.05 with Bonferroni correction. Additionally, we computed Spearman’s correlation between previously identified abnormal IVFC/IVFC and the total score of HAMD.
Association between distributions of abnormal IVSC and IVFC and neurotransmitter receptors/transporters
Then, we examined the relationship between distributions of abnormal IVSC and IVFC and neurotransmitter receptors/transporters profiles to explore potential molecular underpinnings. Neurotransmitter receptors/transporters profiles were obtained from a PET-derived neurotransmitter receptors/transporters atlas (Hansen & Shafiei, Reference Hansen and Shafiei2022b), which includes serotonin (5HT1A (Savli et al., Reference Savli, Bauer, Mitterhauser, Ding, Hahn, Kroll, Neumeister, Haeusler, Ungersboeck, Henry, Isfahani, Rattay, Wadsak, Kasper and Lanzenberger2012), 5HT1B (Baldassarri et al., Reference Baldassarri, Hillmer, Anderson, Jatlow, Nabulsi, Labaree, Cosgrove, O’Malley, Eissenberg, Krishnan-Sarin and Esterlis2018; Gallezot et al., Reference Gallezot, Nabulsi, Neumeister, Planeta-Wilson, Williams, Singhal, Kim, Maguire, McCarthy, Frost, Huang, Ding and Carson2010; Matuskey et al., Reference Matuskey, Bhagwagar, Planeta, Pittman, Gallezot, Chen, Wanyiri, Najafzadeh, Ropchan, Geha, Huang, Potenza, Neumeister, Carson and Malison2014; Murrough et al., Reference Murrough, Czermak, Henry, Nabulsi, Gallezot, Gueorguieva, Planeta-Wilson, Krystal, Neumaier, Huang, Ding, Carson and Neumeister2011; Murrough et al., Reference Murrough, Henry, Hu, Gallezot, Planeta-Wilson, Neumaier and Neumeister2011; Saricicek et al., Reference Saricicek, Chen, Planeta, Ruf, Subramanyam, Maloney, Matuskey, Labaree, Deserno, Neumeister, Krystal, Gallezot, Huang, Carson and Bhagwagar2015; Savli et al., Reference Savli, Bauer, Mitterhauser, Ding, Hahn, Kroll, Neumeister, Haeusler, Ungersboeck, Henry, Isfahani, Rattay, Wadsak, Kasper and Lanzenberger2012), 5HT2A (Beliveau & Ganz, Reference Beliveau and Ganz2017, 5HT4 (Beliveau & Ganz, Reference Beliveau and Ganz2017, 5HT6 (Radhakrishnan et al., Reference Radhakrishnan, Nabulsi, Gaiser, Gallezot, Henry, Planeta, Lin, Ropchan, Williams, Morris, D’Souza, Huang, Carson and Matuskey2018; Radhakrishnan et al., Reference Radhakrishnan, Matuskey, Nabulsi, Gaiser, Gallezot, Henry, Planeta, Lin, Ropchan, Huang, Carson and D’Souza2020), 5HTT (Beliveau & Ganz, Reference Beliveau and Ganz2017), norepinephrine (α4β2 (Baldassarri et al., Reference Baldassarri, Hillmer, Anderson, Jatlow, Nabulsi, Labaree, Cosgrove, O’Malley, Eissenberg, Krishnan-Sarin and Esterlis2018; Hillmer et al., Reference Hillmer, Esterlis, Gallezot, Bois, Zheng, Nabulsi, Lin, Papke, Huang, Sabri, Carson and Cosgrove2016), M1 (Naganawa et al., Reference Naganawa, Nabulsi, Henry, Matuskey, Lin, Slieker, Schwarz, Kant, Jesudason, Ruley, Navarro, Gao, Ropchan, Labaree, Carson and Huang2021), VAChT (Aghourian, Legault-Denis, Soucy, & Rosa-Neto, Reference Aghourian, Legault-Denis, Soucy and Rosa-Neto2017; Bedard et al., Reference Bedard, Aghourian, Legault-Denis, Postuma, Soucy, Gagnon, Pelletier and Montplaisir2019), cannabinoid (CB1 (D’Souza et al., Reference D’Souza, Cortes-Briones, Ranganathan, Thurnauer, Creatura, Surti, Planeta, Neumeister, Pittman, Normandin, Kapinos, Ropchan, Huang, Carson and Skosnik2016, Neumeister et al., Reference Neumeister, Normandin, Murrough, Henry, Bailey, Luckenbaugh, Tuit, Zheng, Galatzer-Levy, Sinha, Carson, Potenza and Huang2012, Normandin et al., Reference Normandin, Zheng, Lin, Mason, Lin, Ropchan, Labaree, Henry, Williams, Carson, Neumeister and Huang2015, Ranganathan et al., Reference Ranganathan, Cortes-Briones, Radhakrishnan, Thurnauer, Planeta, Skosnik, Gao, Labaree, Neumeister, Pittman, Surti, Huang, Carson and D’Souza2016), dopamine (D1 (Kaller et al., Reference Kaller, Rullmann, Patt, Becker, Luthardt, Girbardt, Meyer, Werner, Barthel, Bresch, Fritz, Hesse and Sabri2017), D2 (Sandiego et al., Reference Sandiego, Gallezot, Lim, Ropchan, Lin, Gao, Morris and Cosgrove2015; Slifstein et al., Reference Slifstein, van de Giessen, Van Snellenberg, Thompson, Narendran, Gil, Hackett, Girgis, Ojeil, Moore, D’Souza, Malison, Huang, Lim, Nabulsi, Carson, Lieberman and Abi-Dargham2015; Smith et al., Reference Smith, Crawford, Dang, Seaman, San Juan, Vijay, Katz, Matuskey, Cowan, Morris, Zald and Samanez-Larkin2019; Zakiniaeiz et al., Reference Zakiniaeiz, Hillmer, Matuskey, Nabulsi, Ropchan, Mazure and Picciotto2019), DAT (Dukart et al., Reference Dukart, Holiga, Chatham, Hawkins, Forsyth, McMillan, Myers, Lingford-Hughes, Nutt, Merlo-Pich, Risterucci, Boak, Umbricht, Schobel, Liu, Mehta, Zelaya, Williams, Brown and Sambataro2018), GABA (GABAa (Nørgaard et al., Reference Nørgaard, Beliveau, Ganz, Svarer, Pinborg, Keller, Jensen, Greve and Knudsen2021), histamine (H3 (Gallezot et al., Reference Gallezot, Planeta, Nabulsi, Palumbo, Li, Liu, Rowinski, Chidsey, Labaree, Ropchan, Lin, Sawant-Basak, McCarthy, Schmidt, Huang and Carson2017), glutamate (mGluR5 (DuBois et al., Reference DuBois, Rousset, Rowley, Porras-Betancourt, Reader, Labbe, Massarweh, Soucy, Rosa-Neto and Kobayashi2016; Smart et al., Reference Smart, Cox, Scala, Tippler, Jaworska, Boivin, Séguin, Benkelfat and Leyton2019), NMDA (Galovic et al., Reference Galovic, Erlandsson, Fryer, Hong, Manavaki, Sari, Chetcuti, Thomas, Fisher, Sephton, Canales, Russell, Sander, Årstad, Aigbirhio, Groves, Duncan, Thielemans and Koepp2021; McGinnity et al., Reference McGinnity, Hammers, Riaño Barros, Luthra, Jones, Trigg, Micallef, Symms, Brooks, Koepp and Duncan2014), opioid (MOR (Kantonen et al., Reference Kantonen, Karjalainen, Isojärvi, Nuutila, Tuisku, Rinne, Hietala, Kaasinen, Kalliokoski, Scheinin, Hirvonen, Vehtari and Nummenmaa2020), and norepinephrine (NET (Belfort-DeAguiar et al., Reference Belfort-DeAguiar, Gallezot, Hwang, Elshafie, Yeckel, Chan, Carson, Ding and Sherwin2018; Ding et al., Reference Ding, Singhal, Planeta-Wilson, Gallezot, Nabulsi, Labaree, Ropchan, Henry, Williams, Carson, Neumeister and Malison2010; Li et al., Reference Li, Potenza, Lee, Planeta, Gallezot, Labaree, Henry, Nabulsi, Sinha, Ding, Carson and Neumeister2014; Sanchez-Rangel et al., Reference Sanchez-Rangel, Gallezot, Yeckel, Lam, Belfort-DeAguiar, Chen, Carson, Sherwin and Hwang2020). PET images were averaged within each study, normalized to the MNI-ICBM 152 non-linear 2009 template, parcellated into 200 brain regions, and Z-scored (Hansen & Shafiei, Reference Hansen and Shafiei2022a, b).
We constructed a multilinear model to assess the association between neurotransmitter receptors/transporters profiles and abnormal IVSC and IVFC patterns (unthresholded t-static vectors) separately. The significance of these multilinear models was evaluated through permutation testing (10000 times) with Bonferroni correction. To further understand the relative importance of each predictor (neurotransmitter receptor/transporter profile) in the model’s overall predictive power, we conducted dominance analysis for each multilinear model. This analysis evaluates the impact of each predictor on the model fit by considering all possible combinations of predictors, resulting in 2n – 1 subset models for n predictors (Budescu, Reference Budescu1993). We used total dominance value to gauge each of the predictor’s relative significance, calculated as the mean increase in explained variance (R2) across all subset models in which it is included (Budescu, Reference Budescu1993).
Reproducibility and disease specificity analysis
To assess the robustness and reproducibility of abnormal IVSC and IVFC patterns, we performed several validation analyses. First, we validated these patterns using an alternative brain atlas, the Automated Anatomical Labeling (AAL) atlas, to ensure results were not dependent on the specific parcellation scheme used (Tzourio-Mazoyer et al., Reference Tzourio-Mazoyer, Landeau, Papathanassiou, Crivello, Etard, Delcroix, Mazoyer and Joliot2002). Second, we employed a split-sample method to evaluate whether abnormal IVSC and IVFC patterns were consistent across two halves of the dataset. We randomly divided patients into two subgroups, computed abnormal patterns (unthresholded t-statistic vectors) for each subgroup separately, and calculated Spearman’s correlation coefficients to measure similarity. This process was repeated 100 times to minimize random selection errors. Third, we conducted leave-one-site-out validation to reduce the influence of specific sites. For each iteration, we excluded one site, computed abnormal IVSC and IVFC patterns using the remaining sites, and calculated Spearman’s correlation coefficients between the new patterns and those from all sites.
Additionally, to investigate whether abnormal IVSC and IVFC patterns are specific to MDD, we compared these patterns with those in schizophrenia using another independent dataset sourced from the Centre for Biomedical Research Excellence (COBRE, http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html) (Collin, Kahn, de Reus, Cahn, & van den Heuvel, Reference Collin, Kahn, de Reus, Cahn and van den Heuvel2014). The COBRE dataset includes 72 patients with schizophrenia and matched 74 HCs. Detailed information about the dataset and data acquisition protocols can be found elsewhere (Collin et al., Reference Collin, Kahn, de Reus, Cahn and van den Heuvel2014).
Results
Demographics and clinical characteristics
Demographics and clinical characteristics are detailed in Table 1.
Table 1. Demographic and clinical information of participants
Note: MDD, major depressive disorder; HCs, healthy controls; HAMD, Hamilton Depression Rating Scale; FD, framewise displacement.
a Two sample t test.
b Chi-square test.
Abnormal IVSC and IVFC patterns in MDD
Compared with HCs, patients with MDD demonstrated significantly higher IVSC in the precuneus and lingual gyrus but lower IVSC in the medial frontal gyrus, calcarine, cuneus, and cerebellum anterior lobe (p < 0.05 with Bonferroni correction). The abnormal pattern of IVSC is depicted in Figure 1 and Supplementary Table S1. Conversely, patients with MDD exhibited increased IVFC across nearly the entire brain (p < 0.05 with Bonferroni correction) including the middle frontal gyrus, anterior cingulate cortex, hippocampus, insula, striatum, and precuneus (Figure 1, Supplementary Table S1).
Figure 1. Abnormal patterns of intersubject variability of structural covariance (IVSC) and functional connectome (IVFC) in patients with MDD and the Spearman’s correlation between them.
Additionally, the spatial distribution of abnormal IVSC pattern (un-thresholded t-static vector) was significantly correlated with the distribution of abnormal IVFC pattern (Spearman’s r = −0.21, p = 0.003, see Figure 1).
More severe gray matter morphological atrophy is associated with higher IVSC
We observed significant spatial correlations between abnormal IVSC patterns and gray matter morphological atrophy, with a Spearman’s correlation coefficient of r = 0.22 (Bonferroni corrected p = 0.006, see Figure 2a). There was also a strong correlation between abnormal IVSC pattern and the regional structural covariance connectome (Spearman’s r = 0.49, Bonferroni corrected p < 0.001). Gray matter morphological atrophy showed a marginal correlation with an abnormal pattern of regional structural covariance connectome (Spearman’s r = 0.13, uncorrected p = 0.07, see Figure 2b).
Figure 2. Associations between abnormal IVSC and gray matter morphological atrophy and abnormal patterns of regional structural covariance connectome (SC differences). (A). The distribution of gray matter morphological atrophy (t statistics of healthy controls vs. patients with MDD) and its Spearman’s correlation with abnormal IVSC. (B). The distribution of SC differences and its Spearman’s correlation with abnormal IVSC. (C). Abnormal IVSC has a mediating effect between gray matter atrophy and abnormal IVFC.
Additionally, mediation analysis result showed that IVSC mediated the effects of gray matter atrophy pathway on IVFC in MDD (total effect (path C) = −0.48, 95% confidence interval (CI) = [−0.93, −0.02], p = 0.033; direct effect (path C′) = −0.34, 95% CI = [−0.56, −0.08], p = 0.143; indirect effect (path AB) = −0.15, 95% CI = [−0.35, −0.03], p = 0.016, Figure 2c).
Abnormal IVFC shows a close relationship with clinical features
No significant effects of episodicity (first episode versus recurrent) and medication status (treated versus untreated) were observed on IVSC (p > 0.05 with Bonferroni correction). Similarly, there were no significant correlations between abnormal IVSC and the total score of HAMD (p > 0.05 with Bonferroni correction).
For IVFC, recurrent patients demonstrated significantly greater variability (p < 0.05 with Bonferroni correction) than first-episode patients in the left striatum, insula, right lingual gyrus, posterior cingulate, and left calcarine (Figure 3). Treated patients showed lower variability (p < 0.05 with Bonferroni correction) than untreated patients in the right insula, superior temporal gyrus, and inferior parietal lobule (Figure 3). No significant correlations between abnormal IVFC and the total score of HAMD were observed (p > 0.05 with Bonferroni correction).
Figure 3. Group comparison results of IVFC between patient subgroups (p < 0.05 with Bonferroni correction).
There were no significant differences between male patients and female patients in either IVSC or IVFC (p > 0.05 with Bonferroni correction).
Association between abnormal IVSC/IVFC pattern and neurotransmitter receptors/transporters
We explored the associations between neurotransmitter receptors/transporters profiles and abnormal IVSC and IVFC patterns using multilinear models. The goodness-of-fit of these models (adjusted R2) was 0.42 (F-statistic (200,180) = 8.68) for IVSC and 0.25 (F-statistic (200,180) = 4.39) for IVFC (see Figure 4a). Permutation testing confirmed that these models performed significantly better than chance, with permutation p < 0.05 after Bonferroni correction (see Figure 4b).
Figure 4. Associations between distributions of neurotransmitter receptors/transporters and abnormal IVSC and IVFC. (a). The goodness of fits (adjusted R2) of multilinear models between distributions of neurotransmitter receptors/transporters and abnormal IVSC and IVFC. The significance is assessed through permutation testing. (b). Permutation testing results indicate that these performances of models (true R2) are significantly better than chance (permutation R2). (c). Dominance analysis is employed to determine the relative importance of predictors for each multilinear model. Total dominance values (d_total), indicative of the predictors’ relative importance, are depicted.
Dominance analysis identified key neurotransmitter receptors associated with each abnormal pattern. Specifically, MOR and CB1 were significant for abnormal IVSC, while GABAa was crucial for abnormal IVFC (see Figure 4c).
Reproducibility and disease specificity analysis
Reproducibility analysis confirmed the robustness and reproducibility of abnormal IVSC and IVFC patterns. Abnormal IVSC/IVFC patterns obtained using the AAL brain atlas were consistent with our main results. Patients with MDD demonstrated both increased and decreased IVSC and overall increased IVFC compared to HCs (see Figure 5a). Additionally, abnormal IVSC patterns remained significantly positively correlated with gray matter atrophy using the AAL brain atlas (see Figure 5a). Second, abnormal IVSC and IVFC patterns showed strong reproducibility between randomly split halves, with mean (± SD) Spearman’s correlations of 0.40 (±0.06) for IVSC and 0.78 (±0.05) for IVFC. The distributions of Spearman’s correlations are shown in Figure 5b. Third, leave-one-site-out validation also supported the consistency of our results, with a mean (± SD) Spearman’s correlations of 0.97 (±0.19) for IVSC and 0.99 (±0.15) for IVFC (see Figure 5c).
Figure 5. Reproducibility analysis results. (a). Abnormal patterns of IVSC and IVFC using AAL and Spearman’s correlation between them. (b). The distribution of Spearman’s correlations between randomly split halves (100 times). (c). Spearman’s correlations between abnormal patterns of IVSC and IVFC using some study sites and those using all sites.
Finally, we investigated abnormal IVSC and IVFC patterns in schizophrenia using an independent dataset from COBRE. Patients with schizophrenia demonstrated decreased IVSC and variable IVFC across brain regions compared to HCs (see Supplementary Figure S1). No significant association was found between abnormal IVSC and IVFC patterns in schizophrenia (Supplementary Figure S1). These results reinforced the MDD specificity of abnormal IVSC and IVFC patterns.
Discussion
This study represents the first systematic investigation of the heterogeneity in structural and functional connectomes in a large sample of patients with MDD. Our findings reveal that patients with MDD exhibit significantly higher IVSC in regions such as the precuneus and calcarine while showing lower IVSC in the cuneus and medial frontal gyrus. Conversely, these patients demonstrate increased IVFC across a broad range of brain regions, including the middle frontal gyrus, anterior cingulate cortex, hippocampus, insula, striatum, and precuneus. Notably, the spatial distribution of abnormal IVSC patterns was negatively correlated with that of abnormal IVFC patterns. Higher IVSC was associated with more severe gray matter morphological atrophy, suggesting that structural covariance connectome abnormalities might be compensatory in response to gray matter loss. Our study also highlights that abnormal IVSC and IVFC patterns are associated with distinct neurotransmitter receptor distributions, indicating different molecular underpinnings. Specifically, MOR and CB1 were linked to abnormal IVSC, while GABAa was associated with abnormal IVFC. This suggests that the variations in these neurotransmitter systems may contribute to the observed neuroimaging heterogeneity in MDD. Additionally, Additionally, our findings underscore the disease-specific nature and reproducibility of these abnormal patterns, which were confirmed through various validation methods. Together, these findings provide systematic insights into the heterogeneity of structural and functional connectomes, underlying observed inconsistent dysconnectivity findings in previous studies and underline the need for addressing the heterogeneity of abnormal structural and functional connectomes to guide individualized clinical diagnosis and treatment evaluations.
In line with the diverse clinical manifestations of psychiatric disorders, neuroimaging studies have consistently observed higher intersubject variability of neuroimaging metrics among psychiatric patients. For instance, patients with schizophrenia demonstrate higher gray matter morphological variations and IVFC compared to healthy controls (Cheng et al., Reference Cheng, Cai, Liu, Yang, Pan, Zhang, Mo, Yu and Zhu2025; Fang, Wen, Niu, Wan, & Zhang, Reference Fang, Wen, Niu, Wan and Zhang2022; Mo et al., Reference Mo, Zhao, Li, Cai, Song, Wang, Yu and Zhu2024; Sun et al., Reference Sun, Liu, Ma, Duan, Wang, Xu, Xu, Xu, Wang, Tang, He and Xia2021). Similarly, previous studies have reported elevated intersubject variability of multi-model connectomes in MDD (Gai et al., Reference Gai, Chu, Li, Guo, Ma, Shi, Che, Zhao, Dong, Li, Xie and Mao2024b; Han et al., Reference Han, Zheng, Li, Zhou, Jiang, Fang, Wei, Pang, Li, Zhang, Chen and Cheng2023; Hou et al., Reference Hou, Jiang, Li, Liu, Hou, Yin, Zhang, Zhang, Xie, Zhang, Kong and Yuan2023). However, these studies have some limitations such as small sample sizes and a focus on single-modal connectomes. Our study, leveraging a large multi-site dataset, represents a significant advancement by systematically investigating abnormal IVSC and IVFC in MDD. Our findings reveal that patients with MDD exhibited higher IVSC in the precuneus and lingual gyrus while showing lower IVSC in the medial frontal gyrus, calcarine, cuneus, and cerebellum anterior lobe. These structural and functional abnormalities are implicated in the pathology of MDD and correlate with diverse clinical manifestations (Gentili et al., Reference Gentili, Ricciardi, Gobbini, Santarelli, Haxby, Pietrini and Guazzelli2009; Smith, Pang, Servatius, Jiao, & Beck, Reference Smith, Pang, Servatius, Jiao and Beck2016; Zhou et al., Reference Zhou, Wang, Ma, Xie, Kang, Xu, Deng, Gong, Nie, Yao, Bu, Wang and Liu2024). Higher IVSC offers insights into previously contradictory findings related to structural and connectome abnormalities in MDD (Han et al., Reference Han, Cui, Zheng, Li, Zhou, Fang, Sheng, Wen, Liu, Wei, Chen, Chen, Cheng and Zhang2023b; Han, Zheng, Li, Zhou, et al., Reference Han, Zheng, Li, Zhou, Jiang, Fang, Wei, Pang, Li, Zhang, Chen and Cheng2023; Hu et al., Reference Hu, Cheng, Tang, Long, Huang, Li, Song, Song, Li, Yin and Chen2024; Wang et al., Reference Wang, Yang, Zhang, Dong, Yu, Yuan and Lei2024; Zeng et al., Reference Zeng, Shen, Liu, Fang, Liu and Hu2015). Notably, we found that higher IVSC correlated positively with gray matter atrophy, suggesting that structural covariance abnormalities might compensate for gray matter loss. It is important to note that lower IVSC in regions such as the medial frontal gyrus does not imply lower intersubject variability in all neuroimaging metrics. For example, despite lower IVSC in the medial frontal gyrus, schizophrenia patients have been reported to exhibit higher intersubject gray matter variability in this region compared to healthy controls (Alnæs et al., Reference Alnæs, Kaufmann, van der Meer, Córdova-Palomera, Rokicki, Moberget, Bettella, Agartz, Barch, Bertolino, Brandt, Cervenka, Djurovic, Doan, Eisenacher, Fatouros-Bergman, Flyckt, Di Giorgio, Haatveit and Karolinska Schizophrenia2019; Brugger & Howes, Reference Brugger and Howes2017; Fang et al., Reference Fang, Wen, Niu, Wan and Zhang2022).
On the functional side, patients with MDD show overall higher IVFC relative to healthy controls. This finding is consistent with a previous study identifying higher intersubject heterogeneity in functional connectome in MDD (Hou et al., Reference Hou, Jiang, Li, Liu, Hou, Yin, Zhang, Zhang, Xie, Zhang, Kong and Yuan2023). However, the authors failed to find significant results using the AAL brain atlas due to the limited sample size (Hou et al., Reference Hou, Jiang, Li, Liu, Hou, Yin, Zhang, Zhang, Xie, Zhang, Kong and Yuan2023). The elevated IVFC in MDD provides an explanation for the inconsistencies observed in prior studies regarding functional connectivity, such as alterations in the default mode network (Javaheripour et al., Reference Javaheripour, Li, Chand, Krug, Kircher, Dannlowski, Nenadic, Hamilton, Sacchet, Gotlib, Walter, Frodl, Grimm, Harrison, Wolf, Olbrich, van Wingen, Pezawas, Parker and Wagner2021; Scalabrini et al., Reference Scalabrini, Vai, Poletti, Damiani, Mucci, Colombo, Zanardi, Benedetti and Northoff2020; Yang et al., Reference Yang, Chen, Chen, Li, Li, Castellanos, Bai, Bo, Cao, Chang, Chen, Chen, Chen, Cheng, Cheng, Cui, Duan, Fang, Gong and Yan2021). The overall elevated IVFC in MDD highlights the need for strategies to address neuroimaging heterogeneity, such as subgroup or individualized analysis. Our mediation analysis indicates that abnormal IVSC mediates the relationship between gray matter atrophy and IVFC. Moreover, IVFC is significantly influenced by illness duration and antidepressant treatment, aligning with the progression of structural dysfunction and offering insights into the progressive nature of MDD (Han et al., Reference Han, Zheng, Li, Liu, Wang, Jiang, Wen, Zhou, Wei, Pang, Li, Zhang, Chen and Cheng2021; Han et al., Reference Han, Zheng, Li, Liu, Wang, Jiang, Wen, Zhou, Wei, Pang, Li, Zhang, Chen and Cheng2023; Morris et al., Reference Morris, Puri, Walker, Maes, Carvalho, Bortolasci, Walder and Berk2019; Verduijn et al., Reference Verduijn, Milaneschi, Schoevers, van Hemert, Beekman and Penninx2015; Wang et al., Reference Wang, Genon, Dong, Zhou, Li, Yu, Yuan, He, Qiu, Feng, Chen and Lei2023). Antidepressant treatment also reduces IVFC, consistent with the known effects of these medications on functional connectivity (Tassone, Nezhad, Demchenko, Rueda, & Bhat, Reference Tassone, Nezhad, Demchenko, Rueda and Bhat2024; Zhang et al., Reference Zhang, Rolls, Wang, Xie, Cheng and Feng2024).
Our study further reveals that abnormal IVSC and IVFC patterns are associated with specific neurotransmitter receptor distributions. Neurotransmitter receptors, which play a crucial role in synaptic plasticity and neural communication, are implicated in the pathophysiology of MDD (Shine, Reference Shine2019). Neurotransmitter dysfunction is thought to drive disorder-specific cortical abnormalities and cross-disorder similarities in psychiatry (Hansen & Shafiei, Reference Hansen and Shafiei2022a, b). The effectiveness of modern antipsychotic and antidepressant drugs depends on the selective manipulation of neurotransmitter function (Hansen & Shafiei, Reference Hansen and Shafiei2022a). This study found that neurotransmitter receptor/transporter distributions account for 42% and 25% of the variability in abnormal IVSC and IVFC, respectively. Dominance analysis highlights MOR and CB1 as significant for abnormal IVSC, while GABAa is crucial for abnormal IVFC. These receptors are involved in mood and reward regulation, with evidence suggesting dysregulation of the endogenous opioid system in depression (Jelen, Stone, Young, & Mehta, Reference Jelen, Stone, Young and Mehta2022). For instance, MOR knockout mice exhibit reduced food-anticipatory activity and diminished motivation (Kas et al., Reference Kas, van den Bos, Baars, Lubbers, Lesscher, Hillebrand, Schuller, Pintar and Spruijt2004; Papaleo, Kieffer, Tabarin, & Contarino, Reference Papaleo, Kieffer, Tabarin and Contarino2007). Although the pathological mechanisms underlying brain connectome alterations in psychiatric disorders remain poorly understood, this study further supports the involvement of receptors in abnormal intersubject variability of brain connectome, providing a new insight into the molecular bias of neuroimaging heterogeneity in MDD.
Despite these insights, our study has limitations. First, the cross-sectional design does not capture changes in IVSC and IVFC over time. Longitudinal studies are needed to validate our findings regarding disease progression. Second, the clinical information available for patients was limited, which restricted our ability to explore potential associations between other clinical features and abnormal IVSC and IVFC patterns in MDD. Third, the influence of comorbidity on abnormal IVSC and IVFC remains unexplored. E.g. whether patients who have comorbidities with other mental disorders would have higher IVSC/IVFC should be investigated in the future.
In conclusion, our study provides a comprehensive analysis of intersubject heterogeneity in multi-modal brain connectomes in MDD using a large, multi-site dataset. We identify distinct patterns of abnormal IVSC and IVFC, with abnormal IVSC positively correlated with gray matter atrophy and mediating the association between gray matter atrophy and abnormal IVFC. The reproducibility and disease specificity of these patterns underscore their relevance, while the association with neurotransmitter receptor distributions bridges neuroimaging heterogeneity and molecular mechanisms. These findings contribute to a deeper understanding of MDD and highlight the need for individualized approaches in clinical diagnosis and treatment.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0033291725000078.
Acknowledgments
The authors thank everyone who participated in this research.
Funding statement
This research study was supported by the Natural Science Foundation of China (81601467, 81871327, 62106229, 62476252), the Medical science and technology research project of Henan province (201701011, SBGJ202102103, SBGJ202101013, SBGJ202302068, S20240045), China Postdoctoral Science Foundation (2022 M712890), and Medical science and technology research project of Henan province (LHGJ20230217).
Competing interest
All authors declared no conflict of interest.
Open access
