Case-control study

We studied participants from the Danish national registry, which is based on record linkage between the Danish Psychiatric Central Register, the Danish Civil Registration System14 and the Danish Newborn Screening Biobank15. Newborn dried blood spots have been systematically collected and stored since May 1, 1981 and stored at −20 degrees C. Cases were randomly selected from all those born in Denmark between 1981–2000 who received a diagnosis of schizophrenia according to International Statistical Classification of Diseases, 10th Revision16 code F20. Most of these cases developed schizophrenia between September 2005 and December 2008. Controls, drawn from the same birth cohort, were individually matched on sex and date of birth, and were alive and free of schizophrenia at the time of onset of the matched case. The Danish Data Protection Agency and the ethics committees of the Central Denmark Region approved this study, and all methods were performed in accordance with relevant guidelines and regulations.

Assessment of vitamin D status

3.2 mm DBS (samples and calibrants) were hydrated by the addition of 50 µL water and shaken for 30 min at room temperature. Vitamin D species were extracted with acetonitrile (500 µL) containing 6,19,19-[2H 3 ]-25OHD 2 and 6,19,19-[2H 3 ]-25OHD 3 (0.5 pmol) as internal standards. After mixing, supernatants (500 µL) were subjected to solid-phase ion exchange using zirconium dioxide (ZrO 2 ) and titanium dioxide (TiO 2 ) (Glygen, USA) to remove ion-suppressing phospholipids. Eluted samples were evaporated to dryness then derivatised with 4-phenyl-1,2,4-triazioline-3,5-dione (PTAD), (50 µL, 0.1 mg/mL in anhydrous ethyl acetate). After 30 min samples were again evaporated to dryness.

Samples were reconstituted in 33% acetonitrile/water (60 µL) on a plate shaker at room temperature (650 rpm, 1 min). Samples were injected onto a Kinetex column XB-C18, 50 × 2.1 mm, 1.5 µm (Phenomenex, CA, USA) by a peak-focusing technique (30 µL sample injection volume +20 µL water from an open 2-mL glass vial) using an integrated pretreatment program. The LC-MS/MS system consisted of a UHPLC (Nexera X2, Shimadzu Ltd, Japan), comprising a SIL-30AC autosampler (50 µL loop), two LC-30AD binary pumps, DGU-20A5 degasser, and CTO-30A oven, and interfaced to a 5500 QTRAP mass spectrometer (SCIEX, Canada) with Atmospheric Pressure Chemical Ionisation (APCI) TurboIonSpray ion source. Reference samples for neonatal dried blood are not available. Therefore we prepared an “in house” dried blood spot sample for inter-assay variance. The inter-assay variation over 2 years was 11.6%17. The sera calibrators supplied by the National Institute of Standards and Technology (NIST 1950) were used in each plate. Using these external calibrants, overall assay inaccuracy was <4%. Previous studies have demonstrated a strong correlation between this dried blood spot assay and cord sera 25OHD concentration (r = 0.85, P < 0.0001)18. Full details of the assay are described elsewhere17.

Total 25OHD was reported as the sum of 25-hydroxyvitamin D 2 (25OHD 2 ) and 25-hydroxyvitamin D 3 (25OHD 3 ) species as previously described, and the inter- and intra-assay coefficients of variability were 6.9% and 11.6% respectively17.

Genotyping and generation of the PRS

The genotyping and quality control details for this sample (henceforth referred to as DK2016) have been published previously19. In summary, DNA was extracted from the neonatal dried blood spots stored at the Statens Serum Institut, whole-genome amplified in triplicate using the Qiagen REPLI-g mini kit [Qiagen, Hilden, Germany] (the three separate reactions were pooled), and then genotyped with Illumina Human 610-Quad BeadChip array (DK2016a) or Illumina Infinium CoreExome beadchip (DK2016b)(Illumina, San Diego, CA). In order to maximizes sample size, we combined our new DK2016 sample with a previously published independent case-control sample (henceforth referred to as DK2010)11. The original samples for DK2010 were 423 cases and 425 controls. The original samples for DK2016a and DK2016b were 387 cases and 391 controls, and 923 cases and 924 controls respectively. 8 subjects were removed for ambiguous sex (genotype versus register values) and 1 subject was removed for a very high concentration of 25OHD 2 . The models that examined the association of 25OHD on risk of schizophrenia was restricted to matched pairs, which are a subset of these samples.

We generated PRS in our Danish sample using allelic effect sizes estimated in the latest GWAS for schizophrenia from the Schizophrenia Working Group of the Psychiatric Genomics Consortium20. The Danish samples were excluded from this ‘discovery sample’, leaving a sample of 34,600 cases and 45,968 controls. Specifically, we selected single nucleotide polymorphisms (SNPs) that were present in both our combined dataset and the PGC data. SNPs were also filtered on PLINK imputation INFO parameter retaining those with values > 0.8 in both data sets. Next, SNPs were pruned based on linkage disequilibrium (LD) r2 (r2 > 0.2 in 500 kb window) using the clumping function in PLINK 1.0721 (i.e. preferentially retaining the most associated SNPs in an LD regions), and we excluded SNPs in the extended Major Histocompatibility Complex (MHC) region (chr6:25–34 Mb; except retaining rs7746199). We created a single PRS for each individual based on SNPs with p-value < 0.05; this decision was justified as the threshold that maximised variance explained in leave-one-out analyses20. These algorithms were applied to datasets within each genomic array and the resultant PRSs were standardised within each chip and sex and then merged.

Statistical analyses

In keeping with our previous analyses of neonatal blood samples11, we derived quintiles for 25OHD based on the control sample, and used conditional logistic regression (Incidence Rate Ratio; IRR, and 95% confidence intervals) to assess the relationship between neonatal 25OHD concentration and risk of schizophrenia. As our previous study identified a non-linear relationship (both the lowest and highest quintiles were associated with increased risks of schizophrenia compared to the fourth quintile), we prespecified the fourth quartile as the reference category. Based on previous research using the same psychiatric case register and based on factors known to be associated with vitamin D status22 we explored the influence of a range of variables on the association between the variables of interest (using likelihood ratio tests). These include maternal, paternal, and sibling history of any mental disorder, sex, age second-generation immigrant status, degree of urbanization at place of birth, maternal and paternal age at the time of the child’s birth, gestational age, birth weight, and birth length. Because the cases and controls were individually matched on sex, date of birth/age, all relative risks were controlled for these variables. Population-attributable fraction was calculated according to the recommendations of Bruzzi et al.23 (equation 10).

For the analyses that examined both 25OHD concentration and PRS scores in the combined samples, the 25OHD concentrations were first standardised within the assay runs (DK2010 and DK2016) and transformed using Van der Waerden-normalised scores24. The range and distribution of 25OHD concentrations were similar between the two samples. A sample of QC-ed and LD-pruned SNPs (approximately 40,000) were used to calculate principal components (PCs) in PLINK 1.925 and ancestral outliers mean +/−4 SD for each of the first two PCs (loss of 414 individuals giving a total sample of 3129). After excluding these outliers, PC axes were recalculated. None of the PC axes were associated with schizophrenia status (up to 10 PCs tested). However, PC1 was significantly associated with both PRS and 25OHD concentration.