A total of 3023 subjects from a community in Nanjing (Jiangsu, China) were enrolled in this case-control study during July to September 2018. The NAFLD cases were recruited from those diagnosed based on “Guideline of prevention and treatment for nonalcoholic fatty liver disease: a 2018 update”. NAFLD can be diagnosed if the following items 1-4 coexist with the fifth or sixth item: (1) no drinking or history of overdose drinking (less than 210g/week ethanol for men and 140g/week for women in the past 12 months); (2) excluding drug hepatitis, hepatitis C virus genotype 3 infection, hepatolenticular degeneration and other specific diseases that could result in fatty liver; (3) serum levels of transaminase and γ-glutamyl transpeptidase (γ-GT) increase mildly to moderately (<5 times above the upper normal limit), usually presenting as an increase of alanine aminotransferase (ALT); (4) metabolic syndrome constituents such as visceral obesity, hyperglycemia, blood lipid disorder and hypertension; (5) the result of liver imaging study meets the imaging diagnostic criteria of diffuse fatty liver; (6) the histological findings of liver biopsy meet the pathological diagnostic criteria of fatty liver disease. The non-NAFLD controls were collected from the same community during the study period and randomly assigned to the control group. The constituent ratios of gender and age between cases and controls were considered similar, according the results of frequency-matching.

Excluded were: (1) less than 18 years old; (2) subjects with infection, acute or chronic gastrointestinal diseases, autoimmune diseases or malignant tumors; (3) subjects with history of other viral hepatitis, alcoholic liver disease or primary liver cancer; (4) subjects of excessive drinking (alcohol consumption ≥30g/d in males and ≥20g/d in females); (5) subjects who received a liver transplant within the previous year, or had complications of advanced liver disease (varicose veins, ascites, etc.); (6) subjects with drug-induced fatty hepatitis; (7) subjects with history of psychiatric disorders.

By reviewing the literature, we assumed the frequency of gene mutation in the general population was 20%, odds ratio (OR) was 1.5, two-sided test α was 0.05, and power of test (1-β) was 90%. Therefore, the minimum sample size was estimated by NCSS-PASS 11 software (Dawson edition; Kaysville, UT) to be 784. This study had a sample size large enough to guarantee the production of reliable results.

The current study protocol was in accordance with the Declaration of Helsinki, and was approved by the Institutional Ethics Review Committee of Nanjing Medical University (Nanjing, China). Written informed consent was obtained before blood test and genetic analysis.

Self-designed questionnaires and an electronic medical record system were used to collect the demographic and clinical characteristics of all participants. 5-mL ethylene diamine tetraacetic acid (EDTA) anticoagulant venous blood was collected from each subject in the morning after an overnight fast. Within 2 hours, the serum and blood cells in each blood sample were separated and frozen at -80°C until further serological tests and genotyping assays.

In two groups (gender and age matched at a ratio of 1:1, 336 subjects in the NAFLD group and 336 in the control group), the serum 25(OH)D₃ level was measured by enzyme-linked immunosorbent assay (ELISA) according to the manufacturer’s instructions (Human 25-Dihydroxy vitamin D3 (25(OH)D3) ELISA Kit; Jin Yibai Biological Technology Co., Ltd.; Nanjing, China).

Target SNPs with potential biological function were selected to better identify the susceptibility sites that affect the development of NAFLD. The selection processes were as follows: (1) Downloading the genotype database of VDR gene in Han Chinese in Beijing (CHB) from the 1000 Genomes Project database (http://asia.ensembl.org/Homo_sapiens/Info/Index); (2) Importing the genotype database into the Haploview software (version 4.2; Broad Institute, Cambridge, MA, USA). Setting parameters: Hardy-Weinberg P-value cutoff = 0.05; minor allele frequency (MAF) = 0.05; r² threshold = 0.8. At this point, 56 tagging SNPs (tagSNPs) were captured; (3) Searching for literature in which these potential tagSNPs were associated with NAFLD, hypertension, type 2 diabetes, hyperlipidemia, metabolic syndrome, inflammatory diseases or immune-related disorders. Finally, eight disease-related target SNPs of VDR gene were selected, including rs3782905 (C>G), rs3847987 (C>A), rs11574129 (T>C), rs2228570 (C>A), rs11568820 (G>A), rs739837 (G>T), rs7975232 (C>A) and rs11168287 (G>A).

Genomic DNA was isolated from EDTA anticoagulated blood samples using magnetic bead method (blood genomic extraction kit; Pangu Genome Nanotechnology Co., Ltd.; Nanjing, China). The eight SNPs were genotyped using TaqMan allelic discrimination assays on the ABI 9700 system (Applied Biosystems, Foster City, California, USA; catalog numbers: C:_3290647_10, C:_2404006_10, C_175992105_10, C:12060045_20, C:_2880808_10, C:_2404007_10, C:28977635_10, C:_2404006_10). The quality control of experimental data was as follows: (1) blinding was adopted in genotyping, so that all laboratory personnel were unclear about the clinical data of the subjects; (2) 10% of the samples were randomly selected for repeated experiments with a repeatability of 100%. The success rates of genotyping all SNPs were higher than 99% in this study.

All statistical analyses were processed using SPSS (version 22.0, SPSS Inc., Chicago, IL, USA) and MedCalc (Version 19.1, Ostend, Belgium). Distributions of demographic and clinical characteristics of two groups were compared using χ ²-test, student’s t test or Mann-Whitney U test, wherever appropriate. Logistic regression analysis was used to calculate odds ratio (OR) and 95% confidence interval (95% CI) for quantifying the association of serum 25(OH)D₃ level with the risk of NAFLD. The correlation of serum 25(OH)D₃ level with VDR SNPs was assessed using general linear regression model adjusted for gender and age. The Hardy-Weinberg equilibrium (HWE) was tested using a goodness of-fit χ ² -test among the control subjects. The relationships between VDR SNPs and the risk of NAFLD were analyzed by dominant model (heterozygote + mutant homozygote VS. wild homozygote), recessive model (mutant homozygote VS. wild homozygote + heterozygote), and additive model (mutant homozygote VS. heterozygote VS. wild homozygote), respectively. For multiple SNPs comparisons, false discovery rate (FDR) correction was used and the P _FDR value ≤ 0.25 was regarded as modest confidence that the correlation represented a positive result (30). Subgroup analysis was performed for positive SNPs, and Q test was performed to calculate the heterogeneity between subgroups. The area under the receiver operating curves (AUROCs) was performed to assess the predictive power of the combination of positive SNPs, environmental behavior factors and mainly related clinical factors for NAFLD risk. A two-tailed test with a P value < 0.05 was regarded as statistically significant in all analyses.

The demographic and clinical characteristics of 1120 NAFLD cases and 1903 controls were summarized in Table 1 . No significant differences were observed in the distribution of gender and age between the two groups (all P>0.05). However, there were significant differences in exercise time, body mass index (BMI), waist circumference (WC), systolic blood pressure (SBP), diastolic blood pressure (DBP), triglyceride (TG), total cholesterol (TC), high density lipoprotein-cholesterol (HDL-C), low density lipoprotein-cholesterol (LDL-C), glucose (GLU), γ-glutamyl transpeptidase (γ-GT), alanine aminotransferase (ALT), aspartate aminotransferase (AST), direct bilirubin (DBIL), total bilirubin (TBIL) and serum 25(OH)D₃ levels (all P<0.05).

According to the level of serum 25(OH)D₃, 672 participants were divided into three categories: VD deficiency (<20ng/mL), VD insufficiency (20ng/mL-30ng/mL) and VD sufficiency (≥30ng/mL) (31). The constituent ratios of three categories were 60.9% (409/672), 21.4% (144/672), and 17.7% (119/672), respectively. The prevalences of NAFLD in the three categories were 59.2% (242/409), 32.6% (47/144) and 39.5% (47/199), respectively, with significant differences (χ ² = 36.366, P<0.05). In addition, multiple comparisons showed that the prevalence of NAFLD in patients with VD deficiency was significantly higher than that in patients with VD insufficient or VD sufficient (all P _Bonferroni<0.017), but there was no significant difference in the prevalence of NAFLD between categories with VD insufficiency and VD sufficiency (P _Bonferroni>0.017) ( Supplementary Figure 1 ).

Further logistic regression analysis suggested that compared to VD deficiency, VD sufficiency (crude OR=0.450, 95%CI=0.297-0.684, P<0.001) and VD insufficiency (crude OR=0.334, 95%CI=0.224-0.499, P<0.001) were associated with a lower risk of NAFLD, and this association was dose-dependent (P _trend<0.001) (model 1). After adjusting for gender, age, BMI and SBP (model 2), the subjects with VD sufficiency (adjusted OR=0.477, 95%CI=0.277-0.821, P=0.008) and VD insufficiency (adjusted OR=0.344, 95%CI=0.206-0.573, P<0.001) all had significantly decreased risk of NAFLD compared with those with VD deficiency, and a significant locus-dosage effect of VD level on NAFLD risk was also observed (P _trend<0.001). When adjusting for the above factors and TG, HDL-C, GLU (model 3), the risk of NAFLD still decreased significantly with the increase of VD level (VD sufficiency vs. VD deficiency: adjusted OR=0.552, 95%CI=0.305-0.998, P=0.049; VD insufficiency vs. VD deficiency: adjusted OR=0.298, 95%CI=0.165-0.539, P<0.001; P _trend=0.005) ( Table 2 ).

Serum 25(OH)D₃ levels were Lg transformed into an approximately normal distribution, and SNPs were coded in an additive genetic model. The general linear regression analysis indicated that there were no significant associations between eight VDR SNPs and serum 25(OH)D₃ levels after adjusting for gender and age (all P>0.05) ( Supplementary Table 1 ).

The genotype distributions of the eight SNPs in both groups were shown in Table 3 . After adjusting for gender and age, logistic regression analyses showed that VDR rs2228570-A variant (AA: adjusted OR=0.782, 95%CI=0.633-0.966, P=0.023; dominant model: adjusted OR=0.837, 95%CI=0.710-0.986, P=0.034; additive model: adjusted OR=0.883, 95%CI=0.794-0.981, P=0.020) and rs11168287-A variant (GA: adjusted OR=0.830, 95%CI=0.707-0.974, P=0.022; dominant model: adjusted OR=0.839, 95%CI=0.721-0.976, P=0.023) all significantly reduced the risk of NAFLD in different models. Under the premise of FDR threshold ≤ 0.25, FDR revealed with modest confidence that the associations between the two SNPs and a low risk of NAFLD were positive (all P _FDR=0.136, in dominant model; Supplementary Table 2 ).

In addition, the combined effects of VDR variants rs2228570 and rs11168287 on the risk of NAFLD were estimated by the number of favorable alleles from the two SNPs, as shown in Supplementary Table 3 . The subjects were first divided into three groups with “0”, “1-2” and “3-4” favorable alleles. Compared with those who had 0 favorable alleles, subjects with “3-4” favorable alleles had significantly decreased risk of NAFLD (adjusted OR=0.750, 95%CI=0.577-0.974, P=0.031) and the more favorable alleles, the lower NAFLD risk, suggesting a significant locus-dosage effect of the combined alleles on NAFLD risk (P _trend=0.039). The subjects were then distributed into two groups with “0” and “1-4” favorable alleles, and we found that the presence of “1-4” alleles was related to a 0.798-fold lower risk of NAFLD (adjusted 95% CI=0.644-0.990, P=0.040).

We further performed the stratification analyses to evaluate the combined effects of VDR rs2228570-A and rs11168287-A alleles on the risk of NAFLD adjusted with gender and age. As shown in Supplementary Figure 2 , the combined protective effects of two alleles were more prominent in subjects ≤ 40 years (adjusted OR=0.696, 95%CI=0.515-0.940, P=0.018), non-hypertension (adjusted OR=0.770, 95%CI=0.602-0.986, P=0.038), non-hyperglycemia (adjusted OR=0.791, 95%CI=0.631-0.991, P=0.042) and non-low HDL-C (adjusted OR=0.710, 95%CI=0.536-0.942, P=0.018). The heterogeneity test discovered no significant differences among all the subgroups (all P>0.05).

A stepwise regression model were performed with gender, age, visceral obesity, hypertension, hyperglycemia, hypertriglyceridemia, Low HDL-C, ALT, exercise time, rs2228570 and rs11168287. The coding of each variable was described in Supplementary Table 4 . The results showed that age ≤40 years (OR=0.716, 95%CI=0.583-0.879, P=0.001), exercise time ≥150 min/week (OR=0.798, 95%CI=0.661-0.963, P=0.019) and rs2228570-A (OR=0.768, 95%CI=0.626-0.942, P=0.011) were independent protective factors of NAFLD. Conversely, visceral obesity (OR=3.653, 95%CI=2.984-4.471, P<0.001), hypertension (OR=1.654, 95%CI=1.290-2.121, P<0.001), hypertriglyceridemia (OR=3.455, 95%CI=2.784-4.288, P<0.001), Low HDL-C (OR=1.879, 95%CI=1.536-2.299, P<0.001) and ALT > 40U/L (OR=2.729, 95%CI=2.060-3.615, P<0.001) were independent risk factors of NAFLD ( Table 4 ).

We further constructed combined factors based on the above factors for assessing NAFLD risk, and the combined factor 1 (a combination of clinical factors) = (-1.588) + (-0.229) × age + 1.333 × visceral obesity + 0.484 × hypertension + 1.233 × hypertriglyceridemia + 0.621 × Low HDL-C + 0.939 × ALT. Similarly, combined factor 2 (a combination of clinical factors and exercise time) = (-1.448) + (-0.254) × age + 1.294 × visceral obesity + 0.512 × hypertension + 1.222 × hypertriglyceridemia + 0.610 × low HDL-C + 0.958 × ALT + (-0.220) × exercise time; combined factor 3 (a combination of clinical factors, exercise time and SNP genetic variant) = (-1.519) + (-0.334) × age + 1.296 × visceral obesity + 0.503 × hypertension + 1.240 × hypertriglyceridemia + 0.631 × low HDL-C + 1.004 × ALT + (-0.225) × exercise time + (-0.264) × rs2228570. In Figure 1 , the AUROCs of these combined factors were 0.770, 0.774, and 0.780, respectively, and the AUROC of combined factor 3 for assessing NAFLD risk was slightly higher than that of combined factor 2 (0.780 vs. 0.774, 95%CI=0.002-0.011, P=0.007) and combined factor 1 (0.780 vs. 0.770, 95%CI=0.005-0.016, P<0.001).