1. Genetics and Genomics
Download icon

GWAS of three molecular traits highlights core genes and pathways alongside a highly polygenic background

  1. Nasa Sinnott-Armstrong  Is a corresponding author
  2. Sahin Naqvi
  3. Manuel Rivas
  4. Jonathan K Pritchard  Is a corresponding author
  1. Department of Genetics, Stanford University, United States
  2. Department of Chemical and Systems Biology, Stanford University, United States
  3. Department of Biomedical Data Sciences, Stanford University, United States
  4. Department of Biology, Stanford University, United States
Research Article
  • Cited 9
  • Views 5,227
  • Annotations
Cite this article as: eLife 2021;10:e58615 doi: 10.7554/eLife.58615

Abstract

Genome-wide association studies (GWAS) have been used to study the genetic basis of a wide variety of complex diseases and other traits. We describe UK Biobank GWAS results for three molecular traits—urate, IGF-1, and testosterone—with better-understood biology than most other complex traits. We find that many of the most significant hits are readily interpretable. We observe huge enrichment of associations near genes involved in the relevant biosynthesis, transport, or signaling pathways. We show how GWAS data illuminate the biology of each trait, including differences in testosterone regulation between females and males. At the same time, even these molecular traits are highly polygenic, with many thousands of variants spread across the genome contributing to trait variance. In summary, for these three molecular traits we identify strong enrichment of signal in putative core gene sets, even while most of the SNP-based heritability is driven by a massively polygenic background.

Introduction

One of the central goals of genetics is to understand how genetic variation, environmental factors, and other sources of variation, map into phenotypes. Understanding the mapping from genotype to phenotype is at the heart of fields as diverse as medical genetics, evolutionary biology, behavioral genetics, and plant and animal breeding.

The nature of the genotype-to-phenotype mapping has been a key motivating question ever since the start of modern genetics in the early 1900s. In those early days, the genetic basis of phenotypic variation was debated between the Mendelians, who were interested in discrete monogenic phenotypes, and the biometricians, who believed that Mendelian genetics were incompatible with the continuous distributions observed for height and many other traits. Those battles were largely resolved by Fisher’s 1918 paper showing that a large number of Mendelian loci, each with proportionally weak effects, can approximate a continuous trait (Fisher, 1918; Provine, 2001; Barton et al., 2017; Visscher and Goddard, 2019). Taken to its extreme, this type of model is referred to as the ‘infinitesimal model’, and it laid the foundations for the growth of quantitative genetics in the 20th century (Lynch and Walsh, 1998).

Despite the importance of the infinitesimal model in the development of the field, for a long time this was mainly a theoretical abstraction. Even though some authors predicted early on that certain human diseases might be polygenic (Penrose, 1953; Gottesman and Shields, 1967), it was recognized that even a few loci (<10) can approximate infinitesimal predictions (Thoday and Thompson, 1976; McGuffin and Huckle, 1990). Thus, prior to the GWAS era it was entirely unclear how many loci would actually affect complex traits in practice (Risch et al., 1999; Visscher and Goddard, 2019). For example, in a 1989 review of quantitative genetics, Barton and Turelli wrote that ‘we still do not know whether the number of loci responsible for most genetic variation is small (5–20) or large (100 or more)’ (Barton and Turelli, 1989). Consistent with this, practioners of human genetics in the pre-GWAS era expected that we might be looking for a small handful of genes per trait; in the 1990s, this motivated hundreds of small studies of complex traits that were only powered to detect large-effect loci. In one typical example, Risch and Merikangas’ foundational 1996 paper on association mapping computed the power for common variants with relative risks between the alternate homozygotes ranging from 2.25 to 16 (Risch and Merikangas, 1996): effect sizes that we now know were unrealistically high.

The advent of GWAS, starting around fifteen years ago, completely transformed our understanding of the genetic basis of a wide variety of human complex traits and diseases (Claussnitzer et al., 2020). While early GWAS studies showed the power of this approach to identify significant and replicable signals, it quickly became clear that the lead variants generally explain only small fractions of the heritability of the corresponding traits (Weedon et al., 2008; Goldstein, 2009). The limited explanatory power of the detected loci became known as the ‘‘mystery of missing heritability’’ (Manolio et al., 2009): a mystery that was largely resolved by work showing that most of the heritability is due to the presence of many sub-significant causal variants (Purcell et al., 2009; Yang et al., 2010). Subsequent work has shown that for most traits the bulk of the SNP-based heritability is spread widely and surprisingly uniformly across the genome (Loh et al., 2015; Shi et al., 2016; O’Connor, 2020), and that most complex traits are in fact enormously polygenic, with various studies estimating >10,000 or even >100,000 causal variants per trait (Zhang et al., 2018; Frei et al., 2019; O'Connor et al., 2019).

Why are complex traits so polygenic?

These findings raise a mechanistic question of how to understand the biological processes that link genotype to phenotype. How should we understand the observations that the lead variants for a typical trait contribute only a small fraction of the heritability, while most of the heritability is driven by tens of thousands of variants, presumably mediated through thousands of genes?

One hypothesis for polygenicity is the observation that many disease or behavioral endpoints are impacted by multiple distinct processes, or endophenotypes (Turkheimer, 2000; Gottesman and Gould, 2003): for example, diabetes is affected by lipids, adiposity and β-cell function (Udler, 2019). For such traits, the genetic basis of the endpoint phenotype is expected to reflect the genetics of all the intermediate phenotypes. While this is surely true, it seems unlikely to fully resolve the question of why specific phenotypes can be affected by tens of thousands of variants, unless the endophenotypes themselves are highly polygenic. Indeed, as we will show in this paper, even relatively ‘simple’ molecular traits such as urate levels can be hugely polygenic, implying that we need additional explanations for high polygenicity.

In recent work, we proposed an alternative model for understanding extreme polygenicity, namely that it may be a consequence of the architecture of gene regulatory networks (Boyle et al., 2017; Liu et al., 2019). Work from several groups has shown that, for an average gene, most of the heritability in gene expression results from large numbers of small trans effects (Price et al., 2008; Liu et al., 2019). Building on this observation, we proposed a conceptual model in which there is a set of ‘core’ genes, defined as genes with a direct effect on the trait that is not mediated through regulation of other genes. Meanwhile, other genes that are expressed in trait-relevant cell types are referred to as ‘peripheral’ genes, and can matter if they affect the expression of core genes. By this definition, transcription factors (TFs) are considered peripheral, but we refer to TFs with coordinated effects on multiple core genes as ‘master regulators’ to acknowledge their special roles. This model primarily denotes different categories of genes (rather than variants) with respect to their roles in trait variation; we discuss how these distinctions may apply to variant effect sizes later in the text.

We proposed that variants near core genes contribute only a small fraction of the heritability, and that instead most trait variance is due to huge numbers of trans-regulatory effects from SNPs with cis-effects on peripheral genes. In what we referred to as the ‘omnigenic’ extreme, potentially any gene expressed in trait-relevant cell types could affect the trait through effects on core gene expression. (Note that this does not mean that every gene would in fact have associated variants, as presumably the distribution of peripheral gene effect sizes would be centered on zero, and in practice not all genes have regulatory variants).

While the omnigenic model is broadly consistent with observations on cis and trans heritability of expression (Liu et al., 2019), it has been difficult to evaluate the model in detail because for most diseases and other traits we know little in advance about which genes are likely to be directly involved in disease biology. Recent efforts to systematically nominate core genes have primarily relied upon associations identified in rare, monogenic disorders (Vuckovic et al., 2020); while promising, such approaches are inherently limited by the ability to discover rare gene-disease associations, which can depend upon a number of factors. Furthermore, we still have highly incomplete information about cellular regulatory networks and trans-eQTLs.

Here, we focus on three molecular traits that are unusually tractable in order to gain insights into the roles of core genes. This work illustrates two key parts of the model: (1) the existence and identity of sets of core genes for each trait and (2) that the core genes contribute only a small fraction of the heritability. We do not directly assess the role of trans-regulatory networks for these traits as well-powered trans-eQTL data do not exist for the relevant cell types.

GWAS of model traits: three vignettes

We investigate the genetic architecture underlying variation in three molecular traits: serum urate, IGF-1, and testosterone levels. For each of these traits, we know a great deal in advance about the key organs, biological processes and genes that might control these traits.

This stands in contrast to many of the traits that have been studied extensively with GWAS, such as schizophrenia (Ripke et al., 2014; Ripke et al., 2020; which is poorly understood at the molecular level) or height (Wood et al., 2014; where we understand more of the underlying biology, but for which a large number of different biological processes contribute variance). We do now know various examples of core genes or master regulators for specific traits (e.g. Sekar et al., 2016; Small et al., 2011; Small et al., 2018), but there are few traits where we understand the roles of more than a few of the lead genes. Among the clearest examples in which a whole suite of core genes have been identified are for plasma lipid levels (e.g. Liu et al., 2017; Lu et al., 2017; Hoffmann et al., 2018, reviewed by Dron and Hegele, 2016; Liu et al., 2019); and for inflammatory bowel disease (de Lange et al., 2017).

As described in more detail below, we performed GWAS for each of these traits in around 300,000 white British individuals from the UK Biobank (Bycroft et al., 2018). For all three traits many of the most significant hits are highly interpretable–a marked difference from GWAS of typical disease traits. While these three molecular traits highlight different types of lead genes and molecular processes, they also have strikingly similar overall architectures: the top hits are generally close to genes with known biological relevance to the trait in question, and all three traits show strong enrichment in relevant gene sets. Most of these genes would be considered core genes (or occasionally master regulators) in the sense of Liu et al., 2019.

At the same time, however, variants near the lead genes and pathways explain only a modest fraction of the heritability. Aside from one major-effect variant for urate, the lead pathways explain ∼10% of the SNP-based heritability. Instead, most of the SNP-based heritability is due to a highly polygenic background, which we conservatively estimate as being due to around 10,000 causal variants per trait.

In summary, these three molecular traits provide points of both contrast and similarity to the architectures of disease phenotypes. From one point of view they are clearly simpler, successfully identifying known biological processes to an extent that is highly unusual for disease GWAS. At the same time, the most significant hits sit on a hugely polygenic background that is reminiscent of GWAS for more-complex traits.

Results

Our analyses make use of GWAS results that we reported previously on blood and urine biomarkers (Sinnott-Armstrong et al., 2021), with minor modifications. In the present paper, we report four primary GWAS analyses: urate, IGF-1, and testosterone in females and males separately. Prior to each GWAS, we adjusted the phenotypes by regressing the measured phenotypes against age, sex (urate and IGF-1 only), self-reported ethnicity, the top 40 principal components of genotype, assessment center and month of assessment, sample dilution and processing batch, as well as relevant pairwise interactions of these variables (Materials and methods).

We then performed GWAS on the phenotype residuals in White British participants. For the GWAS we used variants imputed using the Haplotype Reference Consortium with MAF >0.1% and INFO >0.3 (Materials and methods), yielding a total of 16M variants. The final sample sizes were 318,526 for urate, 317,114 for IGF-1, 142,778 for female testosterone, and 146,339 for male testosterone. One important goal of our paper is to identify the genes and pathways that contribute most to variation in each trait. For gene set-enrichment analyses, we annotated gene sets using a combination of KEGG (Kanehisa and Goto, 2000) and previous trait-specific reviews, as noted in the text. We considered a gene to be ‘close’ to a genome-wide significant signal if it was within 100 kb of at least one lead SNP with p<5e-8. The annotations of lead signals on the Manhattan plots were generally guided by identifying nearby genes within the above-described enriched gene sets, or occasionally other strong nearby candidates.

Genetics of serum urate levels

Urate is a small molecule (C5H4N4O3) that arises as a metabolic by-product of purine metabolism and is released into the blood serum. Serum urate levels are regulated by the kidneys, where a set of transporters shuttle urate between the blood and urine; excess urate is excreted via urine. Urate is used as a clinical biomarker due to its associations with several diseases. Excessively high levels of urate can result in the formation of needle-like crystals of urate in the joints, a condition known as gout. High urate levels are also linked to diabetes, cardiovascular disease, and kidney stones.

The genetics of urate have been examined previously by several groups (Woodward et al., 2009; Köttgen et al., 2013; Nakayama et al., 2017; Nakatochi et al., 2019; Boocock et al., 2019; Tin et al., 2019 and recently reviewed by Major et al., 2018). The three strongest signals for urate lie in solute carrier genes: SLC2A9, ABCG2, and SLC22A11/SLC22A12. A recent trans-ancestry analysis of 457 k individuals identified 183 genome-wide significant loci (Tin et al., 2019); their primary analysis did not include UK Biobank. Among other results, this study highlighted genetic correlations of urate with gout and various metabolic traits; tissue enrichment signals in kidney and liver; and genetic signals at the master regulators for kidney and liver development HNF1A and HNF4A.

Performing GWAS of urate in the UK Biobank data set, we identified 222 independent genome-wide significant signals, summarized in Figure 1A (further details in Supplementary file 1). Remarkably, six of the 10 most significant signals are located within 100 kb of a urate solute transport gene. A recent review identified 10 genes that are involved in urate solute transport in the kidneys (Wright et al., 2010; Anzai et al., 2007); in addition to the six transporters with extremely strong signals, two additional transporters have weaker, yet still genome-wide significant signals (Figure 1B). Hence, GWAS highlights eight out of 10 annotated urate transporters, although some transporters were originally identified using early GWAS for urate levels. The two genes in the pathway that do not have hits (SMCT1 and SMCT2; also known as SLC5A8 and SLC5A12) do not directly transport urate, but instead transport monocarboxylate substrates for URAT1 to increase reabsorption rate (Bobulescu and Moe, 2012) and thus may be less direct regulators of urate levels.

Figure 1 with 1 supplement see all
Genetic basis of serum urate variation.

(A) Genome-wide associations with serum urate levels in the UK Biobank. Candidate genes that may drive the most significant signals are indicated; in most cases in the paper, the indicated genes are within 100 kb of the corresponding lead SNPs. (B) Eight out of 10 genes that were previously annotated as being involved in urate transport (Wright et al., 2010; Anzai et al., 2007) are within 100 kb of a genome-wide significant signal. The signal at MCT9 is excluded from figure and enrichment due to its uncertain position in the pathway (Fisel et al., 2018). (C) Urate SNP-based heritability is highly enriched in kidney regulatory regions compared to the genome-wide background (analysis using stratified LD Score regression). Other tissues show little or no enrichment after removing regions that are active in kidney. See Figure 1—figure supplement 1 for the uncorrected analysis.

Among the other top hits, five are close to transcription factors involved in kidney and liver development (HNF4G, HNF1A, HNF4A, HLF, and MAF). These are not part of a globally enriched gene set, but recent functional work has shown that the associated missense variant in HNF4A results in differential regulation of the urate solute carrier ABCG2 (Tin et al., 2019), while the MAF association has been shown to regulate SLC5A8 (Leask et al., 2018). Finally, two other loci show large signals: a missense variant in INHBC, a TGF-family hormone, and a variant in/near GCKR, a glucose-enzyme regulator. Both variants have highly pleiotropic effects on many biomarkers, although the mechanisms pertaining to urate levels are unclear.

While most of the top hits are likely associated with kidney function, we wanted to test whether other tissues contribute to the overall SNP-based heritability (Figure 1C). To this end, we used stratified LD Score regression to estimate the polygenic contribution of regulatory regions in 10 previously defined tissue groupings (Finucane et al., 2015). Serum urate SNP-based heritability was most-highly enriched in kidney regulatory regions (29-fold compared to the genome-wide average SNP, p=1.9e-13), while other cell types were enriched around 8-fold (Figure 1—figure supplement 1; see also Tin et al., 2019). We hypothesized that the enrichment for other tissues might be driven by elements shared between kidney and other cell types. Indeed, when we removed active kidney regions from the regulatory annotations for other tissues, this eliminated most of the signal found in other cell types (Figure 1C). Thus, our analysis supports the inference that most serum urate heritability is driven by kidney regulatory variation.

Finally, while these signals emphasize the role of the kidneys in setting urate levels, we wanted to test specifically for a role of urate synthesis (similar to recent work on glycine [Wittemans et al., 2019]). The urate molecule is the final step of purine breakdown; most purines are present in tri- and monophosphates of adenosine and guanosine, where they act as signaling molecules, energy sources for cells, and nucleic acid precursors. The breakdown pathways are well known, including the genes that catalyze these steps (Figure 2A).

Modest enrichment of signals among genes involved in urate biosynthesis.

(A) Urate is a byproduct of the purine biosynthesis pathway. The urate component of each molecule is highlighted. (B) The same pathway indicating genes that catalyze each step. Genes with a genome-wide significant signal within 100 kb are indicated in red; numbers in gray indicate the presence of additional genes without signals. Pathway adapted from KEGG.

Overall, we found that genes in the urate metabolic pathway show a modest enrichment for GWAS hits relative to all annotated, protein coding genes as a background (2.1-fold, p=0.017; Figure 2B). XDH, which catalyzes the last step of urate synthesis, has an adjacent GWAS hit, as do a number of upstream regulators of urate synthesis. Nonetheless, the overall level of signal in the synthesis pathway is modest compared to that seen for kidney urate transporters, suggesting that synthesis, while it plays a role in common variation in urate levels, is secondary to the secretion pathway. In contrast, remarkably, nearly all of the kidney urate transporter genes are close to genomewide significant signals; there are additional strong signals in kidney transcription factors, as well as a strong polygenic background in kidney regulatory regions.

Genetics of IGF-1 levels

Our second vignette considers the genetic basis of IGF-1 (insulin-like growth factor 1) levels. The IGF-1 protein is a key component of a signaling cascade that connects the release of growth hormone to anabolic effects on cell growth in peripheral tissues (Laron, 2001). Growth hormone is produced in the pituitary gland and circulated around the body; in the liver, growth hormone triggers the JAK-STAT pathway leading, among other things, to IGF-1 secretion. IGF-1 binding to IGF-1 receptor, in turn, activates the RAS and AKT signaling cascades in peripheral tissues. IGF-1 is used as a clinical biomarker of growth hormone levels and pituitary function, as it has substantially more stable levels and a longer half-life than growth hormone itself. The growth hormone–IGF axis is a conserved regulator of longevity in diverse invertebrates and possibly mammals (van Heemst, 2010). In humans, both low and high levels of IGF-1 have been associated with increased mortality from cancer and cardiovascular disease (Burgers et al., 2011). We note that while IGF-1 rises sharply in puberty, our analyses are focused on middle-aged individuals. IGF-1 is a major effect locus for body size in dogs (Sutter et al., 2007), and IGF-1 levels are positively associated with height in UK Biobank (Figure 3—figure supplement 1).

Previous GWAS for IGF-1, using up to 31,000 individuals, identified around half a dozen genome-wide significant loci (Kaplan et al., 2011; Teumer et al., 2016). The significant loci included IGF-1 itself and a signal close to its binding partner IGFBP3.

In our GWAS of serum IGF-1 levels in 317,000 unrelated White British individuals, we found a total of 354 distinct association signals at genome-wide significance (Figure 3, further details in Supplementary file 2). Eight of the most significant signals are key parts of the IGF-1 pathway (Figure 4). The top hit is an intergenic SNP between IGFBP3 and another gene, TNS3 (Supplementary file 2; p=1e-837). IGFBP3 encodes the main transport protein for IGF-1 and IGF-2 in the bloodstream (Firth and Baxter, 2002). The next most significant hits are at the IGF-1 locus itself and at its paralog IGF-2. Two other lead hits are associated with the IGF transport complex IGFBP: IGFALS, which is an IGFBP cofactor that also binds IGF-1 in serum (Baxter et al., 1989), and PAPPA2, a protease which cleaves and negatively regulates IGFBPs (Overgaard et al., 2001). Three other lead hits lie elsewhere in the growth hormone–IGF axis: GHSR is a pituitary-expressed receptor for the signaling protein ghrelin which negatively regulates the growth hormone (GH) signaling pathway upstream of IGF-1 (Laron, 2001); and FOXO3 and RIN2 lie in downstream signaling pathways (Stitt et al., 2004).

Figure 3 with 5 supplements see all
Genetic basis of IGF-1 variation.

Manhattan plot showing the locations of major genes associated with IGF-1 levels in the IGF-1 pathway (yellow), transcription factor (blue), pleiotropic gene (red), or unknown function (black) genes sets.

GWAS hits in the IGF-1 pathway.

Bolded and colored gene names indicate that the gene is within 100 kb of a genome-wide signficant hit. Gray names indicate absence of a genome-wide signficant hit; gray numbers indicate that multiple genes in the same part of the pathway with no hit. Superscript numbers indicate that multiple genes are located within the same locus and hence may not have independent hits. (A) Upstream pathway that controls regulation of IGF-1 secretion into the bloodstream. (B) Downstream pathway that controls regulation of IGF-1 response.

Additional highly significant hits that are not directly involved in the growth hormone–IGF pathway include the liver transcription factor HNF1A (also associated with urate [Tin et al., 2019]); variants near two genes–GCKR and KLF14–that are involved in many biomarkers, although to our knowledge the mechanism is unclear; and variants at two additional genes CENPW and ZNF644.

Given the numerous lead signals in the IGF-1 signaling cascade, we sought to comprehensively annotate all GWAS hits within the cascade and its sub-pathways. We compiled lists of the genes from KEGG and relevant reviews from five major pathways in the growth hormone–IGF axis (Figure 4, Materials and methods). Four of the five pathways show extremely strong enrichment of GWAS signals. The first pathway regulates growth hormone secretion, acting in the pituitary to integrate ghrelin and growth hormone releasing hormone signals and produce growth hormone. This pathway shows strong enrichment, with 14 out of 32 genes within 100 kb of a genome-wide significant signal (7.3-fold enrichment, Fisher’s exact p=5.4e-7). The second pathway, IGF-1 secretion, acts in the liver, where growth hormone triggers JAK-STAT signaling, leading to IGF-1 production and secretion (Dehkhoda et al., 2018). This pathway again shows very strong enrichment of GWAS signals (10/14 genes, 23-fold enrichment, p=4.9e-8). The third pathway, serum balance of IGF, relates to IGF-1 itself, and its paralogs, as well as other binding partners and their regulators in the serum. Here 10/18 genes have GWAS hits (11.7-fold enrichment, p=1.5e-6).

We also considered two downstream signaling pathways that transmit the IGF signal into peripheral tissues. Most notably, many of the genes in the AKT branch of the IGF-1 signaling cascade were close to a genome-wide significant association including FOXO3 (9/31 genes; 3.8-fold enrichment, p=0.002). In contrast, the RAB/MAPK/RAS pathway was not enriched overall (p=0.59), although one key signaling molecule (RIN2) in this pathway was located at one of the strongest hits genome-wide. The observation of strong signals downstream of IGF-1 suggests the presence of feedback loops contributing to IGF-1 regulation. This is consistent with work proposing negative feedback from downstream pathways including AKT and MAPK to growth hormone activity (Li et al., 2009).

Lastly, given that most of the strongest hits lie in the same pathway, we were curious whether there might be evidence for epistatic or non-additive interactions. Experiments in molecular and model organism biology regularly find interaction effects between genes that are close together in pathways (Tong et al., 2004; Scanga et al., 2000; Bassik et al., 2013; Fischer et al., 2015; Wang et al., 2014), but evidence for epistatic interactions between GWAS variants is extraordinarily rare (Ritchie and Van Steen, 2018), potentially due to GWAS hits lying in unrelated pathways, having modest effect sizes, or most often not representing the causal variant. We found no signal of epistasis among the 77 most significant (p<1e-20) lead SNPs for IGF-1 (Figure 3—figure supplement 2), and weak enrichment of signal among the top 38 urate lead SNPs (Figure 3—figure supplement 2 inset; see Materials and methods). Similarly, IGF-1 lead SNPs (p<5e-8) showed weak, global inflation of test statistics for departures from additivity (e.g. domininance or recessivity) (Figure 3—figure supplement 3), while the two most significant urate hits showed significant minor dominant (SLC2A9) and minor recessive (ABCG2) effects that were nevertheless substantially smaller than the additive effects (Figure 3—figure supplement 3 inset). Genome-wide paired difference tests for the SLC2A9 variant showed no signal (Figure 3—figure supplement 4). Building upon previous studies (Wei et al., 2014; Zaitlen et al., 2013) that have found little evidence of epistasis or dominance in human populations, these results indicate that non-linear genotype effects, while likely present to some degree, are substantially weaker than additive components, even when the strongest effects are concentrated in the same biological pathways and would thus be more likely to show epistasis.

In summary for IGF-1, we found 354 distinct associations that surpass genome-wide significance. The lead variants show strong enrichment across most components of the growth hormone-IGF axis, including the downstream AKT signaling arm, suggesting regulatory feedback. Among the strongest hits we also find involvement of one transcription factor (HNF1A) and two other genes of unclear functions (GCKR and KLF14) that have pleiotropic effects on multiple biomarkers, perhaps due to overall effects on liver and kidney development.

Testosterone

Our third vignette describes the genetic basis of testosterone levels. Testosterone is a four carbon-ring molecule (C19H28O2) that functions as an anabolic steroid and is the primary male sex hormone. Testosterone is crucial for the development of male reproductive organs and secondary sex characteristics, while also having important functions in muscle mass and bone growth and density in both females and males (Tracz et al., 2006; Herbst and Bhasin, 2004). Circulating testosterone levels range from about 0.3 to 2 nmol/L in females and 8 to 33 nmol/L in males (Figure 5—figure supplement 1).

Testosterone is synthesized from cholesterol as one possible product of the steroid biosynthesis pathway. Synthesis occurs primarily in the testis in males, and in the ovary and adrenal glands in females. Testosterone production is stimulated by the hypothalmic-pituitary-gonadal (HPG) axis: gonadotropin-releasing hormone (GnRH) signals from the hypothalamus to the pituitary to cause production and secretion of luteinizing hormone (LH); LH in turn signals to the gonads to produce testosterone. The HPG axis is subject to a negative feedback loop as testosterone inhibits production of GnRH and LH by the hypothalamus and pituitary to ensure tight control of testosterone levels (Javorsky et al., 2017). Testosterone acts on target tissues via binding to the androgen receptor (AR) which in turn regulates downstream genes. Approximately half of the circulating testosterone (∼40% in males, ∼60% in females [Dunn et al., 1981]) is bound to sex hormone binding globulin (SHBG) and is generally considered non-bioavailable. Testosterone breakdown occurs primarily in the liver in both females and males.

Previous GWAS for serum testosterone levels studied up to 9000 males, together finding three genome-wide significant loci, the most significant of which was at the SHBG gene (Ohlsson et al., 2011; Jin et al., 2012). While this paper was in preparation, two studies reported large-scale GWAS of testosterone levels in UKBB individuals, finding significant sex-specific genetic effects (Flynn et al., 2021; Ruth et al., 2020). Previous studies of young adults found minimal correlation of salivary testosterone levels between opposite-sex dizygotic twins (Grotzinger et al., 2018). In our preliminary analysis, we found that testosterone shows minimal genetic correlation between the sexes, in contrast to other biomarkers including urate and IGF-1 (Figure 7—figure supplement 1). We therefore performed sex-stratified GWAS of testosterone, in contrast to the combined analysis used for urate and IGF-1.

Here, we performed testosterone GWAS in UKBB females (N = 142,778) and males (N = 146,339) separately. We discovered 79 and 127 independent genome-wide significant signals in females and males, respectively (Figure 5, further details in Supplementary file 34). We note that a recent paper reported larger numbers of independent genome-wide significant signals (245 and 231 in females and males, respectively); this was likely due to the inclusion of individuals with broader European ancestry, as well as a less stringent definition of independence used by Ruth et al (Ruth et al., 2020).

Figure 5 with 1 supplement see all
Manhattan plots for testosterone.

(A) Females. (B) Males. Notice the low overlap of lead signals between females and males. FAM9A and FAM9B have been previously proposed as the genes underlying the KAL1 locus (Ohlsson et al., 2011).

In females, six of the most significant signals are close to genes involved in testosterone biosythesis (Figure 5A); together these results suggest that the steroid biosynthesis pathway is the primary controller of female testosterone levels. Among these, the top hit is at a locus containing three genes involved in hydroxylation of testosterone and estrone, CYP3A4, CYP3A5, and CYP3A7 (Kandel et al., 2017; Lee et al., 2003; Kuehl et al., 2001). Two other lead hits (MCM9 and FGF9) are involved in gonad development (Lutzmann et al., 2012; Wood-Trageser et al., 2014; Colvin et al., 2001).

Strikingly, and in agreement with recent studies and in agreement with recent studies (Flynn et al., 2021; Ruth et al., 2020), the lead hits in males are largely non-overlapping with those from females. Overall, the male hits affect a larger number of distinct processes. Three of the most significant signals affect the steroid biosynthesis pathway (SRD5A2, UGT2B15, and AKR1C); three are involved in either upstream activation (NR0B2) (Vega et al., 2015) or downstream signaling (the androgen receptor, AR, and its co-chaperone FKBP4), respectively; and two have been implicated in the development of the GnRH-releasing function of the hypothalamus (KAL1) (Franco et al., 1991) or the gonads (NR2F2) (Qin et al., 2008). However, the largest category, including the most significant hit overall, is for a group of eight distinct variants previously shown to affect sex hormone binding globulin (SHBG) levels (Coviello et al., 2012). SHBG is one of the main binding partners for testosterone–we will discuss the significance of SHBG below.

Steroid biosynthesis

Given our observation of numerous lead hits near steroid hormone biosynthesis genes, we curated the female and male hits in the KEGG pathway (Figure 6). We observed that nearly all major steps of the pathway contained a gene near a genome-wide significant SNP in either females or males: 31 out of 61 genes are within 100 kb of a genome-wide significant signal in males, females or both. Indeed, the KEGG steroid hormone pathway shows strong enrichment for signals in both females and males (26-fold enrichment, p=2.5e-8 in females; 11-fold enrichment, p=1.2e-4 in males; Figure 6—figure supplement 1). While this pathway shows clear enrichment in both females and males, the major hits do not overlap. At two loci, AKR1C and PDE2A, female and male hits co-occur at the same locus, but are localized to different SNPs (Figure 7—figure supplement 1). More broadly, male hits and female hits tend to occur in different parts of the steroid hormone biosynthesis pathway: catalytic steps involved in progestagen and corticosteroid synthesis and metabolism only showed hits in females, while most male hits were concentrated within androgen synthesis, either upstream or downstream of testosterone itself (Figure 6).

Figure 6 with 2 supplements see all
Pathway diagram for steroid hormone biosynthesis showing GWAS hits for females and males.

The text color indicates genes within 100 kb of a genome-wide significant hit for females (orange), males (blue), or both females and males (black). Gray gene names or numbers indicates genes with no hits. Colored superscripts indicate multiple genes from the same locus (and hence may reflect a single signal). ‘S*’ indicates that an additional, sulfonated metabolite, along with the catalytic step and enzymes leading to it, is not shown. Pathway from KEGG; simplified based on a similar diagram in Wikipedia, 2012.

Genetics of testosterone regulation in males versus females

One remarkable feature of the testosterone data is the lack of sharing of signals between females and males. This is true for genome-wide significant hits, for which there is no correlation in the effect sizes among lead SNPs (Figure 7A), as well as genome-wide, as the global genetic correlation between females and males is approximately zero (Figure 7—figure supplement 1).

Figure 7 with 5 supplements see all
Sex differences in genetic variation in testosterone.

(A) When comparing lead SNPs (p<5e-8 ascertained in either females or males), the effects are nearly non-overlapping between females and males. Other traits show high correlations for the same analysis (see urate and SHBG in inset). (B) Schematic of HPG axis signaling within the hypothalamus and pituitary, with male GWAS hits highlighted. These variants are not significant in females. (C) Global genetic correlations, between indicated traits (estimated by LD Score regression). Thickness of line indicates strength of correlation, and significant (p<0.05) correlations are in bold. Note that LH genetic correlations are not sex-stratified due to small sample size in the UKBB primary care data (N = 10,255 individuals). (D) Proposed model in which the HPG axis and SHBG-mediated regulation of testosterone feedback loop is primarily active in males. Abbreviations for all panels: SHBG, sex hormone-binding globulin; CBAT, calculated bioavailable testosterone; LH, luteinizing hormone.

As we show below, two aspects of testosterone biology can explain these extreme sex differences in genetic architecture. First, the hypothalmic-pituitary-gonadal (HPG) axis plays a more significant role in regulating testosterone production in males than in females. This is due to sex differences in both endocrine signaling within the HPG axis and the tissue sources of testosterone production. Second, SHBG plays an important role in mediating the negative feedback portion of the HPG axis in males but not in females.

To assess the role of HPG signaling, we searched for testosterone GWAS hits involved in the transmission of feedback signals through the hypothalamus and pituitary (Figure 7B, genes reviewed in Skorupskaite et al., 2014). We also considered hits from GWAS of calculated bioavailable testosterone (CBAT), which refers to the non-SHBG-bound fraction of total teststerone that is free or albumin-bound, and can be inferred given levels of SHBG, testosterone, and albumin and assuming experimentally determined rate constants for binding (Vermeulen et al., 1999). CBAT GWAS thus controls for genetic effects on total testosterone that are mediated by SHBG production.

We found hits for both male testosterone and male CBAT throughout the HPG signaling cascade (Figure 7B). These include genes involved in the direct response of the hypothalamus to testosterone (AR, FKBP4) (Smith et al., 2005); modulation of the signal by either autoregulation (TAC3, TACR3) (Skorupskaite et al., 2014) or additional extrinsic endocrine signals (LEPR) (Ahima et al., 1996; Barash et al., 1996); downstream propagation (KISS1) (Messager et al., 2005) and the development of GnRH-releasing neurons in the hypothalamus (KAL1, CHD7) (Cariboni et al., 2004; Layman et al., 2011); and LH-releasing gonadotropes in the pituitary (GREB1) (Li et al., 2017). All these hits showed more significant effects on CBAT as compared to total testosterone (Figure 7—figure supplement 3), suggesting that their primary role is in regulating bioavailable testosterone.

Importantly, these HPG signaling hits do not show signals in females. To further investigate the different roles of the HPG axis in males versus females, we performed GWAS of LH levels using UKBB primary care data (N = 10,255 individuals). (Recall that LH produced by the pituitary signals to the gonads to promote sex hormone production.) If HPG signaling is important for testosterone production in males but not females, variants affecting LH levels should affect testosterone levels in males but not females. Consistent with this, we found significant positive genetic correlation between LH and male but not female testosterone (male rg=0.27, p=0.026; female rg=0.084, p=0.49; Figure 7C). These results were similar when considering measured testosterone and LH levels rather than genetic components thereof (Supplementary file 5).

Two known features of the HPG axis can explain the lack of association in females. First, the adrenal gland, which is not subject to control by HPG signaling, produces ∼50% of serum testosterone in females. Consistent with this idea, GWAS hits for female testosterone cluster in steroid hormone pathways involving progestagen and corticosteroid synthesis (Figure 6), processes known to occur largely in the adrenal. Female testosterone hits are also specifically enriched for high expression in the adrenal gland relative to male testosterone hits (Figure 7—figure supplement 4).

Second, for the ovaries, which produce the remaining ∼50% of serum testosterone in females, the net effect of increased LH secretion on testosterone production is expected to be diminished. This is because the pituitary also secretes follicle-stimulating hormone (FSH), which in females stimulates aromatization of androgens (including testosterone) into estrogens (Ulloa-Aguirre and Michael Conn, 2014). In males, FSH does not stimulate androgen aromatization but is instead required for sperm production. Consistent with differential roles of FSH, a previously described GWAS hit for menstrual cycle length at FSHB (Laisk et al., 2018) shows suggestive association with testosterone in females but not males (Supplementary file 6).

In addition to the role of HPG signaling, the presence of many SHBG-associated variants among the top hits in male testosterone suggests that SHBG also underlies many of the sex-specific genetic effects (Figure 5B). We found high positive genetic correlation between female and male SHBG, as well as between SHBG and total testosterone in males but not females (Figure 7C). Additionally, we found a significant negative genetic correlation between SHBG and CBAT in both females and males, but of a far larger magnitude in females than males (Figure 7C). Together, these observations suggest that while SHBG regulates the bioavailable fraction of testosterone in the expected manner in both females and males, there is subsequent feedback in males only, where decreased CBAT leads to increased total testosterone.

Combining these observations, we propose that increased SHBG leads to decreased bioavailable testosterone in both females and males, and in males this relieves the negative feedback from testosterone on the hypothalamus and pituitary gland, ultimately allowing LH production and increased testosterone production (Figure 7D). The lack of SHBG-mediated negative feedback in females is likely due in part to the overall weaker action of the HPG axis, as well as the fact that female testosterone levels are too low to effectively inhibit the HPG axis. This idea is supported experimental manipulations of female testosterone, which result in significant reductions of LH only when increasing testosterone levels to within the range typically found in males (Serafini et al., 1986).

In summary, we find that many of the top signals for female testosterone are in the steroid biosynthesis pathway, and a smaller number relate to gonadal development. In contrast, the lead hits for male testosterone reflect a larger number of processes, including especially SHBG levels and signaling components of the HPG axis, in addition to biosynthesis and gonadal development. These differences in the genetic architecture of female and male testosterone are so extreme that these can be considered unrelated traits from a genetic point of view.

Polygenic architecture of the three traits

We have shown that the lead signals for all three traits are highly concentrated near core genes and core pathways. As an additional confirmation of these enrichments, we found that (i) random sets of SNPs matched to urate, IGF-1, or testosterone GWAS hits showed far lower overlap with the corresponding core pathways (Figure 7—figure supplement 5), and (ii) an alternative approach (de Leeuw et al., 2015) showed highly similar gene-set enrichments to those we observed (Supplementary file 7). Given this observation, we wondered whether these traits might be genetically simpler than typical complex diseases–most of which are highly polygenic, and for which the lead pathways contribute relatively little heritability (Boyle et al., 2017; Shi et al., 2016).

To address this, we first estimated how much of the SNP-based heritability is explained by variation at genes in enriched pathways (see Supplementary files 810 for pathways and genes used). We used HESS to estimate the SNP-based heritability in each of 1701 approximately-independent LD blocks spanning the genome (Shi et al., 2016; Berisa and Pickrell, 2016). Plotting the cumulative distribution of SNP-based heritability across the genome revealed that, across all four traits, most of the genetic variance is distributed nearly uniformly across the genome (Figure 8A).

Figure 8 with 14 supplements see all
Despite clear enrichment of core genes and pathways, most SNP-based heritability for these traits is due to the polygenic background.

(A) Cumulative distribution of SNP-based heritability for each trait across the genome (estimated by HESS). The locations of the most significant genes are indicated. Insets show the fractions of SNP-based heritability explained by the most important genes or pathways for each trait. (B) Estimated fractions of SNPs with non-null associations, in bins of LD Score (estimated by ashR). Each point shows the ashR estimate in a bin representing 0.1% of all SNPs. The inset text indicates the estimated fraction of variants with a non-null marginal effect, that is, the fraction of variants that are in LD with a causal variant. (C) Simulated fits to the data from (B). X-axis truncated for visualization as higher LD Score bins are noisier. Simulations assume that π1 of SNPs have causal effects drawn from a normal distribution centered at zero (see Materials and methods). The simulations include a degree of spurious inflation of the test statistic based on the LD Score intercept. Other plausible assumptions, including clumpiness of causal variants, or a fatter-tailed effect distribution would increase the estimated fractions of causal sites above the numbers shown here.

In aggregate, core genes contribute modest fractions of SNP-based heritability, with the exception of the SLC2A9 locus, which HESS estimates is responsible for 20% of the SNP-based heritability for urate. Aside from this outlier gene, the core pathways contribute between approximately 1–11 percent of the SNP-based heritability.

Numbers of causal variants

We next sought to estimate how many causal variants are likely to contribute to each trait (Zhang et al., 2018; Frei et al., 2019; O'Connor et al., 2019). This is fundamentally a challenging problem, as most causal loci have effect sizes too small to be confidently detected. As a starting point we used ashR, which is an empirical Bayes method that estimates the fraction of non-null test statistics in large-scale experiments (Stephens, 2017). As described previously, we stratified SNPs from across the genome into bins of similar LD Score; we then used ashR to estimate the fraction of non-null associations within each bin (Boyle et al., 2017). (For this analysis, we used the 2.8M SNPs with MAF >5%.) We interpret this procedure as estimating the fraction of all SNPs in a bin that are in LD with a causal variant.

For each trait, the fraction of non-null tests increases from low levels in the lowest LD Score bins to above 50% in the highest LD Score bins. Overall we estimate that around 45–50% of SNPs are linked to a non-zero effect variant for urate, IGF-1 and male testosterone, and 30% for female testosterone (Figure 8B). These estimates were robust to halving the sample size of the input GWAS, and were substantially higher than for randomized traits (simulated by permuting the IGF-1 and urate phenotypes) (Figure 8—figure supplement 1).

We next conducted simulations to understand how these observations relate to the numbers of causal variants (Figure 8C). To make this identifiable, we assume that a fraction 1-π1 of all SNPs have an effect size that is exactly zero, while the remainer (π1) draw their effect size from a single normal distribution with mean zero. Our goal is to estimate π1. We simulated phenotypes for the UK Biobank individuals assuming a range of values of π1 (Materials and methods). Causal variants were chosen uniformly at random from among the 4.4M SNPs with MAF >1%; effect sizes were simulated from a normal distribution with mean zero, and variances set to produce the observed SNP heritabilities (0.3 for urate, IGF-1, and male testosterone, and 0.2 for female testosterone). We also allowed for a degree of over-inflation of the test statistics (i.e. allowing for an inflation factor as in Genomic Control [Devlin and Roeder, 1999])–this was important for fitting the positive ashR estimates at low LD Scores. We then matched the simulations to the observed ashR results to approximate the numbers of causal variants.

Overall, our estimates range from 0.1% of all 4.4M variants with MAF >1% in female and male testosterone (∼4000 causal sites) to 0.3% of variants for urate (∼12,000 causal sites). These results imply that all four traits are highly polygenic, though considerably less so than height (for which we estimate 2%, or 80,000 causal sites in UK Biobank; Figure 8—figure supplements 2 and 4).

Furthermore, there are three reasons to suspect that these numbers may be underestimates. First, causal variants are likely to be clumped in the genome instead of being uniformly distributed; simulations with clumping require a larger number of causal variants to match the data (Figure 8—figure supplement 5). Second, if the distribution of effect sizes has more weight near zero and fatter tails than a normal distribution, this would imply a larger number of causal variants (see analysis assuming a T-distribution, Figure 8—figure supplement 6). Third, stratified LD Score analysis of the data suggests that some of the apparent evidence for overinflation of the test statistics (Supplementary file 11) may in fact be due to a higher proportion of causal variants occurring in lower LD Score bins (Gazal et al., 2017) rather than population stratification, as the annotation-adjusted intercepts for all traits but height are consistent with 1 (no population stratification).

We note that the proportion of causal variants estimated by ashR is substantially lower in low-MAF bins, even in infinitesimal models, presumably due to lower power (Figure 8—figure supplements 7 and 8). We overcame this by using a parametric fit, which is robust to inflation of test statistics (Figure 8—figure supplements 9 and 10); the resulting estimates were relatively similar, albeit slightly higher, than when using the simulation-matching method (Figure 8—figure supplement 4). We note that it is still critical to match samples by heritability and sample size, as in the simulation method (Figure 8—figure supplement 11), and to use correct covariates in the GWAS (Figure 8—figure supplement 12).

As an alternative approach, we used the program GENESIS, which uses a likelihood model to fit a mixture of effect sizes using 1–2 normal components, and a null component (Zhang et al., 2018; Supplementary file 12). Assuming a single normal distribution, the results for the molecular traits were very similar to our results: male testosterone 0.1%; female testosterone 0.2%; urate 0.3%; IGF-1 0.4%. The GENESIS results for a mixture of two normal distributions resulted in a significantly higher overall likelihood, and estimates roughly threefold higher than our estimates: male testosterone 0.6%; female testosterone 0.7%; urate 1.1%; IGF-1 1.1%. GENESIS estimates for height were lower than ours (0.6% and 1.2%, respectively); it is possible that there is a downward bias at high polygenicity as GENESIS estimates for a simulated fully infinitesimal model were 2.7%.

In summary this analysis indicates that for these molecular traits, around 10–15% of the SNP-based heritability is due to variants in core pathways (and in the case of urate, SLC2A9 is a major outlier, contributing 20% on its own). However, most of the SNP-based heritability is due to a much larger number of variants spread widely across the genome, conservatively estimated at 4000–12,000 common variants for the biomarkers and 80,000 for height.

Discussion

In this study, we examined the genetic basis of three molecular traits measured in blood serum: a metabolic byproduct (urate), a signaling protein (IGF-1), and a steroid hormone (testosterone). We showed that unlike most disease traits, these three biomolecules have strong enrichments of genome-wide significant signals in core genes and related pathways. At the same time, other aspects of the data are reminiscent of patterns for complex common diseases, including high polygenicity, little indication of allelic dominance or epistasis, and enrichment of signals in tissue-specific regulatory elements spread across the genome.

Our main results are as follows.

  • Urate: The largest hits for urate are in solute carrier genes in the kidneys that shuttle urate in and out of the blood and urine. Remarkably, eight out of ten annotated urate transporters have genomewide significant signals. A single locus, containing SLC2A9, is responsible for 20% of the SNP-based heritability. While the urate transport pathway was previously known to be enriched in GWAS hits (Tin et al., 2019), we further demonstrate that the purine biosynthetic pathway, from which urate is produced as a byproduct, is modestly enriched for signals (2.1-fold). Several master regulators for kidney and liver development are among the most significant hits. Aside from SLC2A9, the overall SNP-based heritability is primarily driven by variants in kidney regulatory regions, both shared across cell types and not.

  • IGF-1: IGF-1 is a key component of a signaling cascade that links growth hormone released from the pituitary to stimulation of cell growth in peripheral tissues. We identified 354 independent genome-wide significant signals. The strongest signals lie in genes that interact directly with IGF-1, including IGFBP3, as well as in the IGF1 gene itself. More generally, we see striking enrichment of hits throughout the growth hormone-IGF cascade–this includes especially the upper parts of the cascade, which regulate IGF-1 release, but also in downstream components of the cascade as well, suggesting a feedback mechanism on IGF-1 levels. These pathway-level enrichments were not identified in previous, less well-powered GWAS of IGF-1 levels (Teumer et al., 2016).

  • Testosterone: In contrast to urate, testosterone shows clear enrichment of signals within the steroid biosynthesis pathway (26-fold in females, 11-fold in males). Remarkably, the genetic basis of testosterone is almost completely independent between females and males, as reported recently (Flynn et al., 2021; Ruth et al., 2020). In females, the lead hits are mostly involved in synthesis. In males, in addition to hits in the synthesis pathway, we see signals throughout the hypothalamic-pituitary-gonadal (HPG) axis which regulates testosterone production in the gonads, as well as in variants that regulate SHBG. Furthermore, in males, increased SHBG reduces negative feedback between testosterone levels and the HPG axis, thereby increasing total serum testosterone. These results provide a mechanistic explanation of the sex differences in testosterone genetics, in addition to showing that GWAS hits can reveal the core biology of a trait even in the context of vastly differing genetic architecture between the sexes.

  • Polygenic background. For each of these traits, the core genes and pathways contribute only a modest fraction of the total SNP-based heritability. Aside from SLC2A9 for urate, the most important core pathways contribute up to about 10% of the total SNP-based heritability. We estimated the numbers of causal variants under a model where causal variants have a normal effect-size distribution. We estimate that there are around 4000–12,000 common variants with non-zero effects on these traits. Using the same method, we estimated about 80,000 causal sites for height. These estimates are likely conservative as several of our assumptions may lead us to underestimate the true values.

Understanding the architecture of complex traits

One of our motivations in this study was to use these three traits as models to extend our understanding of the genetic architecture and types of genes underlying complex traits.

Many of the advances of 20th century genetics resulted from focused study of the functions of major-effect mutations; this principle has been extended in the GWAS era into interpreting the impact of lead signals. And yet, at the same time, most heritability is driven by the polygenic background of small effects at genes that are not directly involved in the trait. The overwhelming importance of the polygenic background of thousands of small effects is a striking discovery of modern GWAS, and demands explanation as it does not fit neatly into standard conceptual models of the relationship between genotype and phenotype.

As discussed in the Introduction, our group recently proposed a simplified conceptual model to understand this phenomenon (Boyle et al., 2017; Liu et al., 2019). We proposed that for any given trait there is a limited set of core genes that are directly involved in the biology of the phenotype, but that most of the heritability is due to SNPs with cis-effects on other (peripheral) genes that are expressed in the same tissues; these in turn affect core genes via trans-regulatory networks. Thus far, it has been difficult to fully test this model because, in general, we do not know many of the genes that may have direct effects in any given trait. We also generally have very limited knowledge of trans-regulatory networks.

The present paper helps to fill part of this gap by studying the genetic basis of three molecular traits where we can a priori identify a large number of core genes and their associated pathways. Thus, our work provides concrete examples of how genetics can affect core biology to an extent that is usually difficult to achieve for disease traits. Furthermore, consistent with the model and our previous analyses of gene expression (Liu et al., 2019), we find that the known core pathways contribute only a modest fraction of the SNP-based heritability, and that the bulk of the heritability is driven by thousands of variants spread across much of the genome.

That said, it remains difficult to test the second part of the model, that is, that most of the heritability passes through trans-regulatory networks. This problem is challenging because we currently lack the sample sizes necessary for inferring trans-regulatory networks in any tissue or cell type, with the possible exception of whole blood. Secondly, its likely that the relevant networks may be active only in particular cell types or at particular times in development, which makes the estimation of trans-regulatory networks even more difficult. However, one promising approach has recently yielded results consistent with the trans-regulatory part of the omnigenic model. Võsa et al have shown that genome-wide polygenic scores for several traits correlate with the expression levels of core genes for those traits, as would be predicted by the model (Võsa et al., 2018). Nonetheless, the trans-regulatory component of the model requires further work.

Another type of explanation for high polygenicity comes from the observation that many traits and diseases are affected by multiple biological processes. Thus, any variants that affect those intermediate processes can potentially be detected in GWAS of the endpoint trait (Turkheimer, 2000; Gottesman and Gould, 2003; Bittante et al., 2012; Pickrell et al., 2016; Udler, 2019). While this process undoubtedly contributes to the polygenicity of many endpoint traits, our data suggest it is unlikely that this type of process drives high polygenicity for these molecular traits. Notably, for urate, we estimated ∼12,000 causal variants, and showed that the vast majority of the SNP-based heritability likely acts through the kidneys. Thus, any explanation for the high polygenicity of urate must presumably depend on the role of genetic variation on kidney function in general, and urate transport in particular.

The huge polygenicity of complex traits also raises questions about how to extract biological insight from GWAS. If there are tens of thousands of associated variants, acting through thousands of genes, then presumably most of these will not be especially helpful for understanding mechanisms of disease (Goldstein, 2009). (In contrast, for constructing polygenic scores, we do in fact care about all variants, as small effects drive most of the phenotypic variance.) This raises the question of how to use GWAS to identify the genes that are actually most proximal to function. This is of course a question that many in the field have wrestled with, for a wide variety of traits (de Leeuw et al., 2015; Pers et al., 2015). Overall, we can expect that the most significant variants will usually point to biologically important genes for the corresponding trait. That said, there are many reasons why significance is not a fully reliable indicator of gene importance: significance depends on both the variant effect size and its allele frequency; the allele frequency is a random outcome of genetic drift and, moreover, selection tends to lower frequencies of the most important variants (Simons et al., 2018; O'Connor et al., 2019); lastly the effect size of the variant depends not only on the importance of the gene for the trait, but also on the magnitude of that variant’s effect on the gene (e.g. as a cis-eQTL). Furthermore, some genes that are biologically important may be entirely missed because they do not happen to have common functional variants. Nonetheless, given all these caveats, we found that for these three molecular traits the lead GWAS hits were indeed highly enriched for core genes, consistent with work for other traits where many of the lead variants are interpretable (Lu et al., 2017; Liu et al., 2017; de Lange et al., 2017).

In summary, we have shown that for three molecular traits, the lead hits illuminate core genes and pathways to a degree that is highly unusual in disease or complex trait GWAS. By doing so they illustrate which processes may be most important for trait variation. For example, for urate, kidney transport is more important than biosynthesis, while for testosterone, biosynthesis is important in both sexes but especially in females. However, in other respects, the GWAS data here are reminiscent of more-complex traits: in particular most trait variance comes from a huge number of small effects at peripheral loci. These vignettes help to illustrate the architecture of complex traits, with lead variants that are directly involved in trait biology alongside a massively polygenic background.

Materials and methods

Population definition

Request a detailed protocol

We defined our GWAS population as a subset of the UK Biobank (Bycroft et al., 2018). We use ∼337,000 unrelated White British individuals as our cohort, filtering based on sample QC characteristics as previously described (Sinnott-Armstrong et al., 2021):

  1. Used to compute principal components (used_in_pca_calculation column).

  2. Not marked as outliers for heterozygosity and missing rates (het_missing_outliers column).

  3. Do not show putative sex chromosome aneuploidy (putative_sex_chromosome_aneuploidy column).

  4. Have at most 10 putative third-degree relatives (excess_relatives column).

  5. Finally, we used the in_white_British_ancestry_subset column in the sample QC file to define the subset of individuals in the White British cohort.

Trait definition

Request a detailed protocol

We perform trait normalization and quality control similarly to previous work (Sinnott-Armstrong et al., 2021). Trait measurements are first log-transformed, then adjusted for genotype principal components, age indicator variables, sex, 5 year age (‘approximate age’) by sex interactions, self-identified ethnicity, self-identified ethnicity by sex interactions, fasting time, estimated sample dilution factor, assessment center, genotyping batch, 20-tile of time of sampling, month of assessment, and day of assay.

Then, individuals were subset to the GWAS population (defined above), separated by sex for testosterone measurements. The final sample sizes were 318,526 for urate, 317,114 for IGF-1, 142,778 for female testosterone, and 146,339 for male testosterone.

GWAS

Request a detailed protocol

We performed GWAS in plink2 alpha using the following command (data loading arguments removed for brevity):

plink2 --glm cols=chrom,pos,ref,alt,alt1,ax,a1count,totallele,a1freq, 
    machr2,firth,test,nobs,beta,se,ci,tz,p omit-ref 
  --covar-variance-standardize 
  --remove [non-White-British, related White British or excluded] 
  --keep [male, female, or all] 
  --geno 0.2 --hwe 1e-50 midp --threads 16

GWAS were then filtered to observed allele frequency greater than 0.001 and INFO score greater than 0.3 for further analyses.

GWAS for paired difference epistasis

Request a detailed protocol

A GWAS was performed in two subsets of individuals – those with two C alleles at rs16890979 (N = 295209) and those with two T alleles at rs16890979 (N = 30184). The following command was used:

plink2 --glm cols=chrom,pos,ref,alt,a1freq,firth,test, 
      nobs,beta,se,ci,tz,p hide-covar 
  --hwe 1e-50 midp --keep [rs16890979 CC or TT individuals] 
  --remove [non-White British] --geno 0.1 --maf 0.001

With covariates including adjusting for age, age squared, genotyping array, and 20 principal components. The residual urate levels, already adjusted for age, sex, global principal components, and technical covariates (Methods) were used as input.

After GWAS completed, SNPs valid in both CC and TT individuals were compared for betas using a paired difference Z test. The test statistic was then converted to a p-value using a standard normal distribution.

LH trait definition

Request a detailed protocol

LH levels were extracted from UK Biobank primary care data using code XM0lv. Separately, LH levels extracted using code XE25I were also included for phenotypic correlation analyses. The median level across observations and log number of observations were recorded for covariate correction below. Individuals with median observations more than 10 times the interquartile range away from the median of medians were discarded. Once these individuals were removed, individuals with observations more than four standard deviations from the resulting mean were also discarded.

For the primary LH code XM0lv, the distribution of raw, cleaned, and covariate-adjusted phenotype values were respectively:

Distribution of raw (left), cleaned (middle), and covariate-adjusted (right) phenotype values for primary luteinizing hormone (LH) code XMOlv.

For the secondary LH code XE25I, the distribution of raw, cleaned, and covariate-adjusted phenotype values were respectively:

Distribution of raw (left), cleaned (middle), and covariate-adjusted (right) phenotype values for secondary LH code XE25I.

For GWAS, the cleaned phenotypes were log-transformed and adjustments were used as covariates.

LH GWAS

Request a detailed protocol

Age, sex, genotyping array, 10 PCs, log number of observations in primary care, and which primary care code produced a given observation were used as covariates.

We performed GWAS in plink2 alpha using the following command (data loading arguments removed for brevity):

plink2 --glm cols=chrom,pos,ref,alt,alt1,ax,a1count,totallele,a1freq, 
    machr2,firth,test,nobs,beta,se,ci,tz,p hide-covar omit-ref 
  --covar-variance-standardize 
  --remove [non-White-British, related White British or excluded] 
  --keep [all White British] 
  --geno 0.2 --hwe 1e-50 midp --maf 0.005 --vif 999

We also performed GWAS of LH code XE25I in a sex stratified fashion using the following command:

plink2 --glm cols=chrom,pos,ref,alt,alt1,ax,a1count,totallele, 
    a1freq,machr2,firth,test,nobs,beta,se,ci,tz,p 
  hide-covar omit-ref --covar-variance-standardize --remove <non-White-British> 
  --geno 0.2 --hwe 1e-50 midp --threads {threads} --maf 0.001 --vif 999;

On genotyped SNPs and imputed variants with a minor allele frequency greater than 1% in the White British as a whole.

GWAS were then filtered to MAF >1% and INFO >0.7. These higher threshold were chosen to reflect the much smaller sample size in the GWAS.

GWAS hit processing

Request a detailed protocol

To evaluate GWAS hits, we took the list of SNPs in the GWAS and ran the following command using plink1.9:

plink --bfile [] --clump [GWAS input file] --clump-p1 1e-4 --clump-p2 1e-4 
  --clump-r2 0.01 --clump-kb 10000 --clump-field P --clump-snp-field ID

We then took the resulting independent GWAS hits and examined them for overlap with genes. In addition, for defining the set of SNPs to use for enrichment analyses, we greedily merged SNPs located within 0.1 cM of each other and took the SNP with the minimum p-value across all merged lead SNPs. In this way, we avoided potential overlapping variants that were driven by the same, extremely large, gene effects.

Gene proximity

Request a detailed protocol

We annotated all genes in any Biocarta, GO, KEGG, or Reactome MSigDB pathway as our full list of putative genes (in order to avoid pseudogenes and genes of unknown function), and included the genes within each corresponding pathway as our target set. This resulted in 17,847 genes. We extended genes by 100 kb (truncating at the chromosome ends) and used the corresponding regions, overlapped with SNP positions, to define SNPs within range of a given gene. Gene positions were defined based on Ensembl 87 gene annotations on the GRCh37 genome build.

Pathway enrichment of GWAS hits

Request a detailed protocol

GWAS hit pathway enrichment was evaluated using Fisher’s exact test. For each pathway for a given trait (Supplementary files 810), genes were divided into those within the pathway and those outside; and separately into genes within 100 kb of a GWAS hit and not. A 2 × 2 Fisher’s exact test was used to estimate the total enrichment for GWAS hits around genes of interest.

For female and male testosterone, we noticed a number of GWAS loci with multiple paralogous enzymes within the synthesis pathway (e.g. AKR1C, UGT2B, CYP3A). To avoid double counting GWAS hits when testing enrichment at such loci, we instead considered the number of GWAS hits (within 100 kb of any pathway gene as above) normalized to the total genomic distance covered by all genes (±100 kb) in the pathway. A Poisson test was used to compare the rate parameter for this GWAS hit/Mb statistic between genes in a given pathway and all genes not in the pathway.

To quantify pathway enrichment expected from random sets of SNPs not associated with a phenotype, we used SNPSnap (Pers et al., 2015) with default settings to obtain 1000 sets of equally-sized random SNPs matched to urate, IGF-1, or testosterone hits in terms of LD, minor allele frequency, and genic distance. For each set of random, matched SNPs, we determined the number of core genes within 100 kb as for the true set of GWAS hits.

To quantify pathway enrichments using an alternative approach, we used MAGMA (de Leeuw et al., 2015) with a 10 kb gene window and with the default competitive mode. We tested enrichment for all gene sets in Biocarta, GO, KEGG, or Reactome MSigDB, as well as additional curated sets of core genes for the three traits.

Partitioned heritability

Request a detailed protocol

Partitioned SNP-based heritability estimates were generated using LD Score regression (Finucane et al., 2015). The BaselineLD version 2.2 was used as a covariate, and the 10 tissue type LD Score annotations were used as previously described (Finucane et al., 2015) in a multiple regression setup with all cell type annotations and the baseline annotations.

Pathway heritability estimation

Request a detailed protocol

We evaluated SNP-based heritability in pathways using two distinct strategies. Initially, we used partitioned LD Score regression (Finucane et al., 2015) but found that the estimates were somewhat noisy, likely because most pathways contain few genes. As such, we used alternative fixed-effect models for which there is increased power.

Next, we calculated the SNP-based heritability in a set of 1701 approximately independent genomic blocks spanning the genome (Berisa and Pickrell, 2016) using HESS (Shi et al., 2016). Next, we overlapped blocks with genes in each pathway. The SNP-based heritability estimates for all blocks containing at least one SNP within 100 kb of a pathway gene were summed to estimate the SNP-based heritability in a given pathway. Pathway definitions were assembled based on a combination of KEGG pathways, Gene Ontology categories, and manual curation based on relevant reviews.

Causal SNP simulations

Request a detailed protocol

All imputed variants with MAF >1% in the White British (4.1M) were used as a starting set of putative causal SNPs. Individual causal variants were chosen at random, with a fraction P of them marked as causal. Each causal variant was assigned an effect size:

βN(0,1)

For our simulations, we used P{0.0001,0.001,0.003,0.01,0.03}.

Next, GCTA was used to simulate phenotypes based on the marked causal variants, using the following command:

gcta64 --simu-qt --simu-causal-loci CausalVariantEffects 
  --simu-hsq 0.3 --bfile UKBBGenotypes"

Producing predicted phenotypes with SNP-based heritability h2=0.3. GWAS were run within both the full set of 337,000 unrelated White British individuals and a randomly downsampled 50%, to approximate the sex-specific GWAS used for Testosterone, across the set of putative causal SNPs. GWAS for the traits, as well as a random permuting across individuals of urate and IGF-1 to act as negative controls, were repeated on this subset of variants as well. In this way, we have a directly comparable set of simulated traits to use, along with the corresponding true traits and negative controls, to ascertain causal sites in the genome.

For the infinitesimal simulations, instead plink was used to generate polygenic scores on the basis of the random assignment of effect sizes to SNPs, and these were then normalized with N(0,σ2) environmental noise such that h2 was the given target SNP-based heritability.

Causal SNP count fitting procedure using ashr

Request a detailed protocol

LD Scores for the 489 unrelated European-ancestry individuals in 1000 Genomes Phase III (Bulik-Sullivan et al., 2015) were merged with the GWAS results along with LD Scores derived from unrelated European ancestry participants with whole genome sequencing in TwinsUK. TwinsUK LD Scores are used for all analyses. Then variants were filtered by minor allele frequency to either greater than 1%, greater than 5%, or between 1% and 5%. Remaining variants were divided into 1000 equal sized bins, along with 5000 and 200 bin sensitivity tests. Within each bin, the ashR estimates of causal variants, as well as the mean χ2 statistics, were calculated using the following line of R:

data %>% filter(pmin(MAF, 1-MAF) > min.af, pmin(MAF, 1-MAF) < max.af) %>% 
    mutate(ldBin = ntile(ldscore, bins)) %>% group_by(ldBin) %>% 
    summarize(mean.ld = mean(ldscore), se.ld=sd(ldscore)/sqrt(n()), 
        mean.chisq = mean(T_STAT**2, na.rm=T), 
        se.chisq=sd(T_STAT**2, na.rm=T)/sqrt(sum(!is.na(T_STAT))), 
        mean.maf=mean(MAF), 
        prop.null = ash(BETA, SE)$fitted_g$pi[1], n=n())

Thus, the within-bin χ2 and proportion of null associations π0 were each ascertained. Next, these fits were plotted as a function of mean.ld to estimate the slope with respect to LD Score, and true traits were compared to simulated traits, described below.

We use two fixed simulated heritabilities, h2=0.3 and h2=0.2, to approximately capture the set of heritabilites observed among our biomarker traits. Traits with true SNP-based heritability among variants with MAF >1% different than their closest simulation might have causal site count over-estimated (for htrue2>hsim2) or under-estimated (for htrue2<hsim2). In addition, most traits in reality have more than zero SNPs with MAF <1% contributing to the SNP-based heritability. Thus, we take these estimates as approximate and conservative.

Effect of population structure on causal SNP estimation

Request a detailed protocol

We expect that population structure might lead to test statistic inflation for causal variant and genetic correlation estimates (Berg et al., 2019). To evaluate this, we performed GWAS for height using no principal components, and evaluated the causal variant count (Figure 8—figure supplement 12).

This suggests that the test statistic inflation is an important parameter in the estimation of causal variants, as is intuitive. As such, we generated estimated SNP counts for five different inflation values (0.9, 1, 1.05, 1.1, and 1.2) and plotted all of them, under the assumption that the best fitting intercept would have the most calibrated estimates. Plots are replicated across these intercepts in the sensitivity analyses shown, as in Figure 8—figure supplement 9.

Evaluating the calibration of causal SNP proportion estimation

Request a detailed protocol

To evaluate calibration of causal SNP estimates, in addition to using simulated traits as the controls, we also generated a randomized control by shuffling the SHBG phenotype values across individuals (Figure 8—figure supplement 3). We performed this analysis using urate and IGF-1 to similar effect (data not shown).

This suggests that the causal variant counts are well calibrated for the randomized traits, even though they lack structure with respect to covariates.

Effect of sample size on causal SNP estimation

Request a detailed protocol

It is important to note that these estimates are still likely power limited even in a study as large as UK Biobank. We make this note on the basis of observed π0 for MAF>5% variants being uniformly higher than 1%<MAF<5% variants in both simulations and observed data for high causal variant counts (Figure 8—figure supplement 8).

As such, we anticipate that future studies with larger samples will yield increased, but asymptotic, estimates of causal SNP percentages among common variants, and treat our estimates as conservative bounds.

Particularly for height (Figure 8—figure supplement 2), while the uncalibrated estimates with the full sample are substantially higher than the half sample, the calibrated estimates are nearly identical. This suggests that trait polygenicity might be an important factor in determining the power of this method at different sample sizes, as height is known to be highly polygenic (Shi et al., 2016).

Effect of binned variant count on causal SNP estimation

Request a detailed protocol

It is possible that the ashR algorithm itself, and not the GWAS, are the power limited step of the analysis. To evaluate this, we ran ashR on 200, 1000, and 5000 equally sized bins along the LD Score axis. We found that increasing bin counts both decrease the standard errors and the intercepts (Figure 8—figure supplement 13) and recommend as many bins as is practical.

Effect of minor allele frequency on causal SNP estimation

Request a detailed protocol

Because we only simulated causal effects among SNPs with MAF >1%, we were concerned that variant effect bins might be biased by the minor allele frequency cutoff. We previously ran with higher MAF cutoffs (25% and 40%) as calibrations on an earlier version of the model, and observed uniformly larger causal SNP percentages. We saw relative robustness to lower thresholds, but overall the fraction of causal variants was lower in the lower MAF bins (Figure 8—figure supplement 7).

Effect of concentrated SNPs on causal SNP estimation

Request a detailed protocol

For each variant, the megabase bin it is contained within was used as a proxy for SNPs in local LD. A within-megabase causal SNP percentage parameter:

PBeta(α,α/ρ)

was chosen such that ρ was the overall expected percentage of causal sites in the genome across a concentration parameter α. For our simulations, we used ρ{0.0001,0.0003,0.001,0.003,0.01,0.03,0.05} and α{10,3,0.3} to represent different degrees of ‘clumpiness’ along the genome.

Genetic correlation between sex-stratified testosterone-related traits

Request a detailed protocol

LD Score regression [Bulik-Sullivan2015-tx] was used to generate genetic correlation estimates. The following command was used:

ldsc.py --rg <traits> --ref-ld-chr eur_ref_ld_chr 
  --w-ld-chr eur_w_ld_chr

where eur_*_ld_chr were downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/.

Residual height comparison with IGF-1

Request a detailed protocol

Height (adjusted for age and sex) and residualized log IGF-1 levels for unrelated White British individuals were plotted against each other, and visualized using geom_smooth.

Pathway diagrams

Request a detailed protocol

Diagrams were drawn using Adobe Illustrator and a Wacom graphics tablet.

PheWAS analysis

Request a detailed protocol

PheWAS were performed using the Oxford Brain Imaging Genetics (BIG) Server (Elliott et al., 2018).

Non-additivity tests

Request a detailed protocol

Residualized trait values were used as the outcome in all models. An ANOVA was performed between a model measuring the effect of genotype dosages versus a model with both genotype dosage effects and indicators for each rounded genotype. In this way, a large number of possible non-additive models are approximated with a single model. Analyses were performed in R 3.4 using lm.

Epistasis tests

Request a detailed protocol

We estimated that for hits with p<1e-20 we would have power to detect interaction components that are at least 10% the magnitude of a main effect (see Materials and methods). Thus, we tested all pairwise interactions among the independent lead SNPs with p<1e-20. Residualized trait values were used as the outcome in all models. An ANOVA was performed between a model measuring the effect of indicators for each rounded genotype (4 degrees of freedom) versus the interaction between the two sets of indicators (8 degrees of freedom). In this way, a large number of possible non-additive models are approximated with a test. Alternative models with dominant-only effect interactions with fewer degrees of freedom were also tested with similar results. Analyses were performed in R 3.4 using lm.

LD score regression for partitioning SNP-based heritability

Request a detailed protocol

We used partitioned LD Score regression (Finucane et al., 2015) to estimate the enrichment of individual tissues. We used the ldsc package and the updated BaselineLD v2.2 annotations with the following command:

ldsc.py --h2 <munged urate summary statistics> \ 
    --ref-ld-chr baselineLD.,<cell type annotations> \ 
    --overlap-annot --frqfile-chr 1000G_frq/1000G.mac5eur. \ 
    --w-ld-chr weights_hm3_no_hla/weights.

where <cell type annotations> were alternative either the default annotations for each of the ten cell type groups (Finucane et al., 2015) or modified versions which were filtered of any regulatory regions overlapping with the kidney cell type, using the following command:

ls 1000G_Phase3_cell_type_groups/*.bed | while read bed; do 
  intersectBed -a $bed -b 1000G_Phase3_cell_type_groups/7.bed -v > 
    1000G_Phase3_cell_type_groups_exclude_kidney/`basename $bed`; 
done

In this way, the cell type exclusive, non-kidney regulatory elements are used.

Data availability

Full raw summary statistics and relevant processed data tables are available on Figshare (https://doi.org/10.6084/m9.figshare.c.5304500.v1), or the lab website (http://web.stanford.edu/group/pritchardlab/dataArchive.html, direct link to google drive https://drive.google.com/drive/u/3/folders/10hCG_Wz8f25E6_sxw6sB8vDtS2OWUW9E).

The following data sets were generated
    1. Naqvi S
    (2021) figshare
    Supplementary Data for Sinnott-Armstrong and Naqvi.
    https://doi.org/10.6084/m9.figshare.c.5304500.v1
The following previously published data sets were used

References

    1. Baxter RC
    2. Martin JL
    3. Beniac VA
    (1989)
    High molecular weight insulin-like growth factor binding protein complex. Purification and properties of the acid-labile subunit from human serum
    The Journal of Biological Chemistry 264:11843–11848.
  1. Book
    1. Fisher R
    (1918)
    The Correlation Between Relatives on the Supposition of Mendelian Inheritance
    Transactions of the Royal Society of Edinburgh press.
    1. Herbst KL
    2. Bhasin S
    (2004) Testosterone action on skeletal muscle
    Current Opinion in Clinical Nutrition and Metabolic Care 7:271–277.
    https://doi.org/10.1097/00075197-200405000-00006
  2. Book
    1. Javorsky BR
    2. Aron DC
    3. Findling JW
    4. Tyrrell JB
    (2017)
    Chapter 4: Hypothalamus and Pituitary Gland
    New York: McGraw-Hill Medical.
    1. Köttgen A
    2. Albrecht E
    3. Teumer A
    4. Vitart V
    5. Krumsiek J
    6. Hundertmark C
    7. Pistis G
    8. Ruggiero D
    9. O'Seaghdha CM
    10. Haller T
    11. Yang Q
    12. Tanaka T
    13. Johnson AD
    14. Kutalik Z
    15. Smith AV
    16. Shi J
    17. Struchalin M
    18. Middelberg RP
    19. Brown MJ
    20. Gaffo AL
    21. Pirastu N
    22. Li G
    23. Hayward C
    24. Zemunik T
    25. Huffman J
    26. Yengo L
    27. Zhao JH
    28. Demirkan A
    29. Feitosa MF
    30. Liu X
    31. Malerba G
    32. Lopez LM
    33. van der Harst P
    34. Li X
    35. Kleber ME
    36. Hicks AA
    37. Nolte IM
    38. Johansson A
    39. Murgia F
    40. Wild SH
    41. Bakker SJ
    42. Peden JF
    43. Dehghan A
    44. Steri M
    45. Tenesa A
    46. Lagou V
    47. Salo P
    48. Mangino M
    49. Rose LM
    50. Lehtimäki T
    51. Woodward OM
    52. Okada Y
    53. Tin A
    54. Müller C
    55. Oldmeadow C
    56. Putku M
    57. Czamara D
    58. Kraft P
    59. Frogheri L
    60. Thun GA
    61. Grotevendt A
    62. Gislason GK
    63. Harris TB
    64. Launer LJ
    65. McArdle P
    66. Shuldiner AR
    67. Boerwinkle E
    68. Coresh J
    69. Schmidt H
    70. Schallert M
    71. Martin NG
    72. Montgomery GW
    73. Kubo M
    74. Nakamura Y
    75. Tanaka T
    76. Munroe PB
    77. Samani NJ
    78. Jacobs DR
    79. Liu K
    80. D'Adamo P
    81. Ulivi S
    82. Rotter JI
    83. Psaty BM
    84. Vollenweider P
    85. Waeber G
    86. Campbell S
    87. Devuyst O
    88. Navarro P
    89. Kolcic I
    90. Hastie N
    91. Balkau B
    92. Froguel P
    93. Esko T
    94. Salumets A
    95. Khaw KT
    96. Langenberg C
    97. Wareham NJ
    98. Isaacs A
    99. Kraja A
    100. Zhang Q
    101. Wild PS
    102. Scott RJ
    103. Holliday EG
    104. Org E
    105. Viigimaa M
    106. Bandinelli S
    107. Metter JE
    108. Lupo A
    109. Trabetti E
    110. Sorice R
    111. Döring A
    112. Lattka E
    113. Strauch K
    114. Theis F
    115. Waldenberger M
    116. Wichmann HE
    117. Davies G
    118. Gow AJ
    119. Bruinenberg M
    120. LifeLines Cohort Study
    121. Stolk RP
    122. Kooner JS
    123. Zhang W
    124. Winkelmann BR
    125. Boehm BO
    126. Lucae S
    127. Penninx BW
    128. Smit JH
    129. Curhan G
    130. Mudgal P
    131. Plenge RM
    132. Portas L
    133. Persico I
    134. Kirin M
    135. Wilson JF
    136. Mateo Leach I
    137. van Gilst WH
    138. Goel A
    139. Ongen H
    140. Hofman A
    141. Rivadeneira F
    142. Uitterlinden AG
    143. Imboden M
    144. von Eckardstein A
    145. Cucca F
    146. Nagaraja R
    147. Piras MG
    148. Nauck M
    149. Schurmann C
    150. Budde K
    151. Ernst F
    152. Farrington SM
    153. Theodoratou E
    154. Prokopenko I
    155. Stumvoll M
    156. Jula A
    157. Perola M
    158. Salomaa V
    159. Shin SY
    160. Spector TD
    161. Sala C
    162. Ridker PM
    163. Kähönen M
    164. Viikari J
    165. Hengstenberg C
    166. Nelson CP
    167. CARDIoGRAM Consortium
    168. DIAGRAM Consortium
    169. ICBP Consortium
    170. MAGIC Consortium
    171. Meschia JF
    172. Nalls MA
    173. Sharma P
    174. Singleton AB
    175. Kamatani N
    176. Zeller T
    177. Burnier M
    178. Attia J
    179. Laan M
    180. Klopp N
    181. Hillege HL
    182. Kloiber S
    183. Choi H
    184. Pirastu M
    185. Tore S
    186. Probst-Hensch NM
    187. Völzke H
    188. Gudnason V
    189. Parsa A
    190. Schmidt R
    191. Whitfield JB
    192. Fornage M
    193. Gasparini P
    194. Siscovick DS
    195. Polašek O
    196. Campbell H
    197. Rudan I
    198. Bouatia-Naji N
    199. Metspalu A
    200. Loos RJ
    201. van Duijn CM
    202. Borecki IB
    203. Ferrucci L
    204. Gambaro G
    205. Deary IJ
    206. Wolffenbuttel BH
    207. Chambers JC
    208. März W
    209. Pramstaller PP
    210. Snieder H
    211. Gyllensten U
    212. Wright AF
    213. Navis G
    214. Watkins H
    215. Witteman JC
    216. Sanna S
    217. Schipf S
    218. Dunlop MG
    219. Tönjes A
    220. Ripatti S
    221. Soranzo N
    222. Toniolo D
    223. Chasman DI
    224. Raitakari O
    225. Kao WH
    226. Ciullo M
    227. Fox CS
    228. Caulfield M
    229. Bochud M
    230. Gieger C
    (2013) Genome-wide association analyses identify 18 new loci associated with serum urate concentrations
    Nature Genetics 45:145–154.
    https://doi.org/10.1038/ng.2500
    1. Liu DJ
    2. Peloso GM
    3. Yu H
    4. Butterworth AS
    5. Wang X
    6. Mahajan A
    7. Saleheen D
    8. Emdin C
    9. Alam D
    10. Alves AC
    11. Amouyel P
    12. Di Angelantonio E
    13. Arveiler D
    14. Assimes TL
    15. Auer PL
    16. Baber U
    17. Ballantyne CM
    18. Bang LE
    19. Benn M
    20. Bis JC
    21. Boehnke M
    22. Boerwinkle E
    23. Bork-Jensen J
    24. Bottinger EP
    25. Brandslund I
    26. Brown M
    27. Busonero F
    28. Caulfield MJ
    29. Chambers JC
    30. Chasman DI
    31. Chen YE
    32. Chen YI
    33. Chowdhury R
    34. Christensen C
    35. Chu AY
    36. Connell JM
    37. Cucca F
    38. Cupples LA
    39. Damrauer SM
    40. Davies G
    41. Deary IJ
    42. Dedoussis G
    43. Denny JC
    44. Dominiczak A
    45. Dubé MP
    46. Ebeling T
    47. Eiriksdottir G
    48. Esko T
    49. Farmaki AE
    50. Feitosa MF
    51. Ferrario M
    52. Ferrieres J
    53. Ford I
    54. Fornage M
    55. Franks PW
    56. Frayling TM
    57. Frikke-Schmidt R
    58. Fritsche LG
    59. Frossard P
    60. Fuster V
    61. Ganesh SK
    62. Gao W
    63. Garcia ME
    64. Gieger C
    65. Giulianini F
    66. Goodarzi MO
    67. Grallert H
    68. Grarup N
    69. Groop L
    70. Grove ML
    71. Gudnason V
    72. Hansen T
    73. Harris TB
    74. Hayward C
    75. Hirschhorn JN
    76. Holmen OL
    77. Huffman J
    78. Huo Y
    79. Hveem K
    80. Jabeen S
    81. Jackson AU
    82. Jakobsdottir J
    83. Jarvelin MR
    84. Jensen GB
    85. Jørgensen ME
    86. Jukema JW
    87. Justesen JM
    88. Kamstrup PR
    89. Kanoni S
    90. Karpe F
    91. Kee F
    92. Khera AV
    93. Klarin D
    94. Koistinen HA
    95. Kooner JS
    96. Kooperberg C
    97. Kuulasmaa K
    98. Kuusisto J
    99. Laakso M
    100. Lakka T
    101. Langenberg C
    102. Langsted A
    103. Launer LJ
    104. Lauritzen T
    105. Liewald DCM
    106. Lin LA
    107. Linneberg A
    108. Loos RJF
    109. Lu Y
    110. Lu X
    111. Mägi R
    112. Malarstig A
    113. Manichaikul A
    114. Manning AK
    115. Mäntyselkä P
    116. Marouli E
    117. Masca NGD
    118. Maschio A
    119. Meigs JB
    120. Melander O
    121. Metspalu A
    122. Morris AP
    123. Morrison AC
    124. Mulas A
    125. Müller-Nurasyid M
    126. Munroe PB
    127. Neville MJ
    128. Nielsen JB
    129. Nielsen SF
    130. Nordestgaard BG
    131. Ordovas JM
    132. Mehran R
    133. O'Donnell CJ
    134. Orho-Melander M
    135. Molony CM
    136. Muntendam P
    137. Padmanabhan S
    138. Palmer CNA
    139. Pasko D
    140. Patel AP
    141. Pedersen O
    142. Perola M
    143. Peters A
    144. Pisinger C
    145. Pistis G
    146. Polasek O
    147. Poulter N
    148. Psaty BM
    149. Rader DJ
    150. Rasheed A
    151. Rauramaa R
    152. Reilly DF
    153. Reiner AP
    154. Renström F
    155. Rich SS
    156. Ridker PM
    157. Rioux JD
    158. Robertson NR
    159. Roden DM
    160. Rotter JI
    161. Rudan I
    162. Salomaa V
    163. Samani NJ
    164. Sanna S
    165. Sattar N
    166. Schmidt EM
    167. Scott RA
    168. Sever P
    169. Sevilla RS
    170. Shaffer CM
    171. Sim X
    172. Sivapalaratnam S
    173. Small KS
    174. Smith AV
    175. Smith BH
    176. Somayajula S
    177. Southam L
    178. Spector TD
    179. Speliotes EK
    180. Starr JM
    181. Stirrups KE
    182. Stitziel N
    183. Strauch K
    184. Stringham HM
    185. Surendran P
    186. Tada H
    187. Tall AR
    188. Tang H
    189. Tardif JC
    190. Taylor KD
    191. Trompet S
    192. Tsao PS
    193. Tuomilehto J
    194. Tybjaerg-Hansen A
    195. van Zuydam NR
    196. Varbo A
    197. Varga TV
    198. Virtamo J
    199. Waldenberger M
    200. Wang N
    201. Wareham NJ
    202. Warren HR
    203. Weeke PE
    204. Weinstock J
    205. Wessel J
    206. Wilson JG
    207. Wilson PWF
    208. Xu M
    209. Yaghootkar H
    210. Young R
    211. Zeggini E
    212. Zhang H
    213. Zheng NS
    214. Zhang W
    215. Zhang Y
    216. Zhou W
    217. Zhou Y
    218. Zoledziewska M
    219. Charge Diabetes Working Group
    220. EPIC-InterAct Consortium
    221. EPIC-CVD Consortium
    222. GOLD Consortium
    223. VA Million Veteran Program
    224. Howson JMM
    225. Danesh J
    226. McCarthy MI
    227. Cowan CA
    228. Abecasis G
    229. Deloukas P
    230. Musunuru K
    231. Willer CJ
    232. Kathiresan S
    (2017) Exome-wide association study of plasma lipids in >300,000 individuals
    Nature Genetics 49:1758–1766.
    https://doi.org/10.1038/ng.3977
  3. Book
    1. Lynch M
    2. Walsh B
    (1998)
    Genetics and Analysis of Quantitative Traits, 1
    Sunderland: Sinauer.
    1. McGuffin P
    2. Huckle P
    (1990)
    Simulation of mendelism revisited: the recessive gene for attending medical school
    American Journal of Human Genetics 46:994–999.
  4. Book
    1. Provine WB
    (2001)
    The Origins of Theoretical Population Genetics, 1971
    Chicago: University of Chicago Press.
    1. Tin A
    2. Marten J
    3. Halperin Kuhns VL
    4. Li Y
    5. Wuttke M
    6. Kirsten H
    7. Sieber KB
    8. Qiu C
    9. Gorski M
    10. Yu Z
    11. Giri A
    12. Sveinbjornsson G
    13. Li M
    14. Chu AY
    15. Hoppmann A
    16. O'Connor LJ
    17. Prins B
    18. Nutile T
    19. Noce D
    20. Akiyama M
    21. Cocca M
    22. Ghasemi S
    23. van der Most PJ
    24. Horn K
    25. Xu Y
    26. Fuchsberger C
    27. Sedaghat S
    28. Afaq S
    29. Amin N
    30. Ärnlöv J
    31. Bakker SJL
    32. Bansal N
    33. Baptista D
    34. Bergmann S
    35. Biggs ML
    36. Biino G
    37. Boerwinkle E
    38. Bottinger EP
    39. Boutin TS
    40. Brumat M
    41. Burkhardt R
    42. Campana E
    43. Campbell A
    44. Campbell H
    45. Carroll RJ
    46. Catamo E
    47. Chambers JC
    48. Ciullo M
    49. Concas MP
    50. Coresh J
    51. Corre T
    52. Cusi D
    53. Felicita SC
    54. de Borst MH
    55. De Grandi A
    56. de Mutsert R
    57. de Vries APJ
    58. Delgado G
    59. Demirkan A
    60. Devuyst O
    61. Dittrich K
    62. Eckardt KU
    63. Ehret G
    64. Endlich K
    65. Evans MK
    66. Gansevoort RT
    67. Gasparini P
    68. Giedraitis V
    69. Gieger C
    70. Girotto G
    71. Gögele M
    72. Gordon SD
    73. Gudbjartsson DF
    74. Gudnason V
    75. German Chronic Kidney Disease Study
    76. Haller T
    77. Hamet P
    78. Harris TB
    79. Hayward C
    80. Hicks AA
    81. Hofer E
    82. Holm H
    83. Huang W
    84. Hutri-Kähönen N
    85. Hwang SJ
    86. Ikram MA
    87. Lewis RM
    88. Ingelsson E
    89. Jakobsdottir J
    90. Jonsdottir I
    91. Jonsson H
    92. Joshi PK
    93. Josyula NS
    94. Jung B
    95. Kähönen M
    96. Kamatani Y
    97. Kanai M
    98. Kerr SM
    99. Kiess W
    100. Kleber ME
    101. Koenig W
    102. Kooner JS
    103. Körner A
    104. Kovacs P
    105. Krämer BK
    106. Kronenberg F
    107. Kubo M
    108. Kühnel B
    109. La Bianca M
    110. Lange LA
    111. Lehne B
    112. Lehtimäki T
    113. Lifelines Cohort Study
    114. Liu J
    115. Loeffler M
    116. Loos RJF
    117. Lyytikäinen LP
    118. Magi R
    119. Mahajan A
    120. Martin NG
    121. März W
    122. Mascalzoni D
    123. Matsuda K
    124. Meisinger C
    125. Meitinger T
    126. Metspalu A
    127. Milaneschi Y
    128. V. A. Million Veteran Program
    129. O'Donnell CJ
    130. Wilson OD
    131. Gaziano JM
    132. Mishra PP
    133. Mohlke KL
    134. Mononen N
    135. Montgomery GW
    136. Mook-Kanamori DO
    137. Müller-Nurasyid M
    138. Nadkarni GN
    139. Nalls MA
    140. Nauck M
    141. Nikus K
    142. Ning B
    143. Nolte IM
    144. Noordam R
    145. O'Connell JR
    146. Olafsson I
    147. Padmanabhan S
    148. Penninx B
    149. Perls T
    150. Peters A
    151. Pirastu M
    152. Pirastu N
    153. Pistis G
    154. Polasek O
    155. Ponte B
    156. Porteous DJ
    157. Poulain T
    158. Preuss MH
    159. Rabelink TJ
    160. Raffield LM
    161. Raitakari OT
    162. Rettig R
    163. Rheinberger M
    164. Rice KM
    165. Rizzi F
    166. Robino A
    167. Rudan I
    168. Krajcoviechova A
    169. Cifkova R
    170. Rueedi R
    171. Ruggiero D
    172. Ryan KA
    173. Saba Y
    174. Salvi E
    175. Schmidt H
    176. Schmidt R
    177. Shaffer CM
    178. Smith AV
    179. Smith BH
    180. Spracklen CN
    181. Strauch K
    182. Stumvoll M
    183. Sulem P
    184. Tajuddin SM
    185. Teren A
    186. Thiery J
    187. Thio CHL
    188. Thorsteinsdottir U
    189. Toniolo D
    190. Tönjes A
    191. Tremblay J
    192. Uitterlinden AG
    193. Vaccargiu S
    194. van der Harst P
    195. van Duijn CM
    196. Verweij N
    197. Völker U
    198. Vollenweider P
    199. Waeber G
    200. Waldenberger M
    201. Whitfield JB
    202. Wild SH
    203. Wilson JF
    204. Yang Q
    205. Zhang W
    206. Zonderman AB
    207. Bochud M
    208. Wilson JG
    209. Pendergrass SA
    210. Ho K
    211. Parsa A
    212. Pramstaller PP
    213. Psaty BM
    214. Böger CA
    215. Snieder H
    216. Butterworth AS
    217. Okada Y
    218. Edwards TL
    219. Stefansson K
    220. Susztak K
    221. Scholz M
    222. Heid IM
    223. Hung AM
    224. Teumer A
    225. Pattaro C
    226. Woodward OM
    227. Vitart V
    228. Köttgen A
    (2019) Target genes, variants, tissues and transcriptional pathways influencing human serum urate levels
    Nature Genetics 51:1459–1474.
    https://doi.org/10.1038/s41588-019-0504-x
    1. van Heemst D
    (2010)
    Insulin, IGF-1 and longevity
    Aging and Disease 1:147.
    1. Wood AR
    2. Esko T
    3. Yang J
    4. Vedantam S
    5. Pers TH
    6. Gustafsson S
    7. Chu AY
    8. Estrada K
    9. Luan J
    10. Kutalik Z
    11. Amin N
    12. Buchkovich ML
    13. Croteau-Chonka DC
    14. Day FR
    15. Duan Y
    16. Fall T
    17. Fehrmann R
    18. Ferreira T
    19. Jackson AU
    20. Karjalainen J
    21. Lo KS
    22. Locke AE
    23. Mägi R
    24. Mihailov E
    25. Porcu E
    26. Randall JC
    27. Scherag A
    28. Vinkhuyzen AA
    29. Westra HJ
    30. Winkler TW
    31. Workalemahu T
    32. Zhao JH
    33. Absher D
    34. Albrecht E
    35. Anderson D
    36. Baron J
    37. Beekman M
    38. Demirkan A
    39. Ehret GB
    40. Feenstra B
    41. Feitosa MF
    42. Fischer K
    43. Fraser RM
    44. Goel A
    45. Gong J
    46. Justice AE
    47. Kanoni S
    48. Kleber ME
    49. Kristiansson K
    50. Lim U
    51. Lotay V
    52. Lui JC
    53. Mangino M
    54. Mateo Leach I
    55. Medina-Gomez C
    56. Nalls MA
    57. Nyholt DR
    58. Palmer CD
    59. Pasko D
    60. Pechlivanis S
    61. Prokopenko I
    62. Ried JS
    63. Ripke S
    64. Shungin D
    65. Stancáková A
    66. Strawbridge RJ
    67. Sung YJ
    68. Tanaka T
    69. Teumer A
    70. Trompet S
    71. van der Laan SW
    72. van Setten J
    73. Van Vliet-Ostaptchouk JV
    74. Wang Z
    75. Yengo L
    76. Zhang W
    77. Afzal U
    78. Arnlöv J
    79. Arscott GM
    80. Bandinelli S
    81. Barrett A
    82. Bellis C
    83. Bennett AJ
    84. Berne C
    85. Blüher M
    86. Bolton JL
    87. Böttcher Y
    88. Boyd HA
    89. Bruinenberg M
    90. Buckley BM
    91. Buyske S
    92. Caspersen IH
    93. Chines PS
    94. Clarke R
    95. Claudi-Boehm S
    96. Cooper M
    97. Daw EW
    98. De Jong PA
    99. Deelen J
    100. Delgado G
    101. Denny JC
    102. Dhonukshe-Rutten R
    103. Dimitriou M
    104. Doney AS
    105. Dörr M
    106. Eklund N
    107. Eury E
    108. Folkersen L
    109. Garcia ME
    110. Geller F
    111. Giedraitis V
    112. Go AS
    113. Grallert H
    114. Grammer TB
    115. Gräßler J
    116. Grönberg H
    117. de Groot LC
    118. Groves CJ
    119. Haessler J
    120. Hall P
    121. Haller T
    122. Hallmans G
    123. Hannemann A
    124. Hartman CA
    125. Hassinen M
    126. Hayward C
    127. Heard-Costa NL
    128. Helmer Q
    129. Hemani G
    130. Henders AK
    131. Hillege HL
    132. Hlatky MA
    133. Hoffmann W
    134. Hoffmann P
    135. Holmen O
    136. Houwing-Duistermaat JJ
    137. Illig T
    138. Isaacs A
    139. James AL
    140. Jeff J
    141. Johansen B
    142. Johansson Å
    143. Jolley J
    144. Juliusdottir T
    145. Junttila J
    146. Kho AN
    147. Kinnunen L
    148. Klopp N
    149. Kocher T
    150. Kratzer W
    151. Lichtner P
    152. Lind L
    153. Lindström J
    154. Lobbens S
    155. Lorentzon M
    156. Lu Y
    157. Lyssenko V
    158. Magnusson PK
    159. Mahajan A
    160. Maillard M
    161. McArdle WL
    162. McKenzie CA
    163. McLachlan S
    164. McLaren PJ
    165. Menni C
    166. Merger S
    167. Milani L
    168. Moayyeri A
    169. Monda KL
    170. Morken MA
    171. Müller G
    172. Müller-Nurasyid M
    173. Musk AW
    174. Narisu N
    175. Nauck M
    176. Nolte IM
    177. Nöthen MM
    178. Oozageer L
    179. Pilz S
    180. Rayner NW
    181. Renstrom F
    182. Robertson NR
    183. Rose LM
    184. Roussel R
    185. Sanna S
    186. Scharnagl H
    187. Scholtens S
    188. Schumacher FR
    189. Schunkert H
    190. Scott RA
    191. Sehmi J
    192. Seufferlein T
    193. Shi J
    194. Silventoinen K
    195. Smit JH
    196. Smith AV
    197. Smolonska J
    198. Stanton AV
    199. Stirrups K
    200. Stott DJ
    201. Stringham HM
    202. Sundström J
    203. Swertz MA
    204. Syvänen AC
    205. Tayo BO
    206. Thorleifsson G
    207. Tyrer JP
    208. van Dijk S
    209. van Schoor NM
    210. van der Velde N
    211. van Heemst D
    212. van Oort FV
    213. Vermeulen SH
    214. Verweij N
    215. Vonk JM
    216. Waite LL
    217. Waldenberger M
    218. Wennauer R
    219. Wilkens LR
    220. Willenborg C
    221. Wilsgaard T
    222. Wojczynski MK
    223. Wong A
    224. Wright AF
    225. Zhang Q
    226. Arveiler D
    227. Bakker SJ
    228. Beilby J
    229. Bergman RN
    230. Bergmann S
    231. Biffar R
    232. Blangero J
    233. Boomsma DI
    234. Bornstein SR
    235. Bovet P
    236. Brambilla P
    237. Brown MJ
    238. Campbell H
    239. Caulfield MJ
    240. Chakravarti A
    241. Collins R
    242. Collins FS
    243. Crawford DC
    244. Cupples LA
    245. Danesh J
    246. de Faire U
    247. den Ruijter HM
    248. Erbel R
    249. Erdmann J
    250. Eriksson JG
    251. Farrall M
    252. Ferrannini E
    253. Ferrières J
    254. Ford I
    255. Forouhi NG
    256. Forrester T
    257. Gansevoort RT
    258. Gejman PV
    259. Gieger C
    260. Golay A
    261. Gottesman O
    262. Gudnason V
    263. Gyllensten U
    264. Haas DW
    265. Hall AS
    266. Harris TB
    267. Hattersley AT
    268. Heath AC
    269. Hengstenberg C
    270. Hicks AA
    271. Hindorff LA
    272. Hingorani AD
    273. Hofman A
    274. Hovingh GK
    275. Humphries SE
    276. Hunt SC
    277. Hypponen E
    278. Jacobs KB
    279. Jarvelin MR
    280. Jousilahti P
    281. Jula AM
    282. Kaprio J
    283. Kastelein JJ
    284. Kayser M
    285. Kee F
    286. Keinanen-Kiukaanniemi SM
    287. Kiemeney LA
    288. Kooner JS
    289. Kooperberg C
    290. Koskinen S
    291. Kovacs P
    292. Kraja AT
    293. Kumari M
    294. Kuusisto J
    295. Lakka TA
    296. Langenberg C
    297. Le Marchand L
    298. Lehtimäki T
    299. Lupoli S
    300. Madden PA
    301. Männistö S
    302. Manunta P
    303. Marette A
    304. Matise TC
    305. McKnight B
    306. Meitinger T
    307. Moll FL
    308. Montgomery GW
    309. Morris AD
    310. Morris AP
    311. Murray JC
    312. Nelis M
    313. Ohlsson C
    314. Oldehinkel AJ
    315. Ong KK
    316. Ouwehand WH
    317. Pasterkamp G
    318. Peters A
    319. Pramstaller PP
    320. Price JF
    321. Qi L
    322. Raitakari OT
    323. Rankinen T
    324. Rao DC
    325. Rice TK
    326. Ritchie M
    327. Rudan I
    328. Salomaa V
    329. Samani NJ
    330. Saramies J
    331. Sarzynski MA
    332. Schwarz PE
    333. Sebert S
    334. Sever P
    335. Shuldiner AR
    336. Sinisalo J
    337. Steinthorsdottir V
    338. Stolk RP
    339. Tardif JC
    340. Tönjes A
    341. Tremblay A
    342. Tremoli E
    343. Virtamo J
    344. Vohl MC
    345. Electronic Medical Records and Genomics (eMEMERGEGE) Consortium
    346. MIGen Consortium
    347. PAGEGE Consortium
    348. LifeLines Cohort Study
    349. Amouyel P
    350. Asselbergs FW
    351. Assimes TL
    352. Bochud M
    353. Boehm BO
    354. Boerwinkle E
    355. Bottinger EP
    356. Bouchard C
    357. Cauchi S
    358. Chambers JC
    359. Chanock SJ
    360. Cooper RS
    361. de Bakker PI
    362. Dedoussis G
    363. Ferrucci L
    364. Franks PW
    365. Froguel P
    366. Groop LC
    367. Haiman CA
    368. Hamsten A
    369. Hayes MG
    370. Hui J
    371. Hunter DJ
    372. Hveem K
    373. Jukema JW
    374. Kaplan RC
    375. Kivimaki M
    376. Kuh D
    377. Laakso M
    378. Liu Y
    379. Martin NG
    380. März W
    381. Melbye M
    382. Moebus S
    383. Munroe PB
    384. Njølstad I
    385. Oostra BA
    386. Palmer CN
    387. Pedersen NL
    388. Perola M
    389. Pérusse L
    390. Peters U
    391. Powell JE
    392. Power C
    393. Quertermous T
    394. Rauramaa R
    395. Reinmaa E
    396. Ridker PM
    397. Rivadeneira F
    398. Rotter JI
    399. Saaristo TE
    400. Saleheen D
    401. Schlessinger D
    402. Slagboom PE
    403. Snieder H
    404. Spector TD
    405. Strauch K
    406. Stumvoll M
    407. Tuomilehto J
    408. Uusitupa M
    409. van der Harst P
    410. Völzke H
    411. Walker M
    412. Wareham NJ
    413. Watkins H
    414. Wichmann HE
    415. Wilson JF
    416. Zanen P
    417. Deloukas P
    418. Heid IM
    419. Lindgren CM
    420. Mohlke KL
    421. Speliotes EK
    422. Thorsteinsdottir U
    423. Barroso I
    424. Fox CS
    425. North KE
    426. Strachan DP
    427. Beckmann JS
    428. Berndt SI
    429. Boehnke M
    430. Borecki IB
    431. McCarthy MI
    432. Metspalu A
    433. Stefansson K
    434. Uitterlinden AG
    435. van Duijn CM
    436. Franke L
    437. Willer CJ
    438. Price AL
    439. Lettre G
    440. Loos RJ
    441. Weedon MN
    442. Ingelsson E
    443. O'Connell JR
    444. Abecasis GR
    445. Chasman DI
    446. Goddard ME
    447. Visscher PM
    448. Hirschhorn JN
    449. Frayling TM
    (2014) Defining the role of common variation in the genomic and biological architecture of adult human height
    Nature Genetics 46:1173–1186.
    https://doi.org/10.1038/ng.3097

Decision letter

  1. Jonathan Flint
    Reviewing Editor; University of California, Los Angeles, United States
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States
  3. Vincent J Lynch
    Reviewer; University at Buffalo, United States
  4. Naomi Wray
    Reviewer
  5. Aravinda Chakravarti
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Thank you for submitting your article "GWAS of three molecular traits highlights core genes and pathways alongside a highly polygenic background" for consideration by eLife. Your article has been reviewed by four peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Patricia Wittkopp as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Vincent J Lynch (Reviewer #2); Naomi Wray (Reviewer #3); Aravinda Chakravarti (Reviewer #4).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.

Our expectation is that the authors will eventually carry out the additional experiments and report on how they affect the relevant conclusions either in a preprint on bioRxiv or medRxiv, or if appropriate, as a Research Advance in eLife, either of which would be linked to the original paper.

Summary:

Pritchard and colleagues use genome-wide association studies (GWAS) to identify genetic variants in the UK Biobank associated with three molecular traits-urate, IGF-1, and testosterone. Elegant and comprehensive analyses clearly demonstrate that the known biology of these traits explain many of the top hits of the GWAS, even when trait biology differs by sex, and that these core signals reside in a sea of polygenic variation.

The authors interpret their results within the framework of the "omnigenetic" model, namely that there is a difference between genetic effects from variants affecting core biological processes directly related to a trait, and a (usually much larger) number of loci that are not directly related to the trait.

Revisions for this paper:

1) The main concern is that the omnigenic hypothesis remains vague and untested. The authors do not address a crucial question: if the major hits are direct and core, are the other hits not core and therefore not important to unravel the biology of the trait?

One test of their hypothesis is that the many other significant associations are at genes that do affect the core genes. Their demonstration, at least for urate, is that the remaining non-core gene heritability is largely from the kidney. This is an important piece of the evidence. But they also need to demonstrate that these peripheral genes affect core gene expression.

They state that "genes that are expressed in trait-relevant cell types are referred to as "peripheral" genes, and can matter if they affect the expression of core genes." This is not demonstrated. The problem with the core versus peripheral definition is that it is imprecise. Is genetic variation in a transcription factor that directly regulates a core gene, core or peripheral?

While this concern does not lessen the value of the GWAS they perform, we'd like to see the authors either address this issue, or re-focus the manuscript so that the context for the GWAS is no longer to test the omnigenic hypothesis.

2) The authors should put their findings, and their interpretation, within the broader context of the literature on quantitative traits. As currently written, the manuscript provides a partial view of complex traits and disease. A better referenced Introduction and Discussion, and acknowledgement that the omnigenic model is consistent with long-published conceptualisations is needed.

For example, they conclude "these vignettes help to illustrate why many diseases are extraordinarily polygenic, as they are usually impacted by multiple biological processes that, like those considered here, are themselves highly polygenic." It needs to be pointed out that this conclusion is consistent with the thinking of the last 50 years.

Some acknowledgement needs to be made that many authors have reported that core biology can be recovered through top hits of a GWAS.

There is a long history of considering molecular phenotypes as endophenotypes or "intermediate" phenotypes. That literatures should be cited and the authors' results related to it. Without additional referencing Figure 9 is presented to the reader as if no others have conceptualised common disease as the endpoint of many contributing polygenic traits ("The point here is that when multiple risk factors-each of which is polygenic-contribute to any given disease, the disease endpoint absorbs the polygenic basis for all of the risk factors together."). For example, Figures 1 and 2 of Gottesman and Gould are conceptually very similar to Figure 9 (https://doi.org/10.1176/appi.ajp.160.4.636). The review of milk coagulation traits in 2012 (http://dx.doi.org/ 10.3168/jds.2012-5507) provides a Figure (#5) conceptually similar to Figure 9. There are other examples (e.g., PLoS Genet 6(9): e1001139. doi:10.1371/journal.pgen.1001139).

The first paragraph of the Discussion states "We showed that unlike most disease traits, these three biomolecules have clear enrichment of genome-wide significant signals in core genes and pathways." The statement reads as if it were a novel finding, but it is expected that as a DNA-trait relationship gets closer SNP associations will be biologically more obvious (and indeed the terms core and peripheral genes have been used long before the advent of the omnigenic model).

Revisions expected in follow-up work:

1) Authors should address the issue of whether their findings can address key elements the omnigenic hypothesis, or re-focus the manuscript so that the context for the GWAS is no longer so tied to the omnigenic hypothesis.

2) Put their findings, and their interpretation, within the broader context of the literature on quantitative traits.

https://doi.org/10.7554/eLife.58615.sa1

Author response

Revisions for this paper:

1) The main concern is that the omnigenic hypothesis remains vague and untested. The authors do not address a crucial question: if the major hits are direct and core, are the other hits not core and therefore not important to unravel the biology of the trait?

We respectfully disagree that the hypothesis remains vague, and in that regard, point to Liu et al., 2019, which lays out our model in detail, elaborating on the verbal model in Boyle et al., 2017. Liu et al. sought to develop a model to understand why it is that so many variants, spread so widely across the genome, can be responsible for the heritability of a typical complex trait. Very briefly, Liu et al. proposed that most heritability is due to variants with cis-effects on peripheral genes. These in turn perturb the regulation of a smaller class of core genes via trans regulatory networks; genetic effects on traits are mediated via core genes. Thus, peripheral genes are important to unraveling the biology of a trait because they explain the trans-regulatory context within which core genes sit.

That said, we acknowledge that at this time the model remains largely conceptual. There are no traits that are understood in sufficient detail to fully evaluate the model. For most traits the core genes are unknown; further, at this time we do not know detailed trans-regulatory networks in any cell type. Nonetheless, we would argue that there is value in articulating specific models for how genetic variation impacts traits. In population genetics and quantitative genetics there is a long history of theoretical models that preceded the relevant data, and – whether these have turned out to be right or wrong – they have been invaluable in guiding and informing future research questions.

In this paper we elucidate one major component of the model, namely the identities and roles of core genes, for three example traits. For most disease traits there is only very incomplete knowledge of likely core genes; here we use these three molecular traits to (1) Provide examples of core genes, and to show huge enrichment of signal around those, and

(2) To show that in these examples, variation at core genes explains only a small fraction of the heritability.

By doing so, our paper starts to connect the dots for certain key aspects of the model that have been difficult to verify in more complicated traits. Clarifying the nature and roles of core genes in specific examples is an important step, as some commentators, including one of the reviewers of this paper, have argued against the concept of core genes. We acknowledge, however, that our paper does not touch on the network part of the model (as now mentioned in the Introduction) – this remains for future work. In summary: we have reframed the manuscript, including the Introduction, to make clear that this study focuses on specific aspects of the omnigenic model that have not, to date, been rigorously tested with real data, while not yet addressing all components of the model.

The authors do not address a crucial question: if the major hits are direct and core, are the other hits not core and therefore not important to unravel the biology of the trait?

We have added a consideration of this and related issues to the Discussion. In brief, if tens of thousands of variants, acting through at least thousands of genes, have nonzero effects, then presumably most of these are not going to be useful for understanding the biological mechanism of a disease.

One test of their hypothesis is that the many other significant associations are at genes that do affect the core genes. Their demonstration, at least for urate, is that the remaining non-core gene heritability is largely from the kidney. This is an important piece of the evidence. But they also need to demonstrate that these peripheral genes affect core gene expression.

We agree that as an additional step to evaluate the omnigenic model, it will be important to study regulation of core gene expression by delineating cellular regulatory networks. While an experimental demonstration of such regulatory connections is beyond the scope of this paper (as we now note in the Introduction), we have analyzed cis- vs. trans-heritability of expression of the sets of curated core genes for each of the model traits in this paper. As a set, we have found that core genes do not have significantly higher expression cis-heritability than all other genes (Figure 8—figure supplement 14). This is consistent with the idea that much of core gene expression is determined in trans, presumably through regulation by peripheral genes. A more comprehensive analysis of this would require enumerating cis vs. trans heritabilities of genes in a variety of tissues, which would be an immense project outside the scope of the current work.

They state that "genes that are expressed in trait-relevant cell types are referred to as "peripheral" genes, and can matter if they affect the expression of core genes." This is not demonstrated. The problem with the core versus peripheral definition is that it is imprecise. Is genetic variation in a transcription factor that directly regulates a core gene, core or peripheral?

In this study, we have defined core genes in line with the definition from Liu et al., 2019. In that light, a transcription factor that regulates other core genes would be a peripheral gene, or perhaps a master regulator if it simultaneously regulates multiple core genes. However, as defined by Liu et al., we consider signaling receptors, such as the androgen receptor for testosterone, to be core even though they are transcription factors, since they directly receive inputs from outside the cell. Transcription factors that act within a cellular regulatory network , rather than directly receiving inputs from outside the cell, would be considered peripheral. Nevertheless, we admit that this is a simplified, conceptual model, such that not all genes will fit neatly into these definitions. We have made these points more clearly in the Introduction.

While this concern does not lessen the value of the GWAS they perform, we'd like to see the authors either address this issue, or re-focus the manuscript so that the context for the GWAS is no longer to test the omnigenic hypothesis.

As noted above, the contribution of this paper to the omnigenic model is in: (1) illustrating the role of core genes for specific traits; (2) showing that these contribute only a small fraction of the heritability; and (3) showing that even these seemingly simpler traits are affected by on the order of 10,000 genes.

2) The authors should put their findings, and their interpretation, within the broader context of the literature on quantitative traits. As currently written, the manuscript provides a partial view of complex traits and disease. A better referenced Introduction and Discussion, and acknowledgement that the omnigenic model is consistent with long-published conceptualisations is needed.

We disagree with the vague claim that “the omnigenic model is consistent with long-published conceptualisations” as well as the implication below that it is “consistent with the thinking of the last 50 years” and believe they represent a misunderstanding of the central question addressed in this manuscript.

Prior to ~2006, researchers in human genetics had almost no idea how many loci might underlie variation in complex traits, and many expected that this might be on the order of tens of loci, and that some loci would have large effect sizes (see references in our Introduction). Only in the last few years has it become clear how many causal variants there are, and that most of the genome contributes to heritability. While we appreciate that theoretical models in quantitative genetics had considered the possibility that many loci might contribute to variation in a trait, there was no empirical evidence that they applied to humans. Most importantly, those studies were not focused on – or even interested in – the mechanisms that led to polygenic inheritance.

In that regard, these new findings pose a new question that would not have seemed relevant prior to the GWAS era: From a mechanistic point of view, how should we interpret the observation that so many variants (and by extension so many genes) contribute to any given trait, and that the lead variants contribute such a small fraction of heritability? The central element of our model rests on the roles of cis and trans-acting expression QTLs. The observation that most expression variance is due to many small trans effects has emerged relatively recently – roughly in the last dozen years – and those results have, arguably, been underappreciated in the field. Thus, beyond the much more generic claim that complex trait variation may be highly polygenic, we’re not aware that anyone else has formulated the specific model previously articulated by Boyle et al., 2017, and Liu et al., 2019.

For example, they conclude "these vignettes help to illustrate why many diseases are extraordinarily polygenic, as they are usually impacted by multiple biological processes that, like those considered here, are themselves highly polygenic." It needs to be pointed out that this conclusion is consistent with the thinking of the last 50 years.

The specific comment here refers to a section in the Discussion noting that when there are multiple biological processes contributing to a disease these can further increase polygenicity in addition to the main themes that we have discussed. This is not really part of the main omnigenic model, which proposes a key role for intracellular networks. We mentioned the multiple-processes point as it presumably also contributes to the polygenicity of many disease endpoints; that said it seems less likely to explain simpler traits such as those studied here, especially urate where the heritability comes almost entirely from the kidneys. Thus, this study provides empirical evidence for one model of why disease polygenicity is so high.

Given that the presence of this section in the Discussion seems to have been distracting, we removed Figure 9 and rewrote the corresponding paragraph to make clearer that we think of this as a complementary model rather than part of the same model.

Some acknowledgement needs to be made that many authors have reported that core biology can be recovered through top hits of a GWAS.

We had reviewed the past literature on core genes in our previous papers on this topic (Boyle et al., 2017, Liu et al., 2019). It seems outside scope to cover that at great length again, but we now point readers to some of this literature with the following:

“We do now know various examples of core genes or master regulators for specific traits (e.g., Sekar et al., 2016, Small et al., 2011, Small et al., 2018), but there are few traits where we understand the roles of more than a few of the lead genes. Among the clearest examples in which a whole suite of core genes have been identified are for plasma lipid levels (e.g., Liu et al., 2017, Lu et al., 2017, Hoffmann et al., 2018), reviewed by Dron et al., 2016, Liu et al., 2019; and for inflammatory bowel disease (de Lange et al., 2017).”.

There is a long history of considering molecular phenotypes as endophenotypes or "intermediate" phenotypes. That literatures should be cited and the authors' results related to it. Without additional referencing Figure 9 is presented to the reader as if no others have conceptualised common disease as the endpoint of many contributing polygenic traits ("The point here is that when multiple risk factors-each of which is polygenic-contribute to any given disease, the disease endpoint absorbs the polygenic basis for all of the risk factors together."). For example, Figures 1 and 2 of Gottesman and Gould are conceptually very similar to Figure 9 (https://doi.org/10.1176/appi.ajp.160.4.636). The review of milk coagulation traits in 2012 (http://dx.doi.org/ 10.3168/jds.2012-5507) provides a Figure (#5) conceptually similar to Figure 9. There are other examples (e.g., PLoS Genet 6(9): e1001139. doi:10.1371/journal.pgen.1001139).

Thank you for these additional citations. We note that our previous version did in fact reference previous papers in this context, including Udler, 2019, who had a similar figure, and Turkheimer, 2000, who has written extensively on this point in the context of behavioral genetics. That said, as noted above, we have removed Figure 9 as this was somewhat extraneous to our main points, shortened the corresponding section in the Discussion, and added these additional references.

The first paragraph of the Discussion states "We showed that unlike most disease traits, these three biomolecules have clear enrichment of genome-wide significant signals in core genes and pathways." The statement reads as if it were a novel finding, but it is expected that as a DNA-trait relationship gets closer SNP associations will be biologically more obvious (and indeed the terms core and peripheral genes have been used long before the advent of the omnigenic model).

Here, use of the word “expected” suggests that perhaps the reviewers are referring to their expectations rather than specific prior data analyses that demonstrate this point.

With regard to the three biomolecules considered here, prior GWAS analysis showed mixed evidence for enrichment in the relevant pathways. While the urate transporters were well-known to be enriched for GWAS hits, enrichment in the urate synthesis pathway was not demonstrated, and no such pathway enrichments were demonstrated for IGF-1. Similarly, Ruth et al., 2020, did not focus on testosterone biology as the focus was mainly on sex differences in genetic effects and their relationship to cardiometabolic traits, and previous papers (e.g. Ohlsson et al., 2011, Ruth et al., 2015, and Prescott et al., 2012) focus almost exclusively on the strong association with SHBG variants in males and have no genome-wide significant associations in females (though Prescott et al., 2012 reports a single CYP4A1 subthreshold association in passing as well).

Revisions expected in follow-up work:

1) Authors should address the issue of whether their findings can address key elements the omnigenic hypothesis, or re-focus the manuscript so that the context for the GWAS is no longer so tied to the omnigenic hypothesis.

We have clarified that in this paper we elucidate one major component of the model, namely the identities and roles of core genes. For most disease traits there is only very incomplete knowledge of likely core genes; here we use these three molecular traits to (1) Provide examples of core genes, and to show huge enrichment of signal around those, and

(2) To show that in these examples, core genes explain only a small fraction of the heritability. We further clarified that the paper does not address the network component of the model.

(2) Put their findings, and their interpretation, within the broader context of the literature on quantitative traits.

We have expanded the Introduction and Discussion to provide more context on the history of studies of complex traits.

https://doi.org/10.7554/eLife.58615.sa2

Article and author information

Author details

  1. Nasa Sinnott-Armstrong

    Department of Genetics, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Contributed equally with
    Sahin Naqvi
    For correspondence
    nasa@stanford.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4490-0601
  2. Sahin Naqvi

    1. Department of Genetics, Stanford University, Stanford, United States
    2. Department of Chemical and Systems Biology, Stanford University, Stanford, United States
    Contribution
    Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Nasa Sinnott-Armstrong
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2635-7967
  3. Manuel Rivas

    Department of Biomedical Data Sciences, Stanford University, Stanford, United States
    Contribution
    Resources, Data curation, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Jonathan K Pritchard

    1. Department of Genetics, Stanford University, Stanford, United States
    2. Department of Biology, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    pritch@stanford.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8828-5236

Funding

American Society for Engineering Education (NDSEG)

  • Nasa Sinnott-Armstrong

Stanford University (Stanford Graduate Fellowship)

  • Nasa Sinnott-Armstrong

National Human Genome Research Institute (HG008140)

  • Jonathan K Pritchard

National Human Genome Research Institute (HG009431)

  • Jonathan K Pritchard

Helen Hay Whitney Foundation (HHWF 2020)

  • Sahin Naqvi

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank members of the Pritchard, Page, Przeworski, Sella, and Bassik labs, as well as Ipsita Agarwal, Evan Boyle, Eric Fauman, Jake Freimer, Rebecca Harris, Yang Li, Xuanyao Liu, Iain Mathieson, Molly Przeworski, Guy Sella, Yuval Simons, and Jeff Spence for helpful discussions or comments; and the UK Biobank and its participants for making this project possible, which we accessed through UK Biobank application number 24983. This work was supported by NIH grants HG008140 and HG009431 (to JKP), a Stanford Graduate Fellowship (to NS-A), a National Defense Science and Engineering Grant (to NS-A), and a Helen Hay Whitney Fellowship (to SN).

Ethics

Human subjects: This research has been conducted using the UK Biobank Resource under Application Number 24983, "Generating effective therapeutic hypotheses from genomic and hospital linkage data" (criteriahttp://www.ukbiobank.ac.uk/wp-content/uploads/2017/06/24983-Dr-Manuel-Rivas.pdf). Based on the information provided in Protocol 44532 the Stanford IRB has determined that the research does not involve human subjects as defined in 45 CFR 46.102(f) or 21 CFR 50.3(g). All participants of UK Biobank provided written informed consent (more information is available at https://www.ukbiobank.ac.uk/2018/02/gdpr/).

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Jonathan Flint, University of California, Los Angeles, United States

Reviewers

  1. Vincent J Lynch, University at Buffalo, United States
  2. Naomi Wray
  3. Aravinda Chakravarti

Publication history

  1. Received: May 5, 2020
  2. Accepted: January 18, 2021
  3. Version of Record published: February 15, 2021 (version 1)

Copyright

© 2021, Sinnott-Armstrong et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 5,227
    Page views
  • 464
    Downloads
  • 9
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Evolutionary Biology
    2. Genetics and Genomics
    Gabriela Santos-Rodriguez et al.
    Research Article

    Many primate genes produce circular RNAs (circRNAs). However, the extent of circRNA conservation between closely related species remains unclear. By comparing tissue-specific transcriptomes across over 70 million years of primate evolution, we identify that within 3 million years circRNA expression profiles diverged such that they are more related to species identity than organ type. However, our analysis also revealed a subset of circRNAs with conserved neural expression