UK Biobank is a prospective cohort study that recruited more than half a million participants aged 40–69 years in the United Kingdom between 2006 and 2010 [18]. At baseline, participants were assessed for their demographic characteristics, lifestyle, medical histories, and women’s reproductive factors, through touchscreen questionnaires and oral interviews. In addition, they underwent various physical examinations such as measurements of weight, height, waist circumference, and hip circumference. Bioelectrical impedance was measured by the Tanita BC418MA body composition analyzer, which uses an eight-electrode setup that estimates body composition in the whole body as well as in the limbs [17]. VAT mass was measured by dual-energy X-ray absorptiometry (DXA) in a subset of 5,109 individuals using the GE Healthcare Lunar iDXA scanner.

Using the methods developed by Karlsson T et al., we estimated the predicted measurements of VAT mass (VAT^) from the multivariable-adjusted model, for men and women separately [17]. The following variables were included in the prediction model: age, menopause status, height, weight, waist circumference, hip circumference, impedance of left arm, impedance of right arm, impedance of left leg, impedance of right leg, and impedance of the whole body. The training dataset included 5,109 participants with VAT mass measured by DXA, and the application dataset included the whole study sample of the baseline. To avoid potential bias due to population structure, only White British participants were included in the study, and those who had ‘Do not know’, ‘Prefer not to answer’ or ‘NA’ variables in the prediction model were set as missing the corresponding variables [21]. After removal, 396,403 participants passed the initial quality control and obtained the VAT^, and the following participants were further excluded from the analysis: (1) individuals with extreme waist (≤ 26 cm) or hip circumference compared to their body size measurements (hip > 2 × BMI + 80 cm); (2) individuals with a difference in impedance between the left and right arm or the left and right leg of more than 120Ω. Detailed information about the exclusion criterion of participants has been published elsewhere [17]. In total, 395,612 participants constituted the final application dataset (Additional File 2: Fig. S1).

Heel BMD was estimated based on ultrasound measurement of the calcaneus during baseline assessment, using a Sahara Clinical Bone Sonometer (Hologic Corporation, Bedford, MA, USA) [19]. After quality control, a total of 483,939 individuals were included in the analysis (Fig. 1). Incident fractures were identified according to the ICD-10 codes at hospital admissions (Additional file 1: Table S3) [19]. Participants with fractures of the skull, face, hands, and feet (usually due to trauma), pathological fractures (due to malignant tumors), atypical femoral fractures, periprosthetic fractures, and fracture healing were excluded from the analysis [19]. Traumatic fractures were not excluded because the cause of the trauma was not well captured by the ICD-10 codes. Fractures at the hip and vertebrae are most likely attributed to low BMD [20, 21]. Therefore, hip, vertebrae, and other fractures were included as outcomes for specific sites of fractures in this study. The date of fracture was determined from hospital admission data. Date of death was determined from death certificates held by the National Health Service (NHS) Information Centre for participants from England and Wales, and the NHS Central Register Scotland for participants from Scotland [18]. In the current analysis, hospital admission data were available until March 31, 2017, and mortality data were available until February 14, 2018. Therefore, we censored follow-up on March 31, 2017, or the date of first fracture or death, whichever occurred first. In total, 451,023 participants without baseline fractures were included in the analysis (Fig. 1).

Flow diagram of study design. VAT^, predicted values of VAT (visceral adipose tissue) mass; BMD bone mineral density; DXA dual-energy x-ray absorptiometry; FN-BMD femoral neck BMD; LS-BMD lumbar spine BMD; FA-BMD forearm BMD; MR mendelian randomization; and IV instrumental variable

Baseline touchscreen questionnaires were used to assess several potential confounding variables, including age, household income, smoking status, alcohol consumption, physical activity, dietary supplements, overall health rating, medical history including the histories of type 2 diabetes (T2D), cardiovascular disease (CVD), or cancer diagnosed by a doctor, and menopausal status and hormone replacement therapy (women only). Lean mass and body weight were measured with the Tanita BC-418 MA body composition analyzer (Tanita Corporation of America, IL). Standing height was measured using the Seca 202 device (SECA, Hamburg, Germany). BMI was calculated as weight in kilograms divided by height in meters squared.

According to gender-specific quartiles of VAT^, baseline characteristics were presented as percentage (%) for categorical variables and mean (standard deviation, SD) for continuous variables, for men and women separately. Heel BMD was normalized using the R ‘scale’ function [22] before analysis. VAT^ was examined as a continuous variable with the unit kg in all models. The regression coefficient (β) and standard error (SE) for the association of BMD with VAT^ were assessed by using the linear regression models. The hazard ratios (HRs) and 95% confidence interval (CI) for fracture risk associated with VAT^ were assessed by using the Cox proportional hazards models. The proportional hazards assumption of Cox models was tested based on the Schoenfeld residuals and almost no evidence of violation was observed (Additional file 1: Table S9). We also used restricted cubic splines with five knots at the 5th, 25th, 50th, 75th, and 95th centiles to flexibly model the association of VAT^ with heel BMD and fracture risk. We tested for potential nonlinearity by using a likelihood ratio test comparing the model with only a linear term against the model with linear and cubic spline terms [23, 24]. All models were adjusted for age (years), annual household income (< £18,000, £18,000-£30,999, £31,000-£51,999, £52,000-£100,000, or > £100,000), lean mass (kg), standing height (cm), smoking status (never, previous, or current), alcohol consumption (less than 3 times/month, 1–2 times/week, or more than 3 times/week), physical activity (metabolic equivalent task [MET] -min/week), intakes of calcium and vitamin D supplements (both as yes or no), overall health rating (excellent, good, fair, or poor), baseline disease history including T2D, CVD, and cancer (all as yes or no), and for women, menopause status (premenopausal or postmenopausal) and use of hormone replacement therapy (yes or no). All analyses were conducted by gender (men and women), and we further grouped men and women by BMI (< 25 kg/m2 and ≥ 25 kg/m2) in these observational analyses. The interactions of sex and BMI with VAT^ on skeletal outcomes were estimated by including a multiplicative factor in all models. To minimize the potential confounding effects and reverse causality brought by diabetes, CVD, and cancer [25,26,27,28,29], we also conducted additional sensitivity analyses by excluding participants who had self-reported physician-diagnosed T2D, CVD, or cancer at baseline. P < 0.05 (two-tailed) was considered to be statistically significant. All analyses were performed using R (version 3.6.1).

In the main analyses, we performed two-sample MR analyses (Fig. 1) using individual-level data from UK Biobank for VAT^ and summary results data for the skeletal outcomes. A genome-wide association (GWA) analysis was performed in the UK Biobank (n = 396,612) to identify instrumental variables of VAT^ (Additional File 3: Method S1). BMD-associated SNPs were identified from 3 separate GWAS summary statistics of European participants’ femoral neck bone mineral density (FN-BMD, n = 32,735), lumbar spine bone mineral density (LS-BMD, n = 28,498), and forearm bone mineral density (FA-BMD, n = 8143) [30]. For fractures, we used summary GWA results of hip fractures (the number of cases = 3,983; the number of controls = 211,460; GWAS ID: 'finn-b-ST19_FRACT_FEMUR') and vertebrae fractures (the number of cases = 2,859; the number of controls = 212,839; GWAS ID: ' finn-b-ST19_FRACT_LUMBAR_SPINE_PELVIS') from the FinnGen Biobank and the largest GWA results of total fractures from a meta-analysis of 25 cohorts [31] (the number of cases = 37,857; the number of controls = 227,116). The random-effect inverse-variance weighted (IVW) method was applied as the primary analysis [32]. We also performed MR Egger [33], weighted median [34], and mode-based regressions as sensitivity analyses. MR analyses were performed using the R ‘TwoSampleMR’ package [35].