GWAS-Informed data integration and non-coding CRISPRi screen illuminate genetic etiology of bone mineral density

Over 1,100 independent signals have been identified with genome-wide association studies (GWAS) for bone mineral density (BMD), a key risk factor for mortality-increasing fragility fractures; however, the effector gene(s) for most remain unknown. Informed by a variant-to- gene mapping strategy implicating 89 non-coding elements predicted to regulate osteoblast gene expression at BMD GWAS loci, we executed a single-cell CRISPRi screen in human fetal osteoblast 1.19 cells (hFOBs). The BMD relevance of hFOBs was supported by heritability enrichment from cross-cell type stratified LD-score regression involving 98 cell types grouped into 15 tissues. 24 genes showed perturbation in the screen, with four (ARID5B, CC2D1B, EIF4G2, and NCOA3) exhibiting consistent effects upon siRNA knockdown on three measures of osteoblast maturation and mineralization. Lastly, additional heritability enrichments, genetic correlations, and multi-trait fine-mapping revealed that many BMD GWAS signals are pleiotropic and likely mediate their effects via non-bone tissues that warrant attention in future screens. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

Low-trauma fragility fractures are a significant and common cause of increased mortality and morbidity in old age1–3. Low bone mineral density (BMD), a highly heritable4–6 and polygenic7,8 trait, is among the most important risk factors for such fractures9. This key trait has been the primary focus of genomic research into fracture etiology. Genome wide association studies (GWAS) have identified over 1,100 signals associated with BMD8, with each representing a possible therapeutic target to treat low BMD. Yet, there remains a fundamental obstacle in the conversion of BMD GWAS discoveries into new treatments, namely the identity of the underlying causal effector genes. Most GWAS loci (~90%)10–12 detect non-coding variant associations that likely confer their effects by altering the expression of nearby genes13–15, though which genes is often less than obvious. Two interrelated issues have impeded broad identification of effector genes at non-coding BMD GWAS loci: 1) the development of a highly-parallelized screening technique capable of linking GWAS variants to their causal genes and 2) the need to determine the cellular and/or tissue context(s) relevant to each locus. Recently, CRISPRi screens targeted to non-coding elements have emerged as a powerful tool to solve the first of these issues. By pairing pooled CRISPRi perturbations of non-coding regulatory elements with single-cell RNA sequencing (scRNA-seq) readouts of gene expression, these screens can scale the identification of causal mediating effector genes to hundreds of loci in parallel without the need for artificial reporter constructs16– 30. However, due to the technical requirement that any cell model used for a CRISPRi screen be easily transfected, these screens have thus far been applied in a limited number of disease contexts17,21,22,26–30, none obviously related to bone biology. Regarding the second issue, while prior efforts have provided clear examples of BMD loci operating in specific cell types, most obviously in the osteoblast lineage responsible for bone deposition31–34, in general the full range of primary cell types that function in BMD pathophysiology remains uncharacterized. The only systematic genomic assessment of this issue was limited to cell types found in scRNA-seq of mouse bone35. Stratified linkage disequilibrium score regression (S-LDSC)36 offers an opportunity to determine BMD-relevant cell types throughout the whole body using human-derived measurements, specifically by finding cell types whose genomic regulatory regions are enriched for trait heritability. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint In this work, we addressed both issues described above, providing a roadmap for elucidating effector genes across non-coding BMD GWAS loci through CRISPRi screening. Specifically, we leveraged S-LDSC to identify cell types and models relevant to BMD etiology, noting among the significant results a transfectable human osteoblast cell model, the human fetal osteoblast 1.19 cell (hFOB). We used this model to conduct a focused, pooled CRISPRi screen of 89 non- coding regions harboring putatively causal variants determined through linkage disequilibrium with GWAS sentinel SNPs and identified 24 perturbed genes. Using short interfering RNA (siRNA) knockdown, we then interrogated the roles of these genes in osteoblast differentiation and function, validating 19 with one or more significant osteoblast effects. Lastly, we corroborated additional S-LDSC heritability enrichments in metabolic and structural-tissue annotations by calculating cross-trait genetic correlations and conducting multi-trait fine- mapping. These analyses revealed that at both the genome-wide and locus-specific level, the genetic etiology of BMD relates to those of other cardiometabolic and anthropometric traits in ways that suggest a substantial proportion of BMD GWAS loci confer their effects in cell types beyond those of bone. This result has critical implications for future functional experiments designed to dissect the architecture of BMD and related disease endpoints. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

S-LDSC identifies cell types relevant to BMD To identify potential model systems in which to execute a CRISPRi screen relevant to the etiology of BMD and fracture, we utilized genome-wide functional genomics to generate evidence supporting potential cell types relevant to BMD. Specifically, we applied stratified linkage disequilibrium score regression (S-LDSC)36 to partition the heritability of BMD within active and open chromatin regions annotated across a range of metabolic and structural cell types. As open and active chromatin generally reflects functional regulatory regions, cell types with enriched epigenetic annotations are more likely to play a role in BMD determination. For our analysis of public-domain and in-house-generated datasets, we employed the largest GWAS to date of BMD8 along with chromatin immunoprecipitation sequencing (ChIP-seq) peaks for activating histone marks (H3K27ac, H3K9ac, H3K4me1, and H3K4me3) and assay for transposase-accessible chromatin using sequencing (ATAC-seq) peaks to assess open chromatin. In total, we performed analysis on 210 genomic annotations across 98 primary cell types and models, which we grouped into 15 tissue categories (Supplementary Table 1). After adjusting for multiple testing, we observed significant BMD heritability enrichment (Bonferroni adj. P < 0.05) across structural tissue annotations, namely those for osteoblasts, connective tissue, and skin (Fig. 1a, Supplementary Table 1). Additionally, annotations for several metabolic tissues – adipose, cardiovascular, central nervous system, gastrointestinal, immune cells, liver, and skeletal muscle – also showed two or more significant enrichments. We then sought to determine whether there were any differences in the tissues relevant to BMD versus fracture by repeating the S-LDSC analysis on a GWAS of fracture incidence8 where we observed a broadly similar, but overall weaker pattern of heritability enrichment across tissue annotations (Supplementary Fig. 1). Unsurprisingly, given the known key role of osteoblasts in BMD determination, among the significant enrichments for BMD were several for primary osteoblasts as well as the hFOB and BMP2-stimulated human mesenchymal stem cell models (hMSC-osteoblasts) previously employed by ourselves31,34,37,38 and others39–41 to interrogate the genetic etiology of BMD in an osteoblast context. Specifically, we observed significant enrichment in primary-osteoblast H3K4me1 (adj. P = 1.15 x 10-12), H3K4me3 (adj. P = 1.01 x 10-11), and H3K27ac (adj. P = 6.81 x 10-9) ChIP-seq peaks in addition to those for H3K27ac in differentiated hFOBs (adj. P = 3.56 x .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 10-8). There was also enrichment in ATAC-seq peaks for differentiated hFOBs (adj. P = 0.013), 3-day differentiated pediatric hMSC-osteoblasts (adj. P = 3.23 x 10-8), 6-day differentiated pediatric hMSC-osteoblasts (adj. P = 2.38 x 10-8), and 3-day differentiated adult hMSC- osteoblasts (adj. P = 2.81 x 10-10). However, we did not observe significant evidence for enrichment in the three annotations for osteoclast models included in our set (min. adj. P = 1) nor in the three annotations for monocytes (min. adj. P = 0.06), the direct precursors of osteoclasts. In fact, the only three annotations among the immune tissues to show enrichment were for hESC-derived CD56+ cultured mesoderm cells (adj. P = 3.27 x 10-6) and primary G-CSF-mobilized hematopoietic stem cells (adj. P = 0.009 and 0.033), both of which represent early progenitors capable of differentiating into many different lineages aside from the monocyte-osteoclast lineage. hFOB CRISPRi screen nominates candidate effector genes at distal BMD GWAS loci Given the enrichment of BMD heritability within the regulatory regions of primary osteoblasts and their corresponding cell models, we designed an osteoblast-focused screen leveraging the easily passaged and highly transfectable hFOB model. Our objective was to elucidate novel BMD effector genes for difficult-to-resolve GWAS signals mediated through distal regulatory effects rather than nonsynonymous coding variation or the disruption of gene promoters. Building on our prior experience with 3D genomic approaches to the elucidation of GWAS signals for other common complex traits42–45, we identified 88 such candidate signals that wholly resided in open chromatin outside of any active gene promoter but which showed physical interactions, as determined by chromatin conformation capture, with the open promoters of expressed genes in the hFOB and/or hMSC-osteoblast models (Fig. 2A, Supplementary Fig 2, Methods). After merging two pairs of closely localized signals that fell within the effective repressive range of CRISPRi and including three additional signals that our group previously reported as associated with pediatric bone accrual37, our final tally of targets screened was 89. We designed three synthetic guide RNAs (sgRNAs) for each of the 89 target sites and combined them into a pool with 27 scrambled non-targeting guides and two positive-control guides that targeted the transcription start sites of RAB1A and SYVN1 (Supplementary Table 2). This pool was incorporated into lentiviral vectors and transfected into permissive hFOBs at a low multiplicity of infection (MOI) optimized to ensure that most viable cells received one sgRNA (Methods). After differentiating and collecting the cells, we used scRNA-seq to determine both .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint the gene expression profile of each cell as well the identities of any gRNAs they contained. Ultimately, we obtained 27,383 high-quality cells following filtering, each containing a single sgRNA (Methods, Supplementary Fig. 3-4). As indication of the quality of the retained cells, we observed expression of key osteoblast marker genes46 in our cell population (Supplementary Fig. 5). All but one sgRNA was successfully transfected into one or more of the 27,383 cells. The one failed transfection was due to the given sgRNA being substantially underrepresented in the lentiviral pool used for the screen (Supplementary Table 2). The remaining sgRNAs, yielded a median of 90 cells per guide with a minimum of 9 and a maximum of 354 (Supplementary Fig. 6). We used SCEPTRE47,48 to test for significant perturbations in our screen which delivered well-calibrated results with minimal statistical inflation (Fig. 2B). Since our target selection

was agnostic to any observed or predicted effect directions, we did not know a priori whether the targets were active enhancers or repressors, and we allowed for both options in our analysis. Beginning with the positive controls, we observed one of two that significantly repressed its target TSS (RAB1A). We lacked power to detect an effect for the other positive control, which targeted SYVN1, given only ~15% of cells expressed the gene. At the 89 target sites, we found 23 instances where perturbation of a target significantly impacted the expression of a gene within 1Mb; these instances spanned 23 distinct genes and 20 of the 89 distinct sites (adj. P < 0.10; Fig. 2C, Supplementary Table 3). All these perturbations were repressive except for one involving RP11-242D8.1, which is consistent with nearly all targets acting as enhancers. Eight of the perturbed genes, for example ARID5B (Fig. 2D), were predicted to be regulated by their target site based on the observed chromatin-conformation capture physical interactions. Conversely, the remaining 15 perturbations, involving genes like NCOA3, whose expression was affected by targeting rs6090584, an intronic variant within EYA2 (Fig 2E), were not reflected in the physical interaction dataset. However, given our experience that physical interactions often indicate distal regulatory activity, we hypothesized there were likely additional true perturbations that we lacked power to detect when examining all genes within 1Mb. As such, we undertook a second targeted analysis focused on the predicted regulatory connections arising from our variant-to-gene mapping (Supplementary Table 4). From this second analysis we observed one additional significant perturbation involving FAM118A (adj. P = 0.076), which brought the total to 24 perturbed genes for 21 targets. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint siRNA Knockdown of perturbed genes reveals effects on osteoblast function (and adipogenesis) We next sought to corroborate the function of the perturbed genes in osteoblasts and validate their status as BMD modulators. To accomplish this, we used siRNA to directly knock down the candidate effector genes resulting from our screen to evaluate their loss of function on osteoblast maturation and mineralization in the hFOB and hMSC-osteoblast models. As a measure of maturation, we assayed for alkaline phosphatase (ALP) staining in both cell models, and to assess mineralization we used alizarin red S (ARS) assays in the hMSC-osteoblast model; hFOBs do not synthesize sufficient mineralized matrix for this second assay. Prior to beginning the assays, we manually reviewed the perturbed genes from the screen settling on a list of 21 to be tested in the assays that did not contain targeted exonic variants and were repressed rather than upregulated by CRISPRi (Methods). After knocking-down expression of the 21 genes in hFOBs, we observed evidence of decreased ALP expression for 17 of the 21 genes relative to scrambled siRNA controls (Benjamini- Hochberg adjusted P < 0.05; Fig. 3 top panel, Supplementary Fig. 7, Supplementary Table 5). When we completed the comparable experiment in hMSC-osteoblasts, we observed that knockdown of six of the 21 genes repressed ALP expression (Fig 3. second panel, Supplementary Fig. 8, Supplementary Table 6). This decrease in the number of significant

is likely due to differences in the experimental model setting, but overall, we observed similar mean fold-changes between the two experiments. Taken together, five of the 21 genes (ARID5B, CC2D1B, EIF4G2, FAM118A, and NCOA3) revealed reduced ALP expression consistently across both cell models, while CXCL12 was the only repressed gene to show effects in the hMSC-osteoblasts but not the hFOBs. Additionally, as a positive control in the hMSC-osteoblast assay, we chose to include EPDR1, a gene we previously validated to have effects on osteoblast function34. As before, we observed that knockdown of EPDR1 resulted in decreased ALP expression. EPDR1 repression also resulted in decreases in the mineralization observed in the ARS assay, as did the repression of six other genes (Fig. 3 third panel, Supplementary Fig. 9, Supplementary Table 7). In total, the repression of four genes (ARID5B, CC2D1B, EIF4G2 and NCOA3) showed consistent effects across all three assays, and 19 genes, all but HOXD10 and CPED1, showed effects on one or more of the osteoblastic phenotypes. These results taken together with those previously published by our team on the CPED1/ING3 locus31, suggest that .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint both HOXD10 and CPED1 represent situations where a gene without clear effects on BMD is co-regulated in osteoblasts by a GWAS-tagged regulatory element alongside a likely BMD- modulating gene, in these cases, HOXD11 and ING3 respectively. To elucidate how knockdown of the genes with observed effects disrupted osteoblast function, we also assessed whether siRNA knockdown would restrict the ability of hMSCs to differentiate along the adipocyte trajectory, an alternative lineage to the osteoblast/osteocyte path. We hypothesized that if the suppression of any genes also disrupted adipogenesis, it would indicate that those genes are involved in upstream pathways that regulate the switch between hMSC proliferation and differentiation. In these adipogenesis assays, the knockdown of ten genes impaired the adipogenic potential of the hMSCs as measured by the formation of intracellular lipid droplets (Fig. 3 bottom panel, Supplementary Fig. 10, Supplementary Table 8). Focusing on the genes that had consistent effects in the osteoblast assays, we observed that the results of their knockdown aligned with the genetic architecture observed at their loci in osteoblast and adipocyte cell models. For example, NCOA3 and CC2D1B repression impaired adipogenesis, and the target site for each intersected an ATAC-seq peak found in both osteoblast and adipocyte cell models (Fig. 2E and Supplementary Fig. 11). In contrast, the ARID5B-linked target resides within osteoblast-specific open chromatin (Fig. 2D), and the EIF4G2-linked target was observed to have an interaction with the EIF4G2 promoter only in the osteoblast hMSC- osteoblast model (Supplementary Fig. 12). However, the lack of an observed interaction in adipocytes may not be a sufficient explanation for EIF4G2’s lack of effect on adipogenesis as a similar hMSC-osteoblast-specific interaction was observed at the TCF7L1 locus (Supplementary Fig. 13), though adipogenic effects were detected for TCF7L1 siRNA knockdown. Genetic correlations and multi-trait fine-mapping indicate pleiotropic genetic architecture for BMD As mentioned above, our primary objective was to validate BMD genes at non-coding GWAS loci using a CRISPRi screening system. To do so, we leveraged the most obvious cell type with BMD heritability enrichment, i.e. the osteoblast lineage. However, given we had successfully implicated effector genes using our strategy in this initial cellular setting, we returned to consider the ramifications of our presented S-LDSC analyses described above, where we observed potential roles for several metabolic and structural cell types in BMD determination. If valid, these observations suggest that to fully characterize the entire genetic architecture at BMD .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint GWAS loci, non-coding CRISPRi screens and other orthogonal functional approaches would need to be applied across a range of cell models, and crucially beyond those traditionally considered directly relevant to bone biology. We reasoned that if these additional tissues are critical to BMD pathophysiology, then BMD should share genetic etiology with phenotypes known to be related to them. To this end, we systematically investigated how the etiology of BMD relates across 37 other anthropometric and cardiometabolic traits using publicly-available GWAS conducted in the UK Biobank8,49 (Supplementary Table 9). We first investigated the cross-trait relationships of BMD at the genome-wide level via LDSC- based genetic correlations. In total, we detected a significant correlation (Benjamini-Hochberg adj. P < 0.05) between BMD and eleven additional traits (Fig 4 and Supplementary Table 10). As expected, we observed a strong positive correlation of the trait with itself and an inverse correlation with the incidence of bone fracture. We also validated a recent result showing an inverse relationship with sex-hormone binding globulin50. Many of the remaining significant correlations were for traits related to body composition and fat including body mass index (rg = 0.07), body fat percentage (rg = 0.05), trunk fat percentage (rg = 0.04), whole-body impedance (rg = -0.06), and whole-body fat mass (rg = 0.05). These correlations align well with the observed heritability enrichments in adipose, but it should be noted that they are all derived from similar anthropometric measurements and are all correlated with one another (Supplementary Fig. 14 and Supplementary Table 11). We further investigated the shared cross-trait etiology of BMD at the locus-specific level using CAFEH (colocalization and fine-mapping in the presence of allelic heterogeneity)51 and an approximate in-sample linkage disequilibrium (LD) reference panel to conservatively fine-map and colocalize 433 independent BMD signals, which we presumed captured at least one BMD- modulating variant each. (Methods and Supplementary Tables 12-14). 123 of the BMD signals were shared with one or more other traits, with 22 of the 37 input traits mapped to one or more signals (Supplementary Fig. 15). The greatest number of signals were shared with height (42 signals); three body-composition and weight-related traits: whole-body water mass (27 signals), whole-body fat-free mass (26 signals), and basal metabolic rate (25 signals); and serum ALP (22 signals). Using hierarchical clustering to group the signals by the traits for which each was mapped, we found a cluster of signals mapped to body-composition traits that suggested a shared genetic .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint etiology (Fig. 5A). Included among this cluster was the single most pleiotropic signal which consisted of a single variant at the CCND2 locus, rs76895963 (MAF = 2.1%), and which mapped to 13 traits including BMD (Supplementary Fig. 16). There were also distinct groups of signals linked to (i) height independent of body-composition, (ii) serum ALP , (iii) blood pressure, and (iv) BMD only. Visualization with uniform manifold approximation projections (UMAP)52 yielded a similar clustering of signals, with observed clusters related to body composition, body- composition-independent height, ALP, blood pressure, and BMD only (Supplementary Fig. 17). Seeking to understand whether these clusters could reflect underlying pathways, we examined the effect directions of the signals across the traits (Supplementary Table 15). In most cases, the signals identified for any trait were split between those that had positively correlated effects on the trait and BMD and those with negatively correlated effects (Fig. 5B). A notable exception was bone fracture incidence for which all six mapped signals had effects that correlated negatively with BMD. When re-clustering the signals accounting for effect directions, we observed that the major clusters all split, indicating there could be multiple pathways relating groups of traits, some with counteracting effects (Supplementary Fig. 18). In summary, these

taken together with the S-LDSC heritability enrichments, support a complex model of BMD genetic architecture that is both pleiotropic across a large subset of GWAS loci and mediated by multiple distinct pathways and cell types. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

In this study, we demonstrate the power of non-coding CRISPRi screens to elucidate unknown biology at BMD GWAS loci implicating 24 putative BMD-modulating genes, 19 of which we later validated in vitro to affect one or more measures of osteoblast maturation or mineralization in hFOBs or hMSC-osteoblasts. Notably, we found four genes (ARID5B, CC2D1B, EIF4G2, and NCOA3) exhibiting directionally-consistent repressive effects in all three osteoblast-focused assays and showed that knockdown of two of them, CC2D1B and NCOA3, also impaired differentiation of hMSCs into terminal adipocytes which aligns with genomic evidence that these genes may impact BMD determination by acting at an earlier time point in hMSC differentiation. Lastly, we observed S-LDSC heritability enrichments, genetic correlations, and multi-trait fine- mapped signals as evidence that a plethora of metabolic and structural cell types and widespread pleiotropic inheritance are critical to BMD etiology. This underlines the future challenge of applying CRISPRi screens and other experimental techniques to the complete resolution of causal effector gene identities at BMD GWAS loci. Regarding our design, we were surprised to observe that only 9 of the 21 significantly-perturbed target sites (42.8%) were observed to regulate one of the predicted genes with whose promoters they interacted in the hFOB and hMSC-osteoblast Capture-C datasets. While this proportion was higher than the 23.9% (32/134) of significantly perturbed targets in a recent screen of K562 cells that were reflected by H3K27ac HiChIP contacts21, it failed to include three genes – KARS, TEAD4, and PRPF38A – which we previously observed to have osteoblastic effects and interact with enhancers at the included pediatric bone accrual loci37. Although this may reflect temporal activities at given loci, this discordance most likely reflects differences in sensitivities between Capture-C and CRISPRi screens and highlights the value of embracing a “confluence of evidence” to nominate putative causal genes. Fortunately, these differences in sensitivity provided the opportunity to nominate CC2D1B, a second well-supported BMD gene at the locus harboring PRPF38A. These two genes provide an interesting example of a GWAS signal tagging multiple co-regulated genes with independent effects on osteoblast function, and likely BMD. In contrast, between this study and our prior work31, we have also observed two loci with co-regulated genes where only one gene has observed osteoblastic effects, namely, HOXD10/HOXD11 and ING3/CPED1. That both patterns of co-regulation exist emphasizes the need to validate putative genes derived from CRISPRi screens or any other nomination schema in orthogonal assays. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Two of the most surprising results from this study were first, the observed metabolic and structural-tissue heritability enrichments; and second, the lack thereof in the tested osteoclast annotations. Of course, we acknowledge that not all cell types with enriched epigenetic annotations are necessarily causal for BMD. Spurious enrichments can result from correlations of epigenetic features with true causal cell types via shared regulatory pathways or linkage disequilibrium with causal variants active in distinct cell types. We are aware of a few methods that leverage expression quantitative trait loci (eQTL) to fine-map causal cell-types53–55 but were cautious to employ them as they are limited by the systematic differences between the discoverability of GWAS and eQTL signals56 and the availability of bone-cell eQTLs, which is restricted to a dataset of primary osteoblasts from surgical explants57 and a dataset of RANKL- stimulated osteoclast-like cells58. Instead, we chose to interpret our S-LDSC results principally at the lower-resolution of tissues, which should be less susceptible to spurious correlations between closely related cell types. Moreover, we have identified genetic correlations and multi-trait fine-mapped signals that support the heritability enrichments. To illustrate, the cardiovascular tissue enrichments could plausibly be linked to the small cluster of blood pressure signals, and the cluster of ALP signals coincides with the observed heritability enrichments in liver, the other major source of serum ALP besides bone59. More obviously, adipose tissue enrichments align well with positive genetic correlations between BMD and body-fat percentage, BMI, and whole-body fat mass. Similarly, the cluster of signals fine-mapped predominantly to basal metabolic rate and whole-body fat- free mass support the observed enrichments in skeletal muscle tissue. Much research has focused on the relationships between body composition and BMD. Briefly, muscle mass, strength, and bone density have generally been found to correlate positively across groups4,60–66 likely sharing a causal relationship67–69, while the association between adiposity and BMD is more complex with some evidence of positive correlations between adiposity and the BMD of weight-bearing bones64,66,67,70,71 that may attenuate or even reverse after accounting for adiposity type65,72–74, at extreme ranges of adiposity66,74,75, or in certain age and sex-based sub- populations64,67,73,76. The complexity of these prior results coincides well with our observation of effect direction heterogeneity across signals fine-mapped for both body-composition traits and BMD. Such heterogeneity could also explain the small magnitudes observed for genetic correlations with body-composition traits though more work would need to be done to confirm this hypothesis and fully elucidate the clearly multi-faceted relationship between body .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint composition, size, and BMD. Given the insight available at this point, it is only obvious that these traits are somehow related, and therefore it is logical that at least a subset of BMD GWAS loci mediate these relationships in cell types related to body composition. In contrast, the lack of BMD heritability enrichment in osteoclast annotations and those of their predecessors, monocytes, seems counterintuitive given osteoclasts’ critical role in remodeling bone and since disrupted osteoclast activity has been proposed to mediate several specific BMD GWAS loci58,77–79. However, as we calculated heritability enrichments on top of the baseline S-LDSC model which accounts for conservation and genomic regulation broadly relevant across cell types, a lack of enrichment in osteoclast annotations does not prohibit isolated BMD GWAS loci from mediating their effects in osteoclasts, it merely suggests that in general, osteoclasts are not as critical to the etiology of BMD as the other enriched cell types. In fact, there may be some orthogonal evidence for this conclusion. A recent study leveraging scRNA-seq of mouse bone and a variant of S-LDSC80 found a similar lack of heritability enrichment in osteoclasts35 and another unrelated study observed that a far larger proportion of BMD-associated SNPs were eQTLs in adipose and skeletal muscle than in osteoclasts81. Of course, other explanations for the lack of enrichment are also possible. It may be that osteoclasts are more critical to BMD in its extreme ranges, such as in osteoporosis patients, but that osteoclast activity is less relevant to BMD variation in healthy individuals, such as those that comprise most of the UK Biobank cohort. Similarly, it could be possible that the RANKL- stimulated osteoclast cell model used for the osteoclast annotations does not match the primary cell closely enough to detect enrichment, perhaps due to an insufficient duration of differentiation. However, detracting from this latter hypothesis, all three osteoclast regulatory annotations were found to strongly overlap the promoters of osteoclast marker genes82 (Methods). There are several key limitations of our work. Firstly, the CRISPRi screen was conducted in an osteoblast cell model, which may not fully reflect primary cell activity. Similarly, targeted variants with significant perturbations cannot be assumed to be causal without validation due to the potential for CRISPRi-induced heterochromatin to span linked, true causal variants or even inhibit unrelated regulatory mechanisms. For example, the targeted variant at the HOXD10/HOXD11 locus, rs847147, lies 350bp upstream of the HOXD11 promoter, and the CRISPRi heterochromatin may have expanded into the promoter region blocking HOXD11 expression independent of rs847147. Additionally, as we focused on osteoblast biology, it is .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint possible that some of the targeted GWAS signals and nominated genes may have additional BMD-relevant effects in other cell types and tissues. Moreover, the S-LDSC enrichments, genetic correlations, and multi-trait fine-mappings depend on the power of their underlying GWAS and are affected by mismeasurement or heterogeneity of the phenotypes. Such phenotypic heterogeneity may explain the fewer enrichments observed in fracture incidence relative to BMD. The computational results may also not be fully transferrable outside of individuals genetically similar to European reference populations83. Lastly, seeking to prioritize precision over recall, we tolerated a high-false negative rate in the CAFEH fine-mapping to ensure that the signals we identified and the traits to which we mapped them were well- supported. For this reason, the 433 signals we report are far fewer than the 1,103 conditionally- independent signals mapped by Morris et al. using the same GWAS8. In summary, we have identified 24 putative causal genes of which we were able to validate 19 with at least one osteoblast effect and provide strongest evidence for four: ARID5B, CC2D1B, EIF4G2, and NCOA3, thereby demonstrating the power of non-coding CRISPRi screens in relevant cell models to elucidate unknown biology and causal genes at BMD GWAS loci. We also characterized the tissues relevant to BMD etiology and corroborated them via genetic correlations and multi-trait fine-mapped signals. Jointly, these results provide a roadmap for how this powerful experimental technique may be applied to the challenging task of resolving effector gene identities at all BMD GWAS loci. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

Acquisition and processing of hFOBs for ATAC-seq and Capture-C hFOBs were purchased from ATCC and maintained in a permissive state at 33.5°C in a 1:1 mixture of Ham's F12 Medium and Dulbecco's Modified Eagle's Medium with 2.5 mM L-glutamine (without phenol red), 10% Fetal Bovine Serum (FBS), and 0.3 mg/ml G418 sulfate solution. All experiments were performed on cells lower than passage 8 and confirmed to be mycoplasma negative. Cells were differentiated by increasing culture temperature to 39.5°C and were harvested for ATAC-seq and Capture-C five days post-differentiation. Matched undifferentiated control cells were also collected at the same time. Three biological replicates of the undifferentiated and differentiated hFOBs were collected for Capture-C with a fourth replicate of the differentiated hFOBs collected for ATAC-seq. Differentiation of primary human mononuclear cells into osteoclasts for ATAC-seq Human bone marrow mononuclear cells purchased from Lonza were utilized to generate and characterize human osteoclasts with minor modifications as described by Cody et al.84 and Susa et al.85 Cells were cultured in alpha-MEM containing 10% FBS supplemented with 33 ng/ml recombinant M-CSF for 2 days before using for differentiation. For differentiation, 2 x 105 cells were seeded onto a well of a 24 well plate and cultured in differentiation medium containing 33 ng/ml M-CSF, 66 ng/ml human RANKL and 1 ng/ml TGF-beta1. Medium containing supplements were re-fed every 3-4 days for a total of 12 days after which the cells were evaluated for morphological changes and stained for tartrate-resistant acid phosphatase (TRAP) with a commercially available leukocyte acid phosphatase kit (SIGMA, Cat. 387-A). Three replicates of differentiated cells were processed to prepare samples for ATAC-seq at 0, 4, 8, and 12 days. Acquisition of pediatric hMSCs and differentiation to osteoblasts for ATAC-seq hMSCs were obtained from the surgical waste of six pediatric patients undergoing ACL reconstruction surgery at the MOTT Children's Hospital, University of Michigan. Samples were processed following the protocol we published previously for adult hMSCs31. Briefly, the bone reamings were digested with collagenase for 3 hours and plated on a 10 cm dish. Cell colonies were lifted with Trypsin-EDTA and cell lines were established. Established lines were characterized by expression of MSC markers, and additionally tested for adipocyte, osteoblast, and chondrocyte differentiation. Validated cells were used for ATAC library generation. Cells .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint were collected at 3 days and 6 days post BMP2-stimulated differentiation and at 3 days post mock stimulation. hFOB, pediatric osteoblast, and osteoclast ATAC-seq library generation Fresh hFOBs, pediatric osteoblasts, and osteoclasts were harvested via Trypsin or TrypLE, followed by a series of DPBS wash steps. 50,000 cells from each sample were pelleted at 550 × g for 5 minutes at 4 °C. The cell pellet was then resuspended in 50 μl cold lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% IGEPAL CA-630) and centrifuged immediately at 550 × g for 10 minutes at 4 °C. The nuclei were resuspended in transposition reaction mix (2x TD Buffer (Illumina Cat #FC-121–1030, Nextera), 2.5 µl Tn5 Transposase (Illumina, 20034197 Cat #FC-121–1030, Nextera) and Nuclease Free H2O) on ice and then incubated for 45 minutes at 37°C. The transposed DNA was then purified using the MinElute Kit (Qiagen), eluted with 10.5 μl elution buffer (EB), frozen and sent to the Center for Spatial and Functional Genomics at CHOP. The transposed DNA was PCR amplified and indexed using the Illumina Nextera Kit (Illumina) and NEBNext High-Fidelity 2x PCR Master Mix (NEB) for 12 cycles to generate each library. The PCR reaction was subsequently purified using AMPureXP beads (Agencourt) and libraries were paired-end sequenced on the Illumina NovaSeq 6000 platform. Human articular chondrocyte isolation Human knee articular cartilage was provided by AlloSource (Centennial, CO) from donors deemed eligible for tissue donation for research purposes. Donor eligibility was determined in accordance with American Association of Tissue Banks (AATB) and Food and Drug Administration (FDA) regulations. Tissue fragments that included subchondral bone and cartilage from both the tibial plateau and femoral condyles were surgically removed from donors (N=3) and immediately placed into pre-chilled (wet ice) wash medium composed of DMEM/F12 medium (Cytiva, #SH30023.01) containing amphotericin (1 ng/ml, Sigma, #50-175-7519), gentamycin (0.05 µg/ml, Gibco, #15750-060), and Pen/Strep (1% v/v, VWR, # K952-100ML) for transport to the laboratory. Articular cartilage was removed from the subchondral bone using a scalpel, minced into 1 mm X 1 mm fragments in a Petrie dish, rinsed three times with phosphate-buffered saline, and placed into a 50 ml conical tube containing 30 ml wash medium supplemented with fetal bovine serum (FBS, 25% v/v, Gibco, #10438-026). Chondrocytes were isolated from the cartilage using a modified version of an established method86. Specifically, after 60 minutes of gentle shaking at 37°C, medium was aspirated and replaced with 30 ml .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint digestion medium composed of DMEM/F12, 25% FBS, 0.05 µg/ml gentamycin, 1% v/v Pen/Strep, and ascorbic acid (100 µg/ml, Sigma, #A4544-25G), and supplemented with pronase (0.3 µg/ml, Roche, #10165921001). Cartilage fragments were digested in this medium for 90 minutes at 37°C with gentle shaking, followed by centrifugation at 300 X g at 4°C for 7 minutes to pellet tissues and cells, followed by resuspension in 30 ml digestion medium supplemented with collagenase II (1.2 mg/ml, Worthington, #LS004177). Cartilage was digested for 18 hours at 37°C with gentle shaking and filtered through a 70 µm strainer, with a single rinse of the tube with wash medium to recover all remaining cells and tissue, which was again strained. Strained

were centrifuged to a pellet at 300 X g for 7 minutes at 4°C, resuspended with fresh wash medium, spun once more at 300 X g for 7 minutes at 4°C to form a pellet, and then finally resuspended in DMEM/F12 containing 20% FBS, 0.05 µg/ml gentamycin, 1% v/v Pen/Strep, and 100 µg/ml ascorbic acid. Isolated suspensions of chondrocytes were counted using a hemocytometer. Human chondrocyte nucleic acid preparation for ATAC-seq Immediately following human articular chondrocyte isolation, nucleic acids were extracted for use in an ATAC-seq experiment. From each human donor (N=3), suspensions of 75,000 cells were centrifuged at 550 X g for 5 min at 4°C, resuspended in 50 µl of ice-cold PBS, and centrifuged once more at 550 X g for 5 minutes at 4 °C. The final cell pellets were resuspended in 50 µl of ice-cold lysis buffer (10 mM Tris-HCL, pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% IGEPAL CA-630) and immediately centrifuged at 550 X g for 10 minutes at 4 °C. The supernatant was discarded, cells were placed on wet ice, and the TDE1 Tagment DNA Enzyme and Buffer Kit (Illumina, #20034197) was utilized according to the manufacturer's instructions. Briefly, chondrocytes were incubated in the TDE1 reaction buffer for 45 minutes at 37°C, followed by addition of 10 µl of 3 M sodium acetate to stop the reaction. Nucleic acid purification was performed using the Qiagen MinElute Kit (#28204) following the manufacturer's instructions. Purified DNA was stored at -20°C until shipped to the CHOP Center for Spatial and Functional Genomics. Libraries were generated and sequenced in the same manner as indicated in the “hFOB, pediatric osteoblast, and osteoclast ATAC-seq library generation”

section. ATAC-seq alignment and peak calling ATAC-seq peaks were called for hFOBs, pediatric osteoblasts, osteoclasts, and chondrocytes using the ENCODE ATAC-seq pipeline87 and default settings. Briefly, this pipeline input pair-end .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint reads from the biological replicates for each cell type and aligned them to GRCh38 using bowtie288, removing any duplicate reads from the alignment. The pipeline then called narrow peaks independently for each replicate using macs289 and removed peaks in ENCODE blacklist regions (ENCFF001TDO). Quality-control metrics were checked against the ENCODE recommended standards, and any data sets found not to meet standards, namely the 0, 8, and 12-day differentiated osteoclasts, were discarded from further analysis. For the S-LDSC analysis, we used the irreproducible discovery rate (IDR)90 optimal peak sets. For target selection using the hFOB and hMSC-osteoblast ATAC-seq results, we re-ran the ENCODE pipeline, aligning this time to hg19, and used the less stringent pooled peak sets. The 4-day differentiated osteoclast peaks were also checked for overlap with the promoters of seven osteoclast marker genes: CALCR, CA2, CTSK, MMP9, SPP1, ACP5, EDNRB82. The promoters were defined as ±1kb from the GRCh38 transcription start sites as obtained from GeneCards91. All promoters except those for CALCR and EDNRB were overlapped by one or more peaks. GWAS summary statistics Summary statistics for BMD (estimated by heel quantitative ultrasound) and bone fracture were obtained from the largest GWAS of each to date8. Each of these GWAS was executed on a population of white British individuals from the UK Biobank92 determined based on genetic similarity to the 1000 Genomes GBR subpopulation83. We also obtained GWAS summary statistics for 36 traits analyzed by the Pan UKBB project in an overlapping population of individuals based on genetic similarity to the 1000 Genomes EUR superpopulation.49,83 (Supplementary Table 9). These metabolic and anthropometric phenotypes were manually selected to represent diverse areas of biology. All summary statistics were downloaded in hg19. Linkage Disequilibrium Score Regression We calculated genetic correlations and heritability enrichment in cell-type specific ATAC-seq and histone ChIP-seq peaks via cross-trait93 and stratified36 linkage disequilibrium (LD)-score regression94 respectively (v1.0.1; https://github.com/bulik/ldsc). We used the 1000 Genomes Project Phase III GRCh38 LD reference95 to calculate LD scores for the regressions and only retained SNPs included in the HapMap Project Phase 3 call set96. All GWAS summary statistics were lifted over to GRCh38 for the analyses using the UCSC LiftOver tool97 (https://genome.ucsc.edu/cgi-bin/hgLiftOver). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Stratified linkage disequilibrium score-regression (S-LDSC) heritability enrichments were calculated on top of the baseline LDSC model using 210 binary genomic annotations across 98 primary cell types and cell models. 198 of the annotations for 86 unique cell types were consolidated ChIP-seq narrow peak calls for four types of activating histone marks: H3K4me1, H3K4me3, H3K9ac, and H3K27ac. These peaks were obtained directly from the Roadmap Epigenomics Project98,99 and were supported by a minimum of 30M reads (Supplementary Table 1). These annotations were downloaded in hg19 and lifted over to GRCh38. We also included in the analysis the newly generated ATAC-seq peaks described above for pediatric hMSCs and hMSC-osteoblasts, chondrocytes, and RANKL-differentiated osteoclast models as well as reprocessed datasets obtained from the public domain. For the public domain datasets, we obtained raw FASTQ files for a H3K27ac ChIP-seq experiment in hFOBs100, a paired ATAC- seq / H3K27ac ChIP-seq experiment in RANKL-differentiated osteoclasts101, and two ATAC-seq studies previously published by our team for hMSC-osteoblasts31 and tonsillar-organoid sorted monocytes43. These FASTQs were reprocessed using standard ENCODE pipelines87 and aligned to GRCh38. The overlapping optimal peak sets were used for the reprocessed ChIP-seq datasets, and optimal IDR peaks were used for the reprocessed ATAC-seq experiments. The osteoclast peaks were also checked for overlap with the marker gene promoters. The ATAC-seq peaks overlapped all promoters except for EDNRB, and the H3K27ac ChIP-seq peaks overlapped 4 of 7 promoters (CA2, CTSK, MMP9, and ACP5). Bonferroni-adjustment102 was used to correct for multiple testing across annotations in the S- LDSC analysis, and all annotations with an adjusted P < 0.05 were deemed to have significant heritability enrichment. Genetic correlations were adjusted for multiple testing via the Benjamini- Hochberg procedure103 as part of two separate analyses. In the first, genetic correlations were calculated for all traits with BMD (Supplementary Table 10), and in the second analysis, genetic correlations were calculated between all unique pairs of the weight and impedance related traits (Supplementary Table 11). Adjusted P < 0.05 were again considered significant for these analyses. Additionally, in the first analysis, two traits (hypoglycemia and type 2 diabetes) were found to have genetic correlations with BMD that could not be estimated by LDSC. These values were ignored when correcting for multiple testing and reporting results. hFOB promoter-focused Capture-C library preparation and sequencing We followed the procedure previously published by our team for the generation and sequencing of the hFOB promoter-focused Capture-C libraries104. For this protocol, 107 fixed hFOB cells .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint were resuspended in dH2O supplemented with protease inhibitor cocktail and incubated on ice for 10 minutes twice. After setting aside 50 μl of cell suspension for pre-digestion QC, the remaining sample was divided into 6 tubes. All incubation reactions were carried out in a Thermomixer (BenchMark) shaking at 1,000 rpm. Samples were pre-digested for 1 hour at 37°C after adding 0.3% SDS, 1x NEB DpnII restriction enzyme buffer and dH2O. We then added a 1.7% solution of Triton X-100 to each tube and continued the incubation an additional hour. In the sample tubes only, we added 10 µL of DpnII (NEB, 50 U/μl) and continued the incubation until the end of the day when another 10 μl DpnII was added to each sample to digest overnight. The next morning, we added another 10 μl DpnII and incubated for a final 2−3 hours. We removed 100 µL of each digestion reaction, pooled them into two 1.5 ml tube, and set them aside for digestion efficiency QC. We heat inactivated the remaining samples at 65°C for 20 minutes, before cooling on ice for 20 minutes. We ligated digested samples overnight at 16°C with T4 DNA ligase (HC ThermoFisher, 30 U/μl) and 1X ligase buffer. The next day, we spiked an additional T4 DNA ligase into each sample and incubated another few hours. We then de-crosslinked the samples overnight at 65°C with Proteinase K (20 mg/ml, Denville Scientific) along with pre-digestion and digestion controls. The following morning, we incubated the controls and ligated samples for 30 minutes. at 37°C with RNase A (Millipore) prior to phenol/chloroform extraction, ethanol precipitation at −20°C, and centrifugation at 4°C and 3000 rpm for 45 minutes to pellet the samples, while controls were pelleted at 14,000 rpm. Pellets were washed in 70% ethanol and centrifuged again as described above. We resuspended the 3C library and control pellets in dH2O and stored both at −20°C. We measured sample concentrations by Qubit and assessed digestion and ligation efficiencies by gel electrophoresis on a 0.9% agarose gel and by quantitative PCR (SYBR green, Thermo Fisher). DNA from each 3C library (10 ug) was sheared to an average fragment size of 350bp using a QSonica Q800R (60% amplitude, 30 seconds on / 30 seconds off, 2-minute intervals for 5 total intervals) at 4°C. After shearing, DNA was purified using AMPureXP beads (Agencourt), DNA size was confirmed on a Bioanalyzer 2100 (Agilent) and DNA concentration measured via Qubit. Libraries were prepared for selection using the SureSelect XT Library Kit (Agilent) following the manufacturer protocol and once again bead purified, then checked for size and concentration as described above. One microgram of adaptor-ligated library was hybridized using the SureSelect XT capture kit (Agilent) and our custom-designed 41K promoter Capture-C .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint probe set31. After amplification and purification, we assessed the quantity and quality of the captured libraries one final time. We paired-end sequenced all promoter-focused capture-C libraries on the Illumina NovaSeq 6000 platform with 51bp read length. Analysis of hFOB Capture-C data We pre-processed paired-end reads from the three hFOB replicates using the HICUP pipeline105 (v0.5.9) aligning reads to the hg19 reference genome with bowtie288. We called significant promoter interactions with genes from the GENCODE Release 19 gene set (GRCh37.p13)106 at 1-DpnII fragment resolution using CHiCAGO107 (v1.1.8) with default parameters except for binsize set to 2500. We also called significant interactions at 4-DpnII fragment resolution by artificially grouping four consecutive DpnII fragments and inputting them into CHiCAGO using default parameters except for removeAdjacent which was set to False. We considered interactions with a CHiCAGO score > 5 at either 1-fragment or 4-fragment resolution to be significant interactions and converted significant interactions to ibed format for use in variant to gene mapping. hFOB RNA-seq Total RNA was isolated from hFOB cells using TRIzol reagent (Invitrogen) following manufacturer instructions, then purified using the Direct-zol RNA Plus Miniprep Kit (Zymol). After measuring concentration (Nanodrop, Invitrogen) and RNA integrity (RIN > 7, Bioanalyzer 2100, Agilent), RNA was depleted of rRNA using the QIAseq FastSelect RNA Removal Kit (Qiagen). RNA-seq libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (NEB) and NEBNext Multiplex Oligos for Illumina (Dual Index Primers, NEB) following standard protocols. Libraries were sequenced on an Illumina NovaSeq 6000, generating ~100 million paired-end 50 bp reads per sample. RNA-seq data were aligned to the hg19 genome with STAR v. 2.6.0c and gene counts were obtained using HTseq, with flags -f bam -r name -s reverse -t exon -m intersection-strict, to count genes from GENCODE Release 19 (GRCh37.p13) annotation plus annotation for lincRNAs and sno/miRNAs from the UCSC Table Browser (downloaded 7/7/2016). Normalized counts for the uniquely mapped read pairs were generated through the transcript per million read method with effective gene length and the resulting values were used in the computation of gene expression percentiles. CRISPRi Target Selection .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint We selected targets for the screen beginning with the list of 1,103 independent BMD signals reported by Morris et al.8 We identified LD-proxies for the signals at an r2 > 0.8 using SNiPA v3.4108 (https://snipa.org/snipa3/) and the 1000 Genomes Phase III hg19 LD reference95. 36 signals could not be mapped through SNiPA and were retained with only themselves as proxies. We mapped proxy variants to candidate effector genes in differentiated hFOBs and hMSC- osteoblasts by identifying promoter-interacting fragments that contained a proxy and overlapped an ATAC-seq peak in the same cell type. We discarded gene nominations for genes with expression of less than 1 transcript per million in the same matched cell types and any in which the implicating variant overlapped the promoter of an expressed gene. Having retained 88 signals with one or more linked gene, we found two pairs of these signals located within 1kb, the distance we assumed as the effective repressive range of CRISPRi109, and we combined them together. Additionally, we added three more targets we previously identified at genomic loci associated with pediatric bone density accrual37 bringing the total to 89 targets. Custom sgRNA Pool Target Design We designed synthetic guide RNAs (sgRNAs) for each target site using CRISPick110,111. We then input the CRISPick recommended sequences into FlashFry112 and discarded any candidate guides FlashFry flagged to have high GC content, polyT sequences, or multiple genomic targets. We then iterated over the ranked list of remaining guides for each target and selected the top three whose binding sites did not overlap any of the previously selected guides. Generation of Helper hFOBs and sgRNA Configuration Optimization. Helper hFOBs expressing the dCas9-CRISPRi-KRAB lentiviral construct under Blasticidin selection (10ug/ml,11 days) were first generated using the Sigma 10X CRISPRi Feature Barcode Optimization Kit (CRISPRI10X). Next, to test for the optimal sgRNA configuration, stocks of the helper hFOB-dCAS9-CRISPRi-KRAB cells were transduced with lentiviral pools containing one of four sgRNA capture sequence configurations and selected with Puromycin (1 ug/ml, 11 days): Capture Sequence One Stem (CS1-STEM), Capture Sequence One Three Prime (CS1-3’), Capture Sequence Two Stem (CS2-STEM), and Capture Sequence Two Three Prime (CS2-3’) (Supplementary Fig. 19). Each lentivirus pool contained sgRNA targeted to the RAB1A TSS and a Negative Control. Cells from all four configurations were subjected to 10X Genomics single cell analysis for both scRNA-seq (GEX library) and Feature Barcoded CRIPSR Capture (CRISPR Capture library). The optimal configuration was determined using both the best fraction of usable guide reads and the best -log2 fold change in RAB1A mRNA expression. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint The CS1-STEM configuration was determined to be optimal for hFOBs (Supplementary Fig. 19). Generation of CRISPRi sgRNA Pool Targeted hFOBs. To generate the hFOBs containing our custom sgRNAs, the same helper hFOB-dCas9- CRISPRi-KRAB cells were used for transduction. Cells were transduced at low MOI (0.2) plus polybrene (8 ug/ml) with a Sigma-Aldrich custom sgRNA lentiviral pool (titer = 5.3 x 108 TU/ml). We selected an MOI of 0.2 to ensure that most viable cells would contain only one sgRNA and determined the optimal titer via the recommended procedure113. Under a Poisson model and perfect selection for transfected cells, ~90% of viable cells were expected to have one sgRNA. The lentiviral vectors followed Sigma-Aldrich’s standard CRISPRi-screen construct design (Supplementary Fig. 20) and were sequenced by the manufacturer to ensure quality prior to shipping (Supplementary Table 2). On day 2 post-transduction, cells were selected with Puromycin (1 ug/ml). Transduction was confirmed at day 8 by blue fluorescent protein (BFP) and frozen for stocks on day 11. Stock hFOB-CRISPRi-KRAB-Pooled-sgRNA cells were grown in 100mm plates under Blasticidin/Puromycin selection for 2 days at 33.5°C, then differentiated for 5 days at 39.5°C. Cells were removed from plates with TrypLE, counted, and diluted to 1000 cells/ul in DPBS+1% FBS. Viability was determined to be around 90% before 160,000 cells (8 lanes of 20K cells each) were processed for both 10X Genomics scRNA-seq (GEX libraries) and Feature Barcoded CRIPSR Capture (CRISPR Capture libraries) at the CHOP Center for Applied Genomics (CAG). Both sets of libraries were sequenced as eight pools on the Illumina Novaseq 6000 system using an S2-100 flow cell. Single-Cell Processing Single-cell FASTQs were initially processed using the CellRanger pipeline (10X Genomics Cell Ranger 3.0.0)114 with default settings. We then used CellBender115 to denoise and filter the raw CellRanger outputs separately for each of the eight pools. The number of droplets and expected number of cells were visually estimated from each pool’s unique molecular identifier (UMI) curve, and we ran Cellbender with a learning rate of 0.00005, 150 training epochs, and a false- positive rate of 1%. After CellBender we used UMAP projections52 and violin plots implemented in Scanpy116 to visualize the remaining droplets for each pool and quickly recognized that all pools but the second contained two clusters of droplets, one with a high number of UMIs and genes per droplet and the other with a low level of cellular complexity. The second pool had an overall low level of complexity indicative of mostly empty droplets and was discarded. For the .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint remaining pools we used the Leiden algorithm117 to cluster the cells and retained only the high complexity cluster. The remaining droplets were then combined across pools and visualized using UMAP and violin plots. We discarded droplets with >10% mitochondrial reads and >90,000 UMI to remove suspected dying cells and doublets respectively and retained 43,144 high-quality cells. We similarly removed 8,822 cells without sgRNAs as these were either non-transfected cells or cells whose guides were not captured in the scRNA-seq. We then tested each non-targeting sgRNA for random assortment against each of the targeting guides using a Fisher Exact Test and a 2 x 2 contingency table. Results were corrected for multiple testing using the Benjamini- Hochberg procedure103 and a significance threshold of 0.05 was used. Finding 15 of 27 non- targeting guides preferentially assorted with one or more targeting guides, we discarded 5,655 remaining droplets containing more than one sgRNA, retaining 27,383 cells. Perturbation Testing and Visualizations We next used the low-MOI version48 of SCEPTRE47 v0.3.0 (https://github.com/katsevich- lab/sceptre) to test all genes within 1Mb of each target site for perturbations. This method pools sgRNAs targeted to the same variants and tests them jointly using a permutation test. We determined which genes to test for each target site by taking the full GENCODE Release 19 gene set (GRCh37.p13)106 and identifying all genes who overlapped or came within 1kb of the binding sites for any of the sgRNAs for the target and were captured in the scRNA-seq. As covariate inputs into SCEPTRE, we included the number of UMIs per cell, the number of unique genes detected per cell, the pool in which each cell was sequenced, the mitochondrial read % per cell, and the top 15 gene expression principal components (PCs). We used Seurat v5.0.1118 to calculate the PCs via the recommended procedure on the 2,000 most variable genes determined via the vst method. We assumed targeted elements could be either enhancers or repressors and allowed for both possibilities by using a two-sided test. Statistical calibration was confirmed visually from the quantile-quantile plot. For the initial perturbation analysis, we considered all tested pairs of genes and targets and corrected for multiple testing using the Benjamini-Hochberg procedure103. In a second analysis focused on predicted variant-to-gene relationships informed by Capture-C, we subset the list of tested gene-target pairs to those reflected in the variant-to-gene mappings, and recorrected .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint using the Benjamini-Hochberg procedure103. For both analyses adjusted P < 0.10 were considered significant. We visualized genomic annotations at targets found to exhibit significant perturbations with pyGenomeTracks v3.8119 (https://pygenometracks.readthedocs.io/). Included among the tracks were the hg19-aligned hFOB ATAC-seq and Capture-C annotations described above as well as previously published ATAC-seq and Capture-C datasets for hMSC-osteoblasts31 and adipocytes differentiated from hMSCs (hMSC-adipocytes)120. For visualization purposes the basic GENCODE Release 19 gene set106 was used. Perturbed genes were plotted in red and all others in blue. Selection of Gene Targets for siRNA Assays Prior to beginning functional assays, we noted we had retained a few targets in our screen that, while not intersecting gene promoters, did overlap exons and were therefore not the focus of our work. We dropped the corresponding perturbed genes – ADAT1, ADCY4, and FBXW4 – from consideration in our siRNA-based assays. For reference, the first two targets resided in occasionally-retained introns within the 5’-UTR of ADAT1 and ADCY4 respectively, and the last was a synonymous coding variant in FBXW4. We also dropped RP11-242D8.1 from inclusion in the assays believing it was a false positive result since it was located at the same locus as the well-established BMD gene, SOST, and was the only gene to show increased expression in the screen, in contrast to the generally repressive nature of CRISPRi. However, after manually reviewing the sub-significant screen results for borderline genes with strong prior biological evidence of causality, we included CALCRL, a gene that showed suggestive evidence of perturbation (P = 0.028) and is closely related to the calcitonin receptor that plays an essential role in bone biology121,122. This brought our total number of assayed genes to 21. hFOB siRNA Treatments and Alkaline Phosphatase Assay Single cell suspensions of hFOBs were seeded into 24-well plates at 45-60K cells per well and allowed to adhere overnight. Transfections were carried out the next day using ON- TARGETplus SMARTRpool siRNA purchased from Horizon Discovery (Supplementary Table 16) and Dharmafect-1 transfection reagent per the manufacturer’s protocol. The next day, growth media was replaced. The plate designated for differentiation into osteoblasts was placed at 39.5oC, while the permissive plate was kept at 33.5oC. Both plates were stained for ALP after 4 days using the Alkaline Phosphatase Staining Kit (Abcam, ab242286) following kit .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint instructions. Plates were photographed and the images were split into 8-bit RGB images using Image J software. Images within the green channel were used to enumerate integrated density values within the cell culture area for each well as previously described31. Assays were repeated six times. hMSC siRNA Treatments and Assays Following our previously published protocol for conducting siRNA assays in hMSC models31,37, we obtained primary bone-marrow derived hMSCs isolated from healthy adult donors (Supplementary Table 17) and characterized them for cell surface expression (CD166 + CD90 + CD105+/CD36-CD34-CD10-CD11b-CD45-) and tri-lineage differentiation (osteoblastic, adipogenic, and chondrogenic) potential. We achieved experimental knockdown of candidate genes using siRNAs as in the hFOB cell models. For osteoblastic differentiation, we plated 15,000 cells/cm2 in alpha-MEM consisting of 16.5% FBS, 25 µg/ml Ascorbic acid-2- phosphate, 5 mM beta-glycerophosphate and 1% insulin-transferrin-selenous acid (osteogenic media) and stimulated them the next day with recombinant human BMP2 (300 ng/ml) (R&D Systems, MN) in serum-free osteogenic media. Cells were harvested at 72 hours following BMP2 treatment for alkaline phosphatase assessment and at 8–10 days for staining with Alizarin red S. ALP and ARS assay plates were scanned on a flatbed scanner and quantified by Image J after splitting the color images into 8-bit RGB images as described above. For differentiation into hMSC-adipocytes, 30,000 cells were seeded on 24 well plates and transfected next day using Dharmafect-1. Cells were allowed to recover for 2 days and adipogenic differentiation was started using 10% FBS alpha-MEM supplemented with Indomethacin, IBMX, and Dexamethasone as described previously31. Media exchange was carried out every 3 days until staining with Oil Red O at 18-21 days. Lipid droplet accumulation was enumerated using Lionheart automated microscope in the Texas Red channel. 4X objective was used to take a montage of 25 different microscopic fields which were then stitched and quantified using the cell count feature. Representative images were taken with a 20x objective with DAPI nuclear staining for reference. Statistical Testing of siRNA Assays The effects of siRNA knockdown on functional measures were assessed by comparing wells with differentiated cells transduced with gene-targeting siRNAs to those targeted with scrambled, control siRNAs. Assays were designed to be at or near saturation in the unperturbed .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint state, so one-sided Student’s t-tests were used to assess the loss of staining upon siRNA knockdown. Prior to statistical testing, all wells were normalized relative to the differentiated control-siRNA well on the plate. For hMSC-based assays, where available, technical replicates reflecting multiple passages of donor lines were averaged together prior to using each donor measurement as an instance for statistical testing. ALP and ARS assays in hMSC-osteoblasts were both conducted in five donors, and adipogenesis assays were conducted in between 2-4 donors due to the availability of matched stocks for the donors used in the osteoblast assays. CAFEH Multi-Trait Fine-Mapping We used the multi-trait fine-mapping algorithm, CAFEH (https://github.com/karltayeb/cafeh), to fine-map and colocalize shared causal BMD signals across the 38 GWAS described above. We first defined 501 BMD-relevant loci by adapting a previously reported approach123. This method involves tiling the genome into 250kb tiles and merging all adjacent tiles with one or more significant variants at a P < 10-6 threshold. Adjacent tiles were padded by 250kb on both sides to form loci and any overlapping loci were merged. Loci with one or more genome-wide significant variants (P < 5 x 10-8) were retained for analysis and the rest discarded. We then identified which non-BMD traits to fine-map at each locus by identifying any that had one or more genome-wide significant variants (P < 5 x 10-8) within the bounds of the BMD-defined locus. For each locus, we input the variants tested in all the mapped GWAS studies into CAFEH. To execute signal fine-mapping, we downloaded published LD reference matrices calculated by the Pan-UKBB team in a subset of 421K participants from the UK Biobank who are genetically similar to individuals in the 1000 Genomes EUR superpopulation.49,83 The UK Biobank-based populations used to create the LD matrices and conduct the GWAS heavily overlap and range in size from 347K to 427K individuals (Supplementary Table 9). We limited ourselves to these datasets, though larger GWAS are available for some of the traits, precisely so that we could minimize mismatch in the LD patterns within the GWAS study cohorts and between the GWAS cohorts and the individuals used to generate the LD matrices. Both types of mismatch are known to produce spurious results and inflated error rates124–126 that can be eliminated by conducting all the GWAS and generating the LD matrices in the exact same sample of individuals, a process known as “in-sample” fine-mapping. Our approach, limited by the summary-level data available in the public domain can be thought of as an approximation that leverages a strongly overlapping set of cohorts and reduces LD mismatch to the greatest extent possible. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint We selected CAFEH for fine-mapping because among a limited number of multi-trait fine- mapping tools51,127–129, CAFEH appeared best suited to the specific task of efficiently mapping a variable selection of traits across a large number of BMD loci. We began applying CAFEH to each locus with the maximum number of signals set to the default, 10. If after initial mapping, at least one signal was detected with a purity value (the absolute correlation between the variants in the signal’s credible set) below 1%, we stopped mapping the locus. If, however, 10 signals were detected with purity values > 1%, we iteratively increased the number of signals by 1 and re-mapped the locus until at least one low-purity signal (< 1%) was detected or a maximum of 30 signals was reached. Seeking to minimize suspicious results and prioritizing precision of reported signals over recall, we post-hoc filtered the CAFEH results by several metrics and reported only signals linked to BMD. First, we reported signals with purity > 50% as signals below this threshold represent instances where the model failed to distinguish between variants in moderately low LD. Second, we linked signals to traits only when the signals’ credible sets had a CAFEH activity score > 0.95 and one or more variants with P < 5 x 10-8 in the trait GWAS. These filters respectively capture signal-trait linkages with strong posterior evidence of trait causality under CAFEH’s Bayesian model and robust frequentist evidence of trait association. Lastly, we reported only signal-trait linkages where a variant in the credible set captures the maximum residual association for the corresponding signal in the given trait. The residual association of each signal represents the remaining GWAS association at each SNP position after removing the effects of the other mapped signals. Mathematically, the residual association is the significance of the residual effect, written as 𝑟−𝑡𝑘 for trait, t, and signal, k. The residual effect can be calculated by subtracting trait t’s first moments of the CAFEH joint model for all signals except k, from the GWAS βs (see equation 63 in CAFEH supplemental methods51) and then dividing the difference by the GWAS standard error. The residual effect approximately follows the unit normal distribution. Under a well-fit CAFEH model, plots of residual association for a given signal should appear similar to LocusZoom GWAS plots130 which either show a null distribution for traits not relevant to the signal or a uni-modal distribution with the credible set at or near the peak for traits linked confidently to the signal. However, in reviewing plots of these residual associations, we observed that often when variants in moderate LD appeared visually in the raw GWAS associations to be causal for .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint different traits, CAFEH would group the variants into a single signal with the credible set driven by the trait with the stronger GWAS association. While these trait linkages could reflect true biology, we doubted that the model was correctly specified for the corresponding signals, so we removed them. Some of these trait linkages were removed by the activity score filter, and we removed the rest by ensuring that the credible set of each signal we reported lay right at the peak of the residual association for all the traits to which it was linked. For reference, we have provided a table of the number of BMD-linked signals detected under different filtering criteria (Supplementary Table 18) and a list of the 1,349 BMD-linked signals identified under the least stringent filtering criteria we considered (only a BMD activity score > 0.5) complete with the information to refilter the signals to any degree of stringency (Supplementary Table 19). After signal filtering, we hierarchically clustered the retained signals by their binarized trait linkages using the complete-linkage method with Euclidean distances implemented in pheatmap131 v1.0.12 (https://CRAN.R-project.org/package=pheatmap). For each signal, we also discretized the CAFEH weight means to -1, 0, and 1. These means reflect the effect-direction relationships between the traits linked to each signal. We standardized the discretized means so that the BMD weight mean would always have a value of 1, and used them to again cluster the signals hierarchically. Additionally, we visualized the binarized trait linkages using the umap package v0.2.10.0 (https://cran.r-project.org/web/packages/umap) and manually annotated the observed clusters by the dominant traits mapped to the signals found in each. DATA AVAILABILITY Upon publication of this manuscript all newly generated data sets will be available on the Gene Expression Omnibus (GEO) at accession number GSE261284. Previously generated raw ATAC- seq and ChIP-seq read files for monocytes43, osteoclasts101, and hFOBs100 are available on GEO at accessions GSE174658, GSE203587, and GSE152942 respectively. Raw capture-C, ATAC-seq, and RNA-seq reads from hMSC-Osteoblasts45 are available on ArrayExpress with the following accession numbers: E-MTAB-6862, E-MTAB-6834, and E-MTAB-6835. GWAS summary statistics obtained from the Genetic Factors for Osteoporosis Consortium (GEFOS) and the Pan-UKBB team are available at the links provided in Supplementary Table 9. CODE AVAILABILITY Public software packages are available at the citations and URLs listed. Custom code for this analysis has been deposited on GitHub (https://github.com/mconery/Grant_hFOB_CRISPRi). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

1. Caliri, A., Filippis, L., Bagnato, G. & Bagnato, G. Osteoporotic fractures: Mortality and quality of life. Panminerva Med. 49, 21–7 (2007). 2. Rizkallah, M. et al. Comparison of morbidity and mortality of hip and vertebral fragility fractures: Which one has the highest burden? Osteoporos. Sarcopenia 6, 146–150 (2020). 3. Sattui, S. E. & Saag, K. G. Fracture mortality: associations with epidemiology and osteoporosis treatment. Nat. Rev. Endocrinol. 10, 592–602 (2014). 4. Krall, E. A. & Dawson-Hughes, B. Heritable and life-style determinants of bone mineral density. J. Bone Miner. Res. 8, 1–9 (1993). 5. Richards, J. B., Zheng, H.-F. & Spector, T. D. Genetics of osteoporosis from genome-wide association studies: advances and challenges. Nat. Rev. Genet. 13, 576–588 (2012). 6. Ng, M. Y . M., Sham, P . C., Paterson, A. D., Chan, V. & Kung, A. W. C. Effect of Environmental Factors and Gender on the Heritability of Bone Mineral Density and Bone Size. Ann. Hum. Genet. 70, 428– 438 (2006). 7. Kim, S. K. Identification of 613 new loci associated with heel bone mineral density and a polygenic risk score for bone mineral density, osteoporosis and fracture. PLOS ONE 13, e0200785 (2018). 8. Morris, J. A. et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat. Genet. 51, 258–266 (2019). 9. Johnell, O. et al. Predictive Value of BMD for Hip and Other Fractures. J. Bone Miner. Res. 20, 1185– 1194 (2005). 10. Tak, Y . G. & Farnham, P . J. Making sense of GWAS: using epigenomics and genome engineering to understand the functional relevance of SNPs in non-coding regions of the human genome. Epigenetics Chromatin 8, 57 (2015). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 11. Sun, Q. et al. From GWAS variant to function: A study of ∼148,000 variants for blood cell traits. Hum. Genet. Genomics Adv. 3, 100063 (2022). 12. Hindorff, L. A. et al. Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proc. Natl. Acad. Sci. 106, 9362–9367 (2009). 13. Zhang, F. & Lupski, J. R. Non-coding genetic variants in human disease. Hum. Mol. Genet. 24, R102– R110 (2015). 14. Trynka, G. et al. Chromatin marks identify critical cell types for fine mapping complex trait variants. Nat. Genet. 45, 124–130 (2013). 15. Degner, J. F. et al. DNase I sensitivity QTLs are a major determinant of human expression variation. Nature 482, 390–394 (2012). 16. Xie, S., Duan, J., Li, B., Zhou, P . & Hon, G. C. Multiplexed Engineering and Analysis of Combinatorial Enhancer Activity in Single Cells. Mol. Cell 66, 285-299.e5 (2017). 17. Xie, S., Armendariz, D., Zhou, P ., Duan, J. & Hon, G. C. Global Analysis of Enhancer Targets Reveals Convergent Enhancer-Driven Regulatory Modules. Cell Rep. 29, 2570-2578.e5 (2019). 18. Gasperini, M. et al. A Genome-wide Framework for Mapping Gene Regulation via Cellular Genetic Screens. Cell 176, 377-390.e19 (2019). 19. Alda-Catalinas, C. et al. Mapping the functional impact of non-coding regulatory elements in primary T cells through single-cell CRISPR screens. Genome Biol. 25, 42 (2024). 20. Fulco, C. P . et al. Activity-by-contact model of enhancer–promoter regulation from thousands of CRISPR perturbations. Nat. Genet. 51, 1664–1669 (2019). 21. Morris, J. A. et al. Discovery of target genes and pathways at GWAS loci by pooled single-cell CRISPR screens. Science 380, eadh7699 (2023). 22. Papalexi, E. et al. Characterizing the molecular regulation of inhibitory immune checkpoints with multimodal single-cell screens. Nat. Genet. 53, 322–331 (2021). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 23. Schraivogel, D. et al. Targeted Perturb-seq enables genome-scale genetic screens in single cells. Nat. Methods 17, 629–635 (2020). 24. Shukla, A. & Huangfu, D. Decoding the noncoding genome via large-scale CRISPR screens. Cell Reprogramming Regen. Repair 52, 70–76 (2018). 25. Cooper, Y . A., Guo, Q. & Geschwind, D. H. Multiplexed functional genomic assays to decipher the noncoding genome. Hum. Mol. Genet. 31, R84–R96 (2022). 26. Wünnemann, F. et al. Multimodal CRISPR perturbations of GWAS loci associated with coronary artery disease in vascular endothelial cells. PLOS Genet. 19, e1010680 (2023). 27. Yihan Wang et al. Enhancer regulatory networks globally connect non-coding breast cancer loci to cancer genes. bioRxiv 2023.11.20.567880 (2023) doi:10.1101/2023.11.20.567880. 28. Armendariz, D. A. et al. CHD-associated enhancers shape human cardiomyocyte lineage commitment. eLife 12, e86206 (2023). 29. Wang, Z. et al. Landscape of enhancer disruption and functional screen in melanoma cells. Genome Biol. 24, 248 (2023). 30. Yang, X. et al. Functional characterization of Alzheimer’s disease genetic variants in microglia. Nat. Genet. 55, 1735–1744 (2023). 31. Chesi, A. et al. Genome-scale Capture C promoter interactions implicate effector genes at GWAS loci for bone mineral density. Nat. Commun. 10, 1260 (2019). 32. Calabrese, G. M. et al. Integrating GWAS and Co-expression Network Data Identifies Bone Mineral Density Genes SPTBN1 and MARK3 and an Osteoblast Functional Module. Cell Syst. 4, 46-59.e4 (2017). 33. Guo, Y . et al. Integrating Epigenomic Elements and GWASs Identifies BDNF Gene Affecting Bone Mineral Density and Osteoporotic Fracture Risk. Sci. Rep. 6, 30558 (2016). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 34. Pippin, J. A. et al. CRISPR-Cas9–Mediated Genome Editing Confirms EPDR1 as an Effector Gene at the BMD GWAS-Implicated ‘STARD3NL’ Locus. JBMR Plus 5, e10531 (2021). 35. Dillard, L. J. et al. Single-Cell Transcriptomics of Bone Marrow Stromal Cells in Diversity Outbred Mice: A Model for Population-Level scRNA-Seq Studies. J. Bone Miner. Res. 38, 1350–1363 (2023). 36. Finucane, H. K. et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat. Genet. 47, 1228–1235 (2015). 37. Cousminer, D. L. et al. Genome-wide association study implicates novel loci and reveals candidate effector genes for longitudinal pediatric bone accrual. Genome Biol. 22, 1 (2021). 38. Medina-Gomez, C. et al. Bone mineral density loci specific to the skull portray potential pleiotropic effects on craniosynostosis. Commun. Biol. 6, 691 (2023). 39. Chen, D. et al. Osteogenic Differentiation Potential of Mesenchymal Stem Cells Using Single Cell Multiomic Analysis. Genes 14, (2023). 40. Abood, A. et al. Identification of Known and Novel Long Noncoding RNAs Potentially Responsible for the Effects of Bone Mineral Density (BMD) Genomewide Association Study (GWAS) Loci. J. Bone Miner. Res. 37, 1500–1510 (2022). 41. He, P . et al. Why SNP rs3755955 is associated with human bone mineral density? A molecular and cellular study in bone cells. Mol. Cell. Biochem. 477, 455–468 (2022). 42. Xia, Q. et al. The type 2 diabetes presumed causal variant within TCF7L2 resides in an element that controls the expression of ACSL5. Diabetologia 59, 2360–2368 (2016). 43. Pahl, M. C. et al. Implicating effector genes at COVID-19 GWAS loci using promoter-focused Capture-C in disease-relevant immune cell types. Genome Biol. 23, 125 (2022). 44. Palermo, J. et al. Variant-to-gene mapping followed by cross-species genetic screening identifies GPI-anchor biosynthesis as a regulator of sleep. Sci. Adv. 9, eabq0844 (2023). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 45. Su, C. et al. Mapping effector genes at lupus GWAS loci using promoter Capture-C in follicular helper T cells. Nat. Commun. 11, 3294 (2020). 46. Sojan, J. M. et al. Bacillus subtilis Modulated the Expression of Osteogenic Markers in a Human Osteoblast Cell Line. Cells 12, (2023). 47. Barry, T., Wang, X., Morris, J. A., Roeder, K. & Katsevich, E. SCEPTRE improves calibration and sensitivity in single-cell CRISPR screen analysis. Genome Biol. 22, 344 (2021). 48. Timothy Barry, Kaishu Mason, Kathryn Roeder, & Eugene Katsevich. Robust differential expression testing for single-cell CRISPR screens at low multiplicity of infection. bioRxiv 2023.05.15.540875 (2023) doi:10.1101/2023.05.15.540875. 49. The Pan UKBB Team. Pan UKBB. https://pan.ukbb.broadinstitute.org. 50. Qu, Y . et al. Genetic Correlation, Shared Loci, and Causal Association Between Sex Hormone- Binding Globulin and Bone Mineral Density: Insights From a Large-Scale Genomewide Cross-Trait Analysis. J. Bone Miner. Res. 38, 1635–1644 (2023). 51. Arvanitis, M., Tayeb, K., Strober, B. J. & Battle, A. Redefining tissue specificity of genetic regulation of gene expression in the presence of allelic heterogeneity. Am. J. Hum. Genet. 109, 223–239 (2022). 52. McInnes, L., Healy, J. & Melville, J. Umap: Uniform manifold approximation and projection for dimension reduction. ArXiv Prepr. ArXiv180203426 (2018). 53. Amariuta, T., Siewert-Rocks, K. & Price, A. L. Modeling tissue co-regulation estimates tissue-specific contributions to disease. Nat. Genet. 55, 1503–1511 (2023). 54. Ongen, H. et al. Estimating the causal tissues for complex traits and diseases. Nat. Genet. 49, 1676– 1683 (2017). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 55. Benjamin J. Strober, Martin Jinye Zhang, Tiffany Amariuta, Jordan Rossen, & Alkes L. Price. Fine- mapping causal tissues and genes at disease-associated loci. medRxiv 2023.11.01.23297909 (2023) doi:10.1101/2023.11.01.23297909. 56. Mostafavi, H., Spence, J. P ., Naqvi, S. & Pritchard, J. K. Systematic differences in discovery of genetic effects on gene expression and complex traits. Nat. Genet. 55, 1866–1875 (2023). 57. Grundberg, E. et al. Population genomics in a disease targeted primary cell model. Genome Res. 19, 1942–1952 (2009). 58. Mullin, B. H. et al. Expression Quantitative Trait Locus Study of Bone Mineral Density GWAS Variants in Human Osteoclasts. J. Bone Miner. Res. 33, 1044–1051 (2018). 59. Magnusson, P ., Degerblad, M., Sääf, M., Larsson, L. & Thorén, M. Different Responses of Bone Alkaline Phosphatase Isoforms During Recombinant Insulin-like Growth Factor-I (IGF-I) and During Growth Hormone Therapy in Adults with Growth Hormone Deficiency. J. Bone Miner. Res. 12, 210– 220 (1997). 60. Bevier, W. C. et al. Relationship of body composition, muscle strength, and aerobic capacity to bone mineral density in older men and women. J. Bone Miner. Res. 4, 421–432 (1989). 61. Sutter, T. et al. Relationships between muscle mass, strength and regional bone mineral density in young men. PLOS ONE 14, e0213681 (2019). 62. Snow-Harter, C., Whalen, R., Myburgh, K., Arnaud, S. & Marcus, R. Bone mineral density, muscle strength, and recreational exercise in men. J. Bone Miner. Res. 7, 1291–1296 (1992). 63. Henderson, N. K., Price, R. I., Cole, J. H., Gutteridge, D. H. & Bhagat, C. I. Bone density in young women is associated with body weight and muscle strength but not dietary intakes. J. Bone Miner. Res. 10, 384–393 (1995). 64. Ho-Pham, L. T., Nguyen, U. D. T. & Nguyen, T. V. Association Between Lean Mass, Fat Mass, and Bone Mineral Density: A Meta-analysis. J. Clin. Endocrinol. Metab. 99, 30–38 (2014). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 65. Katzmarzyk, P . T. et al. Relationship between abdominal fat and bone mineral density in white and African American adults. Interact. Bone Adipose Tissue Metab. 50, 576–579 (2012). 66. Kim, W. et al. The relationship between body fat and bone mineral density in Korean men and women. J. Bone Miner. Metab. 32, 709–717 (2014). 67. Lee, S. J., Lee, J.-Y . & Sung, J. Obesity and Bone Health Revisited: A Mendelian Randomization Study for Koreans. J. Bone Miner. Res. 34, 1058–1067 (2019). 68. Song, J. et al. Causal associations of hand grip strength with bone mineral density and fracture risk: A mendelian randomization study. Front. Endocrinol. 13, (2022). 69. Liu, C. et al. Osteoporosis and sarcopenia-related traits: A bi-directional Mendelian randomization study. Front. Endocrinol. 13, (2022). 70. Felson, D. T., Zhang, Y ., Hannan, M. T. & Anderson, J. J. Effects of weight and body mass index on bone mineral density in men and women: The framingham study. J. Bone Miner. Res. 8, 567–573 (1993). 71. Ma, B. et al. Causal Associations of Anthropometric Measurements With Fracture Risk and Bone Mineral Density: A Mendelian Randomization Study. J. Bone Miner. Res. 36, 1281–1287 (2021). 72. Zhu, K. et al. Relationship between visceral adipose tissue and bone mineral density in Australian baby boomers. Osteoporos. Int. 31, 2439–2448 (2020). 73. Bland, V. L. et al. Metabolically favorable adiposity and bone mineral density: a Mendelian randomization analysis. Obesity 31, 267–278 (2023). 74. Hu, J. et al. Associations of visceral adipose tissue with bone mineral density and fracture: observational and Mendelian randomization studies. Nutr. Metab. 19, 45 (2022). 75. Liu, P .-Y ., Ilich, J. Z., Brummel-Smith, K. & Ghosh, S. New Insight into Fat, Muscle and Bone Relationship in Women: Determining the Threshold at Which Body Fat Assumes Negative Relationship with Bone Mineral Density. Int. J. Prev. Med. 5, 1452–1463 (2014). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 76. Kemp, J. P ., Sayers, A., Smith, G. D., Tobias, J. H. & Evans, D. M. Using Mendelian randomization to investigate a possible causal relationship between adiposity and increased bone mineral density at different skeletal sites in children. Int. J. Epidemiol. 45, 1560–1572 (2016). 77. Mullin, B. H. et al. Characterisation of genetic regulatory effects for osteoporosis risk variants in human osteoclasts. Genome Biol. 21, 80 (2020). 78. He, D. et al. A longitudinal genome-wide association study of bone mineral density mean and variability in the UK Biobank. Osteoporos. Int. 34, 1907–1916 (2023). 79. Dong, H. et al. Comprehensive Analysis of the Genetic and Epigenetic Mechanisms of Osteoporosis and Bone Mineral Density. Front. Cell Dev. Biol. 8, (2020). 80. Timshel, P . N., Thompson, J. J. & Pers, T. H. Genetic mapping of etiologic brain cell types for obesity. eLife 9, e55851 (2020). 81. Greenbaum, J. et al. A multiethnic whole genome sequencing study to identify novel loci for bone mineral density. Hum. Mol. Genet. 31, 1067–1081 (2022). 82. Takeshita, S., Kaji, K. & Kudo, A. Identification and Characterization of the New Osteoclast Progenitor with Macrophage Phenotypes Being Able to Differentiate into Mature Osteoclasts. J. Bone Miner. Res. 15, 1477–1488 (2000). 83. Fairley, S., Lowy-Gallego, E., Perry, E. & Flicek, P . The International Genome Sample Resource (IGSR) collection of open human genomic variation resources. Nucleic Acids Res. 48, D941–D947 (2020). 84. Cody, J. J. et al. A simplified method for the generation of human osteoclasts in vitro. Int. J. Biochem. Mol. Biol. 2, 183–189 (2011). 85. Susa, M., Luong-Nguyen, N.-H., Cappellen, D., Zamurovic, N. & Gamse, R. Human primary osteoclasts: in vitro generation and applications as pharmacological and clinical assay. J. Transl. Med. 2, 6 (2004). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 86. Shen, J. et al. DNA methyltransferase 3b regulates articular cartilage homeostasis by altering metabolism. JCI Insight 2, e93612. 87. Benjamin C. Hitz et al. The ENCODE Uniform Analysis Pipelines. bioRxiv 2023.04.04.535623 (2023) doi:10.1101/2023.04.04.535623. 88. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012). 89. John M. Gaspar. Improved peak-calling with MACS2. bioRxiv 496521 (2018) doi:10.1101/496521. 90. Qunhua Li, James B. Brown, Haiyan Huang, & Peter J. Bickel. Measuring reproducibility of high- throughput experiments. Ann. Appl. Stat. 5, 1752–1779 (2011). 91. Stelzer, G. et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr. Protoc. Bioinforma. 54, 1.30.1-1.30.33 (2016). 92. Sudlow, C. et al. UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age. PLOS Med. 12, e1001779 (2015). 93. Bulik-Sullivan, B. et al. An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47, 1236–1241 (2015). 94. Bulik-Sullivan, B. K. et al. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat. Genet. 47, 291–295 (2015). 95. Auton, A. et al. A global reference for human genetic variation. Nature 526, 68–74 (2015). 96. Altshuler, D. M. et al. Integrating common and rare genetic variation in diverse human populations. Nature 467, 52–58 (2010). 97. Hinrichs, A. S. et al. The UCSC Genome Browser Database: update 2006. Nucleic Acids Res. 34, D590–D598 (2006). 98. Bernstein, B. E. et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat. Biotechnol. 28, 1045–1048 (2010). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 99. Kundaje, A. et al. Integrative analysis of 111 reference human epigenomes. Nature 518, 317–30 (2015). 100. Cottone, L. et al. Aberrant paracrine signalling for bone remodelling underlies the mutant histone- driven giant cell tumour of bone. Cell Death Differ. 29, 2459–2471 (2022). 101. Bae, S. et al. RANKL-responsive epigenetic mechanism reprograms macrophages into bone- resorbing osteoclasts. Cell. Mol. Immunol. 20, 94–109 (2023). 102. Bonferroni, C. Teoria statistica delle classi e calcolo delle probabilita. Pubblicazioni R Ist. Super. Sci. Econ. E Commericiali Firenze 8, 3–62 (1936). 103. Benjamini, Y . & Hochberg, Y . Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B Methodol. 57, 289–300 (1995). 104. Su, C. et al. 3D promoter architecture re-organization during iPSC-derived neuronal cell differentiation implicates target genes for neurodevelopmental disorders. Prog. Neurobiol. 201, 102000 (2021). 105. Wingett, S. W. et al. HiCUP: pipeline for mapping and processing Hi-C data. Preprint at https://doi.org/10.12688/f1000research.7334.1 (2015). 106. Frankish, A. et al. GENCODE 2021. Nucleic Acids Res. 49, D916–D923 (2021). 107. Cairns, J. et al. CHiCAGO: robust detection of DNA looping interactions in Capture Hi-C data. Genome Biol. 17, 127 (2016). 108. Arnold, M., Raffler, J., Pfeufer, A., Suhre, K. & Kastenmüller, G. SNiPA: an interactive, genetic variant- centered annotation browser. Bioinformatics 31, 1334–1336 (2015). 109. Thakore, P . I. et al. Highly specific epigenome editing by CRISPR-Cas9 repressors for silencing of distal regulatory elements. Nat. Methods 12, 1143–1149 (2015). 110. Doench, J. G. et al. Optimized sgRNA design to maximize activity and minimize off-target effects of CRISPR-Cas9. Nat. Biotechnol. 34, 184–191 (2016). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 111. Sanson, K. R. et al. Optimized libraries for CRISPR-Cas9 genetic screens with multiple modalities. Nat. Commun. 9, 5416 (2018). 112. McKenna, A. & Shendure, J. FlashFry: a fast and flexible tool for large-scale CRISPR target design. BMC Biol. 16, 74 (2018). 113. Sigma Aldrich. CRISPRi Human Whole Genome and Long Non-Coding Library Screening User Manual. (2021). 114. Zheng, G. X. Y . et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017). 115. Fleming, S. J. et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat. Methods 20, 1323–1335 (2023). 116. Wolf, F. A., Angerer, P . & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018). 117. Traag, V. A., Waltman, L. & van Eck, N. J. From Louvain to Leiden: guaranteeing well-connected communities. Sci. Rep. 9, 5233 (2019). 118. Hao, Y . et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. (2023) doi:10.1038/s41587-023-01767-y. 119. Lopez-Delisle, L. et al. pyGenomeTracks: reproducible plots for multivariate genomic datasets. Bioinformatics 37, 422–423 (2021). 120. Hammond, R. K. et al. Biological constraints on GWAS SNPs at suggestive significance thresholds reveal additional BMI loci. eLife 10, e62206 (2021). 121. Xie, J. et al. Calcitonin and Bone Physiology: In Vitro, In Vivo, and Clinical Investigations. Int. J. Endocrinol. 2020, 3236828 (2020). 122. Naot, D., Musson, D. S. & Cornish, J. The Activity of Peptides of the Calcitonin Family in Bone. Physiol. Rev. 99, 781–805 (2019). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint 123. Anurag Verma et al. Diversity and Scale: Genetic Architecture of 2,068 Traits in the VA Million Veteran Program. medRxiv 2023.06.28.23291975 (2023) doi:10.1101/2023.06.28.23291975. 124. Chen, W. et al. Improved analyses of GWAS summary statistics by reducing data heterogeneity and errors. Nat. Commun. 12, 7117 (2021). 125. Kanai, M. et al. Meta-analysis fine-mapping is often miscalibrated at single-variant resolution. Cell Genomics 2, (2022). 126. Yang, Z. et al. CARMA is a new Bayesian model for fine-mapping in genome-wide association meta- analyses. Nat. Genet. 55, 1057–1065 (2023). 127. Wallace, C. A more accurate method for colocalisation analysis allowing for multiple causal variants. PLOS Genet. 17, e1009440 (2021). 128. Zhou, F. et al. Leveraging information between multiple population groups and traits improves fine- mapping resolution. Nat. Commun. 14, 7279 (2023). 129. Yuxin Zou, Peter Carbonetto, Dongyue Xie, Gao Wang, & Matthew Stephens. Fast and flexible joint fine-mapping of multiple traits via the Sum of Single Effects model. bioRxiv 2023.04.14.536893 (2024) doi:10.1101/2023.04.14.536893. 130. Pruim, R. J. et al. LocusZoom: regional visualization of genome-wide association s36can results. Bioinformatics 26, 2336–2337 (2010). 131. Kolde, R. & Kolde, M. R. Package ‘pheatmap’. R Package 1, 790 (2015). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint

The authors would like to thank AlloSource and the University of Colorado Interdisciplinary Joint Biology Program Biorepository for providing healthy human articular cartilage tissue samples for the chondrocyte ATAC-seq experiments and also Ms. Samantha Landgrave for isolating chondrocytes from this tissue to support sequencing. Additionally, the authors would also like to acknowledge the Center for Applied Genomics (CAG) at CHOP for their assistance with the 10X scRNA-seq library generation. FUNDING MJZ is supported by the University of Colorado Gates Grubstake Award; EK by the NSF (DMS 2113072 and DMS 2310654); ADW by NIAID (R01AI154773); ADW and SFAG by NIDDK (R01DK122586); BSZ by NCATS (UL1 TR001878); SFAG, BSZ, and KDH by the NICHD (R01 HD100406); SFAG and KDH by the NIA (R01 AG072705); SFAG and BV by the NIDDK (UM1 DK126194); KDH by the Henry Ruppenthal Family Professorship for Bioengineering and Orthopaedic Surgery; and SFAG by the Daniel B. Burke Endowed Chair for Diabetes Research. AUTHOR CONTRIBUTIONS Conceptualization: MC, JAP, KT, ADW, BFV, AC, SFAG; Methodology: MC, JAP, YW, KT, MCP , EK, BFV, AC, SFAG; Investigation: MC, JAP, YW, KT, DAV, LJF; Visualization: MC, JAP , YW; Funding acquisition: MJZ, EK, ADW, BSZ, KDH, SFAG; Supervision: BFV, KDH, AC, SFAG; Writing – original draft: MC, JAP, YW, DAV; Writing – review & editing: MC, JAP , YW, KT, MCP, DAV, CLA, MJZ, EK, ADW, BSZ, BFV, KDH, AC, SFAG COMPETING INTERESTS The authors declare no competing interests. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint FIGURES Fig. 1 Bone Mineral Density Partitioned Heritability Enrichments across 98 Cell Types. Each bar in the figure represents a particular genomic annotation (H3K27ac, H3K9ac, H3K4me1, H3K4me3, or open-chromatin) measured in a specific primary cell type or cell model. Negative log10 Bonferroni-adjusted p-values are plotted along the y-axis. The dashed line reflects a Bonferroni-adjusted significance cutoff of 0.05. Coloring reflects manually curated tissue categories for each cell type. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Fig. 2 hFOB CRISPRi Screen Targets and Perturbation Results. (a) Breakdown of 89 screen targets by origin. (b) Quantile-quantile plot of targeting and non-targeting (negative-control) sgRNA tests. (c) Volcano plot of targeting sgRNA screen results. The 23 genes that exhibited significant perturbation at a target site are labeled. pyGenomeTracks plots of the (d) ARID5B and (e) EYA2 loci showing ATAC-seq read and Capture-C chromatin loops measured in hMSC- Osteoblasts, hFOBs, and hMSC-adipocytes. Targeted SNPs and genes are plotted below. All isoforms of genes perturbed in the CRISPRi screen are colored red. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Fig. 3 Assays of siRNA Knockdown on Osteoblast and Adipocyte Maturation and Function. Alkaline phosphatase (ALP) assay in hFOBs (top panel), ALP assay in hMSC- Osteoblasts (second from top), alizarin Red S (ARS) assay in hMSC-Osteoblasts (third from top), and adipogenesis assay in hMSCs (bottom). siRNA targets are listed along the x-axis and normalized assay measurements, along the y-axis. All measurements were normalized relative to their plate’s respective control-siRNA well such that the controls have a value of one (dashed red line). hMSC-based assays were measured and normalized in the treated state (with BMP2 or adipogenic induction media). Any siRNAs resulting in significant decreases of assayed measurements are displayed in blue. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Fig. 4 Genetic Correlations between BMD and 38 Traits Estimated in White British Individuals. Genetic correlations are displayed the vertical axis and traits along the horizontal. The 38 traits include the correlation of BMD with itself. Significant correlations (Benjamini- Hochberg Adj. P-Value < 0.05) are labeled with an asterisk. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint Fig. 5. Cross-Trait Signal Sharing and Effect Directions of 433 BMD Signals across 22 Traits. (a) Heatmap and dendrogram of trait mappings. Traits are plotted down the y-axis with individual signals plotted along the x-axis. Black coloring indicates when a signal modulates the corresponding trait. Signals and traits are clustered hierarchically using the complete-linkage

with Euclidean distances. (b) Breakdown of signals mapped per trait split into signals with positively correlated effects on BMD and given trait (red) and signals with negatively correlated effects (blue). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 20, 2024. ; https://doi.org/10.1101/2024.03.19.585778doi: bioRxiv preprint