Integrated analysis of environmental and genetic influences on cord blood DNA methylation in new-borns

ARTICLE Integrated analysis of environmental and genetic influences on cord blood DNA methylation in new-borns Darina Czamara et al. # Epigenetic processes, including DNA methylation (DNAm), are among the mechanisms allowing integration of genetic and environmental factors to shape cellular function. While many studies have investigated either environmental or genetic contributions to DNAm, few have assessed their integrated effects. Here we examine the relative contributions of prenatal environmental factors and genotype on DNA methylation in neonatal blood at variably methylated regions (VMRs) in 4 independent cohorts (overall n=2365). We use Akaike’s information criterion to test which factors best explain variability of methylation in the cohort-specific VMRs: several prenatal environmental factors (E), genotypes in cis (G), or their additive (G +E) or interaction (GxE) effects. Genetic and environmental factors in combination best explain DNAm at the majority of VMRs. The CpGs best explained by either G, G +E or GxE are functionally distinct. The enrichment of genetic variants from GxE models in GWAS for complex disorders supports their importance for disease risk. https://doi.org/10.1038/s41467-019-10461-0 OPEN Correspondence and requests for materials should be addressed to E.B.B. (email: [email protected]). # A full list of authors and their affiliations appears at the end of the paper. NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 1 1234567890():,; Foetal or prenatal programming describes the process by which environmental events during pregnancy influence the development of the embryo with on-going implications for future health and disease. Several studies have shown that the in utero environment is associated with disease risk, including coronary heart disease1,2, type 2 diabetes3, childhood obesity4,5as well as psychiatric problems6and disorders7–9. Environmental effects on the epigenome, for example, via DNA methylation, could lead to sustained changes in gene transcription and thus provide a molecular mechanism for the enduring influences of the early environment on later health10. Smoking during pregnancy influences widespread and highly reproducible differences in DNA methylation at birth11. Less dramatic effects have been reported for maternal body mass index (BMI)12, preeclampsia and gestational diabetes13,14. Possible epigenetic changes as a consequence of prenatal stress are less well established15. Some of these early differences in DNA methylation persist, although attenuated, through childhood11,16 and might be related to later symptoms and indicators of disease risk, including BMI during childhood17,18 or substance use in adolescence19. These data emphasise the potential importance of the prenatal environment for the establishment of inter-individual variation in the methylome as a predictor or even mediator of disease risk trajectories. In addition to the environment, the genome plays an important role in the regulation of DNA methylation. To this end, the impact of genetic variation, especially of single nucleotide polymorphisms (SNPs) on DNA methylation in different tissues, has resulted in the discovery of a large number of methylation quantitative trait loci (meQTLs, i.e., SNPs significantly associated with DNA methylation status20). These variants are primarily in cis, i.e., at most 1 million base pairs away from the DNA methylation site20–22 and often co-occur with expression QTLs or other regulatory QTLs23–25. The association of meQTLs with DNA methylation is relatively stable throughout the life course21.In addition, SNPs within meQTLs are strongly enriched for genetic variants associated with common disease in large genome-wide association studies (GWAS) such as BMI, inflammatory bowel disease, type 2 diabetes or major depressive disorder21,23,24,26. Environmental and genetic factors may act in an additive or multiplicative manner to shape the epigenome to modulate phenotype presentation and disease risk27. However, very few studies have so far investigated the joint effects of environment and genotype on DNA methylation, especially in a genome-wide context. Klengel et al.28, for instance, reported an interaction of the FK506 binding protein 5 gene (FKBP5) SNP genotype and childhood trauma on FKBP5 methylation levels in peripheral blood cells, with trauma associated changes only observed in carriers of the rare allele. The most comprehensive study of integrated genetic and environmental contributions to DNA methylation so far was performed by Teh et al.29. This study examined variably methylated regions (VMRs), defined as regions of consecutive CpG-sites showing the highest variability across all methylation sites assessed on the Illumina Infinium HumanMethylation450 BeadChip array. In a study of 237 neonate methylomes derived from umbilical cord tissue, the authors explored the proportions of the influence of genotype vs. prenatal environmental factors such as maternal BMI, maternal glucose tolerance and maternal smoking on DNA methylation at VMRs. They found that 75% of the VMRs were best explained by the interaction between genotype and environmental factors (GxE) whereas 25% were best explained by SNP genotype and none by environmental factors alone. Collectively, these studies highlight the importance of investigating the combination of environmental and genetic contributions to DNA methylation and not only their individual contribution. The main objective of the present study is to extend our knowledge of combined effects of prenatal environment and genetic factors on DNA methylation at VMRs. Specifically, this is addressed by: (1) assessing the stability of the best explanatory factors across different cohorts and whether this extends to all environmental factors, (2) dissecting differences between additive and interactive effects of gene and environment not explored in Teh et al., (3) testing whether VMRs influenced by genetic and/or environmental factors might have a different predicted impact on gene regulation and (4) evaluating the relevance of genetic variants that interact with the environment to shape the methylome for their contribution to genetic disease risk. Our results show that across cohorts genetic variants in combination with prenatal environment are the best predictors of variance in DNA methylation. We observe functional differences of both the genetic variants and the methylation sites best explained by genetic or additive and interactive effects of genes and environment. Finally, the enrichment of genetic variants within additive as well as interactive models in GWAS for complex disorders supports the importance of these environmentally modified methylation quantitative trait loci for disease risk. Results Cohorts and analysis plan. We investigated the influence of the prenatal environment and genotype on VMRs in the DNA of 2365 newborns within 4 different cohorts: Prediction and Prevention of Pre-eclampsia and Intrauterine Growth Restrictions (PREDO, cordblood)30, the UCI cohort (refs. 31–33, heel prick), the Drakenstein Child Health Study (DCHS, cordblood)34,35 and the Norwegian Mother and Child Cohort Study (MoBa, cordblood36). A description of the workflow of this manuscript is given in Fig. 1and the details for each of the cohorts are given in Table 1. We analysed 963 cord blood samples from the PREDO cohort with available genome-wide DNA methylation and genotype data. Of these samples, 817 had data on the Illumina 450k array (PREDO I) and 146 on the Illumina EPIC array (PREDO II). The main analyses are reported for PREDO I, and replication and extension of the results is shown for PREDO II as well as for three independent cohorts including 121 heel prick samples (UCI cohort, EPIC array) as well as 258 (DCHS, 450 K and EPIC array) and 1023 cord blood samples (MoBa, 450 K array). We tested 10 different prenatal environmental factors covering a broad spectrum of prenatal phenotypes (see Table 1) (referred to as E), as well as cis SNP genotype (referred to as G), i.e., SNPs located in at most 1MB distance to the specific CpG, additive effects of cis SNP genotype and prenatal environment (G +E) and cis SNP×environment interactions (GxE) for association with DNA methylation levels (see Fig. 1). We then assessed for each VMR independently which model described the variance of DNAm best using Akaike’s information criterion (AIC)37. In all models, we corrected for child’s gender, ethnicity (using MDScomponents), gestational age as well as estimated cell proportions to account for cellular heterogeneity. Variably methylated regions.Wefirst identified candidate VMRs, defined as regions of CpG-sites showing the highest variability across all methylation sites. In PREDO I, we identified 10,452 variable CpGs that clustered into 3982 VMRs (see Supplementary Data 1). Most VMRs (n=2683) include 2 CpGs. As detailed in Supplementary Note 1, the distribution of methylation levels of CpGs within these VMRs is unimodal, (see Supplementary Fig. 1A), VMRs are enriched in specific functional regions of the genome, correlate with differences in gene ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 2NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications expression, and overlap with sites associated with specific prenatal environmental factors. To examine the factors that best explain the variance in methylation in these functionally relevant sites, we chose the CpG-site with the highest MAD score as representative of the VMR. These CpGs are named tagCpGs. The correlation between methylation levels of tagCpG and average methylation of the respective VMR was high (mean r=0.85, sd r=0.08), suggesting that the tag CpGs are valid representatives of their VMRs. Furthermore, tagCpGs are mainly uncorrelated with each other (mean r=0.03, sd =0.12). Which models explain methylation of tagCpGs best? We next compared the fit of four models for each of the 3,982 tagCpGs (see Fig. 1): best SNP (G model), best environment (E model), SNP+environment (G +E model) and SNP× environment (GxE model). Association results for each model are listed in Supplementary Data 2–5. For each tagCpG, the model with the lowest AIC was chosen as the best model (see Methods section). In total, 40.6% of tagCpGs were best explained by GxE (n=1616), followed by G (30%, n=1, 194) and G +E (29%, n=1171) (Fig. 2a). E explained most variance in one tagCpG. All tag CpGs and the respective SNPs and environments from the best model are listed in Supplementary Data 6–8 and Supplementary Table 1. With regard to environmental factors, 27.0% of tagCpGs best explained by the G +E model were associated with environmental factors related with stress or, in particular, glucocorticoids (i.e., maternal betamethasone treatment), 40.8% with general maternal factors (mostly maternal age) and 32.20 % with factors related to metabolism (pre-pregnancy BMI, hypertension, gestational diabetes). For best model GxE tagCpGs, the proportions of environmental factors were similar with 22.2, 44.1 and 33.7%, respectively (see Fig. 2b). We next looked into the delta AIC, i.e., the difference between the AIC of the best model to the AIC of the next best model (see Supplementary Note 2). GxE models appear to be winning by a significantly larger AIC margin over the next best model, when compared to the other types of winning models (see Fig. 2c). DeepSEA prediction of SNP function. We were next interested in understanding the functionality of both the VMRs as well as the associated SNPs in the G, GxE and G +E models. For this we restricted the analyses only to potentially functional relevant SNPs using DeepSEA38 and not all linkage disequilibrium (LD)- pruned SNPs as described above. DeepSEA, a deep neural network pretrained with DNase-seq and ChIP-seq data from the ENCODE39 project, predicts the presence of histone marks, DNase hypersensitive regions (DHS) or TF binding for a given 1 kb sequence. The likelihood that a specific genetic variant influences regulatory chromatin features is estimated by comparing predicted probabilities of two sequences where the bases at the central position are the reference and alternative alleles of a given variant. We reran the four models now restricting the cis-SNPs to those 36,241 predicted DeepSEA variants that were available in our imputed, quality-controlled genotype dataset. Top results for models including G, GxE and G +E are depicted in Supplementary Data 9–12. Results were comparable to what we observed before: 1195 (30.09%) of tagCpGs presented with best model G, 1193 CpGs (30.04%) with best model G +E, 1510 CpGs (38.02%) with best Determine variably methylated regions (VMRs): CpG-sites with MAD-score > ninetieth percentile and at least 2 consecutive CpGs with at most 1 kb distance Model E: tagCpG ~ environmental phenotypes Keep model with lowest AIC across all E models tagCpG: choose CpG-site with highest MADscore within each VMR as representative Model G: tagCpG ~ cis DeepSEA variants Keep model with lowest AIC across all G models Model G+E: tagCpG ~ cis DeepSEA variants + environmental phenotypes Keep model with lowest AIC across all G+E models Model GxE: tagCpG ~ cis DeepSEA variants x environmental phenotypes Keep model with lowest AIC across all GxE models Determine model with lowest AIC across E, G, G+E and GxE models as best model for each tagCpG Functional annotation of tagCpGs/DeepSEA variants stratified by best model E, G, G+E, GxE Replication of partition in best model E, G, G+E and GxE in independent cohorts For each tagCpG For all DeepSEA SNPs in 1 MB cis distance to tagCpGs For ten prenatal E For ten prenatal E x DeepSEA SNPs in 1 MB cis of tag CpG Fig. 1 Flow diagram of VMR analysis NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 ARTICLE NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 3 model GxE and 74 CpGs (1.86%) with best model E (Fig. 3a) and also showed similar differences in delta-AIC and proportions of E categories (see Supplementary Note 3). Only 10 tagCpGs did not present with any DeepSEA variant within 1MB distance in cis and were therefore not further considered. All respective CpGenvironment-DeepSea SNP combinations are depicted in Supplementary Data 13–16. The distribution of best models was not influenced by the degree of variability of DNA methylation, but was comparable across the whole range of DNA methylation variation (see Supplementary Note 4 and Supplementary Fig. 2). A slight enrichment for G +E models was observed in longer VMRs with at least 3 CpGs (p=9.00 × 10−06,OR=1.39, Fisher-test, see Supplementary Fig. 3). In conclusion, also when we focus on potentially functionally relevant SNPs, it is the combination of genotype and environment which best explains VMRs. We observed that, as expected, different types of exposures or maternal factors have different relative impact on DNA methylation (see Supplementary Note 5). However, even for those exposures with the highest fraction of VMRs best explained by E alone, combined models of G +E and GxE remain the best models in even higher fractions of VMRs (see Supplementary Fig. 4B). Functional annotation of different best models. Focusing on combinations between tagCpGs, environmental factors and DeepSEA variants, we found functional differences for both the SNPs as well as the tagCpGs (see Supplementary Note 6) within the different models. Overall, 895 DeepSEA variants were uniquely involved in best G models, 905 were uniquely in best G+E models and 1162 uniquely in best GxE models. As a DeepSEA variant can be in multiple 1 MB-cis windows around the tagCpGs, several DeepSEA variants were involved in multiple best models: 138 DeepSEA variants overlapped between G and GxE, 118 between G and G +E and 147 between GxE and G +E VMRs. We observed no significant differences with regard to gene-centric location for DeepSEA variants involved only in G models, only in G +E models or in multiple models. However, DeepSEA variants involved only in GxE models were significantly depleted for promoter locations (p=3.92 × 10−02,OR=0.79, Fisher-test, see Supplementary Fig. 5A). Although no significant differences were present, DeepSEA SNPs involved in the G and G +E model were located in closer proximity to the specific CpG (model G: mean absolute distance =256.8 kb, sd =291.2 kb, model G +E: mean absolute distance =244.8 kb, sd =284.0 kb, Supplementary Fig. 5B) whereas DeepSEA SNPs involved in GxE models (mean absolute distance =352.6 kb, sd =305.3 kb) showed broader peaks around the CpGs. With regards to histone marks, DeepSEA variants in general were enriched across multiple histone marks indicative of active transcriptional regulation (Fig. 4c). DeepSEA variants involved in best model G +E showed further enrichment for strong transcription (p=7.19 × 10−03,OR=1.34, Fisher-test) as well as depletion for quiescent loci (p=7.17 × 10−03,OR=0.78, Fisher-test). In contrast, GxE DeepSEA variants were significantly enriched in these regions (p=2.62 × 10−02,OR=1.22, Fishertest, Fig. 4d). Taken together, these analyses indicate that both the genetic variants and the VMRs in the different best models (G, GxE and G+E) preferentially annotate to functionally distinct genomics regions. Table 1 Overview of investigated cohorts Cohort PREDO I PREDO II DCHS I DCHS II UCI MoBa Sample size 817 146 107 151 121 1023 Methylation array Illumina 450 K Illumina EPIC Illumina 450 K Illumina EPIC Illumina EPIC Illumina 450 K Methylation data processing Funnorm and Combat Funnorm and Combat SWAN and Combat BMIQ and Combat Funnorm and Combat BMIQ and Combat SNP genotyping Illumina Human Omni Express Exome Illumina Human Omni Express Exome Illumina PsychArray Illumina GSA Illumina Human Omni Express Illumina HumanExome Core Infant gender male 433 (53.0%) 75 (51.4%) 63 (58.8%) 83 (55.0%) 65 (53.7%) 478 (46.7%) Maternal age mean (sd) 33.28 (5.79) 32.25 (4.92) 26.27 (5.87) 27.42 (5.93) 28.47 (4.91) 29.92 (4.32) Partity mean (sd) 1.05 (1.02) 0.87 (1.03) 0.98 (1.12) 1.09 (1.07) 1.11 (1.15) 0.83 (0.88) Caesarian section 169 (20.7%) 36 (24.7%) 19 (17.6%) 35 (23.2%) 37 (30.6%) 228 (22.3%) Pre-pregnancy BMI mean (sd) 27.42 (6.40) 25.37 (5.79) Not available Not available 27.90 (6.44) 24.05 (4.19) Maternal smoking yes Exclusion criterion Exclusion criterion 7.40 (10.52)a4.94 (9.43)a10 (8.2%) 148 (14.4%) Gestational diabetes yes 183 (22.4%) 19 (13.0%) No cases available No cases available 9 (7.4%) 15 (1.5%) Hypertension yes 275 (33.7%) 31 (21.2%) 2 (0.19%) 2 (1.3%) 7 (5.8%) 50 (4.9%) Betamethasone treatment yes 35 (4.3%) 2 (1.5%) Not available Not available No cases available Not available Anxiety score mean (sd) 33.93 (7.90)b34.43 (8.38)b5.70 (4.15)c5.32 (3.91)c1.67 (0.41)d4.79 (1.36)e Depression score mean (sd) 11.34 (6.47)f11.53 (6.98)f17.64 (12.10)g12.52 (11.55)g0.68 (0.41)h5.24 (1.57)e aBased on ASSIST Tobacco Score bSTAI sum scores cSRQ-20 dSTAI average scores eBased on Hopkins Symptom Checklist fCESD sum scores gBDI-II hCESD average score ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 4NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications Replication of best models in independent cohorts. To assess whether the relative distribution of the best models for VMRs and DeepSEA variants was stable across different samples, we assessed the relative distribution of these models in 3 additional samples (DCHS I and DCHS II, UCI and PREDOII) with VMR data both from the Illumina 450 K as well as the IlluminaHumanEPIC arrays. Information on these cohorts is summarised in Table 1 and the number of VMRs, the distribution of VMR methylation levels, VMR length and specific SNP information are given in Supplementary Note 7 and Supplementary Fig. 6. While major maternal factors overlapped among the cohorts - such as maternal age, delivery method, parity and depression during pregnancy - there were also differences, as the nonPREDO cohorts did not include betamethasone treatment but additionally included maternal smoking (see Table 1). Despites these differences and differences in the total number of VMRs, the overall pattern remained stable: in all 4 analyses, DCHS I, DCHS II, UCI and PREDO II, we replicated that E alone models almost never explained most of the variances, while G alone models explained the most variance in up to 15% of the VMRs; G+E in up to 32%; and GxE models in up to 60% (see Fig. 5and Table 2). The importance of including G for a best model fit could also be observed for maternal smoking, described as one of the most highly replicated factors shaping the newborns’methylome11 and present in the replication but not the discovery cohort PREDO I. These analyses are detailed in Supplementary Note 8. We were also able to replicate our finding showing that GxE VMRs were enriched for OpenSea positions with a trend on the 450 K array (DCHS I, OR =1.11, p=5.03 × 10−02, Fisher-test) and significantly for the EPIC array data (PREDOII: p=2.96 × 10−06,OR=1.29, UCI: p=3.79 × 10−02,OR=1.09, DCHSII: p=2.91 × 10−04,OR=1.16, Fisher-tests). For all additional cohorts, the delta AIC for best model GxE to the next best E G G+E GxE 0 5 10 15 ab c p = 4.78×10–80 p = 2.22×10–96 40.58% 29.41% 29.98% 0.00 0.25 0.50 0.75 1.00 PREDOI pruned genotypes Percentage of CpGs Type E G G+E GxE 7.47% 11.82% 12.66% 24.64% 9.30% 9.05% 24.81% 8.54% 10.59% 12.65% 28.41% 9.67% 8.47% 21.59% 6.62% 7.96% 10.67% 30.64% 6.29% 6.90% 30.56% 0.00 G+E GxE tagCpGs 0.25 0.50 0.75 1.00 Percentage of CpGs Type Anxiety score Betamethasone intake Delivery mode Depression score Gestational diabetes Maternal age Maternal hypertension Maternal pre–pregnancy BMI Ogtt Parity Fig. 2 VMR analysis in pruned PREDO I dataset. aPercentage of models (G, E, GxE or G +E) with the lowest AIC explaining variable DNA methylation using the PREDO I dataset with pruned SNPs. bDistribution of the different types of prenatal environment included in the E model with the lowest AIC (right), in the combinations yielding the best model GxE (middle), or the best model G +E models (left). To increase readability all counts <3% have been omitted. cDeltaAIC, i.e, the difference in AIC, between best model and next best model, stratified by the best model. Y-axis denotes the delta AIC and the X-axis the different models. The median is depicted by a black line, the rectangle spans the first quartile to the third quartile, whiskers above and below the box show the location of minimum and maximum beta-values. P-values are based on Wilcoxon-tests NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 ARTICLE NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 5 model was also significantly higher as compared to CpGs with G, EorG+E as the best model. Overall, 387 tag CpGs overlapped between PREDO I, PREDO II, DCHS I and DCHS II (see Supplementary Fig. 7), which allowed us to test the consistency of the best models for specific VMRs across the different cohorts. Over 70% of the overlapping tagCPGs showed consistent best models in at least 3 cohorts (see Fig. 6) with GxE being the most consistent model (for over 60% of consistent models, see Supplementary Fig. 8). Focusing only on EPIC data (PREDO II, DCHSII and UCI), we identified more, namely 2091, tag CpGs that overlap across the three cohorts and here 86% show a consistent best model in at least two of the three cohorts, despite differences in study design, prenatal phenotypes and ethnicity. Thus, the additional cohorts not only showed a consistent replication of the proportion of the models best explaining variance of VMRs but also consistency of the best model for specific VMRs. Within this context, we observed the GxE models were the most consistent models across the cohorts (see Supplementary Fig. 8), with 85% of the CpGs with consistent models across 5 cohorts having GxE as the best model. Furthermore, we could validate specific GxE combinations between PREDO I and MoBa as shown as in the Supplementary Note 9, in Supplementary Data 17 and 18 and in Supplementary Fig. 9. Disease relevance. Finally, we tested whether functional DeepSEA SNPs involved in only G, only GxE and only G +E models in PREDO I for their enrichment in GWAS hits. We used all functional SNPs and their LD proxies (defined as r2of at least 0.8 in the PREDO cohort and in maximal distance of 1MB to the target SNP) and performed enrichment analysis with the overlap of nominal significant GWAS hits. We selected for a broad spectrum of GWAS, including GWAS for complex disorders for which differences in prenatal environment are established as risk factors, but also including GWAS on other complex diseases. For psychiatric disorders, we used summary statistics of the Psychiatric Genomics Consortium (PGC) including association studies for autism40, attention-deficit-hyperactivity disorder41, ab c 38.02% 30.04% 30.09% 1.86% 0.00 0.25 0.50 0.75 1.00 Percentage of CpGs Type E G G+E GxE 36.49% 6.76% 44.59% 12.16% 32.97% 5.61% 40.08% 21.34% 31.6% 5.28% 37.55% 25.57% 30.13% 5.50% 41.65% 22.72% 31.55% 5.10% 37.55% 25.80% 0.00 0.25 0.50 0.75 1.00 E G G+E GxE VMRs Percentage of CpGs Type Island OpenSea Shelf Shore 35.13% 18.92% 31.60% 27.03% 8.11% 31.72% 21.34% 5.02% 28.2% 3.01% 5.61% 27% 22.04% 5.28% 30.93% 3.52% 6.37% 30.27% 18.34% 6.42% 29.47% 3.78% 7.28% 31.51% 19.52% 5.45% 29.08% 3.28% 6.44% 0 25 50 75 100 E G G+E GxE VMRs Percentage of CpGs Type 1st Exon 1st Intron 3′ UTR 5′ UTR Distal intergenic Downstream (<=3 kb) Other exon Other intron Promoter PREDO I DeepSEA SNPs Fig. 3 VMR analysis in DeepSEA annotated SNPs in PREDO I dataset. aPercentage of models (G, E, GxE or G +E) with the lowest AIC explaining variable DNA methylation using the PREDO I dataset with DeepSEA annotated SNPs. bDistribution of the locations of all VMRs and tagVMRs with best model E, G, G+E and GxE on the 450k array using only DeepSEA variants in relationship to CpG-Islands based on the Illumina 450 K annotation. cDistribution of gene-centric locations of all VMRs and tagVMRs with best model E, G, G +E and GxE on the 450k array using only DeepSEA variants ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 6NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications bipolar disorder42, major depressive disorder43, schizophrenia44 and the cross-disorder associations including all five of these disorders45. Additionally, we included GWAS of inflammatory bowel disease46, type 2 diabetes47 and for BMI48. Nominal significant GWAS findings were enriched for DeepSEA variants and their LD proxies per se across psychiatric as well as nonpsychiatric diseases (Fig. 7a). However, G, GxE and G +E DeepSEA variants showed a differential enrichment pattern above all DeepSEA variants (Fig. 7b), with strongest enrichments of GxE DeepSEA variants in GWAS of autism spectrum disorder (p< 2.20 × 10−16,OR=2.07 above DeepSEA, Fisher-test), attention-deficit-hyperactivity disorder (p< 2.20 × 10−16,OR= 1.71, Fisher-test) and inflammatory bowel disease (p< 2.20 × 10−16,OR=1.71, Fisher-test) and G +E DeepSEA variants in GWAS for attention-deficit-hyperactivity disorder (p=9.54 × 10−36,OR=1.23, Fisher-test) and inflammatory bowel disease (p=1.85 × 10−52,OR=1.30, Fisher-test). While SNPs with strong main meQTL effects such as those within G and G +E models have been reported to be enriched in GWAS for common disease, we now also show this for SNPs within GxE models that often have non-significant main G effects. Discussion We evaluated the effects of prenatal environmental factors and genotype on DNA methylation at VMRs identified in neonatal blood samples. We found that most variable methylation sites were best explained by either genotype and prenatal environment a b c 0 1 2 3 Odds ratio E G G+E GxE Active TSS Bivalent enhancer Bivalent/poised TSS Enhancers Flanking active TSS Flanking bivalent TSS/Enh Genic enhancers Heterochromatin Quiescent/low Repressed PolyComb Strong transcription Transcr. at gene 5′ and 3′ Weak repressed PolyComb Weak transcription ZNF genes & repeats 0.0 2.5 5.0 7.5 10.0 12.5 Odds ratio G G+E GxE Active TSS Bivalent enhancer Bivalent/poised TSS Enhancers Flanking active TSS Flanking bivalent TSS/Enh Genic enhancers Heterochromatin Quiescent/low Repressed PolyComb Strong transcription Transcr. at gene 5′ and 3′ Weak repressed PolyComb Weak transcription ZNF genes & repeats 0.8 0.9 1.0 1.1 1.2 1.3 OR d Significant Non significant Fig. 4 Functional annotation of VMR-mapping in DeepSEA annotated SNPs in PREDO I dataset. aHistone mark enrichment for all VMRs. The Y-axis denotes the fold enrichment/depletion as compared to no-VMRs. Blue bars indicate significant enrichment/depletion, grey bars non-significant differences based on Fisher-tests. bHistone mark enrichment for tagVMRs with best model E, G, G +E and GxE relative to all VMRs. Green colour indicates depletion, red colour indicates enrichment. Thick black lines around the rectangles indicate significant enrichment/depletion based on Fisher-tests. cHistone mark enrichment for all DeepSEA variants in the dataset. Blue bars indicate significant enrichment/depletion based on Fisher-tests. dHistone mark enrichment for all DeepSEA variants involved in models where either G, G +E or GxE is the best model as compared to all tested DeepSEA variants. Green colour indicates depletion, red colour indicates enrichment. Thick black lines around the rectangles indicate significant enrichment/depletion based on Fishertests NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 ARTICLE NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 7 interactions (GxE) or additive effects (G +E) of these factors, followed by main genotype effects. This pattern was replicated in independent cohorts and underscores the need to consider genotype in the study of environmental effects on DNA methylation. In fact, VMRs best explained by G, G +E or GxE and their associated functional genetic variants were located in distinct genomic regions, suggesting that different combinatorial effects of G and E may impact VMRs with distinct downstream regulatory effects and thus possibly context-dependent impact on cellular function. We also observed that functional variants with best models G, G +E or GxE, all showed significant enrichment within GWAS signals for complex disorders beyond the enrichment of the functional variants themselves. While this was expected for G and G +E models based on results from previous studies21,23,24,26, it was surprising for GxE SNPs, as these often do not have highly significant main genetic effects. Their specific enrichment in GWAS for common disorders supports the importance of these genetic variants that moderate environmental impact both at the level of DNA methylation but also, potentially, for disease risk. The fact that GxE and G +E best explained the majority of VMRs (see Fig. 5) and that GxE models were selected by a larger margin than the other models (see Fig. 2c) was consistently found across all tested cohorts. These findings are in line with a previous report by Teh et al.29 who performed a similar analysis based on AIC in umbilical cord tissue. Differences to the findings by Teh et al. are discussed in the Supplemental Discussion. Using data from four different cohorts, we not only saw comparable proportions of VMRs best explained by the different models, but also saw in the VMRs common across cohorts that specific VMRs had consistent best models (see Fig. 6). This is in line with the fact that VMRs best explained by G, GxE or G +E show functional differences and may differentially impact gene regulation. In addition to consistent findings using AIC-based approaches, we also observed some indication for validation of individual GxE and G +E combinations on selected VMRs using p-value based criteria, with a small number of specificG+E and GxE effects on VMRs replicating between the PREDO I and the MoBa cohort. The low number of specific replications could be due to lack of overall power as well as larger differences in prenatal factors between these two cohorts (see Table 1). As shown in Supplementary Fig. 4B, which specific G and E combinations best explain VMRs is also dependent on the specific prenatal factors. Larger and more homogenous cohorts regarding exposures will be needed for such analyses to be more conclusive. While E alone was rarely the best model, it should be pointed out that main environmental effects on DNA methylation were observed (see Supplementary Data 3), and consistent with previous large meta-analyses such as in the case of maternal smoking (see Supplementary Note 7). Within the MoBa cohort, the cohort with the largest proportion of maternal smoking, 10% of all tagCpGs were best explained by maternal smoking alone. However, in all other cohorts, where smoking was less prevalent, the inclusion of genotypic effects in addition to maternal smoking explained more of the variance. This supports that while main E effects on the newborn methylome are present, genotype is an important factor that, in combination with E, may explain even more of the variance in DNA methylation. VMRs best explained by either E, G, G +E or GxE and their associated functional SNPs were enriched for distinct genomics locations and chromatin states (see Fig. 4), suggesting that VMRs moderated by different combinations of G and E may in fact have 38.02% 30.03% 30.09% 53.97% 30.78% 15.12% 55.21% 28.94% 15.30% 60.08% 23.69% 12.22% 4.01% 56.57% 32.20% 11.11% 0.00 PREDOI_450K DCHSI_450K PREDOII_EPIC UCI_EPIC DCHSII_EPIC 0.25 0.50 0.75 1.00 Type E G G+E GxE Fig. 5 VMR analysis in PREDO I and replication datasets. Percentage of models (G, E, GxE or G +E) with the lowest AIC explaining variable DNA methylation in PREDO I (450 K), DCHS I (450 K), PREDO II (EPIC), UCI (EPIC) and DCHS II (EPIC) Table 2 VMRs and best models across cohorts Cohort PREDO I PREDO II DCHS I DCHS II UCI Sample-size 817 146 107 151 121 Methylation array Illumina 450 K Illumina EPIC Illumina 450 K Illumina EPIC Illumina EPIC # VMRs 3972 8547 6072 10,005 9525 Proportion: best model E (%) 2.0 <1 <1 <1 4.1 Best model G (%) 30.0 15.0 15.8 11.5 12.8 Best model G +E (%) 30.0 29.0 29.8 32.1 24.1 Best model GxE (%) 38.0 56.0 54.3 56.3 59.0 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 8NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications distinct functional roles in gene regulation. Overall, VMRs best explained by GxE were consistently enriched for regions annotated to the OpenSea regions with lower CpG density and located farthest from CpG Islands49. Open Sea regions have been reported to be enriched for environmentally-associated CpGs with for example exposure to childhood trauma50 and may harbour more long-range enhancers. In addition to their position relative to CpG islands and their CpG content, G, GxE and G +E VMRs and their associated functional SNPs also showed distinct enrichments for chromatin marks. Compared to 450 K VMRs in general, VMRs with GxE as the best models were relatively depleted in regions surrounding the TSS, while VMRs with G +E were relatively enriched in these regions (see Fig. 4), suggesting that GxE VMRs are located at more distance from the TSS than G +E VMRs. To better map the potential functional variants in these models and to compare methylation-associated SNPs from a regulatory perspective, we used DeepSEA38, a machine learning algorithm that predicts SNP functionality from the sequence context based on sequencing data for different regulatory elements in different cell lines using ENCODE data39. We identified the SNPs with putatively functional consequences on regulatory marks by DeepSEA and compared putative regulatory effects of G, G +E and GxE hits. Relative to the imputed non-DeepSEA SNPs contained in our dataset, these predicted functional DeepSEA SNPs were enriched for TSS and enhancer regions and depleted for quiescent regions, supporting their relevance in regulatory processes (see Fig. 4). Compared to DeepSEA SNPs overall, DeepSEA SNPs within the three different best models also showed distinct enrichment or depletion patterns. Similar to GxE VMRs, likely functional GxE SNPs also showed a relative depletion in TSS regions while G +E SNPs showed enrichment in genic enhancers. Overall, both the VMRs as well as the associated functional SNPs appear to be in distinct regulatory regions, depending on their best model. In addition, GxE functional SNP and tagCpGs were located farther apart than SNP/tagCpG pairs within G or G +E models (see Supplementary Fig. 5B), supporting a more long-range type of regulation in GxE interactions on molecular traits as compared to all genes; a similar relationship has been reported previously for GxE with regard to gene expression in C. elegans51,52. SNPs associated with differences in gene expression but also DNA methylation have consistently been shown to be enriched among SNPs associated with common disorders in GWAS21,24,26,53. The functional genetic variants that were within G, GxE or G +E models predicting variable DNA methylation were even enriched in GWAS association results (beyond the baseline enrichment of DeepSea SNPs per se). The fact that such enrichment was observed for not only G and G +E SNPs, with strong main genetic effects, but also for GxE SNPs, with smaller to sometimes no main genetic effect on DNA methylation underscores the importance of also including SNPs within GxE models in the functional annotation of GWAS. A detailed catalogue of meQTLs that are responsive to environmental factors could support a better pathophysiological understanding of diseases for which risk is shaped by a combination of environment and genetic factors. Finally, we want to note the limitations of this study. First, we restricted our analyses to specific DNA methylation array contents that are inherently biased as compared to genome-wide bisulfite sequencing, for example. In addition, we restricted our analysis to VMRs, which also limits the generalisability of the a b 0.0 0.5 1.0 1.5 Odds ratio Fishertest Significant G G+E GxE ADHD ASD BMI BP CrossDisorder IBD MDD SCZ T2D 1.00 1.25 1.50 1.75 2.00 OR Fig. 7 Enrichment of DeepSEA variants for GWAS associations. aEnrichment for nominal significant GWAS associations for all tested DeepSEA variants and their LD proxies for GWAS for ADHD (attentiondeficit hyperactivity disorder), ASD (autism spectrum disorder), BMI (body mass index), BP (bipolar disorder), CrossDisorder, IBD (inflammatory bowel disease), MDD (major depressive disorder), SCZ (schizophrenia) and T2D (Type 2 diabetes). The Y-axis denotes the fold enrichment with regard to non-DeepSEAvariants. Blue bars indicate significant enrichment/ depletion based on Fisher-tests. bEnrichment for nominal significant GWAS hits for DeepSEA variants and their LD proxies involved in best models with G, G +E or GxE as compared to all tested DeepSEA variants. Green colour indicates depletion, red colour indicates enrichment. Thick black lines around the rectangles indicate significant enrichment/depletion based on Fisher-tests 3.62% 21.96% 49.1% 25.58% 0.00 0.25 0.50 0.75 1.00 Overlapping tagCpGs Percentage of CpGs Type 2 or less consistent cohorts 3 consistent cohorts 4 consistent cohorts 5 consistent cohorts Fig. 6 Consistency of best models across cohorts. Percentage of consistent best models in overlapping tag CpGs of PREDO I (450 K), DCHS I (450 K), PREDO II (EPIC), UCI (EPIC) and DCHS II (EPIC). Overlapping VMRs included significantly more CpGs as compared to all VMRs (p< 2.2 × 10−16, Wilcoxon-test, mean =4.43) NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 ARTICLE NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 9 Health, University of Cape Town, Cape Town 7925, South Africa. 27 South African Medical Research Council (SAMRC), Unit on Risk and Resilience in Mental Disorders, Cape Town 7505, South Africa. 28 Department of Paediatrics & Child Health and SAMRC Unit on Child and Adolescent Health, University of Cape Town, Cape Town 7505, South Africa. 29 Department of Epidemiology, Harvard T. H. Chan School of Public Health, Boston, MA 02115, USA. 30 Department of Mathematics, Technische Universität München, Munich 85748, Germany. 31 Department of Psychiatry and Behavioral Sciences, Emory University School of Medicine, Atlanta 30329, USA. † A full list of consortium members appears at the end of the paper. Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium Naomi R. Wray32,33, Stephan Ripke34,35,36, Manuel Mattheisen37,38,39,40, Maciej Trzaskowski32, Enda M. Byrne32, Abdel Abdellaoui41, Mark J. Adams42, Esben Agerbo40,43,44, Tracy M. Air45, Till F.M. Andlauer1,46, Silviu-Alin Bacanu47, Marie Bækvad-Hansen40,48, Aartjan T.F. Beekman49, Tim B. Bigdeli47,50, Douglas H.R. Blackwood42, Julien Bryois51, Henriette N. Buttenschøn39,40,52, Jonas Bybjerg-Grauholm40,48, Na Cai53,54, Enrique Castelao55, Jane Hvarregaard Christensen38,39,40, Toni-Kim Clarke42, Jonathan R.I. Coleman56, Lucía Colodro-Conde57, Baptiste Couvy-Duchesne58,59, Nick Craddock60, Gregory E. Crawford61,62, Gail Davies63, Ian J. Deary63, Franziska Degenhardt64,65, Eske M. Derks57, Nese Direk66,67, Conor V. Dolan41, Erin C. Dunn68,69,70, Thalia C. Eley56, Valentina Escott-Price71, Farnush Farhadi Hassan Kiadeh72, Hilary K. Finucane58,73, Andreas J. Forstner64,65,74,75, Josef Frank76, Héléna A. Gaspar56, Michael Gill77, Fernando S. Goes78, Scott D. Gordon79, Jakob Grove38,39,40,80, Lynsey S. Hall42,81, Christine Søholm Hansen40,48, Thomas F. Hansen82,83,84, Stefan Herms64,65,75, Ian B. Hickie85, Per Hoffmann47,64,65, Georg Homuth86, Carsten Horn87, Jouke-Jan Hottenga41, David M. Hougaard40,48, Marcus Ising88, Rick Jansen49, Eric Jorgenson89, James A. Knowles90, Isaac S. Kohane91,92,93, Julia Kraft35, Warren W. Kretzschmar94, Jesper Krogh95, Zoltán Kutalik96,97, Yihan Li94, Penelope A. Lind57, Donald J. MacIntyre98,99, Dean F. MacKinnon78, Robert M. Maier33, Wolfgang Maier100, Jonathan Marchini101, Hamdi Mbarek41, Patrick McGrath102, Peter McGuffin56, Sarah E. Medland57, Divya Mehta33,103, Christel M. Middeldorp41,104,105, Evelin Mihailov106, Yuri Milaneschi49, Lili Milani106, Francis M. Mondimore78, Grant W. Montgomery33, Sara Mostafavi107,108, Niamh Mullins56, Matthias Nauck109,110, Bernard Ng108, Michel G. Nivard41, Dale R. Nyholt111, Paul F. O’Reilly56, Hogni Oskarsson112, Michael J. Owen113, Jodie N. Painter57, Carsten Bøcker Pedersen40,43,44, Marianne Giørtz Pedersen40,43,44, Roseann E. Peterson47,114, Erik Pettersson51, Wouter J. Peyrot49, Giorgio Pistis55, Danielle Posthuma115,116, Jorge A. Quiroz117, Per Qvist38,39,40, John P. Rice118, Brien P. Riley47, Margarita Rivera56,119, Saira Saeed Mirza66, Robert Schoevers120, Eva C. Schulte121,122, Ling Shen89, Jianxin Shi123, Stanley I. Shyn124, Engilbert Sigurdsson125, Grant C.B. Sinnamon126, Johannes H. Smit49, Daniel J. Smith127, Hreinn Stefansson128, Stacy Steinberg128, Fabian Streit76, Jana Strohmaier76, Katherine E. Tansey129, Henning Teismann130, Alexander Teumer131, Wesley Thompson40,83,132,133, Pippa A. Thomson132, Thorgeir E. Thorgeirsson128, Matthew Traylor134, Jens Treutlein76, Vassily Trubetskoy35, André G. Uitterlinden135, Daniel Umbricht136, Sandra Van der Auwera137, Albert M. van Hemert138, Alexander Viktorin51, Peter M. Visscher32,33, Yunpeng Wang40,83,133, Bradley T. Webb139, Shantel Marie Weinsheimer40,83, Jürgen Wellmann130, Gonneke Willemsen41, Stephanie H. Witt76, Yang Wu32, Hualin S. Xi140, Jian Yang33,141, Futao Zhang32, Volker Arolt142, Bernhard T. Baune45, Klaus Berger130, Dorret I. Boomsma41, Sven Cichon64,75,143,144, Udo Dannlowski142, E.J.C. de Geus10,145, J. Raymond DePaulo78, Enrico Domenici146, Katharina Domschke147, Tõnu Esko36,106, Hans J. Grabe137, Steven P. Hamilton148, Caroline Hayward149, Andrew C. Heath118, Kenneth S. Kendler47, Stefan Kloiber88,150,151, Glyn Lewis152, Qingqin S. Li153, Susanne Lucae88, Pamela A.F. Madden118, Patrik K. Magnusson51, Nicholas G. Martin79, Andrew M. McIntosh42,63, Andres Metspalu106,154, Ole Mors40,155, Preben Bo Mortensen39,40,43,44, Bertram Müller-Myhsok1,46,156, Merete Nordentoft40,157, Markus M. Nöthen64,65, Michael C. O’Donovan113, Sara A. Paciga158, Nancy L. Pedersen51, ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 16 NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications Brenda W.J.H. Penninx49, Roy H. Perlis68,159, David J. Porteous160, James B. Potash161, Martin Preisig55, Marcella Rietschel76, Catherine Schaefer89, Thomas G. Schulze76,122,162,163,164, Jordan W. Smoller68,69,70, Kari Stefansson155,165, Henning Tiemeier66,166,167, Rudolf Uher168, Henry Völzke131, Myrna M. Weissman102,169, Thomas Werge40,83,170, Cathryn M. Lewis56,171, Douglas F. Levinson172, Gerome Breen56,173, Anders D. Børglum38,39,40 & Patrick F. Sullivan51,174,175 32 Institute for Molecular Bioscience, The University of Queensland, Brisbane 4072 QLD, Australia. 33 Queensland Brain Institute, The University of Queensland, Brisbane 4072 QLD, Australia. 34 Analytic and Translational Genetics Unit, Massachusetts General Hospital, Boston, MA 02114, USA. 35 Department of Psychiatry and Psychotherapy, Universitätsmedizin Berlin Campus Charité Mitte, Berlin 14129, Germany. 36 Medical and Population Genetics, Broad Institute, Cambridge, MA 02142, USA. 37 Centre for Psychiatry Research, Department of Clinical Neuroscience, Karolinska Institutet, Stockholm 17177 SE, Sweden. 38 Department of Biomedicine, Aarhus University, Aarhus 8000, Denmark. 39 iSEQ, Centre for Integrative Sequencing, Aarhus University, Aarhus 8000, Denmark. 40 iPSYCH, The Lundbeck Foundation Initiative for Integrative Psychiatric Research, Aarhus 8000, Denmark. 41 Department of Biological Psychology & EMGO+Institute for Health and Care Research, Vrije Universiteit Amsterdam, Amsterdam 1081 BT, Netherlands. 42 Division of Psychiatry, University of Edinburgh, Edinburgh EH10 5HF, UK. 43 Centre for Integrated Registerbased Research, Aarhus University, Aarhus 8210, Denmark. 44 National Centre for Register-Based Research, Aarhus University, Aarhus 8210, Denmark. 45 Discipline of Psychiatry, University of Adelaide, Adelaide 5000 SA, Australia. 46 Munich Cluster for Systems Neurology (SyNergy), Munich 81377, Germany. 47 Department of Psychiatry, Virginia Commonwealth University, Richmond, VA 22903, USA. 48 Center for Neonatal Screening, Department for Congenital Disorders, Statens Serum Institut, Copenhagen 2300, Denmark. 49 Department of Psychiatry, Vrije Universiteit Medical Center and GGZ inGeest, Amsterdam 1081 NL, Netherlands. 50 Virginia Institute for Psychiatric and Behavior Genetics, Richmond, VA 23298, USA. 51 Department of Medical Epidemiology and Biostatistics, Karolinska Institutet, Stockholm 17177 SE, Sweden. 52 Department of Clinical Medicine, Translational Neuropsychiatry Unit, Aarhus University, Aarhus 8240, Denmark. 53 Human Genetics, Wellcome Trust Sanger Institute, Cambridge, CB10 1SA, UK. 54 Statistical genomics and systems genetics, European Bioinformatics Institute (EMBL-EBI), Cambridge, CB10 1 SD, UK. 55 Department of Psychiatry, University Hospital of Lausanne, Prilly, Vaud 1004, Switzerland. 56 MRC Social Genetic and Developmental Psychiatry Centre, King’s College London, London WC2R 2LS, UK. 57 Genetics and Computational Biology, QIMR Berghofer Medical Research Institute, Herston 4006 QLD, Australia. 58 Centre for Advanced Imaging, The University of Queensland, Saint Lucia 4072 QLD, Australia. 59 Queensland Brain Institute, The University of Queensland, Saint Lucia 4072 QLD, Australia. 60 Psychological Medicine, Cardiff University, Cardiff CF14 4XN, UK. 61 Center for Genomic and Computational Biology, Duke University, Durham, NC 27705, USA. 62 Division of Medical Genetics, Department of Pediatrics, Duke University, Durham, NC 27708, USA. 63 Centre for Cognitive Ageing and Cognitive Epidemiology, University of Edinburgh, Edinburgh EH8 9JZ, UK. 64 Institute of Human Genetics, University of Bonn, Bonn 53127 DE, Germany. 65 Life & Brain Center, Department of Genomics, University of Bonn, Bonn 53127, Germany. 66 Epidemiology, Erasmus MC, Rotterdam 3015 Zuid-Holland, Netherlands. 67 Psychiatry, Dokuz Eylul University School Of Medicine, Izmir 35220, Turkey. 68 Department of Psychiatry, Massachusetts General Hospital, Boston, MA 02114, USA. 69 Psychiatric and Neurodevelopmental Genetics Unit (PNGU), Massachusetts General Hospital, Boston, MA 02114, USA. 70 Stanley Center for Psychiatric Research, Broad Institute, Cambridge, MA 02142, USA. 71 Neuroscience and Mental Health, Cardiff University, Cardiff CF24 4HQ, UK. 72 Bioinformatics, University of British Columbia, Vancouver V5Z 4S6 BC, Canada. 73 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02142, USA. 74 Department of Psychiatry (UPK), University of Basel, Basel 4002, Switzerland. 75 Human Genomics Research Group, Department of Biomedicine, University of Basel, Basel 4031, Switzerland. 76 Department of Genetic Epidemiology in Psychiatry, Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University, Mannheim 68159 Baden-Württemberg, Germany. 77 Department of Psychiatry, Trinity College Dublin, Dublin 8, Ireland. 78 Psychiatry & Behavioral Sciences, Johns Hopkins University, Baltimore, MD 21287, USA. 79 Genetics and Computational Biology, QIMR Berghofer Medical Research Institute, Brisbane 4006 QLD, Australia. 80 Bioinformatics Research Centre, Aarhus University, Aarhus 8000, Denmark. 81 Institute of Genetic Medicine, Newcastle University, Newcastle upon Tyne NE1 3BZ, England. 82 Danish Headache Centre, Department of Neurology, Rigshospitalet, Glostrup 2600, Denmark. 83 Institute of Biological Psychiatry, Mental Health Center Sct. Hans, Mental Health Services Capital Region of Denmark, Copenhagen 4000, Denmark. 84 iPSYCH, The Lundbeck Foundation Initiative for Psychiatric Research, Copenhagen 8000, Denmark. 85 Brain and Mind Centre, University of Sydney, Sydney 2050 NSW, Australia. 86 Interfaculty Institute for Genetics and Functional Genomics, Department of Functional Genomics, University Medicine and Ernst Moritz Arndt University Greifswald, Greifswald 17489 Mecklenburg-Vorpommern, Germany. 87 Roche Pharmaceutical Research and Early Development, Pharmaceutical Sciences, Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd, Basel 4070, Switzerland. 88 Max Planck Institute of Psychiatry, Munich 80804, Germany. 89 Division of Research, Kaiser Permanente Northern California, Oakland, CA 94612, USA. 90 Psychiatry & The Behavioral Sciences, University of Southern California, Los Angeles, CA 90033, USA. 91 Department of Biomedical Informatics, Harvard Medical School, Boston, MA 02115, USA. 92 Department of Medicine, Brigham and Women’s Hospital, Boston, MA 02115, USA. 93 Informatics Program, Boston Children’s Hospital, Boston, MA 02115, USA. 94 Wellcome Trust Centre for Human Genetics, University of Oxford, Oxford OX3 7BN, UK. 95 Department of Endocrinology at Herlev University Hospital, University of Copenhagen, Copenhagen 2730, Denmark. 96 Institute of Social and Preventive Medicine (IUMSP), University Hospital of Lausanne, Lausanne, VD 1010, Switzerland. 97 Swiss Institute of Bioinformatics, Lausanne, VD 1015, Switzerland. 98 Division of Psychiatry, Centre for Clinical Brain Sciences, University of Edinburgh, Edinburgh EH16 4SB, UK. 99 Mental Health, NHS 24, Glasgow G12 0XH, UK. 100 Department of Psychiatry and Psychotherapy, University of Bonn, Bonn 53105, Germany. 101 Statistics, University of Oxford, Oxford OX1 3LB, UK. 102 Psychiatry, Columbia University College of Physicians and Surgeons, New York, NY 10032, USA. 103 School of Psychology and Counseling, Queensland University of Technology, Brisbane, QLD 4059, Australia. 104 Child and Youth Mental Health Service, Children’s Health Queensland Hospital and Health Service, South Brisbane, QLD 4000, Australia. 105 Child Health Research Centre, University of Queensland, Brisbane, QLD 4101, Australia. 106 Estonian Genome Center, University of Tartu, Tartu 51005, Estonia. 107 Medical Genetics, University of British Columbia, Vancouver, BC V6H 3N1, Canada. 108 Statistics, University of British Columbia, Vancouver, BC V6T 1Z4, Canada. 109 DZHK (German Centre for Cardiovascular Research), Partner Site Greifswald, University Medicine, University Medicine Greifswald, Greifswald, Mecklenburg-Vorpommern 17489, Germany. 110 Institute of Clinical Chemistry and Laboratory Medicine, University Medicine Greifswald, Greifswald, Mecklenburg-Vorpommern 17489, Germany. 111 Institute of Health and Biomedical Innovation, Queensland University of Technology, Brisbane, QLD 4059, Australia. 112 Humus, Reykjavik 101, Iceland. 113 MRC Centre for Neuropsychiatric Genetics and Genomics, Cardiff University, Cardiff CF24 4HQ, UK. 114 Virginia Institute for Psychiatric & Behavioral Genetics, Virginia Commonwealth University, Richmond, VA 23298, USA. 115 Clinical Genetics, Vrije Universiteit Medical Center, Amsterdam 1081HV, Netherlands. 116 Complex Trait Genetics, Vrije Universiteit Amsterdam, Amsterdam 1081 HV, Netherlands. 117 Solid Biosciences, Boston, MA 02139, USA. 118 Department of Psychiatry, Washington University in Saint Louis NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 ARTICLE NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications 17 School of Medicine, Saint Louis, MO 63110, USA. 119 Department of Biochemistry and Molecular Biology II, Institute of Neurosciences, Center for Biomedical Research, University of Granada, Granada CP 18100, Spain. 120 Department of Psychiatry, University of Groningen, University Medical Center Groningen, Groningen 9700 RB, Netherlands. 121 Department of Psychiatry and Psychotherapy, Medical Center of the University of Munich, Campus Innenstadt, Munich 80336, Germany. 122 Institute of Psychiatric Phenomics and Genomics (IPPG), Medical Center of the University of Munich, Campus Innenstadt, Munich 80336, Germany. 123 Division of Cancer Epidemiology and Genetics, National Cancer Institute, Bethesda, MD 20892, USA. 124 Behavioral Health Services, Kaiser Permanente Washington, Seattle, WA 98112, USA. 125 Faculty of Medicine, Department of Psychiatry, University of Iceland, Reykjavik 101, Iceland. 126 School of Medicine and Dentistry, James Cook University, Townsville, QLD 4811, Australia. 127 Institute of Health and Wellbeing, University of Glasgow, Glasgow G12 8RZ, UK. 128 deCODE Genetics/Amgen, Reykjavik 101, Iceland. 129 College of Biomedical and Life Sciences, Cardiff University, Cardiff CF14 4EP, UK. 130 Institute of Epidemiology and Social Medicine, University of Münster, Münster, Nordrhein-Westfalen 48149, Germany. 131 Institute for Community Medicine, University Medicine Greifswald, Greifswald, Mecklenburg-Vorpommern 17489, Germany. 132 Department of Psychiatry, University of California, San Diego, San Diego, CA 92093, USA. 133 KG Jebsen Centre for Psychosis Research, Norway Division of Mental Health and Addiction, Oslo University Hospital, Oslo 0407, Norway. 134 Clinical Neurosciences, University of Cambridge, Cambridge CB2 1QW, UK. 135 Internal Medicine, Erasmus MC, Rotterdam, Zuid-Holland 3015, Netherlands. 136 Roche Pharmaceutical Research and Early Development, Neuroscience, Ophthalmology and Rare Diseases Discovery & Translational Medicine Area, Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd, Basel 4070, Switzerland. 137 Department of Psychiatry and Psychotherapy, University Medicine Greifswald, Greifswald, Mecklenburg-Vorpommern 17475, Germany. 138 Department of Psychiatry, Leiden University Medical Center, Leiden 2333 ZA, Netherlands. 139 Virginia Institute of Psychiatric & Behavioral Genetics, Virginia Commonwealth University, Richmond, VA 23298, USA. 140 Computational Sciences Center of Emphasis, Pfizer Global Research and Development, Cambridge, MA 02139, USA. 141 Institute for Molecular Bioscience; Queensland Brain Institute, The University of Queensland, Brisbane, QLD 4072, Australia. 142 Department of Psychiatry, University of Münster, Münster, Nordrhein-Westfalen 48149, Germany. 143 Institute of Medical Genetics and Pathology, University Hospital Basel, University of Basel, Basel 4031, Switzerland. 144 Institute of Neuroscience and Medicine (INM-1), Research Center Juelich, Juelich 52425, Germany. 145 Amsterdam Public Health Institute, Vrije Universiteit Medical Center, Amsterdam 1081 BT, Netherlands. 146 Centre for Integrative Biology, Università degli Studi di Trento, Trento, Trentino-Alto Adige 38123, Italy. 147 Department of Psychiatry and Psychotherapy, Medical Center, University of Freiburg, Faculty of Medicine, University of Freiburg, Freiburg 79104, Germany. 148 Psychiatry, Kaiser Permanente Northern California, San Francisco, CA 94115, USA. 149 Medical Research Council Human Genetics Unit, Institute of Genetics and Molecular Medicine, University of Edinburgh, Edinburgh EH4 2XU, UK. 150 Department of Psychiatry, University of Toronto, Toronto, ON M5T 1R8, Canada. 151 Centre for Addiction and Mental Health, Toronto, ON M6J 1H4, Canada. 152 Division of Psychiatry, University College London, London W1T 7NF, UK. 153 Neuroscience Therapeutic Area, Janssen Research and Development, LLC, Titusville, NJ 08560, USA. 154 Institute of Molecular and Cell Biology, University of Tartu, Tartu 51010, Estonia. 155 Psychosis Research Unit, Aarhus University Hospital, Risskov, Aarhus 8200, Denmark. 156 University of Liverpool, Liverpool L69 3BX, UK. 157 Mental Health Center Copenhagen, Copenhagen Universtity Hospital, Copenhagen 2100, Denmark. 158 Human Genetics and Computational Biomedicine, Pfizer Global Research and Development, Groton, CT 06340, USA. 159 Psychiatry, Harvard Medical School, Boston, MA 02215, USA. 160 Medical Genetics Section, CGEM, IGMM, University of Edinburgh, Edinburgh EH4 2XU, UK. 161 Psychiatry, University of Iowa, Iowa City, IA 52246, USA. 162 Department of Psychiatry and Behavioral Sciences, Johns Hopkins University, Baltimore, MD 21287, USA. 163 Department of Psychiatry and Psychotherapy, University Medical Center Göttingen, Goettingen, Niedersachsen 37075, Germany. 164 Human Genetics Branch, NIMH Division of Intramural Research Programs, Bethesda, MD 20892-9663, USA. 165 Faculty of Medicine, University of Iceland, Reykjavik 101, Iceland. 166 Child and Adolescent Psychiatry, Erasmus MC, Rotterdam, Zuid-Holland 3015, Netherlands. 167 Psychiatry, Erasmus MC, Rotterdam, Zuid-Holland 3015, Netherlands. 168 Psychiatry, Dalhousie University, Halifax, NS B3H 2E2, Canada. 169 Division of Epidemiology, New York State Psychiatric Institute, New York, NY 10032, USA. 170 Department of Clinical Medicine, University of Copenhagen, Copenhagen 2200, Denmark. 171 Department of Medical & Molecular Genetics, King’s College London, London WC2R 2LS, UK. 172 Psychiatry & Behavioral Sciences, Stanford University, Stanford, CA 94305-5717, USA. 173 NIHR BRC for Mental Health, King’s College London, London SE5 8AF, UK. 174 Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC 27514, USA. 175 Psychiatry, University of North Carolina at Chapel Hill, Chapel Hill, NC 27514, USA ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10461-0 18 NATURE COMMUNICATIONS | (2019) 10:2548 | https://doi.org/10.1038/s41467-019-10461-0 | www.nature.com/naturecommunications