Introduction
The psychosis spectrum ranges from serious, enduring, and disabling illness to transient, sub-threshold psychotic experiences in non-clinical populations (Guloksuz and van Os Reference Guloksuz and van Os2018). It represents a wide range of symptoms including aberrant thinking and reasoning, perceptual abnormalities, cognitive disturbance, as well as motivational and social deficits. Consistent with the extended psychosis phenotype model, prevalence is estimated at 5–8% for psychotic experiences in the general population, 3% for clinical psychotic disorders, and 0.5% for arguably the most severe end of the spectrum meeting diagnostic criteria for schizophrenia (van Os et al., Reference van Os, Linscott, Myin-Germeys, Delespaul and Krabbendam2009).
The aetiological and pathophysiological theories of psychosis spectrum have evolved to encompass genetic and environmental factors and their interaction (EUGEI investigators, 2014). The concordance rates between twin pairs suggest the presence of genetic factors with heritability estimates of up to 80% for schizophrenia and 73% for the wider phenotype (Hilker et al., Reference Hilker, Helenius, Fagerlund, Skytthe, Christensen, Werge, Nordentoft and Glenthoj2018). More recent molecular genetic studies have confirmed that schizophrenia spectrum disorder, as a common complex trait, has a polygenic architecture, which is mainly shaped by many common allele variants with small effect sizes that are normally distributed among the general population (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). With the advent of the genome-wide association study approach, the Psychiatric Genomics Consortium has identified 145 significant loci associated with schizophrenia (Pardinas et al., Reference Pardinas, Holmans, Pocklington, Escott-Price, Ripke, Carrera, Legge, Bishop, Cameron, Hamshere, Han, Hubbard, Lynham, Mantripragada, Rees, MacCabe, McCarroll, Baune, Breen, Byrne, Dannlowski, Eley, Hayward, Martin, McIntosh, Plomin, Porteous, Wray, Caballero, Geschwind, Huckins, Ruderfer, Santiago, Sklar, Stahl, Won, Agerbo, Als, Andreassen, Baekvad-Hansen, Mortensen, Pedersen, Borglum, Bybjerg-Grauholm, Djurovic, Durmishi, Pedersen, Golimbet, Grove, Hougaard, Mattheisen, Molden, Mors, Nordentoft, Pejovic-Milovancevic, Sigurdsson, Silagadze, Hansen, Stefansson, Stefansson, Steinberg, Tosato, Werge, Collier, Rujescu, Kirov, Owen, O'Donovan and Walters2018). It is now possible to calculate an individual score summarising the level of genetic risk for schizophrenia, known as polygenic risk score for schizophrenia (PRS-SCZ) (Pardinas et al., Reference Pardinas, Holmans, Pocklington, Escott-Price, Ripke, Carrera, Legge, Bishop, Cameron, Hamshere, Han, Hubbard, Lynham, Mantripragada, Rees, MacCabe, McCarroll, Baune, Breen, Byrne, Dannlowski, Eley, Hayward, Martin, McIntosh, Plomin, Porteous, Wray, Caballero, Geschwind, Huckins, Ruderfer, Santiago, Sklar, Stahl, Won, Agerbo, Als, Andreassen, Baekvad-Hansen, Mortensen, Pedersen, Borglum, Bybjerg-Grauholm, Djurovic, Durmishi, Pedersen, Golimbet, Grove, Hougaard, Mattheisen, Molden, Mors, Nordentoft, Pejovic-Milovancevic, Sigurdsson, Silagadze, Hansen, Stefansson, Stefansson, Steinberg, Tosato, Werge, Collier, Rujescu, Kirov, Owen, O'Donovan and Walters2018).
Similarly, several environmental exposures have been associated with a schizophrenia spectrum disorder, such as childhood adversities, cannabis use, urbanicity, migration, ethnic minorities, hearing impairment, and perinatal factors (Linszen et al., Reference Linszen, Brouwer, Heringa and Sommer2016; Radua et al., Reference Radua, Ramella-Cravaro, Ioannidis, Reichenberg, Phiphopthatsanee, Amir, Yenn Thoo, Oliver, Davies, Morgan, McGuire, Murray and Fusar-Poli2018; Stilo and Murray, Reference Stilo and Murray2019). In accordance with the diathesis-stress model, there is evidence supporting gene–environment interaction in the aetiology of schizophrenia (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019) and mood disorders (Geoffroy et al., Reference Geoffroy, Etain and Houenou2013; Colodro-Conde et al., Reference Colodro-Conde, Couvy-Duchesne, Zhu, Coventry, Byrne, Gordon, Wright, Montgomery, Madden and Ripke2018; Arnau-Soler et al., Reference Arnau-Soler, Adams, Clarke, MacIntyre, Milburn, Navrady, Hayward, McIntosh and Thomson2019a, Reference Arnau-Soler, Macdonald-Dunlop, Adams, Clarke, MacIntyre, Milburn, Navrady, Hayward, McIntosh and Thomson2019b). A recent case-control study found evidence for additive interactions between molecular genetic liability for schizophrenia (i.e. PRS-SCZ) and emotional abuse, emotional neglect, sexual abuse, bullying, and regular cannabis use, suggesting that a multitude of environmental factors and PRS-SCZ are independently and jointly associated with schizophrenia (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019).
To better accommodate the multiplicity of exposures associated with schizophrenia (Guloksuz et al., Reference Guloksuz, Rutten, Pries, Ten Have, de Graaf, van Dorsselaer, Klingenberg, van Os and Ioannidis2018), a cumulative environmental exposure score for schizophrenia – exposome score for schizophrenia (ES-SCZ) – was recently designed and validated through predictive modelling approaches in training and validation data sets of two independent cohorts that followed identical measurement methods for environmental exposures (Pries et al., Reference Pries, Lage-Castellanos, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasoglu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten, van Os and Guloksuz2019). This summary measure is generated using weighted coefficients derived from a single model to take into account the interdependency of exposures. Therefore, ES-SCZ prevents overestimation of the weights per exposure that are likely to occur when correlations between exposures are ignored, e.g. weighted estimates of individual exposures from meta-analyses or simple summation of exposures. Recent studies indicate that the ES-SCZ is associated with psychosis risk states (Guloksuz et al., Reference Guloksuz, Pries, Ten Have, de Graaf, van Dorsselaer, Klingenberg, Bak, Lin, van Eijk, Delespaul, van Amelsvoort, Luykx, Rutten and van Os2020) as well as mental and physical health (Pries et al., Reference Pries, van Os, Ten Have, de Graaf, van Dorsselaer, Bak, Lin, van Eijk, Kenis, Richards, O'Donovan, Luykx, Rutten and Guloksuz2020b) in the general population.
By leveraging aggregate scores of genetic (PRS-SCZ) and environmental liability (ES-SCZ) in the current study and in accordance with a previous study (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019), we aimed to test gene–environment interaction across the psychosis spectrum in a multinational multicentre sample of patients diagnosed with schizophrenia spectrum disorder, their siblings, and healthy comparison participants.
Methods
Study population
Data were derived from the Workpackage 6 (WP6) of the European Network of National Networks studying Gene–Environment Interactions in Schizophrenia (EUGEI) and the Genetic Risk and Outcome for Psychosis (GROUP) studies, collected using uniform assessment schedules between 2010 and 2015 in the Netherlands, Turkey, Spain, and Serbia (Korver et al., Reference Korver, Quee, Boos, Simons and de Haan2012). Both projects were approved by the Medical Ethics Committees of all participating sites and conducted in accordance with the Declaration of Helsinki. All respondents provided written informed consent and, in the case of minors, such a consent was also obtained from parents or legal guardians. Patients were diagnosed with schizophrenia spectrum disorders according to the DSM-IV-TR (average duration of illness since the age of the first contact with mental health services = 9.9 years). Unrelated controls with no lifetime psychotic disorder were recruited from the same population as the cases. Exclusion criteria for all participants were a diagnosis of psychotic disorder due to another medical condition, a history of head injury with loss of consciousness, and an intelligence quotient <70.
EUGEI WP6 (‘vulnerability and severity’) was a cross-sectional study specifically conducted to investigate the role of gene–environment interaction of the vulnerability and severity of schizophrenia spectrum disorder and its intermediate phenotypes in a family-based setting.
GROUP is a naturalistic longitudinal cohort study that started in 2004 in the Netherlands and Dutch-speaking part of Belgium and collected data at baseline, 3 and 6 years follow-ups over an approximate 10-year period, with the aim of studying the interplay of genetic and environmental factors impacting vulnerability and resilience in psychotic disorders. Individuals in the sibling group who manifested lifetime psychotic disorder over the study period were reassigned to the patient group.
Further details of the GROUP and EUGEI projects are provided elsewhere (Korver et al., Reference Korver, Quee, Boos, Simons and de Haan2012; EUGEI investigators 2014). The current analyses used a merged data set of GROUP baseline data and EUGEI WP6 cross-sectional data including 1699 patients, 1753 siblings, and 1542 unrelated healthy comparison participants who were of Caucasian white ethnic origin and had available genotype data.
Outcomes
Diagnosis of schizophrenia spectrum disorder
Patients were diagnosed with schizophrenia spectrum disorders according to the DSM-IV-TR. The diagnosis was confirmed by the Operational Criteria Checklist for Psychotic and Affective Illness (McGuffin et al., Reference McGuffin, Farmer and Harvey1991) in EUGEI WP6, and by the Schedules for Clinical Assessment in Neuropsychiatry (Wing et al., Reference Wing, Babor, Brugha, Burke, Cooper, Giel, Jablenski, Regier and Sartorius1990) and the Comprehensive Assessment of Symptoms and History (Andreasen et al., Reference Andreasen, Flaum and Arndt1992) in GROUP.
Schizotypy trait
In both GROUP and EUGEI, the Structured Interview for Schizotypy-Revised (SIS-R) was administered to siblings and healthy comparison participants. The SIS-R is a semi-structured interview containing 20 schizotypal symptoms and 11 schizotypal signs rated on a four-point scale (Kendler et al., Reference Kendler, Lieberman and Walsh1989; Vollema and Ormel, Reference Vollema and Ormel2000). Symptoms are defined as verbal responses to standardised questions concerning, for example, magical ideation, illusions, and referential thinking. Signs refer to behaviours that are rated by the interviewer such as goal-directedness of thinking and flatness of effect. Questions and rating procedures are standardised. Guided by previous research, 31 item scores were reduced a priori to two-dimensional scores representing the means of seven positive schizotypy items (covering the areas of referential thinking, psychotic phenomena, derealisation, magical ideation, illusions, and suspiciousness) and eight negative/disorganised schizotypy items (covering the areas of social isolation, sensitivity, introversion, restricted affect, disturbances in associative and goal-directed thinking, poverty of speech, and eccentric behaviour) (van Os et al., Reference van Os, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasoglu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and Guloksuz2020).
Genetic and environmental liability measures
Exposome score for schizophrenia
The exposome score in the current analyses was calculated based on our previously validated estimates (Pries et al., Reference Pries, Lage-Castellanos, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasoglu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten, van Os and Guloksuz2019) for constructing cumulative environmental load in this data set. Using the log odds from our previous report, we generated the ES-SCZ by summing log-odds weighted environmental exposures (each exposure defined as absent = ‘0’ and present = ‘1’) including cannabis use, hearing impairment, winter-birth, and childhood adversity domains (emotional and physical neglect, emotional, physical and sexual abuse, and bullying). The definition of each exposure conformed to previous work in this data set.
Childhood adversities were assessed using the Childhood Trauma Questionnaire (CTQ) Short Form (Bernstein et al., Reference Bernstein, Stein, Newcomb, Walker, Pogge, Ahluvalia, Stokes, Handelsman, Medrano, Desmond and Zule2003). This consists of 28 items, rated on a five-point Likert scale, measuring five domains of maltreatment (emotional and physical neglect; emotional, physical, and sexual abuse). The psychometric characteristics of the translated versions (Spanish, Turkish, Dutch, and Serbian) of the CTQ have been comprehensively studied (Sar et al., Reference Sar, Akyuz, Kundakci, Kiziltan and Dogan2004; Thombs et al., Reference Thombs, Bernstein, Lobbestael and Arntz2009; Hernandez et al., Reference Hernandez, Gallardo-Pujol, Pereda, Arntz, Bernstein, Gaviria, Labad, Valero and Gutierrez-Zotes2013). To dichotomise each childhood adversity domain (0 = ‘absent’ and 1 = ‘present’), consistent with previous work in the EUGEI (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019), we used the following cut-off scores for each domain: ⩾9 for emotional abuse; ⩾8 for physical abuse; ⩾6 for sexual abuse; ⩾10 for emotional neglect; and ⩾8 for physical neglect.
Cannabis use was assessed by a modified version of the Cannabis Experiences Questionnaire (Barkus et al., Reference Barkus, Stirling, Hopkins and Lewis2006) in the EUGEI WP6 (0 = ‘none’; 1 = ‘only once or twice’; 2 = ‘a few times a year’; 3 = ‘a few times a month’; 4 = ‘once or more a week’; 5 = ‘everyday’), and by the L section of the Composite International Diagnostic Interview (Robins et al., Reference Robins, Wing, Wittchen, Helzer, Babor, Burke, Farmer, Jablenski, Pickens, Regier, Sartorius and Towle1988) in the GROUP (0 = ‘none’; 1 = ‘less than weekly’; 2 = ‘weekly’; 3 = ‘daily’). Consistent with previous work (van Winkel et al., Reference van Winkel, van Beveren and Simons2011; Pries et al., Reference Pries, Guloksuz, Ten Have, de Graaf, van Dorsselaer, Gunther, Rauschenberg, Reininghaus, Radhakrishnan, Bak, Rutten and van Os2018; Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019; Radhakrishnan et al., Reference Radhakrishnan, Guloksuz, Ten Have, de Graaf, van Dorsselaer, Gunther, Rauschenberg, Reininghaus, Pries, Bak and van Os2019), a binary regular cannabis use variable was constructed by using the cut-off value of one or more per week during the lifetime period of most frequent use.
In accordance with previous studies investigating the association between season of birth and schizophrenia in the Northern hemisphere sites (Davies et al., Reference Davies, Welham, Chant, Torrey and McGrath2003), the high-risk birth period was defined based on the winter solstice (December–March), and a binary winter-birth exposure was constructed. Hearing impairment was defined based on self-reported hearing impairment in the last 12 months (0 = ‘absent’ and 1 = ‘present’).
The history of bullying by peers (emotional, psychological or physical violence) before 17 years of age was assessed using the short version of the Retrospective Bullying Questionnaire (Hunter et al., Reference Hunter, Mora-Merchan and Ortega2004; Schäfer et al., Reference Schäfer, Korn, Smith, Hunter, Mora-Merchán, Singer and Meulen2004) that measures the severity of the bullying experience: 0 = ‘none’; 1 = ‘some (no physical injuries)’; 2 = ‘moderate (minor injuries or transient emotional reactions)’; 3 = ‘marked (severe and frequent physical or psychological harm)’. Exposure to childhood bullying was dichotomised using ⩾1 as the cut-off point (0 = ‘absent’ and ⩾1 = ‘present’).
Polygenic risk score for schizophrenia
Samples of all individuals were genotyped at Cardiff University Institute of Psychological Medicine and Clinical Neurology, using a custom Illumina HumanCoreExome-24 BeadChip genotyping arrays containing probes for 5 70 038 genetic variants (Illumina, San Diego, CA). Genotype data were called using the GenomeStudio package and transferred into PLINK format for further analysis. Quality control was conducted in PLINK v1.07 (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007) or with custom Perl scripts. Variants with a call rate <98% were excluded from the data set. Hardy–Weinberg equilibrium p value was calculated separately in Turkish, northern European, and southern European samples. Variants with Hardy–Weinberg equilibrium p value <1 × 10−6 in any of these three regions were excluded from the data set. After QC, 5 59 505 variants remained. Samples with a call rate <98% were excluded from the data set. A linkage disequilibrium pruned set of variants was calculated using the – indep-pairwise command in PLINK (maximum r 2 = 0.25, window size = 500 single nucleotide polymorphisms (SNPs), window step size = 50 SNPs) and used for further analyses. Homozygosity F values were calculated using the – het command in PLINK, and outlier samples (F < −0.11 or F > 0.15) were excluded. The genotypic sex of samples was calculated from X chromosome data using the check-sex command in PLINK, and samples with different genotypic sex to their database sex were excluded. Identity-by-descent values were calculated for the sample in PLINK. Samples with one or more siblings among the genotyped samples according to the database but no identified genotypic siblings (defined as p̂ >0.35 and <0.65) were excluded. After these were removed from consideration, samples with two or more siblings in the database that were not supported by the genotypic data were also excluded. After visually observing the clustering of errors by genotyping chips, we decided to exclude chips with a high proportion of errors. All samples on chips with five or more sample exclusions due to heterozygosity or call rate (out of 12 possible samples) were excluded. All samples on chips with four or more sample exclusions due to sex or relative checks were also excluded unless their identity was corroborated by concordance between database and genotype relatedness data with a sample on another chip. Principal components (PCs) were calculated in PLINK using linkage disequilibrium (LD) pruned variants after combining the data set with the Thousand Genomes reference data set. After quality control, genotypes were imputed on the Michigan Imputation Server using the Haplotype Reference Consortium reference panel (version 1.1) and the programmes Eagle for haplotype phasing and Minimac3 for imputation (Das et al., Reference Das, Forer, Schonherr, Sidore, Locke, Kwong, Vrieze, Chew, Levy, McGue, Schlessinger, Stambolian, Loh, Iacono, Swaroop, Scott, Cucca, Kronenberg, Boehnke, Abecasis and Fuchsberger2016; Loh et al., Reference Loh, Danecek, Palamara, Fuchsberger, Reshef, Finucane, Schoenherr, Forer, McCarthy, Abecasis, Durbin and Price2016). After imputation, variants with an imputation r 2 > 0.6, minor allele frequency (MAF) > 0.1% and call rate >99% were retained (82 77 535 variants). Best-guess genotypes were generated from genotype probabilities using PLINK.
PRS-SCZ was constructed using summary statistics from the Psychiatric Genomics Consortium (PGC2) genome-wide association study in both samples (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). There was no overlap between the PGC2 and the current data sets. Clumping was performed in imputed best-guess genotypes for each data set using PLINK (maximum r 2 = 0.2, window size = 500 kb, minimum MAF = 10%, minimum INFO score = 0.7), and variants within regions of long-range LD around the genome (including the major histocompatibility complex) excluded (Price et al., Reference Price, Weale, Patterson, Myers, Need, Shianna, Ge, Rotter, Torres, Taylor, Goldstein and Reich2008). PRS-SCZ were then constructed from best-guess genotypes using PLINK at 10 different p-value thresholds (PT = 1, 0.5, 0.3, 0.2, 0.1, 0.05, 0.01, 1 × 10−4, 1 × 10−6, 5 × 10−8). Consistent with previous research in the field (Allardyce et al., Reference Allardyce, Leonenko, Hamshere, Pardinas, Forty, Knott, Gordon-Smith, Porteous, Haywood, Di Florio, Jones, McIntosh, Owen, Holmans, Walters, Craddock, Jones, O'Donovan and Escott-Price2018; Sorensen et al., Reference Sorensen, Debost, Agerbo, Benros, McGrath, Mortensen, Ranning, Hjorthoj, Mors, Nordentoft and Petersen2018) and previous work in this data set, we used PT = 0.05 for our primary analysis, as this threshold optimally captures liability to the disorder in the Psychiatric Genomics Consortium analysis (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014).
Statistical analysis
Stata software version 15.0 was used for the analysis (StataCorp, 2017). Supplementary Table S1 reports missing data. The analyses were conducted on both multiple imputed data and raw data. Under the assumption of missing at random, the multiple imputation chained equation (Royston and White, Reference Royston and White2011) was applied with 20 imputations restricted to in-range values (relative efficiency ⩾ 99%). ES-SCZ was calculated after imputing missing values of the environmental exposures (cannabis use, hearing impairment, winter-birth, and childhood adversity domains). All the analyses were run on multiple imputed data and pooled using Rubin's rules (Rubin, Reference Rubin2004). To test gene–environment interaction, additive models were chosen over multiplicative models prior to data collection (EUGEI consortium meeting, 14 December 2013), consistent with previous work (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019), and given that they provide a superior representation of biological synergy (Rothman, Reference Rothman1976) and inform public health decisions within the sufficient cause framework (Rothman et al., Reference Rothman, Greenland and Walker1980; Kendler and Gardner, Reference Kendler and Gardner2010). For all analyses, random intercept multilevel mixed regression models, taking into account the clustering of participants within countries, were applied. Models including PRS-SCZ were a priori adjusted for ancestry using 10 PCs and adjusted models included age and sex as covariates. The nominal significance threshold was set to p = 0.05.
For the case-control analyses, as utilised in previous studies (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019; Guloksuz et al., Reference Guloksuz, Pries, Ten Have, de Graaf, van Dorsselaer, Klingenberg, Bak, Lin, van Eijk, Delespaul, van Amelsvoort, Luykx, Rutten and van Os2020), ES-SCZ and PRS-SCZ were dichotomised at the quartile cut-off points based on the control distribution within each country (to account for differences between countries that may arise due to ethnic and geographical variation). The highest quartile was considered the binary risk state for schizophrenia (hereafter PRS-SCZ75 and ES-SCZ75). Multilevel logistic regression models were applied to test the independent and joint effects of PRS-SCZ75 and ES-SCZ75 (independent variables) on the diagnosis of schizophrenia (i.e. case-control status; dependent variable). Departure from additivity was tested using the relative excess risk due to interaction (RERI) (Knol and VanderWeele, Reference Knol and VanderWeele2012; VanderWeele and Knol, Reference VanderWeele and Knol2014). RERI greater than zero was defined as a positive deviation from additivity and considered significant when the 95% confidence interval (CI) did not contain zero. Conforming to early work in this sample, we applied the delta method to calculate the RERI using the odds ratios derived from the model (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019). Furthermore, sensitivity analyses were conducted using the bootstrap percentile method to estimate additive interaction between continuous PRS-SCZ and ES-SCZ in unimputed data (N = 1000 bootstrap replications) (Richardson and Kaufman, Reference Richardson and Kaufman2009).
In unaffected siblings and healthy comparison participants, the effects of continuous measures of PRS-SCZ, ES-SCZ, and their interaction on continuous measures of schizotypy dimensions (total, positive, and negative) as dependent variables were tested with multilevel linear regression models, where the coefficient of the product term (PRS-SCZ×ES-SCZ) reflects the departure from additivity (Knol et al., Reference Knol, van der Tweel, Grobbee, Numans and Geerlings2007).
Previous analyses did not indicate a gene–environment correlation between the individual environmental exposures and PRS-SCZ75 in the control sample (Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019). Furthermore, for the current analyses, we tested gene–environment correlation between the continuous (ES-SCZ and PRS-SCZ) and dichotomised (ES-SCZ75 and PRS-SCZ75) exposome and genetic risk scores applying multilevel linear and logistic regression, respectively. Nagelkerke's R2 was calculated based on logistic regression with case-control status as the dependent variable.
Results
Sample demographic data, SIS-R scores, PRS-SCZ75 and ES-SCZ75 distributions are reported in Table 1. Missing data are reported in the Supplementary material (Table S1).
Sample characteristics
ES-SCZ75, exposome score for schizophrenia (75% cut-point); n, number of individuals; PRS-SCZ75, polygenic risk score for schizophrenia (75% cut-point); s.d., standard deviation; SIS-R, the structured interview for schizotypy – revised.
a Within siblings and control groups (3295 individuals).
PRS-SCZ explained 15% of the variance in case-control status (OR = 1.30 [95% CI 1.25, 1.34], p < 0.001) and 20% after adjusting for age, sex, and country (OR = 1.30 [95% CI 1.26, 1.35], p < 0.001). ES-SCZ explained 28% of the variance in case-control status (OR = 2.52 [95% CI 2.29, 2.78], p < 0.001) and 33% after adjusting for age, sex, and country (OR = 2.40 [95% CI 2.17, 2.66], p < 0.001).
There was no evidence for gene–environment correlation, as PRS-SCZ75 was not strongly or significantly associated with ES-SCZ75 in the control group (OR = 1.08 [95% CI 0.78, 1.51], p = 0.635), neither after adjusting for age and sex (OR = 1.08 [95% CI 0.78, 1.51], p = 0.638) nor when using the continuous scores; PRS-SCZ and ES-SCZ (B = −0.008 [95% CI −0.028, 0.013], p = 0.478; adjusted B = −0.008 [95% CI −0.029, 0.012], p = 0.429).
Main and joint effects of PRS-SCZ75 and ES-SCZ75 on case-control status
PRS-SCZ75 was associated with case status (OR = 2.91 [95% CI 2.48, 3.40], p < 0.001; adjusted for age and sex: OR = 2.85 [95% CI 2.43, 3.35], p < 0.001); and ES-SCZ75 was associated with case status (OR = 4.99 [95% CI 4.22, 5.90], p < 0.001, adjusted for age and sex: OR = 4.90 [95% CI 4.14, 5.81], p < 0.001). There was evidence for a positive additive interaction between PRS-SCZ75 and ES-SCZ75 (RERI = 7.29 [95% CI 3.73, 10.85], p < 0.001), also after adjusting for age and sex (Table 2 and Fig. 1). Sensitivity analyses using continuous PRS-SCZ and ES-SCZ confirmed G × E interaction (RERI = 1.77 [95% CI 1.00, 3.24], p = 0.003; adjusted RERI = 1.80 [95% CI 1.01, 3.32], p = 0.004). Results from the analyses using unimputed data corroborated these results and are reported in the Supplementary material (Table S2 and Figure S3).
Additive effect of the polygenic risk score for schizophrenia, 75% cut-point (PRS-SCZ75), and the exposome score for schizophrenia, 75% cut-point (ES-SCZ75) on case-control status, adjusted for age, sex, and ten PCs; RERI: relative excess risk due to interaction.
Interaction of PRS-SCZ75 and ES-SCZ75 on case-control status
CI, confidence interval; ES-SCZ75, exposome score for schizophrenia (75% cut-point); PRS-SCZ75, polygenic risk score for schizophrenia (75% cut-point); RERI, relative excess risk due to interaction.
Adjusted for sex, age, and ten PCs.
Main and joint effects of continuous PRS-SCZ and ES-SCZ on SIS-R dimensions
PRS-SCZ was significantly associated with the SIS-R dimensions in the unaffected sibling/healthy comparison participants sample (total: B = 0.011 [95% CI 0.006, 0.015], p < 0.001; positive: B = 0.012 [95% CI 0.007, 0.018], p < 0.001; negative: B = 0.010 [95% CI 0.005, 0.014], p < 0.001) also after adjusting for age and sex (Table 3). ES-SCZ was also significantly associated with the SIS-R dimensions (total: B = 0.088 [95% CI 0.078, 0.098], p < 0.001; positive: B = 0.103 [95% CI 0.090, 0.116], p < 0.001; negative: B = 0.074 [95% CI 0.064, 0.085], p < 0.001), also after adjusting for age and sex (Table 3). There was evidence for a significant interaction between ES-SCZ and PRS-SCZ on the SIS-R dimensions (total: B = 0.006 [95% CI 0.003, 0.009], p < 0.001; positive: B = 0.005 [95% CI 0.002, 0.009], p = 0.002; and negative: B = 0.006 [95% CI 0.003, 0.009], p < 0.001), also after adjusting for age and sex (Table 3). Results from the analyses in unimputed data confirmed the results in imputed data and are reported in the Supplementary material (Table S4).
Main and joint effects of PRS-SCZ and ES-SCZ on SIS-R scores
B, regression coefficient from the multilevel model; CI, confidence interval; ES-SCZ, exposome score for schizophrenia; PRS-SCZ, polygenic risk score for schizophrenia; SIS-R, the structured interview for schizotypy – revised.
All analyses were adjusted for age and sex.
a Additionally adjusted for ten PCs.
Discussion
To the best of our knowledge, this is the first study testing the role of gene–environment interaction using aggregate scores of environmental and genetic liability across the spectrum of psychosis expression. In the case-control design, we found evidence for additive interaction between PRS-SCZ and ES-SCZ increasing the odds for schizophrenia. Similarly, evidence emerged for interaction between PRS-SCZ and ES-SCZ on schizotypal traits when investigating G×E interaction in the group of unaffected siblings and healthy comparison participants.
By using aggregate scores for genetic and environmental liability for schizophrenia, we provided further support for the role of gene–environment interaction in schizophrenia spectrum disorder (Bernardo et al., Reference Bernardo, Bioque, Cabrera, Lobo, Gonzalez-Pinto, Pina, Corripio, Sanjuan, Mane, Castro-Fornieles, Vieta, Arango, Mezquida, Gasso, Parellada, Saiz-Ruiz, Cuesta and Mas2017; Guloksuz et al., Reference Guloksuz, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasog Lu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and van Os2019) and replicated recent findings of suggestive, but not nominally statistically significant, additive interaction between PRS-SCZ and environmental risk score for schizophrenia in a first episode psychosis cohort (Mas et al., Reference Mas, Boloc, Rodriguez, Mezquida, Amoretti, Cuesta, Gonzalez-Penas, Garcia-Alcon, Lobo, Gonzalez-Pinto, Corripio, Vieta, Castro-Fornieles, Mane, Saiz-Ruiz, Gasso, Bioque and Bernardo2020). When PRS-SCZ75 and ES-SCZ75 were analysed as binary modes of risk factors, the relative excess risk due to the interaction was 6.79 and the corresponding 95% CI was above 2, suggesting a ‘mechanistic’ interaction, which means that the risk of developing schizophrenia for some individuals exists only when both genetic and environmental risks are present together but not when either genetic or environmental risk is present alone. The results further suggest that the PRS-SCZ and ES-SCZ explain 15 and 28% of the variance in case-control samples, respectively.
In a previous study, we demonstrated that the extent of sub-threshold phenotypic expression of schizophrenia polygenic risk is contingent on having a sibling with a psychotic disorder, suggesting a gene–environment interaction underlying schizotypy expression (van Os et al., Reference van Os, Pries, Delespaul, Kenis, Luykx, Lin, Richards, Akdede, Binbay, Altinyazar, Yalincetin, Gumus-Akay, Cihan, Soygur, Ulas, Cankurtaran, Kaymak, Mihaljevic, Petrovic, Mirjanic, Bernardo, Cabrera, Bobes, Saiz, Garcia-Portilla, Sanjuan, Aguilar, Santos, Jimenez-Lopez, Arrojo, Carracedo, Lopez, Gonzalez-Penas, Parellada, Maric, Atbasoglu, Ucok, Alptekin, Saka, Arango, O'Donovan, Rutten and Guloksuz2020). In the light of this new evidence, we tested for the first time the putative role of gene–environment interaction in schizotypy traits. In line with our previous inference, we have now demonstrated that molecular genetic liability for schizophrenia moderates the effect of environmental liability for schizophrenia on phenotypic expression of overall, positive, and negative schizotypy traits in unaffected participants. Although much research has investigated the role of familial sensitivity to individual environmental exposures (e.g. cannabis use and childhood adversities) underlying subclinical psychosis expression (Modinos et al., Reference Modinos, Iyegbe, Prata, Rivera, Kempton, Valmaggia, Sham, van Os and McGuire2013; EUGEI investigators, 2014), only a few studies have utilised PRS to investigate the role of G×E in intermediate psychotic phenotypes (Ronald and Pain, Reference Ronald and Pain2018). A recent study from the 1966 Northern Finland Birth Cohort showed that high birth weight, a risk factor for familial schizophrenia in this cohort, increased the association between PRS-SCZ and social anhedonia, suggesting a gene–environment interaction (Liuhanen et al., Reference Liuhanen, Suvisaari, Kajantie, Miettunen, Sarin, Jarvelin, Lonnqvist, Veijola and Paunio2018). Similarly, a study conducted in a general population twin cohort demonstrated that while PRS-SCZ was not independently associated with affective dysregulation and psychosis proneness, PRS-SCZ increased sensitivity to the effect of childhood adversities on affective dysregulation and psychosis proneness (Pries et al., Reference Pries, Klingenberg, Menne-Lothmann, Decoster, van Winkel, Collip, Delespaul, De Hert, Derom, Thiery, Jacobs, Wichers, Cinar, Lin, Luykx, Rutten, van Os and Guloksuz2020a). Although not a direct test of gene–environment interaction, a study of healthy young males assessed during their compulsory military service showed that there was a negative association between PRS-SCZ and positive schizotypy at military induction (stressful condition) but not at follow-up, providing further support for the key role of environment in the phenotypic expression of schizotypy traits (Hatzimanolis et al., Reference Hatzimanolis, Avramopoulos, Arking, Moes, Bhatnagar, Lencz, Malhotra, Giakoumaki, Roussos, Smyrnis, Bitsios and Stefanis2018). Taken together, although warranting further replication in independent cohorts, these findings imply that the phenotypic expression of schizotypical traits involves underlying genomic liability for schizophrenia that operates, at least in part, through sensitising individuals to the exposome.
The major strengths related to the study population were threefold: sufficient sample size to detect gene–environment interactions, access to comprehensive genotype, phenotype, and exposure data collected through validated interviews conducted by trained psychiatrists, psychologists or research assistants, and the geographical and cultural diversity of the sample that may increase the variation of environmental exposures and thereby provide increased power and replicability to detect interaction effects across populations (Ritz et al., Reference Ritz, Chatterjee, Garcia-Closas, Gauderman, Pierce, Kraft, Tanner, Mechanic and McAllister2017).
The ES-SCZ was constructed using predictive modelling that mutually adjusted for the interdependency of exposures to prevent overestimation of the weights per exposure. ES-SCZ was fully compatible with this study population and clearly outperformed other aggregate scores that were based on meta-analytical estimates or simple summation of exposures as shown previously. Notwithstanding, ES-SCZ was limited by the degree to which exposures were available in the data set, and therefore did not include other exposures that might be of importance, such as obstetric and pregnancy complications (Garcia-Rizo and Bitanihirwe, Reference Garcia-Rizo and Bitanihirwe2020). Furthermore, childhood adversities and cannabis use were retrospectively assessed and may be affected by recall bias (Baldwin et al., Reference Baldwin, Reuben, Newbury and Danese2019). Although population stratification was controlled using PCs and no gene–environment correlation was detected, unmeasured environmental confounding might still be present. The cross-sectional design did not allow for investigating the dynamic nature of gene–environment interaction over time.
In conclusion, we have shown that the interplay between exposome load and genetic liability for schizophrenia contributes to the phenotypical expression of psychosis across the extended phenotype. Our findings provide further empirical support to the notion of aetiological continuity of psychosis spectrum and pave the way for future longitudinal studies of schizotypy traits in the general population. As some individuals may only develop schizophrenia spectrum disorder if both the genetic and environmental vulnerabilities are present, the current findings highlight the importance of the combined effect of genomic and exposomic liability for clinical practice. Furthermore, the results suggest that health care strategies may benefit from focusing on modifiable environmental factors.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S2045796020000943.
Data
The data that support the findings of this study are available from the corresponding author upon reasonable request under the condition of the approval of the EUGEI and GROUP steering committees.
Acknowledgements
The authors are grateful to the patients and their families for participating in the project. They also thank all research personnel involved in the GROUP project, in particular J. van Baaren, E. Veermans, G. Driessen, T. Driesen, E. van't Hag and J. de Nijs. All the DNA samples from Turkey were provided by the Ankara University Brain Research Center Biobank, which was supported by Ankara University Scientific Research Projects Coordination Unit (project no. 10A6055003, 2010).
Author contributions
SG and L-KP conceived the idea of this study and developed the plan for analysis. SG, L-KP, and GAdF performed the statistical analysis and wrote the first draft. SG, JvO, and BPFR provided supervision and expert knowledge on gene–environment interaction. ALR and MoD provided expert knowledge on psychiatric genetics and processed genotyping data. SG and L-KP contributed to data cleaning and database management for initial use and later reuse. All authors contributed to the collection of data and interpretation of the results, revised the manuscript and approved the final version.
Financial support
The EUGEI project was supported by the European Community's Seventh Framework Program under grant agreement no. HEALTH-F2-2009-241909 (Project EU-GEI). Dr O'Donovan is supported by MRC programme grant (G08005009) and an MRC Centre grant (MR/L010305/1). Dr Rutten was funded by a VIDI award number 91718336 from the Netherlands Scientific Organisation. Drs Guloksuz and van Os are supported by the Ophelia research project, ZonMw grant number: 636340001. Dr Arango was supported by the Spanish Ministry of Science and Innovation; Instituto de Salud Carlos III (SAM16PE07CP1, PI16/02012, PI19/024); CIBERSAM; Madrid Regional Government (B2017/BMD-3740 AGES-CM-2); Fundación Familia Alonso and Fundación Alicia Koplowitz.
Conflict of interest
Celso Arango has been a consultant to or has received honoraria or grants from Acadia, Angelini, Gedeon Richter, Janssen Cilag, Lundbeck, Minerva, Otsuka, Roche, Sage, Servier, Shire, Schering Plough, Sumitomo Dainippon Pharma, Sunovion, and Takeda. Michael O'Donovan is supported by a collaborative research grant from Takeda Pharmaceuticals.
Ethical standards
The projects were approved by the Medical Ethics Committees of all participating sites and conducted in accordance with the Declaration of Helsinki. All respondents provided written informed consent and, in the case of minors, such consent was also obtained from parents or legal guardians.
Appendix
Genetic Risk and Outcome of Psychosis (GROUP) Investigators in EUGEI (GROUP-EUGEI) investigators are: Behrooz Z. Alizadeha, Therese van Amelsvoortb, Richard Bruggemana, Wiepke Cahnc,d, Lieuwe de Haane, Bart P. F. Ruttenb, Jurjen J. Luykxc,f,g, Jim van Osc,b,h and Ruud van Winkelb, i
aUniversity of Groningen, University Medical Center Groningen, University Center for Psychiatry, Rob Giel Research Center, Groningen, The Netherlands; bMaastricht University Medical Center, Department of Psychiatry and Neuropsychology, School for Mental Health and Neuroscience, Maastricht, The Netherlands; cUniversity Medical Center Utrecht, Department of Psychiatry, UMC Utrecht Brain Centre, Utrecht University, Utrecht, The Netherlands; dAltrecht, General Menthal Health Care, Utrecht, The Netherlands; eAmsterdam UMC, University of Amsterdam, Department of Psychiatry, Amsterdam, The Netherlands; fGGNet Mental Health, Apeldoorn, The Netherlands; gDepartment of Translational Neuroscience, UMC Utrecht Brain Center, University Medical Center Utrecht, Utrecht University, Utrecht, The Netherlands; hKing's College London, King's Health Partners, Department of Psychosis Studies, Institute of Psychiatry, London, UK and iKU Leuven, Department of Neuroscience, Research Group Psychiatry, Leuven, Belgium
Open access
