The Polygenic and Monogenic Basis of Blood Traits and Diseases

1A major aspiration in human genetics is to understand how genetic variation impacts complex traits and diseases. Recent genome-wide association studies (GWAS) have identified thousands of genetic variants associated with complex phenotypes and provided insights into their genetic architecture. This has led to the recognition that complex trait heritability is polygenic, resulting from the cumulative effects of many genetic loci throughout the genome, each of modest effect size.

2Hematopoiesis is a valuable paradigm for studying complex trait genetic architecture, since blood cell phenotypes are commonly measured in large population-based studies and the production of blood cells is a highly regulated, hierarchical, and intrinsic process that can be readily studied. While there have been advances in understanding genetic loci associated with blood cell production, the spectrum of human genetic variation impacting hematopoiesis remains incompletely defined.

3Most variants contributing to complex trait heritability are noncoding and located in genomic regulatory regions within relevant cell types. The availability of epigenomic and transcriptomic profiles for hematopoietic stem and progenitor and lineage-committed cells enable mechanistic dissection of the roles that different classes of genes have in hematopoiesis. Prior studies of blood cell traits have suggested that master transcription factors (TFs) may be impacted by genetic variation, and it is likely that further studies may uncover additional roles for, and variation of, key hematopoietic regulators. Another priority is to advance understanding of network connectivity between trait-associated genes and variants, and this understanding can be informed by theoretical models. Recently, an “omnigenic” model has been proposed in which two types of genes (“core” versus “peripheral”) differentially contribute to complex trait heritability. However, the extent to which the omnigenic model applies to various complex traits and diseases remains unclear and controversial.

4Finally, although rare variants with large effects generally do not individually contribute substantially to overall complex trait variance, they can often highlight important biologic mechanisms and contribute to rare hematologic disorders, many of which are characterized by variable penetrance or expressivity. In addition, polygenic contributions of many variants with small effects can yield disease risk odds ratios comparable to or larger than that of known monogenic variants. Therefore, large population-based datasets can help to both reclassify the pathogenicity and penetrance of disease-associated variants, as well as understand the contribution of polygenic variation to the risk of blood diseases or as modifiers of rare variants that contribute to presumed monogenic blood disorders.

5We leveraged the power of the UK Biobank cohort to perform a genome-wide discovery analysis in N = 408,112 participants of European ancestry, investigating 29 blood cell phenotypes (Table S1). In parallel, we also performed tests for genetic associations with a subset of 15 phenotypes available in an additional 154,973 European ancestry participants from the Blood Cell Consortium (BCX) (Figure 1A, Table S2). A separate analysis of non-European participants is reported in a companion paper. Overall, this discovery effort identified 16,643 autosomal and 257 X-linked conditionally independent (Method Details) trait-variant associations from the first stage discovery and an additional 141 from the BCX meta-analysis (Tables S3 and S4). The 16,900 associations were assigned to 7,122 genomic loci (5,106 not described before) using a linkage disequilibrium (LD) clumping approach. Each locus was represented by a unique tag variant (between-tag pairwise LD r2 ≤ 0.8), and for simplicity, throughout the paper we use the term “sentinel variant” to refer to either a clump tag variant or a trait-specific conditionally independent signal. Overall, we nearly tripled the number of loci reported prior to this study. We assessed replication rates across three exemplar phenotypes (platelet count [PLT], lymphocyte count, and red blood cell count) for 210 variants on chromosome 1 in the Million Veteran Program (MVP, N = 271,280). We found that nearly all of them had directionally concordant effect size estimates (Pearson’s R2 = 0.94; Figure S1A), and 196 (93%) variants replicated at a nominal significance threshold (p < 0.05). The non-replicating ones exhibited similar effect sizes as in the discovery cohort but lacked power due to MVP having less than half the sample size of the discovery cohort (Figure S1A, zoom-in panel). Using a Bayesian method that accounts for multiple independent signals (Method Details, Figure 1B), we fine-mapped 3,100 (19% of 16,643 autosomal) associations to a single putative causative variant (> 95% posterior probability [PPFM]) (Table S5), and more than half of the associated signals (n = 9,149, 55%) to fewer than 10 variants (Figure 1C). As expected, rare signals are more likely to be fine-mapped to smaller credible sets (Figure 1C). We assigned sentinels to genes using a stringent variant effect predictor (VEP) worst-consequence annotation to obtain a distribution of functional categories. Overall, 8,866 sentinels (83%) were annotated to a gene using this approach, of which 69% were intronic, 24% were in regulatory regions, and 7% were in protein-coding regions (5.5% non-synonymous and 1.5% synonymous; Figure 1D). The credible set size distribution (number of variants per credible set) was consistent across traits (Figure 1E).

6Hematopoiesis is a finely tuned process involving coordinated expression of hundreds of genes, and it is likely that a subset of the variants associated with peripheral blood cell counts and indices acts upon master regulators of this process. To identify whether genes discovered by GWAS identify networks of coregulated genes, we accessed a published coexpression network of 7,509 protein-coding genes expressed in whole blood (Figures 2A and 2B; Method Details). Under the stringent VEP worst-consequence criteria used earlier, 25% of network genes (n = 1,874 genes) were annotated to a GWAS signal. A more permissive VEP any-consequence criteria annotated an additional 2.5% (27.5%, n = 2,070) genes. When all genes in the fine-mapping regions were considered (± 250-kb window), 78% of network genes could be linked to a GWAS locus, and 88% of sentinels were in proximity (< 250 kb) to a network gene, suggesting that genes linked to association signals are likely to be coregulated. Where possible, gene assignments were also validated using colocalization with (expression quantitative trait loci) cis-eQTLs derived from six trait-matched blood cell types (platelets n = 424; CD19+ B cells, CD8+ T cells, CD4+ T cells and CD15+ neutrophils n = 300; CD14+ monocytes n = 1,490). Across 667 colocalizing cis-eQTLs, eGenes matched VEP worst-consequence genes in 65% of the cases (Figure S1B) and were contained in fine-mapping regions in 97% of the cases.

7Biological networks are organized hierarchically. The recently proposed “omnigenic” model postulates that a small number of genes at the center (or “core”) of the network are directly implicated in diseases or phenotypes of interest, but the variants in these genes contribute only a small proportion of the overall trait heritability. Most of the trait heritability is attributable to a much larger number of “peripheral” gene variants with small effect sizes that contribute to subtler physiological perturbations of phenotypes through trans-regulatory effects on core genes. We thus sought to empirically test the main assumptions of the omnigenic model, compared to a more continuous “infinitesimal” model of disease heritability in order to inform its utility for disease gene discovery. We accessed a manually curated list of genes causative for stem cell and myeloid disorders (SMD, 206 genes); bleeding, thrombotic, and platelet disorders (BPD, 104 genes); and bone-marrow failure (BMF) syndromes (80 genes; Table S7). GWAS loci for blood cell indices tended to be strongly enriched in and near Mendelian blood disorder genes (by 2.1-fold, p = 1.9×10−22), a phenomenon already described for many complex traits. We then asked whether these Mendelian genes had properties expected of core genes.

8A first assumption of the model is that core genes are strongly enriched at the center of biological networks (Figure 2C). Overall, we observed strong enrichments of both GWAS (fold enrichment [FE] = 1.86, permutation p < 10−4) and Mendelian (i.e., core, FE = 3.86, p < 10−4) genes in the full blood coexpression network compared to permuted sets of protein-coding genes of similar size (Figure 2D; Table S6). Importantly, Mendelian genes had more connections in the coexpression network compared to other (non-Mendelian) genes, consistent with a centrality scenario (valid for coexpression cut-offs at 0.4–0.8, p ranging from 4×10−4 to 0.02, Wilcoxon test; Figure 2E). Finally, the expression of Mendelian genes was more correlated with other Mendelian genes (median coexpression coefficient = 0.11) than random sets of genes (median = 0.095, p = 0.007 permutation test). A second assumption is that variants assigned to core genes have larger effect sizes than peripheral genes. When compared to variants of comparable minor allele frequency (MAF) assigned to other genes, variants assigned to Mendelian genes (including previously unreported ones) showed significantly higher absolute effect sizes across all functional categories tested (fold change ranging from 1.87- to 2.73-fold increase; Figure S1C). Third, core genes should be more phenotype-specific as opposed to peripheral associations which act as regulators and could be shared across different phenotypes. We show by quartile-quartile (Q-Q) plot enrichments that this pattern holds true for Mendelian versus peripheral blood traits associations in a selection of eight non-blood related traits (Figure S2).

9The model also predicts that peripheral variants explain a large proportion of trait heritability through trans-regulation of core genes. To test this hypothesis, we accessed a large set of recently reported blood trans-eQTLs. Mendelian genes were strongly enriched as targets of trans-eQTLs, compared to other GWAS genes (2.11-fold, Wilcoxon test, p = 4.7×10−5), after matching for expression levels and trans-eQTL Z-scores to account for differences in detection power, with the caveat that there may be other unaccounted factors involved. At a correlation cut-off of 0.8, a coexpression subnetwork of 26 GWAS-associated genes was centered on three known Mendelian genes causative for spherocytosis (SLC4A1, EPB42) and congenital anemia (KLF1; Figure S3A). Interestingly, these factors all play key roles in red blood cell cytoskeleton formation, a process regulated by KLF1. Another example is a subnetwork containing known platelet specific genes GP9, GP1BB, and ITGA2B, and eight other strongly coexpressed genes (Figure 2F). All of these genes are trans-regulated by the ARHGEF3 gene, a known master regulator of megakaryopoiesis. While these results are broadly compatible with expectations of the omnigenic model, first- and second-degree coexpression network neighbors of Mendelian genes were also enriched for GWAS associations (p < 1×10−3, permutation test) and thus had properties attributable to both core and peripheral genes. This indicates either that these loci may fit a more continuous infinitesimal model, or that our current proposed set of core genes is incomplete.

10We next sought to delineate relevant cell states impacted by core and peripheral gene networks. To this end, we integrated all fine-mapped (FM) variants (PPFM > 0.1%) with chromatin accessibility profiles (ATAC-seq) of 18 human hematopoietic progenitor populations. First, we noted that FM variants falling within hematopoietic open chromatin were strongly enriched in gene targets (assigned by VEP worst consequence) compared to non-accessible variants (OR = 1.4, Fisher’s p < 2.2×10−16), consistent with variants acting via trans-regulation of genes in hematopoietic cell states. Next, we used g-chromVAR, a high-resolution cell type enrichment method, to determine the hematopoietic populations most enriched for chromatin accessibility containing FM variants for 22 blood cell traits, including 6 new traits compared to a previous study in a smaller subset of the UK Biobank. There were 43 lineage-specific enrichments surpassing experiment-wide significance (corrected for 18 cell types ×22 traits, p < 1.26×10−4) (Figure 3A), of which 20 were new, including novel enrichments in granulocyte-monocyte progenitor (GMP) cell subsets for variants regulating monocyte, eosinophil, and neutrophil counts.

11We then wondered whether certain trait-cell type enrichments would strengthen when restricting to core genes for corresponding blood diseases. To this end, we calculated enrichments for four platelet traits, considering only the variants mapping to core genes for BPD. Whereas the gene-agnostic analysis produced significant enrichments in both megakaryocytes (n = 4) and its less differentiated myeloid precursors (n = 7), the core-gene restricted approach led to a strong signal for megakaryocytes (n = 4) but a lack of enrichment in any other population (Figure 3B). This suggests that in addition to their roles in Mendelian disease, core genes are also enriched for trans-regulatory variants acting specifically in their causal cell type.

12Next, we sought to predict nucleotide-specific effects of variants on chromatin accessibility. We used deltaSVM, a support-vector machine classifier, to train genomic sequence features of the ATAC-seq from 18 hematopoietic cell populations (Figures S3B and S3C), and then applied the model to predict the allele-specific, cell type-specific impact of FM variants on chromatin accessibility. Out of 215,694 variants with PPFM > 0.001 for one or more blood traits, we identified 22,152 variants with an absolute deltaSVM score above the 99th percentile for at least one hematopoietic cell type. Absolute deltaSVM score was negatively associated with MAF (linear regression p < 2.2×10−16) and positively associated with FM PPFM (linear regression p = 1.0×10−3) (Figure S3D). Variants assigned to a gene by VEP worst consequence had stronger predicted effects on chromatin accessibility compared to intergenic variants (Student’s t test, p = 1.2×10−3); however, there was no significant difference in deltaSVM between variants assigned to “core” versus “peripheral” genes, suggesting that variant-mediated modulation of hematopoietic transcription occurs across the entire gene regulatory network rather than disproportionately impacting core genes.

13To further characterize the regulatory effects of these variants, we predicted the potential for FM variants to disrupt 426 human TF motifs. Across motif-disrupting variants, alternative alleles predicted to increase chromatin accessibility (deltaSVM score > 99th percentile) had a significantly higher motif matching score compared to the reference allele (Method Details). The reverse was also true, indicating that deltaSVM scores track with the potential for variants to break or create TF motifs (Figure 3C). Moreover, this trend was cell type specific, as evidenced by the fact that variants affecting lineage-determining TFs had a higher deltaSVM score within lineage-specific cell types, such as GATA1-disrupting variants within erythroid progenitors, compared to other hematopoietic populations (Figures 3D–3F). We next sought to integrate these functional annotations in order to gain novel insights into biologically relevant variants. For example, variant rs72928038, previously identified in a locus associated with lymphocyte count, was fine-mapped here as the likely causal variant (PPFM = 0.78), with the minor allele A (MAF = 18%) corresponding to decreased lymphocyte count. The variant maps to intron 1 of lymphoid TF BACH2 colocalizes with a H3K27ac histone QTL in CD4+ T cells and has high chromatin accessibility in the CD4+ and CD8+ T lymphoid populations. Interestingly, the lymphocyte count-decreasing minor allele has strongly negative deltaSVM scores (i.e., predicted to decrease chromatin accessibility) in the lymphoid lineage and is predicted to disrupt multiple TF motifs at the BACH2 locus, including those with known roles in lymphocyte development, such as STAT3 and ETS1 (Figures 3G and 3H). This variant has been previously implicated in risk for several autoimmune conditions including rheumatoid arthritis and vitiligo. These lines of evidence suggest that rs72928038 may affect lymphocyte count by altering the binding of specific lymphoid TFs within T cell progenitors. Altogether, our functional characterization of non-coding blood trait variants highlights the value of incorporating lineage-specific chromatin accessibility profiles and motif disruption analyses to nominate high-confidence mechanisms.

14The large sample size and dense imputation in this study gave us unprecedented statistical power to discover variants with low MAF and to assess their impact on human disease. First, we identified 574 rare (minor allele count [MAC] > 20, MAF < 1%) blood trait variants which were either conditionally independent lead variants and/or strongly fine-mapped (PPFM>0.5), of which 512 (89.2%) were previously unreported. These variants had larger effect sizes (p < 2×10−16, t test) on blood traits as expected and were enriched for protein-coding consequences compared to other variants with similar PPFM and/or lead conditional independence (27.2% versus 4.86%, χ2-test p < 2.2×10−16; Figure 4A, Figure S3E). Remarkably, these rare variants were strongly enriched for assignment to Mendelian blood genes (OR = 3.2, Fisher’s p = 4.22×10−14), even after excluding known pathogenic variants (Table 1; OR = 2.9, Fisher’s p = 4.46×10−11), but were not enriched for non-Mendelian genes (OR = 1.2, Fisher’s p = 0.18). These data support the hypothesis that a small group of high-effect rare variants disproportionately affect core genes for a complex trait.

15Given the large effects of these variants on blood traits, we next sought to test their pleiotropic associations with other clinical and disease traits. Thus, we performed a phenome-wide association study (PheWAS) for 456/574 rare variants using summary statistics for 529 well-represented clinical phenotypes from the UK Biobank cohort (https://www.leelabsg.org/resources). There were 112 significant associations involving 27 variants (Bonferroni-corrected p threshold of 9.45×10−5; Figure 4B, Table S8), of which 110 (98.2%) are currently unreported in the GWAS Catalog. Several biologically coherent associations stand out, including the missense variant rs78534766 in ADCY7, associated with autoimmune conditions (hypothyroidism, inflammatory bowel disease) and eosinophil traits; several variants near PIEZO1 associated with varicose veins and erythroid traits; and a variant (rs45611741) in the 5′ UTR of APOA5 associated with hypercholesterolemia and Mean Corpuscular Volume (MCV), as well as Mean Corpuscular Hemoglobin Concentration (MCHC). Altogether, the PheWAS analysis revealed a variety of novel and relevant disease associations for rare blood trait variants and could point toward common mechanistic roles for these pleiotropic loci.

16Splice-altering genetic variants are a prevalent and under-recognized class of variation underlying genetic disorders and complex trait regulation. We hypothesized that a subset of blood trait variants, especially those that are rare with large effect sizes, may be mediated by splice alterations. We utilized a state-of-the-art neural net classifier, SpliceAI, to predict FM variants with splice-altering consequences. The delta score has been shown to closely track with the validation rate of cryptic splice variants, thus approximating its splice-altering probability. Across 215,694 FM variants (PPFM>0.1%), we identified 109 variants with a putative splicing consequence in 106 unique genes (delta score > 0.2) (Figure S3F). Of these, 11 (10%) were rare (MAF < 1%) and confidently fine-mapped (PPFM>0.5; Table S9). Strikingly, 85% (93/109) of the variants, including 9/16 with delta score > 0.8, fell in non-canonical splice sites, meaning they lie outside the essential GT and AG splice junction dinucleotides. In addition, putative splice variants had lower MAF (Mann-Whitney U p = 5.08×10−8) and higher PPFM (Mann-Whitney U p = 9.89×10−6) compared to other FM variants (Figures S3G and S3H). Even when matched by MAF and PPFM, splice variants also had a 1.5-fold higher GWAS effect size (Mann-Whitney U p = 0.007). To validate these in silico predictions, we examined isoform variation in RNA-sequencing data of 465 participants from the Geuvadis project and used the LeafCutter tool to identify splicing quantitative trait loci (sQTLs). After excluding variants with insufficient statistical power in GEUVADIS, LeafCutter quantified differential splicing effects for 28/109 (26%) putative splice variants. Of these, 23/28 (82%) were identified as sQTLs at a 5% false discovery rate. For example, two common variants falling within 95% credible sets for PLT and lymphocyte count (rs8113779, MAF = 16% and rs12898397, MAF = 37%) were predicted to produce donor gain splice alterations in CD3EAP and ULK3 respectively. These effects were validated by LeafCutter (rs8113779, adjusted p = 3.54×10−47; rs12898397, adjusted p = 3.39×10−104), with alternative splice sites produced by these variants in Geuvadis (Figures 4C and 4D). Finally, we highlight a previously unreported splice variant which was too rare to be quantified by LeafCutter but has interesting biological connections. rs139178017 (MAF = 0.53%) is a strongly FM variant in a novel association locus for red cell distribution width (RDW) (PPFM = 0.73) and MCV (PPFM = 0.4). It is predicted to induce a donor loss splice alteration for transferrin receptor 2 (TFR2), a partner of the erythropoietin receptor and a known regulator of erythropoiesis. Compared to non-carriers, the 4 carriers of rs139178017 harbored substantially increased transcripts with intron retention adjacent to this variant (Figure 4E). These findings support the idea that large GWAS are well powered to identify splice variants with large phenotypic effects, and these splice variants represent a currently under-appreciated mechanism of trait regulation in GWAS loci.

17Our study has identified the largest number of variants ever associated with a single group of correlated phenotypes. While each common (MAF ≥1%) variant accounts for a small effect, their joint effect may be substantial. We used different variant selection criteria to build weighted polygenic scores (PGSs) based on the UK Biobank study and selected the one yielding most predictive power for the 29 blood measurements in an independent cohort (INTERVAL study; Method Details; Table S10). Remarkably, PGS based on hundreds of common sentinel variants (135–689 depending on trait) were shown to be more predictive than larger SNP sets employing more liberal significance thresholds, in line with findings for autoimmune diseases but in contrast to other common human diseases. The proportion of phenotypic variance explained (R2) by the PGS ranged between 2.5% for basophil count to 27.3% for mean platelet volume. Estimates obtained for the same score in an independent cohort of 2,314 French Canadians (CARTaGENE) for 15 available traits were broadly comparable, confirming portability of the PGS between European-ancestry groups (Figure 5A). The causal relationship between genetic variants determining eosinophil count and asthma risk has been previously demonstrated. Focusing on this exemplar disease, we can show that the eosinophil count PGS was also significantly associated with asthma incidence in UK Biobank (odds ratio [OR] = 1.17, 95% confidence interval [CI] = 1.13–1.21, p = 1.02×10−19), suggesting the potential utility of PGSs for blood biomarkers in the clinic.

18Intriguingly, the behavior of the PGSs suggests that the current discovery sample sizes may have achieved saturation of biological signals for blood cell traits. To begin to test this hypothesis, we modeled different discovery measures (total number of variants, loci, genes, and heritability explained) as a function of increasing discovery sample sizes. The best-fitting model shows a quadratic rate of discovery decrease across all tested measures and traits (Figures S4A–S4D). However, while the total numbers of associations detected does not seem to reach a plateau, the heritability explained does (Figures 5B and 5C), suggesting that GWAS with larger sample sizes will provide new discoveries, but of smaller and smaller effects, with the exception of unobserved rare variants. In line with the fact that variants assigned to Mendelian genes have higher effect sizes, these showed a faster saturation curve compared to other genes (Figure S4D). Larger independent discovery datasets will be required to conclusively validate this observation.

19Finally, we wondered if multiple sentinels at a single locus could underlie associations with complex diseases and help define an allelic series at pharmacologically relevant genes. 20% of blood trait loci had ≥2 sentinels, including some unusually large sets (Figure 5D). We overlapped these regions with colocalization results for 18 common human diseases (Figure 5E). Figures 5F–5J show 6 instances of such conditionally independent variant sets, of which 3 involve a known drug target. For example, the type I diabetes (T1D) locus tagged by rs5845323 on chromosome 9 contains one rare and six common variants, all associated with eosinophil percentage (Figure 5I). While the colocalizing T1D variant is intronic in C1QTNF6 gene, the coding-synonymous one from the series is in IL2RB (interleukin 2 receptor subunit beta). It has recently been proposed that the cancer drug Aldesleukin (recombinant IL-2, which binds IL2RB) may be repurposed to treat T1D at low doses, and the drug is currently in phase II clinical trial for this therapeutic application (ClinicalTrials.gov ID: NCT01862120). There were 3 colocalizing loci between asthma and eosinophil count and/or percentage and a further three novel PheWAS associations of rare non-coding variants near known asthma genes (GATA3, RAD50, and IL33). One of the rare variants is part of a 270-kb set of sentinels on chromosome 5 associated with eosinophil count, including another rare variant and 5 common signals (Figure 5K). The genes implicated are C5orf56 (IRF1-AS1 or IRF1 antisense RNA 1), IRF1, IL5, RAD50, IL13, KIF3A, and IL4. Interestingly, both IL5 and IL4 are current therapeutic targets for treating a number of allergic diseases. Overall, this large set of conditionally independent variants informs future efforts to define allelic series to study genes of pharmacological importance.

20Mendelian blood disorders display considerable heterogeneity in penetrance and expressivity. Furthermore, estimates of effect size and penetrance of pathogenic variants tend to be inflated when ascertained from patient populations. While the PGSs defined by the common variants discovered in this study explain a substantial proportion of variance of respective phenotypes, the extent to which polygenic variation contributes to the manifestation of rare diseases remains to be determined. To address this question, we first explored the genetic landscape of classical blood disorders in UK Biobank. We annotated each protein-coding sentinel variant using (1) ClinVar, (2) Human Gene Mutation Database (HGMD), and (3) a recently curated list of variants for rare blood disorders from the Rare Disease Pilot for the 100,000 Genomes Project (NIHR-RD). Overall, 101 sentinels were included in one or more databases above, of which 80% were coding, 10% were annotated to 3′ or 5′ UTRs, and the remaining were splice or intronic variants. 16/101 (16%) were annotated to be pathogenic in either ClinVar or HGMD using strict criteria (Method Details), involving 11 genes (PKLR, HFE, HBB, PIEZO1, TMPRSS6, JAK2, MPO, CSF3R, MPL, GP9, and WAS; Table 1). Only 5/16 variants satisfied the pathogenicity criteria in both ClinVar and HGMD. Of these five, two variants previously reported as pathogenic for autosomal recessive diseases (rs116100695 in PKLR for pyruvate kinase deficiency of red cells and rs1800730 in HFE for hemochromatosis) were found in apparently healthy homozygous UK Biobank participants. Similarly, we found apparently healthy homozygous carriers for other four recessive variants, reported as pathogenic in HGMD, but not in ClinVar (rs61755431 in PKLR for pyruvate kinase deficiency or red cells, rs138156467 in CSF3R for neutropenia, rs61745086 in PIEZO1 for dehydrated stomatocytosis and rs41316003 in JAK2 for erythrocytosis and thrombocytosis). This lack of disease phenotype may be indicative of low penetrance, missing health record data, misannotation of the pathogenicity, or undiscovered compensatory effects either by rare variants or polygenic variation. For two additional recessive variants (rs137853120 in TMPRSS6 for iron-refractory iron deficiency anemia and rs5030764 in GP9 for Bernard-Soulier Syndrome) we observed no homozygous carriers, but heterozygous carriers were around 3 times more likely to have blood indices outside the normal range (hemoglobin < 12 g/dl, PLT < 150×109/l), demonstrating previously unreported dosage-dependent effects (OR = 3.25, 95% CI = 1.85, 5.37, p = 5×10−5 [Figures 6A and 6B] and OR = 3.79, 95% CI = 2.40–5.68, p = 1.1×10−9, respectively). Data were inconclusive for the remaining 8 variants, either because there were no homozygous carriers in UK Biobank (rs113403872 in PKLR, rs28928907 in MPL, rs33946267 in HBB, rs146220228 in WAS), or because the disease presented mild symptoms that are not easily detectable (rs35897051 and rs119468010 in MPO for myeloperoxidase deficiency).

21We next compared the effects of PGS and rare monogenic variants. The average effect of each standard deviation of PGS ranged from 0.16 to 0.47 SD (depending on trait) and was thus comparable to that of a rare pathogenic variant carried in heterozygosity (Figure 6C). We hypothesized that in the previous examples, a low penetrance in a rare disease could be explained by background polygenic variation, for instance in cases where the rare disease mutation carriers have a polygenic effect in the opposite direction that compensated for a high-impact rare mutation. However, the PGS of identified homozygous rare variant carriers was not different from the population mean (defined arbitrarily as PGS > 2SDPGS or tested by logistic regression for variants with more than 10 homozygotes). Hence, the polygenic effects alone were not sufficiently extreme to explain the lower disease prevalence in homozygous rare variant carriers, at least with our current PGS definition.

22Comparing effect sizes can also be used to screen for potential new pathogenic mutations. Variants of uncertain significance (VUS) and missense variants showed a broad distribution of effect sizes, with tails approaching the range of pathogenic ones, and could harbor putative new pathogenic variants. Among the 16 missense variants with the largest effect sizes, two previously uncharacterized ones were in known Mendelian genes (rs139473150 in TUBB1, associated with platelet count and rs201514157 in SPTA1 associated with immature reticulocytes). Platelet count was also associated with rs149254521 in PEAR1, a gene previously identified by an intronic variant in a platelet GWAS. Subsequent functional studies showed that this gene is involved in platelet aggregation, which is consistent with the phenotypes observed here. Two missense variants associated with monocyte count (rs140221307 and rs149771513) were in IL17RA, which has been implicated in monocyte function in mice. Similarly, E2F4 is known to be essential in mouse erythropoiesis, and here its missense variant rs61735430 was strongly associated with mean reticulocyte volume. The other variants were in TIE1 and PLEKHO2 associated to platelets (rs140190628 and rs143331139 respectively), IFRD2 (rs200622087) associated with reticulocytes, TF (rs8177318 and rs150854910) with mean corpuscular hemoglobin (MCH) and MCV, CXCR2 (rs61733609 and rs55799208) with white blood cells, C10orf54 (rs201859625, in VSIR encoding V-set immunoregulatory receptor) with monocytes, and FAM46C (rs148397151) with MCH. The high predictive value of the PGS, alongside future catalogs of additional rare and private variants from whole genome sequencing (WGS), will enable robust modeling and could help explain the heterogeneity in many monogenic blood disorders.

23As shown earlier, polygenic contributions can yield effect sizes on blood traits comparable to or larger than that of known monogenic variants. We therefore sought to explore whether they affect predisposition to rare diseases of the blood. In UK Biobank, we extracted ICD10 codes for a total of 29,080 patients and controls, with 423 blood diseases. We then considered the subset of these participants that were excluded from the discovery GWAS and we estimated their weighted sentinel-based PGS (as detailed earlier). We then fit logistic regression models to test the associations of the PGS with rare disorders of the blood. For the first time, we showed that PGSs derived for blood parameters can influence risk for several rare blood disorders (Figure 6D). For instance, we showed that a higher PGS for red cell count was positively associated with incidence of secondary polycythemia, a disorder characterized by elevated hematocrit (OR = 1.55, 95% CI = 1.21–2.00, p = 6.5×10−4). A high PGS for MCH was protective for iron deficiency anemia (OR = 0.86, 95% CI = 0.79–0.94, p = 9.1×10−4). A high PGS for neutrophil count decreased the risk of aplastic anemia (OR = 0.74, 95% CI = 0.64–0.87, p = 2.9×10−4), which manifests as cytopenias due to depletion of hematopoietic stem cells and failure of blood cell production. Finally, the PLT PGS was negatively associated with thrombocytopenia or low PLT (OR = 0.83, 95% CI = 0.76–0.91, p = 3.8×10−5). We replicated findings for platelet related PGSs in an independent cohort of 1,199 BPD patients and 7,308 controls with WGS (Figure 6D). We showed that an increased PGS for PLT resulted in a protective effect against BPD disorders, including thrombocytopenia (OR = 0.92, 95% CI = 0.86–0.98, p = 0.007). These results refine our understanding of rare blood disease heterogeneity and the contribution of the polygenic background of an individual to the manifestations of a rare disease known to be caused by high-impact pathogenic variants.

24Hematopoiesis is a highly regulated hierarchical process. Genetic variation leading to alteration of blood cell counts can teach us fundamental lessons about this process and serve as a paradigm for studying complex trait genetic architecture. Here we present a large set of GWAS association results for a set of traits that illuminate numerous aspects of hematopoiesis. The magnitude of this discovery set enabled unprecedented statistical power to explore current paradigms in complex trait genetics as well as build a bridge between GWAS in general population cohorts and existing knowledge of monogenic blood disorders.

25While an omnigenic model has recently been proposed to explain complex disease or trait architecture, this framework has been met with some skepticism. Demonstrating the validity of the omnigenic model has relevance for disease gene discovery for differentiating genes with a potential to underpin pathological variation in human traits from those that control variation within healthy physiological ranges. Here we leveraged the knowledge accrued on the hematopoietic system to carry out a first empirical assessment of this model in the context of blood cell trait variation. By defining core genes as those that are found in Mendelian blood disorders and using a coexpression network, we describe properties consistent with this omnigenic model. Specifically, when compared to other GWAS-associated genes, Mendelian genes (1) were enriched among GWAS hits, (2) harbored variants with larger effect sizes, (3) had predominantly blood-specific effects, (4) were coregulated, (5) had central properties in the coexpression network, and (6) were enriched among trans-eQTL targets. Other observations were inconclusive or suggestive of a more continuous (infinitesimal) pattern of inheritance: (1) there was an apparent continuum of effect sizes of variants on disease odds, (2) there was no significant difference in functional scores (chromatin accessibility, deltaSVM) between variants assigned to “core” versus “peripheral” genes, and (3) first- and second-degree neighbors of Mendelian genes were compatible with “core” and “peripheral” functions. The relatively incomplete ascertainment of rare and private DNA sequence variation in both population cohorts and rare disease cases implies that our knowledge of core genes is likely incomplete. Despite these uncertainties, the emerging picture of underlying network connectivity regulating blood traits harbors potential for discovering new pathogenic genes and drug targets. As an example, we identified a subset of 11 closely coexpressed genes (including three known platelet genes (GP9, ITGA2B, and GP1BB) that is coregulated by the same trans-acting eQTL in the ARHGEF3 gene. The use of a large dataset with concurrent genetic and gene expression data in different cell states will be necessary for further quantitative validation of this model.

26Polygenic variation has a substantial contribution to variation in complex quantitative traits and disease risk, sometimes yielding effects comparable to those of rare pathogenic variants. Using only the sentinel signals from our discovery GWAS, we built PGSs explaining up to 28% of phenotypic variance. We explored the polygenic effects jointly with pathogenic variants and as phenotype modulators in patients with rare blood disorders. While we found that 16 known monogenic variants were each associated with quantitative blood traits, 52 participants homozygous for five rare recessive pathogenic variants appeared to be healthy with normal blood count and indices. This suggests that the penetrance of pathogenic variants may be overestimated in many instances, as was recently shown. Differences in PGS could not explain the reduced penetrance, but our analysis may be limited by the diseases we had adequate statistical power to assess. Conversely, we observed strong allele dosage-dependent effect sizes for two heterozygous variants (previously reported as recessive), that could lead to disease especially if coinherited with an adverse PGS. For example, heterozygous carriers of the GP9 variant rs5030764 were three times more likely to have a PLT below the normal range (< 150K/ul). Finally, we observed a significant association between phenotype-relevant PGSs and rare blood disorders for thrombocytopenia, secondary polycythemia, anemia, and aplastic anemia, regardless of the presence or absence of known rare variants in patients. This highlights a substantial polygenic modulating effect on presumably monogenic disorders and lays the groundwork for future studies aiming to define the impact of genetic background on the variable penetrance and expressivity in blood disorders.

27In summary, through the largest study of blood cell trait variation to date, we provide new insights into the regulation of blood cell parameters and how genetic variation may contribute to the variability observed in rare blood disorders that are presumed to have a monogenic etiology. Our findings provide a novel framework for considering an individual’s genetic background and how this may impact the presentation of blood diseases. Finally, the lessons learned from this study of hematopoiesis will likely be more broadly applicable to a wide range of other complex diseases and traits.

28Further information and requests may be directed to the Lead Contact, Nicole Soranzo (ns6@sanger.ac.uk).

29This study did not generate new unique reagents.

30Summary statistics are available to download from: ftp://ftp.sanger.ac.uk/pub/project/humgen/summary_statistics/UKBB_blood_cell_traits/ for UK Biobank and http://www.mhi-humangenetics.org/en/resources for the meta-analysis. The accession numbers for the UK Biobank summary statistics reported in this paper are GWAS Catalog: GCST90002379–GCST90002407.

31The code generated during this study is publicly available at GitHub https://github.com/bloodcellgwas/manuscript_code/.

32Following the success of the Blood Cell Consortium ( and), the Blood Cell Consortium Phase 2 (BCX2) continues to identify novel common and rare variants associated with blood cell traits using imputed genotype data based on Haplotype Reference Consortium (HRC) or the 1000 Genomes Project (Phase 3, version 5) for European ancestry cohorts and non-European ancestry cohorts, respectively. BCX2 comprises 746,667 participants from 40 discovery cohorts and five ancestries: European, African American, Hispanic, East Asian, and South Asian. BCX2 is divided into two working groups: European that consists of 563,946 participants from 26 cohorts and that is the focus of this study, and trans-ethnic. Detailed descriptions of the participating cohorts are summarized in Table S1. All participants provided written informed consent, and local research ethics committees and institutional review boards approved the individual studies.

33Genotyping array and pre-imputation quality control (QC) for each participating cohort is provided in Table S2. Genotype QC metrics included MAF (> 0), call rate (> 98%) and Hardy-Weinberg equilibrium p (> 10−6). The pre-imputation sample exclusion criteria (Table S2) included call rate (> 95%), heterozygosity rate (> median+3∗IQR), gender mismatches, duplicates, and outliers from principal component analysis with reference samples from 1000 Genomes Project. All genotypes were on Genome Reference Consortium Human Build 37 (GRCh37) forward strand (https://www.well.ox.ac.uk/∼wrayner/strand/). All the cohorts checked strand and allele orientation in the variant call format files prior to imputation using checkVCF (https://genome.sph.umich.edu/wiki/CheckVCF.py). Finally the imputation was performed using servers available at https://imputation.sanger.ac.uk/ or http://imputationserver.sph.umich.edu/index.html with requesting the HRCr1.1 2016 reference panel, EUR population and Quality Control & Imputation Mode. All cohorts followed this procedure for the imputation of autosomal variants except for UK Biobank (UKBB) and INTERVAL that had their genotype imputation described elsewhere.

34When possible, we excluded samples with any of the following: pregnancy (when complete blood count (CBC) done), acute medical/surgical illness (when CBC done), blood cancer, leukemia, lymphoma, chemotherapy, myelodysplastic syndrome, bone marrow transplant, congenital or hereditary anemia (e.g., hemoglobinopathy such as sickle cell anemia or thalassemia), HIV, end-stage kidney disease, dialysis, EPO treatment, splenectomy, cirrhosis and those with any of the following extreme measurements: WBC count > 100∗109/L with > 5% immature cell or blasts, WBC > 200∗109/L, Hemoglobin > 20 g/dL, Hematocrit > 60%, Platelet > 1000∗109/L. For the WBC subtypes (e.g., basophils count) we used the relative count, i.e., the total WBC count multiplied by the proportion for each cell type (e.g., basophils percentage). Raw phenotypes were regressed on age, age-squared, sex, principal components and cohort specific covariates (e.g., study center, cohort, etc) if needed, WBC related traits were log10 transformed before regression modeling. Residuals from the modeling were obtained and then inverse normalized for cohort level association analysis or GWAS. All cohorts followed the same exclusions and phenotype modeling except for UKBB and INTERVAL that had their procedure described elsewhere. The cohort level association analyses were then conducted using a linear mixed effects model in order to account for known or cryptic relatedness (e.g., BOLT-LMM, EPACTS https://github.com/statgen/EPACTS and rvtests with the additive genetic model. Linear mixed effects models have been shown to effectively account for both population structure and inter-individual relatedness within the UK Biobank cohort, along with having increased discovery power over simple linear regression with principal components.

35Cohort level association analysis results went through a standard QC procedure using EasyQC R package (https://www.uni-regensburg.de/medizin/epidemiologie-praeventivmedizin/genetische-epidemiologie/software/). The mapping file and allele frequency reference data (GRCh37/hg19) from HRC were used to harmonize variant names across cohorts and to check allele frequency discrepancies between cohorts and the HRC reference panel, respectively. We generated a unique ID for each variant using the form of chromosome:position_allele1_allele2 where alleles were ordered lexicographically or based on indel length as tri-allelic variants and/or indels of the same chromosome:position were observed. In addition to allele frequency plots, quantile-quantile (Q-Q) plots and SE-N (i.e., inverse of the median standard error versus square root of sample size) plots were also checked to detect systematic inflation, and different phenotypic variances due to mis-specified phenotype transformation or regression model, different study design or different study population, etc.

36Post-QC’ed and pre-processed European cohort results were then meta-analyzed by GWAMA using inverse variance weighted fixed effects approach. We applied an imputation quality filter of INFO score ≤0.4 and a minor allele count (MAC) filter of MAC ≤5 for each variant in the meta-analysis, except for three large cohorts UKBB (N = 487,409), WHI (N = 17,682) and GERA (N = 53,822), where a more stringent MAC filter of MAC ≤ 20 was applied to exclude extremely rare variants with extreme effects prior to meta-analysis.

37We performed an exact conditional analysis using a stepwise multiple linear regression approach in UKBB. Stepwise multiple linear regression aims to identify a parsimonious subset of variants which explain the significant associations identified by univariable GWAS. For each blood phenotype, the set of genome wide significant variants was partitioned into the largest number of blocks such that no pair of blocks are separated by fewer than 5Mb, subject to the restriction that no block contained more than 2,500 variants. Blocks were generated independently for each phenotype. For each phenotype and each block, we identified a parsimonious set of variants explaining the signal in that block using a stepwise conditional linear regression algorithm. Each iteration of the algorithm had two stages: i) addition of variants to the model and ii) removal of variants from the model. Convergence occurred when neither addition nor removal of any variant improved the model fit sufficiently for a t test p < 8.31x10−9. The variants in the model at convergence represent a parsimonious set for the block.

38Let “M” represent the ‘current model,’ The algorithm is initialised with M as the empty model, containing just an intercept term, and develops with the following steps:

391. Update M by inserting the variant in the block with the lowest univariable association p value, into the model.

402. In turn for each variant∗ in the block not in M, compare M to the model generated by augmenting M with the variant using a t test. Record the p value from each comparison.

413. If the least p value recorded in step 2 is greater than 8.31x10−9 terminate the algorithm.

424. Update M by adding the variant with the least p value recorded in step 2.

435. In turn for each variant in M, compare M to the model generated by removing the variant from M, using a t test. Record the p values from each comparison.

446. If the greatest p value recorded in step 5 is smaller than 8.31x10−9 go to step 2.

457. Update M by removing the variant with the greatest p value recorded in step 5.

468. Go to step 5.

47∗ When comparing M to M augmented by a variant, we are testing to see if this new variant represents a genetic signal independent of the variants in M. In the situation where a potential new variant is in high LD (r2 > 0.9) with a variant already in M we assume that the this variant cannot represent an independent signal and we do not proceed to calculate its P value.

48All linear regression was performed using the fastLM from the R package RcppEigen.

49For each phenotype, following identification of conditionally significant variants in each block, all conditionally significant variants within each chromosome were put into a single linear model and tested with the same multiple stepwise linear regression algorithm described above, but starting at step 5. The union across chromosomes of the resulting sets of variants is the ‘conditionally significant’ list of variants for the blood cell phenotype, also referred to as “sentinel variants” throughout the text.

50Using conditional and joint analysis as implemented in GCTA (https://gump.qimr.edu.au/gcta), we identified independent association results in the meta-analyses at p < 5x10−9. To define novel associations, we tested these variants using the same exact multivariate approach as above, in the UKBB, while conditioning on the variants identified by the previous step.

51We checked replication in the Million Veteran Project Cohort, for chromosome 1 and three different traits, one per each major cell type (platelet counts, lymphocyte counts and red cell counts). The replication significance threshold was set to a nominal level (p < 0.05) with the same direction of effect.

52Statistical fine-mapping was performed in the UKBB cohort, using FINEMAP v1.3.1 (http://www.christianbenner.com/). Input windows were defined as +- 250 kb from a conditionally independent signal. In case of multiple sentinels generating overlapping windows these were merged together, resulting in window size ranging from 500kb to 1.38Mb. The number of conditionally independent signals in each window was used as prior knowledge for the maximum number of causative variants to be searched (–n-causal-snps option) and the prior standard deviation for effect sizes was set to 0.08 (–prior-std option). The LD structure was computed from the same samples included in the GWAS analysis. 95% credible sets were defined as minimal sets of variants jointly covering at least 95% of the posterior probability of including the true causative signals.

53We performed colocalization using gwas-pw between GWAS of 10 hematological traits and transcriptomic profiling of Platelets, CD4+, CD8+, CD14+, CD15+, and CD19+ cells. Where MPV, PDW, PLT#, and PCT were colocalized with eQTLs from Platelets, NEUT# and NEUT% were colocalized with eQTLs from CD15+ cells, LYMPH# and LYMPH% were colocalized with eQTLs from CD4+, CD8+, and CD19+ cells, and MONO# and MONO% were colocalized with eQTLs from CD14+ cells. Loci of colocalization were defined by the recombination regions identified by. Colocalization was only performed if the locus contained a conditionally independent variant in LD r2 > 0.8 with the eQTL sentinel. Results were filtered to include only those with posterior probability for colocalization higher than 80% resulting in a set of colocalized loci which are considered ‘highly likely’ to be colocalized.

54The list of Mendelian genes was retrieved from a manually curated list, compiled as part of the NIHR BioResource rare disease sequencing project. It includes genes causative for stem cell and myeloid disorders (SMDs, 206 genes), bleeding, thrombotic and platelet disorders (BPD, 104 genes) and bone-marrow failure syndromes (BMF 80 genes). We refer to Mendelian SNPs as those assigned by VEP (worst consequence) to one of the Mendelian genes. To test for differences in absolute effect sizes we matched for MAF between Mendelian and other SNPs using the R package “MatchIt”, separately for each functional annotation, with at least 10 variants per group. We then tested the absolute effect size distribution shift using the Wilcoxon test as implemented in R.

55A co-expression matrix computed from the whole blood of 2,168 participants was used. The matrix quantifies correlations between genes, replicated across 2 different cohorts. The edges between genes were defined by imposing variable hard cut-offs on co-expression coefficients, e.g., two genes are linked in the network if their co-expression is higher than the cut-off. The following cut-offs were used (0.05, 0.1, 0.2, …, 0.8) but the results did not generally depend on the specific cut-off (unless otherwise stated). The overlap enrichment between GWAS genes and network genes was computed by random permutations of gene sets, in particular we used the following steps:

561. Annotate GWAS associations by VEP - this gives us a gene annotation for 83% of GWAS variants, which we refer to as “GWAS genes”

572. Iterations:

58- Randomly select a set of genes from the Ensembl v99 list of protein coding genes (using the R function “sample”). The size of the set matches the size of the test set

59- Overlap the random set of genes with GWAS genes to calculate how many of the random genes are among the GWAS associations

60- Repeat 10,000 times

613. Compare the observed overlap with the background distribution

62The numbers of links per gene were compared between Mendelian genes and all other genes by Wilcoxon test. For co-expression among Mendelian genes, median absolute co-expression coefficients were computed for equal sized random draws of genes.

63The enrichment of trans-eQTLs targeting Mendelian genes was computed similarly to the above comparison between effect sizes of Mendelian and other GWAS genes. First, we downloaded median gene expression levels in whole blood from one of the cohorts included in eQTLGen, BIOS (http://bbmri.researchlumc.nl/atlas/#data). We rank-inverse normalized the median expression levels. Then, to account for differences in power detection due to higher expression levels of Mendelian genes and higher Z-scores of Mendelian-targeting trans-eQTLs, we used the R package “MatchIt” to select a matched subset. This included trans-eQTLs targeting non-Mendelian genes with matched Mendelian expression levels and Z-scores for trans-eQTLs targeting Mendelian genes. We included the maximum possible number of trans-eQTL-gene pairs with these characteristics (N = 9,258). Finally, we compared the number of trans-eQTLs per gene in the two groups by Wilcoxon test.

64The enrichment of GWAS genes among first and second degree neighbors to Mendelian genes was computed as follows: i) determine the list of neighboring genes based on the specific cut-off, ii) intersect with 1000 random permutations of gene sets of the same size as the GWAS list, iii) compare to the actual intersection. The second degree neighbors were defined as neighboring genes to all first-degree neighbors.

65Bias-corrected enrichment of blood trait variants for chromatin accessibility of 18 hematopoietic populations was performed using g-chromVAR, whose methodology has been previously described in detail. In brief, this method weights chromatin features by fine-mapped variant posterior probabilities and computes the enrichment for each cell type versus an empirical background matched for GC content and feature intensity. For chromatin feature input, we used a consensus peak set for all hematopoietic cell types with a uniform width of 500 bp centered at the summit. For variant input, we included all variants with fine-mapped PPFM > 0.1%.

66DeltaSVM is a machine learning model which uses sequence composition to predict cell type-specific open chromatin. It then uses this sequence vocabulary to predict the change in chromatin accessibility from each variant. We trained on two ATAC-Seq datasets: 1) 18 hematopoietic populations sorted from bone marrow, and 2) 8 stages of primary erythroid differentiation. For each dataset, we trained on strong ATAC peaks in the > 80th percentile of counts matrix from each cell type. Standard 5-fold cross-validation was used to calculate AUROC. We then scored each variant with a posterior probability of association greater than 0.001 for all populations to determine variants predicted to alter chromatin accessibility. Here, a positive deltaSVM score is interpreted as a prediction where the variant increases chromatin accessibility whereas a negative score would reduce chromatin accessibility.

67Prediction of the effects of fine-mapped variants on transcription factor binding sites (TFBS) was performed by using the motifbreakR package and a collection of 426 human TFBS models (HOCOMOCO). For 115,609 fine-mapped variants with PPFM > 0.1%, we applied the ‘information content’ scoring algorithm and used a p cutoff of 1 × 10−4 for TFBS matches; all other parameters were kept at default settings.

68To identify associations between blood trait variants and clinical phenotypes, we conducted a phenome-wide association study (PheWAS) using summary statistics of 1,403 clinical phenotypes analyzed from the UK Biobank (https://www.leelabsg.org/resources). As input, we started with 574 rare variants with 0.00005 < MAF < 0.01 which were either conditionally independent lead signals or had fine-mapped PPFM > 0.50. To avoid studying phenotypes with too few cases to capture these low allele frequencies, we only included phenotype-variant results for which the expected_case_minor_AC, calculated as 2 ∗ variant_MAF ∗ num_cases, was greater than 25. This resulted in the final inclusion of 529 clinical phenotypes (case numbers ranging from 1,236 - 77,977) across 456 variants which had pheWAS data. The Bonferroni-corrected significance threshold for pheWAS was calculated as 0.05 / 529 phenotypes = 9.45x10−5.

69To predict splice variants, we used SpliceAI, a deep neural network that accurately predicts splice junctions from genomic sequence. We obtained prediction scores for all possible single nucleotide variants in the reference genome, which were released along with the SpliceAI tool, and extracted scores for all variants with fine-mapped PPFM > 0.001 in one or more blood traits from the UK Biobank GWAS. We considered a variant to have a putative splicing consequence if it had a delta score > 0.2 for one or more splicing consequences (acceptor gain, acceptor loss, donor gain, donor loss); this threshold was shown to be enriched for splice variants and have high sensitivity.

70Validation of SpliceAI predictions was performed using RNA-seq data on lymphoblastoid cell lines (LCL) of 465 participants from the Geuvadis project 23 (http://www.geuvadis.org). We aligned paired-end reads to the hg19 reference genome with STAR, allowing for novel splice junctions. To systematically evaluate predicted splice-altering variants, we processed junction files for all 465 samples using the LeafCutter workflow and evaluated changes in splicing clusters that overlapped the variant. For each variant, read alignments were merged into two groups, all carriers versus all non-carriers, and visualized in the form of sashimi plots using the ggsashimi tool. The threshold for the minimum number of reads supporting a junction to be drawn was 100 for rs8113779 and rs12898397, and reduced to 15 for rs139178017 given the low number of carriers.

71The polygenic scores (PGSs) were computed as weighted sums of genotypes, weighted by their effect size on the phenotype (beta coefficient), using the PLINK score function. Beta coefficient estimates were computed in UK Biobank and PGS scores were tested in an independent cohort (INTERVAL). Positive effect alleles were included in order to get a positive contribution for each carried allele and consequently a positive correlation with the phenotype. The following SNP inclusion criteria were compared: a) all genome-wide SNPs after LD pruning with PLINK at 0.8 cut-off; b) LD-pruned SNPs with GWAS p < (0.05, 5x10−4, 5x10−6, 5x10−8); c) conditionally independent variants and fine-mapped variants with posterior probability > 0.5; d) conditionally independent variants. Resulting PGSs were then standardized and Pearson’s R coefficients of correlation between each PGS and its relevant trait were computed for comparison. The phenotypic variance explained was computed as R2. Validation in a further independent cohort of French-Canadians (European ancestry) called CARTaGENE was performed using the same protocol, for the best performing PGS (conditionally independent variants). A linear regression model between the adjusted phenotypes and the PGS, adjusted by sex, age and principal components, was used to determine the PGS’s effect sizes per SD.

72To explore discovery saturation we chose 4 large GWAS analyses in cohorts of increasing sample size: INTERVAL (N∼35k), UK BiLEVE (N∼43,5k), UK Biobank 1st release (N∼83k) and UKBiobank full cohort (N∼400k). For all of these we had conditionally independent associations identified by the same method as described above. For each trait and cohort we determined the number of conditionally independent variants detected by GWAS; the number of genes identified by these variants (using VEP worst consequence annotation); the number of associated loci and we further subset the genes as Mendelian or others. Associated loci were defined based on LD-blocks computed in Pickrell et al. which had at least one conditionally independent signal. First, a linear projection of the first 3 data points was visually inspected to determine in the 4th data point fitted the expected. Then 4 different regression models were tested to determine which one best described the full dataset: i) y∼x; ii) y∼sqrt(x); iii) y∼sqrt(x)+x; (iv) y∼log(x). Here y represents the counts (number of variants/genes/loci associated) and x represents the cohort size. Similarly we computed the heritability explained by the set of variants identified by each cohort and searched for the best fitting model. The heritability was computed as R2 of the multivariate model including the relevant variants in the full UK Biobank cohort. The model fitting the Mendelian genes versus others was computed across all pooled phenotypes.

73We performed pairwise colocalization analysis between GWAS studies of 29 hematological parameters from the UK Biobank cohort and 18 different autoimmune and inflammatory related disorders. Our analysis was performed using summary statistics collected following GWAS of the respective studies. An inner merge was performed with variants tested for each hematological parameter and each respective disease risk GWAS. Colocalization analysis was then performed following the same protocol described above for eQTLs. Allelic series were defined as fine-mapping blocks including 2 or more associated sentinels.

74We used UK Biobank participants with ICD10 codes in Chapter III (“Diseases of the blood and blood-forming organs,” D500-D899 codes), who were excluded from the GWAS discovery. After computing the PGS, we performed a logistic regression for the disease status, including sex, age, 10 principal components and any other co-occurring blood disorders as covariates. P values were corrected by Bonferroni correction for the number of diseases (i.e., ICD10 codes) tested. We included only ICD10 codes with at least 40 cases for a total of 49 disorders.

75Conditionally independent variants for each trait, as well as fine-mapped variants with posterior probability of being causative greater than 50% were pulled together for the pathogenicity annotation. Genes were assigned to variants by VEP worst consequence [release 84]. The set of Mendelian genes was manually curated by the NIHR-RD project. We focused on three different sources of pathogenicity annotations: the ClinVar database (https://www.ncbi.nlm.nih.gov/clinvar/ - accessed 2018-11-05), the Human Gene Mutation Database, version pro 2018.4 (HGMD) (http://www.hgmd.cf.ac.uk/ac/index.php - accessed 2019-01-17) and a manually curated list of novel pathogenic variants, produced by the NIHR-RD sequencing project. Variants were matched by chromosome, position and alleles, in GRCh37. The following parameters were considered: a) ClinVar: categorical pathogenicity assignment (yes/no/unknown), the star rating (1-4, 1 being the most uncertain and 4 the most certain) and the Review Status; b) HGMD: categorial pathogenicity assignment and Rank Score indicating the pathogenicity confidence on a continuous scale from 0 to 1, 1 being certain pathogenicity assignment and 0 being very uncertain. The set of pathogenic variants was defined with high confidence, imposing pathogenicity in ClinVar with at least 2 stars or pathogenicity in HGMD (“DM”) with rank score greater than 0.1. Variants reported by NIHR-RD for the first time were assigned to the “variants of uncertain significance” category (VUS). To assess the effects of such pathogenic variants and their penetrance, two types of data were considered: full blood count diagnostic cut-offs as used in the clinics and ICD-10 codes for blood disorders (Chapter III), as recorded by UKBB. Participants with full blood counts in the normal range and no ICD-10 code were considered healthy. The joint modeling of rare variants and PGS was performed only for variants with more than 10 homozygotes, using logistic regression and relevant covariates, as above. For variants with less than 10 homozygotes, we checked if these for systematic PGS deviation from the population mean (defined as PGS > 2∗SDPGS).