Introduction
Patients with mood disorders, including bipolar disorder (BD) and major depressive disorder (MDD), exhibit overlapping dominant psychopathology, particularly in the clinical presentation of depressive episodes, which poses ongoing challenges for early diagnosis and intervention (Cassano et al., Reference Cassano, Rucci, Frank, Fagiolini, Dell'Osso, Shear and Kupfer2004; Kendler & Gardner, Reference Kendler and Gardner1998; Krystal, Reference Krystal2014; Kupfer, Frank, & Phillips, Reference Kupfer, Frank and Phillips2012; Peralta & Cuesta, Reference Peralta and Cuesta1998). Despite decades of effort, there is still a notable absence of definitive biomarkers for mood disorders. This is not completely surprising because studies have mainly used nosology that differentiates mood disorders according to clinical symptoms in the lack of any objective biomarker, although the clinical diagnostic system, Diagnostic and Statistical Manual of Mental Disorders, has already dramatically revolutionized the psychiatry field. While long proposed as distinguishable diagnostic categories, MDD and BD share substantial key attributes as implicated by complementary sources of evidence from neuroimaging, histological, and genetic studies (Brambilla, Perez, Barale, Schettini, & Soares, Reference Brambilla, Perez, Barale, Schettini and Soares2003; O'Donovan & Owen, Reference O'Donovan and Owen2016; Sheline, Reference Sheline2003; Wei et al., Reference Wei, Duan, Womer, Zhu, Yin, Cui and Jiang2020, Reference Wei, Womer, Sun, Zhu, Sun, Duan and Zhang2023; Xia et al., Reference Xia, Womer, Chang, Zhu, Zhou, Edmiston and Xu2019). This suggests a broader continuum between MDD and BD than previously assumed. Thus, it is imperative to comprehend the multidimensional neurobiological intermediate phenotypes underlying psychopathology of mood disorders, transcending traditional diagnostic categories for MDD and BD. Such understanding is essential for deciphering the diverse pathways that contribute to the heterogeneous clinical symptoms observed across mood disorders.
Recent advancements in characterizing distinct imaging or genetic biotypes have shown promise in revealing biological heterogeneity of mental illness (Marquand, Rezek, Buitelaar, & Beckmann, Reference Marquand, Rezek, Buitelaar and Beckmann2016; Sarrazin et al., Reference Sarrazin, Cachia, Hozer, McDonald, Emsell, Cannon and Hamdani2018; Zhang, Sweeney, Bishop, Gong, & Lui, Reference Zhang, Sweeney, Bishop, Gong and Lui2023). Subtyping studies have been advocated as a vital progression toward a more neurologically grounded understanding of heterogeneity in mood disorders (Drysdale et al., Reference Drysdale, Grosenick, Downar, Dunlop, Mansouri, Meng and Etkin2017; Sun et al., Reference Sun, Sun, Lu, Dong, Zhang, Wang and Wei2023). Numerous investigations have employed various clustering methods to detect neuroimaging-based biotypes in transdiagnostic or diagnosis-specific groups (Chang et al., Reference Chang, Womer, Gong, Chen, Tang, Feng and Zhang2021; Ge, Sassi, Yatham, & Frangou, Reference Ge, Sassi, Yatham and Frangou2022; Sun et al., Reference Sun, Sun, Lu, Dong, Zhang, Wang and Wei2023). Recent depression studying studies mainly characterizing the distinct emotional symptoms between subgroups (Drysdale et al., Reference Drysdale, Grosenick, Downar, Dunlop, Mansouri, Meng and Etkin2017; Han et al., Reference Han, Zheng, Li, Zhou, Jiang, Fang and Zhang2022). Ge et al., also identified brain maturational subtypes using neuroimaging profiling and characterized the differences in psychopathology and cognition behaviors among youth with mood and anxiety disorders (Ge et al., Reference Ge, Sassi, Yatham and Frangou2022). Remarkably, they also delineated the biopsychosocial context in which subgroup patients arise (Ge et al., Reference Ge, Sassi, Yatham and Frangou2022). Distinct neuropathological mechanisms may underlie heterogeneity in the presentation and progression of the clinical phenotype. As we mentioned earlier, assessing the validity of subtyping or clustering results necessitates gathering additional evidence from multidimensional biological data sources that extend beyond clustering algorithm selection and just clinical symptom differences (Chang et al., Reference Chang, Womer, Gong, Chen, Tang, Feng and Zhang2021; Zhang, Wang, & Zhang, Reference Zhang, Wang and Zhang2022). Furthermore, the extent to which genetic heterogeneity influences or interacts with phenotypic expression has barely been explored and individual-level variability, including environment, genetic, or other factors, may lead to different levels of disease liability.
Neuroimaging has yielded a profusion of potential biomarkers for mood disorders (Sheline, Reference Sheline2003). Magnetic resonance imaging studies have revealed the shared and distinct gray matter volume (GMV) reductions in bilateral anterior cingulate and medial frontal cortices, insula in MDD and BD (Drevets, Reference Drevets2000; Jiang et al., Reference Jiang, Wang, Jia, Sun, Kang, Zhou and Tang2021). A potential explanation for the intricate findings is the presence of distinct underlying pathophysiology, which resulted in varied patterns of gray matter volume (GMV) abnormalities among different subtypes of patients with mood disorders. Notably, case–control designs, which are supremely prevailing in psychiatry research, yet overlook inter-individual variances that play a critical role in mapping the heterogeneous disease phenotype (Lv et al., Reference Lv, Di Biase, Cash, Cocchi, Cropley, Klauser and Cetin-Karayumak2021). Studies benchmarked each individual scan in the context of normative age-related GMV trends and computed brain charts using individualized centile scores or deviations of neuroanatomical maps (Bethlehem et al., Reference Bethlehem, Seidlitz, White, Vogel, Anderson, Adamson and Areces-Gonzalez2022; Chen, Holmes, Zuo, & Dong, Reference Chen, Holmes, Zuo and Dong2021). Individual deviations from normative ranges in brain mappings have effectively revealed subgroups of patients with major psychiatric disorders in recent studies (Jiang et al., Reference Jiang, Wang, Zhou, Palaniyappan, Luo, Ji and Huang2023; Sun et al., Reference Sun, Sun, Lu, Dong, Zhang, Wang and Wei2023). However, the brain morphometric heterogeneity characterized by individual deviations was supposed to have distinct associations with cell-type specific functions and genetic susceptibility in mental illness (Di Biase et al., Reference Di Biase, Geaghan, Reay, Seidlitz, Weickert, Pebay and Zalesky2022; Wen et al., Reference Wen, Fu, Tosun, Veturi, Yang, Abdulkadir and Singh2022), which remain elusive in mood disorders.
Our current study intended to quantify the brain morphological heterogeneity in mood disorders by mapping region-level changes of gray matter volume (GMV) at the level of individual patients. We chart individual GMV deviations utilizing a normative modeling method that maps inter-individual variances in reference to the healthy control range (Marquand et al., Reference Marquand, Kia, Zabihi, Wolfers, Buitelaar and Beckmann2019), and use unsupervised clustering method to identify subtypes in patients with MDD and BD. The resulting subtypes were then characterized by using clinical behaviors, cell-specific transcriptomic profiles, and polygenic risk scores (PRS). To minimize the potential influence of medication confounders, we initially conducted our study in a discovery cohort comprising drug-naïve and drug-free patients and then validated the imaging-based subtyping in a replication cohort (Voineskos et al., Reference Voineskos, Mulsant, Dickie, Neufeld, Rothschild, Whyte and Lerch2020). It was hypothesized that neuroimaging-based individual deviation biomarkers coupling with multi-dimensional distinct features including clinical behaviors, and genetic risks, can comprehensively describe the underlying heterogeneity in mood disorders.
Methods and materials
Participants
We recruited 114 MDD patients, 60 BD patients and 404 healthy controls (HC). All the patients were recruited from the inpatient department of the Shenyang Mental Health Center and the outpatient clinic of the Department of Psychiatry of the First Affiliated Hospital of China Medical University, Shenyang, China. All the patients were diagnosed by two experienced psychiatrists using the Structured Clinical Interview for the Diagnostic and Statistical Manual of Mental Disorders, 4th Edition, Text Revision (DSM-IV-TR) in patients aged 18 years and older and the Schedule for Affective Disorders and Schizophrenia for School-Age Children-present and Lifetime Version (K-SADS-PL) in those patients under the age of 18 years. The patients were either drug-naïve or had been drug-free for more than two months from oral psychotropic medications before enrollment. The severity of depressive, anxious, and psychopathological symptoms was respectively assessed by using the 17-item Hamilton Depression Rating Scale (HAMD), 14-item the Hamilton Anxiety Rating Scale (HAMA), and Brief Psychiatric Rating Scale (BPRS). We also evaluated participants’ cognitive function performance by using Wisconsin Card Sorting Test (WCST). Patients and HC were excluded if they had any major medical condition, neurological disorder, and MRI contraindications (see online Supplementary Material). HC participants were recruited from the local community through advertisement. The participants over 18 years old signed a written consent form themselves. If the participants age were <18 years, their parental/legal guardian provided written informed consent. This study was approved by the Ethics Committee of the First affiliated Hospital of China Medical University (Shenyang, China; Approved number: 2015-27-2). To validate the imaging phenotypes, we involved 268 patients who received medication treatment as the replication cohort. These patients were also recruited and scanned at the First Affiliated Hospital of China Medical University. The details of patients in the replication cohort are described in the online Supplementary Material.
Image acquisition and MRI processing
Structural MRI scanning was conducted on a Signa HDx 3.0 T superconductive MRI system (GE Healthcare, Little Chalfont, UK) at the First Affiliated Hospital, China Medical University. The details of MRI scanning parameters are presented in online Supplementary Material. T1-weighted images were preprocessed by using the Computation Anatomy Toolbox (CAT) implemented in Statistical Parametric Mapping (SPM 12) for voxel-based morphometry (VBM) calculation (http://www.neuro.uni-jena.de/cat/) (Gaser & Dahnke, Reference Gaser and Dahnke2016). The quality of images was assessed by using the automated weighted average image quality rating (IQR). We applied the cut-off (> = 80%, > = grad B) to ensure high-quality images for analysis as low-quality images can lead to GM underestimations in preprocessing. Voxel wised gray matter density was calculated and obtained applying CAT12 default parameters. The detailed preprocessing information is presented in online Supplementary Material. The Desikan–Killiany (DK) parcellation atlas partitioned the cortex into 68 cortical regions was used as the regions of interest (ROI) template. ROI-wise gray matter density for each brain region was used for investigating between-group differences.
GMV normative deviation calculation and subtypes clustering
Quantile regression was used to obtain a normative range of regional GMV variation as a function of age and sex descripted in a previous study (Lv et al., Reference Lv, Di Biase, Cash, Cocchi, Cropley, Klauser and Cetin-Karayumak2021). We positioned individuals with MDD and BD on the normative percentile charts based on HC and expressed three kinds of continuous measurement of deviation from the generated normative range including the 5th percentile (z5) quantile regression predictor, the 50th percentile (z 50) quantile regression predictor and the 95th percentile (z95) quantile regression predictor, as individual deviation z-scores for each brain region, representing the difference from normative GMV calculated across all HC individuals. We then obtained individual-specific profiles of regional GMV deviations (z 50 maps) to perform subtyping analysis for all patients. Ward's linkage measurement and hierarchical clustering algorithm were used to identify clusters of patients based on the GMV deviations maps (Fig. 1). We assessed the stability of clustering based on scores of the adjusted rand index (ARI) for 2–5 clusters which were divided by using the hierarchical cutoffs (Hubert & Arabie, Reference Hubert and Arabie1985). We then chose the clustering solution with the highest averaged ARI following 10-fold cross-validation with 1000 times permutations, as the optimal clustering results. In this study, the 174 drug-free patients constitute the discovery dataset (dataset 1), while the 268 medicated patients comprise the validated dataset (dataset 2). We validate the distinct GMV deviations between subtypes in dataset 1 and assess the correlations of regional deviation findings between dataset 1 and dataset 2.
Figure 1. Schematic overview of the workflow in this study. (a) The neuroimaging data were from 174 drug-naïve and drug-free patients with MDD/BD and 404 participants as HC; we used quantile regression model and created normative curves with age based on GMV in HC group and calculated individual deviations for each region; patients were clustered in two subtypes using hierarchical clustering; we then validated the subtyping stability using ARI scores and examined GMV deviations reproducibility in a replication cohort; (b) we compared clinical symptoms and cognitive symptoms differences between subtypes; we also created symptoms network using HAMD, HAMA and BPRS factors and compared global and nodal network properties between subtype groups; (c) we then utilized AHBA brain-wide gene expression data to selected GMV deviations related genes, and identified their biological process and cell type components; and (d) we use PRS-AD and PRS-MD scores to compare AD genetic risk and MDD genetic risk between subtype groups.
To further analyze the significantly altered GMV deviations in patients, for each cortical region, we categorized patients as either: (1) within the HC's normative range of variation, labeled as normal; (2) significantly exceeding the HC's normative range, labeled as supra-normal; or (3) significantly below the HC's normative range, labeled as infra-normal. We used z scores to quantify individual GMV deviation from the 5% and 95%, which was then utilized to guide the above mutually exclusive classification. We acquired the standard deviations for z scores from the bootstrapped confidence intervals (CI). We defined the supra-normal as any individual exceeding the 95% CI for the 95th percentile (z 95 > 1.96), and infra-normal as any individual below the 5% CI for the 5th percentile (z 5 < − 1.96). Each patient was represented with a GMV deviation encoding map in which brain areas were numerically encoded with either −1 (infra-normal), 0 (normal), or +1 (supra-normal), and a GMV deviation z map with matched z 5 or z 95 values.
The average of the above encoded number across all patients in each subtype group produced a whole-brain summary measurement defined as the average abnormality rate, which quantified the overall percentage of GMV deviation. The average of the GMV deviation z map across all patients in each subtype group defined as the average abnormality extent (mean values of z 95 or z 5), which quantified the overall extent of GMV deviation.
Clinical information validation in subtype groups
We used two-sample t test to compare the between-subtype differences of clinical measures including total scores of HAMD, HAMA, and BPRS. We also computed the between-group differences in WCST performance among groups. We included five indices of WCST, i.e., correct response (CR), completed categories (CC), total errors (TE), perseverative errors (PE), and non-perseverative errors (NPE). The significant level was set as p FDR < 0.05, false discovery rate (FDR) correction – was used for controlling the false positive rate.
We created clinical symptoms network by using correlations between symptoms factors in each subtype and compared the network topological properties based on graph theory (Rydin et al., Reference Rydin, Milaneschi, Quax, Li, Bosch, Schoevers and Lamers2023; Ye et al., Reference Ye, Zalesky, Lv, Loi, Cetin-Karayumak, Rathi and Di Biase2021). Previous studies suggested psychotic symptoms are considered as the distinct subtype of affective disorders and frequently occur accordingly with depressive and anxious symptoms in all stages of illness (Tonna, De Panfilis, & Marchesi, Reference Tonna, De Panfilis and Marchesi2012). In this study, the overall nodes of the symptom network included 17 items of HAMD, 14 items of HAMA, and 5 factors of BPRS (Chang et al., Reference Chang, Womer, Gong, Chen, Tang, Feng and Zhang2021). We conducted paired t tests to examine differences in network strength ranges between subtypes. The details are presented in online Supplementary Materials.
Imaging transcriptomics and virtual histology analysis in subtype groups
Transcriptional data were acquired from the open access Allen Human Brain Atlas (AHBA) database (https://human.brain-map.org/). The AHBA dataset was preprocessed according to previous practical guide proposed by Arnatkevic et al. (Arnatkeviciute, Fulcher, & Fornito, Reference Arnatkeviciute, Fulcher and Fornito2019). The five steps of data preprocessing are shown in online Supplementary Materials. Then, we calculated a mean of all tissue samples in a brain area and obtained the matrix (34 regions × 10 027 genes) of gene expression values. We conducted the multivariate regression approach of partial least squares (PLS) to investigate the transcriptome-imaging associations (Li et al., Reference Li, Seidlitz, Suckling, Fan, Ji, Meng and Chen2021; Morgan et al., Reference Morgan, Seidlitz, Whitaker, Romero-Garcia, Clifton, Scarpazza and Donohoe2019; Romero-Garcia et al., Reference Romero-Garcia, Seidlitz, Whitaker, Morgan, Fonagy, Dolan and Bullmore2020). The details of GMV deviations related to gene expression analysis are presented in online Supplementary Materials. Then we obtained two PLS1 gene lists respectively for subtype 1 (PLS1-subtype 1) and subtype 2 (PLS1-subtype 2). Finally, we separately performed biological process enrichment for PLS1-subtype 1 and PLS1-subtype 2 using Metascape (https://metascape.org/) with false discovery rates correction (p FDR < 0.05) (Zhou et al., Reference Zhou, Zhou, Pache, Chang, Khodabakhshi, Tanaseichuk and Chanda2019).
Following the procedure of Seidlitz et al. (Reference Seidlitz, Nadig, Liu, Bethlehem, Vértes, Morgan and Clasen2020), we obtained corresponding gene sets of seven canonical cell classes (Habib et al., Reference Habib, Avraham-Davidi, Basu, Burks, Shekhar, Hofree and Regev2017; Lake et al., Reference Lake, Chen, Sos, Fan, Kaeser, Yung and Zhang2018; Li et al., Reference Li, Santpere, Imamura Kawasawa, Evgrafov, Gulden, Pochareddy and Sestan2018), including excitatory (Neuro.ex), and inhibitory neurons (Neuro.in), astrocytes (Astro), microglia (Micro), endothelial cells (Endo), oligodendrocyte precursors (OPC), oligodendrocytes (Oligo), and performed virtual histology analysis (Li et al., Reference Li, Seidlitz, Suckling, Fan, Ji, Meng and Chen2021; Zong et al., Reference Zong, Zhang, Li, Yao, Ma, Kang and Hu2023). The details of cell-type-specific gene lists selection are presented in online Supplementary Materials. We separately overlapped the gene lists of the seven cell types with the PLS1-subtype 1 and PLS1-subtype 2. First, we tested whether genes of PLS1-subtype 1 and PLS1-subtype 2 significantly overlapped with these cell-type specific genes. The p value of the ratio of overlapped gene number to the cell type gene number was calculated by a permutation test (with p FDR < 0.05) (Li et al., Reference Li, Seidlitz, Suckling, Fan, Ji, Meng and Chen2021). Second, we compared between-subtype differences in the number of overlapped genes of seven cell types with the PLS1-subtype 1 and PLS1-subtype 2, using chi-squared test. The significant level was set as false discovery rates correction p FDR < 0.05.
Polygenic risk scores analysis
Details about genotyping and quality control, imputation, and calculation of polygenic risk scores are shown in the online Supplementary Materials. Associations of polygenic risk scores (PRS) including PRS for MDD (PRS-MDD) and PRS for Alzheimer's disease (PRS-AD) were analyzed respectively using logistic regression in both subtype 1 and subtype 2. We calculated Nagelkerke's pseudo-R 2 as a measurement to represent the variance explained by the logistic regression model. We estimated PRS-MDD and PRS-AD at six different levels of p-value thresholds (ranging from 1.0e−06 to 0.1) in both subtypes 1 and 2 (Euesden, Lewis, & O'reilly, Reference Euesden, Lewis and O'reilly2015). The significance level of polygenic risks logistic regression model was set as p FDR < 0.05.
Results
Clinical and demographic data in the two identified subtypes
The clinical, demographic, and cognitive profiles based on clinical diagnosis of the dataset 1 are shown in online Supplementary Table S1. We used GMV deviations (z 50 maps) and hierarchical clustering method and identified two subtypes in the MDD and BD samples (n = 174), subtype 1 (74 MDD and 38 BD patients), and subtype 2 (40 MDD and 22 BD patients) (online Supplementary Table S2). The clustering stability analysis showed that two clustered subtypes had the highest ARI score among k = 2–5 clusters (online Supplementary Figure S3). The distribution of clinical diagnosis (MDD and BD) did not vary between subtype 1 and subtype 2 (x2 = 0.04, p = 0.83, online Supplementary Table S2). There were no differences in age, gender, or education between subtypes 1 and 2 (online Supplementary Table S2), indicating that the subtype discrimination was not by demographic characteristics. The detailed demographic information of subtypes 1 and 2 identified in the validation cohort are presented in online Supplementary Table S12.
GMV deviations differences between the two subtypes
In subtype 1, patients showed supra-normal deviations in frontal cortex (orbitofrontal and rostral middle frontal regions, cingulated cortex, and paracentral cortex) and infra-normal deviations in temporoparietal joint area, compared to HC (z 95 > 1.96, z 5 < −1.96, p < 0.05). In subtype 2, patients showed widespread infra-norm deviations in frontal, temporal, parietal, cingulate cortex, and supra-norm deviation in occipital cortex, compared to HC (z95 > 1.96, z5 < −1.96, p < 0.05) (Fig. 2A). In subtype 1, there were more patients labeled as supra-normal in more than half (57.4%) of the brain regions, particularly in prefrontal, cingulate, and paracentral cortex. In contrast, there were more patients labeled as infra-normal in subtype 2 in 90% of the cortical regions (Fig. 2B and 2C). The group averaged regional GMV deviations values (z 50 maps) were mapped and found to show significantly correlations between datasets 1 and 2 in both subtypes 1 (r = 0.29, p < 0.05) and 2 (r = 0.45, p < 0.0001) (online Supplementary Fig. S4). The top 10 regions with supra-normal deviations in frontal cortex in subtype 1 and regions with infra-normal deviations in frontal cortex and cingulated cortex in subtype 2 were found consistent in both datasets 1 and 2 (online Supplementary Table S9, Table S10).
Figure 2. Regions had significant GMV deviations and subject percentages with supra/infra norm regional deviations in subtypes 1 and 2. (a) regions with group mean GMV deviations values z95 > 1.96 or z5 < −1.96 were presented in both subtypes 1 and subtype 2; (b) we presented the supra and infra GMV deviations individual percentages for each region in both subtypes 1 and 2; (c) the details of supra and infra GMV deviations individual percentages in all regions were presented both subtypes 1 and 2; the red color: supra deviations; the blue color: infra deviations; the color bar: individual percentage.
Clinical and demographic profiles in the two subtypes
Of the two neuroimaging subtypes we identified, subtype 2 had significantly higher total scores of HAMD than that of subtype 1 (t = 2.67, p FDR = 0.025, Fig. 3A), while there were no significant between-subtype differences in the total scores of HAMA, BPRS, or illness duration (p > 0.05, online Supplementary Table S2). We further explored between-subtype differences in scores for each HAMD item and found that subtype 2 had more prominent guilt (t = 2.16, p = 0.032), early insomnia (t = 2.369, p = 0.019), insight (t = 2.16, p = 0.032), and general somatic symptoms (t = 2.219, p = 0.028) compared with that of subtype 1 (online Supplementary Fig. S6).
Figure 3. Clinical profiles in the two subtypes. (a) Subtype 2 had significantly higher total scores of HAMD than that of subtype 1; (b) subtype 1 had poorer WCST performance than that of HC; (c) we created symptom networks for subtype 1 and 2 based on graph theory; (d) subtype 2 had higher network connectivity strength relative to that of subtype 1.
In contrast, as to the between-group comparisons of cognitive function, subtype 1 had poorer performance in the WCST-PE (t = 4.00, p FDR = 0.001), WCST-NPE (t = 2.19, p FDR = 0.08), WCST-CR (t = −3.65, p FDR = 0.001), WCST-CC (t = −3.63, p FDR = 0.001), and WCST-TE (t = 3.61, p FDR = 0.001) compared with that of HC, while the WCST performance of subtype 2 did not differ from that of HC (Fig. 3B).
Furthermore, based on graph theory, we created symptom networks for subtypes 1 and 2 using all the items of HAMD, HAMA, BPRS-five factors, and found the global network density and strength were higher in subtype 2 (density = 0.51, strength = 6.76) comparing to that in subtype 1 (density = 0.42, strength = 5.24) (Fig. 3C). Compared with subtype 1, subtype 2 also had higher global network strength under the same network density threshold (t = 10.17, p = 1.0e−06, Fig. 3D) and nodal network strength and degree (online Supplementary Fig. S5).
Transcriptomics signatures differences between subtype groups
We utilized a multivariate PLS regression approach to clarify the transcriptional characteristics associated with changes in GMV deviations in subtypes 1 and 2. The PLS1-subtype 1 explained 26.38% of the covariance of GMV deviations across the whole cortex. The PLS1-subtype 1 demonstrated similar spatial gene expression patterns with the deviations map of subtype 1, with positive scores in the frontal, cingulate, precentral cortex, and negative scores in the temporal and occipital cortex (Fig. 4A). Moreover, we detected a significantly positive association between the PLS1-subtype 1 scores and changes of GMV deviations across the whole cortex (PLS1-subtype 1: r = 0.52, p = 0.0017). The PLS1-subtype 1 was enriched in genes involved in Alzheimer's disease, cell or blood morphogenesis, and cellular component pathways (p FDR < 0.05; Fig. 4c, online Supplementary Table S11). The intersection of endothelial cell-related genes with the PLS1-subtype 1 gene list was significantly larger than the intersection of endothelial cell-related genes with the PLS1-subtype 2 gene list. The intersection of oligodendrocytes-related genes with the PLS1-subtype 1 gene list was significantly larger than the intersection of oligodendrocytes-related genes with the PLS1-subtype 2 gene list. (Fig. 4e, online Supplementary Table S3, S4).
Figure 4. Transcriptomics features and virtual histology of GMV deviations differences in subtypes 1 and 2. (a) Group mean GMV deviations scores in subtype 1 and PLS1-subtype 1 scores in in left hemisphere; (b) group mean GMV deviations scores in subtype 2 and PLS1-subtype 2 scores in in left hemisphere; (c) top 17(p FDR < 0.05) GO biological process and KEGG pathways enriched using PLS1-subtype 1 genes; (d) top 17(p FDR < 0.05) GO biological process and KEGG pathways enriched using PLS1-subtype 2 genes; (e) the number of overlapped genes between PLS1-subtype 1 genes and cell-specific genes of seven cell types; the black star: significant different between subtypes 1 and 2 using chi-square test at p FDR < 0.05; the between-subtype differences of percentages of the overlapped genes in Endo and Oligo; the percentages of overlapped genes mean the number of overlapped genes divide the number of PLS1 genes; the black star: significant different between subtypes 1 and 2 using chi-square test at p FDR < 0.05; (f) the number of overlapped genes between PLS1-subtype 1 genes and cell-specific genes of seven cell types; the red star: significant enriched in Micro using permutation test at p FDR < 0.05; the distribution of permutation of the number of overlapped genes between random PLS1-subtype 2 genes and micro genes; red dash line present real number of overlapped genes.
The PLS1-subtype 2 demonstrated similar spatial gene expression patterns with the deviations map of subtype 2, with positive scores in the occipital cortex, and negative scores in the frontal, cingulate, and temporal cortex (Fig. 4b). Moreover, we detected a significantly positive association between the PLS1-subtype 2 scores and changes of GMV deviations across the whole cortex (PLS1-subtype 2: r = 0.48, p = 0.0039). The PLS1-subtype 2 was enriched in genes involved in ion transport, calcium signaling, trans synaptic signaling, and response to stress biological process (p FDR < 0.05; Fig. 4d). The PLS1-subtype 2 genes were significantly enriched in microglia-related genes (the overlapped genes number = 117, permutation p FDR = 0.005; Fig. 4F), while marginally enriched in inhibitory neurons-related genes (the overlapped genes number = 157, permutation pFDR = 0.057; online Supplementary Table S5).
Genetic risk differences between subtype groups
Four PRS-AD scores at thresholds of p values of 1.0e−06 (N SNPs = 38), 1.0e−03 (N SNPs = 626), 1.0e−02 (N SNPs = 3007), and 0.1 (N SNPs = 8797) indicated significant differences between subtype 1 (n = 19) and HC (n = 160) (p FDR < 0.05), explaining 4.2%, 5.9%, 4.1%, and 5.4%, respectively, of the variation in subtype 1. Compared to HC, we observed no significant difference in PRS-AD in subtype 2 (n = 29). The PRS-MDD score at thresholds of p values of 1.0e−03 (N SNPs = 1357) indicated significant differences between subtype 2 (n = 29) and HC (n = 160) (p FDR < 0.05), explaining 10.8% of the variation in subtype 2. Compared to HC, there was no significant difference of PRS-MDD in subtype 1. The fitted PRS-AD and PRS-MDD scores for subtypes 1 and 2 are presented in Fig. 5.
Figure 5. PRS-AD and PRS-MDD scores in subtypes 1 and 2. (a) PRS-AD model were significantly fitted in subtype 1 at SNP threshold p < 1.0e−06, p < 1.0e−03, p < 1.0e−02, p < 1.0e−01, while PRS-AD model was not fitted in subtype 2; (b) PRS-MDD model were significantly fitted in subtype 2 at SNP threshold p < 1.0e−03, while PRS-MDD model was not fitted in subtype1; red bar presents the PRS-AD or PRS-MDD fit model had a significant level of p FDR < 0.05.
Discussion
In this study, we first used a novel hierarchical clustering method based on GMV deviations to identify two subtypes in mood disorders, supra-normal (subtype 1) and infra-normal (subtype 2) dominant subtypes, which differed in clinical behaviors, cell-specific transcriptomic profiles, and PRS. The supra-normal dominant subtype had significantly increased GMV deviations in frontal, cingulate, primary motor cortex, significantly impaired cognitive function (executive function) in the WCST performance, and significantly higher genetic risk for Alzheimer's disease. Furthermore, the genes of which transcriptional levels showed spatial associations with variations of GMV deviations in the supra-normal subtype enriched in Alzheimer's disease pathways, as well as biological processes such as cell, blood morphogenesis, and cellular component. In addition, the supra-normal subtype-related genes were significantly enriched in endothelial cells. In contrast, infra-normal dominant subtype demonstrated significantly decreased GMV deviations in frontal, temporal, parietal, and cingulate cortex, significantly severe depressive symptoms, and significantly higher genetic vulnerability to MDD. The genes from the spatial correlation and virtual histology analyses of infra-normal subtype significantly enriched in ion transport, calcium signaling, trans synaptic signaling, and response to stress biological process and inhibitory neurons. Collectively, our findings indicated opposite brain developmental patterns as a core and distinctive characteristic across the mood disorder continuum. Furthermore, the supra-normal and infra-normal dominant subtypes, delineated by this feature exhibited distinct associations with clinical behaviors, cell-specific transcriptomic profiles, and genetic risks. These imaging phenotypes between subtypes replicated in an independent dataset. Our current results provide insight into the heterogeneity of mood disorders in biological and behavioral terms and indicate a potential advancement toward categorical subtyping and, potentially, the development of precise individualized treatment for patients with MDD and BD.
A primary finding from this study is the identification of the supra-normal dominant subtype that had significantly increased GMV deviations in frontal cortex. Although the reasons for this abnormality are currently uncertain, one potential explanation for this phenomenon is that it may be associated with an inflammatory response usually occurred in the early phase of mood disorder (Chen et al., Reference Chen, Hsu, Huang, Tsai, Tu and Bai2023; Gritti, Delvecchio, Ferro, Bressi, & Brambilla, Reference Gritti, Delvecchio, Ferro, Bressi and Brambilla2021; Qiu et al., Reference Qiu, Lui, Kuang, Huang, Li, Zhang and Gong2014). In the initial stage of inflammation, it is known that endothelial cells, regulating inflammation through the expression of adhesion molecules and the release of cytokines and chemokines (Pober & Sessa, Reference Pober and Sessa2015), can activate leukocytes and promote their recruitment to sites of inflammation, which could increase cortical morphology (Kaplanski, Reference Kaplanski2018). Interestingly, our current study revealed that the variations of GMV deviations in the supra-normal subtype were spatially associated with the transcriptional levels of genes, which were mainly enriched in endothelial and oligodendrocytes cells as well as Alzheimer disease pathways. There is growing evidence suggesting that both oligodendrocytes and endothelial cells play important roles in cognitive function and Alzheimer's disease (Desai et al., Reference Desai, Mastrangelo, Ryan, Sudol, Narrow and Bowers2010). Oligodendrocytes are responsible for myelination of axons, which is essential for proper neuronal communication and plasticity (Fields, Reference Fields2015; Frühbeis et al., Reference Frühbeis, Fröhlich, Kuo, Amphornrat, Thilemann, Saab and Nave2013). Recent studies have demonstrated that disruptions in oligodendrocyte function and myelination can lead to cognitive impairments (Henn et al., Reference Henn, Noureldein, Elzinga, Kim, Savelieff and Feldman2022). Similarly, endothelial cells play a crucial role in maintaining the integrity of the blood-brain barrier (BBB), which is important for regulating the exchange of substances between the blood and the brain (Zlokovic, Reference Zlokovic2011). Dysfunction of the BBB has been associated with various neurological disorders and has been implicated in the pathogenesis of neurocognitive disorders such as Alzheimer's disease (Sweeney, Sagare, & Zlokovic, Reference Sweeney, Sagare and Zlokovic2018). More importantly, the validation findings, presented in turn, demonstrated that patients in this supra-normal subtype had significantly higher genetic risk for Alzheimer's disease. These findings potentially explain that depression with cognitive impairment leading to a ‘pseudodementia’ presentation may serve as risk factor of Alzheimer's disease (Byers & Yaffe, Reference Byers and Yaffe2011; Green et al., Reference Green, Cupples, Kurz, Auerbach, Go, Sadovnick and Edeki2003; Kessing & Andersen, Reference Kessing and Andersen2004). Consistently, as to the between-group comparisons of cognitive function, this subtype had poorer performance in the executive cognition functions such as WCST.
In contrast, subtype 2 showed evidence of infra-normal development pattern involving almost globally lower GMV deviations. Our findings can be considered as potential evidence of developmental abnormalities in packing density and cell size which resonates with anomalies in cerebral maturation observed in neurodevelopmental disorders (Ge et al., Reference Ge, Sassi, Yatham and Frangou2022). In addition, the PLS1-subtype 2 genes were enriched in biological processes such as trans-synaptic signaling and regulation of cellular response to stress, and the genes from the virtual histology analyses of this subtype significantly enriched in microglia and inhibitory neurons, the most abundant neurons in GABAergic system. Microglia are key immune cells in the brain that play a role in regulating synaptic pruning and neuroinflammation (Inta, Lang, Borgwardt, Meyer-Lindenberg, & Gass, Reference Inta, Lang, Borgwardt, Meyer-Lindenberg and Gass2017). Recent studies have suggested that dysregulation of microglial function may contribute to the development of depressive symptoms (Calcia et al., Reference Calcia, Bonsall, Bloomfield, Selvaraj, Barichello and Howes2016). In addition, inhibitory neurons, which utilize the neurotransmitter GABA, are crucial in modulating the excitatory activity of neural circuits involved in mood regulation (Brambilla et al., Reference Brambilla, Perez, Barale, Schettini and Soares2003; Croarkin, Levinson, & Daskalakis, Reference Croarkin, Levinson and Daskalakis2011; Kalueff & Nutt, Reference Kalueff and Nutt2007). Dysregulation of GABAergic signaling has also been implicated in the pathophysiology of depression. Together, dysfunction of microglia, inhibitory neurons, and GABAergic signaling may contribute to the development of depressive symptoms. Consistently, our validation of PRS analysis demonstrated that patients in this subtype had significantly higher genetic risk for MDD. Of more importance, this subtype had higher total scores of depression severity than subtype 1. As to the comparison of each HAMD subitem, subtype 2 had more prominent guilt, early insomnia, insight, and general somatic symptoms. These findings indicated that our neuroimaging subtype delineation not only separates mildly and severely ill patients with mood disorder, but also identifies their subitem heterogeneity regarding guilt, early insomnia, general somatic, and insight features, which have long been considered in mood disorder research (Peralta & Cuesta, Reference Peralta and Cuesta1998; Peterson & Benca, Reference Peterson and Benca2006; Zhao et al., Reference Zhao, Wu, Zhang, Mellor, Ding, Wu and Peng2018). In addition, the symptom network analysis revealed that patients in subtype 2 had a higher network connectivity strength than subtype 1. The greater symptom network strength may reflect a pattern of denser symptom interactions possibly owing to the globally lower GMV deviations in the infra-normal patients. Findings in this study advance the clinical conceptualization by relating them to diverse regional changes in structural brain systems.
Mood disorders are heterogeneous diseases genetically, neurobiologically, and clinically, posing a considerable challenge for treatment and management of patients. Although current development of medication has paid attention to the complex clinical syndromes, advances in psychopharmacology for mood disorders have been restricted for decades. It may be important for nosology and treatment development that neurobiological and genetic heterogeneity in these syndromes would be further elucidated to establish targeted interventions for discrete subtypes similar to many fields in medicine. Our findings indicate a move forward in that direction. Although there has been a growing interest in the heterogeneity of mood disorders in recent years, most previous work has addressed this issue by sorting patients separately according to their behavioral, imaging, or genetic features. This approach may limit our ability to delineate the degree to which genetic heterogeneity impacts or interacts with neuroimaging and behavioral phenotypes. Interestingly, imaging transcriptomics analysis promoted the heuristic attempts and findings of cognitive impairment subtypes may have a genetic risk of AD and subtypes with more severe depressive symptoms may have genetic vulnerability to MDD. Uniting neuroimaging with data from genetics as well as clinical data and integrative computational strategies would enable imaging phenotypes subtyping to traverse the knowledge gap between genetic heterogeneity and clinical observations (Brennand, Landek-Salgado, & Sawa, Reference Brennand, Landek-Salgado and Sawa2014; Krystal, Reference Krystal2014). Thus, establishing multi-dimensional characteristics of biological heterogeneity may provide a more prospective approach for subtype identification in terms of treatment development and clinical utility.
Studies that use trans-diagnostic methods are emerging as converging evidence implicated core features across mood disorders, with an increasing focus on the neuroimaging biomarkers and genetic risks from a systems perspective. Our current findings defined two distinctive subtypes across traditional clinical diagnostic boundaries. In each subtype reported herein, both MDD and BD were represented. The mismatch between clinical diagnoses and subtypes may partially explain frequent inconsistent findings among previous studies according to clinical diagnosis. The constraints of present diagnostic framework are apparent (Clementz et al., Reference Clementz, Sweeney, Hamm, Ivleva, Ethridge, Pearlson and Tamminga2016; Ivleva et al., Reference Ivleva, Clementz, Dutcher, Arnold, Jeon-Slaughter, Aslan and Meda2017; Meda et al., Reference Meda, Clementz, Sweeney, Keshavan, Tamminga, Ivleva and Pearlson2016; Pearlson, Clementz, Sweeney, Keshavan, & Tamminga, Reference Pearlson, Clementz, Sweeney, Keshavan and Tamminga2016). Refining the present diagnostic system with corresponding objective biological measures (e.g. supra-normal and infra-normal GMV) would generate more biologically homogeneous categories, which play a critical role in the development of more personalized and effective treatments.
There are several limitations to this work. First, the use of gene expression profiles from the healthy human brain in AHBA to explain GMV deviations changes is limited to the extent that transcription in patients could be different from those in healthy brains. Additionally, the AHBA gene expression data including the right hemisphere of the brain for only two participants, restricts the representation of the spatial correlation between entire brain transcriptional level and GMV deviations. Future studies need to incorporate more extensive datasets to comprehensively explore imaging transcriptomic relationships. Second, the genetic data size in our study was small. Future work with a larger sample size of genetic data is desirable. Third, we only used WSCT scores to test cognitive functions between subtype groups. Some other cognitive batteries would be validated in future work. Fourth, we did not obtain specific data about illness severity of the first episode in medication-free patients and were not able to detect the associations between GMV alterations and prior depressive/anxious/psychotic symptoms.
In this study, mood disorders were characterized into two subtypes associated with brain morphology, behavioral phenotypes, and genetic risk profiles. Our findings provide vital links between MRI-derived phenotypes, spatial transcriptome, and genetic vulnerability, illustrating the feasibility of integrating multi-omics information across multi-dimensional biological scales. The identification and validation of the two subtypes provide a potential framework for future studies about the underlying neuropathological mechanisms and genetic heterogeneity as well as precise clinical care in mood disorders.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291724000886.
Data and code availability
Human gene expression data that support the findings of this study are available in the Allen Brain Atlas (‘Complete normalized microarray datasets’, https://human.brainmap.org/static/download). The code for gene expression preprocessing can be found in previous study and can be downloaded at https://github.com/BMHLab/AHBAprocessing. The probe-to-gene annotations were obtained by the Re-annotator toolkit (v1.0.0, https://sourceforge.net/projects/reannotator/). All data supporting the findings of this study are provided within the paper and its supplementary information. Additional datasets are available from the corresponding author upon request. All the brain imaging, clinical behavior, and genetic results and the analysis codes in this study are available at https://github.com/zhengjunjie1234/Mood_Disorders_subtypeing.git.
Acknowledgements
This study was supported by grants from National Science Fund for Distinguished Young Scholars (81725005), National Natural Science Foundation of China-Guangdong Joint Fund (NSFC-Guangdong Joint Fund) (U20A6005), National Natural Science Foundation of China (62176129), National Key Research and Development Program (2022YFC2405603) and China Postdoctoral Science Foundation (2022M721681).
Author contributions
Authors Junjie Zheng, Lili Tang, Huiling Guo, Pengfei Zhao, Xizhe Zhang, Yanqing Tang, and Fei Wang were involved in participants recruitment and data collection. Junjie Zheng and Lili Tang executed the neuroimaging and multi-omic analysis. Junjie Zheng, Fay Y. Womer, Fei Wang, and Xiaofen Zong wrote the first draft of the manuscript. Yanqing Tang and Fei Wang guided the study design. Fei Wang supervised the whole study and revised the manuscript. All the authors contributed to the final version of the paper.
Competing interests
The authors report no biomedical financial interests or potential conflicts of interest.
Open access
