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)(4)(5)(6)(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 R 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.

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.
"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).
Step : CT biomarkers related to schizophrenia and treatment T 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 T W 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 : 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 .
/fpsyt. . 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 tstatistic 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 Ztests. 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.

DNAm of PLS 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: 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 highdimensional datasets, i.e., DNAm values of PLS1 genes. Laplacian score was performed as follows: where L is Laplacian Matrix, D g is the degree matrix of L, andx 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).

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).

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).
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 PLS 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 imagingtranscriptomic 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).

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 .

FIGURE
Imaging-transcriptomic analysis treatment-related biomarkers. (A) T-statistic map was constructed for the paired comparison of cortical thickness between Patients_ W and Patients_ W; (B) Partial least squares (PLS) regression analysis showed that expression levels of genes in the PLS 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 PLS genes.
found that patients with major depressive disorder exposed to olanzapine had cortical thinning throughout all lobes across a 36week period compared with those taking a placebo (31). Crosssectional 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)(4)(5)(6)(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 . /fpsyt. .  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 Frontiers in Psychiatry frontiersin.org . /fpsyt. . 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 "threestep" 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.