Introduction

In 2019, 32 million Europeans were estimated to have osteoporosis with the highest number of affected persons being estimated for Germany [1]. Osteoporosis is marked by reduced bone mineral density (BMD) and an increased risk of fractures. Affected patients suffer from increased morbidity and mortality [2] which is accompanied by reduced quality of life [3].

According to the World Health Organization (WHO), osteoporosis is defined by a dual X-ray absorptiometry (DXA) T-score of −2.5 or lower [4]. German [5] and international guidelines [1], therefore, recommend DXA measurements for diagnosis of osteoporosis. An alternative method to BMD measurement by DXA is the quantitative ultrasound (QUS) technique. Although QUS measurements are not directly comparable to DXA, they are free of ionising radiation, allow assessing physical bone properties, and predict fracture risk at different sites [6].

Osteoporosis risk is determined by both, genetic and non-genetic, factors. Non-genetic factors include lifestyle-related conditions like being underweight and smoking as well as unmodifiable factors including age and intake of certain medications [7]. A parental history of hip fracture is another important risk factor that points to the impact of the individual genetic disposition to bone phenotypes [8]. Indeed, not only BMD (h2 = 50–80%) but also ultrasound bone properties (h2 = 40–50%) possess a high heritability (h2) [8]. This resulted in efforts to find genetic variants contributing to BMD, bone microarchitecture, and fractures [8]. Causative mutations for rare but severe bone diseases such as osteogenesis imperfecta and sclerosteosis were successfully identified [9]. These important discoveries not only broaden our understanding of the physiology of bone but may also contribute to the identification of new pharmacotherapeutic targets [10].

Next to monogenic mutations causing rare diseases, it is of high interest to assess common genetic variants that contribute to the disease risk in the general population. Genome-wide association studies (GWAS) in large samples provide a powerful tool for the identification of disease-related genes in complex and highly prevalent diseases like osteoporosis [8]. So far, several GWAS have been performed to identify genes associated with BMD, QUS measures, bone microarchitecture, and fractures to finally trace genetic variants involved in the development of osteoporosis (for a review see [11]). Meta-analyses of GWAS identified and replicated several loci associated with BMD including single nucleotide polymorphisms (SNPs) related to WNT16, RSPO3, TMEM135, and many more [12]. In the latest and largest GWAS in the field, including more than 426,000 individuals from the UK Biobank, Morris and Kemp et al. [13] identified 1103 conditional independent SNPs associated with QUS-derived estimated BMD (eBMD) at a genome-wide significant level. These were mapped to 518 loci, among which DAAM2 was selected, in-depth characterised, and suggested as a promising target for further investigation [13].

Here, we followed up on the results of Morris and Kemp et al. [13]. We assessed the associations of the 1103 SNPs with the QUS-derived stiffness index in two cohorts of the population-based Study of Health in Pomerania (SHIP). We aimed to examine which associations replicate in our precisely phenotyped population. Moreover, we investigated our results in two cohorts of patients evaluated for osteoporosis with DXA measurements at the spine and femoral neck.

Materials and Methods

Study Populations

SHIP

SHIP was established to collect and analyse data on health and disease in Northeast Germany. It consists of two non-overlapping, population-based cohorts, SHIP-START and SHIP-TREND. Both cohorts are based on representative samples of the adult inhabitants of the study region. Details on study design and sampling can be found elsewhere [14]. The study was approved by the ethics committee of the University of Greifswald and is conducted in line with the Declaration of Helsinki including obtainment of written informed consent from all participants.

In the present study, data from the second follow-up of the SHIP-START cohort (SHIP-START-2, n = 2333) and the baseline examination of the SHIP-TREND cohort (SHIP-TREND-0, n = 4420) were analysed as only in these study waves, quantitative ultrasound measurements (QUS) at the heel were performed. Data collection in SHIP-START-2 and SHIP-TREND-0 was conducted in parallel between 2008 and 2012 with similar methods and protocols [14]. From the total of 6753 SHIP-START-2 and SHIP-TREND-0 participants, we excluded all subjects with missing QUS or genotyping data, and all subjects that were treated with systemic glucocorticoids (ATC classification: H02AB), bisphosphonates (ATC classification: M05BA, M05BB), or other drugs affecting bone structure and mineralisation (ATC classification: M05BX). The final study population comprised 2108 SHIP-START-2 and 3557 SHIP-TREND-0 participants.

All SHIP participants in both cohorts underwent an extensive computer-assisted personal interview on lifestyle, medical history, and socio-demographic characteristics, and a large range of medical tests (for details see [14]). Standardised measurements of body height and weight were performed with calibrated scales, and body mass index (BMI) was calculated as weight (kg)/height2 (m2). All participants were offered whole body magnetic resonance imaging (MRI). From the images, the amount of abdominal subcutaneous (SAT) and visceral (VAT) adipose tissue was quantified [15]. Women aged 60 years or older and women aged between 40 and 60 years without menstrual cycling were classified as postmenopausal, all further women as premenopausal. Regular medication intake was categorised according to the anatomical-therapeutic-chemical (ATC) classification system. Information on secondary causes of osteoporosis was not collected.

OsteoGene

OsteoGene (DRKS ID: DRKS00016601) is a prospective study recruiting patients evaluated for osteoporosis at the community health centre MVZ endokrinologikum Göttingen (Germany). Enrolled patients were aged between 18 and 88 years and had a 20% increased 10-year fracture risk for vertebral or hip fractures. According to German guidelines [5], these patients underwent several diagnostic measures including DXA measurements. 98.3% of the patients were therapy naïve. Intake of inhalative or oral glucocorticoids was defined as exclusion criteria. In contrast, intake of calcium or vitamin D supplements, or hormone replacement therapy was no exclusion criteria. Additionally, information on secondary causes of osteoporosis was collected. The final study population included 232 patients that were recruited between December 2017 and October 2020. The study was approved by the Ethics Review Committee of the University Medical Center Göttingen. All participants provided written informed consent.

HSD

HSD was performed retrospectively in 452 German subjects. In short, patients that were to be evaluated for osteoporosis in the endocrine outpatient clinic of the University Medical Center in Göttingen were enrolled. The study population included in the present analyses comprised 399 patients with complete information on age and information on BMD of at least one location (spine, femoral neck). Additionally, information on secondary causes of osteoporosis was collected. Further details on the HSD cohort have been published previously [16].

Assessment of Bone Properties and Fractures

SHIP

QUS measurements were performed at the heel of both feet using an Achilles InSight System (GE Medical Systems Ultrasound, GE Healthcare, Chalfont St Giles, UK). In short, two ultrasound parameters, the broadband ultrasound attenuation (BUA) and the speed of sound (SOS), were measured. These measures were combined to form the stiffness index (SI) according to the following formula: SI = (0.67 × BUA) + (0.28 × SOS)-420. The stiffness index serves as an indicator of the osteoporotic fracture risk. Statistical analyses were performed with data from the foot with the lower stiffness index. QUS measurements were not performed when the participant had implants, prostheses, or amputations in or below the knee, wounds, or infections distal to the knee, or oedema. Data from participants who reported an injury or surgery below the knee within twelve months prior to the measurement, who used a wheelchair or could not correctly place the feet into the device, were excluded from the statistical analyses. Data on self-reported incident fractures since the baseline examination were collected in SHIP-START-2 and data on selected lifetime fractures (proximal humerus, vertebral, hip, or femoral neck fractures) were collected in SHIP-TREND-0.

OsteoGene

Areal BMD (g/cm2) was measured by dual-energy X-ray absorptiometry (DXA) at the lumbar spine (L1-L4), total femur, and femoral neck of both legs using a LUNAR Prodigy instrument (GE Healthcare, Chicago, IL, USA). T-Score and Z-Score were automatically determined by the instrument. For analyses, no less than two vertebrae and only vertebrae without fractures were included. Vertebral fractures were also assessed by DXA scan. Peripheral fracture rate was assessed by already available X-ray or MRI scans.

HSD

BMD was determined by DXA measurements at the lumbar spine and the left femoral neck. Fractures were self-reported and partially cross-checked against radiology reports and fracture clinic attendance.

Genotyping

SHIP

SHIP-START participants were genotyped applying the Affymetrix Genome-Wide Human SNP Array 6.0 (Santa Clara, CA, USA). SHIP-TREND-0 participants were genotyped applying either the Illumina Infinium® HumanOmni2.5 BeadChip or the Illumina Infinium® Global Screening Array (San Diego, CA, USA). Genotyping was performed according to the manufacturer’s protocol. Whole-genome imputation was performed on the Michigan Imputation Server using the HRC reference panel (version r1.1 2016).

OsteoGene and HSD

DNA was isolated from blood samples with the QIAamp DNA Blood Mini Kit (Qiagen, Hilden, Germany). The rhAmp SNP Genotyping System (Integrated DNA Technologies, Carolville, IA, USA) was used to genotype rs2707518 (CPED1/WNT16; Assay ID: CD.GT.FSGQ5187.1), rs3779381 (WNT16; CD.GT.PBLY8533.1), rs603424 (PKD2L1; Hs.GT.rs603424.A.1), rs10239787 (JAZF1; Hs.GT.rs10239787.T.1), and rs6968704 (JAZF1; Hs.GT.rs6968704.T.1) in 5 µl reactions in 384-well plates according to the manufacturer’s protocol but using undiluted DNA. The PCR, data collection, and analysis were conducted in a QuantStudio 12k Flex Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA).

Statistical analyses

Characteristics of the SHIP participants and the patients evaluated for osteoporosis are reported as means with standard deviation or proportions.

In SHIP, associations between the SNPs and stiffness index were determined separately for the two cohorts using multivariate linear regression models implemented in EPACTS version 3.2.6 patched (http://csg.sph.umich.edu//kang/epacts/download/). Sex and age were defined as covariates. As genotyping in the SHIP-TREND cohort was performed with two arrays, three further covariates were defined for this cohort: genotyping array and the first two genetic principal components. The individual results were combined by fixed effects inverse-variance weighted meta-analysis using METAL [17]. The false discovery rate (FDR) at 5% using the Benjamini–Hochberg procedure was calculated to account for multiple testing [18]. Results were called significant when the FDR was < 0.05. We report effect estimates with standard error, p-value, and FDR from these models. The results of the meta-analyses were further illustrated in a plot depicting the absolute effect size in relation to the minor allele frequency (MAF).

We then selected the five SNPs with the lowest p-values from the meta-analysis for genotyping in the OsteoGene and HSD study cohorts: rs2707518, rs3779381, rs115242848, rs10239787, and rs603424. These SNPs are located in CPED1/WNT16, WNT16, LOC101927709/ EN1, JAZF1, and PKD2L1, respectively. In addition, rs6968704 (JAZF1) which was also significantly associated with stiffness index, was selected for genotyping. Thus, a total of six SNPs were genotyped in the patient cohorts. Associations between the SNPs and BMD at the femoral neck or spine were assessed with linear regression models adjusted for sex and age (IBM SPSS Statistics v.26, IBM, Armonk, NY, USA). Next to a combined model with pooled data from OsteoGene and HSD, also separate models for the two cohorts were calculated. Finally, we compared the MAFs of the six selected SNPs between SHIP participants and patients evaluated for osteoporosis.

Linkage disequilibrium was analysed using SNiPA [19] with the following settings: Genome assembly: GRCh37, Variant set: 1000 Genomes, Phase 3 v5, Population: European. Data plotting was performed with GraphPad Prism v.5.01.

Co-localisation analyses were conducted to assess effects of genetically predicted gene expression mRNA levels from 49 tissues obtained via eQTLs from the GTEx v8 database (EUR sample, https://gtexportal.org/) on stiffness index. To increase the robustness of these analyses, two different co-localisation methods were applied, focussing on the intersection of the significant results.

For both methods, the associations with the stiffness index of all SNPs within 1.1 Mb around rs603424, as well as the eQTLs of the corresponding regions per tissue, were extracted. First, Bayesian co-localisation analyses were conducted using the R-package “gtx” version 2.1.6 (https://github.com/tobyjohnson/gtx, ‘coloc.fast’ function with 100 kb SNP window and default parameters and prior definitions), which implemented the co-localisation method of Giambartolomei et al. [20]. For all co-localisation analyses, a posterior probability (PP) of ≥ 0.80 of the H4 test (both trait and expression data are associated and share the same single causal variant) was applied to identify significant results.

Second, the SNP rs603424 was tested and plotted for co-localisation with the tissue-specific mRNA levels by applying the summary-data-based Mendelian randomisation (SMR) method [21]. The method includes a test whether the effect on expression observed at a SNP is independent of the signal observed in the trait association (SMR test) and a second test that evaluates if the eQTL and trait associations can be attributable to the same causative variant by performing a heterogeneity test (HEIDI test). Significance for co-localisation of the gene expression and the trait signals was defined by pSMR < 0.001, where additionally a pHEIDI ≥ 0.05 indicates the same underlying causal variant.

Finally, we assessed the associations between rs603424 and the amount of SAT, VAT, and the ratio of VAT/SAT in SHIP-START-2 and SHIP-TREND-0. Cohort-specific linear regression analyses with log-transformed adipose tissue markers as outcome and rs603424 as exposure were calculated. The adjustment of the models followed the adjustment in the genome-wide association study. Following this, the results were combined by a fixed effects inverse-variance weighted meta-analysis analogue to the GWAS meta-analysis.

Results

General characteristics of the SHIP participants and the HSD and OsteoGene patients are listed in Table 1. In the two patient cohorts, women were overrepresented, while the sex ratio was balanced in the SHIP cohorts. Patients in the OsteoGene cohort were older (average age 66.0 years) than patients in the HSD cohort (56.2 years), in SHIP-START-2 (56.8 years) and SHIP-TREND-0 (50.9 years). Among HSD and OsteoGene patients, secondary osteoporosis was diagnosed in 48.9% and 51.3%, respectively. Moreover, patients evaluated for osteoporosis had on average a lower BMI than SHIP participants. Fractures were reported by less than 10% and intake of vitamin D or calcium supplements by less than 2% of SHIP participants. These values were expectedly higher in the patients (Table 1). The QUS-based stiffness index was comparable between SHIP-START-2 and SHIP-TREND-0, while BMD at the spine and femoral neck was lower in patients from the HSD cohort than in OsteoGene patients.

Table 1 Characteristics of SHIP participants and patients evaluated for osteoporosis
Full size table

The associations between the 1103 SNPs reported by Morris and Kemp et al. [13] and the QUS-based stiffness index were examined in the two SHIP cohorts, and the results were combined in a meta-analysis. This analysis, including 2108 SHIP-START-2 and 3557 SHIP-TREND-0 participants (total of 5665), yielded 45 significant associations after correction for multiple testing (Table 2). Among the associated SNPs, there were two in the WNT16 locus (rs2707519 and rs3779381), as well as SNPs in EN1, JAZF1, PKD2L1, SPTBN1, GPC6, TMEM135, and further loci. From this list, we selected the five SNPs with the lowest p-values for further analyses in patients evaluated for osteoporosis: rs2707518 (CPED1/WNT16), rs3779381 (WNT16, intron 1), rs115242848 (LOC101927709/EN1), rs10239787 (JAZF1, intron 2), and rs603424 (PKD2L1, intron 2). We additionally chose rs6968704 (JAZF1, intron 2) to receive more information about the JAZF1 locus.

Table 2 SNPs associated with stiffness index after correction for multiple testing. Results from the meta-analysis combining SHIP-START-2 and SHIP-TREND-0 (n = 5665)
Full size table

The meta-analysis combined the effect of the minor alleles of both SHIP cohorts of each of the 1,103 SNPs on the stiffness index. A positive effect size indicates an increase in stiffness index, and a negative effect size a decrease in stiffness index per minor allele. SNPs whose minor alleles were carried by less than 5% of the SHIP participants (MAF below 0.05) had the greatest effects on stiffness index (effect size below −3.0 or above +3.0), while SNPs with more frequent minor alleles had lower effects on stiffness index (MAF above 0.05 and effect sizes between −2.2 and +2.5) (Fig. 1). The selected SNPs (rs2707518, rs3779381, rs115242848, rs10239787, rs603424, rs6968704) had MAFs between 0.013 and 0.370 and effect sizes between −1.96 and +6.70 (Fig. 1, Table 1).

Fig. 1
figure 1

Effect size according to minor allele frequency (MAF) of all analysed 1103 SNPs in the meta-analysis (n = 5665 SHIP-START-2 and SHIP-TREND-0 participants). Orange dots represent SNPs that were significantly associated with the stiffness index after correction for multiple testing. Orange dots that were named represent the SNPs that were selected for further analysis in the patients evaluated for osteoporosis. Blue dots represent SNPs that were not significantly associated with the stiffness index

Full size image

Subsequent analysis in patients of the HSD or OsteoGene studies revealed no statistically significant associations between the six genotyped SNPs and femoral neck or spine BMD. Also in analyses combining HSD and OsteoGene, no statistically significant results were obtained (Supplemental Table 1). Yet we observed trends toward differences in MAF between SHIP participants (representative of the general population) and the patients of the HSD and OsteoGene studies (Fig. 2 and Supplemental Table 2). While 24.0% of the SHIP participants carried the minor allele of rs3779381, there were only 22.9% minor allele carriers among the patients. The minor allele was positively associated with stiffness index in our meta-analysis [β =  +2.22, standard error (stderr) = 0.37] indicating a protective effect on BMD. Comparable observations were made for the protective minor alleles of rs2707518 and rs698704, although the differences in MAF between SHIP participants and patients evaluated for osteoporosis were below 0.3%. Also, rs115242848 was associated with an increased stiffness index in our meta-analysis. Its MAF was very low in SHIP (1.30%) but even lower in the patients (0.49%). Rs603424 on the other side was inversely associated with stiffness index (β = −1.96, stderr = 0.42), and fewer minor allele carriers were observed among the SHIP participants (17.8%) than among the HSD and OsteoGene patients (19.5%). Only rs10239787, which was inversely related to stiffness index (β = − 1.59, stderr = 0.42) had a higher MAF in SHIP participants (35.8%) than in the patients (35.4%). When assessing the two patient cohorts separately, we observed that all trends for differences in MAF were present in the OsteoGene cohort, while in the HSD cohort, only the results for rs3779381, rs115242848, and rs603424 were stable (Supplemental Table 2 and Supplemental Fig. 1).

Fig. 2
figure 2

Comparison of the minor allele frequency (MAF) of the selected SNPs in the general population in SHIP-START-2 and SHIP-TREND-0 (denoted as SHIP), and the MAF of patients evaluated for osteoporosis in HSD and OsteoGene (denoted as patients). We hypothesised that a minor allele with protective effect on BMD (as inferred from the meta-analysis) is less represented in patients evaluated for osteoporosis than in individuals from the general population, whereas a minor allele with a negative effect on BMD is more represented. Green colour indicates confirmation of the hypothesis, grey colour indicates no confirmation

Full size image

To test whether the genetic effect of a locus on BMD was mediated by gene expression of a nearby gene, we applied two co-localisation methods. These methods test for an association of gene expression and stiffness index estimated by SNP-mRNA and SNP-BMD association results. Both co-localisation methods indicated that higher mRNA expression of stearoyl-CoA desaturase (SCD) in visceral adipose tissue is associated a higher stiffness index (Fig. 3, Supplementary Tables 3 and 4). Our analyses further revealed that rs603424 affects SCD expression. Rs603424 is located in intron 2 of PKD2L1, a neighbouring gene of SCD. Its minor allele was associated with a decreased SCD mRNA expression in visceral adipose tissue. Thus, the minor allele of rs603424 may tag for effects on ultrasound-based heel stiffness index via changes in SCD gene expression. To follow-up on this, we examined whether rs603424 was associated with the amount of SAT, VAT or the ratio of VAT/SAT. Respective analyses yielded, however, no statistically significant results (Supplemental Table 5).

Fig. 3
figure 3

Co-localisation results. Illustration of the summary data-based Mendelian randomisation (SMR) test for the stiffness index and expression quantitative trait loci (eQTL) in adipose visceral tissue at the rs603424 locus. The upper box shows the regional association plot of the genetic associations with the stiffness index, with level of significance of the SMR test (y-axis) for each mRNA transcript in the locus indicated by a diamond positioned at the corresponding gene. A significant SMR test represented by a purple diamond indicates an association of the transcript level of the respective genes (purple label) with the trait. The filled purple diamond indicates a HEIDI test p-value > 0.05, thus, a likely co-localisation. The lower box shows the regional association distribution with changes in expression of the highlighted (purple) gene transcript. In both boxes, the x-axis refers to GRCh37/hg19 genomic coordinates

Full size image

Discussion

In the present analyses, we confirmed several loci associated with QUS-based bone properties in 5665 adults from the general population. Six of these SNPs, mapping to WNT16, EN1, JAZF1, and PKD2L1, were genotyped among 631 patients evaluated for osteoporosis. Tendencies towards differences in MAF between individuals from the general population and patients evaluated for osteoporosis support effects of five of these SNPs on bone substance. Moreover, co-localisation analyses pointed to a causal effect of gene expression of SCD on heel bone stiffness index (for a more detailed description of the respective genes see Supplemental Table 6).

While previous studies described in total more than a thousand genetic variants related to BMD [8, 11], eBMD [13], or fractures [22], it largely remains unclear what the causal mechanisms underlying the observed associations are. Methods to uncover causality include the identification and co-localisation of eQTL [23]. These approaches generate insights into the (patho-)physiology of bone metabolism. We complemented results obtained from a GWAS meta-analysis with eQTL co-localisation results and assessed MAFs of the promising variants in individuals from the general population and two cohorts of patients evaluated for osteoporosis. While the most prominent association with bone stiffness in our data was observed for two known SNPs in the WNT16 locus, eQTL co-localisation analyses pointed to a causal role of rs603424 in bone integrity.

WNT16 (rs2707518, rs3779381)

The Wnt-signalling pathway is crucial for the maintenance of bone homeostasis [24]. WNT16 encodes Wnt-16 and is part of the WNT gene family, which was identified to positively regulate osteoblast [25] and negatively regulate osteoclast differentiation [26]. It is highly expressed in osteoblasts of cortical bone [26] and was shown to be crucial for the preservation of cortical [26] and trabecular bone mass [27]. Wnt-16 has been identified as a target of glucocorticoid action with Wnt-16 suppression resulting in decreased bone formation [28]. Genetic variations in WNT16 were previously associated with BMD and fracture risk at genome-wide significance [29, 30]. Among these variants, two missense variants rs2707466 (Thr > Ile) and rs2908004 (Gly > Arg) were repeatedly shown to be related to BMD [31, 32].

In the present work, we focused on two non-coding SNPs, rs2707518 (CPED1/WNT16) and rs3779381 (WNT16), which were highly associated with eBMD in the UK Biobank [13] and with stiffness index in the SHIP cohorts. Rs2707518 is highly genetically linked (r2 = 0.99) to rs2536195 (CPED1/WNT16) which previously showed an association with stiffness index in a GWAS meta-analysis [33]. Rs3779381, in turn, was observed to be related to an increased osteoporosis risk in postmenopausal, overweight Chinese women [29]. Among the patients evaluated for osteoporosis, both SNPs demonstrated a lower MAF (0.366 for rs2707518 and 0.229 for rs3779381). The observed difference could, thus, represent an indicator of patient selection. Also, the effect direction observed in SHIP is in line with the results from the UK Biobank [13]. The minor T (rs2707518) and G (rs3779381) alleles were associated with increased stiffness index or eBMD. The minor alleles may, thus, confer a protective effect on bone properties, and carriers of these protective minor alleles may be less often found among patients than among the general population. Therefore, our data infer that these non-coding variants in the WNT16 locus may play an important role in maintaining bone integrity.

PKD2L1 (rs603424)

Morris and Kemp et al. [13] were the first to demonstrate a genome-wide significant association of rs603424 with adult eBMD. Their analyses suggest a potentially harmful effect of that SNP on eBMD [13], which was confirmed by our finding of an inverse association with the QUS-based stiffness index. The observation of a tendency towards fewer minor allele carriers among the SHIP participants (17.8%) than among the patients evaluated for osteoporosis (19.5%) further provides support for a potentially harmful effect of this SNP.

Co-localisation analyses suggest that rs603424 may impact on bone stiffness via modification of stearoyl-CoA desaturase (SCD) expression in adipose tissue. SCD is a key enzyme in lipogenesis that catalyses the synthesis of saturated fatty acids to monounsaturated fatty acids [34]. In humans, two SCD isoforms have been identified, SCD1 and SCD5. While SCD5 is mainly expressed in the brain and pancreas, SCD1 is more ubiquitously expressed, e.g. in adipose tissue, liver, brain, heart, and pancreas [34]. Over the last years, several studies reported SCD1 activity to be related to a disturbed lipid metabolism in obesity and non-alcoholic fatty liver disease as well as tumour malignancy [34, 35]. Yet the actions of SCD1 are complex and still not fully elucidated. Thus, SCD1 was described to exert beneficial effects by acting anti-inflammatory while inhibition of SCD1 increased saturated fatty acid levels and inflammation [36]. More important for the present work, SCD1 was reported to be involved in bone homeostasis, by promoting osteogenic differentiation of bone marrow mesenchymal stem cells [37]. In line with this, a connection between SCD and the Wnt pathway was identified as SCD1 and SCD2 were shown to provide a necessary lipid modification for Wnt biogenesis and pathway activation [38].

While our data point to an association of rs603424 with bone properties via affecting SCD1 expression, an association between rs603424 and SAT, VAT or the ratio of VAT/SAT was not significant (all p-values > 0.58). However, it must be noted that only about half of the SHIP participants underwent MRI and provided data for the analyses, which strongly limits statistical power. Moreover, it might be awarding to assess further measures of quantity or quality of adipose tissue and to unravel mechanistic insights of SCD1. Our results presented here, will further need follow-up in functional analyses.

EN1 (rs115242848)

Within our analyses, rs115242848 showed the largest effect on the stiffness index (β =  +6.702). This low-frequency SNP is an intergenic variant close to EN1 which encodes the homeobox protein engrailed-1 [39]. Styrkarsdottir and colleagues [40] previously observed an association of rs115242848 with hip and spine BMD. In the SHIP data, the minor T allele of rs115242848 was associated with higher stiffness index values, indicating a protective effect on bone health. Quite similar, genome-wide significant effects on total body BMD were reported from the Life-Course GWAS meta-analysis [41]. Moreover, associations of rs115242848 with increased lumbar spine BMD [42], lumbar spine area [40], forearm and femoral neck BMD [42], and inverse associations with osteoporosis [43] and fractures [13] were reported but missed genome-wide significance. The MAF of rs115242848 was low among the SHIP participants (0.013) but higher than in other populations (0.010 in 1000G Europe, 0.008 in TWINSUK). The MAF in patients evaluated for osteoporosis (0.005) was even lower. Taking this low MAF and, thus, the low number of individuals carrying this variant in the HSD and OsteoGene cohort into account, no definite conclusions on the association of rs115242848 with BMD can be made.

JAZF1 (rs10239787, rs6968704)

JAZF1 encodes the transcriptional corepressor JAZF Zinc Finger 1 [44, 45]. Previously, associations of JAZF1 with 124 different traits have been reported [46]. Among these, associations with anthropometric measures are dominating, but associations with diabetes, prostate cancer, and heel BMD have also been found [46]. The two, in our data, replicated SNPs had an opposite effect on the QUS-based stiffness index. While rs10239787 was related to lower values (β = −1.59), rs6968704 was associated with higher values (β =  +1.43). The potentially deleterious impact of rs10239787 on eBMD is in line with the results of Morris and Kemp et al. [13]. An inverse association of rs10239787, but without genome-wide significance, was further reported with total body BMD [41] and a corresponding positive association with fractures [47]. Regarding rs6968704, we are the first to independently confirm genome-wide significance with bone properties. Moreover, our results suggest a tendency towards higher MAF of the protective minor allele in the general population than in patients evaluated for osteoporosis. This implies that individuals lacking this protective variant may be overrepresented among the patients. This notion is, however, quite speculative, especially as previous studies demonstrated positive associations of rs6968704 with total body BMD [41] and inverse associations with lumbar spine area [48] but both without genome-wide significance.

Strengths and Limitations

The present study stands out due to its sample of well-characterised individuals from a general population of European ancestry, and the simultaneous complementary evaluation of two samples of patients evaluated for osteoporosis. One limitation is that in the population-based SHIP cohorts, dual-energy X-ray absorptiometry (DXA) scans, using ionising radiation, were impossible due to ethical objections. This prohibited a direct comparison of BMD between SHIP participants and HSD or OsteoGene patients. Moreover, due to differences in bone structure throughout the skeleton, the QUS-based measurements at the heel are not directly comparable to measurements obtained at the femoral neck or spine. Nevertheless, QUS measurements allow osteoporotic fracture risk prediction [6]. The main results from our genetic analysis, i.e. the six selected SNPs, were fostered by the comparison of the MAFs of the minor alleles between SHIP participants and HSD and OsteoGene patients. Yet we cannot exclude that the use of the two methods to evaluate bone traits prevented a full confirmation in the patients evaluated for osteoporosis, i.e. significant associations in the regression models. We can further not rule out that differences in patient characteristics between HSD and OsteoGene contributed to more stable MAF differences in the OsteoGene cohort, than in the HSD cohort. Although there are differences, the majority of results for the HSD and OsteoGene cohorts points in the same direction and allows drawing similar conclusions. Above this, our study has the limitation that, in contrast to the study by Morris and Kemp et al. [13], it may not be large enough to estimate the effect of rare variants and might be underpowered to identify respective associations. It is further not large enough to conduct analyses stratified by sex, menopausal status, or aetiology of osteoporosis. Analyses in stratified, less heterogeneous samples, generally have a higher statistical power and yield more precise effect estimates due to the reduction of data variability. This advantage is, however, outweighed by the reduction of sample size, which substantially decreases statistical power. Despite these limitations, our analyses among the heterogeneous patient cohorts generally support the results from the genome-wide association study. This, in turn, strongly argues for the robustness of our results. The genetic effects obtained from the population-based cohorts assuming a commonly used additive genetic association model were reflected in MAF differences when compared to the patient cohorts, although the analyses were underpowered to reach statistical significance. Furthermore, we could not distinguish between dominant or recessive effects.

Conclusion

Taken together, our results confirm 45 genetic variants to be associated with ultrasound-based heel stiffness index in adult men and women from the general population. Among these, two non-coding variants in the WNT16 locus (rs2707518, rs3779381) as well as a novel, possibly causal association of rs603424 mediated via SCD, were demonstrated to impact the examined bone parameter. Our results, thus, highlight the effect of the Wnt-16 pathway in the regulation of bone properties and indicate a role of SCD expression in adipose tissue on bone substance.