**1. Introduction**

Type 1 diabetes (T1D) is a chronic autoimmune disease, resulting from autoimmune degradation of pancreatic ß-cells leading to the inability to produce and/or use insulin [1]. T1D patients carry a genetic susceptibility to autoimmune disease development, with first-degree relatives of those affected also carrying an increased risk of developing the disease [2,3]. Undiagnosed or untreated T1D can result in hyperglycaemia, increasing the risk of developing microvascular and macrovascular injuries/health complications, such as nephropathy, ischemic heart disease and stroke [4]. Estimates of those with T1D below age 20 had risen to over a million in 2017, with evidence of increasing incidence worldwide [5]. Presently, there are no established treatments identified for the prevention of T1D and the search for genetic and environmental triggers remains ongoing.

Emerging evidence suggests low vitamin D status may play a role in T1Dpredisposition. Vitamin D is a steroid prohormone, with nutrition status approximated via serum 25- hydroxyvitamin D [25(OH)D] concentrations [6]. Notably, 25(OH)D deficiency is strongly associated with skeletal pathology, however, in the advent of vitamin D receptors being discovered throughout the body, there now is a greater acknowledgment of broader disorders associated with deficiency, including autoimmune issues, such as T1D and multiple sclerosis [7,8]. Recent evidence indicates an important role for active vitamin D [1,25(OH)2D] in

**Citation:** Najjar, L.; Sutherland, J.; Zhou, A.; Hyppönen, E. Vitamin D and Type 1 Diabetes Risk: A Systematic Review and Meta-Analysis of Genetic Evidence. *Nutrients* **2021**, *13*, 4260. https:// doi.org/10.3390/nu13124260

Academic Editors: Daniel-Antonio de Luis Roman and Ana B. Crujeiras

Received: 27 October 2021 Accepted: 25 November 2021 Published: 26 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

immune regulation [9]. Mechanistic explanations for 1,25(OH)2D include immunomodulatory action leading to cytokine regulation, reducing the likelihood of destruction of pancreatic ß-cells [10]. Another potential mechanism is through direct protection of pancreatic ß-cells, serving to preserve barrier exclusion of pathogens, likely significant in the prevention of autoimmune disorders [11]. Such mechanistic insight has underpinned novel immune-modulatory concepts for the prevention of T1D.

Association between serum 25(OH)D concentrations and T1D risk is supported by evidence from in vitro and animal experiments [12–14], as well as human observational studies [15–18] and ecological correlation [19]. In animal studies, oral administration of the activated form of vitamin D was found to protect nonobese diabetic mice from T1D [12–14], while human observational studies have shown reduced levels of serum 25(OH)D are associated with increased risk of T1D [15,17]. In the aetiology of T1D observational studies have also shown support of vitamin D supplementation in being inversely associated with T1D [16,18,20]. Animal experimental data, therefore, indicate low 25(OH)D concentrations may be involved in T1D predisposition, however, a causal role of impaired vitamin D metabolism in the aetiology of T1D in humans is ye<sup>t</sup> to be implicated, and stronger forms of evidence—less effected by confounding or reverse causation—are required.

Using selected vitamin D related genetic variants, it is possible in a genetic epidemiological setting to establish evidence of an etiological role of 25(OH)D in T1D pathophysiology. Since 25(OH)D synthesis is regulated by genes, single nucleotide polymorphisms (SNPs) may alter the bioavailability and target effects of vitamin D metabolites. Large-scale genome-wide association studies (GWAS) have identified several SNPs from genes influencing 25(OH)D levels; *CYP2R1*, *DHCR7/NADSY1*, *GC*, *CYP24A1*, *AMDHD1* and *SEC23A*, which have been used as genetic instrumental variants in this study [21,22].

As individual studies may not have enough statistical power to identify an association between selected genetic variants affecting serum 25(OH)D concentrations and T1D, a metaanalysis is a useful statistical tool to pool data from published studies, where increasing the statistical power can give more accurate estimates of effect sizes. In this study, we perform a systematic review and meta-analysis of all existing studies reporting an association between selected 25(OH)D related genetic variants (exposure) and T1D risk (outcome) in humans (population). This topic provides a further scientific understanding of T1D pathophysiology and the potentiality of preventing T1D through increases in 25(OH)D concentrations.

### **2. Materials and Methods**

This systematic review and meta-analysis followed the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines [23]. Registration: PROS-PERO (ID CRD42021224844), https://www.crd.york.ac.uk/prospero/ (accessed on 10 January 2021).

### *2.1. Search Strategy*

A search was conducted in four databases: Ovid Medline (1964-present), Ovid Embase (1947-present), Web of Science (1975-present), IEU OpenGWAS (2020-present) from inception to April 2021. The primary search terms were as follows: humans, single nucleotide polymorphism, genetic variation, type 1 diabetes mellitus and vitamin D. The selection of articles in Medline and Web of Science was performed using Medical Subject Headings (MeSH) to define these descriptors. The selection of articles in Embase was performed using Emtree (Embase subject headings) to define these descriptors. Boolean operators (e.g., OR, AND, NOT) were also combined with keywords and subject headings. An initial pilot search was undertaken to improve inclusion clarity of study inclusion and exclusion, improving accuracy and consistency. The strategy was developed by one reviewer (L.N.) and proofread for syntax, spelling and overall structure by two reviewers (E.H. and J.S.). As part of the development process, we used two relevant, existing studies [24,25] for validation purposes, testing if our search strategy could identify them. The set of search terms

was slightly modified between databases due to different system procedural limitations, however, the overall approach remained as consistent as possible across each database. The selection of studies through OpenGWAS, as well as the UK Biobank, was prepared using R 4.0.2 software, conducting an SNP-based search for the selected genetic variants and their proxies (r2 > 0.8), locating any additional studies fitting the inclusion criteria. Full search strategies are presented in Supplementary Tables S1–S4.

### *2.2. Inclusion and Exclusion Criteria*

Studies testing exposure of selected genetic variants or their proxies with r2 > 0.8 influencing 25(OH)D pathways for association with T1D status and 25(OH)D concentrations, were of interest. Eligible studies met the population, exposure, outcome (PEO) approach [26] as follows:


The following exclusion criteria were also used:


No language, publication status, or publication date limitations were imposed.

### *2.3. Study Selection and Data Extraction*

Literature was searched in duplicate independently by two authors (L.N. and J.S.), and approved by a third author (E.H.). After excluding duplicates, article selection was carried out in two passes. In the first pass, title and abstract screening occurred for the selection of relevant papers meeting the eligibility criteria. In the second pass, proposed articles from the first pass were screened in full text for compliance with inclusion criteria. To ensure literature saturation, reference lists of obtained studies from original database searches were manually scanned for potential unidentified additional studies by one author

(L.N.), with eligibility confirmed by a second author (J.S.). Furthermore, OpenGWAS was used to identify unpublished studies, locating one FinnGen cohort study sharing summarylevel data fitting the search parameters. Datasets were also identified in the UK Biobank, a large-scale prospective cohort study.

Data were extracted independently by two authors (L.N. and J.S.) using a predetermined data extraction template. The following data were extracted from the articles included in this systematic review: first author; region/demographic information; publication year; study design characteristics; participant characteristics, including gender and ethnicity if reported; the number of cases and controls studied; mean age (or range) at the onset of T1D in cases; outcome measure, diagnostic criteria of T1D; mean age (or range) of the control group; how the controls were selected; genotyping methods, genotype distribution, and allele frequency in cases and controls; all reported patient outcome measures; key findings; protocol availability and funding sources. Corresponding authors were contacted by e-mail for missing or unreported data a maximum of three attempts, to avoid any assumptions made from unclear information. All disagreements were resolved by consensus, or with the input of a third author (E.H. or A.Z.).

### *2.4. Statistical Analysis*

All mentioned statistical analyses were performed with STATA 16.0 software (Stata Corporation, College Station, TX, USA) and R 4.0.2 software by two authors (L.N. and A.Z.). For each variant, OR per vitamin D-increasing allele was extracted from individual studies for the meta-analysis, as per the SUNLIGHT consortium [21]. If a study did not contain the selected vitamin D variant, the result of its proxy (r2 > 0.8) was extracted and used to estimate the related effect. In studies where the OR per vitamin-D-increasing allele was not reported, we estimated the allelic effect from the contingency table of T1D distribution by SNP genotypes, where OR was computed by dividing the odds of T1D in the heterozygotes (i.e., with 1 25-hydroxyvitamin D increasing allele) by that in the homozygotes (i.e., with 0 25-hydroxyvitamin D increasing allele). Meta-analysis was performed using the randomeffects model (REM, restricted maximum likelihood method) [28]. Heterogeneity between studies was assessed using Cochran's Q test the I2 statistic, with heterogeneity considered to be substantial if the *p*-value for the Cochran's Q test < 0.05 or I2 > 50%. All *p*-values were for two-tailed tests, and <0.05 was considered statistically significant.

We conducted sensitivity analyses by removing a single study at a time, evaluating the integrity of the results. Subgroup analysis was performed by restricting the sample to the Caucasian population, to examine the possible effects of population stratification. Initial protocol pre-specified plan for further MR analyses, which were not conducted as it was considered redundant given clear results.

#### *2.5. Risk of Bias and Credibility of the Evidence Assessment*

The methodological quality of eligible studies was evaluated using Critical Appraisal Skills Program tools for cohort and case-control studies [29]. Two authors (L.N. and J.S.) independently completed risk of bias assessment and recorded supporting justification and information for each domain to optimise the tool's value (met; partially met; not met; unclear). The domains were: Are the results of the study valid? Were the cases recruited in an acceptable way? Was the exposure accurately measured to minimise bias? Have authors taken account of the potential confounding factors in the design and/or in their analysis? How precise are the results? (size of confidence intervals). Results were compared by categorising each study for study quality (risk of bias) judgement (low, some concerns, high). Articles were judged as 'low' when five or more domains were met. Conversely, articles were judged as 'high' when three or more domains were unmet. Disagreements were resolved by a third author (E.H.). Outcome reporting bias was assessed by comparing outcomes specified in protocols, with outcomes reported in corresponding publications. Where protocols were not available, outcomes specified in the methods and results sections of publications were compared.

Two reviewers assessed the risk of bias due to missing results in a synthesis (L.N. and A.Z.). Potential publication bias was assessed by examining for asymmetry using Begg's funnel plot for each SNP [30]. If publication bias was present, the plot would be asymmetric, indicating a deficiency in publications with negative results. No further formal assessment of publication bias, such as Egger's test was performed, due to insufficient studies [31].
