Genetic influences on circulating retinol and its relationship to human health

The common and rare variant genetic architecture of circulating retinol

We integrated common and rare variant data from up to 22,274 individuals of European ancestry in our discovery meta-analyses to estimate genetic effects on circulating retinol (Figure 1A-B, Online Methods). Firstly, variants from the INTERVAL and METSIM studies were meta-analysed (N_Meta = 17,268 – termed METSIM+INTERVAL, Figure 1, Online Methods). After harmonisation, there were 8,173,975 common and 5,091,050 rare [minor allele frequency (MAF) < 1%] overlapping variants between INTERVAL and METSIM, respectively. As retinol effect sizes were estimated in the same units in both studies (plasma SD units, quantified by the same instrument), we conducted both a fixed-effects inverse-variance weighted (IVW) meta-analysis and a sample-size weighted meta-analysis of Z scores (Stouffer’s method). Secondly, we conducted an additional meta-analysis (N_Meta = 22,274) which also included data from two other studies (ATBC and PLCO), termed the METSIM+INTERVAL+ATBC+PLCO meta-analysis. However, there were markedly fewer variants available in this meta-analysis after imputing the ATBC and PLCO summary statistics and harmonising with METSIM+INTERVAL (N_Var = 3,896,351, Online Methods). As a result, we focused on the METSIM+INTERVAL meta-analysis as the primary discovery dataset.

• Download figure
• Open in new tab

Figure 1. Common variant influences on circulating serum retinol.

(a) Manhattan plot of the meta-analysis of common variants shared between the INTERVAL and METSIM cohorts (Stouffer’s sample size weighted meta-analysis). Variant-wise -log₁₀ P-values for association are plotted, with the dotted red line denoting genome-wide significance. The closest genic transcription start site is labelled for each lead SNP. (b) Manhattan plot, as above, for the larger sample size meta-analysis that includes the ATBC and PLCO cohorts, but with fewer variants available for meta-analysis (c) Estimates of SNP heritability of retinol (h²), with the error bars denoting the standard errors of the estimates. The first two panels denote estimates using the BLD-LDAK model and the LDAK-thin model, respectively, both using LD tagging files from the Great British ancestry participants in the UK Biobank. The last panel estimates heritability using the LDSR model with LD from the 1000 genomes European participants. Estimates were for the METSIM+INTERVAL meta-analyses (Stouffer’s and IVW), as well as the larger meta-analysis including ATBC/PLCO. (d) Empirical Bayes’ estimation of non-null effects on retinol genome-wide, stratified by bins of ascendingly sorted LD score by magnitude. The LD score bins were different for each panel – 1000 bins, 1000 genomes European LD scores (top left); 5000 bins, 1000 genomes European LD scores (top right); 1000 bins, UKBB white British LD scores (bottom left); 5000 bins, white British LD scores (bottom right). Each point represents the proportion of non-null effect sizes for that bin, with the trendline estimated using a generalised additive model for the relationship between the LD score bin and the proportion of non-null effects.

Considering common variants from the HapMap3 panel (MAF > 0.05, outside of the major histocompatibility complex (MHC) region), we observed a relatively subtle inflation of retinol signals across the genome as indexed by the mean X² statistic, with mean X² values around 1.03-1.04, regardless of the meta-analysis considered (Supplementary Table 1). Common variant SNP heritability (h²_SNP) estimated using the linkage disequilibrium (LD) score regression (LDSR) approach and the 1000 genomes European reference panel was between 6%-7% and nominally statistically significant, although with somewhat large standard errors (2.6%-3.2%) (Figure 1C, Supplementary Table 2). SNP heritability estimates of circulating retinol increased to between 10%-13%, but were still noisy, using two alternate models to estimate h²_SNP and the UK Biobank (UKBB) as the LD reference (Figure 1C, Online Methods, Supplementary Table 2). Partitioned h²_SNP across tissues and cell-types demonstrated nominal enrichment in biologically logical contexts such as liver, adipose, pancreas, and blood (Supplementary Figure 1). The somewhat large h²_SNP standard errors are likely a product of sample size; however, we then conducted further analysis to explore the extent of the polygenic signal associated with retinol across the genome using an Empirical Bayes’ method (Online Methods). This method was utilised to model the number of non-null effects on retinol genome-wide stratified by bins of LD scores that index the extent of LD a variant exhibits with other variants (Figure 1D). Across all LD score bins, we estimated that the mean fraction of common variants across the genome with non-null effects on retinol was between 1.4% to 2.4%, depending on the modelling parameters used. In line with expectation, the proportion of non-null retinol effects was very high (> 50%) when considering the variants that display the most extensive LD (highest LD score bins). The application of these analyses to two GWAS of another vitamin (25-hydroxyvitamin D₃) with either comparable or much larger sample sizes^{23, 25}, suggested that the less-polygenic architecture of retinol observed in this study may become more diffuse across the genome with greater sample sizes, in line with many other quantitative traits (Supplementary Figure 2). However, we caution that further investigation will be required as these data become available to confirm this.

Next, we processed the common variant results of both meta-analyses to identify genome-wide significant loci associated with circulating retinol (P_GWAS < 5 × 10^-8). In the primary discovery meta-analysis (METSIM+INTERVAL, Stouffer’s method), we uncovered eight genome-wide significant loci, six of which were not reported in the previous Mondul et al. retinol GWAS (Table 1). The absolute effect sizes of these lead SNPs were between 0.066 and 0.172 SD in circulating retinol per effect allele, as derived from the IVW meta-analysis. We observed minimal heterogeneity between the two cohorts for these lead SNPs, although heterogeneity was slightly more marked for rs6601299. Replication was then attempted in the TwinsUK cohort (N up to 1621, Online methods). Considering the mean association across all timepoints retinol was measured, as well as the twin pairs separately, 7 out of the 8 lead SNPs in the loci had effect sizes in the same direction, which was greater than expected by chance alone (Binomial P = 0.035, Supplementary Table 3).

We also investigated the effect of retinol associated lead SNPs on factors associated with dietary intake of retinoids. Using a GWAS of retinol intake derived from a self-reported 24-hour dietary recall in the UK Biobank (UKBB, N = 62,991)²⁶, we found no evidence to suggest that any of the effect of these genetic signals on circulating retinol is mediated through influencing dietary intake behaviours (Supplementary Table 4, Supplementary Figure 3).

In the larger sample-size meta-analysis with fewer variants available across all input datasets (METSIM+INTERVAL+ATBC+PLCO), six of the eight genome-wide significant loci from the smaller meta-analysis were available to test for association. Five of the six loci available became more statistically significant in this larger meta-analysis, whilst the chromosome 8 locus obtained a very similar level of statistical significance in both meta-analyses (Table 2). It should be noted that there was relatively large heterogeneity for the lead SNP at the TTR locus on chromosome 18 locus between cohorts. This was due to a more significant association in the ATBC+PLCO GWAS, although this region was still very strongly associated in METSIM and INTERVAL in the same direction.

We also estimated rare variant (frequency < 1%) effects on circulating retinol using variants available in the METSIM+INTERVAL meta-analysis. Despite relatively low power to detect rare variant association, we identified a genome-wide significant rare variant signal on chromosome five (chr5:86765041:T:C, dbSNP ID: rs138675130) associated with reduced circulating plasma retinol per C allele (-0.441 SD, SE = 0.0709, P = 6.37 × 10^-9) with no significant heterogeneity between the contributing studies. The frequency of this C allele in Europeans (gnomAD v.3.1.2) ranges from 0.5% in non-Finnish Europeans to 0.8% in Finns, and whilst it is marginally rarer in Africans and South Asians, it is entirely absent in the East Asian and Middle Eastern populations in that database. The variant is intergenic and in FinnGen release 8, the C allele was associated at phenome-wide significance (P < 1 × 10^-5) with increased odds of benign neoplasm of the eye and adnexa, as well as whooping cough. The closest canonical transcription start site to this variant is that of COX7C, which encodes a subunit of a terminal component of the mitochondrial respiratory chain. We also uncovered several suggestively significant rare variant signals (P < 1 × 10^-5), including three non-synonymous variants in the genes FREM2, NAXD, and CHD1 (Supplementary Table 5). Of these, the rare NAXD non-synonymous allele suggestively associated with lower retinol (rs3742192) had some in silico evidence to suggest deleteriousness (Supplementary Table 5), although it is classed as benign in ClinVar. Finally, to boost power, we also statistically aggregated rare variants to genes (Online Methods). There were no significant retinol genes after Bonferroni correction of these gene level association results (P < 3.02 × 10^-7); however, there were two novel suggestively associated genes (P < 3.02 × 10^-5), GALM and ZDHHC18 (Supplementary Table 6).

Prioritisation of retinol-associated genes reveals novel mechanisms influencing circulating retinol

In common variant loci, prioritising causal genes can be difficult due to confounding factors like linkage disequilibrium. For circulating retinol, we employed a multi-faceted approach to prioritise genes that are confidently associated. Firstly, we sought to interrogate the eight genome-wide significant loci uncovered in the main discovery meta-analysis (METSIM+INTERVAL) to uncover plausible causal genes. This was achieved by adapting an integrative pipeline developed in previous work that considers annotation, probabilistic finemapping, integrative scoring, and in silico prediction (Online Methods, Supplementary Table 7). We describe these results further for each locus in the supplementary text but summarise the prioritisation evidence forthwith. There was quite consistent evidence in four loci for a likely causal gene (GCKR, FOXP2, RBP4, and TTR). RBP4 (chromosome 10 locus) and TTR (chromosome 18 locus) were previously associated with retinol in the only other dedicated GWAS of this trait and form the complex that transports retinol in serum²², thereby, having a direct biological link to retinol abundance. GCKR, a gene encoding a protein that binds to and regulates the key metabolic enzyme glucokinase, was confidently the causal gene for the locus on chromsome 2 with the rs1260326 lead SNP. This gene is known to have a large and varied metabolic association profile due to the role of glucokinase in glycaemic and lipid related processes, amongst others^27–29. Interestingly, the lead SNP, and most likely causal variant derived from probabilistic finemapping, was a common missense allele (rs1260326), whereby the retinol increasing T allele corresponds to a substitution of leucine for proline in the GCKR protein. Previous experimental work suggests that this variant impacts GCKR affinity for glucokinase^{27, 30}.

The transcription factor gene FOXP2 was also strongly supported by multiple lines of evidence as a causal gene for the locus on chromosome seven. The role of this gene in the brain and in relation to neurological phenotypes like language has been extensively studied³¹, but less so in the periphery despite relatively high expression across many different systemic tissues. Therefore, we analysed RNA-sequencing data of FOXP2 overexpression in a cell-line not derived from the central nervous system (human osteosarcoma epithelial cell line) and revealed the transcriptional correlates of FOXP2 overexpression were enriched for a broad range of pathways related to factors including extracellular matrix biology, glycosylation, and interleukin signalling, amongst many others (Supplementary Tables 8-11, Online Methods).

The remaining four loci exhibited less clear evidence of which gene to prioritise, although all point to potentially interesting functional mechanisms. On chromosome eight, there is some evidence to support PPP1R3B as a gene that influences retinol, which encodes a catalytic subunit of the phosphatase PP1 that is implicated in relevant metabolic processes like glycogen synthesis³². However, other lines of evidence point to a role of long-noncoding RNA in this locus. The loci on chromosomes 16 and 20 are noteworthy as the closest transcription start sites to the respective lead SNPs are two transcription factors (TF) from the Maf family (MAF and MAFB). Interestingly, MAFB has been shown to regulate both TTR and RBP4 expression in various tissue contexts from human or murine studies^{33, 34}. As only some lines of evidence support these two TF, further functional characterisation of these two loci is warranted. Interestingly, in the locus on chromosome 16, two of the other genes with some evidence for a retinol-related function (MAFTRR and LINC01229) have been shown experimentally to regulate MAF expression and are also associated with other biochemical traits like urate, further supporting the role of the Maf family on retinol abundance in serum³⁵. Finally, the remaining locus on chromosome 2 (2:122078406-122084285) had the least interpretable functional prioritisation results. The closest transcription start site to the lead SNP was another TF (TFCP2L1) that has broad physiological roles including in the kidney³⁶.

We then sought to expand our scope for gene discovery beyond genome-wide significant retinol-associated loci through further integration of genetics with transcriptomics and proteomics (Online Methods, Supplementary Tables 12-13). Firstly, we leveraged multivariate models of genetically regulated expression (GReX) to perform a transcriptome-wide (liver, whole blood, adipose, small intestine, pancreas, and breast mammary tissue) and proteome-wide (plasma) association study (TWAS/PWAS) of circulating retinol. Tissues for the TWAS were selected based on prior knowledge of retinol biology and the results of the partitioned heritability analyses (Online Methods), whilst plasma was the only relevant tissue available for PWAS. After applying multiple-testing correction to the TWAS and PWAS individually (FDR < 0.05) and testing whether there was a shared causal variant via colocalisation [Posterior probability (PP) of a shared causal variant (H4), PP_H4 > 0.8], we identified strong evidence of four additional retinol-associated genes outside of genome-wide significant loci (at least +/-1 megabase away). These were as follows: MLXIPL, which binds to carbohydrate response elements to regulate triglycerides^37-39; GSK3B, a gene that encodes a member of the glycogen synthase kinase family involved in metabolism and glycaemic homeostasis⁴⁰; the tankyrase gene (TNKS) implicated in processes like Wnt signalling⁴¹; and INHBC, part of the inhibin family of proteins with important endocrine functionality⁴². Genetically predicted mRNA expression of MLXIPL in adipose, pancreas, and breast mammary tissue was inversely associated with circulating retinol levels. Conversely, TWAS analyses revealed that genetically predicted expression of GSK3B and TNKS was positively associated with circulating retinol levels. Finally, genetically predicted plasma protein expression of INHBC showed a positive correlation with retinol levels (Z_PWAS = 4.72). We then used a more conservative approach whereby finemapped variants that influence protein expression (pQTLs) were used as instrumental variables (IV) to estimate the causal effect of plasma proteins on retinol using Mendelian randomisation (MR, Online Methods, Supplementary Table 14). We applied the same filters to the results (FDR < 0.05 and PP_H4 > 0.8). These analyses further supported that upregulated INHBC likely increases serum retinol, with each SD increase in plasma protein expression associated with a small but highly statistically significant [0.05 per SD in expression, 95% CI: 0.03, 0.07] impact on circulating retinol. In line with expectation, pQTL-MR, and subsequent colocalisation, further genetically validates that elevated RBP4 protein abundance correspondingly increases serum retinol with somewhat large effect (0.6 [95% CI: 0.48, 0.72] SD in retinol per SD in plasma RBP4 expression). There was also evidence that RBP4 protein expression and retinol colocalise under the hypothesis of a single causal variant (PP_H4 = 1). Considering the eight genes prioritised in this and the previous section (RBP4, GCKR, FOXP2, TTR, MLXIPL, GSK3B, TNKS, INHBC), we found that these genes exhibited upregulated expression in the liver (P_Adjusted < 0.05) upon analysing data from 54 GTEx tissues. This further consolidates the salience of hepatic processing to genetic influences on circulating retinol abundance. Pathway analyses of these genes demonstrated that they were enriched amongst factors involved in carbohydrate metabolism (Supplementary Table 15).

Wide-ranging evidence of causal effects of retinol across the human clinical phenome

The role of retinol in human health and disease has been of long-standing interest. However, most evidence has been observational, limiting the ability for causal inference. Moreover, randomised controlled trials (RCT) of interventions like retinol supplementation and synthetic retinoids have only been performed for a small fraction of the traits implicated through observational studies. We sought to increase our understanding of causal effects of retinol on human health by leveraging genetic variants associated with retinol uncovered in this study as IVs. Given certain assumptions are met, these genetic proxies of retinol can be utilised to estimate causal effects of circulating retinol at scale using Mendelian randomisation (MR) (Online Methods). Firstly, we utilised a single IV in RBP4 (rs10882283), as this gene has a clear and unambiguous association with circulating retinol levels, and therefore, is less likely to be prone to horizontal pleiotropy than other retinol-associated loci. We do caution, however, that RBP4 does exert some other functionality that may not be directly related to retinol transport⁴³, although this gene is still likely the best available single IV associated with circulating retinol at genome-wide significance. We utilised this RBP4 IV (METSIM+INTERVAL meta-analysis effect size) to estimate the effect of retinol on over 19,500 outcomes in the IEUGWAS database (IEUGWASdb, Online Methods, Supplementary Table 16). After multiple-testing correction (FDR < 0.01), retinol was found to putatively exert causal effects on outcomes including several lipid traits, leukocyte counts (total leukocytes and neutrophils), reticulocytes, optic disc area, and resting-state connectivity of a functional MRI (fMRI) derived network edge. For instance, each SD in circulating retinol was associated with a -0.1 SD [95% CI: -0.14, -0.06] decrease in leukocyte count, whilst this unit retinol increase was estimated to increase optic disc area by 0.22 SD [95% CI: 0.12, 0.32]. We further interrogated a representative subset of these associations using colocalisation to test if the signals are driven by a shared causal variant in RBP4 (Supplementary Table 17). Colocalisation strongly supported that the effect of circulating retinol on leukocyte count, optic disc area and the edge in the fMRI network arises from a shared causal variant in RBP4. In contrast, the effect of retinol on lipids through RBP4 was shown to likely arise from linkage due to the proximal gene FFAR4 that encodes a free-fatty acid receptor. Therefore, the lipid findings most likely do not represent a causal impact of circulating retinol, at least through RBP4, but rather the influence of FFAR4.

To boost power for discovery of causal retinol effects, we then utilised all independent (LD r² < 0.001) genome-wide significant SNPs as IVs (Online Methods). Firstly, we used the inverse-variance weighted estimator with multiplicative random effects (IVW-MRE) to estimate the effect of retinol on IEUGWASdb outcomes for which at least six of the IVs were available (> 17,000 outcome phenotypes). There was moderate positive correlation between the IVW-MRE and single RBP4 estimates across all traits tested (r = 0.43, Supplementary Figure 4). Whilst well powered, the IVW-MRE assumes all IVs are valid, which is unlikely in practice. Therefore, we developed a pipeline to prioritise the most confident causal relationships that survive multiple-testing correction considering the IVW-MRE estimates (FDR < 0.01, Online Methods, Figure 2A, Supplementary Tables 18-20). This was comprised of three tiers, with Tier #1 being the highest level of evidence. Retinol/outcome pairings that were assigned a tier had to exhibit no significant heterogeneity between IV-outcome effects, a non-significant MR-Egger intercept (which screens for unbalanced pleiotropy), and not be driven by a single IV. Four additional MR methods with differing assumptions were then applied in this study (Online Methods). From the trait pairings that passed the above heterogeneity and pleiotropy filters, Tier #1 traits were those for which all five methods were statistically significant, whilst Tier #2 traits had 4/5 significant methods, and Tier #3 3/5 methods significant. There were no Tier #1 retinol/outcome trait pairings, but several Tier #2 and Tier #3 trait pairings (Figure 2B-C), with all of them directionally consistent with the estimates from the single RBP4 IV, supporting their validity. Broadly, we found that retinol increased body fat related measures, resting-state fMRI connectivity of several network edges, as well as food consumption phenotypes related to carbohydrates. Retinol also exhibited evidence of a relationship with the cortical thickness and surface area of several brain regions, as well as microbiome composition and keratometry measurements. These results were dominated by continuous traits, for which we are better powered, however, there were some binary traits assigned Tier #2 or Tier #3 evidence. Specifically, retinol was associated with decreased odds of tuberculosis sequalae and coxarthrosis (arthrosis of the hip), whilst it was associated with increased odds of non-specific skin eruptions, adverse asthma/COPD medication effects, and dental problems. We caution that all these inferred associations require further investigation, and should be treated with requisite caution, as reviewed elsewhere^{19, 44, 45}. A consideration of this approach that leverages multiplicative random effects with relatively few IVs (< 10), is the potential influence of residual standard errors below 1 on the estimation of the IVW standard errors. We see some evidence of this impact on the IVW-MRE estimates for Tier #2 and Tier #3 traits – as these traits exhibit no significant heterogeneity between IV estimates. Specifically, whilst all fixed effect IVW results are still highly statistically significant for these traits, they have larger standard errors than the IVW-MRE, indicative of residual standard errors < 1. This is a function of the MRE not scaling the standard error of the IVW by the model’s residual standard error like in the fixed effects model. We discuss this further in the supplementary text and in supplementary figure 5. However, these issues only impact the P-value of the IVW-MRE relative to that of the fixed effects, and Tier #2/Tier #3 traits still exhibit non-zero evidence across multiple methods and no indication of a single IV driving the association. Moreover, using instead a fixed effects IVW estimator as the primary test for trait pairings with no heterogeneity (Cochran’s Q P < 0.05) yields similar outcomes being prioritised as most statistically significant after FDR correction (Supplementary Text). It also important to consider when interpreting these estimates that MR approaches are only interpretable under the assumptions they make, and as a result, any potential causal relationship reported here requires further validation, ideally using a randomised control trial design.

• Download figure
• Open in new tab

Figure 2. Estimated causal effects of circulating retinol across the human clinical phenome.

(a) Prioritisation pipeline overview for retinol causal estimates [inverse-variance weighted estimator with multiplicative random effects (IVW-MRE)] that survive multiple-testing correction (FDR < 0.01). These estimates are subjected to tests for heterogeneity and pleiotropy (Online Methods), with a tier then assigned based on how many of the five Mendelian randomisation (MR) methods applied are at least nominally statistically significant. In panels (b) and (c), the left-hand plot denotes the Z-score (beta/SE) from the MR IVW-MRE estimates. Positive Z scores denote a positive IVW-MRE estimate of the effect of circulating retinol on that trait, and vice vera. The traits are coloured by their broad phenotypic category. The right-hand plot visualises the Z score using the IVW-MRE model verses that of the MR estimate using the RBP4 IV alone (Wald Ratio). The dotted lines approximately represent nominal statistical significance (P < 0.05). In panel (b), just tier #2 traits are plotted (Online Methods), whilst panel (c) plots both tier #2 and tier #3 traits.

We hypothesised that the putative effect of retinol on body fat could be one explanation for its relationship in this study to brain phenotypes beyond a direct effect of retinoid signalling, particularly given that obesity and adiposity have been linked with MRI related indices⁴⁶. However, using IVs for body fat percentage, we did not find any strong evidence that it is causally related to any of the retinol-associated brain regions (Supplementary Table 21), suggesting a direct effect of retinol on these regions/networks or a relationship induced through some other unobserved confounder. Reverse causality for these Tier #2 and Tier #3 exposures was then also considered, although for binary traits this should be treated purely as a test of the null hypothesis given the difficulties in using binary traits as IVs⁴⁷. There was no strong evidence for reverse causality of any of these traits (Supplementary Table 22). One exception to this was in relation to expression of the protein PEAR1, for which there was very nominal evidence of bidirectional effects.

A limitation of the above phenome-wide analyses is that the multiple-testing burden that arises from the inclusion of over 17,000 traits may obscure retinol effects on binary disease phenotypes, as these are usually less powered than GWAS of continuous traits. To overcome this, we also applied the above pipeline using all retinol IVs to 1141 binary endpoints with at least 1000 cases in FinnGen release 8 (not featured in IEUGWASdb), allowing a phenome-wide analysis of electronic health record derived binary outcomes. There were eight disease phenotypes that retinol was associated with after multiple testing correction (Figure 3A, FDR < 0.01), which increased to 19 with an exploratory FDR < 0.1 threshold (Supplementary Table 23). After applying the above pipeline to these results that considers heterogeneity, pleiotropy, and consistency across MR methods, there were four disease endpoints with Tier #3 evidence (Figure 3B, Supplementary Table 24). Specifically, retinol was estimated to increase the odds of congenital malformations of the heart and great arteries, whilst it was protective for type 2 diabetes with coma and inflammatory liver disease. As these traits had tier #3 evidence, there was some inconsistency in the strength of the results across different MR methods, and therefore, these relationships should be interpreted cautiously. One of the most active areas of research in retinol epidemiology is the relationship between retinol and cancer risk⁴⁸. The estimated effect of circulating retinol on the odds of any malignant neoplasm was not significantly different than one - OR = 0.97 [95% CI: 0.91, 1.04], P = 0.423 (Supplementary Figure 6). However, there was some indication of a protective effect of retinol on squamous non-small cell lung cancer, which approached the threshold for statistical significance after FDR correction - OR = 0.64 [95% CI: 0.51, 0.80], P = 8.19 × 10^-5, q = 0.01. There was also some data to support retinol having effects on other respiratory neoplasm endpoints (Supplementary Figure 6). Given previous observational evidence that retinol is protective for lung cancer⁴⁸, as well as some data supporting the use of synthetic retinoids like bexarotene in lung neoplasms⁴⁹, this relationship warrants further exploration.

• Download figure
• Open in new tab

Figure 3. Genetic estimates of the relationship between circulating retinol and binary disease endpoints in the Finnish population.

(a) Volcano plot denoting the phenome-wide MR estimates of retinol (SD units) and the log odds of each binary endpoint in FinnGen release 8. The x-axis denotes the log odds (IVW-MRE), whilst the y-axis denotes the -log₁₀ P-value of the MR estimate. The dotted horizontal line is approximately equivalent to the magnitude of P-value that is estimated at an FDR of 1%. Categories highlighted by a blue box on the legend indicate that an endpoint in that category survived multiple-testing correction. (b) Retinol effects on binary FinnGen endpoints after multiple-testing correction (FDR < 0.01) with the strongest evidence from all phenotypes tested based on heterogeneity, pleiotropy, and consistency across MR estimators (Tier #3, Online Methods). Each panel denotes the effect sizes (odds ratio with 95% confidence intervals per SD in circulating retinol) for the IVW-MRE, the IVW with fixed effects (IVW-FE), the MR-Egger, Weighted median, and Weighted Mode, respectively.

We also investigated evidence for bidirectional effects involving these Tier #3 traits that retinol is genetically predicted to influence. As described above, using binary traits as exposures in MR should be treated with suitable caution and are often underpowered. The use of either genome-wide significant or suggestively significant SNPs (P < 1 × 10^-5) as IVs did not indicate evidence of reverse causality of these diseases to retinol (Supplementary Table 25). However, given some of the statistical limitations of these analyses, such effects warrant further consideration.

Genetic evidence that lipids and kidney function influence circulating retinol

It is also clinically valuable to understand exposures and diseases that impact circulating retinol abundance. To explore this in greater detail, we leveraged retinol as an outcome trait in MR analyses. We utilised a diverse range of thousands of continuous and ordinal phenotypes from IEUGWASdb as exposures in a similar pipeline described above (Online Methods). Several lipid species were demonstrated to putatively influence retinol abundance after multiple-testing correction (FDR < 0.01, Figure 4A, Supplementary Table 26); for example, triglycerides were implicated to increase circulating retinol whilst cholesteryl ester related traits decreased circulating retinol. We also saw a positive effect of the frequency of solarium and sun lamp use on retinol, which may arise from behavioural related mechanisms. Furthermore, our findings suggest that increased retinol levels are associated with a susceptibility-weighted MRI measure in the left putamen, called T2*. T2*, reflecting magnetic susceptibility relative to tissue water, can be influenced by factors such as iron and calcium content. This association may also be influenced by behavioural or other pathways that warrant further investigation.

• Download figure
• Open in new tab

Figure 4. Exploring the causal effects of continuous exposures on circulating retinol.

(a) Exposure traits that demonstrated a significant causal estimate (IVW-MRE) on circulating retinol after multiple-testing correction (FDR < 0.01). Traits are coloured relative to their broad phenotypic category. (b) Multivariable MR (MVMR) models investigating the effect of creatinine and major lipid species on circulating retinol. Each panel represents the results from a different MVMR model (each with different underlying assumptions (Online Methods)). The exposure – retinol relationship plotted is conditional on the three other traits in the model.

After applying the same tiering system used for retinol as an exposure, we observed Tier #1 evidence strongly supporting a causal effect of creatinine on circulating retinol. This relationship is biologically plausible and likely represents an association with kidney function^{50, 51}. Due to the biological complexity of lipid traits, we observed significant heterogeneity between IV effects, and as a result, they were not assigned a tier in our analysis. However, the effect of triglycerides on increasing circulating retinol levels is consistent with established knowledge of retinol biology, as well as this study implicating two genes that are mechanistically confirmed to impact triglycerides (GCKR and MLXIPL). We then leveraged the CAUSE model to distinguish causal effects of creatinine on circulating retinol from correlated pleiotropy that may arise between these two traits due to the extensive polygenicity of creatinine⁵². We found that a model that includes a causal effect of creatinine on retinol was more parsimonious than a model of pleiotropy alone (‘sharing model’) through comparison of these models using the Bayesian expected log pointwise posterior density (ELPD) method (Supplementary Figure 7, ΔELPD_{Sharing vs Causal} = -4.34, P = 8.9 × 10^-3). Given that IVs for creatinine could plausibly act through lipid species like triglycerides to influence circulating retinol, we then constructed multivariable MR (MVMR) models that estimated the creatinine to retinol relationship conditioned on high density lipoprotein (HDL), low density lipoprotein (LDL), and triglycerides (Online Methods). While there was some evidence that the effect of creatinine on retinol could arise due to triglycerides, there was also evidence to suggest an independent effect of both triglycerides and creatinine on increasing circulating retinol, depending on the modelling parameters used (Figure 4B, Supplementary Text). In addition to investigating the causal effects of exposures on retinol, we also examined the possibility of bidirectional effects. We found very weak evidence that retinol has a negative effect on creatinine (P = 0.023), but there was no significant evidence to suggest bidirectional relationships between retinol and the other implicated exposures.

Finally, we explored pharmacological agents and molecular perturbagens that may influence circulating retinol (Online Methods). Considering the novel genes prioritised in this study with an assigned direction of expression (TWAS/PWAS, pQTL MR), it was found that GSK3B is a drug-target known to be inhibited by lithium and related compounds. This may be of clinical interest as it suggests that lithium, utilised as a therapy in mood disorders, may decrease circulating retinol via its inhibition of GSK3B given that genetically predicted expression of this gene was positively associated with retinol. We then employed computational signature mapping to further characterise pharmacological agents related to retinol (Online Methods). However, these analyses did not yield any compounds for which the in vitro transcriptomic signature significantly matched or opposed genetically predicted expression associated with retinol after multiple-testing correction (Supplementary Table 27). We then considered perturbagen signatures aggregated to biological pathways or overall mechanisms of action (MOA) groups of compounds (Online Methods). After multiple-testing correction (FDR < 0.05), there were 13 gene-set based perturbagen signatures that were significantly similar to the directionality of genetically predicted expression associated with retinol (Supplementary Table 28). For example, the expression signature of compounds in the HDAC inhibitor MOA opposed expression genetically predicted to increase serum retinol. This can be interpreted as while no single HDAC inhibitor was significantly associated with retinol, there was at least some evidence for the overall relationship with this MOA. This accords with the suggested effect of HDAC inhibitors like valproic acid on downregulating expression of RBP4^53–55, a gene not included in the signature mapping analyses given the large effect of its encoded protein on retinol.

Genetically proxied retinol can identify individuals outside of the normative range of circulating retinol for a given age

We were also interested in evaluating the performance of a genetically proxied index of circulating retinol, that is, a circulating retinol polygenic score (PGS). The independent TwinsUK cohort was utilised to tune and evaluate retinol PGS (Online Methods). We used several methods to evaluate the performance of a retinol PGS in this cohort. Firstly, we test different retinol PGS configurations in a random selected training subset of the model (70% of cohort), using a linear mixed model to account for relatedness between the twin pairs. The best performing retinol PGS configuration in the training subset explained approximately 2.12% of the phenotypic variance of retinol when applied to the remaining 30% of the cohort (mean variance explained across three retinol measurement timepoints). We also found similar performance when applied to each subset of twin pairs (Supplementary Table 29). A limitation of this approach for using the same cohort for tuning and testing the retinol PGS is that the estimated effect sizes may not be representative. As a result, we employed an approach to tune weights for the PGS using the summary statistics alone. This was achieved through leveraging the principles of probabilistic finemapping to update variant weights by their posterior probability of association (Online Methods). Due to the modest polygenicity of retinol, this method upweights a small number of variants. Despite this, these scores were still significantly associated with circulating retinol in TwinsUK (Supplementary Table 29).

Like most micronutrients, circulating retinol has been shown previously to have a complex relationship with age⁵⁶. However, population-level approaches investigating these effects do not account for inter-individual variability. We hypothesised that normative modelling could be used to characterise individual patterns of circulating retinol, and to evaluate the contribution of genetics to these individualised effects. Normative modelling, derived from the application of growth charts in paediatric medicine, aims to estimate normative reference ranges of variation in the population (e.g., of circulating retinol) based on age and/or other relevant variables. Here, we established reference ranges for circulating retinol as a function of age using generalised-additive models for location, scale, and shape (GAMLSS) frameworks in the TwinsUK dataset (Online Methods, Supplementary Text). These models were constructed in the full cohort, as well as in one twin subset only for comparison. Briefly, this involved identifying the optimal distribution for the GAMLSS model using all samples, followed by splitting the data into two partitions. One partition was utilised to estimate centiles (5^th, 25^th, 50^th, 75^th and 95^th), and the left-out subjects were benchmarked against this reference chart to determine the position of their individual retinol measurement (Online Methods). This process was then repeated using the opposite subject allocation. An individual’s retinol level was classified as infra-normal if their measured retinol fell below the 5^th percentile for their age, and as supra-normal if their retinol value exceeded the 95^th percentile for their age.

We then investigated the extent to which retinol PGS was associated with these individual profiles of circulating retinol with respect to age (Supplementary Table 30). By way of example, we report results forthwith from the more conservative modelling approach which only used one half of the twins. Retinol PGS was at least nominally significantly associated with supra and infra-normal deviations except for supra-normal deviations at the first (youngest visit) for which there was only a trend observed (Supplementary Figures 8-9). For example, at the second visit each SD in retinol PGS was associated with an approximately 70% [95% CI: 23%, 136%] increase in the odds of exhibiting supra-normal retinol levels for a given age relative to all remaining participants, whilst conversely reducing the odds of displaying infra-normal levels by approximately 37% [95% CI: 12%, 55%]. These data suggest that genetics is a non-zero contributor to circulating retinol levels that fall outside the normative range for a given age. In future, a normative modelling approach could be utilized to examine additional factors, including dietary intake, as well as the interplay between genetics and other influences, that might contribute to individualized deviations in retinol levels relative to the population benchmarks.

Study cohorts

The proceeding section outlines the datasets included in the genome-wide meta-analysis of circulating retinol, as well as the replication cohort.

Genome-wide meta-analyses

We conducted genome-wide meta-analysis of common variants (MAF ≥ 0.01) using METAL (version March 2011), followed by a rare variant (MAF < 0.01) meta-analysis also with METAL. The METSIM and INTERVAL cohorts were integrated for the primary meta-analysis as they both had expansive genome-wide coverage of common and rare variants. A sample size weighted meta-analysis of Z scores (Stouffer’s method) was utilised for this purpose. We also conducted the METSIM+INTERVAL meta-analysis via an inverse-variance weighted estimator with fixed effects to estimate the effect sizes of effect alleles in SD units of plasma retinol given this was the unit of both the METSIM and INTERVAL GWAS, as well as both studies using the same Metabolon Inc. platform for metabolite quantification. Heterogeneity between the studies was assessed using Cochran’s Q test. We then conducted a meta-analysis with fewer available variants that also included ATBC+PLCO using Stouffer’s method (as the unit for this GWAS differed from METSIM and INTERVAL). This larger sample-size meta-analysis was also restricted to common variants as there was very limited rare variant coverage in ATBC+PLCO.

The FUMA v1.4.1 (Functional Mapping and Annotation of Genome-Wide Association Studies) platform was utilised to annotate variants, define lead SNPs, and infer loci boundaries for genome-wide significant signals (P < 5 × 10^-8)⁸⁴. We utilised the default settings for defining independent significant SNPs (r² ≤ 0.6) and lead SNPs (r² ≤ 0.1). LD estimation was achieved using the 1000 genomes phase 3 European reference panel, with LD blocks within 250 kb of each other merged into a single locus. We then attempted to replicate the lead SNPs from the eight genome-wide significant loci (METSIM+INTERVAL) in TwinsUK, as described in the previous section. Given the small sample size of TwinsUK, we sought to ascertain if the lead SNPs were directionally consistent with the meta-analysis. This was achieved by taking the mean SNP Z score across the three timepoints for the two sub-cohorts of unrelated participants, with a binomial test utilised to infer whether the number of lead SNPs (mean Z) that were directionally consistent was greater than chance alone (Binomial P < 0.05).

In the rare-variant meta-analysis, we annotated variants using the Functional Annotation of Variants (FAVOR) online resource⁸⁵. Phenome-wide association profiles of selected variants were also investigated using the pheweb browser collated from FinnGen release 8 (https://r8.finngen.fi/)86. Rare variants were then aggregated to genes through leveraging the characteristics of the Cauchy distribution^{87, 88}. In this approach, gene-wise P values are summed and then transformed to approximate a Cauchy distribution, which due to its heavy tail is insensitive to correlations amongst the P values. This behaviour of the Cauchy distribution is important as covariance amongst rare variants is difficult to estimate, and therefore, this approach guards against inflated type I error due to potential unknown covariance/LD between rare variants. Code for implementing the Cauchy aggregation was adapted from https://github.com/yaowuliu/ACAT.

SNP heritability estimation

Summary statistics for the meta-analyses were ‘munged’ using the munge_sumstats.py script from the ldsc repository of scripts (https://github.com/bulik/ldsc) and only common (MAF > 0.05) HapMap3 variants outside of the MHC retained. SNP heritability was then estimated using the LDSR model and the 1000 genomes phase 3 reference panel⁸⁹. By way of comparison, we then estimated SNP heritability using the LDAK model via the SumHer implemented in LDAK v.5.2 (https://dougspeed.com/)90. Pre-computed tagging files derived from 2000 white British individuals in the UKBB for HapMap3 SNPs were utilised to calculate SNP heritability using the LDAK-thin and BLD-LDAK models, as described elsewhere^{90, 91}. We also estimated partitioned SNP heritability via LDSR using a multi-tissue and cell-type panel⁹².

Empirical Bayes’ modelling of the genetic architecture of retinol

We investigated the polygenicity of the genetic architecture of circulating retinol using an Empirical Bayes’ adaptive shrinkage method termed ashR⁹³. Functions to perform this method were implemented via the ashR R package v2.2-54 (https://github.com/stephens999/ashr). Briefly, this approach models effect sizes, along with their standard error, as a mixture of zero and non-zero effects. Empirical Bayes’ inference is performed under the assumption that the distribution of these variant effects is unimodal, which is a realistic assumption for genetic effects on complex traits. The METSIM+INTERVAL IVW meta-analysis was utilised for this as it provides an interpretable effect size and standard error for each variant (plasma SD units). In line with previous work⁹⁴, we annotated each HapMap3 variant with its corresponding LD score from the 1000 genomes phase 3 European reference panel, as well LD scores from the UKBB White Great British samples for comparison, and sorted these into bins of similar LD scores (N_Bins=1000 and N_Bins=5000). The ashR Empirical Bayes’ inference of the proportion of non-zero effects was undertaken in each LD score bin, followed by calculating the mean across all bins. We utilised a generalised additive model to plot a smoothed trend line of the relationship between increasing LD score bin (higher LD score) and the proportion of non-zero effects.

Gene prioritisation

Gene prioritisation was performed within genome-wide significant loci, as well as outside of loci that obtained genome-wide significance. Due to the better coverage of common variants, we focused on the METSIM+INTERVAL meta-analysis for this analysis. The pipeline for prioritising putative causal genes within the eight genome-wide significant loci was adapted from a previous GWAS performed by our group⁹⁵. The following criteria were utilised in this study: 1) closest transcription start site (TSS) to the lead SNP, 2) closest gene (any) to the lead SNP, 3) gene encoding retinoid transporter or enzyme in locus, 4) non-synonymous variant in locus, 5) the most statistically significant GTEx eGene [expression quantitative trait loci (eQTL) signal] in locus, 6) most significant GTEx eGene in finemapped credible set for eQTL signal [posterior inclusion probability (PIP) > 0.1, DAP-G method]^{96, 97}, 7) strongest plasma pGene [protein quantitative trait loci (pQTL) signal] drawn from finemapped pQTLs (PIP > 0.5)⁹⁸, 8) highest CADD score⁹⁹, 9) lowest RegulomeDB score¹⁰⁰, 10) the OpenTargets V2G predicted gene for the lead SNP¹⁰¹, and 11) genes physically mapped to SNPs in the 95% credible set derived from probabilistic finemapping (assuming a single causal variant such that LD did not have to be modelled)¹⁰². We used a prior variance of 0.15 to approximate Bayes’ factors from variant-wise effect sizes (METSIM+INTERVAL IVW meta-analysis), in line with previous work finemapping association signals with quantitative traits¹⁰³. We characterised which genes satisfied the greatest number of the above criteria on a per locus basis.

FOXP2 was one of our confidently prioritised genes but its biological significance outside of the brain is less well understood. To investigate this, we analysed RNA sequencing data from an in vitro experiment that overexpressed FOXP2 in a human osteosarcoma epithelial cell line (U2OS) via transfection of wild-type FOXP2 expressing plasmids. Raw read counts from five control cell line replicates versus five plasmid transfected FOXP2 overexpression replicates were downloaded from the Gene Expression Omnibus (GEO) resource (GEO Accession: GSE138938). Data normalisation, filtration, and differential expression analyses were performed using the edgeR package version 3.34.0¹⁰⁴. Specifically, raw counts were firstly normalised to library size and lowly expressed genes with fewer than 10 raw counts in the smallest library removed via a counts-per-million thresholding approach. Data were inspected before and after the filtration step via coefficient of variation (BCV) and multidimensional scaling (MDS) plots (Supplementary Figure 10). Differential expression for each gene that survived quality control was then performed using exact tests for differences in the means between two groups of negative-binomially distributed counts. We defined a differentially expressed gene as those which survived multiple-testing correction using the Bonferroni method (P_Corrected < 0.05), with three different absolute log₂ fold change (FC) cut-offs considered: |log₂FC| > 1.5, |log₂FC| > 2, and |log₂FC| > 5. The use of Bonferroni correction and large absolute FC thresholds is very conservative; however, given the volume of differentially expressed genes, and a relatively large number of replicates for a cell line experiment boosting power, we believe these strict parameters are warranted to priortise the most salient FOXP2 associated signals. The overrepresentation of each set of candidate genes amongst biological pathways and other ontology sets was tested using g:Profiler¹⁰⁵.

To identify potential causal genes that have not reached genome-wide significance at our current sample size, we integrated the circulating retinol GWAS with genetic effects on mRNA and protein expression. Firstly, we conducted a transcriptome and proteome-wide association study (TWAS/PWAS) of circulating retinol using the FUSION approach¹⁰⁶. As outlined previously^107–110, FUSION leverages models of cis-acting genetically regulated expression (GReX) that exhibit statistically significant non-zero heritability. Variant weights from GReX models are integrated with the effect of those same variants on retinol to estimate the direction of genetically regulated expression associated with increasing circulating retinol. TWAS GReX (mRNA) were estimated previously using GTEx v8 (http://gusevlab.org/projects/fusion/). We selected the following biologically informative tissues to perform TWAS based on known retinol biology or tissues that exhibited at least nominally significant (P < 0.01) enrichment of SNP heritability in the partitioned-LDSR model. The selected tissues were: small intestine terminal ileum, pancreas, liver, adipose (visceral omentum), adipose (subcutaneous), breast (mammary tissue), and whole blood. It has been suggested previously that TWAS signals from tissues that are less directly trait relevant can induce spurious associations, which is why we limited our hypothesis space to these tissues¹¹¹. Protein GReX were derived from plasma (ARIC study), as outlined elsewhere⁹⁸. We applied Benjamini-Hochberg false discovery rate (FDR) correction across all TWAS Z, followed by all PWAS Z. Colocalisation between GReX models and retinol was performed for all TWAS/PWAS signals that were at least nominally significant via the coloc package as implemented by FUSION (single shared variant hypothesis)¹⁰³. We then considered a more conservative approach to priortise proteins for whom expression could be causally linked to circulating retinol through leveraging finemapped plasma pQTLs (PIP > 0.5, ARIC) as instrumental variables (IV) for Mendelian randomisation^{19, 44, 98}. The Wald ratio method was implemented for proteins with single IVs, whilst an inverse-variance weighted estimator with fixed effects was utilised for proteins with more than 1 IV. Fixed effects were used for the IVW rather than multiplicative-random effects in this instance as no protein had > 4 IVs. The use of IV based approaches for identifying trait-associated genes versus GReX has been discussed extensively elsewhere¹⁰⁹. Colocalisation was also performed for proteins that survived FDR correction. Mendelian randomisation was performed using the TwoSampleMR package v0.5.6, with IVs clumped through leveraging LD from the 1000 genomes phase 3 European reference panel to retain only independent pQTLs (r² < 0.001).

Finally, we investigated the tissue specificity of prioritised genes from GWAS loci and the TWAS/PWAS/MR approach using FUMA. To do this, we compared the expression of these prioritised genes using a t-test (one-sided and two-sided) against all other available genes in GTEx v8 on a per tissue basis (54 tissues), followed by applying Bonferroni correction to these P-values⁸⁴. We also conducted pathway analyses of these genes using g:Profiler with default parameters¹⁰⁵.

Causal inference

We developed and implemented a comprehensive pipeline to leverage this retinol GWAS to identify putative causal effects of retinol on traits across the human clinical phenome, as well as traits that causally influence retinol in the reverse direction. This was achieved using Mendelian randomisation (MR), which has been reviewed extensively elsewhere^{19, 45, 69}. Firstly, we considered circulating retinol as the MR exposure. The lead SNP in RBP4 was first chosen as a single IV to proxy circulating retinol as out of all the genes implicated in genome-wide significant loci, RBP4 has exhibits the most specificity in terms of its relationship with serum retinol. We utilised the Wald ratio method to estimate the effect of the circulating retinol increasing rs10882283-A allele (METSIM+INTERVAL IVW effect size) on over 19,000 outcomes in the IEUGWASdb v6.9.2 resource via the ieugwasr package version 0.1.5²⁶. To correct for multiple testing, we applied the false discovery rate (FDR) method, and we retained only those retinol/outcome pairs with estimates that were significant below the 1% FDR threshold (q < 0.01). Subsequently, we investigated whether significant trait pairs colocalised using coloc v5.1.0 with default priors.

Although the single IV approach is more conservative, power to detect causal effects can be boosted by using all available independent (r² < 0.001, 1000 genomes phase 3 European reference panel) genome-wide significant SNPs as IVs. The inverse-variance weighted estimator with multiplicative random effects (IVW-MRE) was utilised to estimate the effect of retinol on outcomes from IEUGWASdb for which at least 6 IVs were available in that GWAS¹¹². The IVW approach has a zero percent breakdown level as it assumes all IVs are valid for use in Mendelian randomisation, which is often unrealistic in practice. We developed pipeline to prioritise the most reliable causal estimates from the IVW-MRE that survived multiple-testing correction (q < 0.01). This involved identifying retinol MR estimates on traits for which the following applied: i) no significant heterogeneity (P < 0.05) between IV exposure-outcome effects tested with Cochran’s Q¹¹³, ii) a non-significant intercept of an MR-Egger model that does not constrain the intercept to pass through the origin¹¹⁴, and iii) no evidence that leaving out any single IV ablates the statistical significance (at least P < 0.05) of the estimate¹¹⁵. Traits satisfying these criteria were then assigned a tier (Tier #1, Tier #2, Tier #3) based on the statistical significance of the retinol causal estimate using other MR methods besides the IVW-MRE that have different assumptions regarding IV validity. These were: the IVW with fixed effects (does not model heterogeneity like with MRE), the weighted median method¹¹⁶, the weighted mode method¹¹⁷, and the MR-Egger method¹¹⁴. The underlying assumptions and methodological considerations of using these methods have discussed extensively elsewhere^{118, 119}. Traits for which the effect of retinol was at least nominally statistically significant using all five methods were assigned as Tier #1, whilst four tests being statistically significant was Tier #2, and three tests being statistically significant was Tier #3. The F-statistic and I² of the IVs was also assessed to ensure they were well powered (F > 10) and suited for MR-Egger (I² > 0.9), respectively¹²⁰. To investigate the effect of body fat percentage on the retinol associated MRI indices, we used IVs from a non-overlapping GWAS of that trait (N=65,831)¹²¹. We then repeated the above process of binary outcomes with at least 1000 cases from FinnGen release 8, which is not included in the current version of IEUGWASdb at time of analysis, to increase power to detect effects on disease endpoints.

We also systematically investigated continuous outcomes that may causally impact circulating retinol. Continuous outcomes were our focus as causal estimates from binary exposures are difficult to interpret and are often less powered in the context of MR⁴⁷. Exposures were filtered from phenotypes available in IEUGWASdb to retain continuous traits, with further filtering to identify traits with ≥ 5 genome-wide significant, independent (r² < 0.001, 1000 genomes phase 3 European reference panel) IVs available in the retinol GWAS. The same pipeline as above was then applied (IVW-MRE FDR < 0.01, followed by sensitivity analyses and tier assignment). The causal estimate of creatinine on retinol was a Tier #1 trait, and due to the pervasive polygenicity of creatinine afforded by its large sample size, we followed up this relationship using the MR model “Causal Analysis Using Summary Effect Estimates” (CAUSE), as described elsewhere using the CAUSE R package v1.2.0⁵². Briefly, this method is a more polygenic approach and seeks to distinguish casual effects from correlated pleiotropy by fitting competing models that account for these terms and comparing them using ELPD.

After using 1 million random variants to estimate nuisance parameters, LD clumping and thresholding was applied to the serum creatinine summary statistics (IEUGWASdb trait ID: met-d-creatinine) in line with the original CAUSE publication (P < 0.001, r² < 0.01, 1000 genomes phase 3 European reference panel). The competing CAUSE models were then fit (null, sharing, and causal), ensuring that all Pareto k estimates were < 0.5 during the model comparison using ELPD.

We then performed a multivariable MR (MVMR) analysis to estimate causal effects of creatinine on retinol conditioned on three major lipid species using GWAS from the global lipids genetics consortium (LDL, HDL, and triglycerides)¹²². MVMR was undertaken in accordance with previous work using the R packages MVMR v0.3 and MendelianRandomization v0.6.0¹¹⁸. Briefly, this entailed identifying variants associated at genome-wide significance with at least one of the four exposures that are independent (r² < 0.001), calculating a conditional F-statistic for multivariable instruments¹²³, and applying four MVMR models (IVW, Egger regression, Weighted Median, and a LASSO based penalised regression approach for selecting the optimal IV configuration)¹²⁴.

Drugs and perturbagens associated with circulating retinol

We also considered drugs that may influence circulating retinol. We searched genes prioritised from our pipeline with an assigned direction of retinol-associated expression using DGIdb v4.2.0 and DrugBank v5 to identify retinol associated genes targeted by drugs^{125, 126}, outside of known drugs that target RBP4 and TTR. Retained drug-gene interactions were restricted to those with known mechanism of action and > 2 lines of supporting evidence. We also utilised computational signature mapping to identify pharmacological agents that may enhance or inhibit the expression of genes associated with circulating retinol. To boost power for signature mapping, we considered all gene that were nominally associated with retinol in the TWAS/PWAS (P < 0.05) that exhibited moderate colocalisation of a shared causal variant (PP_H4 > 0.4). These genes were uploaded to the Connectivity Map Query online tool to quantify the similarity, termed connectivity, with drug perturbagen associated expression profiles, as outlined in detail elsewhere¹²⁷.

Polygenic scoring

We used the independent TwinsUK replication cohort, described in a preceding section, to investigate retinol polygenic scores (PGS). PGS were applied to genotyped variants and high confidence imputed variants (R² > 0.8) in TwinsUK. There were two different methodologies implemented to construct PGS: LD clumping and thresholding (LD C+T) and a probabilistic finemapping based method that scales variant effect sizes based on their posterior probability of causality, thereby upweighting signals more likely to be causal (RápidoPGS)^{128, 129}. In the LD C+T approach, variants were clumped using the within sample LD of TwinsUK at the following P-value thresholds: 5×10^-8, 1×10^-5, 1×10^-3, 0.01, 0.05, 0.1, 0.5, and 1. Additive PGS were then profiled using PRSice2 v2.3.5¹³⁰. The RápidoPGS applies probabilistic finemapping (with Bayes’ factors approximated using Wakefield’s method) to independent LD blocks genome-wide such that variant-wise posterior probabilities of causality can be estimated. This in essence is a ‘shrinkage’ approach to PGS to account for double-counting effects that arise due to correlated effect sizes induced by LD; however, does not inherently require an independent genotyped sample from the GWAS for tuning. Approximate Bayes’ factors in each LD block were derived assuming a prior variance of 0.15, conventionally used for quantitative traits. This parameter choice was compared to a data driven approach to estimating the prior variance based on SNP heritability, as outlined elsewhere¹²⁹. Variant-wise effect sizes are multiplied by their posterior probabilities before PGS calculation, as above.

As we are using a twin cohort, there were two different approaches we utilised to tune the optimal PGS configuration from the LD C+T approach. Firstly, we randomly split the cohort into a training (70% of participants) and test (30% of participants) partition. A linear mixed model was then fit in the training partition with fixed effects of PGS, age, and metabolomics batch and a random effect of family ID to account for twin relatedness. This was applied for retinol measured at each of the three visits for all the P-value thresholds. The marginal R² from a null model with no PGS was subtracted from the full model to infer the best performing PGS in the training partition (mean marginal R² across three visits). The variance explained of the best performing P-value threshold was then estimated using the same approach in the test partition. By way of comparison, we also split the twins into separate unrelated cohorts and used one set as training and one as testing. Fixed effects instead of mixed effects linear regression was then implemented in a similar fashion to above, additionally covaried for five SNP derived principal components. We then repeated all the above for the two probabilistic finemapping weighted PGS (prior variance = 0.15 and data driven prior variance).

Normative modelling

We built a normative model of retinol as a function of age per study visit in TwinsUK. This was achieved using a generalised additive model for location (μ), scale (σ), and shape (GAMLSS)^{131, 132}, implemented in R v4.4.1. The GAMLSS approach is useful in this application as it is semi-parametric and able to account for factors such as heteroskedasticity and non-Gaussian distributions. Our modelling approach can be summarised as follows: we firstly fit a GAMLSS model in the full sample for a variety of GAMLSS distribution families implemented by the gamlss R package v5.4.12. The model on retinol at visit i, i ∈ {1,2,3} for each of the GAMLSS families set the μ term as the first order fractional polynomial of age, along with metabolomics batch as an additional covariate, with the term σ just the first order fractional polynomial of age to model the scale of the distribution. The model fit of each of the tested families was then assessed using the Bayesian information criterion (BIC) and the Akaike information criterion (AIC). We repeated the above also including measured body mass index (BMI) at that timepoint in the μ term by way of comparison. We chose the GAMLSS family (Box-Cox t distribution) through considering the model performance (minimum AIC and BIC) over all three visits (Supplementary Text, Supplementary Figures 11-12). We then split the cohort in half, separating the twins, and fit normative centile curves using the selected GAMLSS family to one half of each batch of retinol measurement, and computed deviations on the other independent sample half on a per batch basis. We repeated this process with the modelling and deviation subsets reversed to compute deviations for all samples. Individuals whose measured retinol was above the model derived 95th percentile were classified as having supra-normal retinol for their age, while those below the 5th percentile were classified as having infra-normal retinol. To guard against overfitting due to the relatedness of twins between the two subsets, we then performed all of the above just using one half of the twins as the full cohort for model fitting, followed by splitting this subset as described above. We tested the relationship between scaled retinol PGS (mean = 0, SD = 1) and infra-normal retinol at each visit was tested in half of the cohort (unrelated) using binomial logistic regression additionally covaried for five SNP derived PCs. The same models were also constructed for supra-normal individuals.

Software and operating systems

The primary analyses in this manuscript were performed either on a MacBook Pro (OS X: Ventura 13.3), an in-house linux cluster (Ubuntu 18.04.5 LTS), or the High-Performance Computing Research Compute Grid of the University of Newcastle [Red Hat Enterprise Linux release 8.1 (Ootpa)]. The primary R version utilised was version 4.1.1 (2021-08-10), with some additional analyses using R version 4.0.3 (2020-10-10) (linux cluster). The Python version utilised was either Python 2.7.17 or Python 3.6.9, depending on the requirements of the analyses.