Genome-wide analysis of health-related biomarkers in the UK Household Longitudinal Study reveals novel associations

Serum biomarker levels are associated with the risk of complex diseases. Here, we aimed to gain insights into the genetic architecture of biomarker traits which can reflect health status. We performed genome-wide association analyses for twenty serum biomarkers involved in organ function and reproductive health. 9,961 individuals from the UK Household Longitudinal Study were genotyped using the Illumina HumanCoreExome array and variants imputed to the 1000 Genomes Project and UK10K haplotypes. We establish a polygenic heritability for all biomarkers, confirm associations of fifty-four established loci, and identify five novel, replicating associations at genome-wide significance. A low-frequency variant, rs28929474, (beta = 0.04, P = 2 × 10−10) was associated with levels of alanine transaminase, an indicator of liver damage. The variant is located in the gene encoding serine protease inhibitor, low levels of which are associated with alpha-1 antitrypsin deficiency which leads to liver disease. We identified novel associations (rs78900934, beta = 0.05, P = 6 × 10−12; rs2911280, beta = 0.09, P = 6 × 10−10) for dihydroepiandrosterone sulphate, a precursor to major sex-hormones, and for glycated haemoglobin (rs12819124, beta = −0.03, P = 4 × 10−9; rs761772, beta = 0.05, P = 5 × 10−9). rs12819124 is nominally associated with risk of type 2 diabetes. Our study offers insights into the genetic architecture of well-known and less well-studied biomarkers.

Scientific RepoRts | 7: 11008 | DOI: 10.1038/s41598-017-10812-1 behavioural, attitudinal, physiological and biomedical variables, including a large panel of the most commonly used clinical biomarkers, have been measured for a representative selection of the sample. This study represents a large sample with very homogenous biomarker measurements, in which recruitment and processing have been carried out consistently and following strict protocols.
Genome-wide association analyses. The genome-wide significance threshold of P < 3.56 × 10 −9 for this study was derived by taking the conventional genome-wide significance threshold (P < 5 × 10 −8 ) divided by the effective number of independent traits analysed (N = 14.05, details in Methods). Across fifteen biomarkers, we observed associations of 54 previously reported loci at this threshold (Fig. 2). This includes a low frequency variant, rs148685782 at 4q31 in the fibrinogen gamma chain gene ( blood clot formation. This variant is a missense mutation and has been previously reported to be associated with fibrinogen levels 7 as well as with hypofibrinogenemia and haemorrhage [8][9][10] . From the discovery phase we carried forward 573 independent (pairwise r 2 < 0.01) variants that were associated with biomarker levels at P < 1 × 10 −5 and were located more than 500 kb away from any known index variant for the respective biomarker. Using data from up to 25,897 samples from 4 independent studies (Supplementary  Table S1), five loci provided evidence of replication and reached P < 3.6 × 10 −9 for the combined analysis of discovery and replication data (   SERPINA1 encodes alpha-1-antitrypsin (AAT), which is a serine protease inhibitor produced in the liver 11 . Low levels of this protein are the hallmark of a genetic disorder called alpha-1 antitrypsin deficiency (A1AD), which leads to liver disease 12 . We identified two novel replicating associations for DHEAS ( Fig. 3B,C). DHEAS is the sulphated form of DHEA, a precursor to major sex-hormones such as testosterone and oestrogen, and is synthesized in the adrenal glands. It is an important marker of adrenal gland function. rs78900934 at chromosome 1p21 (WEAF = 30.9%, beta[SE] = 0.05[0.01], P = 5.9 × 10 −12 ) is located 1 kb upstream of a pseudogene, peptidylprolyl isomerase A pseudogene 7 (PPIAP7). This gene shows a high degree of similarity to cyclophilin A (PPIA), the product of which is involved in a number of biological processes including signal transduction 13 , inflammation 14 and apoptosis 15 . At the second novel locus associated with DHEAS the index variant, rs2911280 at 16q13 (WEAF = 7.5%, beta[SE] = 0.09[0.01], P = 6.0 × 10 −10 ), is located in an intron of the gene encoding c-Maf inducing protein (CMIP), thought to play a role in the T-cell signalling pathway. 16 Two novel replicating associations with HbA1c levels were identified (Fig. 3D,E). HbA1c represents the three-month average plasma glucose concentration and is used to diagnose as well as manage type 2 diabetes. The index variant at 12q13, rs12819124 (WEAF = 46.7%, beta[SE] = −0.03[0.01], P = 4.2 × 10 −9 ) lies in an intron of RP1-228P16.4, a long non-coding RNA. The index variant of the second novel locus, rs761772 at 17q25 (WEAF = 12.4%, beta[SE] = 0.05[0.01], P = 4.9 × 10 −9 ), lies within a non-coding exon in the transmembrane channel-like 6 (TMC6) gene and has been shown to affect the expression of TMC6, as well as TNRC6C antisense RNA 1 (TNRC6C-AS1) and transmembrane channel like 8 (TMC8), in cardiac, thyroid, and vascular tissue, as well as whole blood in the GTEx database 17 .

Discussion
We identify five new biomarker loci, across common and low frequency variants, associated with DHEAS, HbA1c and ALT. We demonstrate polygenic heritability of the majority of biomarkers included in this study and observe large differences in their polygenic component. To our knowledge this is the first report of SNP-based heritability estimates for DHEAS, insulin-like growth factor 1, testosterone and urea. The large sample set with homogeneous biomarker measurements afforded by UKHLS enables reliable estimation for this population. We also identify genetic correlations between several of the biomarkers. Genetic correlation between two traits is an indicator of shared genetic factors and consequently genome-wide pleiotropy. The patterns of heritability and genetic correlations we observe for lipid biomarkers are consistent with previous reports in independent samples 18 . For total and LDL cholesterol, the SNP-based heritability is less than 10% whilst for HDL it is higher at 23.2%. All these estimates represent a lower bound for the narrow sense heritability. Our estimate of the negative genetic correlation between levels of HDL-cholesterol with triglycerides of rg = −0.67 is similar to the estimate derived from an independent study (rg = −0.61) 19 . High levels of triglycerides are mechanistically related to low levels of HDL 20,21 , which could explain the reverse influence of the shared genetic factors on the biomarkers. We show for the first time that polygenic factors for triglyceride are also negatively correlated with DHEAS. There is a statistically significant genetic correlation between CRP and fibrinogen levels, which could be due to shared inflammation pathways. Finally, the genetic correlation we observe between creatinine and urea is a previously unreported  and highly biologically plausible finding as both markers are increased in blood when glomerular filtration rate declines, reflecting impaired kidney function. Characterising the genetic architecture of health-related biomarkers in this way is informative with respect to their biology as well as the design of future association studies. While each known locus individually explains only a small proportion of the variance in biomarker levels, these analyses demonstrate that the joint effect of many variants can be much larger.
We examined less-well studied health-related biomarkers in addition to routine blood measures used in clinical practice. This made it possible to identify novel associations of common and low frequency variants with DHEAS, HbA1c and ALT. These associations could provide novel biological insights. rs2911280, which we found to be associated with DHEAS, is located in an intron of the gene encoding c-Maf inducing protein (CMIP). CMIP is a highly pleiotropic gene, and is associated with several metabolic traits such as adiponectin and HDL cholesterol levels. Cholesterol is a precursor of DHEA in its synthesis process 22 .
rs28929474 at 14q32 is associated with levels of ALT, which is used in clinical practice to assess liver damage. This variant is located in SERPINA1, encoding the serine protease inhibitor alpha-1-antitrypsin (AAT), which is largely produced in the liver. Associations of variants in this gene were previously found for cortisol 23 and height 24 . Mutations of this gene can cause alpha-1 antitrypsin deficiency (A1AD) which can lead to an accumulation of aberrant proteins in hepatocytes causing liver damage 25 . This in turn may elevate levels of ALT, warranting future assessment of the association between this signal and liver-related clinical endpoints.
We identify two novel associations with HbA1c levels. In a lookup using published data from an independent large-scale meta-analysis by the MAGIC consortium 26 , rs12819124 was associated with HbA1c levels with P = 1.8 × 10 −6 . The direction of effect was consistent with our findings. rs12819124 was also nominally associated with risk of type 2 diabetes at P = 0.025 using data from the DIAGRAM study 27 . Moreover, association results from published cohorts suggest a possible pleiotropic association with mental disorders and wellbeing (P = 9.0 × 10 −6 for bipolar disorder and schizophrenia 28 , P = 6.4 × 10 −5 for subjective wellbeing 29 ). No HbA1c association results were available for rs761772 in MAGIC. For a proxy SNP, rs429216 (r 2 = 0.75), the p-value for the association with HbA1c was in the same direction and reached P = 2.7 × 10 −3 .
The UKHLS sample size is modest compared to some of the previous large-scale GWAS meta-analysis efforts (e.g., >45,000 individuals for HbA1c levels 26 ). The relative gain in power leading to novel locus identification in this study can be attributed to several factors. Two of the newly reported signals have relatively low allele frequency (2% and 7.5%, respectively). These were captured here through use of the Illumina HumanCoreExome array and imputation to a comprehensive reference panel consisting of 1000 Genomes combined with the UK10K haplotypes 30 . A further power advantage was afforded by the homogeneous measurement of biomarkers in UKHLS and in two of the replication studies. Each biomarker was measured using the same assay for each sample, and processed on the same machine, avoiding loss of information due to diverse biomarker assays with different sensitivity, dynamic range and detection limit, potentially leading to power reductions 31 .
Larger-scale homogeneous studies and synthesis in massive-scale meta-analyses will help further elucidate the genetic architecture of medically-relevant biomarker traits. Insights into the genetic determinants of population variation in biomarker levels can help us to understand basic processes involved in maintaining health.

Methods
Ethics. Participants gave informed written consent for their blood to be taken and stored for future scientific analysis. The UKHLS has been approved by the University of Essex Ethics Committee and the nurse data collection by the National Research Ethics Service (10/H0604/2). Study population. The Table 3. Association results of replicating novel signals. function: variant functional consequence; nearest gene: gene nearest to lead variant with 500Kb from either side; chr: chromosome; EA/NEA: effect allele/non-effect allele; EAF; effect allele frequency; beta(SE), p-value: effect size (standard error) and p-value; N: total number of individuals analysed for this variant; r 2 : imputation accuracy.
interview. As recruitment was household based, the study contains related individuals. The study includes phenotypical data for a representative sample of participants for a wide range of social and economic indicators as well as a biological sample collection encompassing biometric, physiological, biochemical, and haematological measurements and self-reported medical history and medication use (https://www.understandingsociety. ac.uk/d/100/7251_User_Guide_Health_Assmt_w2_w3.pdf?1392855567). For each participant non-fasting blood samples were collected through venepuncture, were centrifuged to separate plasma and serum, aliquoted and frozen at −80 °C. DNA has been extracted and stored for genetic analyses. For replication, data were available for 5533 individuals from ELSA 32 , 9888 from Fenland 33 (Supplemental Table 1), 7621 from HRS (http://hrsonline.isr.umich.edu) 34 , 2859 from NCDS 35 . These studies have been described in detail elsewhere. Sample collection were carried out consistently and analysed by the same laboratories for UKHLS, ELSA and NCDS.
Biomarker measurements. In total, biomarker data were successfully obtained from 13,107 eligible individuals who gave consent to give blood samples to be stored for future analysis (https://www.understandingsociety.ac.uk/d/154/7251-UnderstandingSociety-Biomarker-UserGuide-2014.pdf?1418057881). All biomarkers were measured from serum (non-fasting), using a variety of suitable assays, and the majority analysed on a single Roche P module analyser 36 . On each machine Internal Quality Controls (IQC) were at regular intervals per day. External Quality Assurance (EQA) systems were in place to monitor all tests. Phenotype transformations and exclusions. The measurements for biomarkers used in the association analyses were prepared according to protocols from the largest genetic association study published for each specific trait at the time when analyses commenced, details for which are available in Supplementary Table S2. Genotyping. In total, 10,484 UKHLS samples have been typed using the Illumina Infinium HumanCoreExome BeadChip Kit ® (12v1-0). This array contains a set of >250,000 highly informative genome-wide tagging single nucleotide polymorphisms as well as a panel of functional (protein-altering) exonic markers, including a large proportion of low-frequency (MAF 1-5%) and rare (MAF <1%) variants. Genotype calling was performed with the gencall algorithm using GenomeStudio (Illumina Inc.). For quality control (QC) we excluded individuals based on the following criteria: sample call rate <98%, autosomal heterozygosity outliers (>3SD), gender mismatches, duplicates as established by identity by descent (IBD) analysis (PI_HAT > 0.9). Individuals with non-European ancestry were also excluded. For this we estimated the genomic kinship between all pairs of individuals along with 1000 Genomes Project data. These were converted to distances and subjected to multidimensional scaling. Prior to variant QC, we first mapped all 538,448 variants to the human reference genome build 37. Variants with Hardy-Weinberg equilibrium p-value < 1 × 10 −4 , call rate below 98% or poor genotype clustering values (<0.4) were excluded, leaving 525,314 variants passing QC. For typed variants in our GWAS analyses that were brought forward for replication we inspected cluster plots manually using Scattershot 0.75 beta (Supplementary Fig. S1). All QC procedures were carried out using PLINK (v1.07) and R. Imputation. We imputed our genotype data using a combined reference panel consisting of 7,562 haplotypes from the UK10K project and 2,184 haplotypes from 1000 Genomes phase 3. Details regarding the creation of this combined imputation panel are described elsewhere 37 . Prior to imputation, we first pre-phased using SHAPEIT (v2.r). Data were then imputed using IMPUTE2 (v2.3.0), resulting in an initial set of 38,310,212 variants. Variants with an IMPUTE info score <0.4, and variants with a Hardy-Weinberg p-value < 1 × 10 −4 were excluded, leaving 26,851,013 variants for analysis. Data availability. The  Statistical analysesh. Heritability analyses and genetic correlations. The proportion of trait variance explained by the genotyped and imputed variants was estimated using the GREML method as implemented in the GCTA software 38, 39 (v1.26). We included all variants with minor allele frequency (MAF) > 0.01. We excluded variants with imputation accuracy less than 0.4. We computed the genetic relationship matrix (GRM) for each autosome and then used GCTA to combine them into one matrix. We excluded relatives from the estimation by filtering based on the GRM using a threshold of 0.1 after inspecting the distribution. This led to the exclusion of 672 individuals for this analysis. We also performed bivariate REML analysis in order to estimate genetic correlations between different biomarkers 40 . We applied a Bonferroni adjusted significance threshold using the effective number of traits for the heritability analyses and using the number of pairs based on the effective number of traits for the genetic correlation analyses.
Association analyses. The association analyses were carried out using a multivariate linear mixed model to account for relatedness as implemented in GEMMA (v0.95). QQ plots and genomic inflation factors, as well as Manhattan plots for traits where we identified novel associations are displayed in Supplementary Fig. S2. Replication analyses were carried out in R and PLINK, following the same trait preparation protocols as used in the discovery stage. The association summary statistics from the replication analyses, as well as the combined discovery and replication stage were meta-analysed using a fixed-effects inverse variance weighted approach implemented in METAL (v2011-03-25). We calculated an adjusted genome-wide significance threshold, for the effective number of traits, as several of our biomarkers have correlated levels. The effective number of traits was derived by computing the eigenvalues for the correlation matrix of the 20 biomarkers (effective number: 14.05). The routinely used GWAS threshold of p < 5 × 10 −8 was then adjusted for this using the Bonferroni approach: 5 × 10 −8 /14.05 = 3.56 × 10 −9 .
Power calculations. We carried out power calculations using Quanto (v1.2.4), for discrete per-variant frequency and (standardised) effect sizes combinations, representative of variants identified in the discovery. Per-trait and per selected variant power analyses showed that we would minimally need 5,000 to 15,000 samples to replicate our variants with P < 3.56 × 10 −9 , two-sided, for testosterone levels, whereas the largest replication sample of 25,000 to 60,000 would be needed for eGFR (Fig. 4).
Selection of replication SNPs, and criteria for novel loci. For replication we selected independent SNPs (LD r2 < 0.1), with MAF > 0.01 and a discovery p-value of P < 1 × 10 −5 and at least > 500 Kb away from the nearest known reported index SNP for a given trait. We also took forward independent rare variants with a MAF ≤ 0.01 that were typed and reached P < 1 × 10 −5 , regardless whether they represented known associations for a given trait. Known index SNPs for all biomarkers analysed in this study were obtained through the GWAS catalog 41 (accessed August 4, 2016) > , supplemented by manual searches in PubMed.