Introduction
Today, while we know much more about the pathophysiology of schizophrenia than we did more than a century ago when Bleuler described the construct of “a group of schizophrenias”, the recovery rates of this illness have not changed much[Reference Jaaskelainen, Juola, Hirvonen, McGrath, Saha and Isohanni1–Reference Norman, MacDougall, Manchanda and Harricharan3]. As Bleuler suspected, converging evidence indicated schizophrenia to be a heterogeneous set, with no single biological underlying process that can be invoked to account for all of the diagnosed patients[Reference Seaton, Goldstein and Allen4]. Nevertheless, we do not yet know how to distinguish one group of schizophrenia from another in an objective manner to enable us to advance in our mechanistic enquiries and for treatment selection––a key step to improve recovery rates in this condition.
In recent years, many studies have been conducted to dissect the heterogeneity of schizophrenia using measures of symptomatology[Reference Dickinson, Pratt, Giangrande, Grunnagle, Orel and Weinberger5, Reference Talpalaru, Bhagwat, Devenyi, Lepage and Chakravarty6], neurocognition[Reference Geisler, Walton, Naylor, Roessner, Lim and Charles Schulz7–Reference Van Rheenen, Lewandowski, Tan, Ospina, Ongur and Neill9], genetics[Reference Di Biase, Geaghan, Reay, Seidlitz, Weickert and Pebay10] and neuroimaging [Reference Clementz, Sweeney, Hamm, Ivleva, Ethridge and Pearlson11–Reference Xiao, Liao, Long, Tao, Zhao and Luo14]. Among them, the use of brain morphology has been the most promising, and the most often replicated subtyping approach to date. Brain morphological features are more stable, with low within-subject fluctuations than symptom rating scores and functional readouts (fMRI, EEG). In addition, MRI-derived measures such as thickness and degree of cortical folding (gyrification) can be quantified objectively in an automatized manner with minimal manual intervention in the quantification process. Thus, brain structure can provide more stable and reliable clustering solutions. Among the various morphological indices, cortical thickness across multiple brain regions has been employed as a feature of interest in many subtyping studies[Reference Xiao, Liao, Long, Tao, Zhao and Luo14–Reference Sugihara, Oishi, Son, Kubota, Takahashi and Murai19], with 2-6 distinguishable subgroups reported across studies of schizophrenia[Reference Xiao, Liao, Long, Tao, Zhao and Luo14–Reference Sugihara, Oishi, Son, Kubota, Takahashi and Murai19]. While the number of subgroups varies among studies, likely due to inherent noise in clinical sampling and variations in analytical approaches, one consistent feature of these studies is the presence of a notable subgroup with widespread cortical thickness reduction (“impoverished cortex”) compared to healthy subjects[Reference Chand, Dwyer, Erus, Sotiras, Varol and Srinivasan12, Reference Xiao, Liao, Long, Tao, Zhao and Luo14] and other patients[Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13, Reference Pan, Pu, Chen, Huang, Cai and Tao20].
We have recently observed inappropriately high cortical glutamate levels in the subgroup with reduced thickness [Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13]. This raises the possibility of a persistent glutamatergic activation at the receptor level––affecting the N-methyl-d-aspartate acid (NMDA) and α-amino-3-hydroxy-5-methylisoxazole propionic acid (AMPA) receptors and voltage-gated calcium channels (VGCC)––leading to calcium overload, oxidative stress and excitotoxic damage to dendritic spines[Reference Hasbani, Schlief, Fisher and Goldberg21–Reference Plitman, Nakajima, de la Fuente-Sandoval, Gerretsen, Chakravarty and Kobylianskii23]. Excessive synaptic elimination with loss of dendritic spines is suspected to underwrite the MRI readout of reduced thickness[Reference Keifer, Hurt, Gutman, Keilholz, Gourley and Ressler24]. This account is consistent with progressive cortical thinning noted in the early phase of psychosis[Reference Hatton, Lagopoulos, Hermens, Scott, Hickie and Bennett25, Reference Buchy, Makowski, Malla, Joober and Lepage26]. Nevertheless, it is important to note that ~80% of patients with schizophrenia do not show deviations (i.e., infra-normal levels) in thickness patterns[Reference Di Biase, Geaghan, Reay, Seidlitz, Weickert and Pebay10, Reference Lv, Di Biase, Cash, Cocchi, Cropley and Klauser27], while some may show supra-normal changes suggestive of adaptive or compensatory response[Reference Guo, Palaniyappan, Liddle and Feng28, Reference Xiao, Lui, Deng, Yao, Zhang and Li29]. Leveraging the variance in thickness alone is unlikely to uncover the full spectrum of heterogeneity in schizophrenia.
Changes in synaptic density can reduce the MRI-based intensity of the grey matter, and affect MR-derived measures based on the location of the gray-white boundary (e.g., cortical thickness, volume). However, contour-based MRI measures such as the gyrification index (ratio of buried to outer cortex) are likely to be more stable in the face of grey matter loss. For example, models of progressive age-related brain atrophy estimate that while 75–164% increase in ventricular volume could occur in adults up to the age of 80, only 2.7% change may occur in the gyrification index. Furthermore, gyrification patterns emerge in utero, and change in response to intrauterine disruptions (e.g., hypoxia[Reference Quezada, Castillo-Melendez, Walker and Tolcos30] and associated excitotoxicity[Reference Kaindl, Favrais and Gressens31]) but remain detectable in later life [Reference Lv, Di Biase, Cash, Cocchi, Cropley and Klauser27]. A large body of evidence now points to the presence of widespread reduction in gyrification (hypogyria) in established cases of schizophrenia [Reference Nesvag, Schaer, Haukvik, Westlye, Rimol and Lange32–Reference Matsuda and Ohi34], in conjunction with poor treatment response patterns [Reference Nelson, White, Kraguljac and Lahti35] and other markers of aberrant neurodevelopment (e.g., Neurological Soft Signs [NSS]) [Reference Hirjak, Kubera, Wolf, Thomann, Hell and Seidl36, Reference Hirjak, Wolf, Paternoga, Kubera, Thomann and Stieltjes37]. Thus, when clustering patients with schizophrenia using neuroanatomical information from multiple brain regions, employing gyrification index as a feature could identify a predominantly “hypogyric” subgroup with developmental aberrations, while thickness could identify a “reduced thickness” subgroup with excitotoxic tissue loss.
One of the major goals of a subtyping exercise in schizophrenia is to test if distinguishable biological mechanisms account for the presence of discrete subgroups of patients. Recently, the genome-wide association studies in schizophrenia have implicated various glutamate receptors and downstream calcium signaling pathways in this illness––such as GRM3[Reference Harrison, Lyon, Sartorius, Burnet and Lane38], GRIN2A[Reference Zhao, Li, Shi, Tang, Chen and Liu39], GRIA1[40] and CACNA1C[Reference Dedic, Pohlmann, Richter, Mehta, Czamara and Metzger41, Reference Liu, Wu, Xia, Yao and Wang42] involved in the functional regulation of metabotropic glutamate receptors (mGluR), NMDA receptors, AMPA receptors, and VGCC. The glutamate receptor system and calcium signaling (regulated by CACNA1C) [Reference Kamijo, Ishii, Horigane, Suzuki, Ohkura and Nakai43–Reference Huang, Mo, Sun, Yu, Li and Wu46] are highly pertinent to the morphological changes in psychosis [Reference Plitman, Patel, Chung, Pipitone, Chavez and Reyes-Madrigal47–Reference Javitt, Zukin, Heresco-Levy and Umbricht50]. While the molecular mechanism influencing the generation of cortical thickness and gyrification are likely to be distinct [Reference Ribeiro, Ventura-Antunes, Gabi, Mota, Grinberg and Farfel51], genetic variations in glutamate transmission could likely affect the development of both cortical gyrification (by influencing neuronal migration and subplate apoptosis [Reference Matsugami, Tanemura, Mieda, Nakatomi, Yamada and Kondo52–Reference Liu, Liu, Yang, Wu, Zhang and Fang55]) and thickness (see Smith and Walsh for a review [Reference Smith and Walsh56]). Furthermore, the ionotropic glutamate receptors (NMDA/AMPA [Reference Guo and Ma57] and VGCC [Reference Michels, Ganjam, Martins, Schratt, Wohr and Schwarting58] play a key role in regulating the glutamate-mediated excitotoxicity. In this context, we hypothesized that the genetic risk variants related to this process will be over-represented in the subgroups with reduced gyrification and thickness, but not in other patients with preserved cortical anatomy. To test this, we studied the polymorphisms pertaining to three main glutamatergic postsynaptic receptors: the NMDA, AMPA, mGluR, and the VGCC.
Morphological features, especially regional thickness, can be affected by several secondary features associated with schizophrenia, e.g., cannabis use, urban living, migration [Reference Neilson, Bois, Gibson, Duff, Watson and Roberts59]. While polygenic risk for schizophrenia has a strong influence on the morphology [Reference Neilson, Shen, Cox, Clarke, Wigmore and Gibson60, Reference Cheng, Frei, van der Meer, Wang, O’Connell and Chu61], the within-group variance among patients appears to reflect these secondary disease factors [Reference Alnaes, Kaufmann, van der Meer, Cordova-Palomera, Rokicki and Moberget62]. As a result, for subgroup-based genetic associations, leveraging within-group variance in healthy groups (especially unaffected siblings) will be important to uncover the suspected genetic associations in morphologically-driven subgroups.
In the current study, we estimate MRI-derived cortical thickness and gyrification across multiple brain regions in 3 groups: patients with schizophrenia, unaffected siblings of patients, and unrelated healthy individuals. Leveraging the within-group heterogeneity of cortical features among healthy as well as clinically and genetically affected groups [Reference Pan, Pu, Chen, Huang, Cai and Tao20, Reference Fischl and Dale63], we identify clusters of patients that have shared morphological features (within-cluster similarity), and study their clinical, cognitive, and genetic profiles (GRM3, GRIN2A, GRIA1, and CACNA1C variants). As cross-sectional phenotyping is insufficient to understand prognostic relevance, we followed up a consenting subsample (43 patients, average of 13.4 months) longitudinally and report prognostic associations. As an exploratory analysis, we also investigated if the cluster membership interacts with genetic polymorphisms in a region-specific manner to affect thickness and gyrification index. To foreshadow the result, we (1) confirm the existence of a “reduced thickness” and a “preserved cortex” subgroups, (2) report for the first time a distinct, predominantly hypogyric subgroup and (3) identify specific glutamate-receptor variations that may influence these anatomical phenotypes in schizophrenia.
Methods and Materials
Participants
Patients (SZ) (n = 114) with a diagnosis of schizophrenia (based on DSM-5 [Reference First64]) and their siblings (n = 42) were recruited from the inpatient and outpatient units at Second Xiangya Hospital of Central South University, Changsha, China from 2016 to 2021. The inclusion criteria of patients: (1) meet the DSM-5 diagnostic criteria for schizophrenia, and the diagnosis was later rechecked after 6 months through face-to-face or real-time video interview; (2) 12< age ≤ 35 years; (3) right-handed; and (4) normal intellectual development. The exclusion criteria included: (1) meet the diagnosis of any mental disorder(s) except schizophrenia in DSM-5; (2) any reported history of substance use, neurological disorder, or serious physical illness in themselves or their first-degree relatives; (3) any contraindication for MRI; (4) left-handedness (as China has a usually low prevalence of left-handedness, and exclusion was more practical than case-control matching) [Reference Kushner65]; (5) history of brain injury or conscious coma; (Reference Talpalaru, Bhagwat, Devenyi, Lepage and Chakravarty6) intellectual disability (IQ<70) and (Reference Geisler, Walton, Naylor, Roessner, Lim and Charles Schulz7) previous electroconvulsive therapy.
We also assessed the longitudinal change in symptoms and functional recovery in a subset of patients (n = 43) who were seeking help in a symptomatic state (81.4% with first-episode of illness, 4.7% for relapse, and 13.9% for persistent symptom burden after first episode) and were fully concordant with follow-up and received continuous antipsychotic treatment based on clinical records for 1-2 years. The average follow-up period was 13.4(±11.5) months. On initial presentation, 58.14% of the follow-up cohort had <1 month of lifetime exposure to antipsychotics; 25.58% of the patients had 3 to 12 months exposure; 16.28% of the patients have been treated for 12 to 30 months. This cohort was primarily recruited to assess subtype-specific differences in clinical and functional improvement rates.
In addition to patients and their siblings, healthy controls (HC) (n = 112) were recruited from the communities and schools in Changsha City. The inclusion criteria of HC are as below: (1) not meeting any diagnostic criteria for any mental disorders; (2) 12< age ≤ 35 years; (3) right-handed; and (4) normal intellectual development. The exclusion criteria were consistent with the criteria of schizophrenia patients, except for diagnosis. The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008. All participants gave written informed consent to the study approved by the local Ethics Committee of Second Xiangya Hospital.
Targeted gene selection and sequencing
In this study, we selected the target regions in the genome to focus on the genes relevant to the receptors and channels on glutamate postsynaptic membrane. Briefly, we included GRM3 (involved in metabotropic glutamate receptors), GRIN2A (involved in NMDA ionotropic glutamate receptors), GRIA1 (involved in KA/AMPA ionotropic glutamate receptors), and CACNA1C (involved in voltage-gated calcium channel). For each of these genes, both coding and non-coding (regulatory) regions were included in the sequencing target. The regulatory genomic regions were comprised of 5’ un-translated region (5’ UTR), 3’ untranslated region (3’ UTR), and intron-exon boundaries (25bp). Custom capture oligos were designed using SureDesign website of Agilent Technologies (Santa Clara, CA)(https://earray.chem.agilent.com/suredesign/).
Blood samples of participants were collected on the day of MRI scan, and then genomic DNA was extracted for sequencing. Genomic DNA (2 μg) was used to target enrichment and to construct a DNA library before targeted sequencing. The genomic DNA was sheared to an average size of 250bp by using of Covaris S220 (Covaris, the USA), and the DNA library preparation and the capture procedure were then performed by using the SureSelect XT Target Enrichment System (Agilent, the USA), following the manual strictly. For all DNA libraries, the Illumina Hiseq2000 sequencing system (Illumina, San Diego, CA, USA) was employed to generate the paired-end 150bp reads raw data. Each sample was sequenced to >80% coverage at a minimum of 30-fold read depth. The Annovar program (dated 2016-02-01) was used for single nucleotide variation (SNV) annotation. Any SNV recorded in dbSNP147, with a minor allele frequency of ≥1% in 1000 genome database, ≥1% in our dataset, and with missing calls in <10% of subjects were considered as single nucleotide polymorphisms (SNP) and included for subsequent individual-variant association analysis (SNPs failing the Hardy-Weinberg equilibrium test at a significance level of 0.0001 were removed). The above bioinformatics analysis was described in Supplementary Material A.
General data collection, cognitive test, and clinical assessment
For all participants, the general intelligence level was evaluated through Wechsler Adult Intelligence Scale (WAIS)-Digital symbol test[Reference Gong66], WAIS-Digit span test (Forward), WAIS-Digit span test (Backward), and WAIS-Arithmetic test. The social functioning level was evaluated through the Social and Occupational Functioning Assessment Scale (SOFAS) [Reference Goldman, Skodol and Lave67]. In terms of cognition, we adopted story retelling and N-back test[Reference Owen, McMillan, Laird and Bullmore68] for memory function assessment, Verbal-Fluency test[Reference Strauss, Sherman, Spreen and Spreen69] for language function assessment, a visual pattern test[Reference Della Sala, Gray, Baddeley, Allamano and Wilson70] for visual perception assessment, and a Wisconsin card sorting test[Reference Heaton, Chelune, Talley, Kay and Curtiss71] for overall cognitive assessment.
For schizophrenia patients, the diagnoses were made by qualified psychiatrists according to DSM-5 criteria. On the same day as the MRI session, the severity of symptoms was evaluated through the Positive and Negative Syndrome Scale (PANSS) [Reference Kay, Fiszbein and Opler72], the Scale for The Assessment of Positive Symptoms (SAPS) [Reference Phillips, Xiong, Wang, Gao, Wang and Zhang73], the Scale for The Assessment of Negative Symptoms (SANS) [Reference Phillips, Xiong, Wang, Gao, Wang and Zhang73], and the Schizophrenia Suicide Risk Scale (SSRS) [Reference Taiminen, Huttunen, Heila, Henriksson, Isometsa and Kahkonen74]. The duration of illness, antipsychotic load (converted into chlorpromazine equivalent per day), and duration of antipsychotic medication were recorded. The PANSS, SAPS, SANS, and SOFAS of patients was assessed after at least 2 months of antipsychotic treatment. The rate of reduction in the scores of PANSS, SAPS, and SANS was calculated as (Baseline-Follow up)/ Baseline, and the improvement rate of SOFAS was calculated as (Follow up-Baseline)/Follow up. Thus, a positive value indicates a better outcome over time in both cases.
Magnetic resonance image acquisition and image processing
The participants were scanned using a Siemens 3.0 Tesla MRI scanner at Second Xiangya Hospital of Central South University at Changsha, China. T1-weighted magnetic resonance imaging data were acquired using a three-dimensional spoiled gradient echo (SPGR) pulse sequence from the sagittal plane, scanning parameter as follow: TR=7.6 ms, TE=3.7 ms,FA =8°, 180 slices, matrix =256*200,and the field of view (FOV)=256×256 mm2,slices were contiguous with a slice thickness of 2 mm. Importantly, during the T1-weighted image acquisition, participants were asked to remain still, and if any motion-related artifacts were detected, the scans were repeated.
A surface-based approach using Free-Surfer (http://surfer.nmr.harvard.edu, version 7.1.1) was used to calculate the cortical thickness and gyrification in the whole brain. Following skull-stripping and intensity correction, the gray–white matter boundary for each cortical hemisphere was determined by tissue intensity and neighborhood constraints. The resulting surface boundary was tessellated to generate multiple vertices across the whole brain before inflating. Using a deformable surface algorithm guided by the gray–CSF intensity gradient, the resulting gray–white interface was expanded to create the pial surface. The inflated surface was then morphed into a sphere followed by registration to an average spherical surface for optimal sulcogyral alignment. Then, the vertex-wise method (advocated by Schaer et al. [Reference Schaer, Cuadra, Tamarit, Lazeyras, Eliez and Thiran75]) was used to continuously assess local gyrification index (LGI) of the entire cortex. This method is an extension of classical two-dimensional GI measurement that calculates the ratio of the pial perimeter over the outer perimeter on coronal sections [Reference Zilles, Armstrong, Schleicher and Kretschmann76]. It provides an LGI for each vertex on cortical surface, which reflects the amount of cortex buried in its immediate locality. After the above procedures, Desikan-Killiany Atlas (68 regions) was used to extract cortical thickness and gyrification of each region using the FreeSurfer software [Reference Desikan, Ségonne, Fischl, Quinn, Dickerson and Blacker77]. Topological defects were corrected manually by two members of the research staff via tktools of freesurfer (https://surfer.nmr.mgh.harvard.edu /fswiki/FsTutorial/TopologicalDefect_tktools).
Statistical analysis
Clustering analysis was conducted in Matlab platform (version x). Before the clustering, the cortical thickness and gyrification of 68 regions were transformed to Z-scores. Then we used the clustering based on Gaussian Mixture Model (GMM) and GAP statistics to identify clusters of participants who shared similar patterns of cortical thickness and gyrification. Gaussian clustering was applied to all participants, including schizophrenia patients, unaffected siblings and HCs. We set cluster numbers from 1 to 6 (6 was the maximum value of the image structural subtype found so far) and GAP statistics to estimate the optimal number of clusters in our data. Then we chose the smallest cluster number that conformed to Gap(k) ≥ Gap(k+1) − Sk+1 as the solution of cluster analysis based on the 1-standard error method suggested by Tibshirani [Reference Tibshirani and Hastie78].
One-way ANOVA (in SPSS 20.0) was used to compare morphological, clinical, demographic, and cognitive indices, as we expected patients to differ from siblings as well as healthy subjects in these phenotypes, with FDR correction to address inflated type 1 error. For data with non-normal distribution (e.g., percentile data on the accuracy of N-back), we used nonparametric Kruskal–Wallis test for statistical analysis. Chi-square analysis was applied for genetic analysis comparing patients and healthy controls. At last, a multivariate generalized linear model with the subgroup based on clusters as the fixed factor was used to test the effect size of all factors including morphological data and phenotypic characteristics.
Results
Demographic, genetic, cognitive, and clinical characteristics of all participants
The sequencing of target gene was conducted in a total of 299 participants (135 SZ, 122 HC, and 42 Siblings), of which 31 participants failed to complete clinical data acquisition and MRI scanning, due to withdrawal of consent for specific procedures or contraindications to or poor quality of MRI acquisition. Finally, 268 participants (114 SZ, 112 HC, and 42 Siblings) completed gene sequencing, general information recording, clinical and cognitive assessment, and MRI scanning (Table 1, with a subset of 43 patients providing follow-up symptom and social functioning scores for prognostic assessment). Significant differences were found in age (P < 0.001) among the three groups. As expected, patients showed significant impairment in all cognitive tasks and social function compared with HC (Table 1). In addition, the siblings also showed significant impairment in all cognitive tasks and social function compared with HC, but had relatively good social function (P < 0.001) and better performance in WCST (P = 0.02 for WCST correct, 0.01 for WCST error) and verbal-fluency test (P = 0.001) compared with schizophrenia.
General information, cognitive performance, genetic information, and clinical characteristics of the participants.
As shown in Table 1, association analysis in target genes showed 12 candidate SNPs were different between schizophrenia patients and healthy controls. These SNPs were located in GRIN2A, GRM3, and CACNA1C, but contrary to our expectations, GRIA1 did not differ between patients and control subjects. The differences between the siblings and HC were seen in GRM3 and CACNA1C, and the differences between the siblings and schizophrenia were limited to CACNA1C. Within 12 SNPs, there were two complete linkages between multiple loci (Link1_rs10491965_10774044_12578811 and Link2_rs2239078_ 74059849). Thus, we obtained 9 candidate SNVs with a diagnostic effect for further investigation.
Across the three groups, a general linear model (GLM) with age and sex as covariates showed no significant difference in gyrification among three groups, with minimal differences in cortical thickness, supporting our expectation of morphometric patterns being driven by smaller subgroups of patients (Supplementary Table S1).
Gaussian clustering and GAP statistics
We explored the possibility of the existence of one to six clusters by gmdistribution (GMM) clustering in Matlab platform, and identified the ideal cluster solution based on GAP statistics to be 3 subgroups (Figure 1A). The numbers of participants composing each subtype were respectively 108 (NSZ/HC/Sb=37/50/21), 54(NSZ/HC/Sb=30/16/8), and 106 (NSZ/HC/Sb= 47/46/13). There are significant differences in the proportion of patients among the three subgroups (P = 0.032). Subgroup 2 (n = 54) which was the smallest of all three subgroups, had a higher (P = 0.009) proportion of patients (55.6% (30/54), n = 30) than subgroup 1 (34.3% (37/108), n = 37), with subgroup 3 (44.3% (47/106), n = 47) being intermediate (Figure 1B). Subgroup 2 also had the lowest proportion of healthy control subjects (Figure 1B). The proportion of siblings in the three subgroups, respectively were 19.4% (n = 21), 14.8% (n = 8), and 12.3% (n = 13), with most siblings clustering together as part of subgroup 1 (50% of all siblings). Most healthy controls also clustered as part of subgroup 1 (41% of all healthy controls).
GAP statistics of the clustering by GMM and the composition of each subgroup. (A) GAP statistics when the cluster number was set from 1 to 6; (B) Percentage of components in each of the three subgroups. Notes: GMM, Gaussian mixture model; HC, healthy controls.
When comparing patients based on their subgroup membership (Supplementary Table S2), those in subgroup 1 had notable hypo-gyrification in widespread regions compared to subgroups 2 and 3, but there were no significant differences between subgroup 2 and 3 (i.e., gyrification trend: subgroup 1< subgroup 3= subgroup 2). With respect to cortical thickness, subgroup 2 had the highest, subgroup 3 the lowest, and subgroup 1 intermediate values across widespread cortical regions (i.e., thickness trend: subgroup 3< subgroup 1< subgroup 2) (Supplementary Table S2).
When comparing patients based on their subgroup membership against all healthy controls as one group, patients from subgroup 1 displayed a generalized “hypogyric” pattern of widespread reduction in gyrification compared with HCs. Subgroup 2 patients had “supra-normal” pattern, characterized by regional increase in cortical thickness and gyrification. The subgroup 3 patients displayed an “impoverished cortex” pattern with regional cortical thinning compared with HCs (Figure 2A).
Characteristics of three patient clusters in morphology, clinic, cognition and candidate SNVs. (A) Differences between each subgroup and healthy controls in regional cortical thickness and gyrification; (B) Differences between each subgroup and healthy controls in candidate SNPs; (C) Phenotypic characteristics of each subgroup. Notes: C1, patient cluster in subgroup 1; C2, patient cluster in subgroup 2; C3, patient cluster in subgroup 3; HC, healthy controls; WCST, Wisconsin card sorting test.
Characteristics of the three patient subgroups in clinical symptoms and cognition
Among the three morphological patient subgroups, there were no significant differences in age and sex. There was no significant difference among the three clusters in estimated total intracranial volume (eTIV), SOFAS, WAIS-Digital symbol test, WAIS-Digit span test (Forward), WAIS-Digit span test (Backward), WAIS-Arithmetic test, N-back, story retelling, onset age, DUP, DoI, DoM, DDD, total PANSS, SAPS, and SRSS, after controlling for the effect of age and sex. However, SANS scores (P = 0.014; hypogyric > supra-normal P = 0.004) and other cognitive tasks, including verbal-fluency test (P = 0.02, hypogyric < supra-normal P = 0.005), visual pattern test (total score, P = 0.001), Wisconsin Card Sorting Test (WCST correct responses, P = 0.013 and errors, P = 0.011) varied as per subgroup membership (Figure 2C). Compared to the “hypogyric” subgroup (cluster 1) and “impoverished cortex” subgroups (cluster 3), the “supra-normal” subgroup (cluster 2) also had better visual pattern test performance (P = 0.002, 0.001), and WCST performance (more correct (P = 0.004, 0.027) and less error (P = 0.002, 0.001) responses).
In the patients with prognostic data (n = 43), a GLM with age and sex as covariates did not find significant effects on follow-up time, symptom reduction rate (PANSS, SAPS, SANS) and the change of all cognitive tasks, but there was a significant difference (P = 0.01) in social function recovery (improvement rate in SOFAS) among the three clusters. The functional recovery of subgroup 3 (“impoverished cortex,” SOFAS change: mean%[SD]= –62%[119], negative sign indicating deterioration over time) was significantly worse than that of the cluster 1 (“hypogyric,” SOFAS improvement = 5%±30%) (P = 0.017) and cluster 2 (“supra-normal,” SOFAS improvement = 13%±25%) (P = 0.004) (Figure 3). Taken together, these results indicated that subgroup 1 was the most hypogyric, cognitively affected group with higher negative symptom burden; subgroup 3 had the most reduction in thickness, with most sluggish functional recovery patterns while the subgroup 2 with supra-normal morphometric indices had the best cognitive, symptomatic and functional recovery patterns.
Difference in improvement rate of SOFAS among three patient subgroups. C1, patient cluster in subgroup 1; C2, patient cluster in subgroup 2; C3, patient cluster in subgroup 3.
Differences in candidate SNVs between the subgroups
In terms of genetic differences, five SNVs within nine candidate SNVs were identified by the comparisons between each heterogeneous cluster and healthy controls (Figure 2B). The prevalence of variants of rs9940680 (P = 0.009) and rs145139281 (P = 0.033) were higher in the hypogyric subgroup compared to healthy controls. The rate of prevalence of risk variants in the supra-normal subgroup was not significantly different compared to healthy controls. The prevalence of variants in the impoverished cortex subgroup (cluster 3) was significantly lower in rs223906 (P = 0.033) and rs556844413 (P = 0.015) compared to healthy controls (Figure 2B). These 5 SNVs with subgroup-level differences were included as candidate SNVs in the following analysis exploring the main effect on morphology (Figure 4) and the interaction between SNVs and morphology (Supplementary Figure S2).
The significant effects of candidate SNVs with differences between subgroups on cortical thickness and gyrification in each significant region. lgi: local gyrification index.
The effect of candidate SNVs with differences between the subgroups on cortical thickness and gyrification
In all participants, Glm with age and sex as covariates showed that five candidate SNPs with effects on heterogeneous groups (Figure 2B) were also associated with regional cortical thickness and gyrification. Concretely, on cortical gyrification, the effects of rs145139281 (in GRM3) were observed in bilateral lingual (P FDR = 0.045, 0.027), bilateral pericalcarine (P FDR = 0.027, 0.036) and left paracentral (P FDR = 0.045); the effects of rs9940680 and rs1420040 (in GRIN2A) were respectively observed in parstriangularis (P FDR = 0.027, 0.018) (Figure 4). On cortical thickness, the effect of rs55684413 (in CACNA1C) was observed in left pericalcarine (P FDR = 0.042); the effect of rs2239016 (in CACNA1C) was observed in right temporal pole (P FDR = 0.018); the effects of rs9940680 and rs1420040 (in GRIN2A) were also respectively observed in left pericalcarine (P FDR = 0.011, 0.02)(Figure 4).
Discussion
Using a data-driven approach that leverages normal, illness-related as well as shared genetic variance in cortical thickness and gyrification in a relatively early stage (<3 years of illness) of schizophrenia, we report three major findings: (1) three morphologically distinguishable clusters of patients are identifiable (reduced gyrification, reduced thickness, and higher thickness subgroups); (2) phenotypically the subgroups are mostly similar but have the varying burden of negative symptoms, cognitive deficits, and functional deterioration, and (3) the subgroup membership influences illness-related variation in glutamatergic receptor polymorphisms. In addition, the unaffected siblings were more often (50% of siblings) clustered with hypogyric patients, but not with patients who had “supranormal” thickness profiles. This indicates that altered thickness profiles may occur in those who experienced overt symptoms. Familial risk per se, while being insufficient to alter cortical thickness, may influence cortical gyrification patterns. Specific illness-related genetic variations (such as glutamate polymorphism) are more likely in patients/siblings with distinct morphological profiles (i.e., hypogyria).
The effort to discover morphologically homogeneous subgroups among schizophrenia has been expanding in recent years [Reference Weinberg, Lenroot, Jacomb, Allen, Bruggemann and Wells8, Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13, Reference Cobia, Csernansky and Wang15, Reference Xie, Zhang, Tang, Zhang, Yu and Gong79]. Despite marked methodological differences (scanner type, metric used, inclusion of healthy variance, etc.), many groups are now reporting at least two broadly consistent subgroups. One with impoverished cortex (reduced thickness or generalized cortical gray matter volume reduction), and the other with preserved or increased cortical grey matter tissue [Reference Weinberg, Lenroot, Jacomb, Allen, Bruggemann and Wells8, Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13, Reference Cobia, Csernansky and Wang15, Reference Pan, Pu, Chen, Huang, Cai and Tao20]. In the current study, we identified a similar solution, with a preserved thickness (cluster 2) group that displayed better cognitive function compared to others, and a subgroup with widespread cortical thinning (cluster 3) in line with our prior work [Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13, Reference Pan, Pu, Chen, Huang, Cai and Tao20] and previous studies [Reference Cobia, Csernansky and Wang15, Reference Xie, Zhang, Tang, Zhang, Yu and Gong79]. Again, in line with other prior studies [Reference Weinberg, Lenroot, Jacomb, Allen, Bruggemann and Wells8, Reference Liang, Heinrichs, Liddle, Jeon, Theberge and Palaniyappan13, Reference Cobia, Csernansky and Wang15, Reference Pan, Pu, Chen, Huang, Cai and Tao20, Reference Xie, Zhang, Tang, Zhang, Yu and Gong79], we found that the impoverished cortex subgroup had a more severe illness burden that we were able to observe using functioning data from follow-up assessments.
To our knowledge, prior studies have not utilized the variations in gyrification patterns alongside thickness for subtyping schizophrenia. This is somewhat surprising given the differences in genetic origins of these two measures [Reference Panizzon, Fennema-Notestine, Eyler, Jernigan, Prom-Wormley and Neale80], the differing spatial distribution of thickness and gyrification deficits in schizophrenia [Reference Rosa, Zugman, Cerqueira, Serpa, de Souza Duran and Zanetti81–Reference Palaniyappan and Liddle83], and the lack of concordance in their progressive trajectories among patients [Reference Pan, Chen, Chen, Lee, Liu and Ouyang84]. Consistent with these distinctions, we noted a specific, globally hypogyric subgroup of patients who had no concurrent thickness abnormalities in comparison with healthy subjects. But these patients had higher negative symptom burden, reduced verbal fluency, and persistent functional deficits, during the early phase of treatment (i.e., 2 months of antipsychotic treatment). These observations are in line with prior studies linking gyrification defects with negative symptoms and cognition [Reference Sallet, Elkis, Alves, Oliveira, Sassi and Campi de Castro85–Reference Takayanagi, Sasabayashi, Takahashi, Komori, Furuichi and Kido87]. These results indicated that the intrinsic mechanisms of cortical thickness decline and gyrification decline in schizophrenia are likely to be different, and also indicated the existence of different subgroups in three-dimensional topological morphology. In addition, we also observed a negative correlation (R = –0.35, P FDR = 0.011) between the cortical thickness of the right medial orbitofrontal and cumulative chlorpromazine equivalent of antipsychotics, indicating that higher dose exposure has a limited effect on the overall cortical morphology.
One of the main goals for subtyping schizophrenia subgroups is to pursue the mechanistic origins of this illness. Our findings relating glutamatergic postsynaptic receptor genes (GRIN2A, GRM3) and CACNA1C gene to subgroups provide an interesting lead in this regard. The subgroup of patients with the most preserved cortical morphology showed no alteration in the glutamate receptor or VGCC risk allele distribution. But the hypogyric subgroup exhibited a higher than expected presence of the schizophrenia-risk SNVs of GRIN2A (rs9940680 and rs1420040) and lower rate of the protective variant of GRM3 (s14513928; a lower rate of this protective variant is also notable in healthy siblings). GRIN2A gene is involved in the synthesis of NMDA receptor complex components, which could directly regulate the permeability of receptors [Reference Swanger, Chen, Wells, Burger, Tankovic and Bhattacharya88]. Increased receptor permeability can lead to an increase in the concentration of cations and glutamate in nerve cells, resulting in excitotoxicity [Reference Plitman, Patel, Chung, Pipitone, Chavez and Reyes-Madrigal47, Reference Shah, Plitman, Iwata, Kim, Nakajima and Chan48]. Shifts in Ca2+ currents through NMDA receptors can lead to notable neuronal apoptosis [Reference Sachser, Santana, Crestani, Lunardi, Pedraza and Quillfeldt89], which if occurring during early development, can lead to altered gyrification patterns. GRM3 gene was involved in the synthesis of metabolic glutamate receptor subunit 3, which can inhibit cyclic adenosine monophosphate [Reference Schaffhauser, Cai, Hubalek, Macek, Pohl and Murphy90] and reduce Ca2+ flow through NMDA receptor, and achieve antioxidant and antineurotoxic effects. Interestingly, animal studies indicate that aberrant gyrification may arise from deficits in astroglial support [Reference Shinmyo, Saito, Hamabe-Horiike, Kameya, Ando and Kawasaki91], which may also relate to glutamate-mediated toxicity during development [Reference Rosenberg and Aizenman92]. Note that the effect of these genetic variations may differ across brain regions (see Supplementary Figure S1).
The subgroup with the most prominent thickness reduction (“impoverished cortex”) had a significantly different mutation rate in CACNA1C (higher in risk SNV rs223906, lower in protective SNV rs556844413). In previous studies, the CACNA1C was generally known as one of the regulators of Ca2+ signaling in the proliferation and survival of neural progenitors [Reference Smedler, Abe, Palsson, Ingvar and Landen45, Reference Malmersjo, Rebellato, Smedler, Planert, Kanatani and Liste93]. Previous studies found Ca2+ signaling plays an important role in glutamate-mediated excitotoxicity [Reference Afanador, Mexhitaj, Diaz, Ordonez, Baker and Angulo94, Reference Ye, Shi and Yin95], which could deduce neuronal apoptosis [Reference Sachser, Santana, Crestani, Lunardi, Pedraza and Quillfeldt89, Reference Croce, Bernardini, Di Cecca, Caltagirone and Angelucci96] resulting in the reduction of cortical thickness [Reference Soeiro-de-Souza, Lafer, Moreno, Nery, Chile and Chaim44, Reference Smedler, Abe, Palsson, Ingvar and Landen45, Reference Wang, Chen, Cao, Dou, Fu and Zhang97]. Nevertheless, the lack of association between the impoverished thickness subgroup and glutamatergic receptor AMPA/NMDA variations suggests a less direct role for NMDA hypofunction, and a more direct role for aberrant glutamate dynamics [Reference Ermolyuk, Alder, Surges, Pavlov, Timofeeva and Kullmann98] in this subgroup. Besides, the characteristics of treatment-resistant patients indeed overlaps with our “impoverished cortex” subgroup, such as more extensive disruption of brain structure [Reference Shah, Plitman, Iwata, Kim, Nakajima and Chan48, Reference Fan, Huang, Tan, Wang, Li and Chen99, Reference Kitajima, Tamura, Sasabayashi, Nakajima, Iwata and Ueno100], worse cognitive performance [Reference Fan, Huang, Tan, Wang, Li and Chen99, Reference Syeda, Wannan, Merritt, Raghava, Jayaram and Velakoulis101] and poor functional recovery. Previous studies found that patients with Treatment-Resistant Schizophrenia have increased glutamate levels in the anterior cingulate cortex, whereas dopamine synthesis in the striatum does not differ from controls in this subgroup of patients [Reference Egerton, Murphy, Donocik, Anton, Barker and Collier102, Reference Mouchlianitis, Bloomfield, Law, Beck, Selvaraj and Rasquinha103]. One predominant model of dopaminergic dysfunction is that it may occur downstream of excitation/inhibition abnormality in the glutamate system [Reference McCutcheon, Krystal and Howes104]. While we do not have striatal positron-emission tomography (PET) data to confirm this, our observations relating to structure and glutamate receptor polymorphisms are broadly consistent with prior evidence implicating glutamatergic/dopaminergic mechanisms in poor outcomes of schizophrenia.
In summary, the regulation of glutamatergic postsynaptic receptor genes and permeability of calcium channels may regulate cell apoptosis and death, affecting brain morphology both during development and in later life. CACNAC1, GRIN2A, and GRM3 have already been targeted in clinical repurposing trials for schizophrenia [Reference Lago and Bahn105, Reference Stansley and Conn106], but these clinical trials also showed heterogeneous outcomes [Reference Lencz and Malhotra107]. The presence of a substantial number of patients (~32.5% in this sample) with supra or near-normal cortical morphology and glutamatergic receptor genes argues for the existence of a non-glutamatergic type of schizophrenia with preserved gyrification, explaining, in part, the heterogeneity of clinical trial outcomes.
The current study has several strengths (multimetric clustering, inclusion of siblings, cluster-based genetic association analysis, and follow-up sample to track functioning); nevertheless, several limitations should also be considered. First, despite a deep sequencing of the target region, we only studied selected genetic locus that we could link a priori with schizophrenia susceptibility and glutamate hypothesis; thus, we did not include all glutamatergic postsynaptic receptor genes. We cannot exclude the regulatory effect of other genes on targeted genes. Second, we lacked transcriptome data to further understand the intermediate process of glutamatergic postsynaptic receptor genes affecting brain morphology. Third, we had limited longitudinal data to study long-term prognostic associations.
To conclude, by linking MRI-derived cortical morphological patterns to glutamatergic and calcium channel variations, we highlight the potential to select patients with certain neuroanatomical features when studying interventions that regulate glutamate/calcium channels in schizophrenia.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1192/j.eurpsy.2023.2408.
Acknowledgments
L.P. acknowledges research support from Monique H. Bourgeois Chair in Developmental Disorders and Graham Boeckh Foundation (Douglas Research Centre, McGill University) and a salary award from the Fonds de recherche du Quebec-Sante ́ (FRQS). This work was also supported by the National Natural Science Foundation of China (Y.P., grant number 82201664; G.W., grant number 81701325), (X.O., grant number 81801353), (Z.L., grant number 82071506), and the Scientific Research Launch Project for new employees of the Second Xiangya Hospital of Central South University (Yunzhi Pan), and Changsha Natural Science Foundation (G.W.).
Author Contributions
Yunzhi Pan conceptualized the article, analysed the data, and critically revised the manuscript; LP contributed to the conception, analytical plan, interpretation of results and writing of the manuscript including its critical revisions. Xuan Ouyang completed the first draft of the article and participated in sample recruitment, scale evaluation, and scanning of magnetic resonance imaging. Xudong Chen, Guowei Wu, Yixin Cheng, Wenjian Tan, Manqi Zhang, Mengjie Deng participated in sample recruitment, clinical evaluation, and scanning of magnetic resonance imaging. Zhening Liu helped conceptualize the article and revise the manuscript.
Competing interest
L.P. reports personal fees from Janssen Canada, Otsuka Canada, SPMM Course Limited, UK, Canadian Psychiatric Association; book royalties from Oxford University Press; investigator-initiated educational grants from Janssen Canada, Sunovion and Otsuka Canada outside the submitted work. No other competing interests were declared. All other authors have declared that there are no conflicts of interest in relation to the subject of this study.
Open access
Comments
No Comments have been published for this article.