**1. Introduction**

Psychiatric disorders are a group of complex diseases, which are mainly characterized by varying degrees of obstacles in mental activities such as cognition, emotion, willpower, and behavior [**?** ]. A meta-analysis showed that overall global prevalence of psychiatric disorders was increased with odds ratio of 1.179 (95% CI: 1.065–1.305) [**?** ]. Additionally, another study has found that the prevalence of the common mental illnesses is continuously rising, particularly in low- and middle-income countries, with many people suffering from depression and anxiety disorders simultaneously [**?** ]. Depression and anxiety disorders are the most common mental disorders in the general population [**?** ]. According to the latest report of World Health Organization, it is estimated that there are 322 million people (4.4% of the world's population) living with depression and more than 260 million people (3.6% of the global population) affected by anxiety disorders [**?** ]. Furthermore, Wang et al. observed that depression and anxiety were significantly associated with higher cancer incidence, cancer-specific mortality and all-cause mortality [**?** ]. Traditional

**Citation:** Zhang, Z.; Yang, X.; Jia, Y.; Wen, Y.; Cheng, S.; Meng, P.; Li, C.; Zhang, H.; Pan, C.; Zhang, J.; et al. Vitamin D and the Risks of Depression and Anxiety: An Observational Analysis and Genome-Wide Environment Interaction Study. *Nutrients* **2021**, *13*, 3343. https://doi.org/10.3390/ nu13103343

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

Received: 25 August 2021 Accepted: 21 September 2021 Published: 24 September 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/).

treatment of depression and anxiety can only relieve rather than cure the condition, causing a tremendous economic burden on individual families and high suicide rate. Therefore, it is urgen<sup>t</sup> to find new ways to treat and prevent depression and anxiety disorders.

Vitamin D (VD) is a member of the steroid hormone family. As a necessary vitamin to human body, it is not only cheap, but also widely available. Besides being ingested from daily diet and VD supplements, it can also be synthesized from 7-dehydrocholesterol in the skin through ultraviolet radiation-b (UVB) [**?** ]. Previous studies on VD have found that genetic and environmental factors can affect the status and metabolism of VD. The data from twin and family studies have suggested that circulating VD concentrations are partly determined by genetics, with heritability ranging from 23% to 80% [**?** ]. Furthermore, genetic studies of VD have found that genetic variations and alterations (e.g., deletions, amplifications, inversions) in genes involved in VD metabolism, catabolism, transport, or binding to its receptors may affect VD levels [**?** ]. Moreover, a study of gene-environment interactions with VD has showed that specific genetic variants associated with VD metabolism may be correlated with prostate cancer risk in VD-deficient patients [**?** ]. As for the biological function of VD, previous studies have demonstrated that VD plays an important role in bone health, reproduction and fertility, and immune function [**?** ]. Growing evidence shows that VD also exerts a grea<sup>t</sup> influence on the development and adult brain, such as maintaining calcium balance and signaling, regulating neurotrophic factors, providing neuroprotection, modulating neurotransmission, and promoting synaptic plasticity [**?** ]. In addition, a recent systematic review has shown that VD deficiency in adulthood may also be associated with adverse brain-related outcomes [**?** ]. The discovery of the functions of VD in the brain and the continuous confirmation of the effects of VD on disease provides a new research direction for the study of psychiatric disorders in recent years.

Although a lot of research has been devoted to exploring the relationship between VD and depression and anxiety, it is still controversial whether VD was associated with the two mental disorders [**?** ]. Some studies have observed an association between VD and depression symptoms [**???** ], while other studies have not found this association [**? ?** ]. Additionally, previous studies on the association between VD and anxiety symptoms observed inconsistent association [**? ?** ]. The discrepant results may be not only related to potential confounding factors, but also to complex etiology of psychiatric disorders. Previous studies on the pathogenesis of psychiatric disorders have found that they are associated with a variety of factors. For example, numerous studies have proved that inflammation is related to the pathogenesis of psychiatric disorders, and that the increase in inflammation will affect the occurrence and development of some psychiatric disorders [**? ?** ]. Other studies have found that both diet and the gu<sup>t</sup> microbiome have a strong influence on emotional behavior and neural processes [**?** ]. Dietary patterns and changes in the microbiome can affect symptoms of depression and anxiety disorder and increase the risk of both disease through the microbiome gut–brain axis (MGBA) [**?** ]. Moreover, some studies have shown that psychiatric disorders are caused by the interplay between genetic susceptibility and environmental risk factors [**? ?** ]. It has been reported that the heritability of major depression [**?** ] and current anxiety symptoms [**?** ] was estimated to be 38% and 31%, respectively. However, previous studies on exploring the associations between VD and depression and anxiety only focused on the influence of environmental or genetic factors, usually not considering the interactions between them, which may underestimate the effects of VD on the risks of depression and anxiety.

Polygenic risk score (PRS) is proposed by running a GWAS on a study sample, selecting SNPs according to the relevant phenotypes, and creating the sum of phenotype related alleles (usually weighted by the SNP-specific coefficients from the GWAS) that can be evaluated in other replication sample [**?** ]. PRS analysis can not only explore the genetic associations between various complex diseases and traits, but also assess the influence of susceptible loci on disease risks [**?** ]. For example, PRS-related studies conducted by Psychiatric Genomics Consortium have found significant associations between PRS of symptom scale score and the risk of schizophrenia [**?** ] as well as the efficacy of antipsychotic medication [**?** ]. Recently, Revez et al. [**?** ] conducted a GWAS of 417,580 UK Biobank study participants and identified 143 genetic loci associated with VD. Using the PRS of VD traits as instrumental variables, we can further explore the associations between the PRS of VD traits and psychiatric disorders.

GWAS study has strong ability in identifying susceptibility genetic loci associated with psychiatric disorders [**?** ]. However, previous studies have shown that most of the phenotypic variation in complex diseases and traits cannot be explained solely by genetic factors, because phenotypic variation can also occur through genetic environment (G × E) interactions, in which the genotypes of different individuals vary in response to environmental stimuli [**?** ]. Therefore, we further adopted the genome-wide environment interaction study (GWEIS), which can not only assess the effects of G × E interactions, but also evaluate the effects of genetic interactions on a genome-wide scale, helping to identify new genetic risk variants and understand the potential biological mechanisms [**?** ]. For example, Rivera et al. found 53 and 34 single nucleotide polymorphisms (SNPS) in additive interactions with smoking in Lofgren's syndrome (LS) and non-LS, respectively, but the association did not persist when assessing the effect of smoking on sarcoidosis without genetic information [**?** ].

Utilizing the individual blood VD levels and calculated VD PRS data in UK Biobank cohort, we conducted regression analysis to evaluate the associations between VD (including blood VD levels and calculated VD PRS data) and depression status and anxiety status in the UK Biobank. Then, based on the results of regression analysis, genome-wide environment interaction studies (GWEIS) were performed to clarify the potential effects of gene × Vitamin D interactions on depression and anxiety.

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

### *2.1. UK Biobank Cohort*

The UK Biobank (UKB) is a large population-based prospective cohort study, with health-linked information, both regarding phenotype and genotype, on approximately 500,000 participants aged between 40 and 69 years from all over the United Kingdom in 2006 and 2010 [**?** ]. All participants were asked to report a series of health status and demographic information through questionnaires and interviews, and approved to use their anonymous data for any health-related research. Informed consent was provided by UKB from all participants. This study was approved by UKB (Application 46478) and obtained participants' health-related records.

### *2.2. UK Biobank Phenotypes of Vitamin D*

A total of 376,803 UKB participants' blood samples were collected for quantitative measurement of 25(OH)D levels via chemiluminescent immunoassay (CLIA), and 343,334 (91.12%) individuals had their vitamin D 25(OH)D levels (UK Biobank data field: 30,890) measured. The analysis was limited to the population of white British individuals (UK Biobank data field: 21,000).

### *2.3. UK Biobank Phenotypes of Depression and Anxiety*

The phenotypes of depression and anxiety were defined according to the previous study [**?** ]. The selection criteria of case group were defined based on self-reports (UK Biobank data fields: 20,002; 20,126; and 20,544). In order to classify participants as much as possible, patient health questionnaire-9 (PHQ-9) [**?** ], general anxiety disorder-7 (GAD-7) [**?** ], and composite international diagnostic interview short-form (CIDI-SF) [**?** ] were used as strict inclusion and exclusion criterion. PHQ-9 and GAD-7 are score scales for depression and anxiety, used to screen and measure the severity of depression and anxiety, respectively. PHQ-9 is a classification scale focusing on nine depression symptoms and signs, with a total score (0–27) [**?** ], and GAD-7 is a classification scale focusing on seven anxious symptoms and signs, with a total score (0–21) [**?** ]. Detailed classification of depression and anxiety are presented in the Supplementary Information. The selection

criteria of the control group were as follows: without symptoms of depression and anxiety defined by CIDI and self -reported, depression PHQ scores and anxiety GAD scores ≤ 5, and without core symptoms.

#### *2.4. UK Biobank Genotyping, Imputation and Quality Control*

A total of 488,377 participants of UKB cohort were genotyped by either the Affymetrix UK BiLEVE Axiom Array or the Affymetrix UKB Axiom arrays [**?** ]. The imputation and the quality control of these genotype results were carried out based on UK10K project reference panel [**?** ] and Haplotype Reference Consortium (HRC) [**?** ] reference panel. Then, we removed the participants who reported inconsistencies between self-reported gender and genetic gender, without ethic consents and imputation data. Additionally, we excluded variants with the Hardy–Weinberg equilibrium test *p* > 1.0 × <sup>10</sup>−5, a minor allele frequency (MAF) of < 0.01, and a genotype missing rate of > 0.05. Ultimately, we used KING software to exclude the genetically related individuals.

### *2.5. GWAS Data of Vitamin D*

The latest large-scale GWAS summary statistics of VD were used here. Briefly, this GWAS dataset detected 18,864 independent SNPs that were statistically associated with VD [**?** ]. The genotype data were quality-controlled and imputed against the HRC and UK10K by the UKB group. Then, a linear mixed model GWAS was implemented in fastGWA to identify the genetic loci associated with 25OHD concentrations. Additionally, a rank-based reverse normal transformation (RINT) was applied to the phenotype, age, and gender. The genotyping batch and the first 40 ancestry PCs were used as covariates in the mixed model. In order to determine the independent associations, a conditional and joint (COJO) analysis [**?** ] was employed on the GWAS results to explain the correlation structure between SNPs in the 10-Mb window (COJO default parameters). Detailed information of the GWAS can be obtained in the published study [**?** ].

### *2.6. PRS Analysis of Vitamin D*

Using the VD associated SNPs from the GWAS (*p* < 5 × <sup>10</sup>−8) [**?** ], the PRS of VD of each individual was calculated as the sum of the risk allele they carried, weighted by the effect size of the risk allele [**?** ]. The PRS of VD was computed by PLINK2.0 [**?** ], according to the formula:

$$\text{PRS} = \sum\_{i=1}^{n} \beta\_i \text{SNIP}\_{\text{im}}$$

PRS denotes the PRS value of VD for UKB subject; n and i, respectively, denote the total number of sample size and genetic markers; βi is the effect parameter of risk allele of the significant SNP associated with VD, which was obtained from the GWAS of VD [**?** ]; and SNPim is the dosage (0, 1, 2) of the risk allele of the SNP associated with VD [**?** ].

### *2.7. Statistical Analysis*

In the UK Biobank cohort, we evaluated associations between vitamin D and depression and anxiety through regression analysis. Specifically, logistic regression analysis was employed to evaluate the associations of self-reported depression and anxiety status with blood VD, VD PRS before COJO adjustment, and VD PRS after COJO adjustment. Linear regression analysis was conducted to test the associations of the PHQ-9 score and the GAD-7 score with blood VD, VD PRS before COJO adjustment, and VD PRS after COJO adjustment. In this regression analysis, blood VD, VD PRS before COJO adjustment, and VD PRS after COJO adjustment were used as independent variables. Self-reported depression, self-reported anxiety, PHQ-9 score, and GAD-7 score were used as outcome variables. Sex, age, and 10 principle components (PCs) of population structure were used as covariates in the regression analysis. A *p* < 0.05 indicated an association. All analyses were conducted by R.

### *2.8. Genome-Wide Environment Interaction Studies (GWEIS)*

The generalized linear regression model of PLINK2.0 [**?** ] was used to estimate the gene × VD interaction effects on the risk of depression and anxiety, using age, gender, and the first 10 PCs as covariates. According to previous research [**?** ], we used PLINK2.0 genetic additive (ADD) models and selected high-quality SNPs through a quality control filters: SNPs with a low call rates (<0.90), low minor allele frequencies (<0.01), or low Hardy–Weinberg equilibrium exact test *p*-values (<0.01) were excluded. *p* < 5.0 × 10−<sup>8</sup> and *p* < 5.0 × 10−<sup>7</sup> were defined as significant and suggestive interactions, respectively. GWEIS results were visualized with the circular Manhattan plots generated by the "CMplot" R script. (https://github.com/YinLiLin/R-CMplot) (accessed on 15 February 2021).

### **3. Results**

### *3.1. General Population Characteristics*

#### 3.1.1. Characteristics of UK Biobank Subjects with Blood Vitamin D Data

For depression traits, with self-reported depression status as the outcome variable, a total of 110,744 participants answered the depression-related questions, and 52,766 were classified into depression group. With the depression PHQ score as the outcome variable, a total of 109,543 participants completed the questionnaire. For anxiety traits, with self-reported anxiety status as the outcome variable, a total of 98,784 participants answered the anxietyrelated questions, and 19,759 were classified into anxiety group. With the anxiety GAD score as the outcome variable, a total of 110,023 participants completed the questionnaire.

#### 3.1.2. Characteristics of UK Biobank Subjects with Vitamin D PRS Data

In the self-reported depression status, a total of 121,685 participants answered depressionrelated questions, of which 58,349 were included in depression group; in the depression PHQ scores, a total of 120,033 participants completed the questionnaire. In the self-reported anxiety status, a total of 108,309 participants answered anxiety-related questions, of which 21,807 were classified into anxiety group. In addition, in the anxiety GAD scores, a total of 120,590 participants completed the questionnaire. Detailed information is shown in Table **??**.


**Table 1.** General population characteristics of this study participants from UK Biobank.

Abbreviations: VD, Vitamin D; VDPRS, polygenic risk score of vitamin D; COJO, conditional and joint analysis; SD, age was described as mean ± standard deviation (SD); PHQ score, patient health questionnaire (PHQ) is used to describe the depression; GAD score, general anxiety disorder (GAD) is used to describe the anxiety of the participants.

### *3.2. Regression Analysis Result*

3.2.1. Associations between Blood Vitamin D and Depression, Anxiety Traits in UK Biobank Cohort

Significant associations of blood VD level with self-reported depression status (odds ratio (OR) = 0.89, *p* = 5.92 × <sup>10</sup>−77) and self-reported anxiety status (OR = 0.92, *p* = 1.46 × <sup>10</sup>−22) were observed. Associations were also observed between the blood VD level, the depression PHQ score (Beta = −0.062, standard error (SE) = 0.003, *p* = 5.95 × <sup>10</sup>−96), and the anxiety GAD score (Beta = −0.030, SE = 0.00, *p* = 1.21 × <sup>10</sup>−21).

3.2.2. Associations between Vitamin D PRS and Depression, Anxiety Traits in UK Biobank Cohort

We observed significant associations of VD PRS before COJO adjustment with selfreported depression status (OR = 0.99, *p* = 3.82 × <sup>10</sup>−2), depression PHQ score (Beta = −0.0060, SE = 0.003, *p* = 3.25 × <sup>10</sup>−2), and anxiety GAD score (Beta = −0.010, SE = 0.00, *p* = 4.36 × <sup>10</sup>−2). In addition, we also observed significant associations of VD PRS after COJO adjustment with a self-reported depression status (OR = 0.99, *p* = 1.84 × <sup>10</sup>−2), a depression PHQ score (Beta = −0.0070, SE = 0.0030, *p* = 9.15 × <sup>10</sup>−3), and an anxiety GAD score (Beta = −0.010, SE = 0.00, *p* = 1.02 × <sup>10</sup>−2). Detailed information is shown in Table **??**.

**Table 2.** The associations between Vitamin D traits and traits of depression and anxiety.


Abbreviations: SE, standard error; T, *t*−test; OR, odd ratios; VD, vitamin D; VDPRS, polygenic risk score of vitamin D; COJO, conditional and joint analysis; PHQ score, patient health questionnaire (PHQ) is used to describe the depression; GAD score, general anxiety disorder (GAD) is used to describe the anxiety of the participants. Note. Logistic regression was used to evaluate the association between blood VD, VD PRS before COJO adjustment, VD PRS after COJO adjustment and self-reported depression and anxiety. Linear regression was used to evaluate the association between blood VD, VD PRS before COJO adjustment, VD PRS after COJO adjustment and PHQ score, GAD score.

### *3.3. GWEIS Analysis Results*

For self-reported depression status, GWEIS identified a significant gene × VD PRS interaction (*p* < 5.0 × <sup>10</sup>−8) at the LRRTM4 gene (rs114086183, *p* = 4.11 × <sup>10</sup>−8). For self-reported anxiety status, significant gene × VD PRS interaction was detected at the GNB5 gene (rs149760119, *p* = 3.88 × <sup>10</sup>−8). For the depression PHQ score, two significant gene × blood VD interactions were identified at SLC11A2 and HIGD1C (rs117102029, *p* = 4.02 × <sup>10</sup>−8). For the anxiety GAD score, we detected multiple significant gene × VD trait interactions, such as SMYD3 (rs142593645, *p* = 2.51 × <sup>10</sup>−8), SEMA3E (rs76440131, *p* = 2.80 × <sup>10</sup>−10), and VTI1A (rs17266687, *p* = 3.09 × <sup>10</sup>−8). Among them, three genes (SEMA3E, DOCK8, TMCO3) were identified by VD PRS before and after COJO adjustment. The visualization of the results is shown in **??????**. Additional detailed results are shown in Table **??**.

**Figure 1.** Genomic regions interacting with blood VD for the PHQ score and the GAD score. (**a**) Depression (PHQ score), A SNP allele was found to significantly interact with blood VD in depression (PHQ score); (**b**) Anxiety (GAD score), seven independent SNP alleles were found to significantly interact with blood VD in anxiety disorder (GAD score). From the center, the first circos depicts the −log10 *p*-values of each variant due to double exposure, i.e., the effect of both SNP allele and blood VD. The second circos shows chromosome density. Red dots represent the *p* < 5 × 10−<sup>8</sup> and green dots represent *p* < 1 × 10−7. The figure was generated using the "CMplot" R script (https://github.com/YinLiLin/R-CMplot) (accessed on 15 February 2021).

(**c**) 

**Figure 2.** Genomic regions interacting with VD PRS after COJO adjustment for depression status, anxiety status, and GAD score. (**a**) Depression status, a SNP allele was found to significantly interact with blood VD in depression status; (**b**) Anxiety status, 2 independent SNP alleles were found to significantly interact with blood VD in anxiety status; (**c**) Anxiety (GAD score), 8 independent SNP alleles interacted significantly with blood VD in anxiety disorders (GAD score). From the center, the first circos depicts the −log10 *p*-values of each variant due to double exposure, i.e., the effect of both SNP allele and VD PRS after COJO adjustment. The second circos shows chromosome density. Red dots represent the *p* < 5 × 10−<sup>8</sup> and green dots represent *p* < 1 × 10−7. The figure was generated using the "CMplot" R script (https://github.com/YinLiLin/R-CMplot) (accessed on 15 February 2021).

**Figure 3.** Genomic regions interacting with VD PRS before COJO adjustment for GAD score. This graph shows 7 independent SNP alleles interacting significantly with blood VD in anxiety disorders (GAD score). From the center, the first circos depicts the −log10 *p*-values of each variant due to double exposure, i.e., the effect of both SNP allele and VD PRS before COJO adjustment. The second circos shows chromosome density. Red dots represent the *p* < 5 × 10−<sup>8</sup> and green dots represent *p* < 1 × 10−7. The figure was generated using the "CMplot" R script (https://github.com/YinLiLin/R-CMplot) (accessed on 15 February 2021).

**Table 3.** Summary of gene −environment interaction analysis among SNP and VD traits for depression and anxiety traits.


Abbreviations: CHR, chromosome; SNP, single nucleotide polymorphism; SE, standard error; ADD, additive effect; VD, vitamin D; VDPRS, polygenic risk score of vitamin D; COJO, conditional and joint analysis; PHQ score, patient health questionnaire (PHQ) is used to describe the depression; GAD score, general anxiety disorder (GAD) is used to describe the anxiety of the participants; *p*−value, estimates of the effect of interaction on depression and anxiety traits by using ADD × VD traits.
