A Causal Web between Chronotype and Metabolic Health Traits

Observational and experimental evidence has linked chronotype to both psychological and cardiometabolic traits. Recent Mendelian randomization (MR) studies have investigated direct links between chronotype and several of these traits, often in isolation of outside potential mediating or moderating traits. We mined the EpiGraphDB MR database for calculated chronotype–trait associations (p-value < 5 × 10−8). We then re-analyzed those relevant to metabolic or mental health and investigated for statistical evidence of horizontal pleiotropy. Analyses passing multiple testing correction were then investigated for confounders, colliders, intermediates, and reverse intermediates using the EpiGraphDB database, creating multiple chronotype–trait interactions among each of the the traits studied. We revealed 10 significant chronotype–exposure associations (false discovery rate < 0.05) exposed to 111 potential previously known confounders, 52 intermediates, 18 reverse intermediates, and 31 colliders. Chronotype–lipid causal associations collided with treatment and diabetes effects; chronotype–bipolar associations were mediated by breast cancer; and chronotype–alcohol intake associations were impacted by confounders and intermediate variables including known zeitgebers and molecular traits. We have reported the influence of chronotype on several cardiometabolic and behavioural traits, and identified potential confounding variables not reported on in studies while discovering new associations to drugs and disease.


Introduction
With the advent of large biobank cohorts to provide genome-wide association study (GWAS), there have been several recent studies looking at the causal influence between chronotype (morningness or eveningness) exposures and both neurobehavioral/psychiatricand cardiometabolic-related outcomes. Lind and colleagues found significant genetic correlations between oversleeping, insomnia, and undersleeping exposures with an outcome of post-traumatic stress disorder [1], but when testing for causality did not find evidence for causal effects of sleep phenotypes on post-traumatic stress. Adams and Neuhausen were interested in the interplay between chronotype and free fatty acid circulation, as well as potential associations between free fatty acid circulation and type 2 diabetes [2]. So as to evaluate this, they conducted two Mendelian randomization (MR) studies using two-sample data, and found that morning chronotype is associated with lower total fatty acid levels (inverse variance-weighted estimator [IVW] β −0.21, p = 0.02) and that elevated fatty acid levels are associated with a decrease in diabetes, granting a protective effect (IVW β −0.23, p = 0.01). They then extended their analysis to include subtypes of free fatty acids and their conclusions held, indicating that a morning chronotype is associated with lower mono-unsaturated fatty acid intake. Richmond and colleagues sought to model sleep traits and risk of breast cancer using chronotype, sleep duration, and insomnia GWAS for instrumental variable selection [3]. They showed a morning chronotype to be protective against breast cancer (odds ratio [OR] 0.85). Additional analyses supported these findings, showing a morning chronotype (IVW OR 0.88) to be protective against breast cancer, while increased sleep duration had a detrimental effect (IVW OR 1. 19). Gibson investigated bi-directional causal effects of smoking on sleep duration and chronotype [4]. They found no clear evidence that smoking initiation influenced sleep behaviors directly, nor evidence for causal effects between chronotype on smoking behavior. However, they did find evidence that insomnia could lead to an increase in smoking behavior (IVW β 1.21, p = 0.02) in an under-powered analysis. Treur modelled caffeine consumption and sleep traits, including chronotype, sleep duration, and history of insomnia [5]. While the association between caffeine consumption and disturbed sleep is well known, their analysis did also show strong genetic correlations between those traits; however, their MR analyses failed to produce significant causal associations. At a wider scale, Lane and colleagues used MR analysis as a follow up to their first GWAS of chronotype using the UK Biobank [6]. They found significant associations between evening phenotype and years of education increasing and self-reported schizophrenia diagnosis, and associations between a morning chronotype and a decreased body mass index.
These studies, each taken in isolation, paint a compelling picture of relationships between chronotype and metabolic traits. However, p-values reported were generally >0.001. If analyzed as a whole set of experiments under multiple testing, results would often fail to achieve statistical significance.
Additionally, most studies fail to address confounding factors which may be identifiable by performing MR tests between exposures and potential confounders. In concurrence with these studies being conducted, large databases of GWAS and MR studies have been built, including the GWAS Catalog [7] and MR-Base [8]. In this study, we have taken advantage of GWAS study statistics and MR analyses surrounding chronotype to identify causal relationships between exposure to chronotype and several cardiometabolic-and mental health-related traits. We have additionally mined MR repositories for potential confounding factors, identifying both non-detrimental intermediate analyses (explanatory mediators and reverse mediators) and confounders which will need to be controlled for in mechanistic studies used to validate MR experiments.

Materials and Methods
A general workflow for this study is depicted in Figure 1. Using the EpiGraphDB R API v1.0 [9], results of pre-computed Mendelian randomization studies were obtained with the following parameters: exposure trait either "Morning/evening person (chronotype)" or "Chronotype", with a p-value threshold of 5 × 10 −8 (GWAS genome-wide significance). Results were filtered to significant random-effects inverse variance-weighted multi-SNP meta-analyses (IVW), with 120 significant associations retained, and additionally manually curated to keep associations between chronotype and mental or cardiometabolic health (see Table S2 for SNP characteristics). Chronotype/chronotype and chronotype/sleep duration studies were discarded, leaving 16 potential studies for investigation ( Figure 1).

Figure 1.
Mendelian randomization (MR) workflow for testing the causal influence of chronotype on traits in the EpiGraphDB database. An unbiased search for GWAS-derived two-sample MR studies was performed with chronotype or moningness/eveningness as reported exposures, and IVW results (p < 5 × 10 −8 ) retained. Potential chronotype IVs were processed with Steiger filtering, LD clumping, and stand/palindromic SNP harmonization. IVW and related MR analyses are then performed on filtered SNP summary statistics, followed by pleiotropy (MR-Egger) and leave-one-out sensitivity analyses. Associations that survive multiple testing correction are then uploaded to EpiGraphDB to investigate potential confounders not present in exposure/trait analyses as reported.

Mendelian Randomization Investigations
Each retained exposure's data were downloaded from MR-Base via the TwoSampleMR package version 0.5.5 [8]. For each instrumental variable (IV) SNP in the exposure study, the following procedure was performed. First, the strandedness of each GWAS was checked to ensure that at each allele, the minor and major alleles were equal. If these were reversed, effect sizes were modified to correct for this. Pallendromic SNPs, which contain alleles represented by the same base pairs on both strands of DNA, were discarded. If SNPs were not present, proxies were found using PLINK with an R 2 of at least 0.8, and strand was checked again [10]. Next, SNPs, in the exposure GWAS set, were clumped by LD to ensure statistical independence. In a window of 10,000 base pairs, an R 2 cutoff of <0.001 was set to obtain haplotype blocks using the European reference panel of the 10,000 Genomes Project [11]. In each exposure/outcome pair, this left a variable number of SNPs for use as valid, independent IVs (see Figure 1 top second step). Effect sizes for each SNP were reported as a β or the transformed log(OR). The Wald ratio was then obtained, giving a measure of the effect of the exposure on the outcome [12]: where β_Y_j is the effect of the IV on the outcome, and β_X_j is the effect of the IV on the exposure is obtained for SNP j, and θ is the effect size. An initial IVW analysis between each exposure/outcome set was performed, with the false discovery rate (FDR) controlled for [13]. An FDR of <0.05 was considered significant. Rather than calculate Wald ratios individually, the outcome GWAS βs or odds ratios are regressed on the exposure in an inverse variance-weighted (IVW) meta-analysis. The slope of the regression line indicates the strength of the effect, as an increase in the unit of outcome per unit of the exposure [14]. In a IVW meta-analysis, the IVW estimate is calculated by: where θ is the inverse variance-weighted average, se is the standard error, is an error term ( Figure 1 third step, right), and other terms are as above. Ten significant IVW studies were then subject to four additional MR methods ( Figure 1, third level left). To address unseen horizontal pleiotropy, the Egger regression was performed [15,16]. The Wald ratios of each SNP are used in meta-regression by taking the inverse variance weights used in the IVW analysis without modeling the intercept. As a result, a causal estimate, similar to IVW, is obtained, adjusted though for horizontal pleiotropy which is necessary so as to ensure IVW validity [16]. MR-Egger regression is an extension of IVW regression. Instead of assuming no intercept term, an intercept is estimated: where θ 0E is the intercept and θ 1E the MR-Egger estimate. If the intercept is equal to zero, then the IVW method and MR-Egger will be equivalent [14]. During the IVW process in MR-Egger, the effect sizes of each SNP must have the same sign, and this decreases the variation between them [16]. Inverse variance-weighted median analyses were performed, using the median of the Wald ratios as an effect size. Unweighted analyses assume that over half of the instruments are valid, while in a weighted median analysis the assumption is that the at least 50% of the weight of the instruments are valid themselves [17]. This approach is robust to directional pleiotropy when compared to a simple IVW meta-analysis.
The mode-based estimator (MBE) clusters Wald ratios before calculating random effects in an IVW meta-analysis [18]. The simple MBE uses unweighted analysis, while the weighted MBE uses inverse variance weighting. First, a smooth empirical density function is calculated for each Wald ratio and are then clustered. The Zero Modal Pleiotropy Assumption states that the biggest cluster with the same ratio estimates will be valid instruments. These were used in each analysis, resulting in fewer SNPs and less power, but an increase in potential robustness to horizontal pleiotropy.
Heterogeneity can occur when individual SNPs do not converge on an estimate; this was estimated by Cochran's Q [19]. In this context, heterogeneity may be a sign of horizontal pleiotropy, wherein SNPs effect the outcome by their influence on other confounding traits [20]. To inspect SNPs for outliers, we performed leave-one-out sensitivity analyses using the IVW method, leaving out one SNP in each analysis.
We used the Steiger test to access the directionality of all causative analyses post hoc [21]. The Steiger test first assesses which variables (exposure or outcome) are influenced by the SNPs used, by testing if the SNPs explain more variance in the exposure than in the outcome with a modified Z statistic. If the p-value of the IVW estimate and the Steiger estimate are both significant, the sign of the Z statistic is used to assign the correct causal direction between exposure and outcome.

Confounder and Intermediate Analysis
To address potential confounders, not seen in pairwise analyses even in the absence of the statistical suggestion of horizontal pleiotropy, existing MR studies were obtained from EpiGraphDB (Figure 1 last step). Each study exposure/outcome pair with significant FDR-corrected IVW results was used to interrogate the database for intermediate, reverse intermediate, collider, and confounding variables. Results from the study database surpassing a p-value threshold of <1 × 10 −5 , using a fixed-effects IVW method only, were retained. The β effect sizes for significant results between exposures and outcomes, exposures and confounders, and outcomes and confounders were used to create a weighted directed graph in Cytoscape [22] version 3.8.2, and the yFiles Organic (force directed) layout to visualize the relationships between traits.
All statistical analyses were performed in R version 4.0.5 [23].

Results
Initially, 120 prospective studies were identified using EpiGraphDB containing chronotype as an exposure,. Of these, there were 28 primary associations, with IVW p-values of <5 × 10 −8 , directly relevant to this study, including measures of alcohol intake, bipolar disorder, T2DM (type 2 diabetes mellitus) as well as testosterone and triglyceride levels. See Table S1. Since the analyses were performed in a high-throughput environment, we re-analyzed each relevant study with the IVW method using summary statistics ( Figure 2). Following filtering and quality control, 10 studies held up to a false discovery rate of p < 0.05. The bi-directional nature of the forest plot exists due the different nature of the chronotype exposure studies (treating morning or evening chronotype as case/control). Each of the 10 studies were further analyzed with a tendency toward eveningness reflecting a high β value in the exposure effect size. Outcome   Figure 2. Causal relationships between exposure to chronotype and various traits. Causal associations mined from EpiGraphDB (IVW p < 5 × 10 −8 ) were analyzed using two chronotype exposure measures. Outcomes are listed with trait and study accession number. Exposure ebi-a-GCST003837 reflects SD increase in chronotype on a continuous morning-evening scale. Exposure ukb-a-11 represents odds of an evening chronotype. Effect sizes are from random-effects inverse variance-weighted analyses with 95% confidence intervals. False discovery rate-corrected p-values are shown.

Chronotype Influences on Diabetes, Alcohol Consumption, and Bipolar Disorder
The strongest associations (β > 0.90) include risk for T2DM and total fatty acid concentration. T2DM included seven independent SNPs, see Table S3. The weakest association was with the MR-Egger method, β = 0.73, which while not highly significant (p = 0.04) did not reveal evidence of horizontal pleiotropy (intercept p-value = 0.98) or heterogeneity (Egger Q p-value 0.99). The IVW analysis and weighted median analysis each concurred (p = 0.002 and 0.009), suggesting little loss of associative signal even if some of the seven SNPs were biased by pleiotropy (Figure 3a). Leave-one-out sensitivity analyses suggests that no one IV dominates the model, and all methods have similar effect sizes as judged by slope, as shown in Figure 3b.
Chronotype has a weaker but more statistically significant association with the amount of alcoholic beer or cider drinks consumed weekly, with the increase reflected in pints per week consumed (mean among all respondents of 3/week). With an MR-Egger intercept p-value of 0.31 and a significant Egger regression coefficient (p-value = 0.04), the initial analysis provides strong evidence that the trait may not be subject to horizontal pleiotropy. The weighted median and IVW method were again strongest (β = 0.089 and 0.79, p = 4.4 × 10 −6 and 9.7 × 10 −4 ), as shown in Figure 4a. There was, however, evidence among the 82 SNPs which passed the threshold for strong heterogeneity (Cochran Q = 333), suggesting pleiotropy or moderating factors for further investigation even in the absence of such suggestions from the MR-Egger intercept significance test. Heterogeneity among SNPs can be observed in Figure 4b, while the leave-one-out analysis shows overlapping confidence intervals in each iteration, suggesting homogeneity on the whole.    The influence of chronotype on bipolar disorder appears substantially varied, with different SNPs with large standard errors apparent in Figure 5. The Cochran's Q statistic (Q = 149, p = 2 × 10 −24 ) echoes the visible heterogeneity of SNPs in the analysis. Nevertheless, there is a consistent β effect size among the IVW, weighted median, and weighted mode analyses (0.188, 0.175, and 0.181) with are highly statistically significant (p = 4 × 10 −4 , 1 × 10 −11 , and 3 × 10 −6 , respectively). This presents a small but powerful signal that an evening chronotype leads to an increase in activation of pathways leading to bipolar dis-

Confounder Case Studies: Bi-Polar Disorder and Alcohol Intake
Each of the 10 exposure/outcome analyses were included in a confounder analysis, searching via existing MR associations for possible intermediate variables. Results in full are available in the Table S4. In Figure 6, a directed graph was created to visually inspect the relationships between types of confounders, outcomes, and potentially causal exposures. Exposures are represented in green, reflecting the two independent GWAS analyses used to access chronotype. Outcomes are reflected in red, and potential intermediate or confounding traits as white notes. As may be evident, there are many potential confounders (111), followed in number by intermediates (52), colliders (35), and reverse intermediates (18) all with p < 1 × 10 −5 in the EpiGraphDB database. The graph shows two clusters sharing potential confounding traits (time to first cigarette and lung cancer). Both are confounders between chronotype and two outcomes: cigarettes confound omega-3 fatty acid concentration and waking too early; lung cancer also confounds omega-3 fatty acid wile also impacting average weekly beer/cider intake.
Not all MR analyses are beset by confounders, however. The relationship between exposure to an evening chronotype and bipolar disorder was not found to be influenced by a statistically significant confounder or collider. In addition to the direct relationship between chronotype and bipolar disorder, a mediator of the odds of ER+ breast cancer was discovered, see Figure 7a. The relationship from chronotype to cancer (β 0.32) was stronger than the unmediated relationship directly to bipolar (β 0.17), but weaker than the possible relationship of breast cancer to bipolar disorder (β 0.42).

Figure 6.
Causal and confounding relationships between chronotype exposures (green nodes) and traits (red nodes). White nodes are potential confounding variables. Edges indicate the direction of causality (arrows). Green edges = causal relationships; purple = colliders, orange = confounders, brown = intermediates, and blue = reverse intermediates.
As a last example study, the relationship between exposure to an evening chronotype and an increase in beer/wine intake was explored. There is a possible bi-directional relationship between chronotype and beer/alcohol intake, mediated by reverse intermediates of a lack of physical activity (types of physical activity in the last 4 weeks: none of the above) and concentration of the UBE2G2 protein. Additionally, direct intermediates include the lifetime number of sexual partners. No colliders were present, but several potential confounders were identified which may need to be conditioned on, depending on study design. These included dietary consumption (milk type and coffee consumed), loneliness and isolation, hip circumference, and other metabolite measures. See Figure 7b. Data underlying the confounder graph in Figures 6 and 7 are in the Supplementary Table S4. (a) (b) Figure 7. (a) An evening chronotype associates with bipolar disorder. In addition to a valid association between chronotype and bipolar disorder, an intermediate association (p < 1 × 10 −5 ) was found for ER+ breast cancer. (b) An evening chronotype associates with increased beer and cider consumption. There is a potential bi-directional relationship between increased alcohol intake and an evening chronotype both directly (blue lines) and mediated by lack of physical activity and elevated UBE2G2. Chronotype may influence the number of sexual partners, which may influence average alcohol intake (brown). There are numerous potential confounders relating body composition, loneliness, and diet.

Discussion
This analysis characterised the influence of an evening chronotype on cardiometabolic and behavioral traits. Chronotype, a measure of circadian rhythm, is a largely endogenous process though there are sustaining queues, zeitgebers, which influence chronotype [24]. The translational importance of understanding circadian biology is not limited to neurobehavioral function; the interplay between metabolism and circadian biology has been highlighted heavily in recent years. The gut microbiome, adipose cytokines, and metabolic hormones from ghrelin to leptin are all strongly regulated by circadian biology [25][26][27]. Recent Mendelian randomization studies have also suggested a strong causal link between chronotype (a gross circadian phenotype) and body composition, free fatty acid circulation, and adiposity [2,28]. The current analysis revealed strong associations between chronotype and 10 other biological traits. In Figure 2, we reproduced causal associations between diabetes and triglycerides. Additionally, we report novel associations between chronotype and metabolomic data (omega-3 fatty acids, total fatty acids, bilirubin). The strongest signals include diabetes and bipolar disorder, highlighting chronotype as a link between metabolic and psychiatric health. While not all IVW MR studies passed multiple testing correction, those that did had complex relationships with other traits, various degrees of heterogeneity and potential interplay in horizontal pleiotropy. All serum triglyceride associations and a link to type 2 diabetes had the same relationship to chronotype: increased risk or concentration of the trait on an exposure to a genetic predisposition to eveningness. Interestingly, there was a significant (albeit small) association between eveningness and a decrease in likelihood to be taking nicorandil, a ventilator used to treat angina (see Figure S4). Nicorandil has been shown in one study to affect the diurnal rhythms in body temperature and heart rate in rats [29]. There is no evidence that the drug affects circadian behavior in humans, although variant angina (treated with nicorandil) has a circadian component, supporting the idea that the influence of circadian biology may contribute to some types of angina that are amenable to nicorandil treatment [30]. The influence of chronotype on alcohol consumption was recently investigated by Hisler and colleagues, who revealed a 24 hour rhythm of alcohol craving in late adolescent adults and showed that an evening chronotype corresponds to a later craving for alcohol [31]. The associations between an evening chronotype and bipolar disorder are well documented [32], and has been proposed as an endophenotype to classify subtypes of bipolar [33], though several interacting traits may play a role in these subtypes. Not surprisingly, there were over 200 potential confounders or mediators observed between chronotype and the ten significant exposures studied. As evidenced by Figure 6, they vary in type from acceptable in an MR analysis (intermediate, possibly reverse intermediate depending on follow-up studies) to those that must be conditioned upon to maintain causality (the most abundant intermediate, confounders) to the more problematic colliders. Colliders are so called because the causal arrows in a directed acyclic graph from exposure and outcome both impact (collide on) the variable. Unlike confounders, conditioning on colliders in regression introduces bias into the association between exposure and outcome [34]. The four associations subject to collider bias are the relationship between chronotype and omega-3 fatty acids, serum total triglycerides, total fatty acids, and lethargy. Possible confounders for the three metabolites include other cardiometabolic parameters, including total fatty acids, fasting insulin, and VLDL concentrations. The release of metabolomic data in the UK Biobank [35] caters for new opportunities to investigate these collisions experimentally since they relate to chronotype and both psychiatric and cardiometabolic diseases. The only collider associated to lethargy was narcolepsy, an intuitive yet potentially complex relationship. The only confounder between an evening chronotype found in this study was ER+ breast cancer. The confounding was as an intermediate, which maintains the causal nature of the analysis while adding an extra explanatory "path" to get from exposure to outcome. While there is no current literature directly linking chronotype to breast cancer and then to the presence of bipolar disorder, this may be investigated in current murine models of breast cancer via knockdown experiments or by manipulating light as an external zeitgeber.
Instrumental variable analysis relies on strong assumptions, and occasionally unverifiable conditions [36]. Subject knowledge is the most reliable method for concluding assumptions are valid, especially when choosing exposures and outcomes correctly. A naive, phenome-wide search for causation between traits is likely to include among many tests non-plausible experimental designs, leading to an artificial increase in multiple testing burden and design flaws. This is evident by the need to prune associations in the first step of this study. In our re-analysis of traits in step 2, we failed to replicate the extremely low p-values present in the EpiGraphDB, even when using the same outcome and exposure studies obtained by a sister database, MR-Base. Differences may be partly explained by analysis choice in the IVW models. We modeled a random effect, assuming that while each SNP contributed overall to the model in a similar direction or magnitude, the Wald ratio generated by each SNP will not be identical. In other words, if each SNP is analogous to a "study" in a meta-analysis, the effect sizes and variances may be similar but a model would not assume they are identical. As a large focus of this study was investigating possible confounders between traits and the possibility of pleiotropy among SNPs, a more conservative random effects approach was chosen as the basis for further analyses. Additionally, the EipGraphDB calculated several models and filtering strategies, and an optimal method for each MR analysis was chosen from a random forest model optimized on artificial data. When re-analyzing data, we chose one strategy, though we included any of their fixed-effect models when incorporating pre-calculated confounder MR analyses as this was the most abundant analysis method in the database. Without access to the underlying filtering methods used in the large repositories, we can speculate that instrumental variable selection, LD clumping, and filtering also played a part. However, this re-analysis has produced associations that stand up to multiple testing while revealing novel associations. This study was limited to viewing causality through the lens of Mendelian randomization alone. In our previous work, we have used gene expression data to identify causal regulatory networks using Bayesian approaches [37]. Opportunities exist to combine these approaches, creating multimodal graphs of gene regulatory networks from that approach with the multitrait networks created in this work.

Conclusions
Our findings have revealed potential causal associations between chronotype and several behavioral and metabolic traits. While using publicly available data to confirm previously known associations (T2DM, cholesterol) we have found novel associations which may be validated experimentally in models. These include associations between chronotype and the vasodilator nicorandil, and a relationship between chronotype and bipolar disorder potentially mediated by estrogen-receptive breast cancer. Lastly, we have revealed potential confounders and colliders impacting the relationship between chronotype and commonly reported cardiometabolic traits which should be addressed through a combination of multivariate and multistage analyses as appropriate.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: