Early treatment response is proposed to be one of the robust predictors of schizophrenia patients’ functional outcomes (Lambert et al., 2008), but it is variable and cannot be precisely predicted by the treating physician. Hence, it is critical to ascertain biomarkers with a predictive potential by investigating their relation to prospective treatment response in the early phase of schizophrenia.

Subcortical structures, including the basal ganglia and some areas of the limbic system, have attracted interest in schizophrenia research (Okada et al., 2016; Van Erp et al., 2016). They contain abundant dopaminergic neurons, the disruption of which is thought to be related to schizophrenia pathology and the formation of psychotic symptoms such as hallucinations and delusions (Howes and Kapur, 2009). Moreover, the remarkable 108 candidate gene loci of schizophrenia identified previously (Schizophrenia Working Group of the Psychiatric Genomics, 2014) also include many genes relevant to dopamine pathways. These studies link both subcortical structure phenotypes and schizophrenia pathology to the dopamine signaling pathways. Moreover, a recent large-scale imaging-genetic study detected the genetic overlap between schizophrenia risk and volumes of subcortical structures including the hippocampus and putamen using genome-wide association study (GWAS) tools (Smeland et al., 2018). The evidence above implies common genetic associations emerging for both subcortical structure phenotypes and schizophrenia pathology. In addition, subcortical structures are the main targets of antipsychotic drugs, and their neuroimaging measurement alterations have also been suggested to be in relation to a lack of response to neuroleptic treatment (Sarpal et al., 2016; Li et al., 2020).

In the domain of neuroimaging, functional connectivity measures of the subcortical areas have been suggested to be an excellent marker for predicting antipsychotic treatment outcomes (Sarpal et al., 2016; Li et al., 2020). The volume of the subcortical structures also shows predictive potential for neuroleptic treatment response (Buchsbaum et al., 2003; Strungas et al., 2003). However, less is known about the predictive role of their structural connectivity patterns. The principal imaging approach available for this purpose is the structural covariance analysis of T1-weighted magnetic resonance imaging (MRI). The “structural covariance,” i.e., morphological correlation, is an organizational pattern in the brain measured across a population and assesses the statistical associations of pairs of brain regions in their anatomical properties such as volume (Alexander-Bloch et al., 2013a,b). Researchers have demonstrated the abnormalities of structural covariance in Parkinson’s disease (Chou et al., 2015; De Schipper et al., 2017; Xu et al., 2017), a neurodegenerative disease caused by dopamine deficiency in the basal ganglia, suggesting the potential associations between structural covariance abnormalities and dopamine signals. Moreover, a positron emission tomography (PET) study revealed the relationship between the baseline dopamine synthesis capacity of subcortical regions and improvements in positive and negative symptoms after treatment (Jauhar et al., 2019). These previous studies led us to speculate that subcortical structural covariance before starting treatment may underlay future variation in response to antipsychotic treatment. Yet, few studies have directly investigated how baseline subcortical morphological covariance relates to treatment outcome in the early phase of schizophrenia.

Another unsolved issue is that less is known about whether the variance of subcortical morphological connectivity is driven by gene expression patterns and epigenetic mechanisms, which are the plausible molecular basis of phenotypic heterogeneity across individuals with schizophrenia. Epigenetic modifications are proposed to mediate between environmental insults and gene expression and alter and stably maintain the expression of genes (Magwai et al., 2021). DNA methylation (DNAm) is the most widely explored epigenetic mechanism and has been suggested to be associated with an antipsychotic drug action mechanism (Ovenden et al., 2018). Although DNAm is to a degree tissue-specific, 10.9% of DNAm sites were moderately robustly associated (r > 0.5) between brain and blood (Braun et al., 2019), implying that DNAm status of peripheral blood cells may act as a surrogate for that of central tissues. Noteworthily, structural covariance is thought to be specifically related to co-expression of a set of genes relevant to neurobiological processes (Romero-Garcia et al., 2018), which suggests the heritability of inter-regional structural covariation. Recent brain expression atlases bridge the gap between epigenetic modifications and brain connectome phenotypes. The Allen Human Brain Atlas (AHBA), a publicly available transcriptomic dataset (Hawrylycz et al., 2012), has been utilized to identify transcriptomic signatures associated with brain network connectivity of individuals with mental disorders, such as schizophrenia and major depressive disorder (MDD) (Morgan et al., 2019; Li et al., 2021), which uncovers the molecular foundation of regional brain vulnerability to major mental disorders. Combining brain network phenotypes, brain gene expression, and DNAm is imperative in the developing of baseline multi-omic biomarkers associated with treatment outcome, as it can integrate multiple omics data to comprehensively understand efficacy heterogeneity across individuals, although few such studies have been conducted.

This current study, therefore, integrates multi-omic measures including subcortical covariance network, AHBA transcriptome data, and peripheral DNAm to comprehensively understand efficacy heterogeneity across individuals. Patients were treated with 8-week risperidone monotherapy to control for the effects of confounders related to multidrug therapy on treatment response. They were divided into responders and non-responders according to the Remission in Schizophrenia Working Group (RSWG) (Andreasen et al., 2005). We applied a seed-based multivariate technique to identify the patterns of subcortical structural covariance in responders and non-responders. Based on previous evidence (Buchsbaum et al., 2003; Strungas et al., 2003; Xu et al., 2017; Ovenden et al., 2018; Romero-Garcia et al., 2018; Yee et al., 2018), it was hypothesized that: (1) patients’ baseline subcortical structural covariance would be associated with patients’ treatment response; (2) variances of baseline subcortical structural covariance between responders and non-responders would be spatially correlated with brain gene expression acquired from the AHBA, and these genes would be enriched in dopaminergic pathways and neurobiological processes; and (3) the DNAm levels of gene of interest (GOI) would also relate to patients’ treatment response. The GOI was defined as the overlapping genes of AHBA genes spatially correlated with baseline brain measures and the 108 schizophrenia candidate loci (Schizophrenia Working Group of the Psychiatric Genomics, 2014), as they were associated with both the schizophrenia pathology and variations in subcortical structural covariance related to treatment response.

This study is among the first to investigate baseline biomarkers related to treatment outcomes in the early phase of schizophrenia using multi-omic measures linking subcortical covariant network, transcriptomic signatures, and peripheral epigenetic modifications. As hypothesized (Figure 4), the subcortical structural covariance before starting treatment may serve as potential predictors of treatment response in psychotic, disorganized, and total symptoms. Phenotypic variations of subcortical connectome related to psychotic symptom response were spatially associated with the expression of genes primarily enriched in neurobiological processes and the dopaminergic pathway. Moreover, the DNAm of GOI demonstrated significant associations with patients’ improvement of psychotic symptoms.

In the psychotic symptom dimension, non-responders had a higher baseline structural covariance in the putamen–hippocampus–pallidum–accumbens pathway compared with responders, implying that a higher level of structural covariance among these subcortical regions may predict a poor treatment response in this symptom dimension. The above subcortical structures, i.e., the hippocampus, pallidum, accumbens, and putamen, have been proposed to be critical for the reward system (Smith et al., 2009; Mizuno et al., 2016; Legates et al., 2018), which are clearly implicated in the pathophysiology of schizophrenia and pharmacological mechanisms of antipsychotic treatment response (Schultz et al., 1997; Kapur and Mamo, 2003). Our approach demonstrates for the first time that a higher baseline structural covariance in the “subcortical reward circuit” may be associated with a poor response to the psychotic symptoms. It is known that antipsychotic agents bind to dopamine receptors primarily in subcortical structures (Eisenberg et al., 2017). Interestingly, a PET study demonstrated that dopamine synthesis capacity of the striatum explained > 40% of the variance in subsequent changes of psychotic symptoms after treatment, suggesting that dopaminergic dysfunction before starting treatment may underlie future variations in treatment response (Jauhar et al., 2019). We speculate that the predictive effect of the putamen—hippocampus–pallidum–accumbens circuit in this study might be mediated by patients’ altered baseline dopamine levels in these subcortical regions. Further studies in combination with PET and MRI are needed to clarify this speculation.

We also found that subcortical circuit biomarkers are associated with treatment efficacy in disorganized symptoms. Specifically, we found consistently decreased baseline structural covariance in the putamen–hippocampus–pallidum–thalamus circuits in non-responders compared with responders. Many recent studies have reported the involvement of the putamen–hippocampus–pallidum–thalamus loop in the working memory (Kalivas et al., 2001; Karlsgodt et al., 2005). Notably, the disorganized symptoms include linguistic symptoms, such as poverty of speech, poverty of content of speech, tangentiality, and derailment (Andreasen et al., 2005). Evidence shows the associations of disorganized symptoms with verbal working memory (Horan et al., 2008; Torniainen et al., 2012). Our preliminary result implies that consistent reductions of baseline structural covariance in the subcortical verbal working memory circuit may predict a poor response to disorganized symptoms.

Subcortical structural covariance is built on the similarity of macrostructural variations. It reflects structural synchronized maturation and similarities in the local micro-architectonic properties (Alexander-Bloch et al., 2013a,b), which is commonly influenced by many certain factors, such as monosynaptic connection in synapse development (Yee et al., 2018) and gene expression (Raznahan et al., 2011). In this study, we detected significant spatial associations between the variation in subcortical volume covariance related to psychotic symptom response and expression of genes enriched for 19 GO biological processes and 1 KEGG pathway. It is known that treatment response to psychotic symptoms has a significant association with baseline dopamine synthesis capacity of subcortical regions (Jauhar et al., 2019). Interestingly, the identified KEGG pathway, i.e., cGMP–PKG signaling pathway, has been reported to mediate the phosphorylation of DARPP-32 in neostriatal neurons induced by glutamate (Nishi et al., 2005). Significantly, mice lacking DARPP-32, a transduction molecule in dopamine signaling that is selectively enriched in the striatum, exhibit profound deficits in their behavioral responses to antipsychotic medication (Fienberg et al., 1998), which suggests the potential role of the cGMP–PKG pathway in dopaminergic neurotransmission and neuroleptic treatment response. Moreover, the detected GO biological processes were primarily involved in neurobiological processes. Interestingly, the schizophrenia risk genes identified previously (Schizophrenia Working Group of the Psychiatric Genomics, 2014) were significantly enriched in the 991 PLS1 genes detected in this study, implying a common genetic basis emerging for both schizophrenia pathology and response-related subcortical covariant network phenotypes. Consistently, a recent large-scale GWAS-based imaging-genetic study detected evidence of genetic overlap between subcortical brain measures and schizophrenia risk (Smeland et al., 2018).

However, we did not detect spatial correlations in the total or disorganized symptom dimension, suggesting that the phenotypic variations of subcortical connectome related to treatment response of the two symptom dimensions may not be transcriptomically underlaid. Future studies are still needed to validate this speculation.

Importantly, we also detected significant correlations between patients’ baseline DNAm of GOI, i.e., HCN1 and MED19, and the longitudinal psychotic symptom alterations. These associations support a functional role for the peripheral DNAm alterations and suggest the validity of patient classification based on subcortical volumetric covariance network markers. These findings, combined with the above transcriptomic basis of variations in subcortical connectome, imply that treatment response to psychotic symptoms may be different from that of other dimensions, as it may have more robust molecular bases. Noteworthily, HCN1, coding for hyperpolarization-activated cyclic nucleotide-gated channels, has been implicated in modulation of dopamine signaling (Rampino et al., 2019). MED19 is proposed to be involved in cancer growth, and its expression inhibits the spread and growth of cancers (Li et al., 2011; Zhu et al., 2013). This study first reveals the associations between baseline DNAm of these two genes and psychotic symptom improvement. Combining brain connectomes, gene transcripts, and DNAm could provide insight into how macroscale brain connectivity impairments are driven by the microscale architecture.

Our results must be interpreted in light of several limitations. First, the sample size of responders and non-responders in each symptom dimension is relatively small due to the challenging requirement of first-episode drug-naïve patients given risperidone monotherapy. Replication in other datasets with a large sample size is needed. Second, we only included tissue samples in the left hemisphere, as only two of the six donors in the AHBA have brain tissue samples in the right hemisphere. Thus, the association between gene expression and treatment response-related alterations in subcortical structural covariance does not represent the condition of bilateral subcortical structures. Third, the DNAm status was derived from peripheral cells rather than brain tissue samples of the subcortical structures, although DNAm sites were moderately and robustly associated between brain and blood (Braun et al., 2019).

This study links structural and molecular alterations relevant to therapeutic response in the early phase of schizophrenia. Subcortical structural covariance and peripheral DNAm could be useful predictors of antipsychotic treatment response, and these results are important for future precision medicine. This study also defines a roadmap for future studies investigating multimodal imaging epigenetic biomarkers for treatment response in schizophrenia.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

The studies involving human participants were reviewed and approved by the Medical Ethics Committee of Henan Mental Hospital and Second Xiangya Hospital of Central South University. The patients/participants provided their written informed consent to participate in this study.

XZ and ML: study design. XZ, TY, MH, and ZL: data acquisition. CH, XH, JX, LL, ML, XD, and JZ: statistical analysis. XZ, CH, and ML: manuscript writing. All authors have read and agreed to the published version of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.