Longitudinal multi-omics alterations response to 8-week risperidone monotherapy: Evidence linking cortical thickness, transcriptomics and epigenetics

Background: Antipsychotic treatment-related alterations of cortical thickness (CT) and clinical symptoms have been previously corroborated, but less is known about whether the changes are driven by gene expression and epigenetic modifications.

Methods: Utilizing a prospective design, we recruited 42 treatment-naive first-episode schizophrenia patients (FESP) and 38 healthy controls. Patients were scanned by TI weighted imaging before and after 8-week risperidone monotherapy. CT estimation was automatically performed with the FreeSurfer software package. Participants' peripheral blood genomic DNA methylation (DNAm) status, quantified by using Infinium^® Human Methylation 450K BeadChip, was examined in parallel with T1 scanning. In total, CT measures from 118 subjects and genomic DNAm status from 114 subjects were finally collected. Partial least squares (PLS) regression was used to detect the spatial associations between longitudinal CT variations after treatment and cortical transcriptomic data acquired from the Allen Human Brain Atlas. Canonical correlation analysis (CCA) was then performed to identify multivariate associations between DNAm of PLS1 genes and patients' clinical improvement.

Results: We detected the significant PLS1 component (2,098 genes) related to longitudinal alterations of CT, and the PLS1 genes were significantly enriched in neurobiological processes, and dopaminergic- and cancer-related pathways. Combining Laplacian score and CCA analysis, we further linked DNAm of 33 representative genes from the 2,098 PLS1 genes with patients' reduction rate of clinical symptoms.

Conclusions: This study firstly revealed that changes of CT and clinical behaviors after treatment may be transcriptionally and epigenetically underlied. We define a “three-step” roadmap which represents a vital step toward the exploration of treatment- and treatment response-related biomarkers on the basis of multiple omics rather than a single omics type as a strategy for advancing precise care.

Introduction

Cortical thickness (CT) relates to the number of neurons and the neuropil within a cohort of cortical neurons originating from a single neuronal progenitor (1). Measurement of CT may therefore allow probing changes in brain maturation and growth, which are increasingly being proposed as neuropathological mechanisms of schizophrenia. CT mapping in this disease has indeed shown abnormalities of the cortical areas, particularly in parts of the temporal and frontal lobes (2). Antipsychotics are commonly used as first-line regimens for the treatment of schizophrenia, as well as bipolar disorders. Noteworthily, accumulating evidence supports that antipsychotic agents may affect CT in schizophrenia patients, inducing excessive thinning widely distributed in most cortical areas (3–7), which has raised concerns about the potential neurotoxic side effects of antipsychotic agents that might need to be taken into consideration when prescribing them.

Although antipsychotic therapy may have potential contribution to cortical morphological alterations, less is known about whether treatment-related changes of CT after treatment are driven by genetic factors. Recent brain expression atlases, bridge the gap between genetic factors and brain morphology phenotypes. The Allen Human Brain Atlas (AHBA), a publicly available transcriptomic dataset (8), has been utilized to identify transcriptomic signatures associated with brain structural alterations of individuals with mental disorders, such as schizophrenia and major depressive disorder (MDD) (9, 10), which uncovers the molecular foundation of regional brain vulnerability to major mental disorders. By using the transcriptomic dataset of AHBA, researchers have uncovered the molecular foundation of regional brain vulnerability to major mental disorders (10). However, less is known about whether CT may be specifically related to co-expression of a set of genes relevant to neurobiological processes.

Epigenetic modifications, some heritable changes that are not due to alterations in DNA sequence, mediate between gene expression and environmental insults, and then alter and stably maintain gene expression levels (11). DNA methylation (DNAm), among the epigenetic modifications, is a stable epigenetic mechanism. Recent research on psychiatric disorders has only begun to investigate the associations between DNAm alterations and antipsychotic therapy outcomes (12). Combining of brain anatomical phenotypes, brain gene expression and DNAm is important for developing treatment- and treatment outcome-related biomarkers, as it can integrate multi-omics data to comprehensively explain heterogeneity of treatment and treatment response across individuals, although few researches have been conducted recently. Noteworthily, although DNAm status is reported to be tissue-specific, 10.9% CpG sites were strongly (r > 0.5) associated between brain tissues and peripheral blood cells (13), implying that DNAm of peripheral blood cells may play as a surrogate for DNAm levels of the brain.

This study firstly integrates multi-omics data including peripheral DNAm, AHBA transcriptomic data and CT to comprehensively explain the transcriptional and epigenetically basis underlying CT alterations and clinical symptom improvement after acute antipsychotic treatment. Patients were given risperidone monotherapy for 8 weeks to control for the confounding effects of multidrug therapy on treatment response. CT estimation was automatically performed with the FreeSurfer software package. Partial least squares (PLS) regression was used to detect the spatial associations between longitudinal CT variations after treatment and cortical transcriptomic data acquired from the AHBA. Patients' peripheral blood genomic DNAm, examined in parallel with T1 scanning at the two time points, was quantified by using Infinium^® Human Methylation 450K BeadChip. We then performed Canonical Correlation Analysis (CCA) to identify multivariate associations between patients' clinical symptom improvement and their longitudinal alterations of DNAm of PLS1 genes. On the basis of previous evidence, it was hypothesized that: (1) alterations of CT and clinical symptoms after antipsychotic treatment would be transcriptionally and epigenetically underlied; (2) phenotypic variations of CT after treatment would be spatially correlated with brain gene expression acquired from the AHBA, and these genes would be enriched in dopaminergic-related pathways and neurobiological processes; (3) the DNAm levels of part of PLS1 genes would relate to patients' clinical response.

Methods

Participants

We recruited 42 drug-naive patients with first-episode schizophrenia (Patients-0W) and their matched healthy controls (n = 38) in Henan Mental Hospital (China) from 12/2012 to 01/2014. This dataset has been utilized by our group (14–16). All patients were diagnosed by experienced psychiatrists based on the Structured Clinical Interview for DSM-IV-TR. Patients' disease duration was <12 months. SCID-non-patient edition was used to scan controls to ensure that they have no history of other mental disorders or neurological diseases.

The study was performed following the Declaration of Helsinki. It was approved by Ethics Committees of the Henan Mental Hospital, Xinxiang, China and the Second Xiangya Hospital of Central South University, Changsha China (protocol code S088, 2012). All subjects in this study provided written informed consent.

Antipsychotic therapy

The starting dose of risperidone was 1 or 2 mg/day, and then slowly titrated after 1 week. Specifically, we increased its dosage at 1-week interval until the patients had obvious clinical improvement. For those individuals who did not show satisfactory improvement after 4 weeks, they were permitted to reach a maximum daily dose of 6 mg/day. All 42 patients were prescribed with risperidone for 8 weeks, and not allowed to take mood stabilizers or antidepressants. The therapeutic safety was assessed weekly through clinical interviews.

Psychopathological assessments

Patients' symptomatic severity was assessed both at baseline and after 8-week treatment by using the Positive and Negative Syndrome Scale (PANSS) (17), which consists of general psychopathology (PANSS-G, 16 items, G1–G16), negative (PANSS-N, 7 items, N1–N7), and positive (PANSS-P, 7 items, P1–P7). The improvement in each symptom dimension was measured by the reduction rate (18). Reduction rates for positive, negative and general psychopathology symptoms were respectively equal to PANSS-P_0W—PANSS-P_8W/PANSS-P_0W-7, PANSS-N_0W—PANSS-N_8W/PANSS-N_0W-7, as well as PANSS-G_0W—PANSS-G_8W/PANSS-G_0W-16.

“Three-steps” analysis detecting longitudinal multi-omics alterations response to treatment

We explored biomarkers correlated to treatment, treatment response as well as schizophrenia pathology using multi-omics measures linking CT measures, transcriptomic signatures and peripheral DNAm values based on the following steps (Figure 1).

FIGURE 1

Figure 1. Study overview of imaging-transcriptomics-epigenetics analysis.

Step 1: CT biomarkers related to schizophrenia and treatment

T1 imaging acquisition

The participants underwent T1 weighted images on a 3.0T Siemens MRI scanner (Vero) at the Magnetic Imaging Center, the Henan Mental Hospital, Xinxiang, China. Patients were scanned both before and after 8-week treatment, while healthy controls only underwent the baseline neuroimaging assessment. Detailed MRI parameters were shown in the Supplementary material. There were 4 patients withdrawing the follow-up T1 scans. We therefore obtained imaging data of 42 Patients_0W, 38 patients after treatment (Patients_8W), and 38 healthy volunteers.

Image analysis

CT estimation was automatically processed with the FreeSurfer software package (v5.3.0; http://surfer.nmr.mgh.harvard.edu/) according to the programme segmentation procedure, the technical details of which were described in the prior publication (19).

Quality control of T1W data

To check the differences of image quality and head motion among the three groups, we computed the Euler number (20) for each MRI image. This quality control method was proposed as a way to quantitatively evaluate image quality (20–22). We then respectively used independent and paired samples t-test to compare the case-control and longitudinal differences of Euler number. We detected no significant case-control or longitudinal differences in the Euler number (Ps > 0.05).

Extraction of CT

The cerebral cortex were subdivided into 68 areas (34 per hemisphere) following the Desikan–Killiany–Tourville (DKT) cortical labeling protocol (23). Each of the 68 regions was defined by the automated FreeSurfer segmentation procedure. We also extracted Total intracranial volume (TIV) to adjust for the confounding factors of brain sizes.

Construction of t-statistic maps

All the CT measures of 68 regions were residualized with respect to TIV, sex and age utilizing generalized linear models before performing between-group comparisons. We respectively used independent samples t-test and paired samples t-test to compare the case-control and longitudinal differences in CT measures, and then obtained case-control and longitudinal t-statistic maps.

Step 2: Transcriptomic basis of variations in CT related to schizophrenia and treatment

Preprocessing of AHBA data

Transcriptional profiles, consisting of expression measures of 20,737 genes (from 58,692 probes), were acquired from the public database AHBA (http://www.brain-map.org) (8). The detailed information about AHBA was shown in Supplementary material. We preprocessed the transcriptional data according to previously published 5 steps (24): (1) annotations from probes to genes; (2) filtering probes; (3) screening representative probes; (4) assigning DKT-68 atlas; (5) normalizing expression levels. The preprocessing related codes are available at github (details see https://github.com/BMHLab/AHBAprocessing respectively). After preprocessing, we obtained 10,027 probes in each brain sample. Not all donors have samples in the right hemisphere (24), we therefore did not analyze tissue samples in the right hemisphere.

PLS regression analysis

We respectively coregistered the two t-statistic maps (left hemisphere) with transcriptional maps (brain tissues of the left hemisphere in the AHBA) according to the DKT-68 atlas. Partial least squares (PLS) regression (9, 21) was then conducted to evaluate the spatial correlations between transcriptional scores and neuroimaging t-statistic maps. In the PLS analysis, the t-statistic map (for 34 brain regions) was set as response variables, and expression data of the 10,027 probes for the 34 brain areas were predictor variables. The first component (PLS1) obtained in the PLS analysis was the linear combination of transcriptional profiles significantly related to the t-statistic map. We conducted Permutation test (1,000 times) to test the null hypothesis that the PLS1 explained no more covariance between t-statistic maps and gene expression levels than expected by chance. We then performed Bootstrapping (500 bootstrap samples) to measure the weight of each PLS1 gene. The ratio of each region's weight to its standard error in the bootstrap was assessed as the Z-values. The normalized weights of each PLS1 gene were ranked through one-sample Z-tests. After FDR correction, we obtained gene lists (with P_FDR < 0.05) contributing to the PLS1.

Enrichment analysis

We performed enrichment analysis for the above detected PLS1 genes with |Z| > 3 (21) (all FDR < 0.05), using an online tool Metascape (https://metascape.org), which provides automated meta-analysis to obtain enriched biological processes (gene ontology, GO) and pathways (Kyoto encyclopedia of genes and genomes, KEGG) (25). We filtered enrichment terms to those pathways and biological processes that met P_FDR < 0.01.

Step 3: Epigenetic biomarkers

DNAm of PLS1 genes

The peripheral blood DNAm status in both Patients_8W (n = 38) and Patients_0W (n = 38) was evaluated in parallel with T1W scanning. Genomic DNAm values were measured in the above samples by using Illumina 450K Genechip. The epigenetic dataset has been used previously by our group (15, 16). The average value of all the CpG sites in each PLS1 gene was used to represent the DNAm level of this gene. We demonstrated detailed information about DNA extraction, quality control (QC), bisulfite conversion, Genechip analysis, microarray data processing in Supplementary material.

CCA analysis

Canonical correlation analysis (CCA) (26), a data-driven analysis among datasets each consisting of multiple variables, was utilized to identify multivariate associations between DNAm of PLS1 genes and clinical data. CCA aims to detect optimal linear combinations of multivariate epigenetic and behavioral measures by maximally relating the DNAm levels to clinical variables. The CCA analysis was performed as follows:

$\begin{array}{l}{U}_i=X{a}_i\end{array}$

$\begin{array}{ll}{V}_i=Y{b}_i,& (1)\end{array}$

where X is the selected gene features by Laplacian score importance, and Y is the clinical behaviors (PANSS scores). Laplacian score (27), a widely used feature selection method, was performed to identify the most representative features from the high-dimensional datasets, i.e., DNAm values of PLS1 genes. Laplacian score was performed as follows:

$\begin{array}{ll}{L}_r=\frac{{\tilde{x}}_r^{T}L{\tilde{x}}_r}{{\tilde{x}}_r^T{D}_g{\tilde{x}}_r},& (2)\end{array}$

where L is Laplacian Matrix, D_g is the degree matrix of L, and ${\tilde{x}}_r^{}$ is the center feature of the r^th original DNAm feature. The DNAm data has N features, and the Laplacian score for the methylation level of the r^th gene is L_r.

The CCA mode is aimed at using a_i and b_i to maximize the Pearson correlation between U_i and V_i, and then the most important DNAm features associated with the clinical behavioral scores can be selected according to the significance of CCA mode.

Here, we utilized CCA analysis to compute the multivariate correlations between patients' baseline DNAm of PLS1 genes and baseline clinical behavioral scores (PANSS_P, _N, _G), as well as between longitudinal alterations of patients' DNAm of PLS1 genes (baseline DNAm_0W–DNAm_8W) and patients' reduction rate of PANSS_P, PANSS_N, and PANSS_G. During performing CCA, the DNAm values of PLS1 genes and symptomatic scores were standardized with z-score transformation.

Taking the longitudinal data as an example, the CCA analysis started with the construction of epigenetic and symptomatic matrices, X_n×s and Y_n×t (n = 38, s = number of representative features of longitudinal alterations of patients' DNAm of PLS1 genes in the Laplacian dimension reduction analysis, t = 3 i.e., PANSS_P, _N, and _G). For each CCA mode, we utilized the permutation test (10,000 times) to test the significance of canonical correlations (26). A CCA mode would be valid only if the permutation test (10,000 times) reached significance (P_FWE < 0.05). If at least one significant CCA mode was identified for the PANSS symptoms, the linear associations between columns of the CCA clinical symptom variate (or estimated clinical projection matrix) with each of the original DNAm variables were estimated using univariate variate-to-variable correlations, i.e., Pearson's correlations. All P-values in the Pearson's correlation analyses were explicitly corrected by false discovery rate (FDR). Code about the CCA analysis is available at the previously published work of other teams (28, 29).

We then further analyzed whether there are overlapping genes of the first significant CCA mode (CCA1) genes and the dopamine-related genes. We used an online tool PathDIP (http://ophid.utoronto.ca/pathDIP) (30) to search dopamine-related gene lists, with the input pathway terms “dopamine,” “dopamine-receptor-mediated,” “dopamine-receptors,” and “dopaminergic.”

Results

Demographics and clinical behaviors

We found no significant between-group differences in the demographic characteristics, see Supplementary Table S1. Patients had significant clinical improvement in positive and general psychopathology symptoms (Ps < 0.001; Supplementary Table S2), while their alterations in negative symptom dimension showed marginal significance (P = 0.060; Supplementary Table S2).

Between-group comparisons of CT

Compared with healthy volunteers, patients_0W had significant increases of CT in the R_pericalcarine (Supplementary Table S3, P = 0.006, FDR corrected), while marginally significant increases in the L_entorhinal, R_entorhinal, and L_pericalcarine cortices (Supplementary Table S3, P = 0.0892, 0.0763, and 0.0794 respectively, FDR corrected).

Significant cortical thinning (Supplementary Table S3, Ps < 0.05, FDR corrected) were observed in bilateral caudalmiddlefrontal, inferiortemporal, lateralorbitofrontal, parsopercularis, parsorbitalis, parstriangularis, rostralmiddlefrontal, superiorfrontal, superiorparietal (R_superiorparietal showed marginal significance), superiortemporal, supramarginal, as well as L_caudalanteriorcingulate, L_fusiform, L_inferiorparietal, L_medialorbitofrontal, L_middletemporal, L_precentral, R_bankssts, R_ posteriorcingulate, R_precuneus and R_insula in Patients_8W relative to Patients_0W.

Spatial associations between gene expression and schizophrenia-related CT alterations

We did not identify significant PLS components associated with the case-control t-statistic map of the CT measures (r = 0.181, P > 0.05).

Spatial associations between longitudinal CT alterations and gene expression

We constructed the t-statistic map representing the longitudinal variations of CT after treatment in the patient group (Figure 2A). As shown in Figure 2B, we detected significant PLS components associated with that t-statistic map (r = 0.558, P < 0.0001), and the PLS1 component had 2,098 genes with normalized PLS1 weights |Z| > 3, including 927 genes with Z > 3, and 1,171 genes with Z < −3 (all P < 0.001 with FDR correction). The PLS1 explained 55.8% of the variance in the longitudinal differences of CT after treatment. We then used Metascape to align the enriched pathways and biological processes for the PLS1 genes. The PLS1 genes were significantly enriched in some neurobiological processes such as modulation of chemical synaptic transmission (with the highest Metascape value), protein phosphorylation and sensory organ development (Figure 2C). They were also enriched in other biological processes such as regulation of cell projection organization, tube morphogenesis, behavior, ion transport, regulation of secretion, regulation of kinase activity, localization within membrane, response to cytokine, et al. (Ps < 0.01 with FDR correction, detailed see Supplementary Table S4). As to the KEGG pathways, the PLS1 genes were significantly enriched in cancer, cAMP signaling and morphine addiction pathways (Ps < 0.01 with FDR correction).

FIGURE 2

Figure 2. Imaging-transcriptomic analysis treatment-related biomarkers. (A) T-statistic map was constructed for the paired comparison of cortical thickness between Patients_0W and Patients_8W; (B) Partial least squares (PLS) regression analysis showed that expression levels of 2098 genes in the PLS1 were spatially associated with the above brain image t-map; (C) We used Metascape to align the GO biological processes as well as KEGG pathways for the PLS1 genes.

We did not detect significant PLS components associated with the case-control t-statistic map of the CT measures (P > 0.05).

CCA analysis: Correlation between DNAm of PLS1 genes and clinical symptoms

Based on the above detected 2,098 PLS1 genes, we further analyzed which genes' methylation levels were related to patients' clinical symptoms. Combining Laplacian score and CCA analysis, we linked DNAm of 33 representative genes (CCA1 loading genes) with three dimensions of PANSS symptoms, i.e., PANSS_P, _N, and _G. For the multivariate correlations between patients' baseline DNAm of PLS1 genes and baseline clinical behaviors scores, we detected no significant CCA modes (P_FWE > 0.05).

For the analysis between longitudinal alterations of patients' DNAm of PLS1 genes (DNAm_0W–DNAm_8W) and patients' reduction rate of PANSS_P, PANSS_N, and PANSS_G, the first CCA mode (CCA1) was significantly correlated (P_FWE = 0.0174, Figure 3A). To quantify which DNAm variables were significantly involved in the first CCA mode, we then analyzed univariate variable-to-variate Pearson's correlations between the DNAm variables and the corresponding PANSS variate, and found 33 significant genes with P_FDR < 0.05 (see Table 1). The corresponding PLS Z-values of the 33 genes computed in the above imaging-transcriptomic spatial association analysis were also shown in Table 1. There are two overlapping genes of 33 CCA loading genes and the 247 dopamine-related genes (obtained from the PathDIP), i.e., PPP2R2C and PRKX (Figure 3B).

FIGURE 3

Figure 3. Multivariate correlations between DNAm of CCA loading genes and clinical symptoms. (A) We identified significant CCA1 loading genes (33 representative features) associated with patients' longitudinal alterations of clinical symptoms; (B) There are two overlapping genes of 33 CCA loading genes and the 247 dopamine-related genes. CCA, canonical correlation analysis. *, important genes overlapping between dopamine related genes and CCA1 loading genes.

TABLE 1

Table 1. Significant genes in the CCA analysis between longitudinal DNAm of PLS1 genes and patients' reduction rate of clinical symptoms.

Discussion

This study firstly investigates multi-omics biomarkers associated with antipsychotic treatment in the first episode schizophrenia using multiple measures linking CT, public transcriptomic signatures and peripheral epigenetic modifications. As hypothesized, alterations of CT and clinical symptoms after antipsychotic treatment may be transcriptomically and epigenetically underlied. Phenotypic variations of CT related to treatment were spatially associated with expression of genes (PLS1 genes) primarily enriched in neurobiological processes, as well as dopaminergic related and cancer pathways. Moreover, the longitudinal alterations of DNAm for the 33 of PLS1 genes demonstrated significant associations with patients' reduction rate of clinical symptoms.

In this study, following 8 weeks of risperidone treatment, brain cortex showed significant thinning primarily in the frontal and temporal lobes, but also in parts of parietal regions. Similarly, one most recent randomized, double-blind, placebo-controlled study found that patients with major depressive disorder exposed to olanzapine had cortical thinning throughout all lobes across a 36-week period compared with those taking a placebo (31). Cross-sectional studies of patients with schizophrenia undergoing chronic (over 5 years) or current treatment suggest excessive thinning in widespread brain regions such as the middle temporal, prefrontal, occipital and parietal regions, and the excessive cortical thinning, progressing across the entire disease course, appears associated with medication intake (3–7). Contrary to our findings, however, a prospective longitudinal study with a sample of 34 unmedicated patients with schizophrenia demonstrated significant increase of CT after shorter term (6 weeks) of acute-phase treatment with risperidone, and the greater increase were correlated with better clinical response (32). Future studies are needed to investigate the potential different effects of short-term 6- and 8-week atypical treatment on human cerebral CT in first-episode patients.

Interestingly, patients' longitudinal variances of CT after treatment showed significant spatial correlations with gene expression of cortical areas derived from AHBA (PLS1 genes). The detected genes were enriched in multiple biological processes including some neurobiological processes, e.g., modulation of chemical synaptic transmission and protein phosphorylation. Consistently, the action of antipsychotic drugs mainly depends on chemical synaptic transmission of dopaminergic neurons (33); protein phosphorylation alterations have been previously suggested to play a vital role in mediating the treatment effects (34). Together with these discoveries, it is plausible to propose that the decreased CT found in this study could be due to atypical induced synaptic remodeling and neural plasticity documented in the frontal, temporal and parietal cortices. Besides, we also found that genes involving the sensory organ development may underly the treatment-related phenotypic variations in CT. Abnormal processing of auditory (35), visional (36), gustative (37), olfactory (38), and tactile (39) stimuli has been previously suggested to be potential markers of mental disorders including schizophrenia. Our current finding would expand the pipeline of ideas for developing new targets of antipsychotic drug action.

As to the KEGG pathways, the PLS1 genes were significantly enriched in cAMP signaling and morphine addiction and cancer pathways. Previous studies have revealed that proteins downstream of cAMP-related pathway, the main messenger for signaling through the dopamine receptors (40), are altered after acute antipsychotic therapy (41); dopamine receptors, the main targets of antipsychotic drugs, could modulate limbic morphine-induced plasticity (42). It would therefore be reasonable to speculate that longitudinal CT alterations after 8-week risperidone treatment may be driven by transcriptional status of dopaminergic signaling associated pathways. Interestingly, the transcriptional state for the cancer pathway may also underlie the longitudinal phenotypic variances in cortical morphology. Recently, antipsychotic agents are extensively being tested for efficacy in patients with various cancers. For example, chlorpromazine may harbor effects on treating human oral cancer (43). In addition, due to their proven ability to cross the blood–brain barrier, antipsychotic medications were also found to have effect on malignant brain tumors (44). These studies reveal the strategies of using the existing antipsychotic agents for cancer, known as “drug repositioning” or “drug repurposing,” and simultaneously imply that the cancer pathway may serve as a possible target for new antipsychotic drug discovery.

DNAm was found to underlie the known heterogeneity of treatment response. Specifically, for the analysis between longitudinal alterations of patients' DNAm of PLS1 genes and patients' reduction rate of PANSS_P, PANSS_N, and PANSS_G, the first CCA mode was significantly correlated. To quantify which DNAm variables were significantly involved in the first CCA mode, we then analyzed univariate variable-to-variate Pearson's correlations between the DNAm variables and the corresponding PANSS variate, and found 33 significant genes, which overlap with two of the dopamine-related genes, i.e., PPP2R2C and PRKX. Evidence from both transcriptional and epigenetic omics simultaneously demonstrated the importance of dopamine signals in antipsychotic treatment and treatment response.

Our findings must be explained in light of some limitations. First, the sample size is small due to the challenging requirement of drug-naïve patients with first-episode schizophrenia given antipsychotic (risperidone) monotherapy. The current findings are still needed to be replicated in other cohorts with large sample size. Second, we did not calculate the gene expression data (from the AHBA) of the right hemisphere, because not all donors have brain tissues of the right hemisphere. Thus, the association between gene expression and treatment-related alterations in CT does not represent the condition of bilateral cortical regions. Third, the methylation levels were derived from peripheral blood rather than brain tissues of the cortical regions, although DNAm values were significantly correlated between brain and blood (13).

Conclusions

This study firstly revealed that changes of CT and clinical behaviors after treatment may be transcriptionally and epigenetically underlied. Phenotypic variations of cortical thickness related to treatment were spatially associated with expression of genes primarily enriched in neurobiological processes, as well as dopaminergic related and cancer pathways. We define a “three-step” roadmap which represents a vital step toward the exploration of treatment- and treatment response-related biomarkers on the basis of multiple omics rather than a single omics type as a strategy for advancing precise care.

Data availability statement

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.

Ethics statement

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

Author contributions

Conceptualization: XZ, MH, GW, SW, QT, and JZ. Participants recruitment: GW, ZN, SM, LK, and NZ. Image analysis: MH and JZ. Multi-omics analysis: MH and QT. Writing—original draft preparation: XZ. Writing—review and editing: XZ, MH, SW, QT, and JZ. Supervision and project administration: MH, SW, QT, and JZ. All authors contributed to the final version of the paper.

Funding

This research was funded by the National Natural Science Foundation of China (No. 81901357), the National Key Research and Development Program of China (STI2030-Major Projects, 2021ZD0202000), and the College Students' Innovative Entrepreneurial Training Plan Program (No. S202210486266).

Acknowledgments

The authors would like to thank all participants for their time and excellent cooperation.

Conflict of interest

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.

Publisher's note

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2023.1127353/full#supplementary-material

References

1. Oertel-Knöchel V, Knöchel C, Rotarska-Jagiela A, Reinke B, Prvulovic D, Haenschel C, et al. Association between psychotic symptoms and cortical thickness reduction across the schizophrenia spectrum. Cereb Cortex. (2013) 23:61–70. doi: 10.1093/cercor/bhr380

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Takayanagi Y, Sasabayashi D, Takahashi T, Furuichi A, Kido M, Nishikawa Y, et al. Reduced cortical thickness in schizophrenia and schizotypal disorder. Schizophr Bull. (2020) 46:387–94. doi: 10.1093/schbul/sbz051

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lesh TA, Tanase C, Geib BR, Niendam TA, Yoon JH, Minzenberg MJ, et al. A multimodal analysis of antipsychotic effects on brain structure and function in first-episode schizophrenia. JAMA Psychiatry. (2015) 72:226–34. doi: 10.1001/jamapsychiatry.2014.2178

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gjerde PB, Jørgensen KN, Steen NE, Melle I, Andreassen OA, Steen VM, et al. Association between olanzapine treatment and brain cortical thickness and gray/white matter contrast is moderated by cholesterol in psychotic disorders. Psychiatry Res Neuroimag. (2018) 282:55–63. doi: 10.1016/j.pscychresns.2018.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. van Erp TGM, Walton E, Hibar DP, Schmaal L, Jiang W, Glahn DC, et al. Cortical brain abnormalities in 4,474 individuals with schizophrenia and 5098 control subjects via the enhancing neuro imaging genetics through meta-analysis (ENIGMA) consortium. Biol Psychiatry. (2018) 84:644–54. doi: 10.1016/j.biopsych.2018.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Liu N, Xiao Y, Zhang W, Tang B, Zeng J, Hu N, et al. Characteristics of gray matter alterations in never-treated and treated chronic schizophrenia patients. Transl Psychiatry. (2020) 10:136. doi: 10.1038/s41398-020-0828-4

PubMed Abstract | CrossRef Full Text | Google Scholar

7. van Haren NE, Schnack HG, Cahn W, van den Heuvel MP, Lepage C, Collins L, et al. Changes in cortical thickness during the course of illness in schizophrenia. Arch Gen Psychiatry. (2011) 68:871–80. doi: 10.1001/archgenpsychiatry.2011.88

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hawrylycz MJ, Lein ES, Guillozet-Bongaarts AL, Shen EH, Ng L, Miller JA, et al. An anatomically comprehensive atlas of the adult human brain transcriptome. Nature. (2012) 489:391–9. doi: 10.1038/nature11405

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li J, Seidlitz J, Suckling J, Fan F, Ji GJ, Meng Y, et al. Cortical structural differences in major depressive disorder correlate with cell type-specific transcriptional signatures. Nat Commun. (2021) 12:1647. doi: 10.1038/s41467-021-21943-5

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li A, Zalesky A, Yue W, Howes O, Yan H, Liu Y, et al. A neuroimaging biomarker for striatal dysfunction in schizophrenia. Nat Med. (2020) 26:558–65. doi: 10.1038/s41591-020-0793-8

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Magwai T, Shangase KB, Oginga FO, Chiliza B, Mpofana T, Xulu KR, et al. methylation and schizophrenia: current literature and future perspective. Cells. (2021) 10:890. doi: 10.3390/cells10112890

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ovenden ES, McGregor NW, Emsley RA, Warnich L, DNA. methylation and antipsychotic treatment mechanisms in schizophrenia: progress and future directions. Prog Neuropsychopharmacol Biol Psychiatry. (2018) 81:38–49. doi: 10.1016/j.pnpbp.2017.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Braun PR, Han S, Hing B, Nagahama Y, Gaul LN, Heinzman JT, et al. Genome-wide DNA methylation comparison between live human brain and peripheral tissues within individuals. Transl Psychiatry. (2019) 9:47. doi: 10.1038/s41398-019-0376-y

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zong X, Hu M, Pantazatos SP, Mann JJ, Wang G, Liao Y, et al. A dissociation in effects of risperidone monotherapy on functional and anatomical connectivity within the default mode network. Schizophr Bull. (2019) 45:1309–18. doi: 10.1093/schbul/sby175

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zong X, Zhang Q, He C, Huang X, Zhang J, Wang G, et al. DNA methylation basis in the effect of white matter integrity deficits on cognitive impairments and psychopathological symptoms in drug-naive first-episode schizophrenia. Front Psychiatry. (2021) 12:777407. doi: 10.3389/fpsyt.2021.777407

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zong X, He C, Huang X, Xiao J, Li L, Li M, et al. Predictive biomarkers for antipsychotic treatment response in early phase of schizophrenia: multi-omic measures linking subcortical covariant network, transcriptomic signatures, and peripheral epigenetics. Front Neurosci. (2022) 16:853186. doi: 10.3389/fnins.2022.853186

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kay SR, Fiszbein A, Opler LA. The positive and negative syndrome scale (PANSS) for schizophrenia. Schizophr Bull. (1987) 13:261–76. doi: 10.1093/schbul/13.2.261

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Leucht S, Kissling W, Davis JM. The PANSS should be rescaled. Schizophr Bull. (2010) 36:461–2. doi: 10.1093/schbul/sbq016

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Fischl B, Salat DH, Busa E, Albert M, Dieterich M, Haselgrove C, et al. Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain. Neuron. (2002) 33:341–55. doi: 10.1016/S0896-6273(02)00569-X

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Rosen AFG, Roalf DR, Ruparel K, Blake J, Seelaus K, Villa LP, et al. Quantitative assessment of structural image quality. Neuroimage. (2018) 169:407–18. doi: 10.1016/j.neuroimage.2017.12.059

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Morgan SE, Seidlitz J, Whitaker KJ, Romero-Garcia R, Clifton NE, Scarpazza C, et al. Cortical patterning of abnormal morphometric similarity in psychosis is associated with brain expression of schizophrenia-related genes. Proc Natl Acad Sci U S A. (2019) 116:9604–9. doi: 10.1073/pnas.1820754116

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Seidlitz J, Váša F, Shinn M, Romero-Garcia R, Whitaker KJ, Vértes PE, et al. Morphometric similarity networks detect microscale cortical organization and predict inter-individual cognitive variation. Neuron. (2018) 97:231–47.e7. doi: 10.1016/j.neuron.2017.11.039

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Desikan RS, Ségonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage. (2006) 31:968–80. doi: 10.1016/j.neuroimage.2006.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Arnatkeviciute A, Fulcher BD, Fornito A, A. practical guide to linking brain-wide gene expression and neuroimaging data. Neuroimage. (2019) 189:353–67. doi: 10.1016/j.neuroimage.2019.01.011

PubMed Abstract | CrossRef Full Text

25. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yu M, Cullen N, Linn KA, Oathes DJ, Seok D, Cook PA, et al. Structural brain measures linked to clinical phenotypes in major depression replicate across clinical centres. Mol Psychiatry. (2021) 26:2764–75. doi: 10.1038/s41380-021-01039-8

PubMed Abstract | CrossRef Full Text | Google Scholar

27. He X, Deng C, Niyogi P. Laplacian score for feature selection. Adv Neural Inform Syst. (2005) 13:507–14.

28. Smith SM, Nichols TE, Vidaurre D, Winkler AM, Behrens TE, Glasser MF, et al. A positive-negative mode of population covariation links brain connectivity, demographics and behavior. Nat Neurosci. (2015) 18:1565–7. doi: 10.1038/nn.4125

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Nichols TE, Holmes AP. Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum Brain Mapp. (2002) 15:1–25. doi: 10.1002/hbm.1058

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Rahmati S, Abovsky M, Pastrello C, Jurisica I. pathDIP: an annotated resource for known and predicted human gene-pathway associations and pathway enrichment analysis. Nucl Acids Res. (2017) 45:D419–D26. doi: 10.1093/nar/gkw1082

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Voineskos AN, Mulsant BH, Dickie EW, Neufeld NH, Rothschild AJ, Whyte EM, et al. Effects of antipsychotic medication on brain structure in patients with major depressive disorder and psychotic features: neuroimaging findings in the context of a randomized placebo-controlled clinical trial. JAMA Psychiatry. (2020) 77:674–83. doi: 10.1001/jamapsychiatry.2020.0036

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Nelson EA, Kraguljac NV, White DM, Jindal RD, Shin AL, Lahti AC, et al. prospective longitudinal investigation of cortical thickness and gyrification in schizophrenia. Can J Psychiatry. (2020) 65:381–91. doi: 10.1177/0706743720904598

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Amato D, Vernon AC, Papaleo F. Dopamine, the antipsychotic molecule: a perspective on mechanisms underlying antipsychotic response variability. Neurosci Biobehav Rev. (2018) 85:146–59. doi: 10.1016/j.neubiorev.2017.09.027

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Jaros JA, Rahmoune H, Wesseling H, Leweke FM, Ozcan S, Guest PC, et al. Effects of olanzapine on serum protein phosphorylation patterns in patients with schizophrenia. Proteom Clin Appl. (2015) 9:907–16. doi: 10.1002/prca.201400148

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Micoulaud-Franchi JA, Hetrick WP, Aramaki M, Bolbecker A, Boyer L, Ystad S, et al. Do schizophrenia patients with low P50-suppression report more perceptual anomalies with the sensory gating inventory? Schizophr Res. (2014) 157:157–62. doi: 10.1016/j.schres.2014.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Güell F, Bernácer J. Anatomical constitution of sense organs as a marker of mental disorders. Front Behav Neurosci. (2015) 9:59. doi: 10.3389/fnbeh.2015.00059

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Compton MT, Ionescu DF, Broussard B, Cristofaro SL, Johnson S, Haggard PJ, et al. An examination of associations between the inability to taste phenylthiocarbamide (PTC) and clinical characteristics and trait markers in first-episode, nonaffective psychotic disorders. Psychiatry Res. (2013) 209:27–31. doi: 10.1016/j.psychres.2013.03.028

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Moberg PJ, Kamath V, Marchetto DM, Calkins ME, Doty RL, Hahn CG, et al. Meta-analysis of olfactory function in schizophrenia, first-degree family members, and youths at-risk for psychosis. Schizophr Bull. (2014) 40:50–9. doi: 10.1093/schbul/sbt049

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ferri F, Costantini M, Salone A, Di Iorio G, Martinotti G, Chiarelli A, et al. Upcoming tactile events and body ownership in schizophrenia. Schizophr Res. (2014) 152:51–7. doi: 10.1016/j.schres.2013.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Boyd KN, Mailman RB. Dopamine receptor signaling and current and future antipsychotic drugs. Handb Exp Pharmacol. (2012) 68:53–86. doi: 10.1007/978-3-642-25761-2_3

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ahmed MR, Gurevich VV, Dalby KN, Benovic JL, Gurevich EV. Haloperidol and clozapine differentially affect the expression of arrestins, receptor kinases, and extracellular signal-regulated kinase activation. J Pharmacol Exp Ther. (2008) 325:276–83. doi: 10.1124/jpet.107.131987

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Rivera A, Suárez-Boomgaard D, Miguelez C, Valderrama-Carvajal A, Baufreton J, Shumilov K, et al. Dopamine D(4) receptor is a regulator of morphine-induced plasticity in the rat dorsal striatum. Cells. (2021) 11:31. doi: 10.3390/cells11010031

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Jhou AJ, Chang HC, Hung CC, Lin HC, Lee YC, Liu WT, et al. Chlorpromazine, an antipsychotic agent, induces G2/M phase arrest and apoptosis via regulation of the PI3K/AKT/mTOR-mediated autophagy pathways in human oral cancer. Biochem Pharmacol. (2021) 184:114403. doi: 10.1016/j.bcp.2020.114403

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Lee JK, Nam DH, Lee J. Repurposing antipsychotics as glioblastoma therapeutics: potentials and challenges. Oncol Lett. (2016) 11:1281–6. doi: 10.3892/ol.2016.4074

PubMed Abstract | CrossRef Full Text | Google Scholar