Bivariate genome-wide association meta-analysis of pediatric musculoskeletal traits reveals pleiotropic effects at the SREBF1/TOM1L2 locus

Bone mineral density is known to be a heritable, polygenic trait whereas genetic variants contributing to lean mass variation remain largely unknown. We estimated the shared SNP heritability and performed a bivariate GWAS meta-analysis of total-body lean mass (TB-LM) and total-body less head bone mineral density (TBLH-BMD) regions in 10,414 children. The estimated SNP heritability is 43% (95% CI: 34–52%) for TBLH-BMD, and 39% (95% CI: 30–48%) for TB-LM, with a shared genetic component of 43% (95% CI: 29–56%). We identify variants with pleiotropic effects in eight loci, including seven established bone mineral density loci: WNT4, GALNT3, MEPE, CPED1/WNT16, TNFSF11, RIN3, and PPP6R3/LRP5. Variants in the TOM1L2/SREBF1 locus exert opposing effects TB-LM and TBLH-BMD, and have a stronger association with the former trait. We show that SREBF1 is expressed in murine and human osteoblasts, as well as in human muscle tissue. This is the first bivariate GWAS meta-analysis to demonstrate genetic factors with pleiotropic effects on bone mineral density and lean mass.


Supplementary
. Prioritization of genes at the 17p11.2 using SMR analysis. Top plot, Brown dots represent the Pvalues for SNPs in our lean mass GWAS meta-analysis, diamonds represent the P values for probes from the SRM test and triangles represent probes without a cis-eQTL at P eQTL < 5x10 -8 . Bottom plot, the eQTL P-values of SNPs from the Westra et al. 1 study for the ILMN_1663035 probe tagging SREBF1. The top and bottom plots include all the SNPs available in the region in the GWAS and eQTL summary data, respectively, rather than only the SNPs common to both data sets. Highlighted in red is the probe tagging SREBF1 that passed the SMR and HEIDI tests.  Table 1. Genes Located within the 17p11.2 locus. We here summarize the information available on the genes located in 17p11.2 related to musculoskeletal traits. OMIM annotations are only considered if they present with a musculoskeletal phenotype. Expression in murine osteoblast are based in clavaria data (see Methods). Expression in skeletal muscle is based in GTEx data [2][3][4] Gene Symbol Gene name and Known function OMIM annotation Mice skeletal phenotype Expression in murine osteoblast Expression in Skeletal muscle

RAI1
Retinoic acid-inducible 1. Although the function of this protein is unknown, it is thought to be involved in nervous system development.

Description of TBLM/TBLH-BMD Associated Loci
Seven of the eight loci associated with TBLH-BMD and TBLM in this pediatric bivariate meta-analysis were previously described in GWAS of BMD at different skeletal sites. In the 1p36.12 locus, the bivariate signal arises from a region in linkage disequilibrium (LD) extending 264 Kb, with GWS variants located within WNT4 and in close vicinity to the 5'region of ZBTB40 (lead SNP rs6684375, P=2.8x10 -9 ). The low LD between the SNPs located at the opposite ends of the described region (r 2 =0.005) are indicative of at least two independent association signals mapping to this locus 14

Genetic Relatedness Matrix Generation for the Generation R Study
As data from Generation R Study comes from and admixed population and the estimation of the relationship coefficient from GCTA under this scenario suffer a bias. One important assumption of this method is that the sampling population is an homogeneous population, the  19 . This issue becomes of importance when high sensitivity is required (i.e., for the assessment of distant relatedness and/or fine pedigree structure), as it is the case of heritability estimates.
Hence, we have modified the genetic relationship matrix (GRM) by using the REAP estimator as follows.
The GRM coefficient from the GCTA methodology as described in 20 is Where A jk is the relatedness coefficient for individuals j and k. N the number of SNPs under study. x i the genotype of individual j or k at marker I (e.g., 0,1,2) and p i the population minor allele frequency of marker i.
REAP condition the population frequency for each individual p m a linear combination of the subpopulation allele frequencies based on individual's ancestry 19 and thus, As originally REAP by default produce kinship coefficients φ jk the coefficients reported should be multiplied by 2 as shown in 19 .