Vitamin D-regulated Gene Expression Profiles: Species-specificity and Cell-specific Effects on Metabolism and Immunity.

Vitamin D has pleiotropic physiological actions including immune system regulation, in addition to its classical role in calcium homeostasis. Hormonal 1,25-dihydroxyvitamin D (1,25D) signals through the nuclear vitamin D receptor, and large-scale expression profiling has provided numerous insights into its diverse physiological roles. To obtain a comprehensive picture of vitamin D signaling, we analyzed raw data from 94 (80 human, 14 mouse) expression profiles of genes regulated by 1,25D or its analogs. This identified several thousand distinct genes directly or indirectly up- or downregulated in a highly cell-specific manner in human cells using a 1.5-fold cut-off. There was significant overlap of biological processes regulated in human and mouse but minimal intersection between genes regulated in each species. Disease ontology clustering confirmed roles for 1,25D in immune homeostasis in several human cell types, and analysis of canonical pathways revealed novel and cell-specific roles of vitamin D in innate immunity. This included cell-specific regulation of several components of Nucleotide-binding Oligomerization Domain-like (NOD-like) pattern recognition receptor signaling, and metabolic events controlling innate immune responses. Notably, 1,25D selectively enhanced catabolism of branched-chain amino acids (BCAAs) in monocytic cells. BCAA levels regulate the major metabolic kinase mammalian Target of Rapamycin (mTOR), and pretreatment with 1,25D suppressed BCAA-dependent activation of mTOR signaling. Furthermore, ablation of BCAT1 expression in monocytic cells blocked 1,25D-induced increases in autophagy marker LAMP1. In conclusion, the data generated represents a powerful tool to further understand the diverse physiological roles of vitamin D signaling and provides multiple insights into mechanisms of innate immune regulation by 1,25D.

Initially identified as the curative agent for nutritional rickets, vitamin D is now known to exert pleiotropic physiological effects (1, 2). Most diets are relatively poor in vitamin D, but it can be produced by photochemical and thermal conversion in skin from the cholesterol precursor 7-dehydrocholesterol in the presence of adequate solar ultraviolet B irradiation. However, in regions outside of the tropics, ultraviolet B radiation is absorbed by the ozone layer during fall and winter, and cutaneous vitamin D synthesis does not occur (vitamin D winter). This can last for 6 months or more at higher latitudes, including many major population centers of Europe and North America (3). Coupled with sun avoidance and, in some regions, conservative dress, this leads to widespread vitamin D deficiency (4). It is therefore important to fully understand the diverse physiological roles of vitamin D and their potential clinical relevance.

Cutaneous or dietary vitamin D enters the blood stream bound to the vitamin D binding protein (DBP), a Gc globulin (5, 6), and is converted in a largely hepatic reaction to the major metabolic form 25-hydroxyvitamin D3 (25D). Subsequent 1α-hydroxylation by CYP27B1 in kidneys and peripheral tissues produces the hormonally active form, 1,25-dihydroxyvitamin D3 (1,25D, calcitriol) (7, 8). While renal CYP27B1 transcription is controlled by calcium homeostatic signals, its extrarenal expression is regulated by distinct physiological inputs and contributes to autocrine/paracrine vitamin D signaling in peripheral tissues (9, 10). Calcitriol binds the vitamin D receptor (VDR), a ligand-activated member of the nuclear receptor family of transcription factors (11-13). 1,25D binding induces VDR heterodimerization with related retinoid X receptors and transcriptional activation by association with specific vitamin D response elements (VDREs) (14, 15), composed of direct 5′-PuGG/TTCA-3′ repeats separated by 3 bp (DR3 elements) (16). Individual gene expression profiling studies have revealed that 1,25D represses expression of approximately 50% of the genes it regulates, and have provided evidence for heterogeneous mechanisms of repression, some of which require interaction with other classes of transcription factors (17).

In line with the pleiotropic actions of vitamin D, the VDR and CYP27B1 are widely expressed, enabling local signaling by 1,25D. Notably, both are present in many cells implicated in innate and adaptive immune responses, consistent with a role of 1,25D in immunity (18-21). CYP27B1 expression is regulated by a complex cytokine network in immune cells (22, 23), and 1,25D boosts antimicrobial innate immune responses in several cell types (19, 21, 24-27) by diverse mechanisms, including induction of expression of antimicrobial peptides such as cathelicidin (CAMP). Vitamin D signaling also inhibits intracellular Mycobacterium tuberculosis growth and strongly enhances infection-induced interleukin (IL)-1β production in human macrophages (27, 28). Induction of IL1B expression by 1,25D is species specific and does not occur in mouse. Indeed, regulation by 1,25D of several genes implicated in human innate immune responses, including antimicrobial peptides, appears to be human/primate specific and is not conserved in rodents (29, 30). Notably, Gombart et al. found that the VDRE in the CAMP gene is embedded in a human/primate-specific Alu repeat transposable element, whose insertion predated the old world–new world primate split (31). Nonetheless, vitamin D signaling does regulate immunity in mice; it has beneficial effects in mouse models of allergies, and autoimmune and inflammatory disorders, including inflammatory bowel disease (32-38), the multiple sclerosis model experimental autoimmune encephalomyelitis (39, 40), diabetes (41-43), systemic lupus erythematosus (44), rheumatoid arthritis (45), and asthma (46).

Although there is growing clinical and epidemiological evidence for the nonclassical physiological roles of vitamin D, in particular its role in immunity, much remains to be learned about the molecular mechanisms underlying these effects (47). As 1,25D signals through the nuclear VDR to directly regulate gene transcription, we performed a large-scale analysis of human and mouse microarray or RNA-seq gene expression profiling studies to understand in depth the tissue- and species-specific mechanisms of vitamin D action. Only datasets with available raw data were included, leading to analysis of 80 human and 14 mouse studies. The results provide numerous insights into the diverse physiological roles of vitamin D signaling, as well as its tissue and species specificity. They show that while there is substantial overlap in the biological processes regulated by 1,25D in human and mouse, the similarities in specific gene-regulatory events are much more limited. Further analysis of the data revealed novel pathways of immune regulation by 1,25D in human cells. These include the cell-specific catabolism of branched-chain amino acids (BCAAs) in macrophages, whose levels serve as signals of nutritional status and control metabolism and activation. In conclusion, the data generated offer several insights into vitamin D action and will provide valuable resources for future mechanistic studies of 1,25D signaling.

Material and Methods

Data collection

The Gene Expression Omnibus public repository was searched for studies using the keyword “vitamin D.” All gene expression studies examining the effects of vitamin D and its metabolites or calcitriol analogs, for which the raw data were available, were included. Studies in VDR knockout mice were excluded from the analysis.

Gene expression analysis

All gene expression analyses were performed using the R statistical package. To get signal intensities from CEL files of Affymetrix microarrays, the oligo package was used and data was normalized and summarized using the Robust Multi-Array Average method, part of the oligo package. Illumina raw input was summarized using Illumina’s BeadStudio or GenomStudio. Illumina, Agilent, Nimblegen, and custom arrays were normalized using the Linear Models for Microarray Data (LIMMA) package in R. RNA-seq raw files (SRA) were downloaded using Aspera and converted to reads in fastq format by the fastq-dump function in the sratoolkit suite from NCBI. A genome index was built based on human/mouse genomes, provided by Calcul Quebec high-performance computing cluster at McGill University, Guillimin, using the Rsubread package. The latter was also employed to align RNA-seq reads in order to obtain read counts based on the “seed-and-vote” paradigm. The resulting data was formatted for downstream analysis using functionalities in the package EdgeR. RNA-seq data manipulation was done on Guillimin. Differential gene expression for microarray and RNA-seq data was assessed using the LIMMA package. LIMMA was chosen for RNA-seq studies as well, rather than the commonly used Cuffdiff method from the Tuxedo suite. The latter appears to perform much worse than the former, particularly in situations featuring small number of replicates (48). Annotation including human orthologs for mouse genes was performed via the biomaRt package, an interface to the BioMart databases at Ensembl.

Quality control criteria

For each dataset, total number of RNA-seq reads sequenced for each condition was similar and so was total GC content, as expected. The Phred offset quality score was 33 and only the single best alignment was kept for downstream analysis. The minimum fragment size for alignment was set to 50. These parameters resulted in 70% to 85% of reads aligning unambiguously to a gene using the RSubread (49) and EdgeR (50) packages in R. Genes with very low read counts were filtered out and only those with counts ≥10 for at least 1 treatment group (in all replicates) were included for downstream analysis. Similarity between replicates in each treatment group was confirmed using a principal components analysis (Fig. S1 (51)). For Illumina microarrays, probe summary files (including control probes) were exported from Illumina’s GenomeStudio (52) without background correction or normalization. For the rest of the microarray platforms, raw files were analyzed directly in R using the affy (53) and limma (54) packages. Signal distribution of the raw and normalized data was visualized using box plots in order to identify potential outliers (Fig. S2 (51)). For all datasets, samples appeared comparable after background correction and normalization.

Enrichment analyses

Enrichment analysis and visualization was performed using the clusterProfiler package (55). Lists of significantly, differentially regulated genes (fold change ≥1.5 and P ≤ .1) for each dataset were generated as the basis for further analysis. Enrichment for disease ontology and for canonical pathways (Kyoto Encyclopedia of Genes and Genomes) terms based on the hypergeometric test was performed using clusterProfiler, which also calls the package DOSE. Enrichment analysis for gene ontology (GO) terms, on the other hand, relied on Fisher’s exact test and the topGO (56) package. Note that the gene universe was taken to be all known genes rather than expressed/measured genes for each dataset. Graphical representation of the enrichment results was done with the aid of clusterProfiler. The package pathview (57) was used for data integration and visualization of the gene expression results. Although the change in expression of certain genes displayed in the Kyoto Encyclopedia of Genes and Genomes pathway is more than 3-fold, the fold change was manually set to 3 for clarity.

Selection and representation of over-represented and overlapping genes

Lists of differentially expressed genes from same-cell–type datasets were combined and filtered to exclude genes that appear as both up- and downregulated. In order to quantify the overlap between human and mouse over-represented genes, mouse differentially expressed genes were annotated with the names and gene IDs for their human orthologues. The overlap was illustrated using Venn diagrams implemented by the VennDiagram (58) package in R. The distribution of the number of appearances of each over-represented gene in same-cell–type datasets is shown by violin plots (displaying the probability density at each value) through the vioplot package in R.

Reagents

1α,25-Dihydroxyvitamin D3 (BML-DM200) was purchased from Enzo Life Sciences and was used at a final concentration of 100 nM. BCAAs were purchased from Sigma Aldrich and used at the following concentrations: L-leucine 115 µM (Cat# L8000), L-valine 200 µM (Cat# V0513), and L-isoleucine 60 µM (Cat# I7403).

Cell culture

SCC25 cells (CRL1628; ATCC) were cultured in Dulbecco’s modified Eagle’s medium/F12 (319-085-CL; Wisent) supplemented with 10% fetal bovine serum (FBS). THP-1 cells were cultured in RPMI-1640 (350-006-CL; Wisent) supplemented with 10% FBS and differentiated into macrophage-like cells by adding 10 nM phorbol 12-myristate 13-acetate for 24 hours. Primary human monocytes were isolated from human blood using Ficoll-Paque density gradient centrifugation (GE Healthcare) and CD14+ cells were selected and cultured and differentiated in RPMI-1640 supplemented with 10% FBS and 20 ng/mL GM-CSF for 7 days. For experiments with 1,25D treatment, media was replaced by RPMI-1640 leucine/valine/isoleucine-free (BCAAs-free) supplemented with 10% charcoal-stripped FBS, BCAAs, and 1,25D or vehicle for 24 hours. For BCAA add-back experiments, after washing with phosphate-buffered saline, cells were starved for 1 hour in BCAA-free and FBS-free RPMI-1640, consistent with several protocols (59-61), followed by addition of BCAAs in 10% charcoal-stripped FBS for 10, 30, or 60 minutes.

RNA extraction, reverse transcription and qPCR

RNA extraction was performed with the FavorPrep™ Tissue Total RNA Mini Kit (Favorgen – FATRK 001) as per the manufacturer’s instructions. cDNA was obtained from 1 µg of RNA using 5× All-in-One RT Mastermix (abm - G485) and diluted 5 times. Quantitative polymerase chain reaction (qPCR) was performed with BrightGreen 2× qPCR MasterMix (abm – MasterMix-LR-XL) on a Roche Applied Science LightCycler 96 machine. Expression of targeted genes was normalized to 18S RNA. All primers are listed in Table S5 (51).

Western Blotting and protein analysis

After differentiation into macrophages, THP-1 cells and primary human monocytes were solubilized in Lysis Buffer (20 mM Tris, pH 8, 150 mM NaCl, 1% Triton X-100, 3,5 mM sodium dodecyl sulfate, 13 mM deoxycholic acid) and proteins extracted were separated on a 4% to 15% Tris/Glycine/sodium dodecyl sulfate gel (Bio Rad). A standard protocol was used for transfer and blotting. All primary antibodies and the antirabbit IgG HRP-linked secondary antibody were purchased from Cell Signaling Technology and used at recommended concentrations (see the list below). Signal was detected using Clarity or Clarity Max ECL chemiluminescent substrates (Bio-Rad) and ChemiDoc Imaging System from the same company was used to quantify the intensity of bands. Changes in protein levels were quantified relative to control using Image Lab software after normalization to β-actin. Western blot images presented are representative of three biological replicates.

List of Western blot antibodies: β-actin: Cat# 4967 (62), Phospho-mTOR (S2448): Cat# 2971 (63), mTOR: Cat# 2972 (64), Phospho-p70 S6 Kinase (Thr389): Cat# 9234 (65), p70 S6 Kinase: Cat# 9202 (66), BCAT1: Cat# 12822 (67), Phospho-BCKDH-E1α (Ser293): Cat# 40368 (68), BCKDH-E1α: Cat# 90198 (69), Antirabbit IgG, HRP-linked Antibody: Cat# 7074 (70).

Generation of BCAT1 knockout THP-1 cell line

CRISPR/Cas9 genome editing was carried out using the protocol of Malina et al. (71). The target sequence, AAGATGCGGCTTGCCAGCTT, was designed using the website http://chopchop.cbu.uib.no/ to target human BCAT1 gene and was cloned in the LeGO-U6LTR-iGCas9-GFP plasmid. The plasmid is a gift from Dr. Jerry Pelletier, McGill University. To produce viral particles containing LeGO-U6LTR-iGCas9-GFP, HEK293T/17 cells were transfected along with packaging plasmids using Calfectin (SL100478, Signagen) reagent. To transduce THP-1 cells, the medium from HEK293T/17 cells containing lentiviral particles was collected, filtered, and added to THP-1 cells. To increase the efficiency of the transduction, the plate containing THP-1 cells in medium with viral particles was centrifuged 1 minute at 1000 RPM and transferred to the incubator. The cells were grown for 3 weeks and then a single-cell green fluorescence protein selection was made by fluorescence-activated cell sorting. Different clonal cell lines were then tested via Western blot for BCAT1 and β-actin expression on the same gel to confirm knockouts.

Immunostaining and particle analysis

Imaging was performed at the Advanced Bioimaging Facility at McGill University. Primary human monocytes or THP-1 cells were cultured and differentiated into macrophages on glass coverslips. Cells were fixed with 4% paraformaldehyde for 15 minutes at room temperature and washed 3 times. After incubation in a permeabilizing/blocking buffer (phosphate-buffered saline + 0,1% Triton + 1% bovine serum albumin), lysosomes were labelled with a primary polyclonal rabbit anti-LAMP1 antibody (1:1000; Abcam Cat# ab24170 (72)), which was detected by a secondary Alexa-488 (for primary macrophages) or Alexa-568 (for differentiated THP-1 cells) conjugated goat antirabbit F(ab′)2 antibody (1:1000; Thermo Fisher Scientific Cat# A-11034 (73); Thermo Fisher Scientific Cat# A-11036 (74)), and nuclei were labelled with DAPI (1/5000; Sigma-Aldrich Cat# MBD0015). Images were acquired on a LSM780 confocal microscope (Zeiss AxioObserver) with a Plan Apochromat 63X/NA=1.4 objective. Images were analyzed using ImageJ software, after the images were thresholded, LAMP1-positive particles were detected, and their number and size were determined using “Analyze particles” tool. Bar graphs are mean ± SD of 3 biological replicates (N). For each replicate, 12 (n) images were taken per condition (N = 3; n = 12).

Statistics

Student’s t-test were performed using GraphPad software.

Results

Reanalysis of 94 available human and mouse profiles of gene expression regulated by 1,25D or its analogs

To investigate in depth the effects of 1,25D or its analogs on gene expression and to define common and unique aspects of vitamin D signaling in human and mouse, we analyzed available microarray or RNA-seq expression profiles. Only datasets from Array Express or the Gene Expression Omnibus public repository with available raw data were included, which yielded 80 human and 14 mouse expression profiles from epithelial cells, fibroblasts, myelomonocytic cells, T cells, B cells, granulocytic cells, and muscle (Table S1; all supplemental tables and figures are available in the following data repository (51)). Differential gene expression analysis was performed using the same pipeline (Fig. 1A) in order to ensure that results were as comparable as possible (see “Materials and Methods” for details). We used the LIMMA Bioconductor package as it performs very well in a variety of settings for both microarray and RNA-seq experiments (48). Lists of differentially expressed genes from same-cell–type datasets were combined and filtered in such a way as to exclude genes that appeared as both up- and downregulated. Over-represented genes for each cell type were selected based on the following criteria: (1) Genes had to be regulated by at least 1.5-fold and P ≤ .10 in individual datasets; (2) upregulated genes could not appear as downregulated in any dataset of the same cell type and vice versa; and (3) because at least 3 datasets were available for all cells types except T cells (1 study only), consistently up- or downregulated genes had to appear in at least 3 datasets. Relaxed parameters were used to partially circumvent limitations imposed by low statistical power in individual datasets arising from small numbers of replicates. Quality control of RNA-seq and microarray datasets was assessed; examples of quality control criteria for human and mouse expression profiles are illustrated in Figs. S1 and 2 (51). Further details on what constitutes an experiment of sufficient quality to include in this analysis are provided in the materials and methods section. Numbers of genes regulated at least 1.5-fold in at least 3 independent studies for each cell type are shown in Fig. 1B; a Venn diagram depiction of data from 5 cell types showed that the majority of these were unique to each cell type, underlying the cell specificity of vitamin D signaling (Fig. 1C and 1D). (Note that the numbers of unique genes in Fig. 1B differs from those in Fig. 1C and 1D because data in Fig. 1B includes all 8 cell types in the table.) Given that many datasets were generated from cells treated for extended times with 1,25D, the genes identified represent a combination of direct and indirect targets of vitamin D signaling. Remarkably, only 2 genes were upregulated at least 1.5-fold in the 5 cell types examined in Fig. 1C. These include, unsurprisingly, CYP24A1, whose product encodes the enzyme that attenuates vitamin D signaling, and CLMN, which encodes calmin, a poorly characterized calponin-like transmembrane domain protein previously identified as part of a 1,25D-responsive signature in breast cancer (75).

Figure 1.

Global analysis of 1,25D-regulated gene expression. (A) Schematic representation of workflow for re-analysis of 80 human and 14 mouse microarray or RNAseq profiling studies of 1,25D-regulated gene expression. (B) List of genes up- or downregulated 1.5-fold, compiled from different cell types. To be included, genes had to be similarly regulated (i.e., up or down) in at least three datasets (except for T cells, for which only one dataset was available). (C, D) Venn diagrams illustrating partial overlap in genes regulated 1.5-fold by 1,25D in epithelial cells, monocytic cells, fibroblasts, PBMCs and granulocytic cells.

Substantial overlap between human and mouse vitamin D-regulated biological processes but not regulated genes

We explored the degree of conservation of vitamin D signaling in human and mouse by examining enrichment for and conservation of biological processes, and through a related but more specific gene-centric method. Significantly regulated genes for each human dataset were used to perform GO over-representation analysis for biological processes (GO_BPs) (Fig. S3, (51)). Some of the GO_BP categories (eg, DNA replication, metabolism, MAPK signaling, etc.) were related to cell growth, proliferation, and regulation of intracellular signaling cascades. Notably, a number of epithelial datasets clustered around extracellular matrix organization-related categories. GO over-representation analysis for cellular components supports this observation, as there was an accumulation of epithelial datasets for categories related to adhesion to and formation of extracellular matrix (Fig. S3 (51)). Leukocyte datasets, on the other hand, clustered around granule formation, which generally occurs following activation (Fig. S3 (51)). Indeed, GO_BP enrichment revealed significant aggregation of these datasets for activation and migration (chemotaxis/homing). Interestingly, enrichment for the chemotaxis was also observed for some fibroblastic and epithelial datasets indicating that these cell types may participate in the regulation of processes typically associated with leukocytes. This can be seen more clearly in Fig. S4 (51), where only biological processes related to immune homeostasis (adhesion, wound healing, extracellular matrix, cell–cell contacts, activation, and migration, chemotaxis, or homing) are displayed for all human datasets. Datasets related to adhesion are boxed in blue, whereas those related to lymphocyte or leukocyte activation are outlined in red and granulocyte/neutrophil activation is boxed in yellow (Fig. S4 (51)).

We employed the same approach for all available mouse datasets. There appeared to be an overlap with the immune homeostatic categories (Fig. S5 (51)) for which the human datasets were enriched, especially for myeloid cells. Quantification of enriched biological processes for human and for mouse datasets (Figs. 2A-C) revealed substantial overlap despite the fact that a number of human datasets (lymphocytes and peripheral blood mononuclear cells [PBMCs]) had no mouse counterparts. This is consistent with vitamin D signaling affecting similar patterns of biological processes in human and mouse. When datasets were paired based on tissue type, the overlap slightly diminished for myeloid cells (Fig. 2B) and decreased more substantially for epithelial cells (Fig. 2C). Generally, epithelial gene expression profiles were quite heterogeneous, consistent with their varied tissues of origin. In addition, the human epithelial datasets are greater in number (21 vs 6) and diversity (breast, corneal, liver, lung, prostate, and pancreatic vs colon and prostate) than mouse (Table S1 (51)). Note that only 16 of the 21 human epithelial datasets (Fig. S4 (51)) had expression profiles that were associated with a number of differentially regulated genes large enough to allow for statistically significant enrichment in GO_BP categories.

Figure 2.

(A-C) Overlap of enriched biological processes between human and mouse. Lists of statistically significantly (P.adjust ≤ .05) enriched biological processes across all datasets for human and for mouse were compiled. The overlaps in biological processes for all cell types (A), myeloid (B), or epithelial (C) cells are illustrated using Venn diagrams. The digits in each diagram indicate the number of biological processes. (D,E) Lists of over-represented up- or downregulated genes for each cell type were compiled for human and mouse. Venn diagrams illustrate the intersection of up- or downregulated over-represented genes between human and mouse dendritic cells (D) or human monocytic and mouse dendritic cells (E).

There was a much greater and somewhat greater overlap in upregulated than downregulated genes in epithelial and monocytic cells, respectively, which is expected as transrepression occurs by more heterogeneous mechanisms than transactivation, and may occur via regulation of tissue-specific transcription factors (17). The same selection criteria were applied to mouse datasets and the intersections of over-represented genes in human/mouse Dendritic Cells (DCs) (Fig. 2D and Table S2 (51)), human monocytic/mouse DCs (Fig. 2E and Table S3 (51)), and human/mouse epithelial cells (Fig. S6A and S6B and Table S4 (51)) were minimal. The distributions of the number of conserved down- and upregulated genes in mouse and human datasets for each cell type are shown using violin plots in Fig. S6C and S6D (51), respectively. In conclusion, while there were substantial similarities between human and mouse cells in terms of biological processes affected by vitamin D signaling, the underlying gene-regulatory events appear quite divergent between the 2 species.

Disease ontology enrichment and pathways analysis reveal extensive regulation of immune function by vitamin D signaling

We performed disease ontology enrichment analysis in all human datasets to identify changes in global gene expression signatures associated with various disease states (Fig. 3). Regulated genes were associated with a number of cancers, infection, and autoimmune and inflammatory disease. The etiologies of the last three categories implicate (altered) vitamin D signaling in compromised immune homeostasis. The molecular basis for the potential roles of 1,25D in immune-related human diseases was probed further by performing pathways analyses on enriched genes from each cell type. We focused primarily on epithelial, monocytic and granulocytic cells because they yielded the greatest number of regulated genes. This analysis revealed extensive regulation of genes implicated in Nucleotide-binding Oligomerization Domain-like (NOD-like) receptor signaling and inflammasome function (Figs. S7-9, (51)). Our previous studies had shown that the gene encoding NOD2 is a direct target of the VDR in epithelial and monocytic cells (21). The current data reveals that, in addition to NOD-like receptors, multiple core components of the inflammasome and its negative regulators are induced by 1,25D, particularly in epithelial cells and to lesser extents in granulocytic and monocytic cells (Figs. S7-9 (51)). We confirmed the regulation by 1,25D of core components of the inflammasome PYCARD (which encodes ASC) and CASP1 in the well-differentiated, 1,25D-sensitive oral squamous epithelial cell line SCC25 and in primary keratinocytes (Fig. S7B and S7C, (51)). We also found that expression of CARD16 (COP1) and PYDC1 (POP1, ASC2), negative regulators of inflammasome signaling (76), was enhanced by 1,25D. However, consistent with the data generated by pathways analysis (compare Figs. S7 and S8, (51)), these regulatory events were not observed in primary human macrophages (Fig. S7D, (51)). These data reveal that 1,25D regulates expression of genes encoding several components, both activators and inhibitors, modulating inflammasome signaling in a cell-specific manner.

Figure 3.

Gene enrichment for disease ontology in human datasets. Lists of significantly regulated genes for each human dataset were used to perform a disease ontology enrichment analysis. The top 5 categories for each dataset are displayed. The datasets (x-axis ticks) are grouped based on cell type, as indicated. p.adjust, Benjamini–Hochberg P value adjusted for multiple testing; count, number of genes per category.

Pathways analysis also provided evidence for monocytic cell-specific induction by 1,25D of BCAA metabolism (Fig. 4, and Figs. S10 and S11 (51)). This is of particular interest because BCAAs function as indicators of metabolic status in myeloid cells, and their elevated levels activate signaling through the central metabolic regulatory kinase mammalian Target of Rapamycin (mTOR) to control autophagy, protein synthesis and other metabolic pathways (77). Induction of BCAA catabolism by 1,25D may thus suppress mTOR signaling. It should also upregulate autophagy, a process stimulated by 1,25D signaling (78, 79), as mTOR activation inhibits autophagy (80). Consistent with pathways analysis, the gene encoding the cytoplasmic enzyme branched-chain aminotransferase BCAT1, which converts BCAAs to ketoacids, is induced in differentiated THP-1 macrophage-like cells and in primary human macrophages (Fig. 5A and 5B). Similarly, expression of genes encoding components of the downstream multicomponent mitochondrial branched ketoacid dehydrogenase complex (BCKDH), which catalyzes the irreversible oxidative decarboxylation of branched-chain ketoacids, are also induced in the same cell types (Fig. 5A and 5B). In addition, expression of the gene encoding the amino acid transporter SLC7A5 is stimulated by 1,25D, suggesting that BCAA uptake is enhanced in 1,25D-treated cells (Fig. 5A and 5B). None of these regulatory events was observed in epithelial SCC25 cells (Fig. S11B (51)), consistent with the cell-specificity observed in pathways analysis.

Figure 4.

Pathways analysis revealing an enrichment for genes regulated by 1,25D implicated in catabolism of branched-chain amino acids leucine, valine and isoleucine in monocytic cells. Red is upregulated and blue represents downregulated genes. Key genes are as follows: 26.1.42, branched-chain aminotransferase 1 (BCAT1); branched-chain ketoacid dehydrogenase components 1.2.4.4, branched-chain ketoacid dehydrogenase E1 alpha subunit (BCKDHA); 1.8.1.4, dehydrolipoate dehydrogenase (DLD); 2.3.1.168 dihydrolipoamide branched-chain transacylase (DBT).

Figure 5.

Validation of regulation by 1,25D of genes implicated in branched-chain amino acids catabolism. Cells in all experiments were treated with 1,25D (100 nM) or vehicle for 24 hours. (A) RT/qPCR analysis of 1,25D-regulated expression in differentiated THP-1 monocytic cells of BCAT1, genes encoding components of the BCKDH complex (BCKDHA, DLD, and DBT), the organic acid transferase SLC7A5, which controls BCAA uptake, as well as BCKDK and PPM1K, which encode the kinase and phosphatase that inactivate and activate the BCKDH complex, respectively. (B) RT/qPCR analysis of 1,25D-regulated expression of gene related to BCAA catabolism in primary human macrophages. (C) Western analysis of expression of BCAT1 protein expression in vehicle- or 1,25D-treated primary human macrophages, as indicated. (D) Western analysis of expression of BCKDHA and its phosphorylated form in primary human macrophages and differentiated THP-1 cells treated with vehicle or 1,25D. Results are average of 2 or 3 biological replicates, bars are mean ± SD *P < .05, **P < .01, ***P < .001, ns = not significant.

Western analysis showed that BCAT1 protein levels are also elevated in 1,25D-treated primary human macrophages (Fig. 5C), although those of the rate-limiting alpha-ketoacid dehydrogenase component (BCKDHA or E1 component) of BCKDH are not affected in spite of its elevated gene expression (Fig. 5D). However, we found that levels of BCKDHA phosphorylation are reduced in 1,25D-treated cells (Fig. 5D), consistent with complex activation. BCKDH activity is inhibited by phosphorylation on S293 by the mitochondrial branched-chain ketoacid dehydrogenase kinase (BCKDK) and activated by dephosphorylation by PPM1K phosphatase (81). Expression of BCKDK transcripts is not inhibited and that of the PPM1K gene is not substantially altered in 1,25D-treated cells (Fig. 5A and 5B), suggesting that 1,25D regulates other pathways that control the enzymatic balance between BCKDK phosphorylation and dephosphorylation.

Collectively, these results suggest that 1,25D signaling in macrophages should enhance BCAA catabolism and thus attenuate the activation of mTOR. To test this, we pretreated primary macrophages with vehicle or 1,25D for 24 hours and then starved cells for 1 hour prior to adding BCAAs back to the media, consistent with several published protocols (59-61). In triplicate biological replicates, pretreatment with 1,25D attenuated the effect of BCAA addition on mTOR phosphorylation at S2448 (Fig. 6A), characteristic of active mTORC1 complexes (82), as well as on induction of P70S6K phosphorylation (Fig. 6B), a target of mTOR signaling. Although there was considerable variation in the fold induction of mTOR and P70S6K phosphorylation between each biological replicate, 1,25D attenuated the induction in each replicate. As mTOR inhibits autophagy, a process important in macrophages for eliminating intracellular pathogens (83, 84), we tested the effects of 1,25D in primary human macrophages on the number and size of LAMP1-positive bodies, markers of autophagosomes. 1,25D treatment significantly increased numbers of these bodies and the total lysosomal surface as measured by LAMP1 staining (Fig. 7A and 7B). These experiments were repeated with wild-type differentiated THP-1 cells and cells in which BCAT1 expression had been ablated by CRISPR-Cas genome editing (Fig. 8A and 8B). While the number of LAMP1-positive bodies increased in 1,25D-treated wild-type THP1 cells, consistent with the data obtained with primary macrophages, no such changes were observed in THP-1 cells depleted for BCAT1 (Fig. 8C and 8D). Collectively, these results are consistent with acceleration of BCAA catabolism by vitamin D signaling inhibiting mTOR activity and contributing to activation of autophagy, a key component of innate immune control of intracellular pathogens.

Figure 6.

Attenuation of mTOR activation upon BCAA stimulation. (A) i. Western analysis of mTOR phosphorylation (S2448) upon stimulation of primary human macrophages with BCAAs. Cells were pretreated for 24 hours with vehicle or 1,25D (100nM) prior to starvation and BCAA supplementation. ii. Quantification of mTOR and mTOR phosphorylation (P-mTOR) at different time points. Controls for each time point were set at 1 and each 1,25D treated condition was plotted relative to the corresponding control. (B) i. Induction of p70S6K phosphorylation after BCAA stimulation in primary macrophages pretreated with vehicle or 1,25D (100nM). ii. Quantification of p70S6K (p70) and p70S6K phosphorylation (P-p70) at different time points. Controls for each time point were set at 1 and each 1,25D treated condition was plotted relative to the corresponding control. A 10-minute time point was not done in all experiments.

Figure 7.

1,25D enhances expression of lysosomal marker LAMP1 expression in human primary macrophages. (A) Representative confocal images of LAMP1 staining in green and nuclei DAPI-stained in blue in macrophages after 24 hours of treatment with vehicle or 100 nM 1,25D. (B) Quantification of the number of LAMP1 positive particles per cell and of the total lysosomal surface in control or 1,25D treated macrophages. Results are average of 3 biological replicates, bars are mean ± SD. *P < .05. Scale bar = 10 µm.

Figure 8.

1,25D-induced increase in lysosomal marker LAMP1 is dependent on BCAT1 expression in differentiated THP-1 cells. (A) Ablation of BCAT1 expression in THP-1 cells by CRISPR/Cas genome editing. Clone 24 was used for subsequent studies. (B) 1,25D-induced BCAT1 expression in wild-type THP-1 cells but not in clone 24. (C) Representative confocal images of LAMP1 staining in red and DAPI-stained nuclei in blue in differentiated THP-1 cells after treatment with vehicle or 100 nM 1,25D (24 hours). Top: wild-type. Bottom: BCAT1-knockout cells. (D) Quantification of the number of LAMP1 positive particles per cell in control or 1,25D treated THP-1 cells wild-type (top), BCAT1 KO (bottom). Results are average of 2 or 3 biological replicates, bars are mean ± SD. *P < .05, **P < .01, ns, not significant. Scale bar = 10 µm.

Discussion

The analysis of the global patterns of 1,25D-regulated gene expression performed above strongly supports the broad and pleiotropic physiological actions of vitamin D signaling. Notably, the vast majority of the gene expression profiling studies to date were performed in cells derived from tissues that are not directly implicated in calcium homeostasis. In addition to epithelial cells of diverse origins, monocytic and granulocytic cells were particularly rich sources of vitamin D-regulated gene expression. Analysis of the datasets derived from multiple cell types showed that there were substantially more commonly (at least three datasets) upregulated than downregulated genes, which is to be expected. ChIPseq studies have revealed a strong enrichment for VDREs in genes upregulated by 1,25D (85, 86), consistent with transactivation being driven largely by the hormone-bound VDR interacting with response elements adjacent to regulated genes. On the other hand, transcriptional repression by the VDR is considerably more heterogeneous (17) and can occur via interactions with several other classes of transcription factors, some of which may be expressed in a cell- and tissue-specific manner.

The data also showed there is considerable overlap in enriched biological processes in human and mouse even though several human datasets, such as lymphocytes and PBMCs had no mouse counterparts. It is thus consistent with 1,25D regulating similar biological processes in human and mouse. However, the intersections of over-represented up- or downregulated genes in several human and mouse datasets were minimal. The relative paucity of gene expression profiling studies in mouse relative to human likely contributed substantially to the lack of overlap. Nonetheless, the underlying gene-regulatory events appear quite divergent between the 2 species. This reinforces our previous results, which revealed that a number of genes implicated in 1,25D-regulated innate immune responses were regulated in human cells, but not in mouse (27, 29, 30).

Collectively, the data show that 1,25D-regulated gene expression controls several fundamental cellular processes related to cell growth, proliferation, and intracellular signaling. The diverse epithelial datasets were enriched for categories related to adhesion and formation of extracellular matrix, although these categories were also enriched in profiles from other cell types such as monocytes. Leukocyte datasets clustered around granule formation and regulatory events associated with activation and migration (chemotaxis/homing). Extracellular biological processes are also regulated by 1,25D as response to extracellular stimulus and positive regulation of cytokine production were enriched in several cell types. The human datasets also reveal the widespread effects of vitamin D signaling on immune-related processes in several cell types, notably innate immune responses to pathogen threat and downstream T cell regulation. While it was known that 1,25D induced the expression of the pattern recognition receptor NOD2 (21), the data reveal widespread effects on NOD-like receptor signaling, particularly in epithelial cells. This includes upregulation of both core components of the inflammasome and some of its negative regulators.

One of the pathways enriched in monocytic cells was branched-chain amino acid catabolism. BCAAs are essential amino acids, and as such can act as sensors of nutritional status. A deficiency in BCAAs leads to downregulation of signaling through mTOR, a central regulator of cellular metabolism. Activation of mTOR will enhance protein synthesis and downregulate autophagy (80). While autophagy contributes to recovery of amino acids through enhanced lysosomal protein degradation under conditions of nutritional deprivation, it is also a critical component of macrophage responses to intracellular infections (83, 84). In this regard, it is noteworthy that inhibition of mTOR signaling in vivo protected mice against lethal infection with Listeria monocytogenes, an intracellular pathogen (87).

The gene expression profiles revealed an upregulation of BCAT1 and components of the branched chain ketoacid dehydrogenase complex in monocytic cells. While BCAT1 catalyzes the reversible production of keto-acids, downstream, the large multisubunit BCKDH complex catalyzes a series of irreversible reactions to catabolize branched-chain keto acids to the corresponding acyl-CoA derivatives. BCAT1 protein levels were increased in 1,25D-treated THP-1 cells and primary macrophages. While the increased levels of mRNA encoding the alpha-ketoacid dehydrogenase (E1) component of the BCKDH complex, BCKDHA, did not translate into altered protein levels, we did observe a diminution in both THP-1 cells and primary macrophages of BCKDHA phosphorylated on S293, which corresponds to activation of the enzyme (81). Together, the above results are consistent with activated BCAA catabolism in 1,25D-treated cells. They are supported by our observations that priming of mTOR signaling by BCAA supplementation is attenuated in 1,25D-treated cells, and the greater number of LAMP1-positive bodies in 1,25D-treated primary human macrophages, consistent with enhanced autophagy. Importantly, the 1,25D-dependent increase in number of LAMP1-positive bodies was also observed in wild-type differentiated THP-1 cells, but not in a line depleted for BCAT1 expression by genome editing. Collectively, these results reveal a mechanism through which 1,25D signaling augments BCAA catabolism, which leads to a downregulation of mTOR activation in macrophages and enhancement of autophagy. Given the role of autophagy in controlling intracellular pathogens, it also represents a novel component of the 1,25D-regulated innate immune response.

The global analysis of 1,25D-regulated gene expression presented here has several strengths and weaknesses. Most importantly, it provides a comprehensive overview of the tissue specificity of 1,25D-regulated gene expression and a comparison of available gene expression data in human and mouse, although the strength of the comparison of mouse and human data are limited considerably by the relative lack of studies in the mouse. In addition, the analysis has enabled us to explore conserved biological pathways for each cell type on a broader scale. From a biological perspective, this is significant given that vitamin D signaling is observed in virtually all cells, especially those involved in the innate and adaptive arms of the immune system (47). Screening of data for genes regulated in specific cell types such as epithelial cells reveals the common genes/pathways regulated in these cells from diverse origins. However, it also omits key genes that may be regulated in a cell-type specific manner, as for some cell types different subtypes exist. For instance, the epithelial cell category includes breast, colon and prostate epithelial cells, all of which have distinct physiological roles that could alter their responses to vitamin D. For this reason, combining different cell subtypes into several broad categories identified common regulated gene but nonetheless presents a limitation in our study.

It is also important to note that several of the human and mouse cell datasets are not derived from identical cell types. For example, dendritic cells in humans are monocyte derived (Table S1 (51)), but in mice, they are bone marrow derived. Likewise, macrophages in mice are derived from the bone marrow; however, in the human datasets, monocytes are taken from peripheral blood and macrophages originate from human alveolar tissue or are differentiated THP-1 cells (Table S1 (51)). This is significant because there are many approaches to generate dendritic cell or macrophage-like cells for research, and it is apparent that even within a species, there are differences in cells produced by different methods.

Biological processes such as granulocyte activation, neutrophil activation, and neutrophil activation involved in immune response were enriched in several cell types, consistent with activation by 1,25D of neutrophil recruitment and function. Notably, the combination of 1,25D and M. tuberculosis infection of macrophages super-induced expression of neutrophil chemokines such as IL-8 (CXCL8) (27), consistent with 1,25D signaling enhancing neutrophil recruitment to sites of infection. In this regard, activation of macrophages induces local 1α-hydroxylation of 25D, suggesting that recruited neutrophils would be exposed to locally produced 1,25D (19). In this context, it is unfortunate that all the profiling studies in granulocytic cells were performed with undifferentiated HL60 cells, which are a poor model for mature neutrophils. It would therefore be of interest to determine the effects of 1,25D on gene expression profiles in primary human neutrophils. Nonetheless, the data collected provides a detailed overview of the tissue-specificity of 1,25D-regulated gene expression. It also provides several insights into the non-classical actions of vitamin D signaling, in particular its role in immune system regulation.