Autocrine vitamin D signaling switches off pro-inflammatory programs of TH1 cells.

The molecular mechanisms governing orderly shutdown and retraction of CD4+ type 1 helper T (TH1) cell responses remain poorly understood. Here we show that complement triggers contraction of TH1 responses by inducing intrinsic expression of the vitamin D (VitD) receptor and the VitD-activating enzyme CYP27B1, permitting T cells to both activate and respond to VitD. VitD then initiated the transition from pro-inflammatory interferon-γ+ TH1 cells to suppressive interleukin-10+ cells. This process was primed by dynamic changes in the epigenetic landscape of CD4+ T cells, generating super-enhancers and recruiting several transcription factors, notably c-JUN, STAT3 and BACH2, which together with VitD receptor shaped the transcriptional response to VitD. Accordingly, VitD did not induce interleukin-10 expression in cells with dysfunctional BACH2 or STAT3. Bronchoalveolar lavage fluid CD4+ T cells of patients with COVID-19 were TH1-skewed and showed de-repression of genes downregulated by VitD, from either lack of substrate (VitD deficiency) and/or abnormal regulation of this system.

Introduction

A significant number of patients with COVID-19 develop severe and life-threatening hyper-inflammation and acute respiratory distress syndrome (ARDS). Mortality from severe COVID-19 remains high, in part due to the limited range of specific immunomodulatory therapies available. Survivors, and those with milder disease, may lose significant tissue function from persistent inflammation and fibrosis, causing chronic lung disease. The efficacy of dexamethasone in reducing mortality indicates the importance of inflammation to disease severity[1]. Improved understanding of the basic mechanisms of COVID-19 will aid rational drug design to reduce morbidity and mortality.

Pro-inflammatory immune responses are necessary for pathogen clearance but cause severe tissue damage if not shut down in a timely manner[2]. The complement system is instrumental in pathogen clearance via recruitment and activation of immune cells[3]. Briefly, complement (C)3, a pro-enzyme, is activated in response to pathogen- or danger-sensing (the lectin pathway), immune complexes (classical pathway) or altered self (alternative pathway) to generate active C3a and C3b fragments, which recruit and activate immune cells and instigate activation of downstream complement components[4]. Complement activation is a pathophysiological feature of ARDS of many etiologies[5] and mediates acute lung injury driven by respiratory viruses[6]. Circulating concentrations of activated complement fragments are high in COVID-19, correlate with severity and are independently associated with mortality[7,8]. Polymorphisms in complement regulators are, likewise, risk factors for poor outcomes[9]. Animal models of other beta-coronaviruses have indicated complement as part of a pathologic signature of lung injury that can be ameliorated by complement inhibition[10]. Emerging clinical trial evidence, from small numbers of treated patients, also points to potential benefits of complement targeting in COVID-19[11].

The complement system is both hepatocyte-derived and serum-effective, but also expressed and biologically active within cells. Notably, activated CD4+ T cells process C3 intracellularly to C3a and C3b via Cathepsin L (CTSL)-mediated cleavage[12]. We have recently shown that SARS-CoV2-infected respiratory epithelial cells express and process C3 intracellularly via a cell-intrinsic enzymatic system to C3a and C3b[13]. This represents a source of local complement within SARS-CoV2-infected lungs, where plasma-derived complement is likely to be absent, and signifies the lung epithelial lining as a complement-rich micro-environment. Excessive complement and IFN-γ-associated responses are both known drivers of tissue injury and immunopathogenesis[14,15]. On CD4+ T cells, C3b binds CD46, its canonical receptor, to sequentially drive T helper (Th)1 differentiation followed by their shut-down, represented by initial production of IFN-γ alone, then IFN-γ together with IL-10, followed by IL-10 alone[16]. Expression of IL-10 by Th1 cells is a critical regulator of Th1 associated inflammation[2]. However, the exact molecular mechanisms governing orderly regulation of Th1 responses culminating in IL-10 expression remain poorly understood and may be critical in the recovery phase of COVID-19 and other Th1-mediated inflammatory diseases.

Vitamin D (VitD) is a fat-soluble pro-hormone carefully regulated by enzymatic activation and inactivation. Most VitD is synthesized in the skin on exposure to ultraviolet B (UVB) radiation from sunlight, then undergoes sequential hydroxylation to 25(OH)VitD and 1,25(OH)2VitD, classically in the liver and kidneys, respectively. VitD has immunomodulatory functions, hence, VitD deficiency is associated with adverse outcomes in both infectious[17] and autoimmune diseases[18]. There are compelling epidemiological associations between incidence and severity of COVID-19 and VitD deficiency/insufficiency[19], but the molecular mechanisms remains unknown.

We found Th1-skewed CD4+ T cell responses in bronchoalveolar lavage fluid (BALF) of COVID-19 patients. As this is a complement-rich microenvironment, we investigated the molecular mechanisms governing orderly shutdown of Th1 responses induced by CD46 engagement. We found that CD46 induces a cell-intrinsic VitD signaling, enabling T cells to both fully activate and respond to VitD. This process was primed by epigenetic remodeling and recruitment of four key transcription factors (TFs), vitamin D receptor (VDR), c-JUN, STAT3 and BACH2. Lastly, we examined these pathways in CD4+ T cells from the BALF of SARS-CoV2-infected patients and found it to be impaired.

Results

COVID-19 CD4+ cells show Th1 and complement signatures

We analyzed single cell RNAseq (scRNAseq) data from the BALF and PBMC of COVID-19 patients and healthy controls (GSE145926, GSE122960 and GSE150728). Because immunity to both SARS-CoV1 and MERS-CoV is mediated by, amongst other cells, IFN-γ-producing CD4+ memory T cells[20] and development of Th1-polarized responses in SARS-CoV2 infection[21] is suspected to contribute to pathogenic hyper-inflammation, we focused our analyses on CD4+ T cells. T cell populations within BALF (Extended Data Fig.1a-b) comprised 5 major sub-clusters, including CD4+ Th cells, according to well-characterized markers (Figs. 1a and Extended Data Fig.1c). Although the proportion of T cells that were CD4+ did not differ between patients and controls (Figs. 1a and Extended Data Fig.1d), 312 genes were upregulated and 134 genes downregulated in patient CD4+ T cells (Fig. 1b and Table S1a). These differentially expressed genes (DEGs) were enriched in noteworthy biological pathways, including interferon-γ response and complement (Fig. 1c and Table S1b). Examination of transcriptional programs by module score indicated that CD4+ T cells in patients were preferentially polarized towards Th1, as opposed to Th2 or Th17 lineages (Fig. 1d). Consistently, expression of core Th1 marker genes were higher in patients (Fig. 1e).

Extended Data Fig. 1

Cellular phenotypes of CD4+ T cells in BALF of patients with COVID-19.

a-b, UMAP representation of scRNAseq showing main clusters of cells from bronchoalveolar lavage fluid (BALF) of patients with COVID-19 and healthy controls (a) and dot plot depicting expression of select marker genes for each cluster (b). Highlighted in both a and b are clusters 2 and 12, which represent T lymphocytes. c, Dot plot showing expression of select marker genes for clusters of cells depicted in Fig. 1a. d, UMAP projection of scRNAseq showing sub-clustering of T cells from bronchoalveolar lavage fluid (BALF) of healthy controls (above) and patients with COVID-19 (below). e, GSEA showing genes more highly expressed in bulk RNA-seq of BALF cells obtained from patients (n=8) with COVID-19 compared to healthy controls (n=20) are enriched in Th1 genes. Box and whisker plots (right) shows the expression of IL10 mRNA in these samples, indicating the median and extending to the minimum, maximum, 25% and 75% quartiles. Data in a-b are from n=9 patients with COVID-19 and n=4 healthy controls, sourced from GSE145926 and GSE122960. Data in c-d are from the same sources but with n=8 patients with COVID-19 and n=3 healthy controls (one sample from each was removed due to too low CD4+ T cell numbers). Data in e are from n=8 patients with COVID-19 and n=20 healthy subjects, obtained from HRA000143. **** p<0.0001 by two-sided Mann-Whitney U-test.

Fig 1

COVID-19 CD4+ T (Th) cells are Th1 skewed.

a, UMAP projection of scRNAseq showing sub-clustering of T cells from bronchoalveolar lavage fluid (BALF) of n=8 patients with COVID-19 and n=3 healthy controls. Stack bars (right) show cumulative cellularities across samples in patients and controls. Dot plot of marker genes for these clusters are shown in Extended Data Fig.1c. b-c, Heatmap showing differentially expressed genes (DEGs; at least 1.5-fold change in either direction at Bonferroni adjusted p-value <0.05 using two-sided Wilcoxon rank sum test) between Th cells of n=8 patients with COVID-19 and n=3 healthy controls (b) and enrichment of those DEGs in Hallmark MSigDB genesets (c). FDR-corrected p-values in c are from hypergeometric tests. Highlighted in red in c are hallmark interferon-γ response and complement pathways. d, Violin plots showing expressions of Th1, Th2 and Th17 specific genes, respectively, summarized as module scores, in BALF Th cells of patients with COVID-19 and healthy controls. Medians are indicated. Exact p-values have been calculated using two-tailed Wilcoxon tests. e, Heatmap showing mean expression of classic Th1 marker genes in BALF Th cells of patients with COVID-19 and healthy controls. Data are sourced from GSE145926 and GSE122960.

Enrichment of complement pathway (Fig. 1c) was notable because a) we recently identified complement as one of the most highly induced pathways in lung CD4+ T cells[22]; b) SARS-CoV2 potently induces complement, especially complement factor 3 (C3), from respiratory epithelial cells[13]; c) COVID-19 lungs are a complement-rich micro-environment[23]; and d) CD4+ T lymphocytes in COVID-19 lungs have a CD46-activated signature[13]. Because CD46 drives both Th1 differentiation and shutdown, characterized by IFN-γ and IL-10 expression, respectively[16], we determined the state of Th1 cells in COVID-19 BALF. IL10 mRNA was dropped out in scRNAseq, but detectable in bulk RNAseq from BALF (Extended Data Fig.1e). Consistently we observed significant enrichment of Th1-related genes in patient cells compared to controls, but ~5-fold lower IL10 (Extended Data Fig.1e). Similar examination within scRNAseq of PBMC (Extended Data Figs.2a-b) did not show meaningful differences in Th1, Th2 or Th17 lineage genes (Extended Data Fig.2c). Collectively, these data indicated the Th1 program and complement signature as features of Th cells at the site of pulmonary inflammation where virus-specific T cells may be concentrated[24], and are consistent the notion that COVID-19 Th1 cells were in the inflammatory phase of their lifecycle compared to healthy controls.

Extended Data Fig. 2

Circulating CD4+ T cells of patients with COVID-19 are not Th1 biased.

a-b, UMAP representation of scRNAseq showing main clusters of cells from peripheral blood mononuclear cells (PBMC) of patients with COVID-19 and healthy controls (a) and dot plot depicting expression of select marker genes for each cluster (b). Highlighted in both a and b are clusters 4, 5 and 6, which represent CD4+ T lymphocytes. c, Violin plots showing expressions of Th1, Th2 and Th17 genes, respectively, summarized as module scores, in PBMC CD4+ T cells of patients with COVID-19 and healthy controls. Data in a-c are from n=6 patients with COVID-19 and n=6 healthy subjects, obtained from GSE150728.

Complement induces an autocrine T cell VitD shutdown program

Prolonged and/or hyper-Th1 activity is pathogenic[14,15]. To discover how shutdown of Th1 cells could be accelerated, we explored how complement regulates Th1 shutdown in healthy cells. CD46, engaged by environmental or intracellularly-generated C3b, works co-operatively with T cell receptor signaling to drive Th1 differentiation then subsequent shutdown[16]. Thus, T cells activated with anti(α)-CD3 and α-CD46 produce IFN-γ, then co-produce IL-10 before shutting down IFN-γ to produce only IL-10[16] (Fig. 2a). T cells secreting neither cytokine (IFN-γ−IL-10− cells) are also appropriately activated as they upregulate activation markers (CD25 and CD69 – Extended Data Figs.3a-b) and proliferate (Extended Data Figs.3c-d). Because this CD46 system is not present in mouse T cells, we explored its function in human CD4+ T lymphocytes. Unless specified otherwise, we used Treg-depleted CD4+ Th cells (CD4+CD25−) throughout. After α-CD3+α-CD46 activation, we flow-sorted cells from each quadrant by surface cytokine capture (Fig. 2a) and performed transcriptome analysis (Extended Data Figs.4a-c). Comparing transcriptomes of IFN-γ+IL-10−, IFN-γ+IL-10+ and IFN-γ−IL-10+ against IFN-γ−IL-10− Th cells, ~2000 DEGs were in common (Fig. 2b and Extended Data Fig.4d and Tables S1c-d). These were enriched for proteins whose molecular function pertained to TF biology (Extended Data Figs.4e-f and 2c and Table S1e), indicating that a key role of CD46 is to regulate TFs. In total, 24 TFs were induced by CD46 in cytokine producing CD4+ cells (Fig. 2d), including VDR (Fig. 2d). VDR was notable for two reasons. First, independent prediction of TFs regulating DEGs of BALF CD4+ T cells and lung biopsies of COVID-19 versus healthy donors returned VDR among the top candidates (Fig. 2e and Table S1f). Second, CYP27B1 was concurrently induced in the transcriptome data (Fig. 2d). CYP27B1 is the 1α-hydroxylase catalyzing the final activation of VitD, converting 25(OH)VitD to biologically active 1,25(OH)2VitD. Inducible expression of CYP27B1 and VDR in Th cells indicated the likely presence of an autocrine/paracrine loop, whereby T cells can both activate and respond to VitD. Although both are described in activated T cells[25,26], the molecular mechanism and consequences are unknown. α-CD3+α-CD28 stimulation of T cells activates C3 processing intracellularly by CTSL, generating autocrine C3b to ligate CD46 on the cell surface[12], which in turn enhances C3 processing to further generate C3b[12]. To establish that CYP27B1 and VDR are induced by complement, we confirmed that α-CD3+α-CD28 and α-CD3+α-CD46 both stimulate these genes in T cells and that this effect was nullified by a cell-permeable CTSL inhibitor, which blocks intracellularly-generated C3b (Fig. 2f). Similarly, T cells from CD46-deficient patients did not upregulate CYP27B1 or VDR with either stimulus (Fig. 2g). Collectively, these data indicated that CD46 ligation by complement induces both the enzyme and receptor to enable Th cells to fully activate and respond to VitD.

Fig. 2

VDR and CYP27B1 are induced by complement and predicted as regulators of the Th1 program in COVID-19.

a, Representative flow cytometry showing IFN-γ and IL-10 in CD4+ T (Th) cells activated with α-CD3+α-CD46 and the four quadrants (A, B, C and D) from which cells were flow-sorted for transcriptome analysis. Live and single cells are pre-gated. b, Venn diagram showing number of DEGs (≥±1.5-fold at unadjusted p-value<0.05 using ANOVA) comparing cells in quadrants B, C and D against A, respectively (n=4 experiments). c, Enrichment of gene ontology molecular function terms in shared DEGs (intersect of Venn diagram in b), ranked by statistical significance. Marked are terms corresponding to transcription factor (TF) activity. d, Heatmap of induced TFs in α-CD3+α-CD46-activated Th cells at each stage of the life-cycle shown in a. Highlighted are VDR and expression of CYP27B1. e, EnrichR-predicted ENCODE and ChEA (ChIP enrichment analysis) TFs regulating the DEGs between COVID-19 vs. healthy donor Th cells (upper panel) and lung biopsies (lower panel). Shown are Benjamini-Hochberg adjusted p-values from hypergeometric tests. f, VDR (left panel) and CYP27B1 (right panel) mRNA in Th cells activated, or not, as indicated, with or without cathepsin L inhibitor (CTSL inh.) (n=5 experiments). g, VDR (left panel) and CYP27B1 (right panel) mRNA in Th cells of a patient with CD46-deficiency, activated, or not, as indicated (n=3 experiments). h-i, Representative flow cytometry (h) and cumulative data from n=6 independent experiments (i) showing IFN-γ and IL-10 in Th cells activated with α-CD3+α-CD46 with, or without, carrier, active [1,25(OH)2D3] or inactive [25(OH)D3] VitD. j, IL10 in Th cells activated with α-CD3+α-CD46 with, or without, inactive [25(OH)D3] VitD, with siRNA targeting VDR, CTSL or CYP27B1, or non-targeting siRNA (NT) (n=5 experiments). Data in a-d are from GSE119416. Data in e, upper panel are from GSE145926 and GSE122960. Data in e, lower panel are from GSE147507. Data in g are from microarrays in[54]. Bars in f-g and i show mean + sem; box plots in j show median value and the range extends from minimum to maximum. All statistical tests are two-sided. *p<0.05 **p<0.01 ***p<0.001 ****p<0.0001 by ANOVA.

Extended Data Fig. 3

Expression of activation markers and proliferation of anti-CD3 + anti-CD46-activated cells.

CD4+ T cells were activated as before using anti-CD3 + anti-CD46 in culture plates. After 48h cells were stained for IFN-γ and IL-10 and co-stained with CD25 or CD69. Separately, CTV-labelled CD4+ T cells were activated in the same manner. a-b, representative flow cytometry histograms (a) and cumulative data from n=3 independent experiments (b) of CD25 and CD69 expression in cells at each stage of cytokine secretion. c-d, representative CTV dilution representing cells that had undergone proliferation (c) and cumulative data from n=3 independent experiments (d). Box and whisker plots in b and d show the medians and extend to the minimum, maximum, 25% and 75% quartiles. * p<0.05 by two-sided ANOVA. All other comparisons were non-significant.

Extended Data Fig. 4

Complement-activated CD4+ T cells are enriched in transcription factors.

a-c, Volcano plots showing differentially expressed genes (DEGs) following activation of CD4+ T cells with α-CD3+α-CD46, comparing IFN-γ+IL-10− cells (a), IFN-γ+IL-10+ cells (b) and IFN-γ−IL-10+ cells (c) to IFN-γ−IL10− cells, respectively. DEGs are defined as at least 1.5-fold change in either direction at unadjusted p-value <0.05 using two-sided ANOVA. Marked in a-c are the IFNG, IL10, VDR and CYP27B1 genes. d, Heatmap showing expression of the 2023 shared DEGs in a-c. e, ClueGo analysis for molecular function terms in the 2023 DEGs shown in d represented as a Cytoscape visualization. Genes are shown in grey, enriched molecular function terms are in red scaled to reflect fold enrichment and edges link genes to molecular function terms. Node sizes reflect enrichment significance. Related terms grouped as families within yellow circles. Four such families represent transcriptional regulation of gene expression and are shown in the inset on the right. f, The top 14 transcription factor molecular function terms are shown, with associated fold enrichments and FDR q-values. Data in a-f are from n=4 experiments.

As proof of principle that this autocrine/paracrine intracellular VitD system is involved in Th1 shut-down, we stimulated CD4+ T cells with α-CD3+α-CD46 with or without active [1,25(OH)2] or inactive [25(OH)] VitD. Both forms of VitD significantly repressed IFN-γ and increased IL-10 (Figs. 2h-i), indicating that the cells had acquired the ability to activate and respond to VitD. We cultured similarly stimulated cells with inactive VitD and silenced CTSL, VDR or CYP27B1 using siRNA and assessed IL10 transcription. A significant reduction in IL10 was evident on silencing any of these components (Figs 2j and Extended Data Fig.5), indicating that they are all required for producing IL-10. Altogether, these data indicated existence of a novel complement-induced intracellular autocrine/paracrine VitD system promoting T cell shutdown.

Extended Data Fig. 5

Efficiency of silencing of indicated targets by siRNA from Fig 2j.

Shown are cumulative data from n=3 independent experiments. Bars show mean + sem. NT, non-targeting. **p<0.01 ***p<0.001 by unpaired two-sided t-test.

VitD induces IL-10 via autocrine/paracrine IL-6/STAT3

We next investigated the effects of VitD on CD4+ T cells. In all subsequent experiments cells were activated with α-CD3 and α-CD28 (as this stimulus also signals through CD46[12]), and we used active 1,25(OH)2VitD. We confirmed that VDR protein was induced by T cell activation (Extended Data Fig.6a), VitD enhanced VDR expression[27] (Extended Data Fig.6a), and VitD-bound VDR translocated to the nucleus (Extended Data Figs.6b-c). VitD upregulated 262 and downregulated 128 genes in Th cells (Fig. 3a and Tables S2a-b), which were not due to alterations in cell proliferation or death (Extended Data Fig.6d). Classical VitD-regulated genes, including CTLA4, CD38 and CYP24A1, were induced and both type 1 (IFNG) and type 17 (IL17A, IL17F, IL22 and IL26) cytokines were repressed (Figs. 3b-c), consistent with previous reports[28]. Genes induced by VitD included two cytokines, IL10 and IL6, and several TFs, including JUN, BACH2 and STAT3 (Figs. 3b-c). As IL-10-producing CD4+FoxP3−type 1 regulatory T (Tr1) cells can be induced from naïve Th cells by VitD[29], we noted that genes induced/repressed by VitD were not enriched in Tr1 signature genes (Extended Data Fig.6e) and did not exhibit the archetypal surface phenotype of Tr1 cells (CD49b+LAG3+30) (Extended Data Fig.6f). Similarly, FOXP3, the master TF of natural Tregs, was not upregulated by VitD (Table S2b). Genes regulated by VitD were most significantly enriched for cytokines (Extended Data Fig.6g). Unexpectedly, IL6, normally a pro-inflammatory cytokine, was the most highly induced gene in the hierarchy of cytokines regulated by VitD (Extended Data Fig.7a). We confirmed repression of IFN-γ and IL-17 and induction of IL-10 and IL-6 proteins in VitD-treated Th cells (Extended Data Fig.7b) and noted a strong dose-response relationship between VitD concentrations and these effects (Fig. 3d). To confirm again the autocrine/paracrine VitD activation system at the protein level (Extended Data Fig.5c), we observed repression of IFN-γ and IL-17 and induction of IL-10 and IL-6 by Th cells treated with only inactive 25(OH)VitD, indicating intracellular conversion of 25(OH)VitD to active 1,25(OH)2VitD (Extended Data Fig.7d).

Extended Data Fig. 6

Phenotype of Vitamin D treated cells.

a, Western blot (left) and cumulative data (right) for VDR at days 1, day 3 and day 5, with Hsp90 as loading control, in both carrier and VitD treated CD4+ T cells. b, Representative immunoblots for VDR and indicated housekeeping proteins in nuclear and cytoplasmic extracts of carrier and VitD treated CD4+ T cells. c, Co-localization of VDR and DAPI in carrier and VitD treated CD4+ T cells, measured on day 2 using ImageStream. Shown are representative frequency histograms indicating overlap between VDR and DAPI in the entire population (left), and cumulative data from n=3 independent experiments (right). d, Cell death assessed by live/dead stain and proliferation assessed by CFSE dilution in CD4+ T cells treated with carrier or VitD after 3 days of culture. e, GSEA showing genes more highly expressed in CD4+ T cells treated with carrier compared to VitD are enriched in Tr1-induced genes (left panel) and genes more highly expressed in CD4+ T cells treated with VitD compared to carrier are enriched in Tr1-repressed genes (right panel; curated from GSE139990). f, Representative flow cytometry plot showing CD49b and LAG-3 expression in CD4+ T cells, with carrier and VitD treatment (left), and quantification of cumulative data (right). g, Top 5 MSigDB canonical pathways enriched in DEGs of VitD vs carrier treated CD4+ T cells (see Fig. 3a). Unless indicated, all cells in Fig. S6 have been activated with α-CD3+α-CD28. Bars represent mean + sem throughout. All experiments have been carried out n=3 times. Shown in e are unadjusted empirical p-values; NES = normalized enrichment score. *p<0.05, ****p<0.0001 by 2-way ANOVA (a, f) and paired t-test (c). All statistical analyses are two-sided.

Fig. 3

Vitamin D induces IL-10 in Th cells by enhancing IL-6-STAT3 signal transduction.

a, Number of DEGs between VitD and Carrier-treated Th cells (≥±1.5-fold change at FDR<0.05). b, Scatter plot showing mRNA expression (RPKM) of genes in Th cells treated with VitD or carrier. VitD-induced and -repressed genes are depicted in red and blue, respectively. Noteworthy genes are annotated (black), including classical VitD-induced genes (orange). n=3 independent biological experiments. c, Heatmap showing expression of select genes from b. d, Dose-response of indicated cytokines from Th cells treated for 72h with VitD. Stars indicate statistically significant changes in comparison to 0nM of VitD; n=3 experiments. e, Pearson correlation between IL-6 and IL-10 concentrations in culture supernatants of VitD-treated Th cells. Shown is the correlation line, plus 95% confidence interval. f, IL-10 concentrations in supernatants of Th cells cultured with VitD, with and without Tocilizumab; n=3 experiments. g, Representative image from n=2 experiments of a phospho-kinase array (array of 43 kinases in duplicate spots) carried out on 3-day lysates of carrier- or VitD-treated Th cells. Location of STAT3 phosphorylated at lysine 705, c-JUN phosphorylated at serine 63 and reference spot (to which all spots are normalized) are indicated. h, Heatmap showing normalized phosphorylation values of differentially phosphorylated proteins following VitD-treatment (please see also Extended Data Fig.7i) in n=2 donors. i, Immunoblots of lysates of Th cells treated with carrier or VitD with, and without, Tocilizumab (Toc) at the concentrations shown. Shown are representative images from n=3 experiments (quantified in Extended Data Fig.7j). j, IL-10 production from Th cells cultured with carrier or VitD, with or without a STAT3 inhibitor (STAT3i). Genotype of cells (WT, STAT3 WT; DN, STAT3 dominant negative) is indicated. n=3 experiments; each dot represents an individual donor. k, IL-10 production from Th cells transfected with control siRNA or siRNA targeting STAT3. n=5 experiments. Unless indicated, all cells in Fig. 3 have been activated with α-CD3+α-CD28. Cumulative data in d, f and j-k depict mean+sem. All statistical tests are two-sided. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001 by one-way (d, j, k) and two-way ANOVA (f).

Extended Data Fig. 7

Hierarchy of cytokines induced by Vitamin D and modulation of IL-6 effects and the phosphor-kinase proteome by Vitamin D.

a, Differentially expressed genes (DEGs) between VitD and carrier treated CD4+ T cells (see Figs. 3a-b) ranked by fold change. Each DEG is marked by a blue dot; each differentially expressed cytokine is marked by an orange dot. Select cytokines have been labelled. b, Cytokine concentrations in supernatants of CD4+ T cell cultures after 5 days of treatment with carrier or VitD. c, Heatmap showing mRNA expressions (log2 TPM) of the 25-hydroxylase enzymes (CYP2R1 and CYP27A1), the 1α-hydroxylase enzyme (CYP27B1) and the 24-hydroxylase enzyme (CYP24A1), responsible for the two steps of Vitamin D activation and its subsequent inactivation, respectively. Data are from CD4+ T cells activated with α-CD3+α-CD28 and cultured with either carrier or VitD, or left unactivated. d, Concentrations of indicated cytokines in culture supernatants of CD4+ T cells treated with escalating doses of 25(OH)VitD for 72h. e, IL6 mRNA, fold change compared to day 1 carrier (above), and IL-6 protein concentration in matched supernatants (below), at days 1, 3 and 5 in carrier and VitD-treated CD4+ T cell cultures. f, Representative flow cytometry plot (left) and cumulative data (right) of intracellular IL-6 expression in T cells treated with carrier or VitD (assay carried out on day 3). Cells gated based on lymphocyte gate (forward scatter, side scatter), singlets, live cells and CD4+ cells. g, IL-10 concentrations in supernatants of CD4+ T cells cultured in the presence of increasing concentrations of IL-6 for 72 hours. h, IL-17 concentrations in supernatants of CD4+ T cells cultured in the presence of increasing concentrations of IL-6, with and without VitD for 72 hours. i, Volcano plot representing changes in protein phosphorylation on phospho-kinase array comparing VitD-treated versus carrier treated cells. Data are from n=2 independent experiments. Thresholds for significance have been set at 1.2 fold change in phosphorylation in either direction at p-value <0.05. Please also see Figs. 3f-g. Marked are phosphoproteins that show significant changes in phosphorylation on VitD treatment. j, Quantification of pY-STAT3, STAT3, p-c-JUN, c-JUN and Hsp90 from immunoblots of lysates of CD4+ T cells treated with carrier or VitD with, and without, Tocilizumab (Toc) at the concentrations shown. Bars show mean + sem from n=3 independent experiments. Please also see Fig. 3h. k-l, Shown are representative Western blots (k) and quantifications (l) of STAT3 and HSP90 in Carrier- and VitD-treated CD4+ T cells in the presence and absence of STAT3 siRNA. Data are representative of n=2 experiments carried out. Unless indicated, all cells in Fig. S7 have been activated with α-CD3+α-CD28. Cumulative data in b, d-h, j depict mean + sem. Unless indicated, all experiments have been carried out n=3 times. All statistical tests are two-sided. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001 by one-way ANOVA (d, j), two-way ANOVA (b, e, h) and paired t-test (f). Statistical comparisons in d and h compare VitD-treated (d) or IL-6-treated (h) cells against untreated cells.

IL-6 is a pleiotropic cytokine expressed by most stromal and immune cells[31]. Yet it is not commonly attributed to CD4+ T cells. We established that IL6 mRNA and protein were produced by T cells and induced by VitD (Extended Data Fig.7e) and that VitD treatment significantly increased intracellular expression of IL-6 (Extended Data Fig.7f). In these experiments, there was a strong correlation between IL-6 and IL-10 produced in response to VitD (Fig. 3e), suggesting a potential causal relationship. Accordingly, we stimulated Th cells with VitD, with or without Tocilizumab, an IL-6 receptor (IL-6R) blocking antibody used clinically to treat IL-6-dependent cytokine release syndrome, including that seen in COVID-19[32]. Tocilizumab significantly impaired IL-10 produced by VitD, indicating that IL-6R signaling induced by autocrine/paracrine IL-6 promotes IL-10 in VitD-treated Th cells (Fig. 3f). IL-6 can cooperate with IL-27[33] or IL-21[34] to promote IL-10 production in mouse T cells. However, both these cytokines were repressed by VitD in our transcriptomic analyses (Table S2b), so it did not appear likely that these cytokines cooperate with IL-6 to induce IL-10 in this setting. Addition of exogenous IL-6 without VitD did not induce IL-10 (Extended Data Fig.7g) but increased pro-inflammatory IL-17 (Extended Data Fig.7h), as has been reported[35]. These data indicate that pro-inflammatory functions of IL-6 may be restricted or averted by the production of anti-inflammatory IL-10 in the presence of VitD in human Th cells (Fig. 3f).

Genes differentially expressed by VitD were enriched in cytokine signaling pathways, which are commonly mediated through phosphorylation of signal-dependent TFs, including JAK-STATs (Extended Data Fig.6g). We therefore carried out a phospho-kinase protein array using lysates of carrier and VitD-treated cells. Five proteins showed significant differences in phosphorylation between carrier and VitD, most notably c-JUN and STAT3, both of which were significantly more phosphorylated (Figs. 3g-h and Extended Data Fig.7i). Immunoblotting confirmed induction of both STAT3 protein and STAT3 phosphorylation by VitD (Fig. 3i and Extended Data Fig.7j). IL-6 potently drives STAT3 activation by phosphorylation[31]. STAT3 phosphorylation induced by VitD was abrogated by blockade of the IL-6R with Tocilizumab, indicating that VitD-induced IL-6 is responsible for STAT3 phosphorylation (Figs. 3i and Extended Data Fig.7j). Conversely, STAT3 protein expression was dependent on VitD but independent of IL-6 signaling, since IL-6R blockade did not reverse its induction by VitD (Figs. 3i and Extended Data Fig.7j). By contrast, both c-JUN expression and phosphorylation were dependent on VitD and mostly independent of IL-6 (Fig. 3i and Extended Data Fig.7j). As VitD-induced STAT3 phosphorylation was mediated by IL-6, we investigated whether STAT3 drives IL-10 produced by VitD. Both a cell permeable STAT3 inhibitor and knock-down of STAT3 by siRNA significantly impaired IL-10 produced in response to VitD (Figs. 3j-k and Extended Data Fig.7k-l). Likewise, VitD failed to produce significant IL-10 from Th cells of patients with hyper-IgE syndrome, that have dominant negative STAT3 mutations and are unable to transduce STAT3 signaling normally (Fig. 3j). Collectively, these data established that VitD induces STAT3 and IL-6, and autocrine/paracrine IL-6R engagement phosphorylates STAT3, which promotes production of IL-10.

VitD alters epigenetics and recruits c-JUN, STAT3 and BACH2

VitD-bound VDR interacts with histone acetyl transferases, transcriptional co-activators, co-repressors, and chromatin remodeling complexes to modulate transcription. We explored the effects of VitD on T cells epigenetic landscapes, using memory CD4+ T cells (these cells express VDR without requiring pre-activation). We profiled genome-wide histone 3 lysine 27 acetylation (H3K27Ac), a marker of active regions of the genome, using cleavage under target and release using nuclease (CUT&RUN) in VitD and carrier-treated cells. VitD induced dynamic changes in histone acetylation genome-wide (Figs. 4a-b), indicating significant changes in enhancer architecture. By 48h after VitD treatment, ~25000 and ~21000 H3K27Ac peaks were induced and repressed, respectively (Figs. 4c-d and Table S3a). Loci of genes transcriptionally induced by VitD showed increased histone acetylation and those repressed by VitD had reduced histone acetylation (Fig. 4c). The size of these peaks was significantly affected by VitD, leading to generation of new super-enhancer (SE) architectures and promotion of existing SEs (Fig. 4e). SEs are complex regulatory domains comprised of collections of enhancers critical for regulating genes of particular importance to cell identity and risk of genetic disease[36,37]. VitD-modified SEs included those associated with BACH2, STAT3, IL10 and other genes induced by VitD (Fig. 4e and Table S3b). To identify potential TFs recruited to these loci, we carried out TF DNA motif finding at H3K27Ac peak loci induced by VitD. The top enriched motifs were VDR, AP-1 family members, notably c-JUN and BACH2, and STAT3 (Fig. 4f). All three of these were TFs transcriptionally induced by VitD (Figs. 3b-c) and two of these, c-JUN and STAT3, were also post-transcriptionally more phosphorylated after VitD (Figs. 3g-h). Thus, we reasoned that they are likely to play an important role in gene-regulation by VitD.

Fig. 4

Vitamin D reshapes epigenetic landscape of Th cells.

a, Genome-wide H3K27Ac CUT&RUN peaks 45mins, 48h and 72h after VitD or carrier-treatment of Th cells. b, differential H3K27Ac peaks (signal ≥0.2, ≥1.5-fold change) after VitD or carrier-treatment of Th cells at the indicated time points. c, Scatterplot showing H3K27Ac CUT&RUN peak signal intensities 48h after VitD or carrier-treatment of Th cells. Indicated are VitD-induced peaks (red) and VitD-repressed peaks (blue). Highlighted are select peaks at loci of interest. Data show a representative example from n=2 independent experiments. d, Heatmaps showing H3K27Ac signal at VitD repressed and VitD-induced peaks (below) and histograms showing normalized signals in carrier and VitD treated cells (48h) above. e, Ranked order of H3K27Ac-loaded enhancers induced by VitD in Th cells after 48h. Super-enhancers (SEs) are indicated. Marked are the relative positions, ranked according to signal intensity (higher = greater signal intensity), of enhancers attributed to selected genes. f, Enriched transcription factor (TF) DNA motifs at H3K27Ac peak loci induced by VitD. Shown are TF families on the left and representative TF members enriched in the data on the right. Unless indicated, all in vitro T cell experiments depicted in Fig. 4 have been activated with α-CD3+α-CD28.

We carried out c-JUN CUT&RUN and VDR, STAT3 and BACH2 cleavage under target and Tagmentation (CUT&Tag) in VitD or carrier treated cells (Extended Data Fig.8a and Table S4) to identify genome-wide distribution of these TFs. Genome-wide binding of all four TFs was increased after VitD (Fig. 5a), and they were recruited to H3K27Ac peak loci induced by VitD (Fig. 5b and Extended Data Fig.8b). VDR, c-JUN, STAT3 and BACH2 each bound ~20-40% of genes differentially expressed by VitD, significantly higher than other loci in the genome (Fig. 5c). Indeed, ~60% of DEGs were bound by at least one of these TFs (Fig. 5d), half of which were bound by more than one (Fig. 5e). Genes bound by two or more TF included STAT3, IL10 and BACH2, (Fig. 5f) indicating the complexity of gene regulation downstream of VitD exposure and interaction of multiple TFs within a gene regulatory network (Fig. 5g). Collectively, these data showed reshaping of the epigenome by VitD, generating new and augmenting existing enhancers, and recruitment of transcriptional regulators to these loci to modify transcriptional output. BACH2, IL10 and STAT3 were exemplars of loci at which VDR binding generated new SEs, to which c-JUN, BACH2 and STAT3 were recruited (Fig. 5h and Extended Data Figs.8c-d) and transcriptional output was increased (Figs. 3b-c).

Extended Data Fig. 8

Vitamin D recruits key transcription factors.

a, histograms showing average cuts per bp in relation to the summits of c-JUN, STAT3, VDR and BACH2 peaks in VitD-treated CD4+ T cells. Shown are data from CUT&RUN (c-JUN) and CUT&Tag (VDR, STAT3 and BACH2) carried out with IgG (cyan) or antibodies specific to each TF (dark blue). p-values from K-S tests are indicated. b, heatmaps showing H3K27Ac and c-JUN signals at VitD-repressed and VitD-induced peaks over time. Time points are indicated. c-d, genome browser tracks at the CTLA4 (c) and STAT3 (d) loci showing H3K27Ac, c-JUN, BACH2, STAT3 and VDR binding in Carrier and VitD-treated cells. Red and blue dots represent peaks in Carrier and VitD-treated cells, respectively. SE denotes super-enhancer regions. Track heights are indicated on the left corner for each track. All cells in Fig. S8 have been activated with α-CD3+α-CD28.

Fig. 5

VitD recruits key transcription factors to shape transcriptional output.

a, Histograms (above) and heatmaps of c-JUN, VDR, STAT3 and BACH2 bound loci in Carrier and VitD-treated Th cells, centered on the peaks. p values by two sample K-S test comparing Carrier to VitD are shown. b, Histograms of c-JUN, VDR, STAT3 and BACH2-bound loci centered on VitD-induced H3K27Ac peaks. p-values by two sample K-S test comparing Carrier to VitD are shown. c, Proportion of genes differentially expressed after VitD treatment (DEGs, from n=3 experiments) bound by each indicated TF. ****p<0.0001 by two-sided Fisher exact test compared to all genes. d, Proportion of DEGs bound by 0, 1, 2, 3 or all 4 of the TFs profiled. e, Venn diagram showing overlap in TF binding between VitD DEGs. f, Heatmap of DEGs showing binding of genes by c-JUN, VDR, STAT3 and BACH2. Select genes have been highlighted on the right. g, Network diagram showing TF binding of genes differentially expressed after VitD treatment. Arrows join TFs to bound genes. Heatmap scale indicates fold change expression after VitD treatment compared to Carrier. h, Genome browser tracks at the BACH2 and IL10 loci showing H3K27Ac, c-JUN, BACH2, STAT3 and VDR binding in Carrier and VitD-treated cells. Red and blue dots represent peaks in Carrier and VitD-treated cells, respectively. SE denotes super-enhancer regions. Track heights are indicated on the left corner for each track. i, ChIP-qPCR for VDR or IgG in VitD-treated Th cells from two donors. Anti-VDR ChIP fragments were probed by qPCR for enrichment of promoters of CYP24A1 (left), STAT3 (middle) and IL6 (right). Data are shown separately for each donor; bars show mean+sem of qPCR. All p-values are from two-sided tests.

We noted that STAT3 is directly bound by VDR (Figs. 5f-g and Extended Data Fig.8d), but we did not find a CUT&Tag VDR peak proximal to the IL6 locus. This may be because IL6 is a lower affinity target, has a distal VDR binding site or because VDR binds this site earlier than the time point at which we carried out CUT&Tag. Thus, because ChIP-qPCR is more sensitive than genome-wide techniques when applied to individual loci, we performed qPCR for STAT3 and IL6 promoters, as well as CYP24A1 (a positive control locus), in anti-VDR ChIP fragments (Fig. 5i). We found strong enrichment of CYP24A1 and STAT3 promoters and moderate enrichment of the IL6 promoter in anti-VDR ChIP fragments, indicating that these loci are all directly bound by VDR, and that VDR binding at IL6 is at lower affinity than the other two loci (Fig. 5i). In summary, these data indicated VitD-induced dynamic changes in genome-wide enhancer architecture and recruitment of several TFs to loci that could explain ~60% of the VitD-dependent variance in gene expression.

BACH2 regulates the VitD response in CD4+ T cells

BACH2 is a critical immunoregulatory TF[38,39]. We confirmed that VitD induces BACH2 protein expression in Th cells (Fig. 6a). Psoriatic skin is rich in CD4+ T cells, psoriasis severity is associated with low active VitD levels[40] and is frequently treated successfully with VitD[41]. We therefore performed confocal imaging on the dermis of patients with psoriasis treated, or not, with VitD. VitD-treatment significantly increased numbers of BACH2+ cells and greater numbers of intranuclear foci of BACH2, compared to untreated skin (Figs. 6b-e), confirming that BACH2 is induced by VitD in vivo. Genes regulated by VitD were ~2-3-fold enriched for BACH2 binding than those not regulated by VitD (Fig. 6f-g). Indeed, BACH2-induced genes and BACH2-repressed genes were more highly expressed in VitD-treated and carrier-treated cells, respectively (Extended Data Figs.9a-b), indicating a BACH2 signature in VitD-regulated genes. Among the most highly enriched in the leading edge of BACH2-induced genes was the IL-6 receptor alpha chain (IL6R) (Extended Data Figs.9a and c). These observations suggested that a significant portion of VitD-driven transcription is BACH2-dependent. Thus, we compared transcriptomes of VitD-treated Th cells from healthy control (BACH2) to those from a BACH2 haploinsufficient patient (BACH2)[36]. VitD-induced genes were significantly enriched in the transcriptomes of wildtype cells (Fig. 6h and Table S5a), indicating that normal BACH2 concentrations are required for appropriate regulation of VitD-induced genes (Extended Data Fig.9e). The same pattern of enrichment was not found for VitD-repressed genes (Extended Data Fig.9d), potentially indicating that half-normal concentrations of BACH2 are sufficient for repression of VitD targets. Thus, some, but not all, type 1 and type 17 inflammatory cytokines (IFNG and IL17F, but not IL17A) were repressed by VitD in BACH2 haploinsufficient cells (Extended Data Fig.9e and Table S5b).

Fig. 6

BACH2 is an important regulator of the Vitamin D response in Th cells.

a, BACH2 Immunoblot in VitD- and carrier-treated Th cell lysates. Shown is a representative example from n=3 experiments. b-e, Representative dermal images of lesional skin from psoriasis patients with (n=2) and without (n=3) VitD supplementation stained for BACH2 (red) and CD3 (green), showing overview (leftmost) and zoomed images (right three images) (b), number of nuclei/image (c), percentage of BACH2+ cells relative to nuclei frequency/image (d) and average number of BACH2 foci/cell (e). For c-e 5-6 images were acquired for each sample and are shown as median values, with minimum, maximum, 25% and 75% quartiles. f-g, Pie-charts comparing percentage of VitD-regulated genes bound and regulated by mouse (f) or human BACH2 (g) against all genes in the genome. The two-sided Fisher exact p-value is shown. h, GSEA showing enrichment in VitD-induced genes comparing transcriptomes of VitD-treated wild-type control (BACH2WT/WT) with BACH2-haploinsufficient Th cells (BACH2WT/L24P). Shown is the empirical p-value; NES = normalized enrichment score. i, Sankey diagram showing the relationship between genes up- and down-regulated by VitD in BACH2 BACH2WT/L24P compared to BACH2WT/WT Th cells, and the binding of those genes by c-JUN, VDR, BACH2 and STAT3. j, Pie-charts comparing percentage of genes up- and down-regulated by VitD in BACH2 sufficient cells that are differently regulated in BACH2-haploinsufficiency. The percentage of those genes bound by BACH2 are shown underneath. k, Top 3 enriched MSigDB canonical pathways pertaining to the genes bound by BACH2 and differently regulated in BACH2-haploinsufficient cells. Indicated are 3 relevant genes contributing to enriched pathways. l, Heatmap showing expression patterns of select VitD-regulated genes in the transcriptomes of carrier- and VitD-treated BACH2WT/WT and BACH2WT/L24P Th cells. m, Genome browser tracks at IL6R locus showing H3K27Ac, c-JUN, BACH2, STAT3 and VDR binding in Carrier and VitD-treated cells. Red and blue dots represent peaks in Carrier and VitD-treated cells, respectively. Track heights are indicated on the left. In vitro T cell experiments depicted in Fig. 6 have been activated with α-CD3+α-CD28. *p<0.01, ****p<0.0001 by unpaired two-tailed t-test.

Extended Data Fig. 9

A subset of Vitamin D-regulated genes are dependent on BACH2.

a-c, GSEA comparing the transcriptomes of carrier and VitD-treated CD4+ T cells against BACH2-bound BACH2-induced genes (a) and -repressed genes (b). Shown in c are the leading edges of the two GSEA enrichment plots in a-b. Marked in a and c is the IL-6 receptor (IL6R). d, GSEA comparing enrichment in VitD-repressed genes of the transcriptomes of VitD-treated BACH2WT/WT haplo-sufficient and BACH2WT/L24P haplo-insufficient CD4+ T cells. e, Circos diagram showing VitD-induced and repressed genes in VitD-treated BACH2WT/WT haplo-sufficient and BACH2WT/L24P haplo-insufficient CD4+ T cells. Cords join shared genes in patient and control. Indicated are the shared VitD-induced genes (CYP24A1, CD38 and IL6) and genes only induced by VitD in the presence of two wild-type copies of BACH2 (IL10 and IL6R). f, Genome browser tracks showing Bach2 ChIP-seq at the IL6ra locus and expression of IL6ra mRNA in CD4+ T cells of Bach2 wild-type (Bach2) and knock-out (Bach2) mice. Track heights are indicated on the left corner for each track. Source data are from GSE45975. Empirical p-values are shown in a-b and d. NES = normalized enrichment score.

To better understand the VitD transcriptional regulatory network, we integrated VitD up- and down-regulated genes in BACH2 and BACH2 CD4+ T cells together with TF binding from CUT&Tag and CUT&RUN (Fig. 6i). 75% and 66% of genes normally up- and down-regulated by VitD, respectively, were not regulated by VitD in BACH2 haploinsufficient cells. 28% and 14% of these, respectively, were directly bound by BACH2 (Fig. 6j) and functionally annotated as cytokine and cytokine-cytokine receptor signaling genes, including IL6R and IL10 (Fig. 6k). As noted, VitD promoted IL-10 via IL-6 signaling through IL-6R and STAT3 (Figs. 3f-j). Despite higher basal mRNA, IL6 was still induced by VitD when BACH2 levels were sub-normal, but IL10 was not, signifying that the IL-6–STAT3–IL-10 axis was disrupted in BACH2 haploinsufficiency (Figs. 6i-l and Extended Data Fig.9e). Indeed, IL6R was not appropriately induced when BACH2 levels were sub-normal (Figs. 6i-l and Extended Data Fig.9e). In animals, the IL6r gene is a direct genomic target of Bach2 and Bach2 knockout status significantly impairs the expression of Il6r (Extended Data Fig.9f). Similar binding of this gene by BACH2 was evident in human cells (Fig. 6m), presumably explaining why BACH2 haploinsufficiency impaired IL6R expression. Since BACH2 also binds IL10 (Fig. 5h), it is likely that BACH2 deficiency also contributes directly to lower VitD-induced IL10 transcription. Collectively, these data indicate that BACH2 is a VitD-induced protein regulating a significant portion of the VitD transcriptome, with a key role in mechanisms that stimulate expression of IL-10.

VitD is predicted to retract the Th1 program in COVID-19

There is compelling epidemiological association between incidence and severity of COVID-19 and VitD deficiency/insufficiency[19], but the molecular mechanisms remain unknown. Given preferential Th1 polarization in COVID-19 BALF (Figs. 1d and Extended Data Fig.1e), we hypothesized that VitD could be mechanistically important for hyper-inflammation in COVID-19 and may be a therapeutic option. We studied expression of VitD-regulated genes (Fig. 3a and Table S2b) in COVID-19 BALF Th cells. Expression of VitD-repressed genes, summarized as module score, was significantly higher in patient Th cells than healthy controls (Fig. 7a). This was further corroborated by Geneset enrichment analysis (GSEA) showing that genes more highly expressed in patient compared to control cells were enriched in VitD-repressed genes (Fig. 7b). On a per cell basis the VitD-repressed module score correlated strongly with the Th1 score in BALF Th cells (Fig. 7c), indicating a reciprocal relationship between Th1 genes and VitD-repressed genes. In contrast, VitD-induced genes were not substantially different (Extended Data Fig.10a), nor were VitD-regulated genes different in PBMC of patients compared to healthy donors (Extended Data Fig.10b). To assess the ability of the Th1 and VitD-repressed gene signatures to distinguish patient from healthy BALF Th cells on a per cell basis, we constructed receiver operating characteristic (ROC) curves of the module scores. The area under the curve (AUC) for the two scores were 0.91 and 0.80, respectively (p<0.00001 for both) (Fig. 7d). Thus, both Th1 and VitD-repressed gene signatures were strong features of COVID-19 BALF Th cells. To contextualize these performances, we calculated AUCs for module scores of every geneset in hallmark and canonical pathways curated by MSigDB (n=2279 genesets). The top performing sets were interferon responses, as expected. The performance of the VitD-repressed geneset was within the top 1% of all genesets (Fig. 7e and Table S6), indicating that this was among the very best performing genesets. We then predicted significant drugs (among n=461) potentially able to counteract genes upregulated in COVID-19 BALF Th cells. Five of the top 10 most significant drugs were steroids and two were VitD (alfacalcidol) (Fig. 7f and Table S7). These were corroborated in independent datasets by GSEA, showing genes more highly expressed in lung biopsies or bulk BALF RNA-seq of COVID-19 compared to healthy control to be enriched in VitD-repressed genes (Figs. 7g-h). Significant numbers of genes modulated by VitD were shared by corticosteroid targets (Extended Data Fig.10c). Steroids, including dexamethasone, increase transcription of VDR so may show therapeutic synergism with VitD. Collectively, these data suggest that, in Th cells from patients with severe COVID-19, the VitD-repressed geneset is de-repressed and that there may be clinical benefit from adding VitD to other immunomodulatory agents (Extended Data Fig.10d).

Fig. 7

VitD is among the top predicted therapeutics to retract the Th1 program in SARS-CoV2-infected lungs.

a, Violin plots showing expressions of VitD-repressed genes, summarized as module scores, in BALF Th cells of patients with COVID-19 and healthy controls. Exact p-values in a have been calculated using two-tailed Wilcoxon tests. b, GSEA showing enrichment in VitD-repressed genes within genes more highly expressed in scRNAseq CD4+ BALF T cells of patients with COVID-19 compared to healthy controls. c, Correlation between module scores of Th1-genes and VitD-repressed genes on a per cell basis in BALF Th cells of patients with COVID-19 and healthy controls. Pearson r and exact p-values are shown. d, Receiver operating characteristic (ROC) curve, evaluating the performance of the Th1 and VitD-repressed module scores to distinguish BALF Th cells of patients with COVID-19 from healthy controls. Shown are the area under the curve (AUC) statistics and p-values. e, Analyses showing the performance of all MSigDB canonical and hallmark genesets to distinguish BALF Th cells of patients with COVID-19 from healthy controls, ranked by AUC values. Marked are the top 2 performing genesets in red and the position of the VitD-repressed geneset within the top 1% of all genesets. f, Top 10 drugs predicted (out of 461 significant drugs) to counteract genes induced in BALF Th cells of COVID-19 patients compared to healthy controls, ordered by adjusted p-value. Highlighted in red is alfacalcidol, an FDA-approved active form of VitD. g-h, GSEA showing enrichment in VitD-repressed genes for genes more highly expressed in bulk RNA-seq lung biopsy specimens (g) and bulk RNA-seq BALF cells (h) of COVID-19 compared to healthy controls. Empirical p-values are shown for GSEA in b, g-h; NES = normalized expression value. p-values in d are from the Mann-Whitney U-statistic. Data in a-f are from n=8 patients with COVID-19 and n=3 healthy controls, sourced from GSE145926 and GSE122960. Data in g-h are from GSE147507 and HRA000143, respectively, and n numbers are indicated.

Extended Data Fig. 10

Vitamin D-induced genes do not distinguish CD4+ BALF T cells of patients with COVID-19 from healthy controls.

a, Violin plots showing expressions of VitD-induced genes, summarized as module scores, of BALF CD4+ T cells of patients with COVID-19 and healthy controls. Data are from n=8 patients with COVID-19 and n=3 healthy controls, sourced from GSE145926 and GSE122960. b, Violin plots showing expressions of VitD-induced and VitD-repressed genes, summarized as module scores, of PBMC CD4+ T cells of patients with COVID-19 and healthy controls. Data are from n=6 patients with COVID-19 and n=6 healthy subjects, obtained from GSE150728. c, Venn diagram showing overlap between COVID-induced genes in CD4+ BALF T cells that are predicted to be normalized by VitD treatment versus those that are predicted to be normalized by steroid drugs. Please also see Table S7. d, Schematic model of autocrine VitD-driven Th1 contraction program and potential intervention of impaired COVID-19 program with VitD and cortico-steroids. e, representative flow sorting strategy and example of post-sort purity obtained for isolation of memory CD4+ T cells in this paper.

Discussion

We showed that cell-intrinsic complement orchestrates an autocrine/paracrine autoregulatory VitD loop to initiate Th1 shutdown. VitD causes genome-wide epigenetic re-modelling, induces and recruits TFs including STAT3, c-JUN and BACH2 that collectively repress Th1 and Th17 programs and induce IL-10 via IL-6-STAT3 signaling. This program is abnormal in lung Th cells of severe COVID-19, which show preferential Th1 skewing, and could be potentially exploited therapeutically by using VitD as an adjunct treatment.

IFN-γ-producing airway Th cells are key components of immunity to coronaviruses, including SARS-CoV1 and MERS-CoV[20]. Th1-polarized responses are also a feature of SARS-CoV2 in humans[24] and severe COVID-19 is accompanied by prolonged, exacerbated, circulating Th1 responses[21]. Complement receptor signaling is a driver of Th1 differentiation and required for effective anti-viral responses[12,42]. C3 cleavage generates C3b, which binds CD46 on T cells. We have previously shown that the lungs in COVID-19 are a complement-rich micro-environment, that local CD4+ T lymphocytes have a CD46-activated signature[13], and show here that these T cells are Th1-polarized. Pro-inflammatory function is important for pathogen clearance, but a switch into IL-10 production is a natural component during successful transition into the Th1 shutdown program and reduces collateral damage[16]. Inability to produce IL-10 results in more efficient clearance of infections but severe tissue damage from uncontrolled Th1 responses results in death[2]. The benefits of remediating inflammatory pathways in severe COVID-19 is demonstrated by successful trials of dexamethasone, an immunosuppressive drug which reduces mortality[1].

VitD has pleiotropic functions in the immune system, including anti-microbial as well as regulatory properties, which are cell- and context-dependent[43,44]. VitD deficiency is associated with higher prevalence and worse outcomes from infections, including influenza, tuberculosis and viral upper respiratory tract illnesses[17], as well as autoimmune diseases, including type 1 diabetes, multiple sclerosis, rheumatoid arthritis and inflammatory bowel disease[18]. Th cells play key roles in all these diseases. Thus, understanding VitD biology in Th cells has potential translational impact.

Our data indicate the complexities of TFs working within networks to regulate sets of genes. After VitD ligates VDR, c-JUN, STAT3 and BACH2 are recruited to acetylated loci, shaping the transcriptional response to VitD. c-JUN is a member of the AP-1 basic leucine zipper family, primarily involved in DNA transcription[45]. This TF was bound adjacent to 40% of VitD-regulated coding loci. BACH2 is a critical immunoregulatory TF[38,39] downregulated in lesional versus non-lesional psoriatic skin[46]. Active VitD concentrations inversely correlate with severity of psoriasis[40]. We found that VitD-treatment of psoriatic lesional skin upregulated BACH2 expression. Both haploinsufficiency and single nucleotide variants of BACH2 associate with monogenic and polygenic autoimmunity, respectively, in humans[36,47]. No BACH2 knock-out humans have yet been identified, suggesting incompatibility of complete BACH2 deficiency with life. Indeed, Bach2 mice succumb to fatal autoimmunity. We found that loss of even 50% of the normal cellular concentration of BACH2 in the haploinsufficient state substantially altered (~70% of) the VitD-regulated transcriptome. Since only a proportion of these genes were directly BACH2-bound, it is probable that BACH2 is a requisite for normal recruitment and function of the other transcriptional regulators.

Both the incidence and severity of COVID-19 are epidemiologically associated with VitD deficiency/insufficiency[19], but the molecular mechanisms remain unknown. We found a link between the inflammatory Th1 program and a VitD-repressed geneset. Attempts to study CD4+ T cells from the site of inflammation were unsuccessful due to the rapid apoptosis of patient cells, but our in silico analyses suggest either dysregulation of the VitD program in COVID-19 or that simple deficiency/insufficiency of substrate (VitD) might explain the epidemiological association.

IL-6 is a pleiotropic, often pro-inflammatory, cytokine. IL-6 is implicated in the COVID-19 “cytokine storm” and targeting of this cytokine specifically have proved beneficial to patients[32]. Our data suggest that pro-inflammatory IL-6 functions may be redirected to production of anti-inflammatory IL-10 by VitD in activated human Th cells. VitD supplementation in children significantly increases serum IL-6 (and non-significantly increases in IL-10)[48] indicating that these observations may also occur in vivo. In the skin, where VitD concentrations are high, IL-6 overexpression protects from injurious stimuli or infection[49] and IL-6-deficiency impairs wound healing[50]. Moreover, an adverse effect of anti-IL-6R for treating inflammatory arthritis is idiosyncratic development of psoriasis[51], indicating a tolerogenic role for IL-6 at this site. Thus, adjunct VitD therapy in severe COVID-19 may divert pro-inflammatory and induce anti-inflammatory effects of IL-6, which may be an alternative to blocking IL-6R signalling.

These data identified the VitD pathway as a potential mechanism to accelerate shutdown of Th1 cells in severe COVID-19. From experience in other diseases, it is likely that VitD will be ineffective as monotherapy. Combination therapy could potentially ameliorate significant adverse effects of other drugs, e.g. high dose corticosteroids, including over-immunosuppression or metabolic side-effects. An important consideration of VitD therapy in COVID-19 is stimulation of IL-6 production from CD4+ T cells. Although autocrine/paracrine IL-6 induces IL-10 in these cells, IL-6 could potentially have pro-inflammatory properties on other cells. These possibilities may be mitigated by adding VitD as an adjunct to other immunomodulators, such as corticosteroids or JAK inhibitors[13]. Of note, two randomized clinical trials with calcifediol, a VitD analogue with high bioavailability not requiring hepatic 25-hydroxylation, comprising >1000 patients together, reported reductions in risk of ICU admission or death when used in addition to standard care (OR of 0.13 and 0.22, respectively[52,53]). These findings are not necessarily specific to COVID-19, since VitD can protect against acute respiratory tract infections in general[17].

In conclusion, we identified an autocrine/paracrine VitD loop permitting Th1 cells to both activate and respond to VitD as part of a shut-down program repressing IFN-γ and enhancing IL-10. These events involved significant epigenetic re-shaping and recruitment of a network of key TFs. These pathways could potentially be exploited therapeutically to accelerate the shut-down program of hyper-inflammatory cells in patients with severe COVID-19.

Methods

Ethics approvals

Human studies, conducted in accordance with the Declaration of Helsinki, were approved by the Institutional Review Board of Guy’s Hospital (reference 09/H0707/86), National Institutes of Health (approval numbers 7458, PACI, 13-H-0065 and 00-I-0159) and Imperial College London (approval number 12/WA/0196 ICHTB HTA license number 12275 to project R14098). All patients provided informed written consent.

Human T cell isolation and culture

Human PBMCs were purified from anonymized leukodepletion cones (Blood Transfusion Service, NHS Blood and Transplantation) or healthy volunteer whole blood packs/buffy coats from the NIH Blood Bank and from fresh patient blood. Leukodepletion cones were diluted 1 in 4 with PBS and layered onto Lymphoprep (Axis-Shield). Whole blood packs/buffy coats were diluted 1 in 2 with PBS onto Lymphocyte Separation Medium (25-072-CV LSM, Corning) and centrifuged.

CD4+CD25− cells were used throughout, unless specified. CD4+ T cells were enriched using RosetteSep Human CD4+ enrichment cocktail (Stem Cell) with CD4+CD25− cells obtained by depletion of CD25+ T cells using CD25 positive selection (CD25 microbeads II, Human, Miltenyi Biotec) followed by unlabeled cell collection. Human memory CD4+ T cells were used for RNA-seq, CUT&RUN and CUT&Tag experiments as these cells express the VDR. These cells were isolated from PBMCs using Miltenyi Memory CD4+ T cell Isolation Kit Human (130091893) or StemCell EasySep Human Memory CD4+ T Cell Enrichment Kit (19157) or by flow sorting.

For flow sorting, bulk CD4+ T cells were enriched as above then stained with antibodies against CD4 (OKT4, Thermo Fisher Scientific), CD45RA (HI100, Biolegend), CD45RO (UCLH1, BD Biosciences) and CD25 (2A3, BD Biosciences) in MACS buffer at 4C for 30 min. CD4+CD25−CD45RO+CD45RA− memory cells were sorted into media to a >99% purity by a FACSAria (BD Biosciences). Sorting strategy and representative post-sort purities are shown in Fig S10e.

Cells were cultured in X-VIVO-15 Serum-free Hematopoietic Cell Medium (04-418Q, Lonza) supplemented with 50 IU/mL penicillin, 50μg/mL streptomycin, and 2mM L-glutamine (Gibco), at 37C 5% CO2 in a non-tissue culture treated 96 well U bottom plate (Greiner) at a density of 106 cells per mL in 200μL. Cells were activated using Human T-Activator CD3/CD28 Dynabeads (11131D, Gibco, cell to bead ratio of 4:1). Where indicated, cells were plate-activated in non-tissue culture treated 48 well plates (Greiner) pre-coated with anti-CD3 (OKT3, from Washington University hybridoma facility), anti-CD3+anti-CD28 (CD28.2, Becton Dickinson) or anti-CD3+anti-CD46 (TRA-2-10, a gift from Dr John P Atkinson, Washington University, USA), with or without 1μM Cathepsin L inhibitor (ALX-260-133-M001, Enzo Life Sciences), with IL-2 included (50 U/mL; PeproTech), in a total volume of 250μL at a density of 5x105/well. Antibody coated plates were prepared by diluting antibodies (2μg/mL) in D-PBS+Ca2+Mg2+ (Gibco) and adding to 48-well plates overnight in a volume of 150μL/well at 4C, subsequently the plates were rinsed with PBS before cell suspension addition and centrifugation (swinging-bucket rotor, 300g, 3min, RT). Cells were additionally cultured with 1α,25-Dihydroxyvitamin D3 (BML-DM200-0050, Enzo Life Sciences) or 25-Dihydroxyvitamin D3 (BML-DM100-0001, Enzo Life Sciences), both reconstituted in 99.8% ethanol (Sigma-Aldrich), used at 10nM unless indicated, with identically diluted ethanol used as carrier control. IL-6 (Biolegend) and Tocilizumab (a gift from Dr Ceri Roberts in Professor Leonie Taams’ lab) were used in functional experiments.

Flow cytometry

Cells were stained 30 min at 4C in a final volume of 100μL staining buffer with antibodies against CD4 (OKT4, Thermo Fisher Scientific), IL-6 (MQ2-13A5, BioLegend), LAG-3 (3DS223H, Thermo Fisher Scientific), CD49b (P1H5, Thermo Fisher Scientific) with DAPI or LIVE/DEAD Fixable Aqua/Violet (Invitrogen) to exclude dead cells. For intracellular staining, cells were treated with the spiked addition of Phorbol 12-myristate 13-acetate (PMA) (50ng/mL, Sigma-Aldrich), ionomycin (1μg/mL, Sigma-Aldrich), GolgiStop (1X, BD Biosciences) and Brefeldin A (1X, BD Biosciences) to media then cultured 5h, subsequently washed in PBS, surface stained, treated with Cytofix/Cytoperm (eBioscience) and subsequent Perm/Wash buffer (eBioscience) steps before incubation with intracellular antibodies for 30 min at 4C. Samples were washed and acquired on a LSRFortessa (BD Biosciences) or the Attune NxT Flow Cytometer (Invitrogen) within 24h and analyzed using FlowJo/v9.9.6. Cell proliferation analysis was carried out by CFSE dilution in polyclonally activated CD4+CD25− T cells with carrier or VitD after 72h.

Cell proliferation and activation markers were assessed in purified CD4+ cells with CellTrace Violet (CTV, C34556, Invitrogen) and surface staining against CD62L-FITC (MEL-14, BioLegend), CD69-PE-Cy7 (FN50, Invitrogen) or CD25-PE-Cy7 (BC96, Invitrogen) and LIVE/DEAD fixable near-IR (Invitrogen).

Cytokine measurement

Supernatants from 96-well plates were aliquoted and stored at -20C. Cytokines were quantified by the Human/Mouse Th1/Th2/Th17 Cytometric Bead Array (CBA) (BD Biosciences) using a FACSCanto II (BD Biosciences), analyzed by FCAP Array v3.0 (BD Biosciences). Alternatively, we used the LEGENDplex Human Inflammation Panel 1 (13plex) (740808, BioLegend) with the Attune NxT Flow Cytometer. Cytokine assays analyzed by FlowJo/v9.

For intracellular cytokine detection, the IFN-γ/IL-10 Secretion Assay Detection Kit (Miltenyi, #130-090-761) was used on live cells after 72h.

Immunoblotting

Cell extracts were prepared by lysis in RIPA (Thermo Fisher Scientific) including 1 in 250 diluted Protease Inhibitor Cocktail Set III (Calbiochem) and 5% 2-Mercaptoethanol, denatured at 95C for 5 min. Protein concentration was quantified using the Quick Start Bradford Protein Assay (BioRad). Proteins resolved by SDS-PAGE on 10% Tris-Glycine (Invitrogen) were transferred onto polyvinylidene fluoride membranes (Millipore) via a XCell II Blot Module (Thermo Fisher Scientific). Immunoblotting was performed with blocking in PBS+10% w/v Blotting-Grade Blocker (1706404, BioRad) for VDR, STAT3 and c-Jun or 3% w/v BSA (Sigma-Aldrich) for pSTAT3 and Phospho-c-Jun made up in PBS+0.1% v/v Tween20 (Sigma-Aldrich) then incubated overnight with antibodies against VDR (1:100, D-6, Santa Cruz), pSTAT3 (1:2000, D3A7, Cell Signaling Technologies), STAT3 (1:1000, 124H6, Cell Signaling Technologies), Phospho-c-Jun (1:1000, S63, R&D systems) or c-Jun (1:2000, L70B11, Cell Signaling Technologies), followed by 2h incubation with HRP conjugated anti-mouse or anti-rabbit antibodies (1:1000, TrueBlot, Rockland), then washed with PBS+0.5% v/v Tween20. Targets were detected using ECL-Plus (Thermo Fisher Scientific), visualized by ImageQuant LAS 4000 mini (GE Healthcare), and quantified (Image Studio Lite, LI-COR).

Nuclear and cytoplasmic extracts were isolated (NE-PER Nuclear and Cytoplasmic Extraction, Thermo Fisher Scientific) and purity determined by immunoblotting against HSP90 (1:2000, C45G5, Cell Signaling Technologies) and H2A.X (1:1000, D17A3, Cell Signaling Technologies).

For BACH2, snap-frozen cell pellets of 6x105 cells were resuspended directly in Laemmli+5% 2-Mercaptoethanol (Bio-Rad) and denatured (95C, 5 min), RT cooled, loaded and resolved at 120V on a Criterion TGX Gel (Bio-Rad). Gel-bound proteins were transferred to nitrocellulose (Trans-Blot Turbo Midi, Bio-Rad) via a Trans-Blot Turbo (Bio-Rad). Membrane blocking occurred (PBS Odyssey Blocking Buffer, LI-COR) prior to the addition of BACH2 (D3T3G, Cell Signaling Technologies) antibody in blocking buffer (1:5000) then rotated overnight, 4C. Membranes washed thrice (PBS+0.1% Tween-20) were exposed to IRDye 800CW Goat anti-Rabbit IgG Secondary Antibody (926-32211, LI-COR) in blocking buffer (1:5000) for 1h, RT. After three washes, antibody-signals were visualized (Odyssey CLx, LI-COR). The membrane was re-probed with anti-β-Tubulin diluted in blocking buffer (2146, Cell Signaling Technologies, 1:1000 dilution) and quantified (Image Studio v5.2, LI-COR).

Human phospho-kinase antibody array

Proteome Profiler Human Phospho-Kinase Array (R&D systems) was carried out on VitD- or carrier-treated CD4+ cells. Cells were lyzed in manufacturer’s buffer, protein concentrations quantified by Quick Start Bradford Protein Assay (BioRad) and samples adjusted to 800μg/mL with lysis buffer. 334μL were loaded per membrane and signals amplified by Chemi-Reagent Mix imaged (ImageQuant LAS 4000 mini) and quantified (Image Studio Lite). Reference spots were used to normalize signals across membranes.

VDR and DAPI co-localization

VDR and DAPI colocalization was performed in CD4+ T cells stained with VDR antibody (1:250, D-6), anti-mouse-AlexaFluor647 (1:400, clone poly4053, BioLegend) and DAPI. Images were acquired via ImageStreamX running Inspire and VDR to DAPI colocalization assessed (IDEAS/v3.0, Amnis). A minimum of 1x105 events were acquired per sample.

Quantitative PCR (qPCR)

2-4x106 cells were lysed in 350μl Trizol (Life Technologies) and RNA extracted by Direct-zol Miniprep (Zymo Research) with on-column genomic DNA digest. RNA was quantified by spectrophotometry (NanoDrop) and reverse transcribed (qPCRBIO cDNA Synthesis, PCR Biosystems). qPCR in 384-well plates with SYBR Green PCR Master Mix (Life Technologies) was used with the ViiA7 real-time PCR system (Life Technologies). Reactions carried out in triplicate with 18s and UBC (geNorm 6 gene kit, ge-SY-6, PrimerDesign UK) as reference alongside IL6 detection (Hs_IL6_1_SG QuantiTect Primer Assay, Qiagen). TaqMan probes (Applied Biosystems) for CYP27B1 (Hs01096154_m1), VDR (Hs01045843_m1), IFNG (Hs00989291_m1), IL10 (Hs00961622_m1), CTSL (Hs00964650_m1), CYP27B1 (Hs01096154_m1), and 18S (Hs99999901_s1) were utilized on reverse-transcribed RNA (iScript, Bio-Rad) with TaqMan Master Mix (Thermo Fisher Scientific). Comparative Ct method was used for analysis using the Viia7 software producing a Ct (relative quantity (RQ) compared to a reference sample).

Inhibition and siRNAs

For STAT3 siRNA, memory CD4+ T cells were sorted and rested overnight in culture medium. 5x106 cells were washed twice with PBS and nucleofected with 10nM STAT3 silencer select siRNA (s745, sense GCACCUUCCUGCUAAGAUUTT and antisense AAUCUUAGCAGGAAGGUGCCT) or Silencer Select Negative Control (Thermo Fisher Scientific) using the Amaxa Human T Cell Nucleofector Kit (Nucleofector 2b Device, U-014 program, Lonza). Nucleofected cells were cultured in media for 6h followed by dead cell removal (Lymphoprep, Axis-Shield). 2x106 recovered cells were activated with anti-CD3/CD28 beads for 72h and knock-down assessed by western.

For chemical STAT3 inhibition, memory CD4+ T cells were sorted and rested overnight. 2x106 cells were activated as before with VitD or carrier and additionally with 150 nM Curcubitacin I, Cucumis sativus L. (Calbiochem) or carrier (ethanol) for 72h.

For VDR, CYP27B1 and CTSL siRNA, total CD4+ T cells were plate-activated with CD3/CD46 antibodies as described above. Dharmacon Accell SMARTPool siRNAs against CYP27B1 (E-009757-00-0005), VDR (E-003448-00-0005), CTSL (E-005841-00-0005), Non-targeting Control (D-001910-01-05) or Green Non-targeting Control (D-001950-01-05) were added to cultures at 1μM for 72h with knock-down assessed by qPCR.

H3K27Ac and c-JUN CUT&RUN

CUT&RUN sequencing was performed on memory CD4+ T cells activated as above in the presence of carrier or VitD for 0.75h, 48h and 72h using the published protocol of Skene et al.[55]. Antibodies against H3K27Ac (ab4729, Abcam), c-JUN (60AB, Cell Signaling Technologies), and non-specific IgG (31235, Thermo Fisher Scientific) were used with pAG-MNase (123461, Addgene) on 5 x 104 cells/target. CUT&RUN (and CUT&Tag below) buffer components were purchased from Sigma-Aldrich and Thermo Fisher Scientific. Post-CUT&RUN, short DNA fragments were prepared for paired-end sequencing (NEB Ultra II, New England Biolabs) with NovaSeq (Illumina). Amplified libraries were quantified by high sensitivity fluorometry (DeNovix) and sized via the 4200 TapeStation (Agilent Technologies).

H3K27Ac reads aligned to the human reference genome (GRCh37; hg19) using Bowtie2[56] with parameters ‘--local --maxins 250’ with sorting and indexing of aligned reads by samtools/1.10[57], reads with mapping quality <30 removed. H3K27Ac peaks were called by MACS2[58] with parameter ‘-f BAMPE --nomodel’. JUN reads were mapped without ‘--local’ parameter. <120 bp fragments were further sorted and indexed by samtools/1.10[57]. JUN induced peaks called by MACS2 parameters ‘--nomodel -g hs -f BAMPE -q 0.01 --SPMR --keep-dup all’. IgG BAM files were included as controls for MACS2. Peaks for H3K27Ac and JUN experiments were screened against previously characterized hypersequencable-regions (Hg19 Blacklist, CutRunTools[59]). H3K27ac and JUN mean binding intensities in peak regions calculated by UCSC bigWigAverageOverBed. Peaks, combined from two conditions, with mean binding intensity greater than 0.2 in at least one condition and fold change of signal>1.5 were differential. Known motifs +/−100bp from differential peak summits were identified (HOMER, findMotifsGenome.pl, parameter ‘-size given’). De novo motif footprinting for c-JUN was performed (CutRunTools[74]). Average cuts per bp proximal to JUN peaks was plotted as a histogram relative to IgG. CUT&RUN tracks and heatmaps were visualized using IGV (Broad Institute) and deepTools[60], respectively.

De novo super-enhancer calling

Super-Enhancers were called using ROSE (Rank Ordering of Super-Enhancers)[37]. H3K27Ac CUT&RUN VitD induced peaks were sorted/indexed (samtools/1.10[57]) and executed with ROSE_main.py with carrier control. SE annotations were screened and corrected, if needed. Hockey-stick plots of rank ordered stitched-enhancers plotted against VitD minus carrier H3K27Ac signal were generated (RStudio/v1.2.5033)

VDR, STAT3, and BACH2 CUT&Tag

VDR, STAT3, and BACH2 genome-wide binding in memory CD4+ T cells activated in the presence of 10nM VitD or carrier for 96h was detected by CUT&Tag as per the published protocol (CUT&Tag V.3, dx.doi.org/10.17504/protocols.io.bcuhiwt6)[61] using fixed nuclei (2min fixation of isolated nuclei, RT) with 75,000 nuclei/target and antibodies targeting either VDR (D2K6W, Cell Signaling Technologies), STAT3 (D3Z2G, Cell Signaling Technologies), BACH2 (D3T3G), or non-specific IgG (31235) diluted 1:50 accordingly and incubated with fixed-nuclei for 1h, RT. Secondary antibody (1:100) (ABIN101961, antibodies-online) was bound for 0.5h, RT. Preloaded pA-Tn5 (C01070001, Diagenode) was tethered to antibodies for 1h, RT. Libraries were sequenced paired-end by Novaseq. Resulting reads were processed and mapped to Hg19 (CutRunTools), with peaks called (SEACR/v1.2, stringent, all mapped fragments, https://github.com/FredHutch/SEACR) and hypersequencable-regions subtracted. Known motif enrichment was determined +/−200bp from peak summits and data visualized as in CUT&RUN.

VDR ChIP-PCR

ChIP-PCR detection of VDR binding at the promoters of CYP24A1, STAT3, and IL-6 was determined from 2x105 CD4+ cells activated with VitD or Carrier for 48h followed by True MicroChIP (C01010132, Diagenode, 3 bioruptions/sample). 4μg of anti-VDR (D-6) or IgG2a isotype control antibody (E5Y6Q, Cell Signaling Technologies) was used during overnight chromatin precipitation at 4C before reverse-crosslinking and purification (MicroChIP DiaPure, Diagenode). qPCR primers were designed to amplify regions within 500bp upstream of the transcriptional start site for each respective gene: CYP24A1 (F: TGACCGGGGCTATGTTCG, R: GGCTTCGCATGACTTCCTG); STAT3 (F: CTGTTCCGACAGTTCGGTGC, R: GCAGGACATTCCGGTCATCTTC); and IL-6 (F: GTAAAACTTCGTGCATGACTTCAGC, R: GGGGGAAAAGTGCAGCTTAGG). Triplicate qPCR using POWRUP SYBR Green (Thermo Fisher) was quantified via CFX384 (Bio-Rad), and results calculated as %Input= 2(-ΔCt [normalized ChIP]).

RNA sequencing and analysis

Memory CD4+ T cells, including from BACH2WT/L24P patients, were activated with VitD or carrier. 6x105 cells were pelleted at 300g for 5 min and RNA extracted by RNAqeous Micro (Thermo Fisher Scientific). 1μg total RNA for each sample was subjected to NEBNext Poly(A) mRNA Magnetic Isolation (E7490) and resulting mRNA prepared for RNA-seq by NEB Ultra II (New England Biolabs) and sequenced by HiSeq (Illumina).

The expression levels of all genes in RNA-seq libraries were quantified by ‘rsem-calculate-expression’ in RSEM/v1.3.1[62] with parameters ‘--bowtie-n 1 --bowtie-m 100 --seed-length 28 --bowtie-chunkmbs 1000’. The bowtie index for RSEM alignment was generated by ‘rsem-prepare-reference’ on all RefSeq genes, downloaded from UCSC table browser in April 2017. EdgeR/v3.26.8[63] was used to normalize gene expression among all libraries and identify DEGs among samples. Microarray analysis sourced from GSE119416 was carried out using Partek Genomics Suite (Partek, Inc.). DEGs were defined using the following criteria: at least 1.5-fold change in either direction at p-value<0.05 for microarray; at least 1.75-fold change in either direction at FDR<0.05 for RNA-seq.

Single cell RNA-sequencing (scRNAseq) and analysis

COVID-19 bronchoalveolar lavage dataset

The pre-processed h5 matrix files for nine COVID-19 patient bronchoalveolar-lavage (BAL) samples and four healthy control BAL samples were from GSE145926[64] and GSE122960[65], and GSM3660650, respectively. Read mapping and basic filtering performed with Cell Ranger (10x Genomics). Further processing using Seurat/v3[66] as follows: Only genes found to be expressed in more than 3 cells and cells with gene number between 200 and 6,000 (UMI count >1,000) retained. Cells with >10% of their unique molecular identifiers (UMIs) mapping to mitochondrial genes or cells with <300 features discarded to eliminate low quality cells. Filtering yielded 66453 cells across 12 samples. Filtered count matrices were normalized by total UMI counts, multiplied by 10,000 and transformed to natural log space. Top 2000 variable features determined by the variance stabilizing transformation function (FindVariableFeatures) using default parameters. All integrated using canonical correlation analysis (FindIntegrationAnchors/IntegrateData) with parameter k.filter=140 [67]. Variants arising from library size and percentage of mitochondrial genes regressed out (ScaleData, Seurat). Principal Component Analysis (PCA) was performed, and the top 50 Principal Components (PCs) were included in a Uniform Manifold Approximation and Projection (UMAP). Clusters identified on a shared nearest neighbor (SNN) modularity graph using the top 50 PCs and the original Louvain algorithm. Cluster annotations are based on canonical marker genes.

Cells identified as “T” and “Cytotoxic T” were subsetted and reprocessed. Samples C146 and C52 were removed due to low T cells (19 and 56 cells, respectively). Integration was performed using FindIntegrationAnchors/IntegrateData with parameters: “dims=1:30 k.filter=124”. Normalization, variable feature detection, scaling, dimensionality reduction, and clustering performed using the top 30 PCs and clustering resolution of 0.3. Annotation guided by marker genes (Extended Data Fig.1c). The contaminating macrophage cluster annotated by canonical marker genes (e.g. CD14) was removed. Statistical differences were assessed by two-tailed Wilcoxon test. DEGs were defined by fold change of >1.5, adjusted p value of <0.05 and expression in at least 10% of cells in either cluster comparison.

Two clusters annotated as CD4+ T cells, one with few cells but expressed similar markers to the main CD4+ cluster. The two were combined and considered together for analysis. DEGs between COVID-19 and Healthy CD4+ T cells were visualized using Morpheus (https://software.broadinstitute.org/morpheus/). DEGs upregulated in COVID-19 CD4+ T cells were subjected to pathway analysis (Fig. 1c) using the hallmark geneset collection from the Molecular Signatures Database (MSigDB/v7.1) or drug prediction from the National Toxicology Program’s DrugMatrix genesets (n=7876) using EnrichR (Fig. 4f). Gene list scores calculated by the AddModuleScore function in Seurat with a control gene set (n of 100). The correlation and the ROC analyses on the module scores comparing patients and healthy controls performed in R or Prism/v8.4.0.

COVID-19 PBMC dataset

The pre-processed R objects for six COVID-19 patient PBMC samples and six healthy control PBMC samples were from GSE150728[68] with read mapping and filtering by Cell Ranger. The exon count matrices were further processed by Seurat as follows: Only genes found to be expressed in more than 10 cells were retained. QC steps for filtering the samples as described[68]. Cells with 1,000-15,000 UMIs and <20% of reads from mitochondrial genes were retained. Cells with >20% of reads mapped to RNA18S5 or RNA28S5, and/or expressed >75 genes per 100 UMIs excluded. SCTransform function normalized the dataset and identified variable genes[68]. PCA was performed and top 50 PCs were included in UMAP dimensionality reduction. Clusters were identified on a SNN modularity graph using the top 50 PCs and the original Louvain algorithm. Cluster annotations were based on canonical marker genes. Gene list module scores were calculated by the AddModuleScore function in Seurat[66]. AUCs for module scores of every geneset in hallmark and canonical pathways curated by MSigDB (n=2279 genesets)[69-71].

Geneset Enrichment Analysis (GSEA) and gene lists

GSEA was carried out as published[69]. Tr1 gene lists were obtained from GSE139990. Data were filtered to include only genes expressed at greater than 0.25 counts per million in at least 2 samples, TMM normalized within the edgeR/v3.28.1, and differential expression performed using the glmQLFit and glmQLFTest edgeR functions. Tr1 signatures defined as DE genes (at least 4-fold change in either direction at FDR<0.05) between Tr1 and Th0 cells. Th1/Th2/Th17 gene lists obtained from transcriptomes of sorted mouse CD4+ T cell subsets[72]. VitD-regulated genes defined as described in RNA-seq methods. BACH2 bound and regulated genes obtained from GSE45975[38]. All gene lists used in this manuscript are in Table S8.

Psoriasis lesional skin staining, imaging, and analysis

Informed consent for all psoriasis patients and data was obtained under the Psoriasis, Atherosclerosis and Cardiometabolic Disease Initiative (PACI, 13-H-0065). Samples included 3 psoriasis patients not on VitD and 2 psoriasis patients on VitD supplementation. 3 mm punch-biopsies were collected from psoriasis lesional skin and placed into 10% formalin. 7μm skin sections were mounted on glass slides, de-waxed, and antigens-retrieved via the citrate buffer method. Sections washed once in PBS and blocked 20 min at RT in 10% normal goat serum then stained with anti-CD3 (1:50, F7.2.38, Thermo Fisher Scientific) and anti-BACH2 (1:40, D3TG) primary antibodies overnight at 4C. The next day washed in 0.01M PBS and Alexa Fluor 488 and Alexa Fluor 594 (1:100, Jackson Immuno Research Laboratories Inc., 115-546-062 and 111-586-047) secondary antibodies added 1h at RT and then washed in 0.01M PBS. Hoechst was added to the samples for 10 min, then washed and cover-slipped with Fluoromount G. Images were collected on a Zeiss 780 inverted confocal microscope at 40X with oil-immersion and analyzed (ImageJ). Representative images prepared in Zen Blue 3.1 (Zeiss). For each sample, 5-6 images were acquired. To quantify the frequency of immune cells in the dermis, images were cropped to 1200 x 1200 pixels to exclude the epidermis and nuclei numbers in the 408 nm channel determined. The number of BACH2 positive cells in the 594 nm channel were calculated per image and the number of BACH2 foci per BACH2 positive cell recorded.

Data presentation and statistical analysis

Figures prepared using Adobe Illustrator (Adobe). Statistical analysis and graphical visualizations carried out in GraphPad Prism (v.8.4.0), XLstat biomed (v2017.4), DataGraph 4.5.1 (Visual Data Tools, Inc.), Cytoscape 3[73] and Circos Table Viewer v.90.63.9[74]. Statistical analyses performed using appropriate paired or unpaired parametric and non-parametric tests, as required. Multiple comparisons performed using ANOVA. p-values <0.05 considered statistically significant throughout.