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.

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.

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.

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.

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

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

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.

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.

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.

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.

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.

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.

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 high-dimensional 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, and 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.”

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

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.

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

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

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