ARTICLE https://doi.org/10.1038/s41467-021-22548-8 OPEN Genetic perturbation of PU.1 binding and chromatin looping at neutrophil enhancers associates with autoimmune disease Stephen Watt1, Louella Vasquez1,17, Klaudia Walter 1,17, Alice L. Mann1,17, Kousik Kundu 1,2, Lu Chen1,2,13, Ying Sims1, Simone Ecker 3, Frances Burden4,5, Samantha Farrow4,5, Ben Farr1, Valentina Iotchkova 1,6,14, Heather Elding1, Daniel Mead 1, Manuel Tardaguila1, Hannes Ponstingl 1, David Richardson 6, Avik Datta6, 1234567890():,; Paul Flicek 6, Laura Clarke 6, Kate Downes4,5, Tomi Pastinen7, Peter Fraser 8,15, Mattia Frontini4,5,9,16,18, Biola-Maria Javierre 8,10,18 ✉, Mikhail Spivakov 8,11,12,18 ✉ & Nicole Soranzo 1,2,18 ✉ Neutrophils play fundamental roles in innate immune response, shape adaptive immunity, and are a potentially causal cell type underpinning genetic associations with immune system traits and diseases. Here, we profile the binding of myeloid master regulator PU.1 in primary neutrophils across nearly a hundred volunteers. We show that variants associated with differential PU.1 binding underlie genetically-driven differences in cell count and susceptibility to autoimmune and inflammatory diseases. We integrate these results with other multi- individual genomic readouts, revealing coordinated effects of PU.1 binding variants on the local chromatin state, enhancer-promoter contacts and downstream gene expression, and providing a functional interpretation for 27 genes underlying immune traits. Collectively, these results demonstrate the functional role of PU.1 and its target enhancers in neutrophil transcriptional control and immune disease susceptibility. 1 Human Genetics, Wellcome Sanger Institute, Genome Campus, Hinxton, UK. 2 School of Clinical Medicine, University of Cambridge, Cambridge, UK. 3 UCL Cancer Institute, London, UK. 4 Department of Haematology, University of Cambridge, Cambridge, UK. 5 National Health Service Blood and Transplant (NHSBT), Cambridge, UK. 6 European Molecular Biology Laboratory, European Bioinformatics Institute, Hinxton, Cambridge, UK. 7 Center for Pediatric Genomic Medicine, Children’s Mercy, Kansas City, MO, USA. 8 Nuclear Dynamics Programme, Babraham Institute, Cambridge, UK. 9 British Heart Foundation Centre of Excellence, Division of Cardiovascular Medicine, Addenbrooke’s Hospital, Cambridge, UK. 10 Josep Carreras Leukaemia Research Institute, Badalona, Barcelona, Spain. 11 Functional Gene Control Group, MRC London Institute of Medical Sciences (LMS), London, UK. 12 Institute of Clinical Sciences, Imperial College Faculty of Medicine, London, UK. 13Present address: Key Laboratory of Birth Defects and Related Diseases of Women and Children, Department of Laboratory Medicine, West China Second University Hospital, State Key Laboratory of Biotherapy, Sichuan University, Chengdu, China. 14 Present address: MRC Weatherall Institute of Molecular Medicine, Radcliffe Department of Medicine, University of Oxford, John Radcliffe Hospital, Headington, Oxford, UK. 15Present address: Department of Biological Science, Florida State University, Tallahassee, FL, USA. 16Present address: Institute of Biomedical & Clinical Science, College of Medicine and Health, University of Exeter Medical School, RILD Building, Exeter, UK. 17These authors contributed equally: Louella Vasquez, Klaudia Walter, Alice L. Mann. 18These authors jointly supervised this work: Mattia Frontini, Biola-Maria Javierre, Mikhail Spivakov, Nicole Soranzo. ✉email: bmjavierre@carrerasresearch.org; mikhail.spivakov@lms.mrc.ac.uk; ns6@sanger.ac.uk NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 N eutrophils are the most abundant circulating leucocytes in Here, to study the effects of PU.1 on immune function in human blood, comprising up to 80% of the white blood human primary neutrophils, we use natural genetic variation cell population. They form the body’s first line of defence across humans as a high-throughput “perturbation experiment” against infection from pathogens1. As inflammatory cells, they to identify PU.1 binding events in neutrophils that are subject to react to chemotactic signals and migrate to a wound or site of genetic control. We show that variants associated with PU.1 infection to protect the organism from pathogenic insults. Neu- differential binding (PU.1 tfQTLs) are enriched for genetic trophils are phagocytes that envelope and consume microbial associations with cell count and susceptibility to multiple auto- invaders and on top of this capacity, they have an arsenal of immune and inflammatory diseases. Integrating this information antimicrobial properties at their disposal. They are highly cyto- with other multi-individual genomic readouts, we show that the toxic cells full of degradative enzymes including proteases, effects of genetically determined variation in PU.1 binding are hydrolases and nucleases2. In addition, neutrophils have the likely propagated through the local chromatin state and 3D ability to produce and release reactive oxygen species, which can enhancer-promoter contacts to control the expression of disease- cross bacterial membranes damaging nucleic acids and proteins3. relevant genes. Jointly, these results provide evidence for a One of the most striking capabilities of these cells is to use their functional role of PU.1 in neutrophils in immune disease aetiol- own DNA as a form of antimicrobial net in a process known as ogy and provide mechanistic insights into the effects of non- NETosis4. This is a programmed form of cell death where the coding genetic variation associated with healthy and pathological neutrophil can release its chromatin to ensnare and kill larger immune traits. microbes that are difficult to phagocytose such as fungal patho- gens. As a collateral effect of their immune function, neutrophils Results secrete highly cytotoxic molecules that can irritate or damage A multi-individual atlas of PU.1 binding and associated surrounding tissues, contributing to inflammatory and auto- genomic readouts. We generated PU.1 chromatin immunopre- immune diseases5. However, it remains unknown to what extent cipitation sequencing (ChIP-seq) data from CD66b+CD16+ aberrations in neutrophil function can drive immune disorders6,7. neutrophils isolated from fully genotyped individuals forming Pinpointing the causal cell types in these disorders is further part of the BLUEPRINT project cohort14,25, with data from a complicated by the known crosstalk between neutrophils and total of 93 donors passing quality control. In addition, we pro- other immune cell types such as B and T lymphocytes through duced multi-individual data for two histone modifications, the cytokine messaging and direct interactions8. active mark H3K4me3 and the Polycomb-associated inhibitory Genome-wide association studies (GWAS) are powerful tools mark H3K27me3 in neutrophils isolated from the same donor for identifying the genetic determinants of complex traits and cohort (n = 107 and n = 106, respectively). We also profiled the diseases, including autoimmunity. However, GWAS are cell type- binding of the key PU.1 partner TF C/EBPβ and of the archi- agnostic and typically implicate non-coding variants that require tectural protein CTCF in neutrophils (n = 22 and n = 30, further interpretation to be actionable. Causal GWAS variants respectively), as well as the binding of the same TFs in classical preferentially localise to DNA regulatory elements such as monocytes (PU.1, n = 10; C/EBPβ, n = 9; CTCF, n = 4) (Fig. 1a, enhancers, which bind transcription factors (TFs) to activate and Supplementary Fig. 1a, b and Supplementary Data 1). We iden- fine-tune target gene expression9,10. Enhancers are often located tified bound regions (peaks) in ChIP-seq data using established large distances (up to megabases) away from their target genes algorithms (see “Methods”) and confirmed data quality by and typically exert their effects via direct DNA looping contacts assessing the fraction of reads in peaks (FRiP), as well as using with their target gene promoters. Therefore, integrating GWAS principal component analysis and other approaches (Supple- results with the readouts of genome conformation and epigenome mentary Fig. 1c–e and Supplementary Fig. 2a, b). Finally, we function is instrumental for identifying the causal cell types and obtained high-resolution chromosomal interaction profiles of molecular mechanisms behind GWAS signals. For example, gene promoters in six individuals in neutrophils and monocytes enhancer activity is associated with specific patterns of chromatin by combining newly generated and published16 Promoter Cap- accessibility and histone modifications (such as histone H3 lysine ture Hi-C (PCHi-C) datasets (Fig. 1a, Supplementary Fig. 3a–c). 4 methylation and lysine 27 acetylation), and we and others have previously shown that GWAS variants can associate with allele- specific imbalances in the magnitude of these readouts that Genetic determinants of PU.1 binding in neutrophils. To study potentially mediate the variants’ expression effects11–14. The genetically determined variation in PU.1 binding to DNA, we allelic effects of non-coding GWAS variants on gene expression15, identified 36,530 high-confidence PU.1 peaks across the 93 as well as their three-dimensional (3D) contacts with gene pro- individuals (Methods) and used normalised read counts at peak moters determined using high-resolution chromosomal archi- regions to determine TF quantitative trait loci associated with tecture profiling, enable the identification of target genes at PU.1 binding (PU.1 tfQTLs; “Methods”). We detected 1868 GWAS traits16,17. independent PU.1 binding QTLs (linkage disequilibrium [LD] PU.1 (encoded by SPI1 gene) is a master TF controlling r2 ≥ 0.8; global false discovery rate [gFDR] <0.05; Supplementary myeloid development18. PU.1 deficiency has profound effects on Fig. 2c–e and Supplementary Data 2). neutrophil maturation and function19. PU.1 is known to bind Lead PU.1 tfQTL SNPs showed a bimodal distribution of multitudes of enhancer elements in myeloid cells20. As a pioneer distances to their respective differential binding peaks (Fig. 1b, TF, PU.1 is able to bind these regions in a repressive chromatin Supplementary Data 3), with just over half of them (55%, 1,036/ state and activate them via the recruitment of other factors21, 1,868) mapping proximally from the peak edge (<2.5 kb; median including core chromatin activators and TFs such as C/EBPβ, distance 264 bp), and the remaining SNPs (45%; 995/1,868) which is expressed throughout neutrophil development22. PU.1 is localising more distally (2.5 kb-1 Mb, median distance 23 kb). required for modulating the response of mouse neutrophils to tfQTL effect sizes were stronger for proximal compared to distal infection23, but the emerging evidence for a role of this factor in variants (Welch two-sided t-test p = 2.2 × 10−16, Fig. 1c) as complex diseases is currently limited to other cell types. For previously observed in cultured lymphoblastoid cells11,12. We example, a common genetic variant within the SPI1 gene was further validated a subset of the detected tfQTLs using allele- found to be associated with Alzheimer’s disease, with microglial specific association analysis26 (“Methods”), which confirmed a cells proposed as the likely causal cell type24. significant allelic imbalance for the majority of the tested peaks 2 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 ARTICLE Fig. 1 Properties PU.1 transcription factor QTLs. a Summary of molecular traits generated as part of this study. b Density plot displaying the distance between sentinel SNPs and their associated PU.1 peaks. The bimodal distribution (grey) can be further subdivided into proximal (solid navy, <2.5 kb) and distal SNP effects (dotted, >2.5 kb). c Boxplot of absolute PU.1 tfQTL effect sizes (beta). Proximal PU.1 tfQTLs exhibit larger effect sizes compared to distal tfQTLs (Welch two-sided t-test). Box plots show the medians (centre lines) and the 25th and 75th percentiles (box edges), with whiskers extending to 1.5 times the interquartile range. n = the number of PU.1 tfQTL in each category. d Proportion of significant tfQTL SNPs with significant allele-specific (A/S) binding. Peaks without suitable heterozygous SNPs were not tested (grey). (98.8% and 95.5% for peaks associated with proximal and distal specific to neutrophils. To this end, we additionally generated variants respectively; Fig. 1d). PU.1-binding maps in primary (CD14+CD16−) monocytes iso- PU.1 tfQTLs associate with differential binding of other TFs in lated from ten BLUEPRINT donors, five of whom were also part a distance-dependent manner. PU.1 binding to DNA leads to of the cohort used to profile PU.1 binding in neutrophils. The extensive chromatin remodelling27 and redistribution of partner absolute majority (93%) of the PU.1-binding peaks associated co-factors28. C/EBPβ is upregulated throughout neutrophil with a tfQTL in neutrophils were also detected in monocytes terminal differentiation22 and has been shown to co-occupy (Supplementary Fig. 4a). The low number of individuals pre- myeloid enhancers at thousands of PU.1 bound sites20. In cluded an independent PU.1 tfQTL analysis in monocytes. contrast, the constitutively expressed CTCF is known to play a Therefore, to assess coordination of genetic effects at PU.1- role in gene regulation by anchoring chromatin interactions29, binding sites across cell types, we assessed the binding signal at but is not known to functionally associate with PU.1. shared PU.1 and C/EBPβ-binding sites in monocytes with indi- We identified 18,862 C/EBPβ and 22,197 CTCF filtered peaks viduals stratified by lead PU.1 tfQTL genotype. We found that from the datasets combined across individuals. We performed monocytes displayed genotype associated changes in the strength tfQTL analysis as before, identifying 427 C/EBPβ and 769 CTCF of binding at proximal SNPs (linear regression p = 3 × 10−9) as putative tfQTLs that reached a nominal p value threshold (p ≤ 1 × observed in neutrophils (p = 2 × 10−13), compatible with shared 10−5; Supplementary Data 2). To assess the shared genetic effects genetic effects between the two cell types. However, the same was on the binding of PU.1 and these two TFs in cis, we focused on not true for distal SNPs (neutrophils p = 4 × 10−7, monocytes proximal PU.1 tfQTLs and assessed the effect size on the nearest p = 0.793; Fig. 2d). This data suggest that proximal tfQTL var- C/EBPβ and CTCF QTLs (p ≤ 1) with respect to the lead SNP. iants in related cell types operate under a similar cis-regulatory PU.1 tfQTLs showed decreasing effects on C/EBPβ binding with mechanism, whereas long-range allelic control decays with dis- increasing distance (linear regression p = 2.2 × 10−16), while no tance and distal QTLs exhibit increased cell type specificity30. distance dependence on CTCF binding was observed (p = 0.113) (Fig. 2a, b). Indeed, significant C/EBPβ QTLs (p ≤ 1 × 10−5, n = PU.1 tfQTLs associate with local chromatin state. Given that 92) were located proximal to the lead PU.1 tfQTL SNP, whereas PU.1 is a pioneer factor, we expect the allelic effects of its binding shared significant CTCF tfQTLs (p ≤ 1 × 10−5, n = 135) were to influence the local chromatin state. To verify this, we detected located distally (Fig. 2c). These results confirm that PU.1 and C/ histone modification-associated QTLs (hQTLs) for the active EBPβ bind collaboratively and suggest long-range genetic effects H3K4me3 and repressive H3K27me3 histone marks in neu- on CTCF binding. trophils and overlapped these data with PU.1 tfQTLs, identifying 621 and 367 shared PU.1 tfQTL/hQTL variants for these two Distal PU.1 binding variants are predominantly cell type spe- marks, respectively (Supplementary Data 2). Variant effects on cific. We asked to what extent the identified PU.1 tfQTLs were PU.1 binding and H3K4me3 tended to be in the same direction, NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 Fig. 2 Effect of PU.1 SNPs on second transcription factor binding. a Mean effect size (95% confidence intervals) for association of proximal sentinel PU.1 SNPs with the nearest C/EBPβ (light blue) and CTCF (red) binding site. The effect size decreases with distance for C/EBPβ (linear model p < 2.2 × 10−16) but not for CTCF (linear model p = 0.113). Beneath: bar chart of number of peaks included in each distance bin. b Genome browser shot of an illustrative example of tfQTL, where SNP rs8057431 (dashed line) alters a PU.1 motif and is associated with a disruption in binding of both PU.1 and C/EBPβ. With signal box plots created from all individuals segregated by genotype. Box plots show the medians (centre lines) and the 25th and 75th percentiles (box edges), with whiskers extending to 1.5 times the interquartile range. p values were obtained by fitting linear mixed models implemented in LIMIX. n = the number of independent donors. c Density plot displaying distance of lead proximal PU.1 SNP to the nearest shared association (p < 10−5) C/EBPβ (light blue) or CTCF (red). d Transcription factor binding intensity for PU.1 (y axis Log RPMs) at shared PU.1 and C/EBPβ tfQTL from five matched individuals in two cell types Neutrophils (left panel) and Monocytes (right panel). Binding sites were segregated by donor genotype (x-axis). tfQTL are categorised as either proximal or distal (33 distal and 43 proximal) and linear models were fitted separately for proximal and distal sites. and in the opposite direction for PU.1 and H3K27me3 (Fig. 3a), functions in lymphoid cells (Supplementary Fig. 4c)32. Notably, consistent with the predominantly activating role of PU.1 at these H3K27ac peaks co-occupied by PU.1 and C/EBPβ displayed loci. However, a significant minority (n = 131) of PU.1 tfQTLs greater hQTL effect sizes compared to peaks bound by PU.1 alone showed the same direction of effect on PU.1 binding and the (Welch two-sided t-test p = 7.28 × 10−7; Fig. 3b), suggesting levels of the H3K27me3 repressive mark, suggesting a role for stronger genetic effects for enhancers at co-occupied sites in PU.1 acting as a repressor23. We additionally took advantage of neutrophils. When compared with monocytes, neutrophil-specific the previously published hQTL data for the enhancer-associated PU.1 and C/EBPβ binding associated with neutrophil-specific histone marks H3K4me1 and H3K27ac14. In total, 808 H3K4me1 chromatin activity (Supplementary Figs. 5 and 6). Collectively, and 946 H3K27ac lead hQTL SNPs overlapped (r2 ≥ 0.8) PU.1 these results show that genetically determined differences in PU.1 tfQTLs. Using the pi1 statistic31, we found evidence of sharing binding associate with differential chromatin states in a tissue- between PU.1 tfQTLs with hQTLs in both neutrophils and specific manner. monocytes (pi1(H3K27ac) = 0.73-0.76, and pi1(H3K4me1) = We next assessed the distance between the PU.1 and histone 0.76-0.80). Sharing between neutrophil PU.1 tfQTLs and hQTLs mark peaks for each shared tfQTL-hQTL genetic association. As detected in CD4 naive T cells was lower (pi1H3K27ac = 0.36-0.72, previously observed in lymphoblastoid cell lines (LCLs)11, there and pi1H3K4me1 = 0.30-0.79; Supplementary Fig. 4b), compatible was a pronounced bimodal distribution of distances between with PU.1 being expressed at a lower level and having distinct PU.1 binding peaks and the locations of H3K27ac and H3K4me3 4 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 ARTICLE Fig. 3 Neutrophil PU.1 tfQTLs and their association with chromatin state. a Bar plot displaying the number of shared associations (r2 ≥ 0.8) between tfQTL and histone QTL. Bar plot is split into two categories, those cases where the beta between phenotypes are positively correlated, i.e., an allele leads to a gain (an increase in signal) in both phenotypes. Or negatively correlated, where an allele associates with a gain in one phenotype and a reduction in the second phenotype. The Pearson correlation between the betas is shown beneath. b Boxplot of absolute beta for H3K27ac neutrophil QTL (no significance threshold) for proximal lead PU.1 SNPs, differentiating H3K27ac regions that are or are not marked by C/EBPβ and/or CTCF binding. Box plots show the medians (centre lines) and the 25th and 75th percentiles (box edges), with whiskers extending to 1.5 times the interquartile range. n = the number of PU.1 tfQTL that intersect H3K27ac hQTL. p value displayed above is a Welch two-sided t-test. c Density distribution plot of log distance between lead PU.1 tfQTL SNPs and shared (r2 ≥ 0.8) histone QTL in neutrophils. marks (Fig. 3c), with around a half of PU.1 peaks localising to less previously implicated in the pathogenesis of inflammatory bowel than 1 kb away from the respective H3K27ac and H3K4me3 disease (IBD)35 (Supplementary Fig. 8a–c). peaks, and others mapping 10–100 kb from them. Given that Distal eQTL signals do not discriminate between direct and H3K4me3 is associated with active promoters, this observation indirect effects. Therefore, to investigate the direct role of PU.1 in highlights the potential long-range regulatory effects of PU.1 long-range gene control, we again turned to our PCHi-C data. binding to distal DNA elements on promoter activity. PIRs enriched in neutrophil specific PU.1-binding sites (pre- sented in Supplementary Fig. 6a) and enhancer-associated H3K4me1/H3K27ac marks were positively correlated with the Allelic effects of PU.1 on enhancer-promoter contacts. To level of expression of the genes they contacted in neutrophils further explore the long-range effects of PU.1 binding, we utilised (Fig. 4b), similarly as shown for other cell types36. In contrast, our high-resolution PCHi-C data in neutrophils and monocytes CTCF binding at PIRs did not correlate with target gene from six donors in total. We detected ~295,000 Promoter Inter- expression (Fig. 4b), as expected given the constitutive nature acting Regions (PIRs) in total across neutrophils and monocytes of many CTCF-mediated chromosomal interactions. In addition, (CHiCAGO score > 5)33, ~110,000 of which were detectable in we observed no association with monocyte-specific PU.1 binding, both cell types. We next explored the extent to which genetic promoter interactions and gene expression levels (Fig. 4b). variation affecting PU.1 binding associates with differential Finally, we compared the effects of distal PU.1 tfQTLs/eQTLs promoter-enhancer contacts. We employed an allele-specific (>25 kb from TSS) mapping within and outside PIRs. PU.1 strategy (“Methods”) to identify heterozygous sites within PIRs tfQTLs mapping to PIRs showed significantly larger effects on the that exhibited allelic imbalance at PCHi-C contacts (Supple- expression of the genes they contacted compared with distance- mentary Fig. 7a-b). We found an enrichment of tfQTLs within matched SNPs that did not map to a PIR (two-sided Fisher’s exact the ~14,000 heterozygous SNPs within PIRs that displayed evi- test p < 2 × 10−16; Fig. 4c), demonstrating that the physical dence of allelic bias in both neutrophils and monocytes (Fig. 4a). enhancer-promoter contacts mediate the transcriptional effects of Notably, the same was true for the hQTLs including the PU.1 binding. Polycomb-associated inhibitory mark H3K27me3, consistent with An example of a PU.1 tfQTL SNP showing allelic imbalance a role of Polycomb repressive complexes in shaping regulatory affecting promoter-enhancer connectivity and gene expres- chromatin architecture34 (Supplementary Fig. 7c–e). sion is rs519989, an eQTL of the LRRC8C gene (Fig. 4d–g). LRRC8C encodes a volume-regulated anion channel subunit upregulated during terminal differentiation of neutrophils37. This Effects of PU.1 tfQTLs on distal gene expression. To assess the and other loci (including Supplementary Fig. 8c) demonstrate effects of PU.1 tfQTLs on distal gene expression, we initially took coordinated genetic influences on PU.1 binding, chromatin advantage of expression QTL (eQTL) data generated from the activity and the formation of promoter interactions in the same individuals14, detecting 609 shared PU.1 tfQTLs/eQTLs regulation of neutrophil gene expression. ([LD] r2 ≥ 0.8), with a median distance from tfQTL to tran- scription start site of 56.5 kb. Notably, for 128 of 905 eQTL genes PU.1 differential binding underpins disease-associated var- PU.1 binding was negatively correlated with gene expression, iants. Genome-wide association studies have been used to iden- reinforcing the role of PU.1 as a dual transcriptional activator and tify regions of the genome associated with phenotypic traits and repressor23. An interesting example of a repressive PU.1 tfQTL disease. In addition, TF-binding maps and epigenetic profiling was rs2149092, associated with RNASET2 gene expression, and has been shown to assist in pinpointing the causal variants within NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 Fig. 4 tfQTLs perturb gene expression through altered chromatin state. a Enrichment of significant tfQTLs (PU.1 and CTCF; Fisher’s exact test p < 1 × 10−5) in PIRs of both neutrophils and monocytes. The bars represent the mean and the error bars the 95% confidence interval. b Heat map showing enrichment of transcription factor or histone modified regions intersecting PIRs, whereby PIRs were ranked into four bins based on the gene expression of the connected baited genes in neutrophils. c Density plot of gene expression QTL Beta value for neutrophil PU.1 SNPs within PIRs (navy) versus distance-matched significant SNPs not in PIRs (grey) (two-sided Fisher’s exact test p < 2 × 10−16), and distribution of beta values for SNPs within <25 kb of transcription start site (purple). The SNPs that are not in PIRs are also significant PU.1 tfQTLs and eQTLs (linear model p < 1 × 10−5 cut off). d CHiCAGO scores for the PIR at tfQTL, segregated by donor genotype for rs519989. n the number of individual donors. e Signal box plots with donors separated by genotype for rs519989 for five molecular traits, f boxplots displaying RNA level for LRRC8C gene segregated by donor genotype for SNP rs519989. n = the number of individual donors. Box plots show the medians (centre lines) and the 25th and 75th percentiles (box edges), with whiskers extending to 1.5 times the interquartile range. p values were obtained by fitting linear mixed models implemented in LIMIX. g Genome browser view of region around LRRC8C gene, QTL regions for each molecular trait are highlighted. Dashed line depicts the location of rs519989. these regions and implicated mechanisms associated with com- surprisingly, not at those associated with infections such Clos- plex traits. To look for evidence of a functional relationship tridium difficile (Supplementary Fig. 9). between our identified transcription factor binding sites and the To formally test the overlap of PU.1 tfQTLs and haematolo- shaping of myeloid traits, we utilised available summary statistics gical traits and diseases we accessed summary statistics from the for monocyte and neutrophil counts38. Using GARFIELD aforementioned GWAS on FBC traits38 and published auto- (GWAS Analysis of Regulatory or Functional Information immune diseases40–45. Applying colocalisation analysis46,47 Enrichment with LD correction)39, we observed enrichment of revealed 26 proximal and 53 distal tfQTLs that shared a genetic variants associated with full blood count (FBC) measure- signal (posterior probability [PP] > 0.9) with at least one GWAS ments within PU.1, C/EBPβ and, interestingly, CTCF bound locus (Fig. 5b, Supplementary Data 4). We detected an over- regions of the neutrophil genome (Fig. 5a). This analysis is sug- lapping signal at 4 loci (rs3784789, rs13089544, rs9262174 and gestive of a relationship between regulatory regions of the neu- rs13035725) between tfQTL, FBC and disease GWAS, 3 of which trophil genome and genetic variants associated with myeloid could be attributed to a lead SNP proximal to a PU.1-binding site. traits. To further investigate whether the identified tfQTLs To determine the putative target genes underpinning PU.1- potentially influence human health, we took advantage of the mediated disease associations, we integrated PCHi-C and eQTL diverse range of reported traits in 361,194 individuals from the data in neutrophils. Overall, 27 high-confidence target genes inter- UK Biobank study. We obtained summary statistics from 24 acting with eQTL loci colocalised with GWAS summary statistics GWAS for autoimmune, allergy and infection associated traits (Supplementary Data 5). generated from the UK Biobank resource (http://www.nealelab.is/ An example of a PU.1 tfQTL associated with myeloid traits is the uk-biobank/). From these we observed inflation of p values for rs791357_C variant associated with decreased neutrophil and PU.1 tfQTLs within GWAS for hay fever, asthma, rheumatoid monocyte cell counts. PCHi-C data show that this region is highly arthritis (RA) and ulcerative colitis (UC), but perhaps somewhat connected to the CPEB4 gene in both neutrophils and monocytes 6 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 ARTICLE Fig. 5 tfQTLs influence cellular phenotype and disease. a Circos plots displaying fold-enrichment of GWAS loci for neutrophil count and monocyte count within neutrophil regions marked by TF binding, modified histones and PU.1 and C/EBPβ specific binding sites in monocytes. Radial grid lines for published GWAS p values, asterisk denotes significance of enrichment for annotation tested at each GWAS p value cut off. b Bar plot displaying the number of GWAS loci that colocalise with PU.1 tfQTL for both full blood count phenotypes (left) and autoimmune disease (right). c Example of colocalised signal for sentinel SNP rs791357 tfQTL with a shared association for both neutrophil and monocyte count traits. Manhattan plots showing −log10(P) value for shared SNPs from published GWAS for neutrophil and monocyte counts obtained by fitting linear models, and PU.1 tfQTL (navy) obtained from LIMIX. With genome browser visualisation of locus for CTCF, C/EBPβ and PU.1 binding shown beneath. The top associated peak is highlighted by the shaded area. In addition, PCHi-C data show that this region is highly connected to the enhancer region in both neutrophils (blue) and monocytes (green). d Boxplot for TF and RNA signal segregated by donor genotype. PU.1 (navy), CEBPβ (light blue), CPEB4 gene expression neutrophil (navy) and CPEB4 gene expression monocyte (olive green). Box plots show the medians (centre lines) and the twenty-fifth and seventy-fifth percentiles (box edges), with whiskers extending to 1.5 times the interquartile range. p values were obtained by fitting linear mixed models implemented in LIMIX. n = the number of individual donors. (Fig. 5c-d). CPEB4 is a cytoplasmic polyadenylation protein which promoter of the RNASET2 gene, with the expression of binds to a recognition sequence in the poly-A tail of mRNAs to this gene itself (p = 5.75 × 10−61) and with a CTCF binding regulate translation48. The SNP rs791357 is a proximal tfQTL for site (p = 8.82 × 10−6) in the gene body. Interestingly, there was an PU.1 (p = 9.05 × 10−21) and C/EBPβ (p = 1.963 × 10−9), an hQTL association with H3K27me3 (p = 3.71 × 10−22) at the promoter for H3K4me3 (p = 1.98 × 10−17) and H3K27ac (p = 1.41×10−33) region and with expression of the RPS6KA2 gene (p = 1.72 × and an eQTL for CPEB4 (p = 1.16 × 10−30) in neutrophils. Similar 10−23), which interacts with the lead SNP. The RPS6KA2 gene sharing of PU.1 and C/EBPβ binding sites was observed in has been implicated as an IBD risk-associated gene in an monocytes with hQTL (H3K27ac p = 8.14 × 10−26) and eQTL for epigenome-wide association study examining differentially CPEB4 (p = 2.55 × 10−18). methylated regions in IBD patients from whole blood49. Here The aforementioned RNASET2 gene (Supplementary Fig. 8) we identify a genetic link suggesting that PU.1 is playing a pivotal has been posited previously as an IBD risk-associated gene in role in the underlying pathogenicity of IBD at this locus. CD4+ T cells35. Moreover, it was identified as a target of TL1A, a We further observed several examples of PU.1-mediated mediator of inflammation in the gut. In our neutrophil QTL data gene regulation that may point to neutrophils as being we identified rs2149092 as the nearest SNP to the PU.1 the effector cell type in the pathogenesis of inflammatory disease. tfQTL (chr6:167371914-167373868), which is in LD(r² = 1) with The GWAS colocalised loci include five shared tfQTL/eQTL in the lead reported GWAS SNP rs1819333. We confirmed a neutrophils without a significant eQTL in monocytes (rs6658646, colocalisation between this PU.1 tfQTL and GWAS for IBD and rs3768421, rs77120970, rs16895831 and rs3817465). The latter is Crohn’s disease (CD) ([PP] = 1). In addition, this tfQTL has a an eQTL for the autoimmune associated gene ILRAP18 (p = significant association with H3K4me3 (p = 6.11 × 10−17) at the 2.16 × 10−27 in neutrophils and p = 0.0457 in monocytes). In NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 addition we identified loci colocalised with genes that have been binding to understand the function of GWAS variants13. Many of previously implicated in neutrophil biology. For example CXCR2 these studies, including multi-individual PU.1-binding is involved in neutrophil migration and has been implicated in profiling11,56, were performed in LCLs. However, using disease- inflammatory disease, through the regulation of neutrophil relevant primary cells provides appropriate cellular context and extracellular trap formation in patients suffering from chronic obviates in vitro artefacts, which is particularly important for obstructive pulmonary disease50. We observed a colocalised neutrophils that are known to have pronounced differences in signal for PU.1 tfQTL, rs13035725 (p = 2.44 × 10−10) with IBD gene expression patterns from the established cell culture and UC. The same SNP is also significantly associated with a models57. Inclusion of a second primary myeloid cell type in our CTCF tfQTL (p = 5.52 × 10−6) and expression of several genes in analysis (monocytes) has enabled pinpointing cell-type-specific neutrophils, including CXCR2 (p = 4.25 × 10−13), ARPC2 (p = effects of genetic variation, and revealing that they preferentially 3.70 × 10−20), AAMP (p = 1.02 × 10−8) and CXCR1 (p = 1.23 × affect distal PU.1 tfQTLs. 10−6), all of which were weakly or not significant in monocytes. The fact that large numbers of GWAS variants tag distal The region also displays high connectivity in the neutrophil enhancers makes the identification of the target genes under- PCHi-C data (CHiCAGO score > 5), but is depleted of significant pinning GWAS associations a challenging task. Insights from interactions in monocytes (Supplementary Fig. 10a, b). expression quantitative trait loci (eQTLs) and 3D chromosomal The PLCL1 locus is another interesting example. The region contacts between enhancers and promoters are proving to be harbours variants associated with susceptibility to several diseases instrumental for GWAS gene assignment16,58. Here we used including CD51, systemic lupus erythematosus52 and allergy53. neutrophil eQTL data from the BLUEPRINT project14 and high- PLCL1 encodes the phospholipase C epsilon or phospholipase C resolution PCHi-C data in primary neutrophils to link PU.1 like I (inactive) signalling protein that has been shown to be tfQTLs and the GWAS variants they underpin with their target involved in receptor turnover but also inhibiting integrin genes. PCHi-C data from six individual donors, which emcom- activity54, suggesting a role in the regulation of cell trafficking. passes the most detailed datasets to date, also enabled us to We detected a colocalised signal with our PU.1 tfQTL and GWAS identify the potential allelic effects of promoter interactions for CD, with our lead tfQTL rs9712275 (p = 2.11 × 10−18) being associated with PU.1 differential binding. This analysis extends in high LD (r2 > 0.9) with the reported GWAS variant rs6738825. existing multi-individual studies using conventional Hi-C to In this case, the tfQTL contacts the PLCL1 gene promoter in both high-resolution promoter interaction analysis59,60. Both our and neutrophils and monocytes and we observed evidence of allelic the previous published analyses, however, used relatively small imbalance in both cell types (Supplementary Fig. 10c, d). The lead numbers of individuals for chromosomal architecture profiling, SNP rs9712275 is an eQTL for the PLCL1 gene in both and further robust associations will likely be detected by future neutrophils (p = 2.86 × 10−27) and monocytes (p = 9.14 × 10−36). higher-powered studies. Finally, we used CATO (Contextual Analysis of TF PU.1 has long been known as the key regulator of myeloid Occupancy)47 to identify collaborating factors that may be development61. However, the role of this factor in mature neu- involved in mediating these traits at shared PU.1 tfQTL/GWAS trophils has remained less well explored. A recent study by loci. Colocalising SNPs were shown to affect binding recognition Fischer et al. has shown that deletion of PU.1 in mouse mature motifs for several PU.1 binding partners, including C/EBP, AP-1, neutrophils affects the expression of a multitude of genes involved ETS, CTF/NF-1, ATF/CREB and RUNX (Supplementary Fig. 11). in immune response and results in a hyperinflammatory phe- These results highlight the likely role of PU.1 and its partners in notype following infection23. Our findings that human neutrophil mediating the functional effects of GWAS variants in neutrophils. PU.1 tfQTLs overlap eQTLs and autoimmune GWAS traits Jointly, these findings demonstrate how genetically determined provide further evidence for the role of this TF in controlling the variation in PU.1 binding, and potentially other partners, transcriptional state of mature neutrophils and point to potential underpins disease traits via altered chromatin states at regulatory clinical implications of its aberrant binding. regions and long-range chromosomal contacts. Our observation that PU.1 binding associates with both active and repressive chromatin marks, as well as both directly and inversely correlates with the effects of overlapping eQTLs, cor- Discussion roborates the evidence from Fischer et al. of the importance of In this study, we have leveraged common genetic variation as PU.1 as a transcriptional repressor in the inflammatory context. “nature’s perturbation screen” to study the functional implica- Our results are also consistent with the established role of PU.1’s tions of the binding of key myeloid TF PU.1 in human primary upstream regulator RUNX3 in autoimmune diseases such as neutrophils. We have identified thousands of DNA variants psoriatic arthritis62. TFs commonly act cooperatively63, and as associated with altered PU.1 binding (both in cis and in trans), expected, many PU.1 binding variants have shared effects on the and show that they also associate with downstream transcrip- binding of its partner TF C/EBP. Cooperativity likely underpins tional regulatory events, including local chromatin state, long- the initially puzzling findings from complementary studies of range enhancer-promoter contacts and gene expression. We other TFs, which showed that the majority of TF-binding variants further show that PU.1-binding variants underpin multiple are located outside the conventional TF motif64. Notably, the GWAS variants for healthy and pathological immune traits, abundance of distal PU.1 binding variants, consistent with pre- including autoimmunity, implicating PU.1 and neutrophils in the vious findings11, suggests that such cooperativity may extend over aetiology of these processes. considerable genomic distances, potentially through DNA GWAS have been deployed successfully over the last 15 years looping30. The shared genetic effects that we observe on the to identify thousands of genetic loci associated with complex binding of PU.1 and the key chromosomal architectural protein traits and disease. To facilitate the utility of this endeavour into CTCF, in many cases with respect to sites mapping distally to tangible improvements in human health there remains a necessity each other, support this view. to identify disease-relevant cell types, biological pathways and Over the recent years, neutrophils have emerged as important target genes. The epigenetic profiling of cell types implicated in regulators of autoimmune response. For example, neutrophils are disease aids in the mapping of disease causative genetic the most abundant cell type found in the synovial fluid between variants10,55. Our study adds to the emerging array of analyses joints of RA patients65. In IBD, neutrophils are known to damage using multi-individual readouts of the chromatin state and TF the mucosal epithelium and release monocyte chemoattractants, 8 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 ARTICLE leading to an extended inflammatory episode in the gut66. Our (Kapa Biosystems KK2602), 18 cycles of PCR followed by a 0.7:1 Ampure XP findings that GWAS variants for autoimmune disease suscept- clean-up. Antibodies for H3K4me3 (C15410003, lot A5051-001D), H3K27me3 (C15410195, lot A1811-001D), CTCF (C15410210, lot A2359-0010) were obtained ibility affect the molecular phenotypes in primary neutrophils from Diagenode, Liege, Belgium. Antibodies for PU.1 (sc-352x, lot B2415 and sc- reinforce the role of this cell type in autoimmunity and suggest 22805x, lot D0609) and C/EBPβ (sc-150x, lot G1814) were obtained from Santa that neutrophils may, at least in some cases, play a driver role in Cruz Biotechnology. the development of these conditions as opposed to merely being “support players” in autoimmune response. Data processing and peak calling. ChIP libraries were sequenced using Illumina Beyond the context of the present study, we release a large HiSeq 2000 and HiSeq 2500 at 50 bp single end reads. Sequenced reads were resource of molecular phenotypes associated with non-coding aligned to reference genome (hg19) using BWA (bwa aln –q 15). Duplicate reads were marked using Picard MarkDuplicates (v1.103). Reads with mapping quality variants in primary immune cells that we envisage to be useful for less than 15 (MAPQ 15) were removed (SAMtools v0.1.18). The fragment size L for the community to further decipher the mechanistic underpinnings each aligned bam was estimated using PhantomPeakQualTools vr18, which uses of immune gene control and variant-disease associations, poten- cross correlation of binned read counts between forward and reverse strands. To tially used in conjunction with the emerging data on in vitro identify highly enriched genomic regions, we used MACS269 (v2.0.10.20131216, standard options) for peak calling with the estimated fragment size from Phan- mutagenesis and CRISPR-targeted chromatin perturbation67. tomPeakQualTools (--shiftsize = half fragment size), with narrow for PU.1, C/ EBPβ, CTCF, H3K4me3 and broad flags set for H3K27me3. For background Methods control ChIP input was created from merging random selected samples. Reads Peripheral adult blood collection. Data generated in this study used donor from 4 pools of 12 individuals for neutrophil input and 2 pools of 6 individuals for samples which were collected as part of the previously described study14. Blood was monocytes (Supplementary Data 1). Significant peaks were selected to be at 1% obtained from donors who were members of the NIHR Cambridge BioResource FDR or less. (http://www.cambridgebioresource.org.uk/) with informed consent (Ethical over- sight by NRES Committee East of England, Hertfordshire, Study title: A Blueprint Data quality. We removed ChIP samples that had a relative strand correlation of Blood Cells, REC 12/EE/0040) at the NHS Blood and Transplant, Cambridge. (RSC) < 0.8 and normalised strand correlation (NSC) < 1.0570. We defined high- Donors were on average 55 years old (range 20–75 years old), with 46% of donors confidence data from those ChIP with RSC > 0.8 and NSC > 1.05. Otherwise, we being male. A unit of whole blood (475 ml) was collected in 3.2% sodium citrate. used genome browser tracks to confirm visually a good ChIP and include it in the An aliquot of this sample was collected in EDTA for genomic DNA purification. A final data set. Supplementary Fig. 1 and Supplementary Data 1 shows quality FBC for all donors was obtained from an EDTA blood sample, collected in parallel control metrics and corresponding principal components, showing no batch effects with the whole blood unit, using a Sysmex Haematological analyser. The level of after PEER correction using K = 10 factors. C-reactive protein (CRP), an inflammatory marker, was also measured in the sera of all individuals. All donors used for the collection had FBC and CRP parameters within the normal healthy range. Blood was processed within 4 h of collection. Normalised read count in the reference peak sets. Consensus peak sets were constructed using dba.peakset function within DiffBind R (ver3.4.0) package. http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/ Isolation of cell subsets. Samples were as those as described14. To obtain pure DiffBind.pdf71. For PU.1, H3K4me3 and H3K27me3, we set the minimum number samples of ‘classical’ monocytes (CD14+CD16−) and neutrophils (CD66b+CD16+) of samples for a peak to be included in consensus to 3, for C/EBPβ, CTCF and we implemented a multi-step purification strategy. Whole blood was diluted 1:1 in a monocyte samples minimum was set to 2. Sex chromosomes were not included in buffer of Dulbecco’s phosphate buffered saline (PBS, Sigma) containing 13 mM the QTL analysis. The reference peak set was filtered further for read counts as sodium citrate tribasic dehydrate (Sigma) and 0.2% human serum albumin (HSA, described below. Next, we generated a quantification signal of ChIP-seq for each PAA) and separated using an isotonic Percoll gradient of 1.078 g/ml (Fisher Scien- donor. Here we only considered read counts under the peaks, as the regions outside tific). Peripheral blood mononuclear cells were collected and washed twice with peaks are more likely to be noise or background signal than true enrichment. For buffer, diluted to 25 million cells/ml and separated into two layers, a monocyte each donor, we generated a vector of log2 reads per million (log2RPM) per peak in rich layer and a lymphocyte rich layer, using a Percoll gradient of 1.066 g/ml. Cells the reference peak set by counting the number of overlapping reads under the from each layer were washed in PBS (13 mM sodium citrate and 0.2% HSA) and peaks (BEDOPS v2.4.14 bedmap-count) and normalised the counts with the total subsets purified using an antibody/magnetic bead strategy. To purify monocytes, number of reads in the library. We further filtered the reference peak set to only CD16+ cells were depleted from the monocyte rich layer using CD16 microbeads consider peaks with log2RPM > 0 in at least 50% of the donors in a given cell type, (Miltenyi) according to the manufacturer’s instructions. Cells were washed in corrected for ten PEER factors and applied quantile normalisation across donors. PBS (13 mM sodium citrate and 0.2% HSA) and CD14+ cells were positively For QTL calling with H3K27me3, two sets of summary statistics are provided on selected using CD14 microbeads (Miltenyi). To purify neutrophils, the dense layer two separate signal matrices. In the first set H3K4me3 peak annotations were used of cells from the 1.078 g/ml Percoll separation was lysed twice using an ammonium in conjunction with H3K27me3 signal to enrich for poised promoter QTLs. In the chloride buffer to remove erythrocytes. The resulting cells (including neutrophils second set broad called H3K27me3 peaks were divided into 2500 bp windows. and eosinophils) were washed and neutrophils positively selected using CD16 microbeads (Miltenyi) according to the manufacturer’s instructions. The purity of each cell preparation was assessed by multicolour FACS using conjugated antibodies Identification of PU.1 and C/EBPβ differential binding sites. We used DiffBind for CD14 (MφP9, BD Biosciences) and CD16 (B73.1 / leu11c, BD Biosciences) for version 1.12.0 with default EdgeR (3.8.3) option to identify peaks which were monocytes, CD16 (VEP13, MACS, Miltenyi) and CD66b (BIRMA 17C, IBGRL- differentially bound between neutrophils and monocytes. For each cell type we NHS) for neutrophils. Purity was on average 95% for monocytes and 98% for used the six best quality samples and their peak sets for this analysis; PU.1: NS1509, neutrophils. NS1510, NS1511, NS1514, NS1516, NS1522, NS1463, NS1464, NS1554, NS1551, NS1437 and NS1490. C/EBPβ; NS1559, NS1565, NS1563, NS1558, NS1562, NS1566, NS585, NS743, NS791, NS729, NS717, and NS793. Quality control plots ChIP-sequencing. Purified cells were fixed with 1% formaldehyde (Sigma) at a generated as part of this analysis include correlation between sample peak sets concentration of approximately 10 million cells/ml. Fixed cell preparations were (Supplementary Fig. 5a) and PCA between samples using normalised counts within washed and stored re-suspended in PBS at 4 °C prior to lysis and sonication. consensus peak set (Supplementary Fig. 5b). Differentially bound peaks were Sonication protocols were performed in a Diagenode PicoRuptor for 8 cycles of considered to those peaks present in at least three individuals with a minimum 30 s on, 30 s off in a 4 °C water cooler. Samples were checked for sonication threefold difference (FDR < 0.05) in binding signal as cut off. The heat map efficiency using the criteria of 150–500 bp, by Agilent DNA bioanalyzer. ChIP-seq visualisation of differentially bound regions Deeptools 2(Galaxy Version was carried out as previously described68 all liquid handling steps were performed 3.3.2.0.1)72. on an Agilent Bravo NGS. Protein A Dynabeads (Invitrogen) were coupled with 2.5 μg of antibody. Sonicated lysate (3–5 million cells) was then added to the bead/ TF enrichments. For determining enrichment of ChIP-seq regions of interest antibody mix and incubated at 4 °C overnight. ChIP-DNA bound beads were within PIRs we used regioneR (1.0.3)73, which performs a statistical evaluation of washed for ten repetitions in cold RIPA solution. Elution of DNA from beads at two sets of genomic regions by permutation testing. We set to 50 permutations the 65 °C for 5 h to reverse the cross linking process. Two microlitres RNase was added randomisation of genomic regions to determine the null. Differentially expressed to ChIP-DNA and incubated at 37 °C for 30 min, followed by 2 μl of Proteinase K genes and gene expression counts. Gene expression counts and list of differentially treatment at 55 °C for 1 h. 1:1.8 ratio of Ampure beads (Beckman Coulter, A63881) expressed genes were available from Ecker et al25. were added to the DNA followed by two cold 70% ethanol washes. ChIP-DNA was eluted in 50 μl elution buffer. Illumina sequencing libraries were prepared on a Beckman Fx liquid handling system. End-repair, A-tailing and paired-end adaptor QTL mapping. Cis-acting QTL mapping was done using the LIMIX package74, ligation were performed using NEBnext reagents from New England Biolabs available from github (https://github.com/PMBio/limix). We considered genetic (E6000S), with purification using a 1:1 ratio of AMPure XP to sample between each variants mapping to within 1 Mb (on each side) of each tested feature (peak), and reaction. Amplification of ChIP-DNA was performed using Kapa HiFi master mix tested their association using linear regression. Models were fit on quantile- NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 normalised PEER residuals, also including a random effect term accounting for Allele-specific analysis of PCHi-C. The genotypes of PCHi-C donors were polygenic signal and sample relatedness (as in the variance component models obtained from Cambridge Bioresource phase 4 (Illumina core exome chip). We above we used the realised relatedness matrix to capture sample relatedness). From phased the genotype using BEAGLE2 (v2.0.5)78 and imputed using Positional the linear regression we obtained the effect size and p value for each tested asso- Burrows–Wheeler Transform and Haplotype Reference Consortium (release 1.1) as ciation. To correct for multiple hypothesis testing, we performed a two-step reference panel, via the Sanger imputation service. We then filtered sites for ≥5% procedure75: first, we corrected for multiple testing across variants for each minor allele frequency, HWE p value ≥ 1 × 10−6, ≤5% sample missingness and molecular outcome using Bonferroni correction and, second, we adjusted the INFO score ≥0.8. We removed indels and only considered biallelic single nucleo- obtained p values for multiple testing across phenotypes within each layer using a tide variants. We used WASP79 to remove PCHi-C reads that are likely to be biased the Q value procedure31, considered QTLs at a significance threshold of 5% FDR. towards the reference allele. We then mapped deduplicated PCHI-C reads on each allele of each SNVs using GATK ASEReadCounter with default parameters, base quality ≥2 and mapping quality ≥15. We then filtered for heterozygous SNVs only PCHi-C. Cells were isolated as described16. One donor was used for preparing each with ≥10 read counts per site and nonzero counts in both alleles. Finally, we only PCHi-C library. In total, 12 PCHi-C libraries were prepared, six using monocytes considered heterozygous sites with allele bias of ≤40% or ≥60%, after removing and six using neutrophils. Approximately 8 × 107 cells per library were re- extreme bias of <1% or >100%. suspended in 30.625 ml of DMEM supplemented with 10% FBS, and 4.375 ml of formaldehyde was added (16% stock solution; 2% final concentration). The fixation Enrichment analysis of tfQTLs and hQTLs in PIRs. Each of these heterozygous reaction continued for 10 min at room temperature with mixing and was then SNVs was annotated based on whether they were located in a PIR and whether they quenched by the addition of 5 ml of 1 M glycine (125 mM final concentration). were significant tfQTLs (PU.1 and CTCF; p < 1 × 10−5) or significant hQTLs Cells were incubated at room temperature for 5 min and then on ice for 15 min. (H3K27me3, H3K4me3, H3K27ac; p < 1 × 10−5). Fisher’s exact tests were carried Cells were pelleted by centrifugation at 400 g for 10 min at 4 °C, and the super- out separately for each sample and for each cell type to test for enrichment of natant was discarded. The pellet was washed briefly in cold PBS, and samples were tfQTLs and hQTLs that fall into PIRs. Finally, the mean and standard deviation centrifuged again to pellet the cells. The supernatant was removed, and the cell were calculated across all samples for each cell type. In another approach, all pellets were flash frozen in liquid nitrogen and stored at −80 °C. HiC library samples were combined across both cell types. SNVs were removed if they were not generation was carried with in-nucleus ligation as described previously76. Chro- observed in at least two samples, or in one sample and in the two cell types, or if the matin was then de-crosslinked and purified by phenol:chloroform extraction. DNA allelic ratio (REF reads/ALT reads) was not consistent across the samples or cell concentration was measured using QuantiT PicoGreen (Life Technologies), and 40 types. Enrichment was tested for SNVs where at least N samples fell into a PIR and μg of DNA was sheared to an average size of 400 bp, using the manufacturer’s at least N samples carried a significant tfQTL or hQTL for increasing number of instructions (Covaris). The sheared DNA was end repaired, adenine tailed and samples N (N = 1, 2, 3, 4). double size-selected using AMPure XP beads to isolate DNA ranging from 250 to 550 bp. Ligation fragments marked by biotin were immobilised using MyOne Streptavidin C1 DynaBeads (Invitrogen) and ligated to paired-end adaptors (Illu- Enrichment of genome-wide association SNPs within ChIP-seq marked mina). The immobilised HiC libraries were amplified using PE PCR 1.0 and PE regions. To test for significant enrichment of trait associated SNPs within regions PCR 2.0 primers (Illumina) with 7 PCR amplification cycles. Capture HiC of of interest, we applied GWAS analysis of regulatory or functional information promoters was carried out with SureSelect target enrichment, using the custom enrichment with LD (GARFIELD)39. H3K27ac and H3K4me1 occupied regions in designed biotinylated RNA bait library and custom paired end blockers according neutrophils were obtained from14. Neutrophil annotations for PU.1, C/EBPβ, to the manufacturer’s instructions (Agilent Technologies). Biotinylated 120mer H3K4me3 and H3K27me3 were generated as described above. With the exception RNA baits were designed to the ends of HindIII restriction fragments overlapping that H3K27me3, regions were not chunked into 2.5 kb bins. Monocyte annotations Ensembl annotated promoters of protein coding, non-coding, antisense, snRNA, are obtained from data shown in Supplementary Fig. 6a for PU.1 and C/EBPβ. miRNA and snoRNA transcripts36. A target sequence was accepted if its GC content ranged between 25 and 65%, the sequence contained no more than two Colocalisation between diseases and molecular trait. To overlap our QTL consecutive Ns and was within 330 bp of the HindIII restriction fragment terminus. results to GWAS catalogue, we calculated the LD information based on our WGS A total of 22,076 HindIII fragments were captured, containing a total of 31,253 data using plink v1.980. For all the QTLs that either directly mapped to the GWAS annotated promoters for 18,202 protein coding and 10,929 nonprotein genes variants or in LD (r2 ≥ 0.8), we considered that the QTL variant overlapped with a according to Ensembl v75 (http://grch37.ensembl.org). After library enrichment, a GWAS signal. For the cases where we further selected six autoimmune diseases, we post capture PCR amplification step was carried out using PE PCR 1.0 and PE PCR took forward the overlapping disease variants with p value ≤5 × 10−8 in six selected 2.0 primers with 4 PCR amplification cycles. PCHi-C libraries were sequenced on studies are coeliac disease40, IBD41, including CD and UC, multiple sclerosis42, the Illumina HiSeq 2500 platform3 50-PE sequencing lanes per PCHi-C library. Type 1 diabetes43, and RA44. The associations of IBD, CD and UC in the European Sequencing reads were processed and mapped with HiCUP (version 0.5.5) and cohorts were used for this study. We also used Type 2 diabetes45 as a negative PCHi-C interactions were called using CHiCAGO (version 0.2.5) with default control. We used a Bayesian colocalization method46,81 to elucidate whether the parameters33. observed overlap between disease and molecular trait may be due to a shared genetic effect. The method calculates the PP, versus the null model of no asso- ciation, for four alternative models: a model where a region or locus contains a Genotyping check of ChIP-Seq and PCHi-C data. Identity matching for each single variant associated with either the molecular trait or disease (models 1 and 2); sample and for each analysis was performed by extracting genotypes from PCHi-C a model where a single causal variant affects association with both traits (model 3); and ChIP-seq and comparing them to SNPs from the WGS data. The first stage of or a model where two distinct associations exist (model 4). The method derives the verifying the sample identity concordance between the ChIP-seq and WGS data PP of each variant in the locus being causal one under different models, and the PP involved pre-processing the BAM files for one autosomal chromosome (chr1) to of a given locus is then the integral sum of the PPs of all variants within, with all remove PCR duplicates and reads with mapping quality score <10. The variants variants under equal prior probability to be causal. The prior for each model is were then called from the resulting BAM file using mpileup from the SAMtools computed to be one that maximises the log-likelihood function46. We acknowledge package77. The variants with QUAL < 20, DP < 5 and GQ < 5 were filtered out. the limitations of the model: it assumes one causal variant in the locus; and in the Then, we compared genotypes of the filtered variants with genotypes generated case of high LD between two causal variants the model has limited power to from WGS and imputation. The genotypes generated were considered to be from distinguish model 4 from model 3. We also note that colocalisation does not imply the same sample if the concordance rate was greater than 90%. a causal relationship between molecular trait and diseases, but may be compatible also with the same variant having independent (‘pleiotropic’) effects on molecular traits and disease. We applied a colocalisation test for each of the 1003 disease- Allele-specific analysis of TF binding. For allele-specific analysis, we used the molecular trait pairs, where the lead SNPs in both traits are in r2 ≥ 0.8. To avoid phased WGS VCF that was also utilised for QTL mapping but here we removed overlapping 2 Mb-wide genetic loci due to features in close proximity (e.g., splicing indels and only considered biallelic single nucleotide variants. We then mapped junctions, genes, histones peaks, CpGs in islands), we tested colocalisation per deduplicated ChIP-seq reads on each allele of each SNVs using GATK ASER- locus, which means that the prior model parameters were estimated using one eadCounter with default parameters, base quality ≥2 and mapping quality ≥15. We locus instead of multiple loci and hence the priors may be overestimated. then filtered for heterozygous SNVs only with ≥10 read counts per site and nonzero counts in both alleles. We required two donors meeting these read counts criteria at each site. To carry out association analysis, we used Rasqual26 with total read Reporting summary. Further information on research design is available in the Nature counts per sample as an offset parameter. Note that Rasqual uses a model that Research Reporting Summary linked to this article. corrects for reference mapping bias and genotyping errors. To correct for non- genetic confounders, we applied PCA with and without permutation on normalised Data availability read counts in log2RPM across all sites and picked the first N components whose Data generated in this study was deposited to the European Genome-phenome Archive explained variances are greater than those from permutation as covariates for under the following accession IDs: transcription factor data: EGAD00001004571; Rasqual. Finally, we only considered SNVs found within peaks to determine the H3K4me3: EGAD00001002711; H3K27me3: EGAD00001002712; PCHi-C: direct allele-specific effect on TF binding of PU.1 and CTCF in neutrophils. 10 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 ARTICLE EGAS00001001911. Data from UK Biobank GWAS graciously provided by the Neale lab, 29. Merkenschlager, M. & Odom, D. T. CTCF and cohesin: linking gene [http://www.nealelab.is/uk-biobank] UKBB GWAS Imputed v3 - File Manifest Release regulatory elements with their targets. Cell 152, 1285–1297 (2013). 20180731. All additional data are available from authors upon request. Source data are 30. Delaneau, O. et al. Chromatin three-dimensional interactions mediate genetic provided with this paper. effects on gene expression. Science 364, eaat8266 (2019). 31. Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. Received: 13 May 2019; Accepted: 17 March 2021; Proc. Natl Acad. Sci. USA 100, 9440–9445 (2003). 32. Nutt, S. L., Metcalf, D., D’Amico, A., Polli, M. & Wu, L. Dynamic regulation of PU.1 expression in multipotent hematopoietic progenitors. J. Exp. Med. 201, 221–231 (2005). 33. Cairns, J. et al. CHiCAGO: robust detection of DNA looping interactions in capture Hi-C data. Genome Biol. 17, 127 (2016). 34. Freire-Pritchett, P. et al. Global reorganisation of cis-regulatory units upon References lineage commitment of human embryonic stem cells. Elife 6, e21926 (2017). 1. Bardoel, B. W., Kenny, E. F., Sollberger, G. & Zychlinsky, A. The balancing act 35. Gonsky, R. et al. Association of Ribonuclease T2 gene polymorphisms with of neutrophils. Cell Host Microbe 15, 526–536 (2014). decreased expression and clinical characteristics of severity in Crohn’s disease. 2. Nauseef, W. M. & Borregaard, N. Neutrophils at work. Nat. Immunol. 15, Gastroenterology 153, 219–232 (2017). 602–611 (2014). 36. Mifsud, B. et al. Mapping long-range promoter contacts in human cells with 3. Nguyen, G. T., Green, E. R. & Mecsas, J. Neutrophils to the ROScue: high-resolution capture Hi-C. Nat. Genet 47, 598–606 (2015). Mechanisms of NADPH Oxidase Activation and Bacterial Resistance. Front. 37. Zhu, Y. et al. Comprehensive characterization of neutrophil genome topology. Cell Infect. Microbiol. 7, 373 (2017). Genes Dev. 31, 141–153 (2017). 4. Kaplan, M. J. & Radic, M. Neutrophil extracellular traps: double-edged swords 38. Astle, W. J. et al. The allelic landscape of human blood cell trait variation and of innate immunity. J. Immunol. 189, 2689–2695 (2012). links to common complex disease. Cell 167, 1415–1429 e19 (2016). 5. Weiss, S. J. Tissue destruction by neutrophils. N. Engl. J. Med. 320, 365–376 39. Iotchkova, V. et al. GARFIELD classifies disease-relevant genomic features (1989). through integration of functional annotations with association signals. Nat. 6. Stockley, R. A. Neutrophils and the pathogenesis of COPD. Chest 121, Genet. 51, 343–353 (2019). 151S–155S (2002). 40. Dubois, P. C. et al. Multiple common variants for celiac disease influencing 7. Segal, A. W. The role of neutrophils in the pathogenesis of Crohn’s disease. immune gene expression. Nat. Genet 42, 295–302 (2010). Eur. J. Clin. Investig. 48, e12983 (2018). 41. Liu, J. Z. et al. Association analyses identify 38 susceptibility loci for 8. Mantovani, A., Cassatella, M. A., Costantini, C. & Jaillon, S. Neutrophils in the inflammatory bowel disease and highlight shared genetic risk across activation and regulation of innate and adaptive immunity. Nat. Rev. populations. Nat. Genet. 47, 979–986 (2015). Immunol. 11, 519–531 (2011). 42. International Multiple Sclerosis Genetics, C. et al. Analysis of immune-related 9. Maurano, M. T. et al. Systematic localization of common disease-associated loci identifies 48 new susceptibility variants for multiple sclerosis. Nat. Genet. variation in regulatory DNA. Science 337, 1190–1195 (2012). 45, 1353–1360 (2013). 10. Farh, K. K. et al. Genetic and epigenetic fine mapping of causal autoimmune 43. Onengut-Gumuscu, S. et al. Fine mapping of type 1 diabetes susceptibility loci disease variants. Nature 518, 337–343 (2015). and evidence for colocalization of causal variants with lymphoid gene 11. Waszak, S. M. et al. Population variation and genetic control of modular enhancers. Nat. Genet. 47, 381–386 (2015). chromatin architecture in humans. Cell 162, 1039–1050 (2015). 44. Okada, Y. et al. Significant impact of miRNA-target gene networks on genetics 12. Grubert, F. et al. Genetic control of chromatin states in humans involves local of human complex traits. Sci. Rep. 6, 22223 (2016). and distal chromosomal interactions. Cell 162, 1051–1065 (2015). 45. Morris, A. P. et al. Large-scale association analysis provides insights into the 13. Tehranchi, A. K. et al. Pooled ChIP-Seq links variation in transcription factor genetic architecture and pathophysiology of type 2 diabetes. Nat. Genet. 44, binding to complex disease risk. Cell 165, 730–741 (2016). 981–990 (2012). 14. Chen, L. et al. Genetic drivers of epigenetic and transcriptional variation in 46. Pickrell, J. K. et al. Detection and interpretation of shared genetic influences human immune cells. Cell 167, 1398–1414 e24 (2016). on 42 human traits. Nat. Genet. 48, 709–717 (2016). 15. Consortium, G. T. et al. Genetic effects on gene expression across human 47. Maurano, M. T. et al. Large-scale identification of sequence variants tissues. Nature 550, 204–213 (2017). influencing human transcription factor occupancy in vivo. Nat. Genet. 47, 16. Javierre, B. M. et al. Lineage-specific genome architecture links enhancers and 1393–1401 (2015). non-coding disease variants to target gene promoters. Cell 167, 1369–1384 e19 48. Richter, J. D. CPEB: a life in translation. Trends Biochem. Sci. 32, 279–285 (2016). (2007). 17. Mumbach, M. R. et al. Enhancer connectome in primary human cells 49. Ventham, N. T. et al. Integrative epigenome-wide analysis demonstrates that identifies target genes of disease-associated DNA elements. Nat. Genet. 49, DNA methylation may mediate genetic risk in inflammatory bowel disease. 1602–1612 (2017). Nat. Commun. 7, 13507 (2016). 18. McKercher, S. R. et al. Targeted disruption of the PU.1 gene results in multiple 50. Pedersen, F. et al. Neutrophil extracellular trap formation is regulated by hematopoietic abnormalities. EMBO J. 15, 5647–5658 (1996). CXCR2 in COPD neutrophils. Eur. Respir. J. 51, 1700970 (2018). 19. Siwaponanan, P. et al. Reduced PU.1 expression underlies aberrant neutrophil 51. Franke, A. et al. Genome-wide meta-analysis increases to 71 the number of maturation and function in beta-thalassemia mice and patients. Blood 129, confirmed Crohn’s disease susceptibility loci. Nat. Genet. 42, 1118–1125 3087–3099 (2017). (2010). 20. Heinz, S. et al. Simple combinations of lineage-determining transcription 52. Langefeld, C. D. et al. Transancestral mapping and genetic load in systemic factors prime cis-regulatory elements required for macrophage and B cell lupus erythematosus. Nat. Commun. 8, 16021 (2017). identities. Mol. Cell 38, 576–589 (2010). 53. Hinds, D. A. et al. A genome-wide association meta-analysis of self-reported 21. Zaret, K. S. & Carroll, J. S. Pioneer transcription factors: establishing allergy identifies shared and allergy-specific susceptibility loci. Nat. Genet. 45, competence for gene expression. Genes Dev. 25, 2227–2241 (2011). 907–911 (2013). 22. Bjerregaard, M. D., Jurlander, J., Klausen, P., Borregaard, N. & Cowland, J. B. 54. Lad, Y. et al. Phospholipase C epsilon suppresses integrin activation. J. Biol. The in vivo profile of transcription factors during neutrophil differentiation in Chem. 281, 29501–29512 (2006). human bone marrow. Blood 101, 4322–4332 (2003). 55. Trynka, G. et al. Chromatin marks identify critical cell types for fine mapping 23. Fischer, J. et al. Safeguard function of PU.1 shapes the inflammatory complex trait variants. Nat. Genet. 45, 124–130 (2013). epigenome of neutrophils. Nat. Immunol. 20, 546–558 (2019). 56. Kilpinen, H. et al. Coordinated effects of sequence variation on DNA binding, 24. Rustenhoven, J. et al. PU.1 regulates Alzheimer’s disease-associated genes in chromatin structure, and transcription. Science 342, 744–747 (2013). primary human microglia. Mol. Neurodegener. 13, 44 (2018). 57. Rincon, E., Rocha-Gregg, B. L. & Collins, S. R. A map of gene expression in 25. Ecker, S. et al. Genome-wide analysis of differential transcriptional and neutrophil-like cell lines. BMC Genom. 19, 573 (2018). epigenetic variability across human immune cell types. Genome Biol. 18, 18 58. Chandra, V. et al. Promoter-interacting expression quantitative trait loci are (2017). enriched for functional genetic variants. Nat. Genet. 53, 110–119 (2021). 26. Kumasaka, N., Knights, A. J. & Gaffney, D. J. Fine-mapping cellular QTLs 59. Greenwald, W. W. et al. Subtle changes in chromatin loop contact propensity with RASQUAL and ATAC-seq. Nat. Genet 48, 206–213 (2016). are associated with differential gene regulation and expression. Nat. Commun. 27. Minderjahn, J. et al. Mechanisms governing the pioneering and redistribution 10, 1054 (2019). capabilities of the non-classical pioneer PU.1. Nat. Commun. 11, 402 (2020). 60. Gorkin, D. U. et al. Common DNA sequence variation influences 3- 28. Hosokawa, H. et al. Transcription factor PU.1 represses and activates gene dimensional conformation of the human genome. Genome Biol. 20, 255 expression in early T cells by redirecting partner transcription factor binding. (2019). Immunity 48, 1119–1134 e7 (2018). NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22548-8 61. Nerlov, C. & Graf, T. PU.1 induces myeloid lineage commitment in the NIHR or the Department of Health and Social Care. M.S. is core-funded by the multipotent hematopoietic progenitors. Genes Dev. 12, 2403–2412 (1998). Medical Research Council of the UK as an MRC Investigator (MC-A652-5QA20). B.M.J. 62. Apel, M. et al. Variants in RUNX3 contribute to susceptibility to psoriatic is funded by the FEDER/Spanish Ministry of Science and Innovation (RTI2018-094788- arthritis, exhibiting further common ground with ankylosing spondylitis. A-I00) and by La Caixa Banking Foundation Junior Leader project (LCF/BQ/PI19/ Arthritis Rheum. 65, 1224–1231 (2013). 11690001). D.R., A.D., L.C., and P. Flicek are supported by the European Molecular 63. Stefflova, K. et al. Cooperativity and rapid evolution of cobound transcription Biology Laboratory. We gratefully acknowledge the participation of all NIHR Cambridge factors in closely related mammals. Cell 154, 530–540 (2013). BioResource volunteers. We thank the Cambridge BioResource staff for their help with 64. Link, V. M. et al. Analysis of genetically diverse macrophages reveals local and volunteer recruitment. We would like to thank members of the EU Blueprint consortium domain-wide mechanisms that control transcription factor binding and including, Dirk Paul, Daniel Rico, Vera Pancaldi, Willem Ouwehand and Henk Stun- function. Cell 173, 1796–1809 e17 (2018). nenberg. In addition, Stephen Wingett from Babraham Institute Bioinformatics team. 65. Wright, H. L., Moots, R. J. & Edwards, S. W. The multifactorial role of Members of the sequencing centres of Wellcome Sanger Institute and Babraham Insti- neutrophils in rheumatoid arthritis. Nat. Rev. Rheumatol. 10, 593–601 (2014). tute. Lorenz Wernisch for advice on statistics. 66. Soehnlein, O. & Lindbom, L. Phagocyte partnership during the onset and resolution of inflammation. Nat. Rev. Immunol. 10, 427–439 (2010). 67. Nasser, J. et al. Genome-wide maps of enhancer regulation connect risk Author contributions Conceived and designed the study, S.W., B.M.J., M.S. and N.S. Performed experiments, variants to disease genes. Preprint at https://www.biorxiv.org/content/ S.W. and B.M.J. Generated experimental resources, F.B., S.F. and B.F. Performed formal 10.1101/2020.09.01.278093v1 (2020). analysis, S.W., L.V., A.L.M., K.K., L.C., Y.Y., S.E., V.I., H.E., M.T., D.R., A.D. and M.S. 68. Aldridge, S. et al. AHT-ChIP-seq: a completely automated robotic protocol for Investigation, S.W., L.V., K.W. and A.L.M. Data curation, L.V., Y.Y., H.P., D.R., A.D. and high-throughput chromatin immunoprecipitation. Genome Biol. 14, R124 L.C. Supervision and study coordination, P. Flicek., L.C., K.D., T.P., P. Fraser., M.F., M.S. (2013). and N.S. Project Administration, D.M., L.C., K.D., P.Fraser., M.F., M.S. and N.S. Per- 69. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, formed primary manuscript writing S.W., A.L.M., M.S. and N.S. R137 (2008). 70. Landt, S. G. et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res. 22, 1813–1831 (2012). Competing interests 71. G, S.R.a.B. DiffBind: differential binding analysis of ChIP-Seq peak data P. Fraser. and M.S. are co-founders of Enhanc3D Genomics Ltd. P. Flicek is a member of (Bioconductor, 2011). the scientific advisory board of Fabric Genomics, Inc., and Eagle Genomics, Ltd. The 72. Ramirez, F. et al. deepTools2: a next generation web server for deep- remaining authors declare no competing interests. sequencing data analysis. Nucleic Acids Res. 44, W160–W165 (2016). 73. Gel, B. et al. regioneR: an R/Bioconductor package for the association analysis of genomic regions based on permutation tests. Bioinformatics 32, 289–291 Additional information (2016). Supplementary information The online version contains supplementary material 74. Lippert, C., Casale, F. P., Rakitsch, B. & Stegle, O. LIMIX: genetic analysis of available at https://doi.org/10.1038/s41467-021-22548-8. multiple traits. Preprint at https://www.biorxiv.org/content/10.1101/003905v2 (2014). Correspondence and requests for materials should be addressed to B.-M.J., M.S. or N.S. 75. Battle, A. et al. Characterizing the genetic basis of transcriptome diversity Peer review information Nature Communications thanks the anonymous reviewers for through RNA-sequencing of 922 individuals. Genome Res. 24, 14–24 (2014). their contribution to the peer review of this work. Peer reviewer reports are available. 76. Nagano, T. et al. Single-cell Hi-C for genome-wide detection of chromatin interactions that occur simultaneously in a single cell. Nat. Protoc. 10, Reprints and permission information is available at http://www.nature.com/reprints 1986–2003 (2015). 77. Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 78. Loh, P. R. et al. Insights into clonal haematopoiesis from 8,342 mosaic published maps and institutional affiliations. chromosomal alterations. Nature 559, 350–355 (2018). 79. van de Geijn, B., McVicker, G., Gilad, Y. & Pritchard, J. K. WASP: allele- specific software for robust molecular quantitative trait locus discovery. Nat. Open Access This article is licensed under a Creative Commons Methods 12, 1061–1063 (2015). Attribution 4.0 International License, which permits use, sharing, 80. Purcell, S. et al. PLINK: a tool set for whole-genome association and adaptation, distribution and reproduction in any medium or format, as long as you give population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007). appropriate credit to the original author(s) and the source, provide a link to the Creative 81. Giambartolomei, C. et al. Bayesian test for colocalisation between pairs of Commons license, and indicate if changes were made. The images or other third party genetic association studies using summary statistics. PLoS Genet. 10, e1004383 material in this article are included in the article’s Creative Commons license, unless (2014). indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Acknowledgements the copyright holder. To view a copy of this license, visit http://creativecommons.org/ This work was supported by the Wellcome Trust grant reference (098051, licenses/by/4.0/. 206194 and 108749/Z/15/Z), EU-FP7 Project BLUEPRINT (HEALTH-F5-2011-282510). This research was supported by the NIHR Cambridge Biomedical Research Centre (BRC- 1215-20014). The views expressed are those of the author(s) and not necessarily those of © The Author(s) 2021 12 NATURE COMMUNICATIONS | (2021)12:2298 | https://doi.org/10.1038/s41467-021-22548-8 | www.nature.com/naturecommunications