Genome-Wide Association Analyses in 128,266 Individuals Identifies New Morningness and Sleep Duration Loci

Using self-reported “morningness”, we generated a binary and a continuous chronotype score. We performed genome-wide association studies on 16,760,980 imputed autosomal variants. Fig 1 presents the overall results for these GWAS. Table 1 presents details of all 16 loci associated at P<5x10^-8.

Analysing UK Biobank data with that from 23andMe provides evidence that at least 13 of the 16 are associated with chronotype. Thirteen of the chronotype signals remained at P<5x10^-8 in a meta-analysis including UK Biobank and 89,283 individuals from 23andMe [14], of which 11 reached P<0.05 in the same direction in 23andMe alone, and 15 of the 16 UK Biobank signals were in the same direction (binomial P = 0.0002) (Table 1). We also attempted to validate the associations in 6,191 European-Ancestry from the Chronogen consortium and 2,532 Korean Ancestry individuals from the Insomnia, Chronotype and sleep EEG (ICE) consortium that used “Gold standard” chronotype questionnaire (Munich Chronotype Questionnaire–MCTQ and Morningness-Eveningness Questionnaire—MEQ). Given the sample size of 5% of the discovery UK Biobank study we assessed directional consistency rather than testing for replication P-values <0.05 or 0.05/16. In the European-Ancestry individuals 11 of the 16 signals were represented. Nine of these 11 variants had the same direction of effect as the discovery UK Biobank cohort (binomial test P = 0.03) and one replicated at Bonferroni significance (rs12140153, P = 0.003). In the Korean study, 9 signals were represented, four of which had the same direction of effect as the discovery UK Biobank cohort (binomial test P = 1.00). The level of directional consistency in these two smaller studies is consistent with what would be expected in cohorts <5% the size of our discovery cohort.

A 23andMe study recently identified 15 loci associated with chronotype [14]. All of the 15 signals were replicated in our study with P<0.05 in the same direction and had meta-analysis P<5x10^-8 (S1 Table). We performed a conditional analysis of our lead chronotype variants by adjusting for the 15 known signals (S2 Table), in order to identify which of our loci coincided with those of Hu et al. [14]. Seven of our 13 replicated signals remained associated at P<5x10^-8 (see Table 1). The addition of these 7 loci brings the number associated with chronotype to 22 (full list in S3 Table).

The variant most strongly associated with chronotype, rs516134 (OR for morningness = 1.21, 95% CI [1.16, 1.26], binary P = 3.7x10^-12, continuous P = 8.9x10^-13) occurs near RGS16, which is a regulator of G-protein signalling and has a known role in circadian rhythms [16] (Table 1 and Fig 2). Another signal occurs near PER2 (lead variant rs75804782, OR = 1.09, 95% CI [1.06, 1.12], binary P = 7.2x10^-10, continuous P = 3.2x10^-7; Fig 3). PER2 is a well-known regulator of circadian rhythms [17–22] and contains a variant, rs75804782, recently shown to be associated with iris formation [23] that is in LD (r² = 0.65, D’ = 0.97) with our reported lead SNP. As there is a reported link between season and reported chronotype [24], we carried out a sensitivity analysis in which we adjusted for month of attendance (to assessment centre); all associations remained genome-wide significant for the reported variants. We tested for enrichment of specific biological and molecular pathways using MAGENTA (Meta-Analysis Gene-set Enrichment of variaNT Associations) [25] but none had a clear link to circadian rhythms (S4 Table).

We performed genome-wide association studies on a binary sleep phenotype and a continuous sleep duration score for 16,761,225 imputed variants. Fig 4 presents the overall results for these GWAS. Three loci reached genome-wide significance. The most strongly associated variant was rs62158211 with an average 2.6 minute (95% CI [1.9, 3.2], P = 5.7x10^-16) per-allele change in sleep duration and occurs at the previously reported association signal near PAX8 [13]. We identified two, novel, conditionally independent, signals that were located ~900kb apart, one upstream and the other downstream of VRK2. The downstream variant, rs17190618, has an average per allele effect of 2.0 minutes (95% CI [1.3, 2.7], P = 1.2x10^-9) on sleep duration. The upstream variant, rs1380703 (which is not correlated with rs17190618, r² = 0.002), has an average per allele effect of 1.6 minutes (95% CI [1.1, 2.2], P = 7.6x10^-9) on sleep duration. On adjusting for month of assessment, we saw similar associations for both rs62158211 (P = 3x10^-16) and rs1380703 (P = 6x10^-9), with no change for rs17190618. Table 2 shows the three sleep duration loci and their lead variants. Fig 5 shows locus zoom plots of the VRK2 association signals. We did not replicate the association of a previously reported variant in ABCC9 [12] with sleep duration (rs11046205, 0.1mins, 95% CI [-0.6, 0.7], P = 0.83).

To replicate the two novel sleep duration hits we used data from 47,180 individuals from a published study [13]. The variant rs17190618 replicated with effect size = 2.1 minutes (95% CI [0.8, 3.3], P = 0.001, meta-analysis P = 5x10^-12). The variant rs1380703 replicated with effect size = 1.3 minutes (95% CI [0.3, 2.2], P = 0.01, meta-analysis P = 3x10^-10).

Using LD-score regression we estimated the heritability of chronotype and sleep duration within UK Biobank to be 0.12 (0.007), and 0.07 (0.007), respectively. There was no evidence of a genetic correlation between sleep duration and chronotype (rG = 0.0177, P = 0.70). Chronotype was nominally genetically correlated with BMI (rG = 0.056, P = 0.048), but not Type 2 diabetes (rG = 0.004, P = 0.99). As the relationship between sleep duration with BMI and risk of T2D is U-shaped (see S1 Fig), we defined two further binary phenotypes; undersleepers (<7 vs. 7–8 hours) and oversleepers (>8 vs. 7–8 hours). There was a strong genetic correlation between undersleeping and BMI (rG = 0.147, P = 1x10^-5), but not T2D (rG = 0.022,P = 0.79). There was also a genetic correlation between oversleeping and both BMI (rG = 0.097, P = 0.039) and T2D (rG = 0.336, P = 0.001). We also performed LD-score regression analyses against a range of other diseases and traits where GWAS summary statistics are publicly available (S5 Table). Schizophrenia was genetically correlated (after adjusting for the number of tests) with hours slept (rG = 0.26, P = 5x10^-10), oversleeping (rG = 0.35, P = 6x10^-8), undersleeping (rG = -0.14, P = 2x10^-3) and chronotype (rG = -0.12, P = 2x10^-4).

The genetic correlations we observed provide general estimates that capture pleiotropic variants (those that affect both traits through different pathways) and associations that are secondary to a variant affecting a trait that causally influences the second trait. Using a genetic risk score of 69 known BMI variants [26] (listed in S6 Table) as an instrumental variable, we next performed Mendelian randomisation analyses in the UK Biobank study to test the potential causal role of BMI in chronotype and sleep. Instrumental variables analyses using variants and their effect sizes identified by previous studies [26] provided no consistent evidence that BMI causally affects self-reported “morningness” (S7 Table). Association statistics of the BMI variants with chronotype are given in S6 Table. We repeated these analyses using a genetic risk score consisting of 55 type 2 diabetes SNPs [27] and did not find any evidence of causality. Performing the reciprocal Mendelian randomization analysis using a genetic risk score of the 13 replicated chronotype variants, with effect sizes obtained from 23andMe, we found no consistent evidence in the UK Biobank data that morningness or eveningness leads to higher BMI (S7 Table). Association of the chronotype-associated variants with BMI are given in S8 Table.

Using the same genetic risk score of 69 known BMI variants as an instrument, we saw no consistent evidence that higher BMI increased an individual’s likelihood of being an undersleeper (IVreg2 P = 0.95, IVW P = 0.04) or an oversleeper (IVreg2 P = 0.29, IVW P = 0.62) in the UK Biobank data (S7 Table). Because there were only three genetic variants of small effect associated with sleep duration, we did not perform any Mendelian Randomisation analyses of sleep on BMI or type 2 diabetes risk.

This study was conducted using the UK Biobank resource. Details of patient and public involvement in the UK Biobank are available online (www.ukbiobank.ac.uk/about-biobank-uk/↗ and https://www.ukbiobank.ac.uk/wp-content/uploads/2011/07/Summary-EGF-consultation.pdf↗). No patients were specifically involved in setting the research question or the outcome measures, nor were they involved in developing plans for recruitment, design, or implementation of this study. No patients were asked to advise on interpretation or writing up of results. There are no specific plans to disseminate the results of the research to study participants, but the UK Biobank disseminates key findings from projects on its website.

We used 128,266 individuals of British descent from the first UK Biobank genetic data release (see http://biobank.ctsu.ox.ac.uk↗). British-descent was defined as individuals who both self-identified as white British and were confirmed as ancestrally Caucasian using principal components analyses (http://biobank.ctsu.ox.ac.uk↗). Of these individuals, 120,286 were classified as unrelated, with a further 7,980 first- to third-degree relatives of these. As the association tests were carried out in BOLT-LMM [31], which adjusts for relationships between individuals and corrects for population structure, we included all 128,266 related white British individuals in the association analyses.

We used imputed variants provided by the UK Biobank. Details of the imputation process are provided at the UK Biobank website (see http://biobank.ctsu.ox.ac.uk↗). For this study we only included the ~16.8M imputed variants with an imputation R² ≥ 0.4, MAF ≥ 0.001 and with a Hardy–Weinberg equilibrium P>1x10^-5.

To perform the association tests, we used BOLT-LMM [31] to perform linear mixed models (LMMs) in the 128,266 individuals. We used BOLT-LMM as it adjusts for population structure and relatedness between individuals whilst performing the association tests with feasible computing resources. As it adjusts for population structure and relatedness, it allowed us to include the additional 7,980 related individuals and therefore improved our power to detect associations. To calculate the relationships between individuals, we provided BOLT-LMM a list of 328,928 genotyped SNPs (MAF>5%; HWE P>1x10^-6; missingness<0.015) for the individuals included in the association analysis and used the 1000 Genomes LD-Score table provided with the software.

As the continuous phenotypes were derived by adjusting for age, gender and study centre, the LMM only included chip (BiLEVE vs. UKBiobank arrays) as a covariate at run-time (see http://www.ukbiobank.ac.uk/wp-content/uploads/2014/04/UKBiobank_genotyping_QC_documentation-web.pdf↗). The binary phenotypes were unadjusted and so included age, gender and chip at run-time. BOLT-LMM reported no improvement of the non-infinitesimal mixed model test over the standard infinitesimal test and so all association results reported in this paper are for the infinitesimal model [31].

Participants (N = 89,283) were from the customer base of 23andMe, Inc. The descriptions of the samples, genotyping and imputation are in [14]. Of the 16 chronotype-associated variants for which we attempted replication, 10 were available from imputation from the 1000 Genomes imputation panel phase 1 pilot. An additional 4 were imputed from the phase 1 version 3 1000 Genomes imputation panel. The final two could not be imputed. We used http://analysistools.nci.nih.gov/LDlink/↗ to find proxies—the best available were rs4729854 for rs372229746 (r² = 0.33), and rs12621152 for rs70944707 (r² = 0.33). We meta-analysed P-values from the discovery and replication samples using sample size weighting implemented in METAL [32].

Genotypes consisting of both directly typed and imputed SNPs were used for the individual GWAS [12]. To avoid over-inflation of test statistics due to population structure or relatedness, we applied genomic control for the independent studies and meta-analysis. Linear regression for associations with normalised chronotype was performed (see [12] for packages used) under an additive model, with SNP allele dosage as predictor and with age, age², gender, normalised sleep duration, season of assessment (dichotomized based on time of the year, and day-light savings time–DST or standard zone time assessments) as covariates. A fixed-effects meta-analysis was conducted with GWAMA (Genome-Wide Association Meta-Analysis) [33] using the inverse-variance-weighted method and low imputation quality (Rsq/proper_info < 0.3) were dropped from the meta-analysis.

Pathway analyses were carried out in MAGENTA[25] using all available libraries provided with the software. We included all imputed variants with association P<1x10^-5 from the continuous chronotype trait. For the results presented in S4 Table, we used gene upstream and downstream limits of 250Kb, excluded the HLA region (default setting) and set the number of GSEA (Gene Set Enrichment Analysis) [34,35] permutations at 10,000 (default). We used HaploReg v4.1[36] to annotate any coding variants within LD r² > 0.8 of the lead variant at each locus.

LocusZoom plots (Figs 2, 3 and 5) were created using the LocusZoom tool [37] (found at http://locuszoom.sph.umich.edu/locuszoom/↗) by uploading summary statistics from the Chronotype and sleep duration GWAS. For the background LD structure, we selected the “1000 Genomes Nov 2014 EUR” panel.

Genetic correlations (see [38] for methodology) between traits were calculated using the LD Score Regression software LDSC (available at https://github.com/bulik/ldsc/↗) [39]. Summary statistics of our traits outputted by BOLT-LMM were first “munged”, a process that converts the summary statistics to a format that LDSC understands and aligns the alleles to the Hapmap 3 reference panel, removing structural variants and multi-allelic and strand-ambiguous SNPs. Genetic correlations were then calculated between our phenotypes and a set of 101 phenotypes for which summary statistics are publicly available (full list in S5 Table). We used precomputed LD structure data files specific to Europeans of HAPMAP 3 reference panel, obtained from (http://www.broadinstitute.org/~bulik/eur_ldscores/↗) as suggested on the LDSC software website.

The 13 variants in Table 1 which reached P<5x10^-8 in combined analyses were used as chronotype instruments in the Mendelian Randomisation analyses. Where binary and continuous traits shared a locus, we selected the top variant of the continuous trait over that of the binary. For loci that reach GW-significance in the binary trait only, we selected the top variant but used the effect size from the continuous trait.

To test for a causal effect of BMI on chronotype and sleep-duration, we selected 69 of 76 common genetic variants that were associated with BMI at genome wide significance in the GIANT consortium in studies of up to 339,224 individuals (S6 Table) [26]. We limited the BMI SNPs to those that were associated with BMI in the analysis of all European ancestry individuals and did not include those that only reached genome-wide levels of statistical confidence in one sex or one stratum only. Variants were also excluded if known to be classified as a secondary signal within a locus. Three variants were excluded from the score due to potential pleiotropy (rs11030104 [BDNF reward phenotypes], rs13107325 [SLC39A8 lipids, blood pressure], rs3888190 [SH2B1 multiple traits]), three due to being out of HWE (rs17001654, rs2075650 and rs9925964) and the last variant due to not being present in the imputed data (rs2033529).

For testing reverse causality of type 2 diabetes on our sleep phenotypes, we used 55 of 65 common variants (listed in S6 Table) known to be associated with type 2 diabetes at genome wide significance in a meta-analysis of 34,840 cases and 114,981 [27], excluding those known or suspected to be pleiotropic.

We performed the Mendelian Randomisation analysis two ways; firstly using instrumental variables (IV) using STATA’s “IVreg2” function [40] and secondly through the inverse-variance weighted (IVW) and MR-Egger methods described in [41]. Analyses were performed in STATA 13.1 (StataCorp. 2013. Stata Statistical Software: Release 13. College Station, TX: StataCorp LP.).

In the instrumental variables method, we generated genetic risk scores (GRS) for BMI and type 2 diabetes using the published list of associated variants and their respective betas. For Chronotype, we generated a GRS using the thirteen replicated variants and their respective betas from 23andMe summary statistics. Using the IVreg2 command, we performed two-stage least squares estimation to calculate the effect of predicted exposure (through the GRS) on the continuous outcome traits. For binary outcomes (type 2 diabetes, undersleeper and oversleeper), we manually carried out the two-stage process by regressing the exposure trait on its GRS and storing both predicted values and residuals. We then used these predicted values and residuals as independent variables in a logistic regression where the dependent variable was the binary outcome.

The inverse-variance weighted (IVW) method is equivalent to a meta-analysis of the associations of the individual instruments and uses associations between the instruments and both the exposure and the outcome to estimate the additive effect of the instruments combined [41]. The MR-Egger method is a modification to the IVW method that allows the inclusion of “invalid” instruments (i.e. those that don't satisfy all three conditions), by performing Egger regression using the summary data of the variants. The IVW and Egger methods operate under the assumption that all instruments are valid, in that they satisfy the three IV conditions: the genetic variants are 1) independent of confounders, 2) associated with the exposure and 3) independent of the outcome. The MR-Egger method, however, accounts for the fact that genetic variants could be pleiotropic and may influence the outcome via pathways other than through the exposure and therefore the resulting association between genetic instruments and the outcome should not be biased by invalid instruments and pleiotropy. The MR-Egger method was used purely as a sensitivity test for the IVW method and so MR-Egger results were not considered if the IVW result did not reach nominal significance.

For the IVW and MR-Egger methods, associations of genetic instruments (variants) with both exposure and outcome phenotypes were generated in STATA by regressing the phenotype against the instrument while adjusting for covariates. As a further sensitivity test, we also repeated these analyses by replacing exposure phenotype-variant associations with their respective published betas and found only slight differences in betas and P-values, though all exposure-outcome associations remained non-significant.

The GWAS summary statistics will be made publically available on our group website http://www.t2diabetesgenes.org/data/↗. The data reported in this paper are available via application directly to the UK Biobank (http://www.ukbiobank.ac.uk/↗).