Introduction

The vitamin D receptor was genetically associated with asthma in 2004 (5) and high intake of vitamin D by pregnant women is associated with about a 50% reduction in asthma risk in the mother’s offspring (6, 7). Together, these exciting findings led to a comprehensive theory about how vitamin D deficiency could influence asthma occurrence through its effects on lung and immune system development and that progressive decreases in vitamin D intake from 1946 onward could be contributing to the epidemic of allergic and autoimmune diseases (1). Vitamin D deficiency is the most common vitamin deficiency in the world today and is particularly prevalent in pregnant women where the fetal lung and immune system are developing (8–10).

To determine the impact of vitamin D supplementation, we performed a clinical trial, the Vitamin D Antenatal Asthma Trial (VDAART) in pregnant women who either had allergies or asthma or had family members with these conditions. The participants were given either 4400 IU of vitamin D3 or 400 IU and were followed throughout their pregnancy and for the first six years of the life of the child (11). The results of the trial were not statistically significant using conventional intent to treat analysis (12, 13) but the reasons were complex, as is often seen when the effects of nutrients are studied. Most importantly, unlike conventional drug trials where one compares drug to no drug or alternative drug, in nutrient trials there is nutrient already present in the placebo group. This creates misclassification that can reduce trial power. When we performed a meta-analysis of the two pregnancy related trials of vitamin D, we got a statistically significant reduction in asthma in the offspring of women who had the higher vitamin D intake during pregnancy and, when we adjusted for the baseline level of vitamin D in the meta-analysis, we got a reduction in asthma risk of 50 percent, exactly what we saw in the observational studies performed previously (14, 15).

To link the results of the VDAART trial directly to the chr17q12 locus we genotyped SNP, rs12936231, in ZPBP2, the gene adjacent to ORMDL3, and stratified the VDAART trial results by maternal genotype at this locus and found that the vitamin D effect in the trial was significantly influenced by genotype with the GG genotype exhibiting a protective effect on asthma risk in the child and the CC or GC genotype not being responsive to vitamin D, additionally this result was related to sphingolipid production with the children protected from asthma having higher sphingolipid production (3).

This SNP (rs12936231) is one of two SNPs (the other being rs4065275) that alter the chromatin state of a regulatory domain, specifically two CTCF binding sites, in the chr17q12-21.1 locus that are correlated with the expression of ORMDL3 (16, 4). 4C-seq assays previously demonstrated that the ORMDL3 promoter interacts with a long-range enhancer in IKZF3 that promotes (or represses) transcription of ORMDL3 in cells expressing both genes and the binding of CTCF per the G allele of rs12936231 in ZPBP2 blocks this interaction resulting in reduced transcription of ORMDL3 on haplotypes with the rs12936231-G allele (4). If rs12936231 is not activated then there is activation of the other CTCF binding site (rs4065275) intronic to ORMDL3 thus favoring the expression of ORMDL3 and increased asthma risk (4).

Mouse models have contributed to our understanding of ORMDL3 function as mice expressing the ORMDL3 transgene exhibit spontaneous increases in airway hyper responsiveness, the essential feature of asthma. This increase in AHR was associated with airway remodeling and peri bronchial fibrosis without airway inflammation (17). ORMDL3 is pleotropic, influencing Ca++ signaling, and the unfolded protein response. Most importantly, increased expression of ORMDL3 inhibits the enzyme serine palmityl-transferase, the rate limiting step in the production of sphingolipids and ORMDL3 TG mice have reduced levels of sphingolipids (18). What has been unclear is how ORMDL3 is controlled.

In the current manuscript we provide evidence for vitamin D regulation of two human chromosomal loci, Chr17q12-21.1 and Chr17q21.2, reliably associated with autoimmune and chronic inflammatory diseases. We demonstrate increased vitamin D receptor (VDR) expression in mouse lung CD4+ Th2 cells, differential expression of these Chr17q12-21.1 and Chr17q21.2 genes in Th2 cells based on vitamin D status and identify the IL-2/Stat5 pathway as a target of vitamin D signaling. Vitamin D induced the expression of the Ikzf3 encoded protein Aiolos to suppress IL-2-signaling and ameliorate cytokine production in Th2 cells thus demonstrating a mechanism for how vitamin D influences autoimmunity.

Materials and methods

In silico analysis of VDR binding sites on Chr17q12-21

For our in-silico analysis, we focused on two genomic regions that are known to harbor variants associated with asthma and autoimmune disease susceptibility. The first one is a 240 kb region that covers parts of the 17q12 and 17q21.1 loci and contains IKZF3, ZPBP2, ORMDL3, GSDMA and GSDMB (chr17: 37,899,254–38,139,253; hg19) and the second one is a 420 kb region in the 17q21.2 locus that includes STAT3, STAT5A and STAT5B (chr17: 40,330,000-40,750,000; hg19). When selecting these genomic ranges, our criterion was to make the window inclusive of any potential enhancer regions targeting the above genes of interest. We used the GeneHancer (19) track in the UCSC Genome Browser (http://genome.ucsc.edu) to determine predicted enhancer-target gene interactions. We obtained VDR, IKZF3, CTCF and Vitamin D-sensitive CTCF binding sites from (20–23), converted their genomic coordinates to hg19 if a different genome build was used in the original publication, and added them as custom tracks on the Genome Browser. We retrieved the known disease and trait associations in these loci from the NHGRI-EBI GWAS Catalog (downloaded on 7/26/2021) (24) and retained all hits that were deemed significant (p<10^{- 5}). We downloaded significant single-tissue eQTLs for IKZF3 and PSMCIP3 in all tissues from the Common Fund (CF) Genotype-Tissue Expression Project (GTEx) Analysis Release V8 (dbGaP Accession phs000424.v8.p2), through the GTEx Portal (https://www.gtexportal.org/home/) in July 2021. We calculated linkage disequilibrium (LD) using the NCI LDmatrix Tool (https://ldlink.nci.nih.gov/?tab=ldmatrix), and by selecting all available populations.

Gene expression analysis in VDAART cohort

To understand the putative mechanism linking IKZF3 expression and Vitamin D to asthma risk we considered umbilical cord blood (CB) RNA-seq data from the VDAART trial. Total RNA was isolated from cord blood using the QIAGEN PAXgene Blood RNA Kit (Qiagen) according to the manufacturer’s protocol. Quality was assessed using the Nanodrop 8000 spectrophotometer (ThermoFisher Scientific). Sequencing libraries were constructed with the TruSeq® Stranded Total RNA Library Prep Globin Kit (Illumina) and the NEXTflex Small RNA sequencing Kit v3 (Boo Scientific). Sequencing was performed using a NextSeq 550 instrument (Illumina). Trimmed reads were mapped to The GRCh38 reference genome using STAR (25). Read counts were computed with htseq (26). The data includes 443 expression profiles of pregnant women with and without asthma undergoing measurements of Vitamin D levels in their CB. Each sample is matched with the observed patterns of wheezing and asthma in the newborn. We filtered all the samples where information about Vitamin D levels, asthma/wheezing patterns or other covariates were missing, obtaining 393 expression profiles. Data was normalized with the voom function of the edgeR R package (27), and Vitamin D levels were binarized using 30 ng/ml as threshold. We performed logistic regression of asthma status as a function of the expression levels of IKZF1 and IKZF3, including their interaction term and the interaction between IKZF3 and Vitamin D. Since IKZF3 forms homodimers and heterodimers with IKZF1, we included IKZF1 in the model. We adjusted the model for maternal asthma status and child race.

Mice

All animal studies were approved by the institutional animal care committee. Mice were maintained in a virus- and parasite-free animal facility under a 12-h light, 12-h dark cycle. C57BL6/J wild type (WT) mice were bred and housed at the MGH animal facility. Vitamin D-deficient mice (C57BL6/J) were bred and maintained in a UV-free environment and weaned onto a similar diet lacking vitamin D metabolites (TD97340; Harlan Teklad) that results in undetectable circulating 25-hydroxyvitamin D levels (28). VDR^-/- mice (Vd^rtm1Mbd/J), on C57BL6/J background, were bred and housed at the MGH animal facility, weaned at 18 d on a diet that maintains normal mineral ion homeostasis in the absence of VDR signaling (2% calcium, 1.25% phosphorus, 20% lactose supplemented diet, TD96348; Harlan Teklad, Madison WI) (29). BALB/c male and female mice were purchased from Charles River Laboratories at 6–8 weeks of age and fed with diets supplemented with varying amounts of vitamin-D3 (TD97340 was used as base diet and supplemented with either 400 IU vitamin D/kg, 1000 IU vitamin D/kg or 4000 IU vitamin D/kg). Mice were bred on the specific diets for two generations to assure stabilization of vitamin D levels. Resulting F2 litters were used for the experiments. Mice were housed under standard conditions with free access to rodent chow and water. All animal experiments were carried out in accordance with NIH Guidelines for the Care and Use of Laboratory Animals as well as guidelines prescribed by the Institutional Animal Care and Use Committee (IACUC) at Brigham and Women’s Hospital (2016N000357), Massachusetts General Hospital (2004N000113) and Harvard Medical Area (05115) (AAALAC 1729).

House dust mite induced-airway inflammation

Mice were exposed to D. pteronyssinus extracts (Greer Laboratories, 10.52 EU/mg endotoxin) via the intranasal application of HDM (25 µg in saline) on days 1,2,3, 8 to 15. Control mice received equal volumes of saline (Hospira Inc., Lake Forest, IL, USA). Twenty-four hours after the final challenge, mice were sacrificed, and specimens sampled.

Bronchoalveolar Lavage

Bronchoalveolar Lavage (BAL) was performed with 2 × 1 mL of PBS/0.02% bovine serum albumin (BSA). The total number of leukocytes was determined using a TC20 Automated Cell Counter (BioRad Systems). To differentiate cell types, cytospins were prepared and stained with Kwik-Diff solutions as per manufacturer’s instructions (Fisher Scientific) and 200 immune cells per slide were counted.

Histological staining

Lungs were perfusion fixed with zinc fixative (BD Pharmingen, San Diego, CA) Tissues were paraffin embedded, divided into sections, and stained with hematoxylin and eosin by MGH DF/HCC Specialized Histopathology Services Core Histology Core. After paraffin embedment, 3-µm sections were stained with hematoxylin and eosin.

Measurements of serum immunoglobulins

Blood samples were collected 24 hours after the final aerosol challenge and serum was prepared. Total IgE levels were measured by ELISA (Chondrex) according to manufacturer’s instructions.

Dissociation of mouse tissues, in vitro polarization of CD4 Th-subsets and flow cytometry analyses

Lungs were perfused via the right heart ventricle with PBS, excised, and digested with collagenase D (2 mg/ml; Roche) in RPM 1640 (Sigma-Aldrich) for 30 minutes at 37°C. The cell suspension was washed twice with PBS. Erythrocytes were removed by hypotonic lysis. Cells were washed twice with PBS and sieved with 30-µm cell strainers (Miltenyi Biotec). The resulting single cells suspension was used to detect leukocyte populations (CD45 BV421 [clone: 104; Biolegend], CD4 PerCPeFluor710 [RM4-5; eBioscience], CD8a APC-Cy7 [53-6.7; Biolegend], CD19 PE-Cy7 [6D5; Biolegend]) and CD4+ Foxp3+ Treg (CD45 BV421, CD4 PerCPeFluor710, Foxp3 PE [FJK-16s; Thermo Fisher]), VDR AlexaFluor647 (D2K6W; Cell Signaling Technologies), Aiolos PE (8B2; Biolegend) and cytokine production by Th-cells (see below). Spleens were excised and mononuclear cells were isolated using a syringe plunger and 100-µm cell strainers. Red blood cells were removed by hypotonic lysis. Tonicity was restored by the addition of RPMI 1640 (Sigma Aldrich). The resulting cell suspension was washed twice with PBS and used for the isolation of naïve CD4+ T cells according to the manufacturer’s instructions (Naive T Cell Isolation Kit II; Miltenyi Biotec). Naive CD4+ T cell purity was assessed by staining for CD4 PerCPeFluor710 (clone: RM4-5), CD44 FITC (IM7), and CD62L PE (mel-14). Naïve T cells were activated with plate bound anti-CD3ε (5 µg/ml, clone: 145-2C11, BD Biosciences) and soluble anti-CD28 (0.5 µg/ml, clone: 37.51, BD Biosciences). Th-subset polarization was induced by exposing naïve CD4+Tcells to following culture conditions for 3 days. Th1 cells: recombinant murine IL-12 (rmIL-12) (10 ng/mL; Peprotech), anti–IL-4 (5 µg/mL; BD Biosciences); Th2 cells: rmIL-4 (20 ng/mL; Peprotech), anti–IFN-γ (5 µg/mL; BD Biosciences), and anti–IL-12 (5 µg/mL; BD Biosciences); Treg: rmTGFβ1 (2 ng/mL; Peprotech), anti–IL-4 (5 µg/mL), anti–IFN-γ (5 µg/mL) for 3 days. Proliferation of Th1 and Th2 cells was induced using rmIL-2 (20 ng/mL; Peprotech) in the presence of Th-polarizing reagents for an additional 2 days. Polarization efficiency was verified by staining for intracellular cytokines (please see below) or Foxp3 expression (see above).

For intracellular cytokine staining, cells were stimulated with phorbol 12-myristate 13-acetate (PMA)/ionomycin and brefeldin A (Cell Stimulation Cocktail; Tonbo Biosiences) for 4 hours. Cells were fixed with Fixation Buffer (Biolegend). Cells were permeabilized using 0.3% saponin in FACS buffer (0.5% FBS in PBS). Antibody staining was performed for 30 minutes in the dark using CD4, IFN-γ (XMG1.2-FITC), IL-10 (JES5-16E3-PE), IL-13 (eBio13A-eFluor660), IL-17A (TC11-18H10.1-PE-Cy7). Cells were washed twice with permeabilization buffer and suspended in FACS buffer. For analysis, 30,000 counts (in the CD4+ gate) were recorded on a FACSFortessa (BD Biosciences) and analyzed with FlowJo software (v10; Tree Star).

Western blot analyses

Cells were lysed in 50 mmol/L Tris/HCl (pH 7.5), 150 mmol/L NaCl, 1% NP-40, and 1× protease and phosphatase inhibitor mix (both Roche) for 15 minutes on ice. Lysates were clarified, and total protein concentrations were measured with BCA (Pierce). Proteins (20 µg) were separated on gradient gels (Mini-PROTEAN^® TGX; BioRad) and transferred on a Immun-Blot PVDF membrane (BioRad). Unspecific binding was blocked using 5 % milk powder in TBST for 30 min. Anti-VDR (D2K6W, 1:1000; Cell Signaling Technologies, CST) was incubated in 0.5 % milk powder in TBST overnight at 4° C. Membranes were washed 3 times with TBST and incubated with secondary anti-rabbit IgG, HRP-linked Antibody (#7074; CST) for 2 hours. To control for equal loading, membranes was probed for β-actin (AC-15; 1:10,000; Sigma-Aldrich) expression for 2 hrs. at RT. Protein bands were visualized with Luminol reagent (Pierce) using the ChemiDoc XRS+system (BioRad).

Immunofluorescence staining for confocal microscopy

In vitro polarized Th2 cells (see above) were stimulated with 100 nm calcitriol (in ethanol) or ethanol [equal v/v] for the final 24 hrs of culture. Cells were washed twice with PBS and spun onto poly-L-lysine coated glass slides. After drying, cells were fixed with 4 % PFA in PBS for 15 min. Cells were permeabilized using 0.2% Triton X-100 in PBS for 20 min. Cells were washed 3 times for 5 min with PBS and incubated with Image-iT FX Signal Enhancer (Thermo Fisher Scientific) for 30 min at RT. After three consecutive washes with PBS, unspecific binding was blocked using blocking buffer (PBS/5 % goat serum/0.3 % Triton X-100) for 60 min at RT. Primary antibodies directed against CD4-FITC (clone; 1:100, BD Biosciences), VDR (D2K6W; 1: 200, CST) were incubated in a humidified chamber over night at 4° C. After 3 washes with PBS, slides were incubated with the secondary reagent, anti-rabbit-Cy3 (1:500; Thermo) for 30 min at RT in the dark. After three more washes with PBS, slides were mounted with ProLong Gold with DAPI according to manufacturer’s instructions (Invitrogen). Images were acquired in a single z-plane using a Zeiss immunofluorescence microscope with a 40x objective.

Total RNA isolation, quality control, and RNA-sequencing analysis

Total RNA was isolated using TRIzol™ Reagent (Invitrogen) according to the manufacturer’s instructions. Quality was assessed using the Nanodrop 8000 spectrophotometer. Sequencing libraries were constructed with the TruSeq® Stranded Total RNA Library Prep Globin Kit (Illumina). Sequencing was performed using a HiSeq. 2500 instrument (Illumina). Trimmed reads were mapped to the GRCm38 reference genome using STAR (25). Read counts were computed with htseq (26). Data were normalized and differential expression was analyzed with the R package DESeq2 (v.1.34.0) (30) using a false discovery rate (FDR) < 0.05.

Pathway enrichment, Overrepresentation and Gene Set Enrichment analysis

KEGG pathway, overrepresentation and gene set enrichment analyses were performed with the “enrichKEGG”, “enrichr” and “GSEA” function in the R package “clusterProfiler” (v.3.10.1) with default settings (31). P-values were adjusted using the BH-FDR correction and enriched terms with a p-adj.<0.05 selected.

Network analysis

RNA-Sequencing data extracted from murine Th2 cells were analyzed using the R package DESeq2 (v.1.22.2) (30). Raw counts were normalized and transformed for variance stabilization using the DESeq2 function rlog (regularized-logarithm transformation). Batch effects were removed using the function removeBatchEffect of the R package LIMMA (32). Gene with counts lower than 5 in at least one sample were removed for downstream analyses. Differential expression analysis was performed between control Th2 cells versus calcitriol stimulated Th2 cells using DeSeq2 (30). The design matrix was built to consider batch effects derived from the RNA isolation process. DESeq2-independent filtering was applied (30) and P values were corrected for multiple hypothesis testing by using the Benjamini-Hochberg (BH)-FDR adjustment (33).

The PPI network used in our network analysis was developed by Silverbush et al. (34). As described in the original paper, this human PPI network was reconstructed by integrating large scale PPI databases, drug response repositories, and cancer genomics data. Only the largest connected component of the network was considered, resulting in an interactome of 15,500 proteins and 234,585 consensus-oriented interactions. Conversion between murine and human gene IDs was performed using the BioMart conversion tool (35). We also assumed a one-to-one correspondence between genes and their protein products in the entire analysis.

To identify the active signaling paths in control and calcitriol stimulated Th2 cells, we modelled the transmission of the molecular signal as a sequential path on the PPI network (called a network path) that starts from a receptor, connects the receptor to multiple intermediate proteins, and leads to a transcription factor. Receptors and transcription factors were selected from a predefined list based on their role in Th2 cell biology and IL-2 expression (Figure 4-table supplement 1). CD3, CD28, IL-4R, IL-13R were selected based on their role in Th2 differentiation. The IL-2R complex was included due to being a target of vitamin D signaling. TGFβRs were selected based on the immunomodulatory activity of this signaling path in CD4+ T cells. role in mediating immunoregulation. Using this subset of receptors and transcription factors as input for the implementation of our network analysis, we computed all the shortest paths connecting each receptor-transcription factor pair on the PPI network. For both control and calcitriol stimulated Th2 cells, we ranked each network path based on the average expression of the path’s nodes. The top selected paths represent the most active downstream processes in Th2 cells in each condition.

mRNA isolation, reverse transcription, and real-time RT-PCR

Total RNA from in vitro polarized Th2 cells was isolated with the RNeasy Mini Kit (Qiagen). First-strand cDNA synthesis was performed with the Omniscript kit (Qiagen). Quantitative real-time PCR was performed with the SsoFast EvaGreen Supermix Kit (BioRad) on a AriaMx Real-time PCR System (Agilent) as per manufacturer’s instructions. Primers sequences used in this study are: Gsdma (sense: 5’-TGC TGT CAA AGG ACG CTA TG-3’; antisense: 5’-TCA ACC AGC TTC AGT TGC AC-3’), IL-2 (sense: 5’-CCC ACT TCA AGC TCC ACT TC-3’; antisense: 5’-ATC CTG GGG AGT TTC AGG TT-3’), Ikfz3 (5’-TGT CAC CTC TGC AAC TAC GC-3’; 5’-TTC CGC AGA ACT CAC ACT TG-3’), Ormdl3 (sense: 5’-CAC ACG GGT GAT GAA C AG TC-3’; antisense: 5’-CTT TGC CTT GGT CTG GAG TC-3’), Stat3 (sense: 5’-GAG GAG CTG CAG CAG AAA GT-3’; antisense: 5’-TCG TGG TAA ACT GGA CAC CA-3’), Stat5a (sense: 5’-TAC ATG GAC CAG GCT CCT TC-3’; antisense: 5’-GTC AAA CTC GCC ATC TTG GT-3’), Stat5b (sense: 5’-TGT GGA TAC AGG CTC AGC AG-3’; antisense: 5’-TGG GTG GCC TTA ATG TTC TC-3’). Relative gene expression levels were calculated with the 2^−ΔΔCt method, with normalization to the murine ribosomal L32 as the housekeeping gene. Data are presented as fold induction over the respective control groups (vehicle or WT). All standard procedures were performed according to the manufacturer’s instructions.

Quantification and statistical analysis

For each figure the number of replicates per experiment is indicated in the corresponding figure legend. In the figures, mean and standard error of the mean (s.e.m.) are presented and error bars represent the mean ± s.e.m. Two group comparisons were analyzed using unpaired two-tailed Student’s t-tests (parametric data) and Mann–Whitney U-test (non-parametric data). Multiple group comparisons were analyzed with Kruskal-Wallis and One-way ANOVA test with Tukey’s post hoc analysis. In vivo experiments were analyzed using Mixed-effect analysis or Two-way ANOVA test with Holm-Šidák’s post-hoc analysis (factors: genotype & exposure, diet & exposure). Prism 9.2 (GraphPad Software) was used to calculate statistics. P values of less than 0.05 were considered statistically significant.

Results

To address the genetic and molecular mechanisms by which vitamin D can influence asthma, allergic and autoimmune disease risk, we studied the genetic loci of two chromosome regions of interest: 1) the 17q12-21.1 region, that includes the well-known asthma-associated ORMDL3 gene (36), and 2) the 17q21.2 region, that includes the STAT5 gene, known to be extensively associated with autoimmune diseases (37, 38). We first examined the entire set of significant disease and trait associations in the Chr17q12-21.1 (181 associations) and Chr17q21.2 (100 associations) loci using the NHGRI-EBI GWAS Catalog (24). Both regions had multiple Th1 and Th2 diseases associations with at least one significant SNP. Chr17q12-21.1 had the highest number of associations related to Th2 diseases (Figure 1A) and 17q21.2 had the highest number of associations related to Th1 diseases (Figure 1-figure supplement 1A). To address the possible link to vitamin D intake, we next investigated a potential overlap of these SNPs with VDR binding sites. Out of the 169 VDR binding sites on Chr17, seven were located in the 17q12-21.1 locus, concentrated in two bands on or near IKZF3 and ZPBP2 (Figure 1B). These VDR binding sites also overlapped with RXRA binding sites (Figure 1B). To assess their functional relevance, we searched for eQTLs in these VDR binding regions. We found that four cis-eQTLs control IKZF3 expression in whole blood and EBV-transformed lymphocytes: rs2941522 and rs12946510 in the enhancer region targeting ORMDL3 and IKZF3, and rs1453559 and rs35564481 in the enhancer region targeting ORMDL3, GSDMA and GSDMB (Figure 1B). Rs1453559 and rs2941522 have previous GWAS associations with both asthma and autoimmune diseases (39–41), rs12946510 has previous GWAS associations with autoimmune diseases (42, 43), and rs35564481 has no known GWAS associations. For the 17q21.2 locus, the only VDR binding site resides near the gene PSMC3IP in an enhancer site targeting STAT5A (Figure 1-figure supplement 1B). PSMC3IP expression is controlled by three eQTLs overlapping with this VDR binding site (rs4793244, rs62078362 and rs111708606), none of which have any known disease associations. These results suggest that vitamin D binding is critical in this genomic region.

Figure 1.

We next examined the linkage disequilibrium (LD) pattern between the Th1/Th2-associated SNPs and the SNPs in the VDR binding sites (Figure 1C). Out of the seven SNPs with VDR binding sites identified in 17q12-21.1 and 17q21.2, the first four in 17q12-21.1 (rs2941522, rs12946510, rs35564481, rs1453559) were in a high-LD block with many other Th1-Th2 associated SNPs in the same region. These same four SNPs in 17q12-21.1 were also in high-LD with two SNPs in ZPBP2 and ORMDL3 (rs12936231 and rs4065275, respectively) that were previously shown to be functionally relevant to asthma (3, 4, 16). Additionally, two of the three SNPs with VDR binding sites in 17q21.2 (rs4793244 and rs62078362) were in strong-LD with three other Th1-Th2 associated SNPs in the17q21.2 region (rs11871801, rs2006141, rs4793090). Overall, these results highlight strong statistical association between vitamin D binding sites and multiple genotypic variations linked to asthma, autoimmune and allergic diseases.

Two SNPs in ZPBP2 and ORMDL3 (rs12936231 and rs4065275, respectively) that were previously shown to be functionally relevant to asthma (3, 4, 16), were also in high linkage disequilibrium with the four SNPs that we identified in the VDR binding regions of the 17q12-21.1 loci (Figure 1C). Importantly, in both the 17q12-21.1 and 17q21.2 regions, the VDR binding sites coincide with CTCF and IKZF3 binding sites, as well as H3K27ac peaks indicating active enhancer sites (Figure 1B, Figure 1-figure supplement 1B). In particular, the two enhancer sites in 17q12-21.1 (GH17J039753 and GH17J039859) were both predicted to interact with the ORMDL3 promoter region (Figure 1B), and the enhancer site in 17q21.2 (GH17J042576) was predicted to interact with the STAT5A promoter region (Figure 1-figure supplement 1B). This suggests a similar mechanism to the one described in (4), whereby SNPs in LD with each other in these enhancer sites affect VDR and CTCF binding, the co-dependence of which in turn could modulate ORMDL3 and STAT5A expression. Although it has been shown that CTCF binding sites can be vitamin D-sensitive (7), no such sites were found in our genomic regions of interest (Figure 1B, Figure 1-figure supplement 1B).

Next, we investigated cord blood gene expression data in the Vitamin D Antenatal Asthma Reduction Trial (VDAART) (see Methods) and looked at the genes of interest in the 17q12 and 17q21 regions and their interaction with cord blood vitamin D levels as predictors of asthma/wheeze risk at age 3 year of age. After adjusting for relevant covariates (see Methods), the baseline expression term of IKZF3 was not found to be significantly associated with asthma risk, while the interaction between vitamin D and IKZF3 was statistically significant (p-val. < 0.05), suggesting that the association of IKZF3 expression with asthma risk was mediated by vitamin D levels in the cord blood (Figure 1-table supplement 1).

Based on our bioinformatic analysis of the Chr17q12 and 17q21 regions, we hypothesize that low vitamin D tissue levels is associated with decreased IKFZ3 expression and increased expression of ORMDL3, STAT3, STAT5A, STAT5B, and IL2.

Vitamin-D receptor (VDR) is expressed in Th2 cells in allergic airway inflammation

To test our hypothesis of vitamin D activity at the Chr17q12-21 locus we utilized a mouse model of allergic airway inflammation and identified vitamin D responsive leukocyte subsets in house dust mite (HDM)-sensitized and challenged mice using flow cytometry (Figure 2A, Figure 2-figure supplement 2A-C). VDR expression was significantly induced in HDM sensitized and challenged mice and almost absent in leukocytes isolated from control (vehicle) lungs (Figure 2-figure supplement 2C) but. Within the CD45+ population, VDR expression was high in CD4+ T cells (p= 0.0357), while only very few CD8+ T cells and CD19+ B cells expressed VDR (Figure 2B, Figure 2-figure supplement 1C HDM). Within the CD4+ T cell population, highest VDR expression was detected in Gata3+ Th2 cells (p= 0.002 vs Tbet+ and p< 0.0001 vs Foxp3+ T cells). VDR expression was very low in Tbet+ Th1 and CD4+ Foxp3+ Treg cells (Fig. 2b). The differential expression of VDR in Th2 cells was further validated by western blot analysis using naïve CD4+ T cells immediately after isolation (0hrs), after activation with CD3/CD28 for 48 hrs and in in vitro polarized Th1, Th2 and iTreg cells with well-defined culture conditions (Figure 2-figure supplement 2A, see Methods section). VDR expression appeared as early as day 1 of Th2 culture conditions and gradually increased during the polarization process (Figure 2-figure supplement 2B). In the absence of the VDR ligand vitamin D, VDR expression localized to the cytoplasm of Th2 cells. Addition of the VDR ligand calcitriol led to a strong translocation of VDR into the nucleus, suggesting functionality of the expressed receptor in the presence of its ligand (Figure 2-figure supplement 2C). These results identify Th2 cells as vitamin D targets in allergic airway inflammation.

Figure 2.

Vitamin D status regulates the expression of key genes on Chromosome17q12-21.1 and Chr17q21.2

Taking cell type specific gene expression and selective expression of VDR in Th2 cells into account, we next analyzed the expression of genes encoded on Chr17q12-21.1 and 17q21.2 loci in WT and vitamin D deficient Th2 cultures as well as in Th2 cells post calcitriol stimulation.

Vitamin D deficiency differentially regulated the expression levels of several Chr17q12-21.1 encoded genes, with lower expression of Ikzf3, Ormdl3 and Gsdma (all p=0.0286) (Figure 2D). Conversely, stimulation with calcitriol induced the expression of the same genes (all p=0.043) (Figure 2E). Expression of the Chr17q21.2 genes Stat3, Stat5a and Stat5b was not affected by vitamin D deficiency (Figure 2F). However, after calcitriol stimulation Stat5a (p=0.0082) and Stat5b (p=0.003) genes were strongly suppressed and Stat3 (p= 0.0048) expression was increased in Th2 cells. (Figure 2G). These results hint to a protective role of vitamin D in allergic airway inflammation, by regulating the expression of the IL-2 downstream signaling molecules Stat5a and Stat5b.

Deficiency in vitamin D signaling augments allergic airway inflammation

To decipher the protective role of vitamin D signaling in allergic airway inflammation, we employed vitamin D deficient and VDR^-/- mice (Figure 3A) (28). HDM sensitization of vitamin D deficient mice led to a stronger reaction to the allergen resulting in more prominent peri-bronchial and perivascular leukocytic infiltrates (Figure 3B). This phenotype was accompanied by lower VDR expression in CD4+ T cells isolated from HDM-exposed lungs compared to WT (Figure 3C). Vitamin D deficiency led to a robust increase in total immunoglobulin (Ig)E levels (Figure 3D). Histological findings were associated with higher total BAL leukocyte (p= 0.0008) (Figure 3E), eosinophil (p= 0.0002) (Figure 3F) and lymphocyte (p= 0.0023) (Figure 3G) numbers in vitamin D deficient mice. Th2 cell numbers (p= 0.0127) were also higher in the HDM exposed mice (Figure 3H) in a manner that was independent of T helper cell skewing (Figure 3-figure supplement 1A, B). Differentiation and recruitment of Th17 and IL-10 producing CD4+ T cells was not affected by vitamin D deficiency (Figure 3-figure supplement 1A, B). To validate these findings from vitamin D deficiency, we next employed VDR^-/- mice (see Methods). As with the vitamin D deficient mice, HDM sensitization and challenge of VDR^-/- mice led to an increased Th2 immune response. The lung phenotype in the VDR^-/- mice was marked by dense leukocytic infiltrates around and mucus plugs in the airways (Figure 3B). The augmented systemic and local immune response increased total IgE titers (p= 0.0489), and total BAL (p< 0.0001 vs WT HDM and p= 0.003 vs Vit D deficient HDM), eosinophil (p< 0.0001 vs WT HDM and Vit D deficient HDM), total lymphocytes (p=0.0194 vs WT HDM and p= 0.0421vs Vit D deficient HDM) and Th2 cells (p= 0.0009 vs WT HDM) (Figure 3D-H). VDR^-/- did not significantly affect recruitment of Treg and Th17 cells (Figure 3-figure supplement 1A, B).

Figure 3.

To explore the mechanisms underlying the augmented Th2-phenotype in the lungs of vitamin D deficient and VDR^-/- mice after HDM allergen challenge, we analyzed the transcriptional profiles of in vitro polarized Th2 cells from WT and vitamin D deficient mice (Figure 3I). Differential expression and enrichment analysis revealed dysregulation of several pathways implicated in Th2-cell activation, cytokine production, proliferation, and survival with prominent changes in the IL-2/ STAT5 - pathway in vitamin D deficient Th2 cells (p= 0.0014) (Figure 3J, Figure 3-table supplements 1, 2). The effects of impaired vitamin D signaling on these pathways was confirmed by qRT-PCR and flow cytometry. Baseline expression of IL-2 was markedly increased at both the RNA (p= 0.0286)) and protein levels (p= 0.0079) in vitamin D deficient and VDR^-/- Th2 cells (Figure 3K, L). In Vitamin D deficient and VDR^-/- Th2 cell in vitro incubation, the elevated IL-2 levels were accompanied by a sharp increase in IL-13 production (VitD deficient: p< 0.0001; VDR^-/-: p= 0.0141) (Figure 3-figure supplement 2).

Taken together, these results indicate a critical role for vitamin D/VDR signaling in regulating type 2 inflammation and identify the IL-2/ Stat5 pathway as a downstream target of vitamin D in Th2 cells.

Vitamin D suppresses the activation of the IL-2/Stat5 pathway and cytokine production in Th2 cells

To ascertain the vitamin D dependent regulation of the IL-2/Stat5 pathway and the impact on the effector program of Th2 cells, we first analyzed the transcriptional profile of WT Th2 cells exposed to calcitriol during differentiation (Figure 4A). Gene ontology enrichment analysis of biological processes revealed regulation of several processes impacting inflammatory responses, including chemotaxis and activation of cells (Figure 4-table supplement 1). Calcitriol stimulation of Th2 cells affected several immune related disease pathways, including asthma (Figure 4-figure supplement 1, Figure 4-table supplement 2). Gene set enrichment analysis (GSEA), performed on the differentially expressed genes, highlighted the negative regulation of IL-2/Stat5 pathway in Th2 cells by calcitriol (Figure 4B, Figure 4-table supplement 3). Calcitriol stimulation suppressed the expression of Il2 (p<0.0001), Stat5a (p<0.0008) and Stat5b (p<0.0008) genes (Figure 4C). Concordantly, levels of activated Itk (phospho-Itk) (p= 0.029) and expression of the IL2Rβ (CD25) (p= 0.0441) were reduced on a per cell basis (Figure 4D).

Figure 4.

Next, the impact of calcitriol stimulation on the effector function of Th2 cells was assessed by flow cytometric analysis for the T2 cytokine IL-13. Activation of the vitamin D/VDR pathway led to significant suppression of IL-13 production in CD4+ T cells (p= 0.0012) (Figure 4E).

To further delineate the molecular impact of vitamin D on Th2 cell biology, we integrated the murine transcriptional profile with a recently developed human Protein-Protein Interaction (PPI) database (34). By modeling the signal transduction process as a sequential path on the PPI network that starts from a receptor and propagates downstream to a transcription factor via multiple intermediate proteins (for details see Methods section, Figure 4-figure supplement 2A), we inferred the active signaling flow in Th2 cells for control (vehicle) and calcitriol stimulated Th2 cells. As input of our network analysis, we selected a subset of receptors and transcription factors that are known to be involved in CD4+ T cell activation, Th2 cell differentiation, and cytokine production (Figure 4-table supplement 4). For clarity of the figure, we display the gene names. Corresponding protein names are attached in Figure 4-table supplement 5. In control Th2 cells, our network analysis revealed the engagement of a variety of receptors including IL4R, IL2RA, CD28, and TNF-receptors. Downstream to these receptors, the most active signaling paths highlighted the activation of STAT3, STAT5A/B, NFkB1, and IKZF1 transcription factors (Figure 4-figure supplement 2B) via multiple adapter proteins of the TRAF and JAK family. Interestingly, these transcription factors are known to be involved in differentiation and activation pathways in Th2 cells Conversely, after calcitriol stimulation, signal from IL4R, CD28 and TNFSF8 propagated over different interaction partners (Figure 4F). Calcitriol specific edges are highlighted in orange color. Molecular mediators included the cytoskeletal proteins VAV1 and RACK1 and emphasized activation of the transcription factors AP-1 (FOS/JUND dimer) and NFATC1.

Our network analysis highlighted a central role of the transcriptional repressor IKZF1 (IKAROS) upon calcitriol stimulation, which was targeted by multiple downstream mediators including IKZF3 (AIOLOS). Notably, the IKZF3 gene encodes for the AIOLOS protein, a direct interaction partner of the IKZF1-encoded IKAROS. Both transcription factors, as homo- and heterodimers, suppress IL-2 expression on the transcriptional levels in mouse and human lymphocytes (44). Calcitriol increased the number of edges connecting intermediate proteins to the transcriptional repressor IKAROS. Among the connecting proteins, AIOLOS interaction with IKAROS was implicated in repression of IL-2 expression in Th2 cells by vitamin D.

Absence of STAT5A/B mediators in the active signaling paths of calcitriol stimulated Th2 cells further suggested repression of the Th2 cell STAT5 pathway with calcitriol compared to control (Figure 4-figure supplement 2B).

Taken together, these results indicate that the specific suppression of the IL-2/Stat5 pathway by vitamin D results in reduced cytokine production by Th2 cells. Our findings suggest that this suppression was mediated by the vitamin D dependent regulation of IKAROS (IKZF1) through the induction of AIOLOS (IKZF3).

Vitamin D supplementation alleviates the allergic phenotype in the lung by suppressing type 2 cytokine production in a dose dependent manner

To test if vitamin D supplementation prevents the development of T2 driven allergic inflammation in the lung, custom rodent diets with select vitamin D doses were used: 400 IU/kg chow (400 IU; low), 1000 IU/kg chow (1000 IU; regular) and 4000 IU/kg chow (4000 IU; high) (Figure 5-figure supplement 1A, see Methods section). The 1000 IU group was used as the control group, as it contains the average amount of vitamin D fortified in standard rodent chow. HDM exposure led to recruitment of leukocytes into the lung tissue and BAL (Figure 5A, B), with a predominant eosinophilic and lymphocytic infiltrate as well as high numbers of IL-13+ Th2 cells (Figure 5C-E). There were trends for higher cellular infiltration in the airways when mice were supplemented with the low vitamin D doses (400 IU), but these inflammatory responses were not significantly different from the phenotype observed with 1000 IU. Lung and airway pathology were similar between these dosing groups, with similar numbers and frequencies of the different leukocyte populations. In sharp contrast, dietary supplementation with the higher vitamin D dose (4000 IU) decreased the inflammatory phenotype in HDM exposed mice compared to the 400 IU and 1000 IU dosed animals. Concordant with the histological findings, BAL leukocytes (p< 0.0001) were decreased with marked reductions in BAL eosinophil (p=0.0016) and lymphocyte numbers (p=0.026) (Figure 5B, D), and the frequency of IL-13+ Th2 cells was significantly reduced in the lungs with 4000 IU of vitamin D (p= 0.0066 vs 400 IU HDM; p=0.0011 vs 1000 IU HDM) (Figure 5E). Of note, dietary supplementation with high vitamin D levels, increased VDR expression in CD4+ T cells from HDM exposed lungs (p< 0.0001) (Figure 5-figure supplement 1B).

Figure 5.

To ascertain that dietary supplementation with vitamin D specifically altered mouse Th2 cell differentiation and the effector program, splenic naïve CD4+ T cells from 400 IU and 4000 IU mice were studied. Of interest, naïve CD4+ T cells did not express VDR until they were exposed to Th2 polarizing conditions. Gata3 expression on per cell basis was lower in 4000 IU Th2 cells (p= 0.0286) (Figure 5-figure supplement 2A). This was accompanied by higher VDR protein expression (p= 0.0286) (Figure 5-figure supplement 2B). The expression of IL-2 (p= 0.0004) and the Th2 cytokines IL-13 (p= 0.0159) and IL-5 (p= 0.0159) were lower in 4000 IU Th2 cultures (Figure 5-figure supplement 2C-E).

To explore the specific induction of Aiolos expression for vitamin D mediated suppression of IL-2 production in Th2 cells, Aiolos expression was monitored throughout Th2 differentiation. Aiolos expression gradually increased under Th2 polarizing conditions in CD4+ T cells, starting as early as day 1 of in vitro culture (Figure 5F). Simultaneous exposure to calcitriol further increased Aiolos protein levels compared to vehicle in Th2 cells (p= 0.0079) (Figure 5G).

To test if vitamin D exerts its inhibitory effect on IL-2 production via Aiolos induction, Ikzf3^-/- mice were employed, exposed to Th2 polarizing conditions and assessed for IL-2 production after calcitriol exposure. IL-2 production in control Ikzf3^-/- Th2 cultures was significantly higher compared to WT Th2 cells. Calcitriol significantly reduced IL-2 levels in both, WT (p= 0.0001) and Ikzf3^-/- Th2 cells (p= 0.0001) (Figure 5H); however, the suppressive effect of calcitriol on IL-2 production was significantly higher in WT Th2 cells (∼ 65% inhibition) compared to Ikzf3^-/- Th2 cells (∼40% inhibition) (p= 0.0001) (Figure 5H).

Taken together, these results indicated that vitamin D can regulate Aiolos expression and suggested a role for Aiolos as a downstream effector for selectively fine-tuning Th2 cell IL-2 production by vitamin D.

Discussion

We utilized human genetics to identify SNPs on chromosome 17q12-21.1 and 17q21.2 that were within GWAS regions for Th1 and Th2 autoimmune diseases. Of note, many of these SNPs were within VDR binding sites. We then determined by eQTL analysis of the genes in these two genomic regions that IKZF3 expression was controlled by four cis-eQTLs in the enhancer region targeting ORMDL3 and IKZF3, and in the enhancer region targeting ORMDL3, GSDMA and GSDMB. Next, we demonstrated that two SNPs in CTCF binding sites were in strong LD with these 4 eQTL SNPs and that the two enhancer sites in 17q12-21.1 were both predicted to interact with the ORMDL3 promoter region, and the enhancer site in 17q21.2 was predicted to interact with the STAT5A promoter region thus, suggesting a plausible human molecular genetic mechanism by which vitamin D could activate IKZF3 to repress ORMDL3 and immune system development. Finally, we leveraged bulk RNAseq in the cord blood of the VDAART trial to show that an interaction for vitamin D level and IKZF3 expression that was associated with reduced asthma risk at 3 years of age.

To confirm our human molecular genetic findings, we chose to utilize mouse models of allergic lung inflammation. We show strong VDR expression in CD4+ Th2 cells and alterations in chr17q12-21.1 and 17q21.2 gene expression induced by vitamin D status. Expression of VDR has been reported for airway epithelial cells and different immune cell types present in the lung (45). In lymphocytes, ex vivo and in vitro studies describe low baseline VDR transcript expression levels that were transiently increased by activating signals and calcitriol stimulation (46). Here, we report a differential VDR expression pattern in CD4+ T cell subsets in allergic inflamed lungs and post lineage specific cytokine stimulation in vitro, with strongest expression in Th2 cells. Activation of VDR suppressed IL-2, IL-5, and IL-13 production by Th2 cells. The expression of VDR in Gata3+Th2 cells in the lung might represent an endogenous control mechanism to curb the pathogenic cell activity and cytokine production after repetitive allergen contact. Ex vivo and in vitro Th1, nTreg and iTreg cells expressed low VDR protein levels at baseline. It remains to be determined whether other Th-subsets described in allergic inflammation express VDR at baseline and if either calcitriol or other secondary signals, specific for each subset, could induce VDR expression during inflammation.

Vitamin D deficiency that had been introduced over two generations and impaired vitamin D signaling led to an augmented allergic phenotype in the lung after allergen exposure with a dominating Th2-signature locally and systemically. This had likely been due to specific deregulation of Th2 cells instead of Th cell skewing, since neither Th17 nor Treg development and recruitment had been significantly affected. In contrast, in mice supplemented with high vitamin D doses over two generations, VDR expression and consequently vitamin D responsiveness was higher in CD4+ T cells from the animals’ lungs, almost abrogating the inflammatory process in the lung after allergen challenge. While the vitamin D dose had been chosen to replicate the VDAART trial in the mouse, the duration of supplementation was adopted to the protocol of generating vitamin D deficient mice. Both factors might influence the outcome here compared to prior rodent studies with inconclusive results.

Differentiation of Th2 cells is driven by IL-4 induced Gata3 expression in naïve CD4+ T cells. This process depends on a positive feedforward loop by IL-2 and its downstream STAT5 signaling pathway for the effective expression of IL-4 (47). Therefore, either interfering with Gata3 expression or the IL-2 pathway represents a promising strategy to finetune Th2-driven immunity. Our data highlight the specific regulation of the IL-2/STAT5 pathway by vitamin D. Vitamin D deficiency and impaired VDR signaling caused elevated IL-2 production by Th2 cells, mounting in high IL-13 production. Elevated IL-13 levels in vivo, could amplify the local T2 responses in the lung. In sharp contrast, exposure of polarizing Th2 cells to calcitriol was able to suppress IL-2 production, expression and activation of STAT5A/B and therefore IL-13 production by Th2 cells. The transcription factor AIOLOS has previously been shown to suppress IL-2 expression on the transcriptional level(44). We report here a vitamin D mediated induction of Aiolos expression in Th2 cells, which is impaired by vitamin D deficiency. High Aiolos expression coincided with lower IL-2 production, ameliorated STAT5A/B expression, and IL-13 production in Th2 cells. In vitro Th2 polarization of naïve CD4+ T cells isolated from mice supplemented with higher vitamin D levels, without further exogenous calcitriol stimulation, developed a weakened Th2 phenotype, characterized by lower IL-2, IL-5, and IL-13 production. Since naïve CD4+ T cells do not express VDR, these results implicate changes in chromosomal accessibility by high vitamin D in the developing immune system, which needs further investigation.

Experiments here with Ikzf3^-/- mice demonstrated significant reductions in vitamin D mediated control of Th2 cell cytokine production, but the vitamin D response was not completely impaired, suggesting additional mechanisms for vitamin D regulation of Th2 cells – several possibilities were uncovered by our PPI analyses. These will be the subject of future research. In addition, we only investigated a Th2 mouse model of asthma and have not extended our mouse work to Th1 models of autoimmunity. We have still to consider other genes that might be linked to these two genetic loci and there are likely other VDR binding sites of importance in both asthma and Th1 autoimmunity. More work to determine the exact levels of vitamin D that confer this protective effect would also be worthy of further investigation.

In summary, we have leveraged clinical trial data, molecular genetic bioinformatic data and mouse model data to outline for the first time, a comprehensive molecular genetic mechanism for how vitamin D influences not only asthma, but both Th1 and Th2 autoimmune disease. The significance of these findings relates to the high prevalence of vitamin D deficiency, especially during pregnancy, and the strong possibility that the epidemic of asthma and autoimmunity might be significantly reduced if vitamin D levels were elevated worldwide.

List of Supplementary Materials

Fig S1 to S11

Tables S1 to S9

Acknowledgements

We wish to thank the VDAART trial participants and clinical staff for their critical work on the VDAART study. We wish to thank also our Channing colleagues for their comments.

Funding

National Heart, Lung, and Blood Institute grant 5P01HL132825

National Heart, Lung, and Blood Institute grant 5UH3OD023268

National Heart, Lung, and Blood Institute grant 5K25HL150336

German Society for Clinical Chemistry and Laboratory Medicine (DGKL) and the German Research Foundation (DFG) grant KI 1868/3-1

Author contributions

A.K. and S.T.W. jointly supervised this work. A.K. and S.T.W. designed the experiments, analyzed data, and wrote the manuscript. A.K., M.G.D., T.R.B., J.R.Q., A.O.S, J.K., T.N., H.Y.P., N.K. and R.E.A performed experiments. A.H., E.M., M.D.M. and R. C. analyzed data. M.D. and K.G. provided critical mouse lines and scientific input. A.A.L. assisted on VDAART related analyses. H.R. and B.D.L. supervised experiments and provided critical scientific input. All authors critically reviewed the manuscript.

Competing interests

Dr. Weiss receives royalties from UpToDate and is an investor in Histolix. All other authors declare that they have no competing interests.

Data and materials availability

All data are available in the main text or the supplementary materials.

Appendix 1 Primer sequences

Appendix 2 Antibodies

Appendix 3