Differential and Overlapping Effects of 20,23(OH)2D3 and 1,25(OH)2D3 on Gene Expression in Human Epidermal Keratinocytes: Identification of AhR as an Alternative Receptor for 20,23(OH)2D3

A novel pathway of vitamin D activation by CYP11A has previously been elucidated. To define the mechanism of action of its major dihydroxy-products, we tested the divergence and overlap between the gene expression profiles of human epidermal keratinocytes treated with either CYP11A1-derived 20,23(OH)2D3 or classical 1,25(OH)2D3. Both secosteroids have significant chemical similarity with the only differences being the positions of the hydroxyl groups. mRNA was isolated and examined by microarray analysis using Illumina’s HumanWG-6 chip/arrays and subsequent bioinformatics analyses. Marked differences in the up- and downregulated genes were observed between 1,25(OH)2D3- and 20,23(OH)2D3-treated cells. Hierarchical clustering identified both distinct, opposite and common (overlapping) gene expression patterns. CYP24A1 was a common gene strongly activated by both compounds, a finding confirmed by qPCR. Ingenuity pathway analysis identified VDR/RXR signaling as the top canonical pathway induced by 1,25(OH)2D3. In contrast, the top canonical pathway induced by 20,23(OH)2D3 was AhR, with VDR/RXR being the second nuclear receptor signaling pathway identified. QPCR analyses validated the former finding by revealing that 20,23(OH)2D3 stimulated CYP1A1 and CYP1B1 gene expression, effects located downstream of AhR. Similar stimulation was observed with 20(OH)D3, the precursor to 20,23(OH)2D3, as well as with its downstream metabolite, 17,20,23(OH)3D3. Using a Human AhR Reporter Assay System we showed marked activation of AhR activity by 20,23(OH)2D3, with weaker stimulation by 20(OH)D3. Finally, molecular modeling using an AhR LBD model predicted vitamin D3 hydroxyderivatives to be good ligands for this receptor. Thus, our microarray, qPCR, functional studies and molecular modeling indicate that AhR is the major receptor target for 20,23(OH)2D3, opening an exciting area of investigation on the interaction of different vitamin D3-hydroxyderivatives with AhR and the subsequent downstream activation of signal transduction pathways in a cell-type-dependent manner.

To better define the signaling pathways and mechanisms underlying the similarities and differences between phenotypic activities of classical 1,25(OH) 2 D3 and the major dihydroxy product of CYP11A1 action on vitamin D3, 20,23(OH) 2 D3, we examined and compared the gene expression profiles of human keratinocytes exposed to these secosteroids. Bioinformatics analysis was performed and differences and similarities in the activities of these structurally similar but distinct dihydroxy-D3 species were compared.
Ingenuity pathway analysis using FC ≥ ±1.5 was performed. The top canonical pathway induced by 1,25(OH) 2 D3 was VDR/RXR signaling (Table 2) (Supplemental Figure S1A), which was expected [3,6,51,81,82]. This was followed by the roles of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis; the role of macrophages, fibroblasts and endothelial cells in rheumatoid arthritis; and Toll-like receptor signaling (Table 2), which is consistent with previously reported functions of 1,25(OH) 2 D3 [3,6,30,31,45,51,54]. Interestingly, the next top nuclear receptor signaling pathway activated by 1,25(OH) 2 D3 was linked to the glucocorticoid receptor (GR) followed by the aryl hydrocarbon receptor [74], PPAR, PPARα/RXRα, LXR/RXR, and RAR ( Table 2). The inclusion of these additional pathways could be secondary to the use of the same dimeric partner, RXR, and communication between receptors, or alternatively by activation by signaling pathways downstream of VDR. For example, it is already known that 1,25(OH) 2 D3 can selectively activate local elements of hypothalamo-pituitary adrenal axis in keratinocytes [71]. The significance of additional nuclear receptor signaling is out of the scope of this paper and is a goal of our future research.  The top canonical nuclear receptor pathway induced by 20,23(OH) 2 D3 was AhR signaling (Supplemental Figure S2A) with VDR/RXR being next (Supplemental Figure S3A) (Table 3). While the identification of the VDR/RXR as the target for 20,23(OH) 2 D3 is consistent with previously reported functional data and molecular modeling [60,65,69,83], identification of the AhR as its primary target was unexpected and hence it was further analyzed in detail as described below. Table 4 shows that for 1,25(OH) 2 D3 the nuclear signaling pathways VDR/RXR, followed by AhR, PPARα/RXRα, RAR and LXR/RXR, were among the top toxicity-related pathways identified. The top signaling pathways for 20,23(OH) 2 D3 were linked to the activation of AhR and VDR/RXR (Table 5). Tables 6 and 7 show certain functional similarities between top diseases and bifunctions affected by both molecules. For example, cancer, and organismal injury and abnormalities, are the top two diseases affected by both molecules. These phenotypic similarities are consistent with previously reported studies comparing the biological effects of 1,25(OH) 2 D3 and CYP11A1-derived D3-hydroxyderivatives, including 20,23(OH) 2 D3, and indicate similarities between the effects on cell proliferation and differentiation, as well as similar anti-inflammatory, photoprotective and anti-cancer actions [23,[60][61][62]64,72,80,84].         Because of the unexpected differences between 1,25(OH) 2 D3 and 20,23(OH) 2 D3, the 6 h incubation experiment was repeated in a similar manner as shown in Figure 2 and microarray analyses were performed using Illumina's HumanWG-6_V2 (Platform GPL13376) chip/array. Average signal values for filtered gene clusters with FC ≥ ±1.5 are shown in Supplemental excel file #2. The heat maps corresponding to relative gene expression and Venn diagrams are shown in Figure 3C. Again, for a 2-fold cut-off value there was only one common gene (CYP24A1) whose expression was stimulated by both 1,25(OH) 2 D3 (80-fold) and 20,23(OH) 2 D3 (2.9-fold). For FC ≥ ±1.5 there were 11 common genes upregulated and 4 downregulated. Again, ingenuity pathway analysis showed that VDR/RXR was the top canonical pathway induced by 1,25(OH) 2 D3, followed by the role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis. As before, other nuclear receptor signaling pathways included LXR/RXR, GR and AhR. For 20,23(OH) 2 D3 the top nuclear receptor signaling pathways were again AhR and VDR/RXR. Of note, this microarray showed that 20,23(OH) 2 D3 upregulated two genes downstream of AhR signaling, CYP1A1 and CYP1B1, by factors of 2.4 and 2.6, respectively. This stimulation was confirmed by qPCR ( Figure 3D). VDR/RXR was identified as the top toxicity pathway for 1,25(OH) 2 D3 and AhR for 20,23(OH) 2 D3.
More robust data were obtained with 24 h of treatment for which the average signal values for filtered gene clusters with FC ≥ ±1.5 are shown in Supplemental excel file #3. Because of the large number of genes affected (Table 1), the heat map of differentially expressed genes and Venn diagrams were generated using the 4-FC value which show three overlapping genes (CYP24A1, MMP3 and SERPINB1) as well as distinct gene expression patterns ( Figure 4). For FC ≥ ±1.5, 93 and 72 common genes were up-and downregulated, respectively. Ingenuity pathway analysis using FC ≥ ±2.0 was consistent with results obtained after 6 h of treatment. Again, the top canonical pathway for 1,25(OH) 2 D3 was VDR/RXR (Supplemental Figure S1B) followed by MIF-related glucocorticoid regulation and regulation of the innate immunity system (Table 8). AhR signaling was also listed. The top canonical pathways induced by 20,23(OH) 2 D3 were AhR signaling (Supplemental Figure S2B) and the cholesterol biosynthesis pathway (Table 9). Interestingly, the involvement of a second nuclear receptor complex was emphasized by VDR/RXR activation (Supplemental Figure S3B), with p53 signaling also being listed. The latter is consistent with the photoprotective properties of 20,23(OH) 2 D3 and activation of p53 by its direct precursor, 20(OH)D3 [23]. The top affected toxicity pathways for 1,25(OH) 2 D3 included VDR/RXR, xenobiotic metabolism, cardiac fibrosis and cytochrome P450s (Table 10). AhR signaling was also listed. For 20,23(OH) 2 D3, AhR signaling was again listed as the top toxicity pathway followed by cholesterol synthesis, p53 signaling and again VDR/RXR (Table 11). The top upstream gene regulation pathways for 1,25(OH) 2 D3 included vitamin D3-VDR-RXR, calcitriol, dexamethasone, progesterone and β-estradiol, while for 20,23(OH) 2 D3, included TP53 (p53 tumor suppressor), β-estradiol, lipopolysaccharide, TNF (tumor necrosis factor) and TGF β1 (transforming growth factor-β1) (not shown).          Tables 12 and 13 show some similarities and differences with dermatological diseases and conditions; with cancer, organismal injury and abnormalities being the main diseases affected by 20,23(OH) 2 D3 and 1,25(OH) 2 D3. With regard to molecular and cellular functions, cellular growth and proliferation, cell death and survival, cellular movement and cell cycle were the major functions for 20,23(OH) 2 D3, and cellular movement, cell signaling, small molecule biochemistry, lipid metabolism and cellular development for 1,25(OH) 2 D3. Among the 25 networks activated by 20,23(OH) 2 D3, the top five included: (1) connective tissue disorders, neurological diseases, organismal injuries and abnormalities, (2) RNA post-transcriptional modification, carbohydrate metabolism and lipid metabolism, (3) connective tissue, developmental, skeletal and muscular disorders, (4) cellular movement, endocrine system disorders, gastrointestinal diseases and (5) nucleic acid metabolism, small molecules biochemistry and dermatological diseases and conditions. Among the 15 networks activated by 1,25(OH) 2 D3, the top five included: (1) cancer, organismal functions, organismal injuries and abnormalities, (2) cell-to-cell signaling and interaction, cellular assembly and organization, cellular development, (3) cellular growth and proliferation, tissue development and cancer, (4) molecular transport, carbohydrate and lipid metabolism and (5) protein degradation, protein synthesis, cellular assembly and organization.   Because of the unexpected finding that AhR signaling represented the top regulatory pathway activated by 20,23(OH) 2 D3, and is validated by qPCR analysis of CYP1A1 and CYP1B1 genes expression ( Figure 3D), we examined whether 20(OH)D3, which is the precursor to 20,23(OH) 2 D3, and 17,20,23(OH) 2 D3 and 1,20(OH) 2 D3, which are downstream metabolites (see Figure 1), also affected the expression of genes linked to AhR in HaCaT keratinocytes. Figure 5A shows that 20(OH)D3 stimulated the expression of CYP1A1 and CYP1B1 in a dose-dependent fashion, with a stimulatory effect also seen for the AhR gene. 17,20,23(OH) 3 D3 (1 µM) could also stimulate CYPA1, CYP1B1 and AhR expression, while 1,20(OH) 2 D3 had only a small effect on CYPB1 and no effect on CYP1A1 and AhR. Finally, we used a Human AhR Reporter Assay System (INDIGO, Biosciences) to analyze the effect of several D3-hydroxyderivatives on AhR-mediated transactivation. The kit contains AhR Reporter Cells that contain the luciferase reporter gene functionally linked to an AhR-responsive promoter, which provides a sensitive surrogate measure of the changes in AhR-mediated activation of luciferase reporter. Figure 6 shows that there was marked activation of AhR activity by 20,23(OH) 2 D3 with weaker but significant activation by 20(OH)D3 or 1,25(OH) 2 D3. Thus, the functional studies support the microarray analysis indicating that hydroxyderivatives of D3 can act on AhR. This finding can be explained by the promiscuous nature of AhR and its activity [85]. An additional mechanistic insight into the above interactions was provided by modeling using the crystal structure of the ligand-binding domain (LBD) of human AhR. The presently available crystal structure of the human AhR (PDB: 5NJ8) is missing the LBD region. A model of the human AhR LBD with bound 20S,23S(OH) 2 D3 was developed as described under Methods. Briefly, the final model was based on the homology modelling template of C-terminal Per-ARNT-Sim domain of Hypoxia-Inducible Factor-2α, PDB entry code 3H82. The sequence identity between human AhR and the modelled sequence is 27%; the alignment is shown in Supplemental Figure S4. Short molecular dynamic simulation runs were performed on selected docked poses of 20S,23R/S(OH) 2 D3 epimers in order to identify binding modes most favorably accommodated in the binding site of the homology model. The selected complex with 20S,23S(OH) 2 D3 was simulated for 100 ns to allow for local structural adjustments of flexible regions to the presence of the vitamin D3 scaffold. The final conformer obtained is referred to as the 'refined AhR model'. Further ligand-induced effects were explored through a 250 ns simulation production run starting with this model. Over the first 130 ns the ligand-induced conformational changes were in the vicinity of F295 and S320. The latter is in a flexible region with two adjacent glycine residues while F295 is part of a loop structure 'covering' the binding pocket. The conformation adopted by 130 ns in these regions were maintained for the rest of the simulation time, likely stabilized by a hydrogen bonding network that formed, involving ligand hydroxyl groups, T289, S320 side chains and the backbone of F295 as shown in a representative simulation snapshot at 230 ns in Figure 7. Interactions of this network link the more rigid beta-sheet structure of the pocket containing T289 with two loop regions. The flexible 'belt' between G309-H326 includes a short helical segment near S320 that also shifted due to the presence of the ligand. This binding mode also changes the preferred orientation of H291 which by 130 ns simulation time forms a stable hydrogen bond with the backbone carbonyl of K292, an interaction not present in the initial or refined AhR models. Alanine mutation of T281, H285 of mouse AhR corresponding to the human residues T289 and H291 was shown to dramatically decrease Hsp90 binding [86]. Figure 7 illustrates that differences in the structural fold of AhR between the homology model, the refined AhR model and the simulation conformer are mainly within loops and the flexible 'belt' region. Data represent means ± SD where * p < 0.05, ** p < 0.01 and **** p < 0.0001 at student t-test; ## p < 0.01 and #### p < 0.0001 at one-way ANOVA test. The proposed binding model of 20S,23S(OH) 2 D3 is shown in Figure 8A through a representative simulation conformer at 230 ns. Figure 8B shows the fraction of simulation time during which interactions are present with each AhR residue, as averaged over 130-250 ns. The most stable polar interactions are hydrogen bonding between 23-OH and T289 at 90% and between 3-OH and S336 at 85%, with S346 also contributing 26% of the simulation time. These interactions anchor the two end regions of the scaffold in the pocket. 20-OH is hydrogen bonding with S320 for 39% of the stimulation time. Due to intra-molecular hydrogen bonding between the ligand hydroxyls, the 20-OH group is positioned to act as a hydrogen bond donor to S320, which allows S320 to interact with F295. This interaction is likely important for these loop conformational changes and their effect on H291. Ligand-protein contacts versus simulation time are shown in Supplemental Figure S5.
Loop conformational changes induced by 20S,23S(OH) 2 D3 are likely specific to this ligand. Therefore, for docking other vitamin D3 analogs the refined AhR model was utilized, applying the Induced Fit method. As shown in Table 14, Glide XP docking scores of three analogs are notably lower compared to other compounds within this set, 20(OH)D3, 1,25(OH) 2 D3 and 1,20(OH) 2 D3. Docked poses of all analogues are very similar, and also closely overlap with 20S,23S(OH) 2 D3 in the refined AhR model. Docked poses for all analogs are displayed in Figure 9, along with the binding mode of 20S,23S(OH) 2 D3 for comparison. Residues from Induced Fit structures contributing to polar interactions with ligands are shown only; all residues in proximity of docked ligands are included in Supplemental Figure S6. Docked Vitamin D3 analogs share similar hydrogen bonding interactions through hydroxyl groups: 1-OH interacts with S365, 3-OH with S336 and possibly S346, 17-OH and 20-OH with S320, 23-OH with T289. Docking results predict that 25-OH interacts with T289. A short, 20 ns molecular dynamic simulation was performed on 20(OH)D3, 1,25(OH) 2 D3, 17,20,23S(OH) 3 D3, starting with docked poses. The ligands maintained the binding mode and predicted interactions during simulation except for 1,25(OH) 2 D3. Therefore, simulation of the latter was extended another 50 ns, during which the pose of 1,25(OH) 2 D3 shifted, disrupting hydrogen bonding between 3-OH and S336 that was only present for 35% of simulation time. In comparison, in the case of 20(OH)D3 and 17,20,23S(OH) 3 D3 the same interaction was present 85% and 66% of time, respectively. Due to mobility of the aliphatic chain in the binding site, 25-OH formed contacts with T289 30% and Y310 37% of the simulation time. 1,25(OH) 2 D3 may have a distinct binding mode and interaction with AhR than the other analogs. Interactions of 17,20,23S(OH) 3 D3 are analogous to those of 20S,23S(OH) 2 D3. However, hydrogen bonding between 17-OH with S320 may interfere with structural changes such as those induced by 20S,23S(OH) 2 D3 during the 250 ns simulation production run. While 20(OH)D3 is also predicted to form analogous contacts, the absence of 23(OH) interactions is likely significant.   −13.0

Modelling Conclusions
Molecular dynamic simulation of the developed AhR-20S,23S(OH) 2 D3 model predicts strong hydrogen bonding interactions between this ligand and T289, S336. A hydrogen bond formed with S320 is also well maintained during simulation. A number of AhR residues have favorable non-polar contacts with the ligand (Figure 8). The simulation trajectory predicts that ligand-specific interactions induce a conformational change in the region in the vicinity of S320 and F295, also leading to a distinct position and interaction of H291. The interaction network that forms during simulation due to the ligand links the beta-sheet structure of the pocket with two loops, restraining the conformation of flexible regions in the binding site. The presented model is also consistent with the observed effect of 20S,23S(OH) 2 D3 on AhR since, in particular, T289 and H291 are essential residues for Hsp90 binding. Docking

Cell Culture
Neonatal foreskins of African American [79] donors were used to isolate neonatal human epidermal keratinocytes (HEK) following standard protocols described previously [69,89]. The use of human tissues were approved both by the IRB at the UTHSC as an exempt protocol #4 and by the IRB at the University of Alabama Birmingham, as they are not subject to FDA regulation and not Human Subject Research. Cells were grown in keratinocyte basal medium (KBM) supplemented with keratinocyte growth factors (KGF) (Lonza, Walkersville, MD, USA) on collagen coated plates [68] and second and third passages were used for the experiments [69]. Human epidermal HaCaT keratinocytes were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with glucose, L-glutamine, pyridoxine hydrochloride (Cell Grow), 5% fetal bovine serum (FBS) (Atlanta Biologicals, Flowery Branch, GA, USA) and 1% penicillin/streptomycin/amphotericin antibiotic solution (Thermo Fisher Scientific, Waltham, MA, USA). Human cells were cultured at 37 • C, with a CO 2 concentration of 5%, 100% humidity, and media were changed every second and/or third day.
Prior to treatment with secosteroids, HaCaT cells were serum deprived for 24 h and the medium was changed to DMEM medium containing 5% charcoal-treated FBS (ctFBS) (Atlanta Biologicals, Flowery Branch, GA, USA) to which D3 hydroxymetabolites from the stock solutions were added. For epidermal neonatal keratinocytes, the KBM with KGF was supplemented with 0.5% bovine serum albumin (BSA) prior to the addition of D3 derivatives.

Microarray Assays
Petri dishes (100 mm in diameter) were seeded with human neonatal keratinocytes that were combined from five different black donors at either passage 2 or 3. After reaching 70-80% of confluence, cells were treated with 10 −7 M of either 20,23(OH) 2 D3 or 1,25(OH) 2 D3, or with 0.1% ethanol (EtOH) as a solvent control for 6 or 24 h. After, these cells were isolated from three plates per each experimental condition and combined for passage 2 and 3, separately ( Figure 2).
The RNA from HEK treated with either 20,23(OH) 2 D3 or 1,25(OH) 2 D3, or 0.1% ethanol control, was isolated using the Absolutely RNA Miniprep Kit (Qiagen, Germantown, MD, USA). High purity RNA samples were subjected to microarray analysis at the Molecular Resources Center at the UTHSC. Expression profiling was accomplished using whole-genome gene expression direct hybridization assay using Illumina's HumanWG-6_V2 (Platform GPL13376) chip/array (Illumina, San Diego, CA, USA). Each array contains full-length 50-mer probes representing more than 22,000 well-annotated RefSeq transcripts, including up-to-date genes derived from the National Center for Biotechnology Information Reference Sequence (NCBI RefSeq) database. Initially, 250 ng total RNA was converted to cDNA, followed by an in vitro transcription step to generate labeled cRNA following the manufacturer's recommendations (Applied Biosystems, Foster City, CA, USA).
The labeled probes were then mixed with hybridization reagents and hybridized overnight to the Human BeadChips. Following washing and staining, the BeadChips were imaged using the Illumina BeadArray Reader to measure fluorescence intensity at each probe. The intensity of the signal corresponds to the quantity of the respective mRNA in the original sample.

Bioinformatics Analysis
For generating networks, a data set containing gene identifiers and corresponding expression values was uploaded into the application. Each identifier was mapped to its corresponding object in Ingenuity's Knowledge Base. A FC of ±2 or ±1.5, where indicated, was set to identify molecules whose expression was significantly differentially regulated. These molecules, called Network Eligible molecules, were overlaid onto a global molecular network developed from information contained in Ingenuity's Knowledge Base. Networks of Network Eligible Molecules were then algorithmically generated based on their connectivity. The Functional Analysis identified the biological functions and/or diseases that were most significant to the entire data set. Molecules from the dataset that met the FC cutoff of ±2 or ±1.5 and were associated with biological functions and/or diseases in Ingenuity's Knowledge Base were considered for the analysis. Right-tailed Fisher's exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that data set is due to chance alone.

Real-Time Reverse Transcription Polymerase Chain Reaction (qRT-PCR)
Semiconfluent cultures of human neonatal keratinocytes or HaCaT cells were treated for 6 h, 24 h or as indicated in the figure legends with vitamin D3 hydroxyderivatives or ethanol, and RNA isolated as described above. Reverse transcription was done using the Transcriptor First Strand cDNA Synthesis Kit (Applied Biosystems, Foster City, CA, USA) with 100 ng RNA per reaction. qRT-PCR was performed using cDNA diluted 10-fold in sterile water and a TaqMan PCR Master Mix. Reactions (in triplicate) were performed at 50 • C for 2 min, 95 • C for 10 min and then 50 cycles of 95 • C for 15 s and 60 • C for 1 min. The primers and probes were designed with the universal probe library (Roche). Data were collected on a Roche Light Cycler 480. The amount of amplified product for each gene was compared to that of Cyclophilin B or GAPDH using a comparative C T (∆∆C T ) method. Supplemental Table S1 lists the primers used for qRT-PCR amplifications.

Interaction of Hydroxyvitamin D Derivatives with AhR
Interaction of 20(OH)D3, 20,23(OH) 2 D3 and 1,25(OH) 2 D3 with AhR was evaluated using the Human AhR Reporter Assay System (INDIGO Biosciences, State College, PA, USA) according to the manufacture's protocol. Briefly, AhR reporter cells were recovered on a 96-well plate frame using the cell recovery medium for 5 h, followed by treatment with vitamin D3 hydroxyderivatives in the compound screening medium for 22 h. After removing the media from the wells, luciferase detection reagent was added to the wells and luminescence was measured using a Cytation 5 Cell Imaging Multi-Mode Reader (Winooski, VT, USA).

Statistical Analyses
Data are presented as means ± SD (n = 3-4), and were analyzed with a Student's t-test (for two groups) or ANOVA using Prism 4.00 (GraphPad Software, San Diego, CA, USA). Statistically significant differences are denoted with asterisks for t-tests or for one way ANOVA with # as indicated in the figure legends.

Data Deposition
The data reported in this paper have been deposited in the Gene Expression Omnibus (GEO) database, https://www.ncbi.nlm.nih.gov/geo (accession no. GEO: GSE117351).

Development of a Human AhR LBD Model Complexed with 20S,23S(OH) 2 D3
The strategy applied utilized tools implemented in the Schrödinger software package, version 2017-4 (Schrödinger, LLC, New York, NY, USA). Homology modelling of the human AhR ligand-binding domain was based on crystal structures of the C-terminal Per-ARNT-Sim domain of Hypoxia-Inducible Factor-2α (HIF-2α PAS-B), PDB entry codes 3H82 and 4XT2. Based on the two templates, two homology models were built using the energy-based homology model building method in Schrödinger. The sequence identity is 27% between human AhR and template sequences in the modelled LBD region. The sequence alignment is shown in Supplemental Figure S4. Residues were numbered according to the human AhR sequence (Uniprot ID P35869). Co-crystallized ligands were included. Modeled loops that contained gaps in the sequence alignment were refined through default loop refinement options. The models were relaxed through restrained energy minimization in Protein Preparation Wizard (OPLS3 force field).
Initial binding mode hypotheses were generated through docking 20S,23R/S(OH) 2 D3 into the two obtained AhR models. Out of the top scoring poses at both AhR models, eight were selected for protein-ligand complex refinement (Prime tool in Schrödinger software), followed by a 10 ns molecular dynamic simulation run for each complex using Desmond. Four poses induced distortions in rigid, AhR beta-sheet/helical backbone structures within 10 ns simulation and were not considered further. Out of the remaining poses, the most favorable contacts were formed by two similar poses of the epimers: 20S,23S/R(OH) 2 D3. Simulation of these poses was extended to 20 ns, which suggested that 20S,23S(OH) 2 D3 is more favorably accommodated than its R epimer. The 20S,23S(OH) 2 D3 complex conformer at 20 ns showed an overall RMSD (root-mean-square deviation) of 1.48 from its homology modelling template (PDB: 3H82). In order to allow flexible regions to adjust to the presence of the bound vitamin D3 scaffold, the model was further simulated for 100 ns or for 230 or 25 ns as indicated. The final conformer was relaxed through restrained energy minimization and is referred to as the 'refined AhR model'. The overall RMSD of this model from its homology modelling template (PDB: 3H82) is 1.58, suggesting stability of the AhR structure over simulation time. The model structure contains only two residues with backbone dihedrals in disallowed regions, both of which are glycines: Gly309 and Gly374.

Docking Method
The induced Fit docking method was used as implemented in Schrödinger, version 2017-4 (Schrödinger, LLC, New York, NY, USA). This method combines Prime tools and Glide docking, taking into account the flexibility of residues in proximity to the ligand. 20,23(OH) 2 D3 was Induced Fit docked into the two human AhR homology models using default parameters except optimization of side chains was extended to 6 Å around the ligand and Glide re-docking was done in extra-precision mode. Docking of vitamin D3 analogs into the refined AhR model also utilized Induced Fit with default options, except that van der Waals scaling parameters for docking were set to 1.0 for both protein and ligand (no scaling), and re-docking was done in extra-precision mode.

Molecular Dynamic Simulation Method
Molecular dynamic simulations were performed using Desmond (Schrödinger, LLC, New York, NY, USA) with the OPLS3 force field. Structures were solved in TIP3P explicit waters with boundary conditions in a 10 Å buffered orthorhombic system. Counter-ions were added. The NPT ensemble was employed with temperature fixed at 300 K and pressure at 1.01 bar. The cutoff radius for Coulombic interactions was set to 10 Å. The trajectory was recorded at 10 ps intervals.

Conclusions
Gene expression profile analysis demonstrated that 20,23(OH) 2 D3 and 1,25(OH) 2 D3 induce distinct and overlapping gene expression patterns in keratinocytes linked to the activation of common (VDR-dependent) and distinct (involving other nuclear receptors) signal transduction pathways. Taking into consideration the strong chemical similarity between 20,23(OH) 2 D3 and 1,25(OH) 2 D3 (Figure 1), the marked differences in gene expression panels (Table 1, Figures 3 and 4) were unexpected. This is because our previous studies predominantly showed phenotypic similarities between the effects of both dihydroxy-vitamin D3 species, such as regulation of cell proliferation and differentiation, and anti-inflammatory, photoprotective and anticancer functions, with only a couple of notable differences where 20,23(OH) 2 D3 displayed no calcemic effects and poor activation of CYP24A1. These predominantly overlapping effects are most likely secondary to the redundancy of downstream phenotypic regulators and intercommunication between distinct transduction pathways at the cell or organ levels. The similarities are likely related to the activation of the VDR. The most significant and unexpected finding is the identification of AhR as the major receptor for 20,23(OH) 2 D3, which also appears to be activated by other CYP11A1-derived vitamin D3 derivatives, and possibly by 1,25(OH) 2 D3 to some extent, as predicted by molecular modeling. The future challenge is to precisely define the interaction of different vitamin D3 hydroxyderivatives with the ligand binding domain of AhR and how it is affected by the location of the OH-group on the side chain or at C1α, and how the activation of downstream signal transduction pathways occurs.