Vitamin D and lumisterol derivatives can act on liver X receptors (LXRs)

The interactions of derivatives of lumisterol (L3) and vitamin D3 (D3) with liver X receptors (LXRs) were investigated. Molecular docking using crystal structures of the ligand binding domains (LBDs) of LXRα and β revealed high docking scores for L3 and D3 hydroxymetabolites, similar to those of the natural ligands, predicting good binding to the receptor. RNA sequencing of murine dermal fibroblasts stimulated with D3-hydroxyderivatives revealed LXR as the second nuclear receptor pathway for several D3-hydroxyderivatives, including 1,25(OH)2D3. This was validated by their induction of genes downstream of LXR. L3 and D3-derivatives activated an LXR-response element (LXRE)-driven reporter in CHO cells and human keratinocytes, and by enhanced expression of LXR target genes. L3 and D3 derivatives showed high affinity binding to the LBD of the LXRα and β in LanthaScreen TR-FRET LXRα and β coactivator assays. The majority of metabolites functioned as LXRα/β agonists; however, 1,20,25(OH)3D3, 1,25(OH)2D3, 1,20(OH)2D3 and 25(OH)D3 acted as inverse agonists of LXRα, but as agonists of LXRβ. Molecular dynamics simulations for the selected compounds, including 1,25(OH)2D3, 1,20(OH)2D3, 25(OH)D3, 20(OH)D3, 20(OH)L3 and 20,22(OH)2L3, showed different but overlapping interactions with LXRs. Identification of D3 and L3 derivatives as ligands for LXRs suggests a new mechanism of action for these compounds.

In addition, Gene Ontology Enrichment Analysis (GO, Panther Platform) has indicated relations between LXR and VDR signaling regulated by 20,23(OH) 2 D3 (one of the most active and promising, non-toxic vitamin D3-hydroxyderivative) in mice fibroblasts (Supplementary Table 10). Similarly to IPA analysis, GO indicated that 20,23(OH) 2 D3 is involved in enrichment of vitamin D response element binding-and vitamin D receptor binding-gene sets, as well as influencing the expression of genes involved in the regulation of metabolic transformations of vitamin D compounds. The main gene sets overrepresented by 20,23(OH) 2 D3 and involved in regulation of the LXR/RXR pathway included low-density lipoprotein particle receptor activity, lipoprotein lipase activity, fatty acid synthase activity, retinoid X receptor binding, S100 protein binding, ABC-type transmembrane transporter activity, retinoic acid-responsive element binding, low-density lipoprotein particle  (B) or upnregulated (C) protein coding genes for absolute fold change ≥ 2 cutoff. The Venn diagrams were prepared using Venny version 2.1.0: https:// bioin fogp. cnb. csic. es/ tools/ venny/ index. html. The heatmap was prepared using ClustVis: https:// biit. cs. ut. ee/ clust vis/. www.nature.com/scientificreports/ binding, high-density lipoprotein particle receptor activity, high-density lipoprotein particle binding, oxidized low-density lipoprotein particle receptor activity, oxysterol binding, cholesterol transfer activity, apolipoprotein A-I binding, peroxisome proliferator activated receptor binding, interleukin-1 type II receptor antagonist activity, hydroxymethylglutaryl-CoA synthase activity, hydroxymethylglutaryl-CoA reductase (NADPH) activity, fatty-acyl-CoA binding, apolipoprotein receptor activity, very-low-density lipoprotein particle receptor activity, farnesyl-diphosphate farnesyltransferase activity, chylomicron binding, bile acid receptor activity, nitric-oxide synthase regulator activity, cholesterol binding, sterol response element binding, nitric-oxide synthase binding and acetyl-CoA binding (Supplementary Table 10).
To test the reliability of our docking studies, we re-docked the co-crystalized ligands with LXRα and LXRβ. All the corresponding ligand pairs (docked VS co-crystalized) displayed root mean square deviation (RMSD) in the range 0.5 Å-0.7 Å (acceptable values are ≤ 2.0 Å), which demonstrated that the docking reproduced the co-crystalized (experimental) ligand poses with high precision (Supplementary Fig. S2).
Overall, all 84 compounds were predicted to bind tightly to the LXRs due to significant hydrophobic interactions and intermolecular hydrogen bond formation. Due to the hydrophobic nature and similar size of the compounds tested, they all fit well into the hydrophobic cavity of the binding site (Supplementary Figs. S3 and S4). In addition to the strong hydrophobic interactions in the central region of the ligand binding domain (LBD), two hydrogen bonding regions on both ends of the LBP contribute to the binding stability of the ligand-LXR complex. Detailed information about the hydrogen bonds and hydrophobic interaction for all 84 ligands is shown in Supplementary Tables 11-15. 3-D representations for the example ligands occupying the hydrophobic cavity of the LXRα LBD (PDBID 5AVI, 3IPQ) and LXRβ LBD (PDBID 5HJP, 1UPV) are shown in Supplementary Figs. S5, S6, S7 and S8 of the Supplementary file. Although all these compounds can bind to the LBD of LXRs, different ligands could have different binding affinities with LXRs and induce different conformational and dynamical motion changes of LXRs to have different effects on gene expression. This is supported by detailed structural, conformation and dynamical motion analyses shown in "Materials and methods". Binding thermodynamics analyses in "Materials and methods" further supports ligand binding specificity with LXRs.
Glide XP docking scores of top ranking poses for the tested ligands for each LXR conformation (Supplementary Table 16) show that novel secosteroidal, 7DHC, L3 and T3 derivatives with a full-length or short side chain have similar docking scores to the natural ligands, 20(OH)C and 22(OH)C (positive controls), for LXRs. Since T3 compounds are unstable 4 , and undergo oxidative modification and 7DHC is reduced to cholesterol by the action of Δ7reductase, we focused our subsequent studies on D3 and L3 hydroxyderivatives.
Functional studies. Activation of transcriptional functions of LXRα and LXRβ. Using the luciferase reporter gene containing the LXR-response element (LXRE), we measured the induction of transcriptional activity in CHO cells and HaCaT keratinocytes by a series of vitamin D3 and L3 compounds (Fig. 4). Figure 4A,B showed a dose-dependent activation of luciferases activity by hydroxyderivatives of D3 and L3, as well as their precursors, with EC 50 values ranging from 10 -7 to 10 -10 M. Figure 4C showed the stimulation of luciferase activity in CHO and HaCaT cells with L3-and D3-derivatives at a concentration of 10 -7 M. While all compounds tested stimulated the transcriptional activity, 1,25(OH) 2 D3, 1,20,23(OH) 3 D3, 1,20,24(OH) 3 D3 and 1,20,25(OH) 3 D3 showed the strongest stimulatory effects among the D3 derivatives. 17,20(OH)2pD, which has a 2C-side chain, displayed the weakest stimulant in CHO cells. These results are in agreement with higher docking scores for the former compounds in comparison to the short side chain derivatives (Supplementary Table 16). L3-hydroxyderivatives showed stronger stimulatory activity in HaCaT cells than in CHO cells, which indicated a degree of cell type specificity.
Investigation of LXRα and LXRβ dynamic interactions with the selected D3 and L3 derivatives using MD simulations and binding free energy analyses. The interactions of D3 hydroxyderivatives with LXR ligand binding domain (LBD) could affect secondary structure, conformation and dynamical motion changes of functionally relevant regions of LXR, further affecting LXR/RXR induced cellular activities. We performed molecular dynamics simulations to provide further insight from conformation, dynamical motion and electrostatic potential perspectives to better understand the activation mechanism of LXRs by D3 and L3 deriva-  . Data are presented as means ± SE, (n = number of assays). Analysis was done using one-way ANOVA for dose responses and the t-test: *p < 0.05, **p < 0.01, ***p < 0.001 or ****p < 0.0001 versus control (ethanol).  . Analysis was done using one-way ANOVA with significance defined as *p < 0.05, **p < 0.01, ***p < 0.001 or ****p < 0.0001.  46 . We also included 20(OH)L3 in these detailed analysis, because of its hydroxyl group at C20, and 20,22(OH) 2 L3 as representatives of a dihydroxylumisterol 20,47 . The comparison of the alignment parameters of the four D3 and two L3 ligands selected with co-crystalized ligands of LXRα (PDBID 5AVI, 3IPQ) and LXRβ (PDBID 5HJP, 1UPV) and the binding energy of the selected ligands with LXRα and LXRβ receptor are shown in Supplementary Tables 17 and 18. Based on alignment parameters and docking score, crystal structures of LXR α (PDBID 5AVI) and LXR β (PDBID 5HJP) were chosen as receptors for the selected four D3 and two L3 compounds for the further investigations. Molecules of interest and the identified LXRα and LXRβ crystal structures are listed in Table 2 Binding thermodynamics analysis. Using the equilibrated last 150 ns MD simulations trajectories, we calculated the binding free energy and binding energy components of LXRα and LXRβ with the selected four D3 derivatives and two L3 derivatives (Table 3). Table 3    1,20(OH) 2 D3 N/A 3.5 × 10 -6 ± 3.1 × 10 -6 (n = 2) 2.7 × 10 -7 ± 2.0 × 10 -7 (n = 2) 1,25(OH) 2 D3 N/A 1.1 × 10 -6 ± 1.3 × 10 -6 (n = 2) 1.4 × 10 -5 ± 1.8 × 10 -5 (n = 2) D3 7.6 × 10 -7 ± 1.0 × 10 -6 (n = 2) N/A 8.9 × 10 -7 ± 3.0 × 10 - Effects of the selected D3 and L3 derivative on conformation, secondary structure, dynamical motion and electrostatic potential of LXRα and LXRβ. Hydrogen bond analyses. While the ligand binding pocket for LXRα and LXRβ is hydrophobic, there are polar or charged residues at the two ends of the cavity 48 . Table 5 shows a hydrogen bond occupancy of no less than 40% between LXRα and D3 or L3 derivatives. More than 80% hydrogen bond occupancy was predicted between 20,22(OH) 2 L3 and THR302 of LXRα. Table 6 shows the hydrogen bond occupancy of no less than 40% between LXRβ and D3 or L3 derivatives. For LXRβ-ligand complexes, more than 80% hydrogen bond occupancy was predicted between 1,20(OH) 2 D3 and SER278, between 1,20(OH) 2 D3 and HID435, between 1,25(OH) 2 D3 and HID435, between 20 (OH)D 3 and SER278, between 20 (OH)D 3 and HID435, between 25 (OH)D 3 and HID435, between 20,22(OH) 2 L3 and THR316, and between 20,22(OH) 2 L3 and GLN438. Hydrogen bond occupancy of more than 10% between D3 or L3 derivatives and LXRα and LXRβ are shown in Supplementary Tables 25 and 26. The results showed that during the dynamic interactions procedure, in addition to the residues shown in Tables 5, 6, additional residues in LBD of LXRα and LXRβ were involved in forming hydrogen bond with D3 or L3 derivatives although the percentage of times that hydrogen bond exist over the equilibrated MD simulation trajectories was less than 40%. Supplementary Figs. S19 and S20 showed the 2D interaction map of D3 and L3 derivatives with ligand binding region of LXRα and LXRβ for the representative complex structure from the clustering analysis of the equilibrated MD simulation trajectories.
The dbscan (density-based spatial clustering of applications with noise) program 49 in Amber 14 was used for the clustering analyses of the equilibrated MD simulation trajectories, and the medoid structure in the largest cluster was chosen as representative structure. The different molecular structure of D3 or L3 derivatives leads to the different binding position and posture in the LXRα and LXRβ binding pocket, which results in different hydrophobic contacts and hydrogen bond occupancy. Thus, while the D3 and L3 derivatives studied were all shown to bind to LXRα and LXRβ (Tables 1 and 3, Fig. 5), the residues in LXRs involved with hydrogen bond formations with the ligands are predicted to be different, contributing to the observed different binding energy component, including electrostatic potential (Table 3). This provides a molecular basis to interpret the experimentally observed different binding characteristics of LXRs with D3 and L3 derivatives ( Table 1, Fig. 5) and support ligand binding specificity with LXRs.
Conformation and secondary structure analyses. The ligand-binding domain of LXRα and LXRβ is a three layered α-helical sandwich which includes two β-sheets (S1 and S2) and 12 helices (h1-h12) 48,[50][51][52] . Activation function-2 (AF-2), in which helix 12 is key, is a ligand-dependent C-terminal sequence that controls LXRα and LXRβ transcriptional activity in response to ligand binding 48,[50][51][52][53] . Conformational fluctuation and secondary structure changes of LXRα and LXRβ by interactions with D3 and L3 derivatives could directly affect the size and shape of the binding pocket and conformation of AF-2 region. This could further affect the binding affinity of D3 and L3 derivatives with LXRα and LXRβ as observed in Tables 1, 3  The root mean squared fluctuation (RMSF) analysis examined the conformation stability of LXRα and LXRβ binding to different D3 and L3 derivatives (Fig. 7B). Overall, RMSF of α-helices in LXR receptors are more stable than β-sheet and the unstructured parts of LXR receptors in the complex with D3 and L3 derivatives. The data in bar graphs (A-C) show significant differences between ligand -treated and control (ethanol treated) cells. Analysis was done using t-test: **p < 0.01, ***p < 0.001 or ****p < 0.0001 versus control (ethanol). For part A the slides were examined using a KEYENCE America BZ-X710 Fluorescence Microscope (Itasca, IL) and captured using KEYENCR BZ-X viewer (version 1.3.0.5, https:// www. keyen ce. com/ produ cts/ micro scope/ fluor escen ce-micro scope/ bz-x700/ index_ pr. jsp). The images were subsequently analyzed using the JACoP plugin (version 2.1.1, https:// image jdocu. tudor. lu/ doku. php? id= plugin: analy sis: jacop_2. 0: just_ anoth er_ coloc aliza tion_ plugin: start) for colocalization analysis 1 with ImageJ (version 1.52a, http:// imageJ. nih. gov/ ij). For part C, images were captured using an Amnis ImageStreamX Mk II Imaging Flow Cytometer (Luminex Corporation) and IDEAS software version 6.2. www.nature.com/scientificreports/ Conformational stability differences for LXRα and LXRβ were observed with the interactions with the different D3 and L3 derivatives. The overall structure of helices in the LBD of LXRα and LXRβ with the interactions with D3 and L3 derivatives was stable. Different degrees of variation in RMSFs for the residues between helices, including the AF-2 region, could directly affect binding affinity of LXRs with D3 and L3 derivatives (Tables 1 and  3, Fig. 5) and the motion of LXRs, further influencing its binding to coactivators for its functions and affecting gene expression (Figs. 1, 2, 3, 4). Analysis of secondary structure by the Define Secondary Structure of Proteins (DSSP) algorithm 54 predicted that there are not significant differences in α-helix, 3 10 -helix and β-strand secondary structure of LXRα and LXRβ between the different D3 and L3 derivatives, and the LBD of both LXRα and LXRβ is mainly formed by α-helices (Fig. 7C). The stable secondary structures are supported by the observation that there is no significant conformational fluctuation difference for α-helix and β-sheet in LXRα and LXRβ bound with different D3 and L3 derivatives (Fig. 7B). Small variations in α-helix occupancy between helices 10-12 for both LXRα and LXRβ when bound with different ligand were observed (Fig. 7C). The position of helix 12 plays a key role in the control of LXR transcriptional activity by determining the recruitment of either coactivators or corepressors 48,[50][51][52][53] . The small secondary structure difference between helices 10 and 12 for LXRα and LXRβ bound with different D3 and L3 derivatives might affect this binding with coactivators as experimentally observed (Fig. 4).

Scientific
Dynamical motion analyses. Since LXRβ has a flexible ligand-binding pocket that can accommodate structurally different ligands 48 , we performed principal component (PCA) analysis for the LBD of LXRs to calculate the distribution of the relative contribution of the first fifty PCA modes of the LXRα (Supplementary Fig. S21) and LXRβ (Supplementary Fig. S22) in the complexes of LXRs with D3 or L3 derivatives. Results showed that the  48 and helix 12 is the key helix for the AF-2 region to bind coactivators required for LXR transcriptional activity 48,[50][51][52][53] . We generated porcupine plots to show the principal dynamic motions of the helix 12/AF-2 region relative to the β-sheet/ helix 6 in the ligand-binding pocket of LXRα and LXRβ by binding with different D3 and L3 derivatives. This includes the first mode of dynamical motion (Fig. 8) and the second and the third modes of dynamical motion ( Supplementary Fig. S23). The porcupine arrows in Fig. 8 and Supplementary Fig. S23 represent the motion  www.nature.com/scientificreports/ direction and the arrow length represents the magnitude of this motion. Results indicate that there is varied motion of the helix 12/AF-2 region relative to β-sheet/helix 6 for LXRs bound with different D3 and L3 derivatives demonstrating the flexible characteristics for ligand binding pocket of LXRα and LXRβ as observed in a previous study 48 . The first three principal modes of the helix 12/AF-2 region in LXRβ showed a large amplitude of motion that is consistent with the high degree of conformational fluctuation for residues in the AF-2 region, as observed in RMSF results shown in Fig. 7B. The moderate motion amplitude of β-sheets for LXRs bound with different ligands is also consistent with the fluctuation in the conformation of β-sheets observed for residues by RMSF (Fig. 7B). In addition to the principal dynamical motion of the helix 12/AF-2 region and β-sheet/helix 6 represented with porcupine arrows (Fig. 8), the dynamical motion of the other helices and unstructured regions that contributed to the binding of LXRs with ligands are shown in the movies: Movie_LXRAlpha_125OH2D3. mpeg and Movie_LXRBeta_125OH2D3.mpeg in the supplemental materials. The movies show the dynamic interaction of 1,25(OH) 2 D3 with LXRα or LXRβ during the equilibrated last 150 ns MD simulations. In the movie of the interaction of 1,25(OH) 2 D3 with LXRβ, a large motion of the unstructured region between helix 1 and helix 3 was observed, which is consistent with the observed large conformation fluctuation for the residues between helix 1 and helix 3 in Fig. 7B (4). Together with the different conformational fluctuation (Fig. 7B), the varied motions for the LBD of LXRα and LXRβ by interactions with D3 and L3 derivatives could directly influence the size and shape of the ligand binding pocket, further resulting in the different binding affinities of ligand with LXRs, as observed in Tables 1 and 3 Table 6. Occupancy of hydrogen bond (no less than 40%) formed between LXRβ and D3 or L3 derivatives.  www.nature.com/scientificreports/ Electrostatic potential analyses. Electrostatic potential analysis of the LBD of LXRα and LXRβ bound with different D3 and L3 derivatives (Fig. 9A,B) could help to understand the effect of D3 and L3 derivatives on the electrostatic potential distribution of LBD, providing a molecular insight into the binding of D3 and L3 derivatives with LXRs and with its function (Table 1, Figs. 4 and 5) from an electrostatic point of view. It is known that the ligand-binding pocket of LXRs extends from helix 12 to the β-sheet lying between helices H6 and H7 48 and helix 12 is the key helix for the AF-2 region to bind coactivators 48,[50][51][52][53] . Varied electrostatic potential distribution in the LBD of LXRα and LXRβ bound with different D3 and L3 derivatives as observed in Fig. 9A,B, including in helix 12 and ligand binding pocket, could directly affect the electrostatic energy and polar solvation energy for the complexes of LXRs with ligands. These varied electrostatic potential distribution could contribute to the different binding free energy between LXRs with ligands, as reported in Table 3, further supporting the binding of LXRs with D3 and L3 derivatives observed experimentally (Table 1 and Fig. 5). Therefore, it is expected that different electrostatic potential distributions in or near helix 12 could affect dimerization of LXRs or binding with the coactivators required for transcriptional activity, as experimentally indicated (Figs. 4, 5, 6).

Concluding remarks
Bioinformatics analyses on deposited microarray and RNAseq data identified the LXR/RXR signaling complex as a target for D3 hydroxyderivatives. This was further substantiated by stimulation of the expression of genes downstream of LXR by secosteroids and by ChIP analysis of chromatin isolated from keratinocytes treated with 1,25(OH) 2 D3 which showed a significant stimulation of LXRα/β binding to the LXRE. The expression of genes downstream of LXR was also stimulated by D3 and L3 hydroxyderivatives. Identification of LXR α/β as receptors for secosteroidal and lumisterol derivatives was confirmed and strengthened by functional assays showing: (i) the activation of the transcriptional functions of LXRα and LXR using the luciferase reporter gene containing the LXRE in CHO and HaCaT cells, (ii) binding in LanthaScreen TR-FRET LXRα and β coactivator assays, and (iii) ligand-induced translocation of LXRα/β to the nucleus. This was further supported by molecular docking with favorable docking scores and favorable binding free energy based on molecular dynamics simulations. Also docked poses for the D3 and L3 derivatives studied aligned well with the crystal structure of ligand-LXRα and LXRβ ligand binding domain complexes and indicated that they share the same ligand binding pocket. These findings allowed us to conclude that LXRα/β are indeed receptors for secosteroids and lumisterol hydroxyderivatives. In addition, molecular docking predicted derivatives of 7DHC, 7DHP, T3 and pT as potential agonists on LXRα/β, opening new possibilities for further experimental testing. Binding thermodynamics analyses based on MD simulations demonstrated that 20,22(OH) 2 L3, 20(OH)L3, 1,20(OH) 2 D3, 1,25(OH) 2 D3, 20(OH)D3, and 25(OH)D3 favorably bind with LXRα and LXRβ, in agreement with the experimental binding assays. Different binding energy components such as van der Waals energy, electrostatics energy, and polar solvation energy were observed between the different ligands in complex with LXRα and LXRβ, which can be attributed to the different dynamic interactions of LBD of LXRs with the four D3 and two L3 derivatives studied that have the different chemical structure and atomic charge parameters. The different D3 and L3 derivatives were predicted to display different hydrogen bonding patterns with the LXRs, varied residue conformation fluctuations and dynamic motion for the LBD of LXRα or LXRβ. These changes could alter the shape, size and electrostatic potential distribution of the LBD pocket, further contributing to the different binding affinities of the D3 and L3 derivatives for LXRα or LXRβ, as experimentally observed. The varied conformation fluctuation, motion, and electrostatic potential distribution of the helix 12/AF-2 region by D3 and L3 derivatives shown as Figs. 7, 8, 9, Supplementary Figs. S21-23 and two movies, could affect the recruitment of transcriptional mediators (e.g. coactivators) and might explain the observed differences in the gene expression profiles between different derivatives. These predictions are consistent with functional and receptor binding analyses and activation of genes with LXRE, and with LXR translocation to the nucleus. However, our study has also some limitations. For example, crystallography, NMR and/or cryogenic electron microscopy studies with secosteroidal or L3 compounds are needed to define the exact nature of the ligand receptor interactions. In addition, in our future studies we plan to include ChiP-Seq and use of LXR knockouts in relation to ligand activated phenotypic activity to better define molecular and biological consequences of these interactions.
The traditional endogenous ligands for LXR are oxysterols including 22R(OH)C, 24S(OH)C, 25(OH)C, 27(OH)C and 20S(OH)C [30][31][32]55 . Considering their structural similarity, it is not surprising that hydroxyderivatives of L3 (9β,10α stereoisomer of 7DHC), which are produced enzymatically by the action of CYP11A1 or CYP27A1 [20][21][22] , are able to bind and upregulate the transcriptional activity of LXRs. Further, molecular modeling results have extended the list of potential ligands for LXR to include 7DHC (photoprecursor of L3) and short side chain derivatives of L3 and 7DHC (pL and 7DHP compounds). The significance of these findings is enhanced by the recently reported biological activity of L3-hydroxyderivatives in the skin and their detection together with L3 itself in the human epidermis and serum 10,21 . Our results predict that these compounds will exert phenotypic effect via activation of LXRs in tissues expressing these receptors. This opens up exciting possibilities for studies on the phenotypic role of novel, endogenously produced LXR ligands at the local and systemic levels and on the nature of the molecular signaling, further supplemented by studies on LXR knock-out mice. We predict that molecular signaling is depend on the cell-type specific expression of nuclear receptors for L3 or related compounds, and the actual crosstalk between LXRs and RORs 10,21,23,36 .
The VDR has long been considered to be the major if not sole nuclear receptor for 1,25(OH) 2 D3 9 . This study not only breaks this dogma but also identifies its precursors and CYP11A1-derived D3-hydroxyderivatives as ligands for LXRs, which is highly significant. The identification of hydroxyderivatives of D3 as ligands for LXRs is unexpected, because of their major structural differences with sterols represented by their 9,10 secosteroidal configuration opposed to the intact B ring in 7DHC and cholesterol 2,46,56 . On the other hand, the binding pocket of LXR is large and potentially other oxysterols could fit rather easily and interact with different regions in the LXR binding pocket. Nevertheless, the A, C and D rings of D3 remain intact and the side chain is identical to cholesterol, differing only with respect to the position(s) of hydroxylation. The ligand binding pocket (LBP) of LXR is lined with mostly non-polar residues and is flexible enough to accommodate D3-hydroxyderivatives to allow the interaction of D3 hydroxyderivatives with LXRs, as predicted by molecular modelling and shown by experimental results. With respect to the VDR, the CYP27B1 mediated addition of an OH to C1α enhances selectivity of classical and non-classical hydroxy D3 compounds towards it 2,9,26,27 . However, this 1α-hydroxylation does not affect their functions as agonists on LXR, except C1α(OH)-derivatives of D3 become inverse agonists on LXRα but not β. This could be due to the different conformational fluctuation, dynamic motion and electrostatic potentials at the AF-2 region, particularly for key helices 12, in LXRα and LXRβ by binding D3 derivatives, which could result in the different transcription activities for LXRα and LXRβ. LXRs binds RXR to form heterodimer for transcriptional activities [30][31][32] . Conformational changes of LXRs by binding with different D3 derivatives could result in the ligand dependent conformational changes of the LXR-RXR heterodimer, which could create a binding surface for interaction with either a co-activator or a co-repressor peptide. Whether LXR-RXR heterodimer recruits a co-activator/repressor peptide also depends on the tissue type and the abundance of the co-activator/ Scientific Reports | (2021) 11:8002 | https://doi.org/10.1038/s41598-021-87061-w www.nature.com/scientificreports/ repressor peptides. Here, we speculate that C1α(OH)-derivatives of D3 influences the protein surface to recruit a co-repressor peptide rather than co-activator peptide. Future crystallography studies, will provide further insight into this phenomenon. Additionally, the identification of tachysterols as potential ligands by molecular modeling opens exciting areas of research for their signaling and functional activities in the human epidermis 10 . Since D3 is produced and activated in the epidermis 2,9,10 , the local concentrations of different D3-hydroxyderivatives will be very high promoting their action on multiple nuclear receptors including VDR, LXR, RORs and AhR 10 . With respect to LXR, it is expressed in the epidermis and plays an important role in the epidermal barrier functions 40,41,57,58 . Similar considerations apply to orally delivered D3 and its metabolites, which will reach high concentration in the gastrointestinal tract (GI) and liver that will include not only production 25(OH)D3 but also of CYP11A1-derived D3-derivatives, since this enzyme is expressed in GI mucosa and by cells in the immune system 59 . The future challenge is to define intra-organ selectivity of nuclear receptors activation by different D3-hydroxyderivatives to achieve the phenotypic effects and the crosstalk between LXR and RORs and perhaps VDR and AhR 10,23,36 .
In conclusion, endogenously produced L3, D3 and most likely also 7DHC and T3 derivatives are identified as ligands for LXRs opening new paradigms on their mechanisms of action, biological functions, and opening diverse opportunities for medicinal chemistry and pharmacology relating to these compounds.

Materials and methods
Source of compounds tested. Vitamin D3, lumisterol, 1,25(OH) 2 D3, 25(OH)D3 and 20(OH)C were purchased from Sigma-Aldrich (St. Louis, MO, USA). Hydroxy-derivatives of D 3 and L 3 were enzymatically synthesized and purified as previously described 20,22,45,[60][61][62][63][64][65] . In addition, 20(OH)D3, 20(OH)L3 and 17,20(OH) 2 pD were chemically synthesized and purified as previously described 18,47,66 . Animal and primary cells protocols. All animal studies followed guidelines outlined by the NIH Guide for the Care and Use of Laboratory Animals and protocols were approved by the Institutional Animal Care and Use Committee at the NIEHS or the UAB. The study was carried out in compliance with the ARRIVE guidelines (http:// www. nc3rs. org. uk/ page. asp? id= 1357). Skin samples from newborn C57BL/6 mice were collected at the NIEHS. Dermal fibroblasts were isolated from the skin of newborn mice, 1-2 day old. Briefly, we followed the protocols described previously 67 . The skin was washed in PBS, cut into smaller pieces and incubated in 2.4 U/ ml of dispase II (Roche) overnight at 4 °C. The dermis was separated from the epidermis and incubated in 0.1% collagenase (Roche) at 4 °C overnight. Next day, it was cut into smaller pieces and incubated for 1 h in 0.25% Trypsin (CellGrow) at 37 °C. Single cell suspension was prepared by pipetting and draining through a cell strainer. After centrifugation, cell pellets were suspended in DMEM medium containing 10% charcoal-stripped fetal bovine serum (cFBS), mouse fibroblasts growth factor (FGF) and 1% antibiotic, and inoculated into tissue culture dishes (Midwest Scientific, TPP). Cells, in early passages, from different mice were pooled together to provide representative samples for experiments. Fibroblasts were cultured in Petri dishes (10 cm in diameter) in DMEM and 10% cFBS and triplicate semi-confluent cultures were treated as described.
HaCaT keratinocytes were cultured as described previously 68 . For RNA isolation the cells were maintained in TPP tissue culture petri dishes (Ø60 mm, 22.1 cm 2 ) in DMEM containing 10% cFBC to reach semi-confluence and then exposed to 10 -7 M compounds or a corresponding concentration of ethanol solvent.
Adult hairless mice SKH-1 (Male, 16 weeks old) and Ptch +/− /SKH-1 (Male, 18 weeks old) were injected subcutaneously with 20 µg and 10 µg/kg of 20(OH)D3, respectively, dissolved in ethanol, further diluted in propylene glycol (PG) and saline (0.9%) sequentially (1:19:80, v/v/v). Controls were injected with the vehicle without the drug. )After 6 h the animals were killed by halothane followed by cervical dislocation and brains were collected for RNA and protein extraction protocols. The procedures followed the guidelines and approvals of the Institute Animal Care and Use Committee of the University of Alabama at Birmingham (IACUC).

RNAseq preparation and analyses.
Triplicate cultures of murine fibroblasts were used for RNA isolation (RNAeasy Micro kit, Qiagen), the RNA samples were combined for each condition and at least 200 ng of RNA for each sample was used to create libraries for RNA sequencing by BGI Americas Corporation Service (Cambridge, MA) 67 . Quality of RNA and its concentration were determined using Agilent 2100.
Pathway analysis using ingenuity. Data was analyzed using the Ingenuity Pathway Analysis (Ingenuity Systems, www. ingen uity. com) 72 . A data set containing gene identifiers and corresponding expression values was uploaded into the application to generate networks. The identifiers were mapped to their corresponding objects in Ingenuity's Knowledge Base. A cutoff of ± 2 of fold change was used to identify molecules whose expression levels were significantly regulated. These molecules, called Network Eligible molecules, were overlaid onto a global molecular network developed from information contained in Ingenuity's Knowledge Base. Based on their connectivity Networks of Network Eligible Molecules were algorithmically generated. The Functional Analysis identify the biological functions or/and diseases that are most significant to the entire data set. Molecules from the dataset that meet the fold change cutoff of ± 2 and are associated with biological functions or/and diseases Quantitative RT-PCR. RNA isolated from either fibroblasts, HaCaT cells or brain was submitted for cDNA synthesis (High Capacity cDNA Reverse Transcription Kit with RNase Inhibitor, Applied Biosystems) following the manufacturers' protocols. RT-PCR was carried out using Cyber green, in triplicates; as previously described 73 .
Approximately, 750 ng of cDNA was used for PCR reactions. GAPDH and CIC-B were used as internal controls. Primer sequences are listed in Supplementary Table 27. Data were analyzed using GraphPad Prism statistical software one way ANOVA or t-test, where appropriate; *p < 0.05; **p < 0.01; ***p < 0.001.
ChIP assay. We followed protocols described in Ref. 40  LXR binding assay. Binding assays of graded concentrations of vitamin D3 or lumisterol hydroxyderivatives to LXRα/β were performed using the LanthaScreen TR-FRET LXRα/β Coactivator kit (Thermo Fisher Scientific, Inc., Waltham, MA) following the manufacture's protocol. Briefly, after addition of LXRα/β-LBD to the derivatives, a mixture of peptide (Fluorescein-TRAP220/DRIP-2 for LXRα or Fluorescein-D22 for LXRβ) and antibody (Tb-anti-GST) was added to the reaction followed by incubation in room temperature for 2 h. TR-FRET ratio was calculated by dividing the emission at 520 nm by the emission at 495 nm using Synergy neo2 (BioTek, Winooski, VT).
Luciferase assays. The effects of vitamin D3 and L3 hydroxyderivatives were tested on the LXRE transcriptional activity using pGreenFire1-LXRE-in-ApoE Lentivirus (System Biosciences, Palo Alto, CA). This plasmid was transduced in CHO and HaCaT cells following the manufacture's protocol and the cells used as a stable LXRE expressed tool. Briefly, the virus was transduced into the cells using TransDux MAX Lentivirus Transduction Enhancer (System Biosciences, Palo Alto, CA) and screened to select the transduced cells using puromycin. The transduced cells were treated with the compounds for 24 h in media containing 5% cFBS (charcoal-treated fetal bovine serum) and the luminescence was read by ONE-Glo Luciferase Assay System (Promega, Madison, WI) using Cytation 5 (BioTek, Winooski, VT).

Ligand induced LXR translocation to the nucleus.
HaCaT cells were plated in DMEM plus 5% cFBS either on a sterile glass cover slip placed in a well of a 6-well plates for immunofluorescence (IF) assays or TPP Petri dishes for image-stream flow cytometry, and treated with 10 -7 M D3-or L3-hydroxyderivatives or 0.1% ethanol (vehicle) as a negative control was added to each well, and the cells were incubated for 12 or 24 h at 37 °C with 5% CO 2 . For IF the cells were fixed in 4% paraformaldehyde for 15 min, incubated in 0.5% Triton X-100 in phosphate-buffered saline (PBS) for 5 min, and washed 3 times in PBS. After blocking with 2% bovine serum albumin (BSA) in PBS for 1 h at room temperature, the cells were treated with an anti-LXRα antibody (Novus Biologicals, LLC, Centennial, CO) at 1:200 in PBS supplemented with 2% BSA overnight in 4 °C. After washing the slides were incubated with a mouse anti-rabbit IgG-FITC (Santa Cruz Biotechnology, Inc, Dallas, TX) at 1:200 dilution for 1.5 h at 37 °C. After washing the slides were mounted using propidium iodide (PI) (Vector Labs, Burlingame, California) as a nuclear counterstain. The slides were examined using a KEYENCE America BZ-X710 Fluorescence Microscope (Itasca, IL). The images were subsequently analyzed using the JACoP plugin for colocalization analysis with ImageJ 74 . Five separate high-power field (40×) images were converted to grayscale and analyzed for each condition. Manders' coefficient (M1), calculated using a preset threshold, was used to quantify the degree of colocalization between PI and LXRα with 0 representing no colocalization, and 1 representing perfect colocalization 75 .
For Image Stream II (Amnis, Seattle, WA, USA) cytometer analyses 76 , HaCaT cells were detached and processed as previously described 77 . The cells were fixed with paraformaldehyde and permeabilized in methanol containing buffer 78,79 and Hoechst stain (r37165, ThermoFisher) and rabbit antibodies to LXRα/β (sc 377260, AF 647, SantaCruz) or VDR (sc-13133, AF 488, SantaCruz) as described previously 76 . Data were analyzed using IDEAS software (Amnis, Seattle, WA, USA). www.nature.com/scientificreports/ LXRβ, a receptor-based approach (molecular docking) was used. Glide in Extra Precision (XP) implemented in Schrödinger (version 2016-1) was used for the docking studies. Taking into account alternative ligand binding modes that could come out of exploring different conformations of the active site, among LXR crystal structures in Protein Data Bank 80 , two LXRα (PDBID 5AVI, 3IPQ) and three LXRβ (PDBID 5HJP, 1PQC, 1UPV) with most structurally diverse conformations were selected for the docking studies, respectively, by structural clustering and visual inspection. The chosen five LXR conformations were energy minimized with OPLS2005 force field. 3D structures of the 84 ligands tested which included a D3 series (full-length side chain), pregnacalciferol (pD)(short side chain), L3, pregnalumisterol (pL) series, 7DHCseries, 7DHP series, T3 and pregnatachysterol (pT) series and natural ligands (controls) including 20(OH)C and 22(OH)C were prepared using LigPrep utility of Schrödinger at pH 7.0 and with OPLS2005 force field. Chiralities were taken into account and assigned correctly based on the chemical structure of each ligand. The 84 ligands were separately docked into each energy minimized structure of the ligand binding domain (LBD) of LXRα and LXRβ. During molecular docking, receptor grids were constructed without assigning any constraints, such as hydrogen bonds or excluded volumes, in order to allow the docking algorithm to explore maximum positions and orientations of the docked ligands. The ligand structures were set as fully flexible while the proteins remained rigid; however, rotatable groups in the residues that could possibly make contacts with the docked ligands were allowed. Docking scores were determined to evaluate the binding of each compound tested with LXRα and LXRβ (see Supplementary Information).

Molecular dynamics (MD) simulation and binding free energy analyses.
With the selected four D3 and two L3 derivative compounds in complex with LXRα or LXRβ, we performed a total of twelve 300 ns MD simulations using the Amber 14 MD simulation package (https:// amber md. org/) 81,82 to determine the binding mode and binding free energy of LXRs with the selected D3 and L3 derivatives, and determine the effect of the D3 and L3 derivative on conformation, secondary structure, dynamical motion and electrostatic potential distribution of LXRα and LXRβ. AMBER force field was used for the simulated systems. A standard MD simulation protocol was performed, as used in our previous studies [83][84][85][86][87][88][89][90][91][92][93] for twelve simulated systems as shown in Table 2: LXRα and LXRβ with four D3 and two L3 derivatives separately. Binding free energies were calculated for the twelve complexes with the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method as described previously 87,[94][95][96] which was implemented using the AMBER 14 MD Software Package.
To determine the simulated systems' equilibration tendencies and its convergence, root mean square deviation (RMSD) of protein backbone atoms over time and binding free energy of the complex structure over the cumulative time were analyzed. Based on equilibrated MD simulation trajectories, we performed a series of analyses to better understand hydrogen bond formation, conformational stability, secondary structure and dynamical motion and electrostatic potential characteristic of LXRα and LXRβ binding with the selected D3 or L3 derivatives, which were performed using the ptraj program of AMBER 14. Electrostatic potential analyses were performed using APBS (v1.2.0, http:// www. poiss onbol tzmann. org/) 97,98 . Details for determining ligand force field including the charge parameters of ligands are described in Supplementary Information section on Molecular Modelling. We also detailed MD simulation protocols, MM-PBSA binding free energy method and the analyses of hydrogen bond formation, conformation, secondary structure, dynamical motion and electrostatic potentials in Supplementary Information: Molecular Modelling.