Vitamin D (VD) is a micronutrient for human health, which is predominantly derived from 7-dehydrocholesterol in the skin under sunlight exposure and is also absorbed from food as a secondary source. Through two steps of hydroxylation, VD is converted to 25-hydroxyvitamin D (25OHD) and then 1,25-dihydroxyvitamin D [1,25(OH)₂D]. The circulating concentration of 25OHD is generally used as a biomarker of VD status, and 1,25(OH)₂D is the functional metabolite working in the manner of a hormone in the osseous and extra-osseous metabolisms^[1]. Previous studies suggested that VD deficiency was associated with abdominal obesity and/or dyslipidemia^[2], and many studies supported the opinion that VD supplementation was beneficial for ameliorating lipid profiles^[3]. However, confounders in various populations led to inconsistency for each of the lipid parameters, namely triglyceride (TG), total cholesterol (TC), low-density lipoprotein cholesterol (LDLC), and high-density lipoprotein cholesterol (HDLC), as suggested by human studies and meta-analyses^{[3, 4]}. Researchers have to conduct a close examination of those potential confounders and pay special attention to the functioning of the VD receptor (VDR), which is the most important mediator of the VD function^[5].

VDR gene variants have the potential to influence biological outcomes^{[5, 6]}. It has been reported that VDR has an effect on skeletal muscles independent of VD, and the decreased expression of VDR protein is related to various disease states and aging^[7]. In the VDR knockout animal model, a decreased size and muscle fiber strength were observed, and VDR expression is also required in myoblasts^{[8, 9]}. Aside from skeletal health status, it was also reported that the single nucleotide polymorphisms (SNPs) of VDR were associated with risks to non-skeletal health^[6], such as myocardial infarction, cancer, and death^[10]. Moreover, these kinds of associations aroused concerns about dyslipidemia since expressions of some lipid metabolism-related genes are regulated by calcitriol-stimulation^[11], and some may even have VD-responsive elements (VDREs) to bind the 1,25(OH)₂D-activated VDR^[12]. Thus, both the VD status and the VDR function have the potential to affect the expressions of lipid metabolic genes, ultimately affecting the physiological activities of cells^[1]. Moreover, by involving biological (in terms of reproductive function, sexual hormones, expressions of heterosome-located genes, etc.) and lifestyle (e.g., smoking, drinking, physical activity, etc.) factors^[13], sexual dimorphism should be taken into careful consideration when conducting medical research^[14].

Located on chromosome 12, the VDR gene spans at least 75 kb in length. From the 5’–3’ direction of its sense (forward) strand, Cdx2 (rs11568820, C > T), Fok1 [rs2228570, A > B (degenerate base standing for C, G, and T, i.e., not A)], Apa1 (rs7975232, C > A), and Taq1 (rs731236, A > G) are the generally the relevant SNPs in various ethnicity groups^[15-17]. These SNPs may affect mRNA expression, amino acid sequence, and protein activity. Particularly, VDR SNPs have been reported to have effects on 25OHD levels in circulation^{[18, 19]}. For example, Santos et al.^[20] studied a group of Brazilian girls aged 7–18 years and found that the VDR wild-types of Bsm1, Apa1, and Taq1 were associated with lower 25OHD levels. Therefore, it is necessary to further explore the relationship between VDR SNPs and specific lipid parameters, as well as VD status in both females and males. In the present study, we focused on the VDR SNPs of Apa1, Cdx2, Fok1, and Taq1 in this regard, which have been studied in genotype-phenotype associations, in order to avoid weak evidence or unknown biological significance for less studied or new SNPs regarding their associations with lipids.

The sample size (N) was estimated with the formula N ≥ (Z_1-α/2/δ)² × p × (1 – p), wherein Z was 1.96 for the two-sided 95% confidential intervals (CIs), P was the dyslipidemia prevalence of 34.64% in Shenzhen adults as known in 2012 but unpublished at the design stage of our study^[21], δ was 2.2%, representing tolerable error, and the calculated N was 1,798.

After signing written consent forms to be subjects in the study, adult volunteers visiting a health examination center in Shenzhen City in Guangdong Province, China, were recruited between July 2013 and January 2014. A questionnaire survey administered via face-to-face interview was conducted to collect the basic health information of each participant. The inclusion criteria were: 1) Han Chinese aged ≥ 20 years old; 2) living in Shenzhen for > 2 years; 3) free of liver diseases, renal diseases, and any cancers in the past six months; and 4) not pregnant (for women). The excluded subjects were those who 1) had severe organic diseases, 2) had acute infection symptoms, allergic diseases, and malignant tumors, 3) had a family history of genetic diseases including familial dyslipidemia, 4) had taken VD supplements in the past six months, 5) had taken medicines to control lipid levels within the past 12 hours, or 6) had taken diuretics, engaged in strenuous exercise, or had overeaten within the 24 hours before the test.

From ulnar veins of the subjects, who had fasted overnight, blood samples were collected with vacuum tubes, and within two hours post-collection, supernatant and blood cells were separated by centrifugation of 3,000 ×g at 4 °C for 10 min. The body mass index (BMI) was calculated as body weight (kg)/ height (m)². The above protocol for the cross-sectional study was approved by the Ethics Committee of the Shenzhen Center for Chronic Disease Control and was in accordance with the Declaration of Helsinki. The study was registered at ClinicalTrials.gov (registration ID number NCT04707612).

Parameters of lipids, glucose, hemogram, etc., were assayed right on the day, and the rest of the samples were aliquoted with EP tubes and stored at -80 °C for later use. The systolic blood pressure (SBP) and diastolic blood pressure (DBP) of the subjects were measured by a standard mercury sphygmomanometer. Glycated hemoglobin (HbA1c) was detected with a Japan Arkray Instruments analyzer (ion chromatography method). The serum levels of TG, TC, LDLC, HDLC, and fasting plasma glucose (FPG) were tested by the Beckman-LX20 automatic biochemical analyzer. The plasma 25OHD concentration was detected with a commercial ELISA kit (Cat. #: AC-57F1, IDS Ltd., UK). According to the Chinese guidelines for dyslipidemia management^[22], TG ≥ 2.3 mmol/L, TC ≥ 6.2 mmol/L, LDLC ≥ 4.1 mmol/L, HDLC < 1.0 mmol/L, or previously diagnosed dyslipidemia, were the criteria used to define the dyslipidemia group (DL), and those subjects who did not match the DL criteria were assigned to the non-dyslipidemia group (ND).

In the logistic analyses, from the viewpoint of preventive medicine, as well as to avoid the small subject numbers for the genotypes of low minor allele frequency (MAF)^[23], the available marginally-elevated cut-off values for TG (≥ 1.7 mmol/L), TC (≥ 5.2 mmol/L), and LDLC (≥ 3.4 mmol/L), as well as the only available cut-off value for HDLC (< 1.0 mmol/L), were used to define the abnormal subgroups for the corresponding lipid parameters.

The protocols for DNA preparation and SNP genotyping described in our previous studies^{[23, 24]} were also adopted in the present study. Briefly, an asymmetric amplification with molecular beacon-based real-time quantitative PCR (MB-qPCR) followed by a melting step was performed to determine the variation of the SNP loci. Before its application to the large sample size analysis, the method was verified with the gold standard of Sanger sequencing (ThermoFisher, Shanghai, China). The information of primers and molecular beacons for the genotyping experiment with qPCR (Roche 480II, Singapore) is presented in Supplementary Table S1, available in www.besjournal.com.

The nucleotide bases used in our study were consistent with those in dbSNP of NCBI, i.e., C > T for rs11568820 (Cdx2), A > B (not A) for rs2228570 (Fok1), C > A for rs7975232 (Apa1), and A > G for rs731236 (Taq1). If the complementary bases [G > A for Cdx2, T > V (standing for A, C, and G) for Fok1, G > T for Apa1, and T > C for Taq1], or the letter standing for the DNA-digestibility with the corresponding restriction enzyme (f > F for Fok1, a > A for Apa1, and T > t for Taq1), were used in the references, the original usages were provided as bracketed annotations in our article.

The statistical analyses were performed using SPSS for Windows version 25 (IBM Corp., Armonk, NY, USA). Clinical data were presented as medians and interquartile ranges (25% to 75%) and compared with rank-sum tests (Wilcoxon rank test or Kruskal-Wallis rank test). The genotypes of the four SNPs were tested with Hardy-Weinberg equilibrium (HWE) analyses for sampling representation. For genotypic comparisons, differences in allele and genotype frequencies were evaluated using the chi-square (χ²) test. The homozygous genotypes of CC for Cdx2, AA for Fok1, CC for Apa1, and AA for Taq1 were used as reference genotypes, respectively. For both genders, the additive, dominant, recessive, homozygous, and allelic models for each of the SNPs were entered into the logistic regression analyses for odds ratios (ORs) and 95% CIs with adjustment for age, BMI, FPG, HbA1c, 25OHD, SBP, and DBP. In order to analyze the interaction between VD nutritional levels and VDR SNPs, the interaction factor calculated by 25OHD concentrations and VDR SNPs were entered into the logistic regression with adjustment for those confounders mentioned above. A simple linear regression was used to assess the association between circulating 25OHD and lipid profiles. A P-value of less than 0.05 was considered to be statistically significant.

A total of 1,987 adults were included in the analysis, aged from 20 to 81 years old. Among them, there were 1,124 females aged 39 (31–49) (median and interquartile range) years old and 863 males aged 36 (30–44) years old. The clinical profiles of the participants were divided into DL and ND for each gender, as summarized in Table 1. It was indicated that several metabolic or metabolism-related parameters were statistically different between DL and ND adults (P < 0.05), such as age, BMI, FPG, HbA1c, SBP, DBP, TG, TC, LDLC, HDLC, etc. In particular, 25OHD concentrations were statistically different between DL and ND groups in men but not in women.

The results of the genotyping experiment on the SNPs with MB-qPCR are shown in Supplementary Figure S1, available in www.besjournal.com. The HWE test showed that the population sampled from the Health Examination Center was representative of the population at large (Supplementary Table S2, available in www.besjournal.com). The genotypic and allelic frequencies of VDR SNPs between DL and ND groups are summarized in Supplementary Table S3, available in www.besjournal.com. In the female participants, both genotypic and allelic frequencies of Fok1 showed differences between DL and ND groups (P < 0.05); while in the male participants, these differences were not observed. The genotypic and allelic frequencies of the other three SNPs, namely Cdx2, Apa1, and Taq1, did not display significant differences between the DL and ND groups for either gender (P ≥ 0.05), and neither did their different genetic models (Supplementary Tables S4–S6, available in www.besjournal.com).

Figure S1.  Melting curve analyses on the polymorphisms of vitamin D receptor gene (VDR). Panel A: Melting curve of VDR Cdx2 (rs11568820, A). Panel B: Melting curve of Fok1 (rs2228570, B). Panel C: Melting curve of VDR Apa1 (rs7975232, C). Panel D: Melting curve of VDR Taq1 (rs731236, D).

Logistic regression analyses were performed with adjustment for age and BMI since the adjustment for 25OHD did not change the overall correlations (see Supplementary Tables S7–S8). For female participants (see Table 2), between the subgroups of LDLC (≥ 3.4 vs. < 3.4 mmol/L), VDR Fok1 presented significance (P < 0.05) in its additive (BB vs. AB vs. AA, P = 0.03, OR = 1.28, 95% CI: 1.03–1.59), recessive (BB vs. AB + AA, P = 0.03, OR = 1.44, 95% CI: 1.04–2.00), homozygous (BB vs. AA, P = 0.02, OR = 2.89, 95% CI: 1.18–7.05), and allelic (B vs. A, P = 0.04, OR = 1.26, 95% CI: 1.01–1.59) models, but between the subgroups of TC (≥ 5.2 vs. < 5.2 mmol/L) and HDLC (< 1.0 vs. ≥ 1.0 mmol/L), no significance (P ≥ 0.05) was found in all models of the four SNPs.