Examining the independent and joint effects of genomic and exposomic liabilities for schizophrenia across the psychosis spectrum.

AIMS: Psychosis spectrum disorder has a complex pathoetiology characterised by interacting environmental and genetic vulnerabilities. The present study aims to investigate the role of gene-environment interaction using aggregate scores of genetic (polygenic risk score for schizophrenia (PRS-SCZ)) and environment liability for schizophrenia (exposome score for schizophrenia (ES-SCZ)) across the psychosis continuum.
METHODS: The sample consisted of 1699 patients, 1753 unaffected siblings, and 1542 healthy comparison participants. The Structured Interview for Schizotypy-Revised (SIS-R) was administered to analyse scores of total, positive, and negative schizotypy in siblings and healthy comparison participants. The PRS-SCZ was trained using the Psychiatric Genomics Consortiums results and the ES-SCZ was calculated guided by the approach validated in a previous report in the current data set. Regression models were applied to test the independent and joint effects of PRS-SCZ and ES-SCZ (adjusted for age, sex, and ancestry using 10 principal components).
RESULTS: Both genetic and environmental vulnerability were associated with case-control status. Furthermore, there was evidence for additive interaction between binary modes of PRS-SCZ and ES-SCZ (above 75% of the control distribution) increasing the odds for schizophrenia spectrum diagnosis (relative excess risk due to interaction = 6.79, [95% confidential interval (CI) 3.32, 10.26], p < 0.001). Sensitivity analyses using continuous PRS-SCZ and ES-SCZ confirmed gene-environment interaction (relative excess risk due to interaction = 1.80 [95% CI 1.01, 3.32], p = 0.004). In siblings and healthy comparison participants, PRS-SCZ and ES-SCZ were associated with all SIS-R dimensions and evidence was found for an interaction between PRS-SCZ and ES-SCZ on the total (B = 0.006 [95% CI 0.003, 0.009], p < 0.001), positive (B = 0.006 [95% CI, 0.002, 0.009], p = 0.002), and negative (B = 0.006, [95% CI 0.004, 0.009], p < 0.001) schizotypy dimensions.
CONCLUSIONS: The interplay between exposome load and schizophrenia genetic liability contributing to psychosis across the spectrum of expression provide further empirical support to the notion of aetiological continuity underlying an extended psychosis phenotype.

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 2018). 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., 2009).

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

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., 2016; Radua et al., 2018; Stilo and Murray, 2019). In accordance with the diathesis-stress model, there is evidence supporting gene–environment interaction in the aetiology of schizophrenia (Guloksuz et al., 2019) and mood disorders (Geoffroy et al., 2013; Colodro-Conde et al., 2018; Arnau-Soler et al., 2019, 2019). 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., 2019).

To better accommodate the multiplicity of exposures associated with schizophrenia (Guloksuz et al., 2018), 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., 2019). 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., 2020) as well as mental and physical health (Pries et al., 2020b) 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., 2019), 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., 2012). 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., 2012; 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., 1991) in EUGEI WP6, and by the Schedules for Clinical Assessment in Neuropsychiatry (Wing et al., 1990) and the Comprehensive Assessment of Symptoms and History (Andreasen et al., 1992) 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., 1989; Vollema and Ormel, 2000). 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., 2020).

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., 2019) 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., 2003). 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., 2004; Thombs et al., 2009; Hernandez et al., 2013). To dichotomise each childhood adversity domain (0 = ‘absent’ and 1 = ‘present’), consistent with previous work in the EUGEI (Guloksuz et al., 2019), 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., 2006) 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., 1988) in the GROUP (0 = ‘none’; 1 = ‘less than weekly’; 2 = ‘weekly’; 3 = ‘daily’). Consistent with previous work (van Winkel et al., 2011; Pries et al., 2018; Guloksuz et al., 2019; Radhakrishnan et al., 2019), 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., 2003), 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., 2004; Schäfer et al., 2004) 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., 2007) 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 r2 = 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., 2016; Loh et al., 2016). After imputation, variants with an imputation r2 > 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 r2 = 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., 2008). 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., 2018; Sorensen et al., 2018) 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, 2011) 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, 2004). 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., 2019), and given that they provide a superior representation of biological synergy (Rothman, 1976) and inform public health decisions within the sufficient cause framework (Rothman et al., 1980; Kendler and Gardner, 2010). 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., 2019; Guloksuz et al., 2020), 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, 2012; VanderWeele and Knol, 2014). 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., 2019). 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, 2009).

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

Previous analyses did not indicate a gene–environment correlation between the individual environmental exposures and PRS-SCZ75 in the control sample (Guloksuz et al., 2019). 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).

Table 1.

Sample characteristics

		Patients n = 1699	Siblings n = 1753	Controls n = 1542	Total N = 4994
Age	Mean (s.d.)	31.49 (8.95)	31.73 (9.62)	33.45 (10.62)	32.18 (9.76)
Sex	% female	30.02	53.17	50.58	44.49
SIS-R total	Mean (s.d.)	–	0.39 (0.33)	0.23 (0.24)	0.32 (0.30)a
SIS-R positive	Mean (s.d.)	–	0.41 (0.42)	0.24 (0.31)	0.33 (0.38)a
SIS-R negative	Mean (s.d.)	–	0.38 (0.34)	0.23 (0.24)	0.31 (0.31)a
PRS-SCZ75	% > 75th	45.44	33.09	23.74	34.40
ES-SCZ75	% > 75th	58.79	36.30	21.88	37.64

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.

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

Table 2.

Interaction of PRS-SCZ75 and ES-SCZ75 on case-control status

	PRS-SCZ75 = 0	PRS-SCZ75 = 1	RERI (95% CI)
	Odds ratio (95% CI)	Odds ratio (95% CI)
ES-SCZ75 = 0	1.0	2.79 (2.24–3.47)P < 0.001	6.79 (3.32–10.26)P < 0.001
ES-SCZ75 = 1	4.86 (3.92–6.02)P < 0.001	13.44 (10.21–17.69)P < 0.001

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.

Fig. 1.

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.

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

Table 3.

Main and joint effects of PRS-SCZ and ES-SCZ on SIS-R scores

Psychopathology measures	Main effect PRS-SCZa			Main effect ES-SCZ			Interactiona
	B	95% CI	p-value	B	95% CI	p-value	B	95% CI	p-value
SIS-R total	0.011	0.007–0.015	<0.001	0.088	0.078–0.098	<0.001	0.006	0.003–0.009	<0.001
SIS-R positive	0.012	0.007–0.018	<0.001	0.103	0.091–0.116	<0.001	0.006	0.002–0.009	0.002
SIS-R negative	0.010	0.005–0.014	<0.001	0.074	0.064–0.085	<0.001	0.006	0.004–0.009	<0.001

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.

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., 2017; Guloksuz et al., 2019) 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., 2020). 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., 2020). 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., 2013; EUGEI investigators, 2014), only a few studies have utilised PRS to investigate the role of G×E in intermediate psychotic phenotypes (Ronald and Pain, 2018). 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., 2018). 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., 2020). 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., 2018). 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., 2017).

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, 2020). Furthermore, childhood adversities and cannabis use were retrospectively assessed and may be affected by recall bias (Baldwin et al., 2019). 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.