Blood samples were collected from a single healthy individual (male, age 56 years, body mass index 25.1, vitamin D status 87.6 nM 25-hydroxyvitamin D₃ in serum), who gave written informed consent to participate in the study. All experiments were performed in accordance with relevant guidelines and regulations related to the VitDbol trial (NCT02063334, ClinicalTrials.gov). The research ethics committee of the Northern Savo Hospital district had approved the study protocol (#9/2014). PBMCs were isolated from freshly collected peripheral blood using Vacutainer CPT Cell Preparation Tubes with sodium citrate (Becton Dickinson) according to manufacturer’s instructions. Deconvolution of RNA-seq data from triplicate solvent-treated samples of each of the three models determined the relative amount of B cells (5.5%), T cells (49.1%), NK cells (19.4%), monocytes/macrophages (23.8%) and other cells (2.2%) within the pool of PBMCs.

PBMCs were washed with phosphate-buffered saline and immediately cultured at a concentration of 0.5 million cells/ml in 5 ml RPMI 1640 medium supplemented with 10% charcoal-depleted fetal calf serum, 2 mM L-glutamine, 0.1 mg/ml streptomycin and 100 U/ml penicillin. Cells were kept at 37 °C in a humidified 95% air/5% CO2 incubator. PBMCs were treated within one hour after taking them into culture with 100 ng/ml LPS (Sigma-Aldrich), 5 µg/ml β-1,3(D)-glucan (BG) (Sigma-Aldrich) or their solvent dimethyl sulfoxide (DMSO) (final concentration 0.1%) and 10 nM 1,25(OH)₂D₃ (Sigma-Aldrich) or its solvent ethanol (EtOH) (final concentration 0.1%) using three different models ( Figure 1A ). In model 1, cells were first exposed for 24 h to LPS, BG or DMSO and then either 1,25(OH)₂D₃ or EtOH were added for another 24 h without a wash-out step. In model 2, cells were first stimulated for 24 h with 1,25(OH)₂D₃ or EtOH and then for additional 24 h with LPS, BG or DMSO. In model 3, cells were incubated for 24 h simultaneously with LPS, BG or DMSO and 1,25(OH)₂D₃ or EtOH. Each in vitro experiment had been performed in three biological repeats within one week with cells from the same donor.

Total RNA was isolated using the High Pure RNA Isolation Kit (Roche) according to manufacturer’s instructions. RNA quality was assessed on an Agilent 2100 Bioanalyzer system (RNA integrity number ≥ 8). rRNA depletion and cDNA library preparation were performed using New England Biolabs kits NEBNext rRNA Depletion Kit, NEBNext Ultra II Directional RNA Library Prep Kit for Illumina and NEBNext Multiplex Oligos for Illumina (Index Primers Sets 1 and 2) according to manufacturer’s protocols. RNA-seq libraries went through quality control with an Agilent 2100 Bioanalyzer and were sequenced on a NextSeq 500 system (Illumina) at 75 bp read length using standard protocols at the Gene Core facility of the EMBL (Heidelberg, Germany).

The single-end, reverse-stranded cDNA sequence reads were aligned (without any trimming) to the reference genome (version GRCh38) and Ensembl annotation (version 93) using STAR (version 2.6.0c) with default parameters. Read quantification was performed within the STAR alignment step (–quantMode GeneCounts). Mapped and unmapped read counts are listed in Table S1 . Ensembl gene identifiers were annotated with gene symbol, description, genomic location and biotype by accessing the Ensembl database (version 101) via the R package BiomaRt (version 2.44.1) (29). Gene identifiers missing external gene name annotation, genomic location or being mitochondrially encoded were removed from the datasets. When a gene name appeared more than once, the entry with the highest average number of counts was kept.

Differential gene expression analysis was computed in R (version 3.6.3) using the tool EdgeR (version 3.28.1) (30) that uses negative binomial distribution to model gene counts. The gene-wise statistical test for differential expression was computed using the generalized linear model quasi-likelihood pipeline (31). In order to mitigate the multiple testing problem, only expressed genes were tested for differential expression. The filtering threshold was adjusted to the expression of the low expressed but highly specific vitamin D responsive gene CYP24A1 (cytochrome P450 family 24 subfamily A member 1). For this purpose, read counts were normalized for differences in sequencing depth to counts per million (CPM). Each gene needed to have an expression of > 0.5 CPM in at least 36 out of 54 samples, in order to be considered. This requirement was fulfilled by 16,861 genes. After filtering, library sizes were recomputed and trimmed mean of M-value normalization applied, in order to eliminate composition bias between libraries. The underlying data structure was explored by visualizing the samples via multidimensional scaling (MDS) ( Figure S1 ). MDS was computed via EdgeR’s function plotMDS() in which distances approximate the typical log2 fold change (FC) between the samples. This distance was calculated as the root mean square deviation (Euclidean distance) of the largest 500 log2FCs between a given pair of samples, i.e., for each pair a different set of top genes was selected. The two principal factors distinguishing the samples’ expression profiles were the type of immune challenge and whether they were treated with 1,25(OH)₂D₃. Thus, the meaningful clustering of samples confirmed the similarity of the triplicates and demonstrates the effects of the treatments. In this line, a design matrix was constructed for the following pairwise comparisons: i) LPS/EtOH (LE) with DMSO/EtOH (DE) reference, ii) BG/EtOH (BE) with DE, iii) DMSO/1,25(OH)₂D₃ (DV) with DE, iv) LPS/1,25(OH)₂D₃ (LV) with LE and v) BG/1,25(OH)₂D₃ (BV) with BE. Trended negative binomial dispersion estimate was calculated using CoxReid profile-adjusted likelihood method and together with empirical Bayes-moderated quasi-likelihood gene-wise dispersion estimates used for generalized linear model fitting. The empirical Bayes shrinkage was robustified against outlier dispersions as recommended (31). Finally, quasi-likelihood F-test was applied to inspect, whether the observed gene counts fit the respective negative binomial model. Only genes with a false discovery rate (FDR) < 0.001 and an absolute FC > 2 were considered. Mean-Difference (MA) plots were generated with vizzy (version 1.0.0), (https://github.com/ATpoint/vizzy) to display the expression profile of each of the 15 comparisons ( Figure S2 ).

Relative cell type composition within the PBMC pool was estimated by deconvolution via the algorithm CIBERSORTx (32) using the default LM22 validated gene-signature matrix and gene expression data of solvent-treated samples of all three models. Estimations are based on 1000 permutations. Venn diagrams were created applying the webtool jvenn (33) (http://jvenn.toulouse.inra.fr) and Manhattan plots were produced in R by using packages ggbio (version 1.36.0) (34) and GenomicRanges (version 1.40.0) (35). Based on transcriptome-wide data pathway analysis was performed via the webtool Enrichr (36, 37) (https://maayanlab.cloud/Enrichr/) utilizing the Kyoto Encyclopedia of Genes and Genomes (KEGG) 2019 Human pathways (38). Adjusted P-values were employed for pathway ranking and the threshold < 0.001 was applied. Integrative database Genecards (https://www.genecards.org) was used for gene product locations and functions.

PBMCs of a single healthy individual were stimulated immediately after isolation with LPS, BG or solvent control (DMSO) in the presence of 1,25(OH)₂D₃ or its solvent (EtOH) ( Figure 1A ). Three different models were applied: in model 1 the cells were first exposed to LPS or BG for 24 h and then for another 24 h to 1,25(OH)₂D₃, in model 2 the sequence was changed, i.e., first 1,25(OH)₂D₃ stimulation for 24 h and then treatment with LPS or BG, and in model 3 immune challenges and 1,25(OH)₂D₃ were applied simultaneously for 24 h. The experiments of each model were performed in three repeats followed by RNA-seq and differential gene expression analysis. When the thresholds FDR < 0.001 and absolute FC > 2 were applied, 1255, 1605 and 1198 differentially expressed genes were detected in models 1, 2 and 3, respectively ( Table S2 and Figure S3A ). For comparison, the influence of cell culture conditions like different treatment times (48 h in models 1 and 2 versus 24 h in model 3) were estimated by differential gene expression analysis of solvent-treated samples of each model ( Figure S3B ). These differences were largely model specific (75.1% of all) and only the five genes ACP5 (acid phosphatase 5, tartrate resistant), ALDH1A1 (aldehyde dehydrogenase 1 family member A1), CCL24 (C-C motif chemokine ligand 24), CD302 (cluster of differentiation 302) and SPARC (secreted protein acidic and cysteine rich) were identified as common genes that are sensitive to cell culture conditions.

In 13 of the 15 single and combined treatments the majority of the responsive genes were downregulated ( Figure 1B ). Within a given model, 23.6 to 33.4% of the responsive genes were downregulated in all treatments, while only 7.4 to 11.1% were exclusively upregulated. Thus, the majority (59.2 to 68.5%) of the responsive genes showed a mixed regulation profile ( Figure S3C ). In total of the three models, 1580 genes responded to LPS, 966 to BG and 1006 to 1,25(OH)₂D₃, from which 503, 388 and 201, respectively, have been previously reported (7, 39) ( Figure S3D ).

In all models, a treatment with LPS alone resulted in the highest count of responsive genes, while lowest numbers were obtained by a combined LPS/1,25(OH)₂D₃ treatment ( Figure 1C ). The number of responsive genes was also reduced by BG/1,25(OH)₂D₃ co-treatment but the effect was less prominent. LPS and BG showed 336, 505 and 375 overlapping genes in models 1, 2 and 3, respectively ( Figure S3E ). For comparison, in the presence of 1,25(OH)₂D₃ there were only 107, 177 and 57 common genes ( Figure S3F ). The count of 1,25(OH)₂D₃-responsive genes was only 288 in model 1, but 645 and 676 in models 2 and 3, respectively. Interestingly, the co-treatment with BG in model 1 increased the number of 1,25(OH)₂D₃-responsive genes, while in models 2 and 3 as well as in combination with LPS the numbers declined, i.e., the count and identity of vitamin D responsive genes was dependent on the co-treatment. The LPS treatment in model 2 is an exception, since in this case the ratio between up- and downregulated genes increased from 0.35 to 1.17 due to pre-treatment with 1,25(OH)₂D₃. The number of genes that are responsive to all three treatments, single and in combination, is rather low: 10 in model 1, 50 in model 2 and 12 in model 3 ( Figure 1C ). In contrast, there are 385, 444 and 298 genes that are in models 1, 2 and 3, respectively, exclusively responsive to the single treatment with LPS. These numbers are significantly higher than the counts for single treatments with BG (140, 49 and 50) or 1,25(OH)₂D₃ (76, 113 and 186).

In summary, the transcriptome of freshly isolated PBMCs shows in a time frame of 1-2 days significant (FDR < 0.001) and prominent (absolute FC > 2) changes in 1580 and 966 genes after immune challenges with LPS and BG, respectively, and in 1006 genes following 1,25(OH)₂D₃ treatment. The counts of the primarily downregulated LPS and BG responsive genes are clearly reduced to a total of 407 and 595, respectively, when the cells are treated 24 h after, 24 h before or in parallel with 1,25(OH)₂D₃. Interestingly, only a pre-treatment of the LPS challenge with 1,25(OH)₂D₃ leads to a majority of upregulated genes, while in the five remaining treatment protocols the proportion of downregulated genes even further increases.

In order to identify key genes responding to either immune challenges by LPS or BG or 1,25(OH)₂D₃ modulation, we focused first on single treatments in all models. From the in total 1580 LPS responsive genes only 24.3% responded in all three models ( Figure 2A ). Similarly, only 27.3% of the 966 BG responsive genes ( Figure 2B ) and 15.5% of 1006 1,25(OH)₂D₃ responsive genes ( Figure 2C ) were common to all models. Thus, most responsive genes have a specificity for one or two models suggesting that the sequence of treatment has a major impact on the responsiveness of the cells.

For understanding the common aspects of the three models, we concentrated on joined responsive genes of the single treatments. Manhattan plots displayed the regular genome-wide distribution of the common responsive genes of LPS ( Figure 2D ), BG ( Figure 2E ) and 1,25(OH)₂D₃ ( Figure 2F ). The number of downregulated responsive genes was at all three treatment conditions higher than the count of upregulated genes. Despite the dominance of downregulation, the most prominent gene expression changes were observed for upregulated genes. Applying an absolute FC > 32 (= 25) threshold highlighted 19 LPS responsive genes (13 up and 6 down), 18 BG responsive genes (16 up and 2 down) and 12 1,25(OH)₂D₃ responsive genes (6 up and 6 down) (named in Figures 2D–F ). The vast majority of these responsive genes are protein coding, but HMGN2P46 is a pseudogene and FAM198B-AS1, AC022509.1 and AC037198.1 are non-coding RNA genes. Interestingly, the top responding genes indicated a number of common responsive genes for LPS and BG treatment [CXCL5 (C-X-C motif chemokine ligand 5), CCL1, CD163, ITGB8 (integrin subunit beta 8), INHBA (inhibin subunit beta A), MMP7 (matrix metallopeptidase 7)] but no overlap with 1,25(OH)₂D₃ stimulation.

We used the transcriptome-wide data for pathway analysis using the webtool Enrichr with the 384, 264 and 156 common responsive genes of LPS, BG and 1,25(OH)₂D₃, respectively, pointed to their top five functions based on KEGG pathways. LPS treatment associated with “Cytokine-cytokine receptor interaction”, “Rheumatoid arthritis”, “NOD-like receptor signaling pathway”, “Salmonella infection” and “Osteoclast differentiation” ( Figure 2G ). The first two functions were also found with BG treatment, in addition to “Toll-like receptor signaling pathway”, “Legionellosis” and “Proteoglycans in cancer” ( Figure 2H ). The latter pathway was also associated with 1,25(OH)₂D₃ treatment alongside “Phagosome”, “Hematopoietic cell lineage”, “ECM-receptor interaction” and “Staphylococcus aureus infection” ( Figure 2I ). When the top five pathways were analyzed for each model separately ( Figure S4 ), LPS treatment resulted for all models in “Rheumatoid arthritis” and “Osteoclast differentiation”, the functions “Cytokine-cytokine receptor interaction” and “NOD-like receptor signaling pathway” were found for models 1 and 3 and “Hematopoietic cell lineage” for models 1 and 2, while “Phagosome”, “Leishmaniasis” and “Influenza A” were model-specific ( Figures S2A–C ). BG treatment highlighted the pathways “Rheumatoid arthritis” and “Cytokine-cytokine receptor interaction” in all models, “Leishmaniasis” in models 2 and 3, while “Proteoglycans in cancer”, “Complement and coagulation cascades”, “ECM-receptor interaction”, “Hematopoietic cell lineage”, “Inflammatory bowel disease”, “Legionellosis” and “Salmonella infection” showed a model-specific fashion ( Figures S2D–F ). In contrast, 1,25(OH)₂D₃ triggered pathways in a more diverse way: “Phagosome”, “Staphylococcus aureus infection”, “Tuberculosis”, “Rheumatoid arthritis” and “Leishmaniasis” associated with two models, while “Hematopoietic cell lineage”, “Toxoplasmosis”, “Cytokine-cytokine receptor interaction”, “Osteoclast differentiation” and “Fluid shear stress and atherosclerosis” were found to be model-specific ( Figures S2G–I ).

Representative responsive genes were selected on the criteria i) being responsive to all treatments in at least one model ii) displaying prominent changes in expression and iii) being involved in the top KEGG pathways. The genes TMEM176A (transmembrane protein 176A), WNT5A (WNT family member 5A), CXCL1, S100A8 (S100 calcium binding protein A8), TNFSF15 (TNF superfamily member 15), CSF1 (colony stimulating factor 1), CD163, INHBA, CCL1, MMP9, CDKN1A (cyclin dependent kinase inhibitor 1A) and TREM1 (triggering receptor expressed on myeloid cells 1) all represent previously reported LPS, BG or 1,25(OH)₂D₃ responsive genes (7, 40, 41) ( Figure S5 ). They represent a 4x3 matrix indicating that the whole group of responsive genes can be classified into 12 categories, such as being primarily responsive only to LPS or BG, both LPS and BG, or only 1,25(OH)₂D₃, as well being all down- or upregulated or showing a mixed response. This highlighted interesting specificities, such as that CCL1 is clearly responsive both immune challenges but it barely responded to treatment with 1,25(OH)₂D₃, whereas TREM1 showed distinct preference for 1,25(OH)₂D₃. The examples of the mixed regulation category indicated that immune challenges led to increased gene expression while 1,25(OH)₂D₃ showed opposite regulation. Furthermore, model-specific differences were observed, where, e.g., TNFSF15 showed distinct responsiveness while CSF1 responded almost the same in all models.

Taken together, the immune challenges LPS and BG display characteristic differences in responsive genes and the respective functions mediated by them, but also reasonable overlap in responding genes and regulated pathways. In contrast, 1,25(OH)₂D₃ primarily regulates a distinct set of genes and in case of joined responsive genes often show opposite direction of gene regulation. Despite these differences, all observed top functions relate to innate and adaptive immunity.

For all models, the effects of either single treatments with LPS or BG and 1,25(OH)₂D₃ were compared with their respective combinations ( Figure 3 ). In model 1, LPS and 1,25(OH)₂D₃ treatments overlapped in 112 genes, only 16 of which responded to the combined treatment of LPS and 1,25(OH)₂D₃ ( Figure 3A ). Individual LPS and 1,25(OH)₂D₃ treatments had in model 2 406 identical genes, 97 of which responded also to the combination of both treatments ( Figure 3B ). In model 3, LPS and 1,25(OH)₂D₃ treatments shared 343 genes, only 23 of which were found with their combination ( Figure 3C ). Similar results were obtained for immune challenge with BG, but compared to LPS the overlaps were larger: in model 1 127 BG and 1,25(OH)₂D₃ responsive genes overlapped, 47 of which in the context of dual stimulation ( Figure 3D ), in model 2 there were 321 identical genes, 123 of which responded to both stimuli ( Figure 3E ), and 320 shared genes in model 3, 89 of which occurred with both treatments ( Figure 3F ).

The combined treatments had reduced the total number of vitamin D responding genes to 407 in presence of LPS ( Figure S6A ) and 595 together with BG ( Figure S6B ). Interestingly, only 23 genes were commonly responding in all models to LPS/1,25(OH)₂D₃, while for BG/1,25(OH)₂D₃ the number was with 166 far higher. Furthermore, the model-specific combined responsive genes were in model 2 with 226 and 191 genes for LPS and BG co-treatment, respectively, clearly higher than in model 1 (66 and 94 genes) and model 3 (15 and 17 genes). Although model 2 had for combined LPS/1,25(OH)₂D₃ treatment a larger responsive gene count than models 1 and 3, only the pathways “ECM-receptor interaction” and “Cytokine-cytokine receptor” passed the threshold ( Figure S6C ). The latter function was also found in model 3, while all five top pathways of model 1 (“Phagosome”, “Proteoglycans in cancer”, “Legionellosis”, “Tuberculosis”, “Amoebiasis”) as well as the remaining four of model 3 (“Allograft rejection”, “Malaria”, “Rheumatoid arthritis” and “Pertussis”) were model-specific. In contrast, for the BG/1,25(OH)₂D₃ combination models 2 and 3 shared the top five pathways “Hematopoietic cell lineage”, “Phagosome”, “Tuberculosis”, “Cytokine-cytokine receptor interaction” and “Osteoclast differentiation” and model 1 at least the first three of them ( Figure S6D ). The two specific pathways of model 1 were “Staphylococcus aureus infection” and “Asthma”. Compared with the pathways highlighted by single treatments, the combined treatments relate more to infectious diseases and their specific pathogens.

Responsive genes serving as representative examples for the effects of combined treatments in comparison with single treatments ( Figure S7 ) were selected by the same criteria as in case of the latter ( Figure S5 ). The combined treatments showed either a boosting, inhibitory or mixed effect on gene expression. Moreover, genes were sorted by being under all conditions downregulated, upregulated or showing a mixed response providing each a 3x3 matrix for LPS and BG. Representative genes for LPS response were FPR3 (formyl peptide receptor 3), TGFBI (transforming growth factor beta induced), ITGB2 (integrin subunit beta 2), CD14, FBP1 (fructose-bisphosphatase 1), SEMA6B (semaphoring 6B), SLC22A23 (solute carrier family 22 member 23), CXCL5 and STAG3 (stromal antigen 3) ( Figure S7A ). The genes TLR4, HLA-DRB5 (major histocompatibility complex, class II, DR beta 5), CCL2, CLMN (calmin), IL1RN (interleukin 1 receptor antagonist), IL1R1 (interleukin 1 receptor type 1), GAL3ST4 (galactose-3-O-sulfotransferase 4), HBEGF (heparin binding EGF like growth factor) and G0S2 (G0/G1 switch 2) represent the BG response ( Figure S7B ). With exception of the genes HLA-DRB5, SLC22A23, STAG3 and GAL3ST4 the example genes are already known as LPS, BG and/or 1,25(OH)₂D₃ responsive genes (7, 39, 42).

In summary, the number of genes responding both to immune challenge and vitamin D, alone and in combination, indicate a descending ranking of models 2, 3 and 1. The joined response to BG and vitamin D shows a far better consensus between the models than that of LPS and vitamin D, both in gene count as well as by pathways. Responsive genes are either boosted or inhibited by dual treatments and often show mixed responses depending on the chosen model.

Integrating the functional consequences of the treatment sequence based on pathway analysis of single ( Figures 2G–I and S2 ) and combined ( Figures S6C, D ) stimulation highlighted the differences of the three models. In model 1, immune challenge with LPS caused chemotaxis and induced cytokine signaling, whereas BG treatment affected proliferation, cell growth and cell migration but also increased cytokine signaling ( Figure 4A ). In contrast, stimulation with 1,25(OH)₂D₃-modulated genes and pathways involved in antigen recognition and phagocytosis. Interestingly, the combined treatment changed the effects of the immune challenges. The modulation of the LPS challenge with 1,25(OH)₂D₃ caused a shift towards phagocytosis, proliferation and cell migration, while the response to BG converted by modulation with 1,25(OH)₂D₃ into differentiation and phagocytosis. In model 2, the effects of all single treatments associated with inflammation, which in case of the immune challenges related to cytokines but with 1,25(OH)₂D₃ linked to pathogen inhibition ( Figure 4B ). Vitamin D modulated both immune challenges so that cytokine signaling was inhibited and in case of BG also phagocytosis was affected. In model 3, single treatment with LPS caused chemoattraction and affected pathogen recognition, while that of BG related to cytokine signaling and inflammation induced by pathogens ( Figure 4C ). In contrast, stimulation with 1,25(OH)₂D₃ alone affected differentiation and caused downregulation of phagocytosis, while in combination with LPS it inhibited cytokine signaling and together BG it initiated differentiation.

Only the genes STAB1 (stabilin 1) and HCAR3 (hydroxycarboxylic acid receptor 3) were in all models responsive to all types of treatments and serve as master examples for monitoring the differences between the models ( Figure 4D ). The STAB1 gene encodes for a highly expressed membrane protein involved in endocytosis, which in every model was downregulated by all types of treatments ( Figure 4E ). The LPS/1,25(OH)₂D₃ co-treatment clearly reduced the change of downregulation being caused by respective single treatments. In contrast, the BG/1,25(OH)₂D₃ treatment resulted in model 1 in an enhanced change in downregulation, in model 2 in no significant effect and in model 3 in a slightly reduced change in downregulation. The HCAR3 gene encodes for a G protein-coupled receptor with low affinity for nicotinic acid. In PBMCs the gene shows a low basal expression, was upregulated by both immune challenges but downregulated by 1,25(OH)₂D₃ and combined treatment. However, the combined treatments led to less change in downregulation than 1,25(OH)₂D₃ alone. Changes in HCAR3 gene expression did not vary much between the three models, although in model 2 LPS had the lowest and BG the highest effect.

Taken together, a co-stimulation with 1,25(OH)₂D₃ is able to change the functional consequences of immune challenges but there are large differences as consequence of treatment sequence, i.e., of the chosen model. The genes STAB1 and HCAR3 are master examples monitoring the complex model-specific response to the modulation of immune challenges by vitamin D.