Systems genetics in diversity outbred mice inform BMD GWAS and identify determinants of bone strength

Abstract

Genome-wide association studies (GWASs) for osteoporotic traits have identified over 1000 associations; however, their impact has been limited by the difficulties of causal gene identification and a strict focus on bone mineral density (BMD). Here, we use Diversity Outbred (DO) mice to directly address these limitations by performing a systems genetics analysis of 55 complex skeletal phenotypes. We apply a network approach to cortical bone RNA-seq data to discover 66 genes likely to be causal for human BMD GWAS associations, including the genes SERTAD4 and GLT8D2. We also perform GWAS in the DO for a wide-range of bone traits and identify Qsox1 as a gene influencing cortical bone accrual and bone strength. In this work, we advance our understanding of the genetics of osteoporosis and highlight the ability of the mouse to inform human genetics.

Fig. 1. Resource overview.

An overview of the resource including data generated and analyses performed. GWAS genome-wide association study, QTL quantitative trait loci, eQTL expression quantitative trait loci, SNP single-nucleotide polymorphism.

Fig. 2. Characterization of the experimental Diversity Outbred cohort.

a Allele frequency per chromosome, across the DO cohort. Intervals represent the eight DO founder strains: A/J (yellow), C57BL/6J (gray), 129S1/SvImJ (beige), NOD/ShiLtJ (dark blue), NZO/HILtJ (light blue), CAST/EiJ (green), PWK/PhJ (red), and WSB/EiJ (purple). b Bone volume fraction and max load across the DO cohort. Insets are microCT images representing low and high bone volume fraction (BV/TV). c Heritability of each bone trait. Phenotypes are colored by phenotypic category: morphology (purple), marrow adiposity (light blue), histomorphometry (green), microarchitecture (olive), and biomechanics (beige). Abbreviations for phenotypes are available in Supplementary Data 1.

Fig. 3. Overview of the network approach used to identify genes potentially responsible for BMD GWAS loci.

Three WGCNA networks (124 total modules) were constructed from RNA-seq data on cortical bone in the DO (N = 192). A Bayesian network was then learned for each module. We performed key driver analysis on each Bayesian network to identify BANs, by identifying nodes (genes) that were more connected to more known bone genes than was expected by chance. We colocalized GTEx human eQTL for each BAN with GWAS BMD SNPs to identify potentially causal genes at BMD GWAS loci. For the key driver analysis, the yellow node indicates the queried gene, red nodes indicate known bone genes, and gray nodes indicate non-bone genes. DO diversity outbred, WGCNA weighted gene co-expression network analysis, BAN bone-associated nodes, eQTL expression quantitative trait loci, BMD bone mineral density, GWAS genome-wide association studies, SNP single-nucleotide polymorphism, PPH4 posterior probability of colocalization, hypothesis 4.

Fig. 4. Identifying SERTAD4 and GLT8D2 as putative regulators of BMD.

a Local 3-step neighborhood around Sertad4. Known bone genes highlighted in green. Sertad4 highlighted in red. b Local 3-step neighborhood around Glt8d2. Known bone genes highlighted in green. Glt8d2 highlighted in red. c Expression of Sertad4 and Glt8d2 in calvarial osteoblasts. For each time point, N = 3 independent biological replicates were examined. Error bars represent the standard error of the mean. TPM transcripts per million. d Single-cell RNA-seq expression data. Each point represents a cell (N = 7092 cells). The top panel shows UMAP clusters and their corresponding cell-type. The bottom two panels show the expression of Sertad4 and Glt8d2. The color scale indicates normalized gene expression value. e Bone mineral density in Glt8d2 knockout mice from the IMPC. N = 7 females and N = 7 males for Glt8d2^−/− mice, N = 1466 females and N = 1477 males for Glt8d2^+/+ mice. Boxplots indicate the median (middle line), the 25th and 75th percentiles (box) and the whiskers extend to 1.5 * IQR. Colors indicate genotype.

Fig. 5. Overview of our approach to QTL fine-mapping.

a Overview of merge analysis. LOD logarithm of the odds, QTL quantitative trait loci, DO diversity outbred, SNP single-nucleotide polymorphism, INDEL insertion–deletion, SV structural variant. b Overview of merge analysis as performed for the identification of missense variants. c Overview of merge analysis as performed for the identification of colocalizing trait QTL/gene eQTL within a locus. The pink columns around the QTL in each association plot represent the QTL 95% confidence intervals. The yellow box in c represents the gene search space for a locus, defined as the region within ±250 Kbp around the outer boundaries of the 95% confidence intervals within a locus. eQTL expression quantitative trait loci.

Fig. 6. QTL (locus 1) on chromosome 1.

a For each plot, the top panel shows allele effects for the DO founders for each of the eight QTL (quantitative trait loci) across an interval on chromosome 1 (Mbp, colors correspond to the founder allele in the legend). Bottom panels show each respective QTL scan. The red horizontal lines represent LOD (logarithm of the odds) score thresholds (genome-wide P < = 0.05). Colors indicate founder mouse strains: A/J (yellow), C57BL/6J (gray), 129S1/SvImJ (beige), NOD/ShiLtJ (dark blue), NZO/HILtJ (light blue), CAST/EiJ (green), PWK/PhJ (red), and WSB/EiJ (purple). b QTL scans across the same interval as a, after conditioning on rs50769082. c QTL scans after conditioning on rs248974780. Phenotype abbreviations: TMD tissue mineral density, ML medial–lateral femoral width, pMOI polar moment of intertia, Imax maximum moment of inertia, Ct.Ar/Tt.Ar bone area fraction, Tt.Ar total area, Ma.Ar medullary area, Ct.Por cortical porosity.

Fig. 7. Characterization of Qsox1.

a The top panel shows allele effects for the DO founders for Ier5 and Qsox1 expression an interval on chromosome 1 (Mbp, colors correspond to the founder allele in the legend). Y-axis units are best linear unbiased predictors (BLUPs). Bottom panels show each respective QTL scan. LOD (logarithm of the odds) score threshold for autosomal eQTL is 10.89 (alpha = 0.05). b Qsox1 expression in calvarial osteoblasts. For each time point, N = 3 independent biological replicates were examined. Error bars represent the standard error of the mean. TPM transcripts per million. c Single-cell RNA-seq expression data. Each point represents a cell (N = 7092 cells). The top panel shows UMAP clusters and their corresponding cell-type. The bottom panel shows the expression of Qsox1. The color scale indicates normalized gene expression value.

Fig. 8. Qsox1 is responsible for several chromosome 1 QTL.

a Representative image of the Qsox1 knockout mutations. b QSOX1 activity assay in serum. Data is grouped by mouse genotype. Boxplots indicate the median (middle line), the 25th and 75th percentiles (box) and the whiskers extend to 1.5 * IQR. Colors indicate mutation type. c–e Femoral morphology in Qsox1 mutant mice. P-values above plots are ANOVA P-values for the genotype term, while P-values in the plots are contrast P-values, adjusted for multiple comparisons. The center points of the plots represent the least-squares mean, while the error bars represent the confidence intervals at a confidence level of 0.95. f–l microCT measurements of chromosome 1 QTL phenotypes in Qsox1 knockout mice. P-values in the plots are contrast P-values, adjusted for multiple comparisons. The center points of the plots represent the least-squares mean, while the error bars represent the confidence intervals at a confidence level of 0.95. m Bone strength (max load, F_max) in the AP orientation, measured via four-point bending. P-values in the plots are contrast P-values, adjusted for multiple comparisons. The center points of the plot represent the least-squares mean, while the error bars represent the confidence intervals at a confidence level of 0.95. n Bone strength (max load, F_max) in the ML orientation, measured via four-point bending. P-values in the plots are contrast P-values, adjusted for multiple comparisons. The center points of the plot represent the least-squares mean, while the error bars represent the confidence intervals at a confidence level of 0.95. o Representative microCT images of the effect of Qsox1 on bone size. All error bars represent confidence intervals at a 95% confidence level. ML medial-lateral femoral width, AP anterior-posterior femoral width, FL femoral length, pMOI polar moment of inertia, Imax maximum moment of inertia, Ct.Ar/Tt.Ar bone area fraction, Tt.Ar total area, Ma.Ar medullary area, TMD tissue mineral density, Ct.Por cortical porosity.