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 (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 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].

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.

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.

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.

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).

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).

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.