Mendelian randomization analysis of lipids for BMD and fracture riskData sources

Single nucleotide polymorphisms (SNP) related to the human lipidome were selected as IVs from our previous genome-wide association study (GWAS) (Fig. 1 and Graphical abstract) [16]. The study was a large-scale GWAS of the human serum lipidome data from 6,057 participants (individuals of European ancestry). The summary-level data for heel BMD from the UK biobank were selected (IEU Open GWAS Project ID: ebi-a-GCST90014022), which included 365,403 individuals of British ancestry [17]. Heel BMD was assessed using quantitative ultrasound (Hologic Sahara system) in the UK Biobank, with standardized T-scores derived from speed of sound and broadband ultrasound attenuation. The summary-level data for osteoporotic fracture was obtained from the UK Biobank (GWAS Catalog ID: GCST006980), which included 53,184 British ancestry cases and 373,611 British ancestry controls [18]. Osteoporotic fracture cases were identified through both hospital records using International Classification of Diseases (ICD)-10 codes (excluding traumatic/pathological fractures) and questionnaire-based self-reports of fractures within the past 5 years. There was no overlap between the populations used for exposure (lipidome) and outcomes (BMD and fracture).

Fig. 1

Flowchart of this study design. MR: Mendelian Randomization; BMD: Bone Mineral Density; GWAS: Genome-Wide Association Study; SNP: Single Nucleotide Polymorphism; IVW: Inverse-Variance Weighted; FDR: False Discovery Rate; TASOAC: Tasmanian Older Adult Cohort; LC–MS/MS: Liquid Chromatography–Tandem Mass Spectrometry

Instrumental variable selection

SNPs associated with lipid class/species at the genome-wide significance threshold (p < 5 × 10–8) were selected as potential IVs. The reference panel for assessing linkage disequilibrium (LD) between these SNPs consisted of the European contingent from the 1000 Genomes Project. The LD threshold for clumping was established at r2 < 0.001, with a clumping window size set to 10,000 kb. SNPs with minor allele frequency (MAF) ≤ 0.01 were excluded, along with palindromic A/T or G/C alleles. If harmonization was not available, or the number of extracted SNPs was fewer than 3, then IVs were re-extracted using a threshold of p < 5 × 10–6.

Observational analysis in Tasmanian Older Adult Cohort (TASOAC) studyStudy population

The TASOAC Study, a prospective population-based cohort study, recruited 1,099 individuals aged 50 to 80 years at baseline, who were selected randomly from the electoral roll in Southern Tasmania [19]. Participants then underwent follow-up assessments at 2.6 years (phase 2), 5.1 years (phase 3), and 10.7 years (phase 4). The number of attendees at each follow-up was 875, 768, and 566, respectively. Blood samples were obtained at the 2.6- and 10.7-year follow-up and used for the lipidomics assay. This study comprised 492 and 310 participants with complete data on BMD, fracture, lipidomics, and covariates at the 2.6- and 10.7-year follow-ups (Fig. 1 and Graphical abstract). Approval for the TASOAC study was granted by the Southern Tasmanian Health and Medical Human Research Ethics Committee (Ref. no: H0006488), and written informed consent was obtained from all participants.

Lipidomics profiling

Lipidomics profiling was measured at both 2.6- and 10.7-year follow-ups. Blood samples were obtained following a minimum 8-h fasting period and subjected to centrifugation at 13,000 xg for 10 min. This fasting duration aligns with standard protocols for lipidomic profiling and is considered sufficient to reduce postprandial variation in circulating lipid levels, as validated in previous studies using liquid chromatography tandem-mass spectrometry (LC–MS/MS) [20, 21]. Serum samples were randomized, and lipid extraction was performed using the butanol/methanol method, as described previously [22]. Lipidomic analysis was conducted using an Agilent 6495C triple quadrupole LC-MS/MS in conjunction with an Agilent 1290 series HPLC system and a ZORBAX eclipse plus C18 column (Agilent, Santa Clara, CA, USA) operating in positive/negative switching polarity, as detailed previously [22, 23]. More details were given in Supplementary file 1.

Hip BMD

The dual-energy X-ray absorptiometry (DXA) (Hologic Delphi densitometer, MA, USA) was employed to conduct the total hip BMD scans at both 2.6- and 10.7-year follow-ups.

Fractures

Fractures were self-reported and documented at both 2.6- and 10.7-year follow-ups. Participants were asked: “List any fracture you may have had since your previous interview for this study. Please list these by the location of the fractures (including vertebral, non-vertebral, hip and major fractures involving the femur, radius, ulnar, vertebral, rib and humerus).”

Covariate measurements

Age, sex, and current smoking status were collected through a questionnaire during an interview at 2.6 years. Standing height and body weight were measured, and body mass index (BMI, kg/m2) was calculated. Physical activity was measured by recording steps per day over seven consecutive days utilizing a pedometer (Omron Healthcare, Kyoto, Japan) at 2.6 years. The detailed criteria for inclusion of pedometer estimates have been described previously [24]. Serum 25-hydroxyvitamin D (25(OH)D) concentrations were measured using liquid-phase radioimmunoassay (Immunodiagnostics Systems Ltd), which detects both 25(OH)D2 and 25(OH)D3. The intra-assay and inter-assay coefficients of variation were 1.8% and 3.3%, respectively [25]. Falls risk scores were derived from the physical profile assessment (PPA) at 2.6 years. The PPA is a validated and reliable tool that assesses five physiological domains: visual contrast sensitivity, reaction time, knee extension strength, proprioception, and postural sway on a foam surface, and standardized Z-scores based upon these five domains, were then calculated [26]. Statin use was self-reported at the 2.6-year follow-up. Energy-adjusted dietary inflammatory index scores were calculated based on validated food frequency questionnaire data collected at baseline, as previously described [27].

Statistical analysisTwo-sample MR and sensitivity analyses

In the two-sample MR analysis, the estimation of the causal effect of exposure on the outcome was carried out using the inverse-variance weighting (IVW) method. This involved combining the ratio estimates for each SNP and effectively translating MR estimates into a weighted regression of SNP-outcome effects on SNP-exposure effects [28]. The fixed-effect IVW was used as the primary analysis method. In cases where there was significant heterogeneity among the IVs, the multiplicative random effect inverse-variance weighted (IVW-MRE) method was utilized to analyze associations. MR-Egger serves as a tool for estimating causal effects through the slope coefficient from Egger regression, and it also can detect small study bias and certain forms of pleiotropy [29]. The weighted median method offers unbiased estimates, even when up to 50% of the information is derived from invalid IVs [30]. Furthermore, the weighted mode approach assumes that the most frequently occurring association estimate remains unaffected by pleiotropy, indicating that it corresponds to the true causal effect [31].

Cochran’s Q statistic, MR-PRESSO analysis, F-statistic and Steiger filtering test method were described in Supplementary file 1.

For all lipid classes, we initially performed MR analysis at the lipid classes-level to explore their associations with BMD/fracture. Subsequently, within the identified significant lipid classes, further MR analyses were performed at lipid species-level, by examining their relationships with the outcomes. The false discovery rate (FDR) correction was used to control multiple testing for analyses at the lipid species-level.

The two-sample MR analyses were performed using the TwoSampleMR (version 0.5.8), Mendelian randomization (version 0.8.0), and MR-PRESSO (version 1.0) in R Software 4.3.2 (https://www.R-project.org). FDR-corrected p-value < 5 × 10–2 was considered statistically significant after controlling multiple testing.

Descriptive analysis in the validation cohort

Participants' characteristics were described using the mean (standard deviation [SD]) for continuous variables and the percentage (number) for categorical variables. The concentrations of lipid classes and species were natural log transformed for normalization and used for analyses.

Associations between lipids and hip BMD over 8 years

Linear mixed-effects models were used to analyze the associations between classes and species of lipids and hip BMD measured at both the phase 2 and phase 4, after adjusting for age, sex, BMI, physical activity, current smoking status, serum levels of vitamin D at phase 2. The interaction terms of lipids and follow-up time in the same linear mixed-effects models were used to estimate the associations between lipid classes and lipid species and changes in hip BMD over 8 years.

Associations between lipids and fracture risk over 8 years

Log binominal mixed-effects models were used to analyze the associations between classes and species of lipids and fracture risk, which were recorded at both the phase 2 and phase 4. These models were adjusted for age, sex, BMI, physical activity, current smoking status, falls risk score, and time-dependent hip BMD. The interaction terms of lipids and follow-up time in the same log binomial mixed-effects models were used to estimate the associations between classes and species of lipids and fracture risk over 8 years.

Sensitivity analyses

Sensitivity analyses were performed by adjusting for statin use and dietary inflammatory index at the 2.6-year follow-up. Additionally, inverse probability weighting method was applied to account for potential attrition bias.

Stata software (V.17) (Stata Corp, College Station, TX, USA) was used to conduct the observational analyses. Statistical significance was defined as a two-tailed p-value < 5 × 10–2.