The study was performed among participants of the prospective population-based Cooperative Health Research in the Region of Augsburg (KORA) study. A total of 4,261 adults, aged 25-74 years, were included at baseline between 1999 and 2001 (S4 visit) with the primary aim to assess health and disease in Southern Bavaria, Germany. Follow-up examinations were conducted after 7 years (F4 visit, 2006 -2008) and after 14 years (FF4 visit, 2013 - 2014) (20–22). All study participants have provided written informed consent. The study was approved by the Ethics Committees of the Bavarian Chamber of Physicians (Ethical Approval Number 06068) adhering to the declaration of Helsinki.

The present analysis includes data from the F4 visit as baseline and FF4 visit as follow-up (average follow-time: 6.5 years). Excluding premenopausal women (n = 602), women with hysterectomy or bleeding due to hormone replacement therapy and younger than 60 years (n = 188), women with missing menopausal status (n = 4), participants without valid FLI information at baseline (n = 47), a total of 2,239 participants (1,456 men and 783 postmenopausal women) were included in the cross-sectional analysis ( Figure 1 ). Due to missing sex hormones information at baseline (n = 60 to 244), the final number of participants for the regression analyses differed by sex hormone (1,328 to 1,417 men; 667 to 762 postmenopausal women) at baseline. For the longitudinal analysis, we further excluded participants lost to follow-up (n = 720) and those without FLI information (n = 14) at the FF4 visit, leaving a sample size of 1,505 participants (941 to 1,003 men; 408 to 468 postmenopausal women).

Details of laboratory, clinical and anthropometric measurements as well as interviews are provided in the Supplementary Materials .

T, DHEA, DHEAS, DHT, progesterone, and 17-OHP were quantified in serum samples which were stored at -80°C until being assayed. The detailed assessment procedure has already been described in detail (23). Samples were prepared and sex hormones were quantified using the AbsoluteIDQ™ Stero17 Kit and electrospray ionization liquid chromatography-mass spectrometry (ESI-LC-MS/MS). The quantification method of the AbsoluteIDQ™ Stero17 Kit has been proved to follow the European Medicines Agency’s Guideline on bioanalytical method validation (July 21^st 2011) (24). Metabolite concentrations were calculated using internal standards and reported in nM or ng/ml. Missing values of sex hormones were imputed (11). Sex hormones were then normalized, and different batches were calibrated (11). SHBG was measured in serum using the chemiluminescent microparticle immunoassay ARCHITECT for the absolute quantification of SHBG (Abbott Diagnostics).

In order to be transported in blood, sex hormones are bounded to SHBG or weakly bounded to albumin. The free circulating sex hormones [e.g. free T (fT), free DHT (fDHT)] represent the bioactive hormones that target tissues. The sum of albumin-bound and free sex hormone is bioavailable sex hormone (e.g. bioavailable T). In the KORA study, fT and fDHT were calculated using mass action equations based on the concentrations of the total hormones and their binding constants to serum SHBG and albumin according to Rinaldi et al. (25).

FLI was calculated from BMI, waist circumference, triglycerides (TG) and gamma-glutamyl-transferase (GGT) with the algorithm developed by Bedogni et al. (14):

FLI = (e ^{0.953*loge (TG) + 0.139*BMI + 0.718*loge (GGT) + 0.053*waist circumference - 15.745})/(1 + e ^{0.953*loge (TG) + 0.139*BMI + 0.718*loge (GGT) + 0.053*waist circumference - 15.745}) * 100, with TG measured in mmol/l, GGT in U/l, and waist circumference in cm, resulting in a score ranging from 0 to 100, with a FLI < 30 ruling out and a FLI ≥ 60 ruling in fatty liver.

We searched the GWAS Catalogue using the full name of the sex hormones and identified the largest GWAS in the European population including total T, bioavailable T (bioT), E2, SHBG (18, 26), DHEAS (17), progesterone and 17-OHP (19). For DHT, the only GWAS in the European population was conducted in a study population at particular high risk of prostate cancer (men only, n= 3225) (27). One GWAS identified a SNP (rs34670419) associated with DHEA in a European population (n=1023); however, the association (p=2e-9) did not reach a genome-wide cut-off of p<5e-11 after multiple-testing adjustment (28). Therefore, we did not include DHT and DHEA in the MR analysis. Summary statistics for total T in men and women, bioT in men and women, E2 in men, SHBG in men and women (Ruth et al., 2020 (18)), DHEAS in men and women combined (Zhai et al., 2011 (17)), and progesterone in men and women, 17-OHP in men and women (Pott et al., 2021) (19) were obtained from the respective publications. Of note, a genetic instrument for E2 in women was not included, as in the GWAS of Ruth et al. most of the women were postmenopausal and showed E2 levels below the limit of detection (78%), which substantially reduced the power of analysis for genetic instruments of E2 and biased the associations towards loci associated with age at menopause (18).

After we included the genome-wide significant SNPs (p<5e-8) for sex hormones and SHBG, we clumped the SNPs if they were in linkage disequilibrium (LD) (LD r²>0.001). The SNP-outcome associations were extracted from the largest GWAS available up to date for MRI measured hepatic proton-density fat fraction (PDFF) in the European population by Parisinos et al. using data from a subsample of UK Biobank (29) ( Supplementary Table 1 ). We chose this study because MRI has been demonstrated as the most definitive non-invasive medical imaging to quantify liver fat content (30). Afterwards, we harmonized the SNP-exposure and SNP-outcome associations and excluded palindromic SNPs ( Supplementary Figure 1 ). Three genetic instruments that could not be matched in the outcome dataset (rs543504257, rs2275560, rs78058190) were excluded from further MR analysis.

Baseline characteristics of the participants were compared among the FLI categories stratified by sex. For continuous variables, the arithmetic mean and standard deviation (SD) are shown if normally distributed or the median and interquartile (IQR) if non-normally distributed. For categorical variables, counts and percentages (%) were displayed. Analysis of variance (ANOVA) was used for continuous variables and chi-square test was used for categorical variables to test the differences between the groups.

Sex-specific correlations of sex hormones were examined by Pearson’s rho. Sex hormone concentrations were sex-specifically z-standardized. The associations between sex hormones and baseline FLI as well as FLI at the follow-up were investigated with linear regression stratified for men and postmenopausal women. Model adjustment was defined a priori. The main model was adjusted for age, conventional lifestyle and cardiometabolic risk factors for sex hormone derangement and ectopic fat accumulation, including smoking (never, ex-smoker, smoker), physical activity (active, inactive), alcohol consumption (no intake, moderate intake, excessive intake), systolic blood pressure (SBP), high-density lipoprotein-cholesterol (HDL-C), low-density lilpoprotein-cholesterol (LDL-C) (all continuous), clinically diagnosed diabetes (yes, no), use of antihypertensive medication (yes, no) and lipid lowering medication (yes, no) ( Supplementary Material ). The model for the longitudinal analysis was additionally adjusted for baseline FLI. In a sensitivity analysis, we further adjusted for continuous C-reactive protein (CRP), thyroid stimulating hormone (TSH), serum albumin and SHBG, which are either closely related to sex hormone derangement or determinant for bioavailable sex hormones. The significance level was set to p<0.0056 to account for multiple tests (9 exposures) using Bonferroni correction.

MR analysis was conducted to investigate the causal relationship between sex hormones/SHBG and hepatic PDFF. More detailed explanation of the methodology is provided in the Supplementary Material . Firstly, we conducted MR analysis with the inverse-variance weighting (IVW) approach. One of the MR assumptions (exclusion restriction) is that the association between the genetic instrument and the outcome goes only through the exposure ( Supplementary Materials ). The IVW approach is only valid if all genetic instruments fulfill the “exclusion restriction” assumption. In case of genetic pleiotropy where the “exclusion restriction” is violated and genetic variants are also associated with other risk factors of the outcome, other robust MR methods provide valid and consistent MR estimates. The weighted-median approach allows up to 50%, the weighted-mode approach 50% - 100% and the MR-Egger approach up to 100% for pleiotropic variants (31, 32). A statistically significant IVW result with directionally consistent MR estimates from all three sensitivity analyses was considered to be a potential causal effect (33). The existence of directional horizontal pleiotropy was defined if the intercept term of the MR-Egger regression significantly differed from zero (p for pleiotropy < 0.05).

For DHEAS, progesterone and 17-OHP, we conducted two-sample MR analysis. Whereas for total T, bioT, E2 and SHBG, MR analysis was carried out in a two-sample setting with population overlap (<10%), since summary statistics for both exposure and outcome were obtained from the UK Biobank. In large scale studies, the precision and bias of MR estimates (except for MR Egger approach) are similar in both two-sample or one-sample (with complete sample overlap) MR settings (34, 35).

In order to investigate the causal effect of T or E2 independent of SHBG, we conducted MR analysis using clusters of genetic instruments with primary effects on specific sex hormone (T or E2 or SHBG) identified previously by Ruth and colleagues ( Supplementary Table 1 ) (18). We also conducted sensitivity analysis excluding the SNPs with larger effects on the metabolic risk factors closely related to fatty liver, including fasting glucose, type 2 diabetes, coronary artery disease, HDL-C, LDL-C, triglycerides, total cholesterol, SBP, DBP, BMI and waist-to-hip ratio adjusted for BMI, than their effects on sex hormones identified by Steiger-Filtering previously (18). The Steiger test filters out SNPs that explain more variance in one phenotype (e.g. outcome/trait closely related to the outcome) than another (e.g. exposure), to reduce potential pleiotropic effects of these SNPs and avoid reverse causality (36).

All analyses were conducted using R statistical software, version 4.2.1, including the MR analyses for which we used the “TwoSampleMR” R package. We used multiple imputation with 5 imputed datasets for covariates with missing values less than 5% for the observational analysis. In the MR analysis, a p<0.0071 (0.05/7 exposure) was considered significant with Bonferroni correction for multiple testing.

Among 2,239 participants eligible for the cross-sectional analysis, the prevalence of FLI ≥ 60 was higher in men (54%, mean age 57 years) than in postmenopausal women (38%, mean age 66 years). For both men and women, participants with higher FLI were significantly older, had higher BMI, larger waist circumference and they were less physically active. They had higher blood pressure, higher blood lipid concentrations, higher CRP levels and higher liver enzyme levels. They also suffered more from diabetes, and were more likely to use antihypertensive or lipid lowering medication ( Table 1 , Supplementary Table 2 ). Among men, lower levels of sex hormones and SHBG were seen with higher FLI. Whereas among postmenopausal women, higher fT and lower DHEA, DHT and SHBG concentrations were observed with higher FLI ( Table 2 ). A correlation matrix between the sex hormones and SHBG is shown in Supplementary Figure 2 .

Multivariable adjusted regression analyses showed that among men, lower T [β, 95%CI: -4.89 (-6.12, -3.66)], DHT [-2.97 (-4.20, -1.73)], progesterone [-2.75 (-4.02, -1.49)], 17-OHP [-3.57 (-4.80, -2.34)] and SHBG [-4.64 (-5.89, -3.39)] were associated with higher FLI at baseline. Among postmenopausal women, higher fT [2.27 (0.77, 3.77)] and lower SHBG [-9.00 (-11.13, -6.87)] were associated with higher FLI at baseline ( Figure 2 and Supplementary Table 3 ). In longitudinal analysis, similar trends followed for both men and women ( Figure 2 and Supplementary Table 4 ). In the sensitivity analysis, additionally adjusting for CRP, TSH, serum albumin and SHBG hardly changed the associations ( Supplementary Table 5 ).

All associations in the longitudinal analysis were attenuated after adjustment for baseline FLI ( Supplementary Table 4 ), possibly due to reverse causation. However, baseline adjustment is only occasionally advantageous, and whether it eliminates or introduces bias depends crucially upon the causal structure relating the variables (37).

For sex hormones, the MR IVW estimates for total T [-0.09 (-0.16, -0.01)] in men and bioavailable T [0.13 (0.03, 0.23)] in women were nominally significant (p<0.05), but they did not pass the significance level of p< 0.0071 after Bonferroni correction. MR analyses with the IVW approach revealed that higher SHBG among women [-0.36 (-0.61, -0.12)] was associated with lower hepatic PDFF. Among men, the estimate was smaller [-0.19 (-0.33, -0.05)], and did not pass the Bonferroni threshold. Sensitivity analyses with weighted median, weighted mode and MR-Egger yielded estimates directionally consistent to the IVW estimates ( Table 3 ). There was no indication of directional horizontal pleiotropy in the above MR analyses (p for pleiotropy from MR-Egger ≥ 0.05) ( Table 3 ).

Due to genetic overlap between T, E2 and SHBG, we used clusters of instrumental SNPs with primary effects on T or E2 or SHBG to investigate the potential causal effect of T or E2 on hepatic PDFF independent of SHBG. MR analysis with clusters of T or E2 showed that there was no association between T and hepatic PDFF independent of SHBG in either men or women. Nor was there any association between E2 and hepatic PDFF independent of SHBG in men ( Supplementary Table 6 ). The IVW estimates for both male SHBG cluster [-0.20 (-0.34, -0.06)] and female SHBG cluster [-0.43 (-0.61, -0.25)] reached statistical significance after Bonferroni correction. All three sensitivity analyses resulted in estimates in the same direction as the IVW estimates ( Supplementary Table 6 ). The male SHBG cluster includes SNPs with primary SHBG increasing effect and secondary increasing effect on total T and decreasing effect on bioT as well as increasing effect on E2, and the female SHBG cluster includes SNPs with primary increasing effect on SHBG and secondary opposing effect on T and bioT. Taken together, this indicated that genetically determined higher SHBG has a decreasing effect on hepatic PDFF in both men and women, probably also through its effect on sex hormones ( Table 3 ).

In order to minimize the pleiotropic effect of SNPs closely associated with metabolic risk factors, we further excluded them from the MR. The association between SHBG and hepatic PDFF attenuated, but maintained the same directionality ( Supplementary Table 7 ).