Vitamin D receptor gene FokI but not TaqI, ApaI, BsmI polymorphism is associated with Hashimoto’s thyroiditis: a meta-analysis

Four VD receptor (VDR) gene polymorphisms (TaqI, ApaI, FokI and BsmI) have been reported to influence Hashimoto’s thyroiditis (HT) risk. However, individual studies have produced inconsistent results. We conducted a comprehensive meta-analysis of eleven case-control studies to better understand roles of the four polymorphisms in HT development. The results showed only FokI polymorphism was significantly associated with the risk of HT (F vs f: OR = 1.44, 95% CI = 1.09–1.91, P = 0.010; FF vs Ff + ff: OR = 1.72, 95% CI = 1.09–2.70, P = 0.019). Subgroup analyses demonstrated the significant effect was only present in Asian population (F vs f: OR = 1.45, 95% CI = 1.07–1.95, P = 0.016; FF vs ff: OR = 1.64, 95% CI = 1.03–2.59, P = 0.036; FF + Ff vs ff: OR = 1.34, 95% CI = 1.00–1.80, P = 0.047; FF vs Ff + ff: OR = 1.64, 95% CI = 1.03–2.64, P = 0.039), but not in Caucasian. For TaqI, ApaI and BsmI polymorphisms, no significant association was found in any model comparison. Based on the current literature, it appears that only VDR FokI polymorphism is associated with HT risk in Asian population, but not in Caucasians; and the TaqI, ApaI and BsmI polymorphisms have not positive association neither in the overall population, nor when stratified by ethnicity. Further well-designed studies with larger sample sizes and different ethnic population are needed to clarify the present findings.


Meta-analysis results.
provides the pooled results regarding the association of the four VDR gene polymorphisms and HT risk under five different genetic models, along with the P-value of Egger's test for publication bias.
FokI polymorphism. Eight studies including 978 cases and 938 controls examined the association of FokI polymorphism and HT risk. Pooled analyses showed a significant association in the allele model (F vs f: OR = 1.44, 95% CI = 1.09-1.91, P = 0.010) and the dominant model (FF vs Ff + ff: OR = 1.72, 95% CI = 1.09-2.70, P = 0.019), but not in the other models (Table 3, Fig. 2). Significant heterogeneity existed in these two models (I 2 = 69.2%, and P = 0.002 for allele model; I 2 = 75.7%, and P = 0.000 for dominant model). Then, Galbraith plot analyses were performed to further explore the sources of heterogeneity. As shown in Fig. 3A and C, the studies performed by Guleryuz et al. 38 and Meng et al. 36 might mainly contribute to the heterogeneity. With exclusion of these studies, the heterogeneity decreased significantly (I 2 = 0% and P = 0.760 for F vs f; I 2 = 0% and P = 0.738 for FF vs Ff + ff) while the overall association remained significant in these two models (F vs f: OR = 1.72, 95% CI = 1.42-2.07, P = 0.000; FF vs Ff + ff: OR = 2.32, 95% CI = 1.79-3.02, P = 0.000) ( Fig. 3B and D). There was one study 35 the genotype distributions in controls departed from HWE. Sensitivity analyses by excluding this study did not change the pooled result of allele model (F vs f: OR = 1.37, 95% CI = 1.03-1.82, P = 0.030), but the P value of the dominant model was borderline (FF vs Ff + ff: OR = 1.54, 95% CI = 0.98-2.43, P = 0.060). Subgroup analyses by ethnicity indicated that the FokI F allele or FF genotype significantly increased the risk of HT in Asians (F vs f: OR = 1.45, 95% CI = 1.07-1.95, P = 0.016; FF vs ff: OR = 1.64, 95% CI = 1.03-2.59, P = 0.036; FF + Ff vs ff: OR = 1.34, 95% CI = 1.00-1.80, P = 0.047; FF vs Ff + ff: OR = 1.64, 95% CI = 1.03-2.64, P = 0.039), but the positive association was not found in Caucasians. However, significant heterogeneity were also detected in two models among studies with Asian population (F vs f: I 2 = 63.4% and P = 0.027; FF vs Ff + ff: I 2 = 65.7% and P = 0.020) ( Table 4). Galbraith plot analyses indicated that Meng et al. 36 might be the source of heterogeneity. With exclusion of this study, the pooled results remain significant (F vs F: OR = 1.64, 95% CI = 1.31-2.04, P = 0.000; FF vs Ff + ff: OR = 2.07, 95% CI = 1.50-2.86, P = 0.000), with no significant heterogeneity (F vs F: I 2 = 0% and P = 0.718; FF vs Ff + ff: I 2 = 0% and P = 0.940). Subgroup analyses by study quality suggested that this positive association only existed in pooled analyses of high-quality studies (F vs f: OR = 1.58, 95% CI = 1.10-2.26, P = 0.013; FF vs Ff + ff: OR = 1.92, 95% CI = 1.09-3.40, P = 0.025).  (Table 3). Similar results were also observed in the subgroup analyses by ethnicity (Table 4). Moreover, sensitivity analyses showed the results did not change meaningfully by excluding two studies 31,34 departed from HWE or one study with low-quality 36 . There was no significant heterogeneity for all models except the allele model (I 2 = 52.1% and P = 0.064). A Galbraith plot analysis suggested that Stefanic et al. 30 might be the source of heterogeneity for the allele model. Omitting this study, the pooled result was still not statistically significant (B vs b: OR = 1.06, 95% CI = 0.85-1.31, P = 0.615), with no significant heterogeneity (I 2 = 0% and P = 0.621).
ApaI polymorphism. Six studies including 766 cases and 813 controls evaluated the association of ApaI polymorphism and HT risk. The meta-analyses demonstrated no positive relationship of ApaI polymorphism and HT risk in the overall population (A vs a: OR = 0.98, 95% CI = 0.82-1.19, P = 0.869; AA vs aa: OR = 0.90, 95% CI = 0.60-1.36, P = 0.615; Aa vs aa: OR = 1.06, 95% CI = 0.82-1.36, P = 0.670; AA + Aa vs aa: OR = 1.01, 95% CI = 0.78-1.32, P = 0.916; AA vs Aa + aa: OR = 0.92, 95% CI = 0.65-1.29, P = 0.620). No significant heterogeneity was found in all the comparisons (all P > 0.05, Table 3). Similar results were found in the subgroup analyses by ethnicity; ApaI polymorphism was not associated with HT risk in Asian or Caucasian populations (Table 3). Sensitivity analyses, by excluding these two studies 33,35 not in HWE or one study with low-quality 36 , suggested that the results were consistent with those of the primary analyses (all P > 0.05).
TaqI polymorphism. A total of 902 cases and 863 controls from seven studies investigated the relationship between TaqI polymorphism and HT risk. The genotype distribution was consistent with HWE in the controls of all studies (all P > 0.05, Table 2). The pooled results showed that the TaqI polymorphism wasn't significantly asso-  Table 3). There was significant heterogeneity for comparison of T vs t and TT vs Tt + tt (I 2 = 70.8%, P = 0.002 and I 2 = 75.4%, P = 0.000, respectively). In the Galbraith plots, two studies 33,35 were outside of the 95%CI from the log OR, causing the heterogeneity in the results. When these two studies were excluded, the heterogeneity decreased significantly, but the pooled results were not changed significantly (T vs t: OR = 1.16, 95% CI = 0.95-1.41, P = 0.147; I 2 = 0% and P = 0.635 for heterogeneity; TT vs Tt + tt: OR = 1.16, 95% CI = 0.90-1.50, P = 0.262; I 2 = 0% and P = 0.788 for heterogeneity). Subgroup analyses by ethnicity found the similar results in Caucasian or in Asian (all P > 0.05) ( Table 4).

Publication bias.
No evidence of publication bias was detected by visual inspections of these funnel plots and Egger's test in all the models regarding the FokI, TaqI and ApaI polymorphism (all P Egger's > 0.05). However, significant publication bias was detected in two models regarding BsmI polymorphism (P Egger's = 0.001 for Bb vs bb and P Egger's = 0.005 for BB + Bb vs bb) ( Table 3, Fig. 4A and C). We used the trim and fill method incorporating the hypothetical studies to recalculate the pooled risk estimate. The pooled analyses continued to show no significant association between BsmI polymorphism and HT risk (Bb vs bb: OR = 0.90, 95% CI = 0.71-1.15, P = 0.397;

Study
Year Country Ethnicity Genotyping method

Discussion
To our knowledge, this is the first meta-analysis specially focused on the association of VDR polymorphism with HT risk. A significant association between the BsmI and TaqI polymorphisms and AITD risk has been reported by a previous meta-analysis 39 . However, in that study, the AITD, including GD and HT, was regarded as an entirety to analyze and only two studies 29,30 concentrated on HT alone among all the contained studies. Although GD and HT shared similar immune-mediated mechanisms characterized by the production of thyroid autoantibodies and by thyroid lymphocytic infiltration, a number of studies has indicated that the two diseases might harbor different susceptibility genes 5,34,40 . Thus, it is necessary to perform a meta-analysis specially focused on HT. Recently, several individual studies [33][34][35][36][37] have been conducted to investigate the association between the VDR gene polymorphisms and HT risk, but results from these studies remain conflictive and inconclusive. The reasons for this discrepancy may be small sample size, extensive geographic variations and difference in lifestyle and ethnicities. Therefore, in order to overcome the potential limitations of individual studies, we performed a meta-analysis and found that VDR FokI but not TaqI, ApaI and BsmI polymorphism was significantly associated with the risk of HT. Furthermore, the positive association of FokI polymorphism was only detected in Asians, not in Caucasians by subgroup analyses based on ethnicity.
Polymorphism FokI (rs2228570), located in the translational initiation site of VDR, which is the only known VDR gene polymorphism that results in the generation of an altered protein [41][42][43] . It can produce two structurally distinct isoforms: a shorter F-VDR or a longer f-VDR protein. The shorter F-VDR protein variant has been reported to be more active than the longer protein variant 44,45 . Transfection experiments showed the presence of short F-VDR resulted in a higher NF-kB-and NFAT-driven transcription capacity compared to the longer f-VDR. Concordantly, human monocytes and dendritic cells with a homozygous FF VDR genotype show higher expression of IL-12 (mRNA and protein) compared to the cells with an ff VDR genotype 46 . Therefore, individual with FF genotype may have a more active immune system and an increased risk to immune-mediated diseases.

Study
Year Genotype distribution in case Genotype distribution in control Eight previous studies investigated the distributional difference of FokI polymorphism in patients with HT and controls, and six found a positive association, but another two studies 36,38 did not. By pooling these results, our meta-analysis demonstrated that the F allele might be a risk factor for susceptibility of HT (OR = 1.44, P = 0.010) and the incidence of HT was significantly higher in FF genotype individuals than that of Ff + ff genotype individuals in overall population (OR = 1.72, P = 0.019). In addition, results from subgroup analyses stratified by ethnicity indicated that HT risk was increased in Asians with FF genotype (OR = 1.64, P = 0.039), but not in Caucasians. This inconsistent result in these two ethnicities may be due to the influence of different genetic backgrounds, lifestyle and environment factors (such as sunlight exposure and diet). In addition, an insufficient number of samples for analysis might lead to unreliable conclusions with deviation in Caucasians. BsmI (rs1544410), ApaI (rs7975232), and TaqI (rs731236) SNPs, located near the 3′ end of the VDR gene, are in strong linkage disequilibrium (LD) with each other. These three SNPs don't change the amino acid sequence of the encoded protein but have been shown to affect gene expression through regulation of mRNA stability 47 . Three studies 30,33,38 indicated TaqI polymorphism was associated with risk of HT in Croatian and Turkish population, but four other studies [34][35][36][37] from China, Japan, Italy and Serbia showed no association. ApaI polymorphism was reported no association with HT risk in previous studies with consistent results. Regarding BsmI polymorphism, the study conduct by Stefanic et al. 30 demonstrated B variant was apparently associated with decreased risk for HT in comparison to the reference b allele, but five other studies didn't find this association. In present meta-analysis, pooled results showed no significant association between HT disease and TaqI, ApaI or BsmI polymorphism. Furthermore, subgroup analyses found similar results, and sensitivity analyses did not change the orientation of pooled results. VDR 3′ -RFLP haplotypes have been positioned within the regulatory element spanning-3′ -untranslated region which contains polymorphic sequences affecting either VDR mRNA stability or VDR transcriptional activity 22,48 . Thus, BsmI, ApaI and TaqI, although functionally most likely anonymous, have been associated with total and allele-specific VDR mRNA expression 22 . Given these three variants strong LD with each other, it is rational to assess the haplotypes effects of VDR polymorphism on HT risk. Meng et al. 36 reported three common haplotypes (ab, Ab and AB) of ApaI-BsmI LD block were not associated with Chinese patients with HT (P > 0.05). Giovinazzo et al. 37 found the distribution of Bat and baT, the two most common BsmI-ApaI-TaqI haplotypes, was not significantly different in HT patients and controls from Italy. In another study 30 conducted in Croatia, the bT and BT of BsmI-TaqI haplotypes were found to be the predisposing and protective haplotypes, respectively. Similarly, common baT as well as the rare BaT haplotypes was associated with increased and decreased risk, respectively. However, we couldn't do meta-analysis due to insufficient published data in these studies. These effects, including effects associated with rare variants or specific stimuli need further research.
Vitamin D, well-known for its role in calcium and bone metabolism, has important effects on immune regulation by binding to the VDR localized in T lymphocytes and macrophages 49,50 . A number of studies 37,38,51-54 have found the serum vitamin D level was lower in subjects with HT than that of healthy controls. This inverse association indicated that vitamin D deficiency might be a causal factor leading to HT. Therefore, vitamin D level might be a significant confounder which should be considered when analyzing the association of VDR and HT risk.

SNPs
Sample size* (case/control) Genetic models  However, a different point of view has also been postulated, which suggested that the low level of serum vitamin D seen in disease is a secondary phenomenon of VDR dysfunction rather than the reason for autoimmunity 55 .
Although vitamin D level is seen as playing an important role, it is VDR dysfunction that is proposed to be the key factor in the autoimmune diseases process 56 . Because VDR is key to innate immune response which is important in the pathogenesis of autoimmune diseases 57,58 , VDR dysregulation greatly compromises the innate immune response. The 25-hydroxyvitamin D3 (25-OHD) level is a reliable parameter reflecting the vitamin D level of the body and usually measured as the level of vitamin D. When VDR dysregulation, the expression of CYP24A1, an enzyme that inactivating 1,25-dihydroxyvitamin D (1,25-OHD) was inhibited. Increased 1,25-OHD will decrease 25-OHD by reducing gene expression and inhibiting expression of CYP27A1 which is an enzyme involved in conversion of vitamin D into 25-OHD 55,59 . Among our included studies, only two studies concurrently provided the information on vitamin D levels and VDR in patients with HT. One study 37 found that the prevalence of vitamin D deficiency in HT patients was significantly higher than that in the control group (70% vs 18.2%; P = 0.0001), but VDR BsmI, ApaI, and TaqI polymorphisms were not associated with HT risk. The other 38 indicated that the prevalence of vitamin D insufficiency in HT cases was significantly higher than controls (P = 0.02) while VDR TaqI, but not FokI polymorphisms is associated with HT. It is unfortunate that neither study analyzed the distributional difference of VDR polymorphisms stratified by vitamin D levels. Therefore, the mechanism and effect for the interaction of vitamin D and VDR in patients with HT need further investigations. Several limitations should be discussed when explaining the results of our meta-analysis. First, lack of adjustments for some factors, such as age, gender, thyroid functional status, circulating vitamin D levels, or dietary vitamin D intake, which may influence the association between VDR variants and risk of HT, might bias the present results. Second, because of unpublished data or limited number of studies, significant publication bias was found in two models regarding BsmI polymorphism, which might have some impact on the final outcome. However, we used trim and fill method to assess the influence of publication bias and found that the results were not significantly changed with or without the addition of hypothetical missing studies. Heterogeneity among studies was also detected in some analyses due to ethnic difference, geographic characteristics and lifestyle. However, our sensitivity analysis showed that studies that contribute to heterogeneity did not significantly alter the conclusions of the overall OR. Third, the statistical power to detect the association may be lower because number of studies included in our meta-analysis is relatively small. However, Ioannidis et al. 60 estimated the median sample size required to detect the observed summary effects in each population addressed in 752 studies is 3,535, which is 13.3-fold more subjects than in each original study. These sample size requirements can be inflated considerably if trying to account for potential bias or heterogeneity. These estimates may be difficult to address even by very large biobanks and observational cohorts. Therefore, meta-analysis is an effective way to explore the truth before the emergence of large sample data. Further studies should be focusing on innovative study designs and strong collaborative efforts.
In conclusion, our meta-analysis suggests that the VDR FokI polymorphism is associated with HT risk in overall population or in Asians, but not in Caucasians. The TaqI, ApaI and BsmI polymorphisms are not associated with HT risk. Further well-designed studies with larger sample sizes and different ethnic population are needed to clarify the present findings. Furthermore, the exact causality and mechanism for the interaction of VDR and HT development need further experimental or animal mechanism studies. The following key words and search terms were used to identify relative publications: "Vitamin D receptor", "VDR", "ApaI", "BsmI", "FokI", "TaqI" and "hashimoto's thyroiditis". The reference lists of identified articles and related reviews were reviewed for additional studies.
Inclusion and exclusion criteria. Studies meeting all of the following inclusion criteria were included: (1) case-control study or cohort study; (2) investigating the association between VDR gene polymorphisms (ApaI, BsmI, FokI and TaqI) and HT risk; and (3) providing the frequencies of the variants in cases and controls or providing sufficient data to calculate the estimation of odds ratios (ORs) with 95% confidence interval (95% CI). Exclusion criteria were as follows: (1) overlapping data; (2) studies without genotype frequency and genotype distribution or insufficient information for data extraction; (3) family-based study design; and (4) abstracts, reviews,  comments or editorial articles lack of necessary raw data. In the case of overlapping data, only the study with the largest population was selected for this meta-analysis.
Data extraction. Two investigators (XF Wang and WL Cheng) extracted data independently. Any disagreement was resolved through discussion. The extracted data included: name of the first author, year of publication, country, ethnicity, number of cases and controls, genotyping method, control sources, and genotype distribution in cases and controls.
Quality Assessment. The quality of included studies was assessed by two independent reviewers (XF Wang and Y Ma) using the Newcastle-Ottawa Scale (NOS) 61 . The NOS judged a study based on three perspectives: selection, comparability and exposure/outcome. The full score was 9 stars. Study that scored above six stars was considered as high quality.

Statistical analysis.
A random-effects model was used to incorporate within-and between-study heterogeneity as this can provide more conservative result than a fixed effects model 62 . Pooled ORs and their respective 95% CIs were calculated to evaluate the association between the four VDR SNPs and HT risk under five genetic models: the allele model (eg, A vs a), the homozygous model (eg, AA vs aa), the heterozygous model (eg, Aa vs aa), the recessive model (eg, AA + Aa vs aa), and the dominant model (eg, AA vs Aa + aa). The Hardy-Weinberg equilibrium (HWE) in controls was tested using the goodness-of-fit χ 2 statistic with one degree of freedom 63 .
Cochrane's Q test and I 2 test were used to assess heterogeneity among trials. Q-test reported a P value < 0.1 or I 2 > 50% was defined as significant heterogeneity 64 . In case of substantial heterogeneity, a Galbraith plot was created to graphically identify the potential outlier studies that might cause the heterogeneity. Then, a meta-analysis was rerun after excluding the outlier studies 65 . Subgroup analyses were performed based on ethnicity and quality of included studies to avoid the potential bias influence. Sensitivity analyses were performed by excluding each individual study or the studies with controls inconsistent with HWE to evaluate the impact of individual study on the pooled risk estimate. Publication bias was evaluated by a visual inspection of funnel plot and Egger's test 66 . If publication bias was indicated, the "trim and fill" method which conservatively imputes hypothetical negative unpublished studies to mirror the positive studies that cause funnel plot asymmetry was performed to further assess the possible effect of publication bias 67 . All P-values were two-tailed. All analyses were performed using Stata 11.0 (Stata Corporation, College Station, TX, USA). This article follows the PRISMA statement 68 and the Cochrane Collaboration guidelines for reporting meta-analysis.