Translational Psychiatry www.nature.com/tp ARTICLE OPEN A longitudinal study of gene expression in first-episode schizophrenia; exploring relapse mechanisms by co-expression analysis in peripheral blood P. Gassó 1,2,41, N. Rodríguez 1,2,41, A. Martínez-Pinteño1,2, G. Mezquida1,2,3,4, M. Ribeiro5,6, J. González-Peñas 4,7, I. Zorrilla4,8,9,10, L. Martínez-Sadurni11, R. Rodriguez-Jimenez4,12,13, I. Corripio4,14,15, S. Sarró 4,16,17, A. Ibáñez 4,18, J. Usall19, A. Lobo 4,20,21, ✉ ✉ C. Moren22,23, M. J. Cuesta5,6, M. Parellada4,7, A. González-Pinto4,8,9,10, E. Berrocoso24,25, M. Bernardo 2,3,4,26 , S. Mas 1,2,4 and 2EPs Group* © The Author(s) 2021 Little is known about the pathophysiological mechanisms of relapse in first-episode schizophrenia, which limits the study of potential biomarkers. To explore relapse mechanisms and identify potential biomarkers for relapse prediction, we analyzed gene expression in peripheral blood in a cohort of first-episode schizophrenia patients with less than 5 years of evolution who had been evaluated over a 3-year follow-up period. A total of 91 participants of the 2EPs project formed the sample for baseline gene 1234567890();,: expression analysis. Of these, 67 provided biological samples at follow-up (36 after 3 years and 31 at relapse). Gene expression was assessed using the Clariom S Human Array. Weighted gene co-expression network analysis was applied to identify modules of co- expressed genes and to analyze their preservation after 3 years of follow-up or at relapse. Among the 25 modules identified, one module was semi-conserved at relapse (DarkTurquoise) and was enriched with risk genes for schizophrenia, showing a dysregulation of the TCF4 gene network in the module. Two modules were semi-conserved both at relapse and after 3 years of follow-up (DarkRed and DarkGrey) and were found to be biologically associated with protein modification and protein location processes. Higher expression of DarkRed genes was associated with higher risk of suffering a relapse and early appearance of relapse (p = 0.045). Our findings suggest that a dysregulation of the TCF4 network could be an important step in the biological process that leads to relapse and suggest that genes related to the ubiquitin proteosome system could be potential biomarkers of relapse. Translational Psychiatry (2021)11:539 ; https://doi.org/10.1038/s41398-021-01645-8 INTRODUCTION been associated with progressive loss of cortical tissue and overall Despite the efficacy of antipsychotic medication in the treatment brain volume [3]. The frequency and the cumulative effects of of schizophrenia-spectrum disorders, this condition is still relapse are related to poor long-term outcomes, a decline in characterized by persistent functional impairment and recurrent quality of life and lower response to subsequent treatment, and it psychotic relapses for most patients. Relapse is characterized by also often requires inpatient hospitalization, increasing healthcare acute psychotic exacerbation and is considered a biologically toxic costs [4]. event that leads to progressive illness course and an overall Around 41–63% of patients will experience a relapse within deterioration of psychosocial functioning [1, 2]. This, in turn, has the first 3 years of a first episode of psychosis (FEP) [5] and up to 1 Department of Clinical Foundations, Pharmacology Unit, University of Barcelona, Barcelona, Spain. 2Institut d’investigacions Biomèdiques August Pi i Sunyer (IDIBAPs), Barcelona, Spain. 3Barcelona Clínic Schizophrenia Unit (BCSU), Neuroscience Institute, Hospital Clínic de Barcelona, Barcelona, Spain. 4Centro de Investigación Biomédica en Red en Salud Mental (CIBERSAM), Madrid, Spain. 5Department of Psychiatry, Complejo Hospitalario de Navarra, Pamplona, Spain. 6IdiSNA, Navarra Institute for Health Research, Pamplona, Spain. 7Department of Child and Adolescent Psychiatry, Institute of Psychiatry and Mental Health, Hospital General Universitario Gregorio Marañón, IiSGM, School of Medicine, Universidad Complutense, Madrid, Spain. 8Department of Psychiatry, Hospital Universitario de Alava, Vitoria, Spain. 9BIOARABA Health Research Institute, Vitoria, Spain. 10 University of the Basque Country, Vitoria, Spain. 11Hospital del Mar Medicar Research Institute (IMIM), Barcelona, Spain. 12Instituto de Investigación Sanitaria Hospital 12 de Octubre (imas12), Madrid, Spain. 13CogPsy Group, Universidad Complutense de Madrid (UCM), Madrid, Spain. 14Psychiatry Department, Institut d’Investigació Biomèdica-Sant Pau (IIB-SANT PAU), Hospital de la Santa Creu i Sant Pau, Barcelona, Spain. 15Universitat Autònoma de Barcelona (UAB), Barcelona, Spain. 16FIDMAG Germanes Hospitalàries Research Foundation, Barcelona, Spain. 17School of Medicine, Universitat Internacional de Catalunya, Barcelona, Spain. 18Department of Psychiatry, Hospital Universitario Ramón y Cajal, IRYCIS, Universidad de Alcalá, Madrid, Spain. 19Etiopatogènia i tractament dels trastorns mentals greus (MERITT) Institut de Recerca Sant Joan de Déu Parc Sanitari Sant Joan de Déu, Barcelona, Spain. 20Department of Medicine and Psychiatry, Universidad de Zaragoza, Zaragoza, Spain. 21Instituto de Investigación Sanitaria Aragón (IIS Aragón), Zaragoza, Spain. 22Cellex, IDIBAPS, University of Barcelona-Hospital Clínic of Barcelona, Barcelona 08036, Spain. 23Centro de Investigación Biomédica en Red (CIBER) de Enfermedades Raras (CIBERER), Madrid 28029, Spain. 24Neuropsychopharmacology and Psychobiology Research Group, Department of Psychology, University of Cádiz, Cádiz, Spain. 25Instituto de Investigación e Innovación Biomédica de Cádiz, INiBICA, Hospital Universitario Puerta del Mar, Cádiz, Spain. 26Department of Medicine, University of Barcelona, Barcelona, Spain. 41These authors contributed equally: P. Gassó, N. Rodríguez. 42A list of authors and their affiliations appears at the end of the paper. ✉email: bernardo@clinic.cat; sergimash@ub.edu Received: 3 May 2021 Revised: 22 September 2021 Accepted: 30 September 2021 P. Gassó et al. 2 80% within the first 5 years of remission from a FEP [6]. The occurred within the last 5 years); (d) not having relapsed after the first main risk factor for relapse is non-adherence to antipsychotic psychotic episode; (e) speaking Spanish correctly; and (f) providing the medication [7, 8]. However, 20–30% of individuals complying signed informed consent form. The exclusion criteria were: (a) having with treatment will develop subsequent relapses [9, 10]. experienced a traumatic brain injury with loss of consciousness; (b) Considering that relapse represents much of the personal and presenting mental retardation understood not only as IQ < 70, but also presenting malfunctioning and problems with adaptive process, and/or social burden of schizophrenia, there is a need to develop more (c) presenting somatic pathology with mental repercussion. efficient interventions or treatment recommendations to prevent From the initial 223 patients recruited in the 2EPs Project, 91 (40.8%) relapse. In this scenario, the search for the potential role of participated in the biological module and provided a biological sample for biomarkers to stratify patients by their individual characteristics gene expression analysis at baseline. Of these, 67 (73.6%) provided according to their risk of relapse is gaining attention. It has been biological samples at follow-up (36 after 3 years and 31 at relapse). proposed that stratification could be undertaken after clinical The study was approved by the investigation ethics committees of all remission with antipsychotics [11]. participating clinical centers. Informed consent was obtained from all Nevertheless, the lack of understanding about the pathophy- participants. For children under the age of 18 years old, parents or legal guardians gave written informed consent before the beginning of their siological mechanisms of relapse hinders research into potential participation in the study, and patients assented to participate. When biomarkers and the development of preventive and therapeutic requested, participants in the study were given a report on the results of interventions. Although different mechanisms have been pro- the tests. This study was conducted in accordance with the Declaration of posed, it appears that dopamine represents the final common Helsinki. pathway in relapse [12]. This agrees with the proposed dopamine super sensitivity theory, wherein the chronic blockade of Clinical assessment dopamine receptors results in compensatory receptor upregula- Demographic data were collected for all patients through semi-structured tion, which in turn leads to an increased sensitivity to endogenous interviews. Diagnoses were determined according to the DSM-IV criteria dopamine that leaves individuals more prone to relapse upon (American Psychiatric Association, 1994), with the SCID-I [26] or the Kiddie- antipsychotic discontinuation [13, 14]. The worsening of long-term SADS [27] depending on age. outcomes after chronic exposure to antipsychotics is generating Clinical symptomatology was assessed using the Spanish validated again an intense debate [15, 16]. Several moderators that version of the Positive and Negative Syndrome Scale (PANSS) [28]. influence neurodevelopment, such as stress [17], may impact this Pharmacological treatment was also recorded during all the visits. The process, and several potentially state-dependent biomarkers have extent of medication non-adherence was assessed using the Morisky Green Levine Medication Adherence Scale [29]. been proposed such as inflammatory markers [18, 19] or neurotrophins [20–22]. To uncover relapse mechanisms, we investigated a molecular Sample collection, RNA isolation, and microarray phenotype such as gene expression as an intermediate measure hybridization between genetic and clinical variation. To this end, we analyzed Ten milliliter of peripheral blood was collected at baseline and at follow- the preservation of blood co-expression modules (clusters of up (3 years or at relapse) in PAXgene Blood RNA tubes (PreAnalytiX genes with highly correlated expression) in a cohort of first- Gmbh, Switzerland). Total RNA was isolated in accordance with the manufacturer’s protocol (PAXgene Blood RNA kit, PreAnalytiX Gmbh, episode schizophrenia patients with less than 5 years of Switzerland). The purity and integrity of RNA was assessed using an evolution evaluated over a 3-year follow-up period. Co- Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). One expression modules were defined at study enrollment with microgram of purified RNA from each of sample was submitted to the patients under remission of their symptoms. Preservation of Kompetenzzentrum für Fluoreszente Bioanalytik Microarray Technology these modules was assessed at relapse or after 3 years of (KFB, BioPark Regensburg GmbH, Regensburg, Germany) for labeling enrollment for those patients who had not experienced relapse. and hybridization to Clariom S Human Array (Affymetrix, Santa Clara, Non-preserved or semi-conserved modules at relapse are CA, USA), following the manufacturer’s protocols. The Clariom S Human expected to be of relevance in its pathophysiology and are Array comprises more than 221,300 probes covering over 337,100 transcripts and variants, which in turn represent 20,800 genes. potential biomarkers of relapse. Genome-wide expression analysis and the WGCNA procedure MATERIAL AND METHODS Microarray data preprocessing was performed using the Oligo R This study is part of the project “Clinical and neurobiological determinants package [30]. The data were standardized using robust multichip of second episodes of schizophrenia. Longitudinal study of first episode of analysis. Multiple probes mapping to the same gene were merged psychosis” (PI11/00325) (2EPs Project), which arose to identify and using the average as the summary of the hybridization values. Co- characterize those clinical, environmental and biological factors that expression modules were identified using the R software package for predict a relapse. weighted gene co-expression network analysis (WGCNA) [31]. Firstly, in order to remove outlier samples, distance-based adjacency matrices of samples were estimated and sample network connectivity according to Study design the distances was standardized. Samples with connectivity less than −5 The 2EPs is a naturalistic, multicenter, coordinated, and multimodal study were considered as outliers and were excluded (Supplementary Fig. of patients with a first psychotic episode of schizophrenia with less than 5 1A). Based on the assumption that non- or semi-preserved modules years of evolution, and takes a 3-year longitudinal-prospective follow-up between remission and relapse may be functionally related to relapse design. The project includes six modules: general, neuroimaging, mechanisms, baseline samples were considered as the reference set for adherence, neurocognition, physical health, and biological. Due to its module detection. The co-expression analysis involved constructing a main goals, the present study was framed within the general and matrix of pairwise correlations between all pairs of genes across all biological modules. The first one assesses the presence or absence of selected samples. Next, the matrix was raised to a soft-thresholding relapses and includes the clinical assessments. The biological module power (β = 6 in this study) to obtain an adjacency matrix (Supplemen- searches for biomarkers potentially involved in second episodes [23]. tary Fig. 1B). To identify modules of co-expressed genes, we constructed the topological overlap-based dissimilarity, which was Subjects then used as input to average linkage hierarchical clustering. This step The inclusion criteria from the 2EPs Project were: (a) age between 16 resulted in a clustering tree (dendrogram) whose branches were and 40 years at the time of first assessment (baseline visit); (b) meeting identified for cutting based on their shape, using the dynamic tree- diagnostic criteria according to DSM-IV for schizophrenia or schizo- cutting algorithm (Supplementary Fig. 1C). The above steps were phreniform disorder [24]; (c) being in remission according to Andreas- performed using the automatic network construction and module sen’s criteria [25] from the first psychotic episode (which should have detection function (blockwiseModules in WGCNA), with the following Translational Psychiatry (2021)11:539 P. Gassó et al. 3 parameters: minModuleSize of 30, reassignThreshold of 0, and RESULTS mergeCutHeight of 0.25. Table 1 shows the demographic, clinical, and pharmacological characteristics of the study participants at baseline (N = 91), at Preservation analysis relapse (N = 31), and after 3 years of follow-up if they had not Module preservation analysis allowed us to evaluate how well the experienced a relapse (N = 36). Non-significant differences in modular structure of the baseline samples is preserved after 36 months baseline variables were observed between the three groups of or at relapse. To do this, according to Langfelder et al. [32], the function patients. As could be expected, significant differences were modulePreservation of WGCNA package was used and a permutation observed in clinical symptomatology at follow-up, with patients test (based on the generation of 200 random permutations) was at relapse exhibiting a significant worsening of their symptoms calculated, which assesses the preservation of the connectivity and with respect to baseline and to patients that had not experienced density of each network module. We used composite module preservation statistics that are constructed to summarize changes in relapse after 3 years of follow-up (post-hoc analysis p < 0.001). module preservation [32]. Zsummary is a composite module preserva- However, baseline symptomatology did not differ between tion statistic that simultaneously assesses whether the genes in a patients who suffer a relapse and patients who had not defined module in the baseline samples remain highly connected after experienced relapse (total PANSS 48.5 ± 14.6 vs. 48.5 ± 14.4, F1 = 36 months (or at relapse) and investigates whether the connectivity 0.935, p = 0.966). Regarding pharmacological treatment at follow- patterns between the genes in the baseline samples remain similar, up, patients who relapsed were treated more frequently with compared with the connectivity after 36 months or relapse samples. olanzapine (p = 0.02), and received anxiolytics (p = 0.007) and The modules with Zsummary > 10 were interpreted as highly pre- lithium (p = 0.016) more often than patients who had not served, if 2 < Zsummary < 10 were defined as semi-conserved and the experienced a relapse. modules with Zsummary < 2 were considered to be non-preserved following Langfelder et al. [32]. The main advantage of this score is that Taking baseline samples (N = 90) as the reference set, we it allowed threshold definition, although it often depends on module identified 25 modules of co-expressed genes (Supplementary Fig. size. We also calculated Zdensity, a density-based preservation statistic 1C). The inferred modules showed different sizes from 41 used to determine whether the genes in a reference module remain (DarkGrey module) to 5627 genes (Turquoise module). Also, highly connected in the test network (Zdensity > 10, density preserva- 1901 genes were assigned to the gray module, which represents tion); and Zconnectivity, a preservation statistic that assesses whether the genes that were not co-expressed based on gene dissimilarity the overall connectivity pattern between genes in a reference module measures and were not assigned to any of the modules. is similar in the reference and test networks (Zconnectivity > 5, Module preservation analysis, using z-score statistics, revealed preserved connectivity) [32]. As a complementary statistic we also two semi-conserved modules after 3 years of follow-up: DarkRed computed medianRank, a composite score based on the ranks of the observed preservation statistics, which shows no dependence on the (53 genes, Zsummary = 5.9, Zdensity = 4.2, Zconnectivity = 7.4) module size. The main disadvantage of medianRank is that it is only and DarkGrey (41 genes, Zsummary = 8.6, Zdensity = 9.1, Zcon- applicable for ranking modules, given relative preservation information nectivity = 7.9) (Fig. 1A). At relapse, module preservation analysis of the modules. identified three semi-conserved modules: DarkRed (53 genes, Network visualization of selected modules showing the correlation Zsummary = 5.9, Zdensity = 4.3, Zconnectivity = 7.5), DarkGrey network adjacencies was performed using Cytoscape 3.8.2 software. (41 genes, Zsummary = 7.4, Zdensity = 6.7, Zconnectivity = 7.9) and DarkTurquoise (48 genes, Zsummary = 8.7, Zdensity = 8.9, Characterization of relevant modules Zconnectivity = 7.5) (Fig. 1B). The analysis of the medianRank To characterize non-preserved or semi-conserved modules, two score confirmed the loss of preservation of the DarkRed and approaches were followed: (1) gene set enrichment analysis of Gene DarkGrey modules in both, non-relapsed (Fig. 1A) and relapsed Ontology terms (Biological Process) and Reactome pathways were patients (Fig. 1B). In contrast, results regarding the semi- assessed using gProfiler R package [33] (“moderate” hierarchical filtering preservation of the DarkTurquoise module in relapsed patients was used and only pathways containing between 10 and 2000 genes should be considered with caution. According to the medianRank, and FDR-adjusted p < 0.05 were accepted); and (2) the genes in each the DarkTurquoise module—not showing high preservation module were used to define a gene set, and each such gene set was tested for overlap with gene sets formed by: compared to other modules with similar size—is far from semi- preserved modules such as DarkRed and DarkGrey, indicating that a. differentially expressed genes in dorsolateral prefrontal cortex Zsummary results for DarkTurquoise module could be influenced (DLPFC) between subjects with schizophrenia (N = 159) versus by module size instead of relapse. control (n = 293) subjects [34]; Since the relapsed and non-relapsed samples at follow-up are b. differentially expressed genes in DLPFC between subjects with subsets of the basal sample, differences in module preservation schizophrenia (N = 258) versus control (n = 279) subjects [35]; may be related to this different sample composition. As a c. genes associated with schizophrenia using gene expression imputa- complementary analysis we performed the preservation analysis tion (transcriptome-wide analysis) across multiple brain regions in using the common samples in the reference and test sets for both, 40,299 schizophrenia cases and 65,264 matched controls [36]; relapsed and non-relapsed sets of patients. This analysis did not d. genes that have association with common genome-wide association study (GWAS) meta-analysis of the CLOZUK and independent shown any difference with the initial analysis (Supplementary Psychiatric Genomic Consortium datasets, excluding related and Fig. 2). overlapping samples (total of 40,675 cases and 64,643 controls) [37]. To investigate the biological functions of the genes included in the semi-conserved modules (Supplementary Table 1), a gene set Gene overlap was assessed using GeneOverlap R package [38]. enrichment analysis was performed. Seventeen biological pro- cesses from Gene Ontology and six Reactome pathways were Statistical procedures found to be enriched in the DarkRed module, mainly related to Differential gene expression of selected genes was assessed using the detection of stimulus and protein ubiquitination processes. The Limma R package [39]. Module eigengenes (the first principal component DarkGrey module was enriched with four Biological processes of each module computed in the WGCNA) at baseline were used to test related to protein location. No enrichment of biological processes the relapse predictive properties of selected modules using nonparametric and Reactome pathways were detected in the DarkTurquoise receiver operating characteristic (ROC) curve analysis, as well as to test the module. Table 2 shows the identified significant terms for each effect of selected modules on the length of time (months from baseline) module. taken to reach relapse using Kaplan–Meier and Cox regression analysis. Among the three semi-conserved modules, only the gene set ROC and survival analysis were assessed using IBM SPSS Statistics 20.0 (statistical analysis software, IBM, Chicago, IL, USA). from the DarkTurquoise module showed significant overlap with gene sets from previous studies (Supplementary Fig. 3). Gene Translational Psychiatry (2021)11:539 P. Gassó et al. 4 Table 1. Demographic, clinical, and pharmacological characteristics of: the 91 participants who provided biological sample at baseline; the 36 participants who provided biological sample at baseline and after 3-year follow-up; and the 31 participants who provided biological sample at baseline and at relapse. Baseline Follow-up Relapse Statistic N 91 36 31 Age, mean ± SD 25.33 ± 5.85 25.47 ± 5.61 24.87 ± 5.40 F2 = 0.1, p = 0.902 Age at first diagnosis, mean ± SD 24.10 ± 5.73 23.83 ± 5.43 23.43 ± 5.27 F2 = 0.0, p = 0.968 Gender, male, n (%) 62 (67.4) 28 (77.7) 20 (64.5) X22 = 1.6, p = 0.447 Ethnicity, Caucasian, N (%) 81 (88.0) 31 (86.1) 30 (96.7) X22 = 2.2, p = 0.323 Symptomatology PANSS positive (mean ± SD) 9.6 ± 3.2 9.8 ± 4.9 20.3 ± 8.1 F2 = 57.0, p < 0.001 PANSS negative (mean ± SD) 13.9 ± 5.4 12.4 ± 5.1 19.3 ± 8.2 F2 = 12.4, p < 0.001 PANSS general (mean ± SD) 25.6 ± 7.9 23.6 ± 8.9 39.3 ± 14.6 F2 = 27.2, p < 0.001 PANSS total (mean ± SD) 49.1 ± 14.7 45.8 ± 17.3 78.8 ± 28.2 F2 = 34.0, p < 0.001 Antipsychotic Aripiprazol, N (%) 34 (40.5) 15 (57.7) 10 (34.5) X22 = 3.3, p = 0.186 Paliperidone, N (%) 24 (28.6) 7 (26.9) 7 (24.1) X22 = 0.2, p = 0.897 Risperidone, N (%) 13 (15.4) 1 (3.8) 5 (17.2) X22 = 2.6, p = 0.262 Olanzapine, N (%) 12 (14.3) 1 (3.8) 9 (31.0) X22 = 7.9, p = 0.02 Clozapine, N (%) 6 (7.1) 6 (23.1) 4 (13.8) X22 = 5.1, p = 0.07 Quetiapine, N (%) 3 (3.6) 0 (0.0) 2 (6.9) X22 = 1.8, p = 0.390 Amisulpride, N (%) 3 (3.6) 0 (0.0) 0 (0.0) X22 = 2.0, p = 0.366 Etamine, N (%) 2 (2.4) 0 (0.0) 0 (0.0) X22 = 1.3, p = 0.514 Asenapine, N (%) 1 (1.2) 0 (0.0) 1 (3.4) X22 = 1.2, p = 0.537 Haloperidol, N (%) 1 (1.2) 0 (0.0) 0 (0.0) X22 = 0.6, p = 0.719 Ziprasidone, N (%) 1 (1.2) 0 (0.0) 1 (3.4) X22 = 1.2, p = 0.537 Co-medication Antidepressant, yes, N (%) 24 (28.6) 7 (26.9) 4 (13.8) X22 = 2.5, p = 0.279 Anxiolytic, yes, N (%) 13 (15.5) 2 (7.7) 11 (37.9) X22 = 9.6, p = 0.007 Lithium, yes, N (%) 2 (2.4) 0 (0.0) 4 (13.8) X22 = 8.2, p = 0.016 Antiepileptic, yes, N (%) 5 (5.9) 2 (7.7) 1 (3.4) X22 = 0.4, p = 0.790 Antiparkinsonian, yes, N (%) 11 (13.1) 1 (3.8) 3 (10.3) X22 = 1.7, p = 0.412 Adherence Morinsky–Green (mean ± SD) 1.3 ± 0.9 1.3 ± 0.7 1.6 ± 1.1 F2 = 1.6, p = 0.259 Non-adherent, N (%) 7 (8.5) 3 (11.1) 6 (20.7) X22 = 3.1, p = 0.218 overlap was detected with genes differentially expressed in DLPFC taken to reach relapse. To this end, baseline module eigengene [35] (OR = 2.7, p = 0.02) and genes associated with schizophrenia values were dichotomized based on percentile distribution, and using gene expression imputation across multiple brain regions the highest 75 percentile was chosen as risk variant for the [36] (OR = 12.8, p = 0.01). Among the overlapped genes (TCF4, analysis. Patients at the highest 75 percentile (N = 20) showed BCL11A, SPRY1) included in the DarkTurquoise, TCF4 was present higher risk of suffering a relapse (OR = 2.10, CI = 1.01–4.33, not only in the studies of Fromer et al. [35] and Huckins et al. [36], beta = 0.742 ± 0.370, p = 0.045), and relapse appeared earlier but also in the GWAS from Pardiñas et al. [37]. Figure 2 shows the (21.60 ± 2.74 vs 27.13 ± 1.85 months, p = 0.045) than the rest of network plots of the direct interactors of TCF4 in the DarkTur- patients (N = 45) (Fig. 3B). quoise module. A loss of density and connectivity could be observed at relapse with respect to baseline, in agreement with the module preservation statistics. DISCUSSION The predictive properties of semi-conserved modules at base- Relapse remains among one of the most challenging aspects of line were tested using ROC curves (Fig. 3A). To this end, module psychosis. In this study, the WGCNA approach was used to identify eigengene values at baseline were used as predictors of relapse. important genes clustered into modules that are likely enriched As could be observed, the DarkRed module showed better for biological pathways that could be dysregulated in relapse. performances (AUC = 0.603, CI = 0.464-0.742), followed by Dark- Identified genes will provide new insights into the molecular Grey (AUC = 0.556, CI = 0.414–0.699) and DarkTurquoise (AUC = mechanism involved in relapse and potential biomarkers for 0.542, CI = 0.395–0.689), although none of the values achieved relapse prediction. significance (p > 0.05). Using a preservation analysis, our study identified two semi- As DarkRed was the module that showed better prediction conserved modules of co-expressed genes (DarkRed and Dark- performance, Kaplan–Meier and Cox regression analysis were used Grey) in both groups after 3 years of follow-up and at relapse. to test its effect on the length of time (months from baseline) Taking into account that these modules could not be considered Translational Psychiatry (2021)11:539 P. Gassó et al. 5 Fig. 1 Module preservation analysis. Scatter plot of the meadianRank and Zsummary composite module preservation statistic and module size in A samples after 3-year follow-up and B samples at relapse. The modules with Zsummary > 10 were interpreted as highly preserved, if 2 < Zsummary < 10 were defined as semi-conserved and the modules with Zsummary < 2 were considered to be non-preserved. state-markers, we hypothesized that genes belonging to these suffering a relapse and earlier appearance of relapse. Interestingly, modules could be influenced by antipsychotic treatment and genes included in the DarkRed module participate in biological course of disease, and therefore could be related to therapeutic processes related to the ubiquitin proteosome system. The response. Considering that therapeutic response could be one of ubiquitin proteosome system regulates protein degradation and the predictors of relapse, we tested the relapse predictive is a master regulator of neuronal development, influencing properties of these modules at baseline. Patients with higher neurogenesis, synaptogenesis, and neurotransmission by deter- expression of genes in the DarkRed module showed higher risk of mining the localization, interaction, and turnover of synaptic Translational Psychiatry (2021)11:539 P. Gassó et al. 6 Table 2. Significant terms (biological processes from Gene Ontology and Reactome pathways) of the gene set enrichment analysis performed in the three semi-conserved modules identified in the present study. Module Term_id Term_name Adjusted p value Term size Term list −9 DarkRed GO:0050911 Detection of chemical stimulus involved in sensory perception 1 × 10 396 13 of smell DarkRed GO:0007600 Sensory perception 8 × 10−7 970 14 DarkRed GO:0007186 G protein-coupled receptor signaling pathway 7 × 10−6 1361 15 DarkRed GO:0016579 Protein deubiquitination 1 × 10−5 291 8 DarkRed GO:0070646 Protein modification by small protein removal 1 × 10−5 309 8 DarkRed GO:0050877 Nervous system process 9 × 10−5 1465 14 DarkRed GO:0003008 System process 1 × 10−2 2327 14 −2 DarkRed GO:0006511 Ubiquitin-dependent protein catabolic process 2 × 10 657 7 DarkRed GO:0019941 Modification-dependent protein catabolic process 2 × 10−2 664 7 DarkRed GO:0043632 Modification-dependent macromolecule catabolic process 2 × 10−2 675 7 DarkRed GO:0051603 Proteolysis involved in cellular protein catabolic process 4 × 10−2 756 7 DarkRed REAC:R-HSA-381753 Olfactory Signaling Pathway 1 × 10−8 380 11 DarkRed REAC:R-HSA-418555 G alpha (s) signaling events 1 × 10−7 522 11 DarkRed REAC:R-HSA-5689880 Ub-specific processing proteases 3 × 10−6 200 7 DarkRed REAC:R-HSA-5688426 Deubiquitination 2 × 10−5 276 7 −5 DarkRed REAC:R-HSA-388396 GPCR downstream signaling 9 × 10 1088 11 DarkRed REAC:R-HSA-372790 Signaling by GPCR 1 × 10−4 1158 11 DarkGrey GO:0042989 Sequestering of actin monomers 2 × 10−3 12 3 DarkGrey GO:0032507 Maintenance of protein location in cell 6 × 10−3 67 4 DarkGrey GO:0045185 Maintenance of protein location 2 × 10−2 98 4 DarkGrey GO:0051651 Maintenance of location in cell 3 × 10−2 226 5 Fig. 2 TCF4 interactome. Network plot of the direct interactor of TCF4 in the DarkTurquoise module at A baseline, B 3-year follow-up, and C relapse. The figure depicts the connectivity patterns (correlation network adjacencies) between genes in the module. proteins [40]. Recently, the ubiquitin proteosome system has been of the leading schizophrenia risk genes [47]—stands out [48]. associated with schizophrenia [41] through gene expression Previous GWAS identified TCF4 polymorphisms linked with studies [42–45] and seems to be affected by antipsychotic schizophrenia and other psychiatric conditions [37, 49–53] as treatment [46]. well as non-neurological genetic diseases [54–56]. Two other Our analysis also identified one module of co-expressed genes genes of the DarkTurquoise module, BLC11A and SPRY1, over- that showed evidence of alteration in intramodular connectivity lapped with genes differentially expressed in previous studies in the relapse samples: the DarkTurquoise module. However, this analyzing brain gene expression [35, 36]. Interestingly, BCL11A, a result could be due to lack of robustness of the modules instead zinc finger protein that regulates transcription, showed decreased of relapse, as seems to be influenced by module size. It is worth co-expression in multiple schizophrenia cohorts, in both per- noticing that the genes included in this module, despite being ipheric and central nervous system, using machine learning to identified in blood samples, showed significant overlap with predict co-expression at the individual level [57]. Moreover, genes associated with schizophrenia in previous studies analyzing genetic variants in this gene were found to be associated with brain gene expression [35, 36]. Among these genes, the member schizophrenia [50] and to co-localize with differential DNA of the helix–loop–helix family transcription factor 4 (TCF4)—one methylation sites [58]. Translational Psychiatry (2021)11:539 P. Gassó et al. 7 Fig. 3 Predictive properties of selected modules. A ROC curve analysis of module eigengenes values at baseline of the three semi-conserved modules (DarkRed, DarkGrey, and DarkTurquoise) identified in the present study. B Cox regression analysis of the length of time (months from baseline) taken to reach relapse using dichotomized DarkRed baseline module eigengene values. Some limitations should be considered in the interpretation 5. Álvarez-Jiménez M, Gleeson JF, Henry LP, Harrigan SM, Harris MG, Killackey E, of our results. Firstly, the sample size that might also limit the et al. Road to full recovery: longitudinal relationship between symptomatic statistical power to detect a difference between groups. remission and psychosocial recovery in first-episode psychosis over 7.5 years. Secondly, although the Morinsky–Green scale was used to Psychol Med. 2012;42:595–606. 6. Bocchio-Chiavetto L, Zanardini R, Tosato S, Ventriglia M, Ferrari C, Bonetto C, et al. assess the adherence of the patients, our sample size is not Immune and metabolic alterations in first episode psychosis (FEP) patients. Brain large enough to stratify patients according to it. Therefore, we Behav Immun. 2018;70:315–24. cannot analyze separately secondary relapse, commonly 7. Kane JM. Improving patient outcomes in schizophrenia: achieving remission, associated with non-adherence, and natural or primary relapse, preventing relapse, and measuring success. J Clin Psychiatry. 2013;74:e18. which represents relapse in the absence of this influencer [12]. 8. García S, et al. Adherence to antipsychotic medication in bipolar disorder and Thirdly, we cannot discard that module preservation changes schizophrenic patients: a systematic review. J Clin Psychopharmacol. during the longitudinal analysis are related to age and normal 2016;36:355–71. development. To rule out the confounding effect of age in our 9. Rubio JM, Schoretsanitis G, John M, Tiihonen J, Taipale H, Guinart D, et al. Psy- results a control group would be necessary, formed by healthy chosis relapse during treatment with long-acting injectable antipsychotics in individuals with schizophrenia-spectrum disorders: an individual participant data volunteers matching the age range of our group with a follow- meta-analysis. Lancet Psychiatry. 2020;7:749–61. up of at least 36 months. Lastly, due to the naturalistic design, 10. Rubio JM, Taipale H, Correll CU, Tanskanen A, Kane JM, Tiihonen J. Psychosis drug treatment was not controlled; the study participants breakthrough on antipsychotic maintenance: results from a nationwide study. maintained their usual treatment. Regarding the slightly Psychol Med. 2020;50:1356–67. differences in drug treatment at follow-up, we cannot rule 11. Rubio JM, Malhotra AK, Kane JM. Towards a framework to develop neuroimaging out the possible effect of these differences in the expression of biomarkers of relapse in schizophrenia. Behav Brain Res. 2021;402:113099. the genes included in semi-preserved modules detected at 12. Remington G, Foussias G, Agid O, Fervaha G, Takeuchi H, Hahn M. The neuro- relapse. Besides these limitations, the strength of this study lies biology of relapse in schizophrenia. Schizophr Res. 2014;152:381–90. in the inclusion of a consistent well-characterized first-episode 13. Muller P, Seeman P. Dopaminergic supersensitivity after neuroleptics: time- course and specificity. Psychopharmacol (Berl). 1978;60:1–11. schizophrenia patient sample in remission due to its natur- 14. Chouinard G, Jones BD. Neuroleptic-induced supersensitivity psychosis: clinical alistic and longitudinal design. and pharmacologic characteristics. Am J Psychiatry. 1980;137:16–21. To our knowledge this is the first study that compares 15. Leucht S, Davis JM. Do antipsychotic drugs lose their efficacy for relapse pre- peripheral gene expression between remission and relapse in a vention over time? Br J Psychiatry. 2017;211:127–9. sample of first-episode schizophrenia patients. Although further 16. Moncrieff J, Steingard S. A critical analysis of recent data on the long-term out- validation in larger samples is still needed, our results provide new come of antipsychotic treatment. Psychol Med. 2019;49:750–3. https://doi.org/ insights into the molecular mechanisms involved in relapse, 10.1017/S0033291718003811. potential candidates for biomarker discovery, and the develop- 17. Myin-Germeys I, Delespaul P, van Os J. Behavioural sensitization to daily life stress ment of preventive and therapeutic interventions. in psychosis. Psychol Med. 2005;35:733–41. 18. Miller BJ, Buckley P, Seabolt W, Mellor A, Kirkpatrick B. Meta-analysis of cytokine alterations in schizophrenia: clinical status and antipsychotic effects. Biol Psy- chiatry. 2011;70:663–71. REFERENCES 19. Miller BJ, Gassama B, Sebastian D, Buckley P, Mellor A. Meta-analysis of lym- 1. Milev P, Ho BC, Arndt S, Andreasen NC. Predictive values of neurocognition and phocytes in schizophrenia: clinical status and antipsychotic effects. Biol Psy- negative symptoms on functional outcome in schizophrenia: a longitudinal first- chiatry. 2013;73:993–9. episode study with 7-year follow-up. Am J Psychiatry. 2005;162:495–506. 20. Buckley PF, Mahadik S, Pillai A, Terry A Jr. Neurotrophins and schizophrenia. 2. Kane JM. Treatment strategies to prevent relapse and encourage remission. J Clin Schizophr Res. 2007;94:1–11. Psychiatry. 2007;68:27–30. 21. Pillai A, Schooler NR, Peter D, Looney SW, Goff DC, Kopelowicz A, et al. Predicting 3. Andreasen NC, Liu D, Ziebell S, Vora A, Ho BC. Relapse duration, treatment relapse in schizophrenia: is BDNF a plausible biological marker? Schizophr Res. intensity, and brain tissue loss in schizophrenia: a prospective longitudinal MRI 2018;193:263–8. study. Am J Psychiatry. 2013;170:609–15. 22. Martinez-Cengotitabengoa M, MacDowell KS, Alberich S, Diaz FJ, Garcia-Bueno B, 4. McCutcheon RA, Pillinger T, Mizuno Y, Montgomery A, Pandian H, Vano L, et al. Rodriguez-Jimenez R, et al. BDNF and NGF signalling in early phases of psychosis: The efficacy and heterogeneity of antipsychotic response in schizophrenia: a relationship with inflammation and response to antipsychotics after 1 year. meta-analysis. Mol Psychiatry. 2021;26:1310–20. Schizophr Bull. 2016;42:142–51. Translational Psychiatry (2021)11:539 P. Gassó et al. 8 23. Bernardo M, Amoretti S, Cuesta MJ, Parellada M, Mezquida G, 2EPs Group. The development and plasticity of cortical and hippocampal neurons. Mol Autism. prevention of relapses in first episodes of schizophrenia: the 2EPs Project, 2018;9:20. background, rationale and study design. Rev Psiquiatr Salud Ment. 49. Schizophrenia Psychiatric Genome-Wide Association Study (GWAS) Consortium. 2021;14:164–76. https://doi.org/10.1016/j.rpsm.2020.07.004. Genome-wide association study identifies five new schizophrenia loci. Nat Genet. 24. American Psychiatric Association. DSM-IV: Diagnostic and Statistical Manual of 2011;43:969–76. Mental Disorders. Washington: American Psychiatric Association; 1994. 50. Schizophrenia Working Group of the Psychiatric Genomics Consortium. Biological 25. Andreasen NC, Carpenter WT Jr, Kane JM, Lasser RA, Marder SR, Weinberger DR. insights from 108 schizophrenia-associated genetic loci. Nature. 2014;511:421–7. Remission in schizophrenia: proposed criteria and rationale for consensus. Am J 51. Wray NR, Ripke S, Mattheisen M, Trzaskowski M, Byrne EM, Abdellaoui A, et al. Psychiatry. 2005;162:441–9. Genome-wide association analyses identify 44 risk variants and refine the genetic 26. Williams JB, Gibbon M, First MB, Spitzer RL, Davies M, Borus J, et al. The Structured architecture of major depression. Nat Genet. 2018;50:668–81. Clinical Interview for DSM-III-R (SCID). II. Multisite test-retest reliability. Arch Gen 52. Cross-Disorder Group of the Psychiatric Genomics Consortium. Identification of Psychiatry. 1992;49:630–6. risk loci with shared effects on five major psychiatric disorders: a genome-wide 27. Kaufman J, Birmaher B, Brent D, Rao U, Flynn C, Moreci P, et al. Schedule for analysis. Lancet. 2013;381:1371–9. affective disorders and schizophrenia for school-age children-present and life- 53. Gelernter J, Sun N, Polimanti R, Pietrzak R, Levey DF, Bryois J, et al. Genome-wide time version (K-SADS-PL): initial reliability and validity data. J Am Acad Child association study of post-traumatic stress disorder reexperiencing symptoms in Adolesc Psychiatry. 1997;36:980–8. >165,000 US veterans. Nat Neurosci. 2019;22:1394–401. 28. Peralta V, Cuesta MJ. Validación de la escala de los síndromes positivo y negativo 54. Baratz KH, Tosakulwong N, Ryu E, Brown WL, Branham K, Chen W, et al. E2-2 (PANSS) en una muestra de esquizofrénicos españoles. Actas Luso Esp Neurol protein and Fuchs’s corneal dystrophy. N Engl J Med. 2010;363:1016–24. Psiquiatr. 1994;22:171–7. 55. Wieben ED, Aleff RA, Tosakulwong N, Butz ML, Highsmith WE, Edwards AO, et al. 29. Unni EJ, Farris KB. Development of a new scale to measure self-reported medi- A common trinucleotide repeat expansion within the transcription factor 4 (TCF4, cation nonadherence. Res Soc Adm Pharm. 2015;11:e133–143. E2-2) gene predicts Fuchs corneal dystrophy. PLoS One. 2012;7:e49083. 30. Carvalho BS, Irizarry RA. A framework for oligonucleotide microarray preproces- 56. Ellinghaus D, Folseraas T, Holm K, Ellinghaus E, Melum E, Balschun T, et al. sing. Bioinformatics. 2010;26:2363–7. Genome-wide association analysis in primary sclerosing cholangitis and ulcera- 31. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network tive colitis identifies risk loci at GPR35 and TCF4. Hepatology. 2013;58:1074–83. analysis. BMC Bioinforma. 2008;9:559. 57. Chen J, Cao H, Kaufmann T, Westlye LT, Tost H, Meyer-Lindenberg A, et al. 32. Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and Identification of reproducible BCL11A alterations in schizophrenia through reproducible? PLoS Comput Biol. 2011;7:e1001057. individual-level prediction of coexpression. Schizophr Bull. 2020;46:1165–71. 33. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. Profiler—a 58. Hannon E, Dempster E, Viana J, Burrage J, Smith AR, Macdonald R, et al. An web server for functional interpretation of gene lists (2016 update). Nucleic Acids integrated genetic-epigenetic analysis of schizophrenia: evidence for co- Res. 2016;44:W83–89. localization of genetic associations and differential DNA methylation. Genome 34. Gandal MJ, Haney JR, Parikshak NN, Leppa V, Ramaswami G, Hartl C, et al. Shared Biol. 2016;17:176. molecular neuropathology across major psychiatric disorders parallels polygenic overlap. Science. 2018;359:693–7. 35. Fromer M, Roussos P, Sieberts SK, Johnson JS, Kavanagh DH, Perumal TM, et al. ACKNOWLEDGEMENTS Gene expression elucidates functional impact of polygenic risk for schizophrenia. This study was supported by the Carlos III Healthcare Institute, the Spanish Ministry of Nat Neurosci. 2016;19:1442–53. Science, Innovation and Universities, the European Regional Development Fund 36. Huckins LM, Dobbyn A, Ruderfer DM, Hoffman G, Wang W, Pardiñas AF, et al. (ERDF/FEDER) (PI08/0208, PI11/00325, PI14/00612); Centro de Investigación Biomé- Gene expression imputation across multiple brain regions provides insights into dica en Red de Salud Mental (CIBERSAM); CERCA Program; Catalan Government, the schizophrenia risk. Nat Genet. 2019;51:659–74. Secretariat of Universities and Research of the Department of Enterprise and 37. Pardiñas AF, Holmans P, Pocklington AJ, Escott-Price V, Ripke S, Carrera N, et al. Knowledge (2017SGR1562 and 2017SGR1355); and Institut de Neurociències, Common schizophrenia alleles are enriched in mutation-intolerant genes and in Universitat de Barcelona. The authors thank the Language Advisory Service at the regions under strong background selection. Nat Genet. 2018;50:381–9. University of Barcelona for manuscript revision. The authors also thank all subjects 38. Shen L, Sinai ISoMaM. GeneOverlap: test and visualize gene overlaps. R package and their families for the time and effort spent on this study as well as Ana Meseguer version 1.26.0. 2020. http://shenlab-sinai.github.io/shenlab-sinai/. for sample collection assistance. 39. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47. 40. Luza S, Opazo CM, Bousman CA, Pantelis C, Bush AI, Everall IP. The ubiquitin AUTHOR CONTRIBUTIONS proteasome system and schizophrenia. Lancet Psychiatry. 2020;7:528–37. The results presented here are part of a broader project, the 2EPs study. MB is the 41. Middleton FA, Mirnics K, Pierri JN, Lewis DA, Levitt P. Gene expression profiling coordinator of the 2EPs study. PG and NF performed the statistical analysis and wrote reveals alterations of specific metabolic pathways in schizophrenia. J Neurosci. the first draft of the manuscript. AM-P performed the sample isolation and 2002;22:2718–29. preparation and participated in the statistical analysis. CM participated in the sample 42. Altar CA, Jurata LW, Charles V, Lemire A, Liu P, Bukhman Y, et al. Deficient preparation. GM participated in the coordination of the sample shipment, the hippocampal neuron expression of proteasome, ubiquitin, and mitochondrial maintenance of the 2EPs database, and in the recruitment and assessment of the genes in multiple schizophrenia cohorts. Biol Psychiatry. 2005;58:85–96. sample. EB is the coordinator of the Biological module of the 2EPs study. MR, JG-P, IZ, 43. Bousman CA, Chana G, Glatt SJ, Chandler SD, May T, Lohr J, et al. Positive LM-S, RR-J, SS, AI, JU, AL, MJC, MP, and AG-P participated in the recruitment and symptoms of psychosis correlate with expression of ubiquitin proteasome genes assessment of the sample. SM designed, supervised, and performed the statistical in peripheral blood. Am J Med Genet B Neuropsychiatr Genet. analysis, performed the interpretation of the results, and wrote the first draft of the 2010;153B:1336–41. manuscript. All the authors, including the 2EPs Group authors listed in the acronym, 44. Arion D, Corradi JP, Tang S, Datta D, Boothe F, He A, et al. Distinctive tran- contributed to the final draft of the manuscript. scriptome alterations of prefrontal pyramidal neurons in schizophrenia and schizoaffective disorder. Mol Psychiatry. 2015;20:1397–405. 45. Hertzberg L, Maggio N, Muler I, Yitzhaky A, Majer M, Haroutunian V, et al. COMPETING INTERESTS Comprehensive gene expression analysis detects global reduction of proteasome A. Ibáñez has received research support from or served as speaker or advisor for subunits in schizophrenia. Schizophr Bull. 2021;47:785–95. https://doi.org/ JanssenCilag, Lundbeck, and Otsuka. J. Saiz-Ruiz has been as speaker for and on the 10.1093/schbul/sbaa160. advisory boards of Adamed, Lundbeck, Servier, Medtronic, Casen Recordati, 46. Seabra G, de Almeida V, Reis-de-Oliveira G, Crunfli F, Antunes ASLM, Martins-de- Neurofarmagen, Otsuka, Indivior, Lilly, Schwabe, Janssen, and Pfizer, outside the Souza D. Ubiquitin-proteasome system, lipid metabolism and DNA damage repair submitted work. M. Bernardo has been a consultant for, received grant/research are triggered by antipsychotic medication in human oligodendrocytes: implica- support and honoraria from, and been on the speakers/advisory board of ABBiotics, tions in schizophrenia. Sci Rep. 2020;10:12655. Adamed, Angelini, Casen Recordati, Janssen-Cilag, Menarini, Rovi, and Takeda. M. 47. Teixeira JR, Szeto RA, Carvalho VMA, Muotri AR, Papes F. Transcription factor 4 Bioque has been a consultant for, received grant/research support and honoraria and its association with psychiatric disorders. Transl Psychiatry. 2021;11:19. from, and been on the speakers/advisory board of has received honoraria from talks 48. Jung M, Häberle BM, Tschaikowsky T, Wittmann MT, Balta EA, Stadler VC, et al. and/or consultancy of Adamed, Angelini, Ferrer, Janssen-Cilag, Lundbeck, Neur- Analysis of the expression pattern of the schizophrenia-risk and intellectual dis- axpharm, Otsuka, Pfizer, and Sanofi, and grants from Spanish Ministry of Health, ability gene TCF4 in the developing and adult brain suggests a role in Instituto de Salud Carlos III. P. A. Saiz has been a consultant to and/or has received Translational Psychiatry (2021)11:539 P. Gassó et al. 9 honoraria or grants from Adamed, CIBERSAM, European Comission, Government of Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims the Principality of Asturias, Instituto de Salud Carlos III, Janssen-Cilag, Lundbeck, in published maps and institutional affiliations. Otsuka, Pfizer, Plan Nacional Sobre Drogas, and Servier. R. Rodriguez-Jimenez has been a consultant for, spoken in activities of, or received grants from Instituto de Salud Carlos III, Fondo de Investigación Sanitaria (FIS), Centro de Investigación Biomédica en Red de Salud Mental (CIBERSAM), Madrid Regional Government (S2010/ BMD-2422 AGES; S2017/BMD-3740), JanssenCilag, Lundbeck, Otsuka, Pfizer, Open Access This article is licensed under a Creative Commons Ferrer, Juste, Takeda, Exeltis, Casen Recordati, and Angelini. Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative ADDITIONAL INFORMATION Commons license, and indicate if changes were made. The images or other third party Supplementary information The online version contains supplementary material material in this article are included in the article’s Creative Commons license, unless available at https://doi.org/10.1038/s41398-021-01645-8. indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory Correspondence and requests for materials should be addressed to 2EPs GroupM. regulation or exceeds the permitted use, you will need to obtain permission directly Bernardo or S. Mas. from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. Reprints and permission information is available at http://www.nature.com/ reprints © The Author(s) 2021 2EPS GROUP M. Bioque2,3,4,26, S. Amoretti1,3,4, A. Andreu-Bernabeu4,7, X. Gurriarán4,7, A. Alonso-Solís4,14,15, E. Grasa4,14,15, P. López4,8,9,10, E. Garcia4,8,9,10, D. Bergé4,27,28, A. Trabsa11, L. Sànchez-Pastor13, O. Jiménez-Rodríguez13, E. Pomarol-Clotet4,16, I. Feria-Raposo16,29, A. Butjosa4,19, M. Pardo30, L. Moreno-Izco5,6, A. M. Sánchez-Torres6, J. Saiz-Ruiz4,18, L. León-Quismondo31, J. Nacher32,33,34, F. Contreras35,36, C. De-la- Cámara4,20,37, M. Gutiérrez4,38,39 and P. A. Sáiz40 27 Hospital del Mar Medical Research Institute (IMIM), Barcelona, Spain. 28Autonomous University of Barcelona, Barcelona, Spain. 29Benito Menni CASM, Sant Boi de Llobregat, Spain. 30Institut de Recerca Sant Joan de Déu, Hospital Sant Joan de Déu, Esplugues de Llobregat, Spain. 31Department of Psychiatry, Hospital Universitario Ramón y Cajal, Madrid, Spain. 32Neurobiology Unit, Program in Neurosciences and Institute of Biotechnology and Biomedicine (BIOTECMED), Universitat de València, Burjassot, Spain. 33 Biomedical Research Networking Centre in Mental Health (CIBERSAM), Madrid, Spain. 34Biomedical Research Institute INCLIVA, Valencia, Spain. 35Bellvitge Biomedical Research Institute IDIBELL, Department of Psychiatry, Bellvitge University Hospital, Hospitalet de Llobregat, Barcelona, Spain. 36Biomedical Research Networking Center for Mental Health Network (CIBERSAM), Barcelona, Spain. 37Hospital Clínico Universitario and Instituto de Investigación Sanitaria (IIS) Aragón, Zaragoza, Spain. 38Araba University Hospital, Bioaraba Research Institute, Psychiatry, Vitoria, Spain. 39Department of Neurosciences, University of the Basque Country (UPV/EHU), Leioa, Spain. 40Department of Psychiatry, School of Medicine, University of Oviedo, Instituto de Investigación Sanitaria del Principado de Asturias, Mental Health Services of Principado de Asturias, Biomedical Research Networking Centre in Mental Health (CIBERSAM), Oviedo, Spain. Translational Psychiatry (2021)11:539