Molecular modeling

Molecular modeling calculations were performed on E4 Server Twin 2 × Dual Xeon-5520, equipped with two nodes. Each node: 2 × Intel® Xeon® QuadCore E5520-2.26Ghz, 36 GB RAM. The molecular modeling graphics were carried out on a personal computer equipped with Intel(R) Core (TM) i7-4790 processor and SGI Octane 2XR12000 workstations.

Analysis of structural properties of Aflatoxin B1

The experimentally determined structures of Aflatoxin B1 (CSD codes: AFLATC and AFLATM) were downloaded from the Cambridge Structural Database (CSD) using the CSDS (Cambridge Structural Database System) software Conquest 1.18. The apparent pKa values of aflatoxin were calculated using ACD/Percepta software. (ACD/Percepta software, version 2017.1.3, Advanced Chemistry Development, Inc., Toronto, ON, Canada, 2017; http://www.acdlabs.com.) The compound was considered neutral in all calculations performed because of the percentage of neutral/ionized forms computed at pH 7.4 (physiological value) using the Handerson–Hasselbalch equation. The compounds were assigned atomic potentials and partial charges using the CVFF force field (Dauber-Osguthorpe et al. 1998).

Structural and bioinformatic analysis

The experimentally determined structures of the LBD of VDR (PDB IDs: 1DB1, 1IE8, 1IE9, 1KB2, 1KB4, 1KB6, 1S0Z, 1S19, 1TXI,1YNW, 2HAM, 2HAR, 2HAS, 2HB7, 2HB8, 3A2I, 3A2J, 3A3Z, 3A40, 3A78, 3AUQ, 3AUR, 3AX8, 3AZ1, 3AZ2, 3AZ3, 3B0T, 3CS4, 3CS6, 3KPZ, 3M7R, 3OGT, 3P8X, 3TKC, 3VHW, 3W0A, 3W0C, 3W0Y, 3WGP, 4G2I, 4ITE, 4ITF, 1RJK, 1RK3, 1RKG, 1RKH, 2O4J, 2O4R, 2ZFX, 2ZL9, 2ZLA, 2ZLC, 2ZMH, 2ZMI, 2ZMJ, 2ZXM, 2ZXN, 3A2H, 3AFR, 3AUN, 3VJS, 3VJT, 3VRT, 3VRU, 3VRV, 3VRW, 3VT3, 3VT4, 3VT5, 3VT6, 3VT7, 3VT8, 3VT9, 3VTB, 3VTC, 3VTD, 3W0G, 3W0H, 3W0I, 3W0J, 3W5P, 3W5Q, 3W5R, 3W5T, 3WT5, 3WT6, 3WT7, 5XPL) and RXRα (PDB IDs: 1BY4, 1FBY, 1FM6, 1FM9, 1G1U, 1G5Y, 1K74, 1LBD, 1MV9, 1MVC, 1MZN, 1R0N, 1RDT, 1XDK, 1XLS, 1XV9, 1XVP, 2ACL, 2P1T, 2P1U, 2P1V, 2ZXZ, 2ZY0, 3DZU, 3DZY, 3E00, 3E94, 3FAL, 3FC6, 3FUG, 3H0A, 3KWY, 3NSP, 3NSQ, 3OAP, 3OZJ, 3PCU, 3R29, 3R2A, 3R5M, 3UVV, 4J5W, 4K4J, 4K6I, 4M8E, 4M8H, 4N5G, 4N8R, 4NQA, 4OC7, 4POH, 4POJ, 4PP3, 4PP5, 3A9E), were downloaded from the Protein Data Bank (PDB; http://www.rcsb.org/pdb/).

All the structures were superimposed by sequence alignment and ligand-induced protein conformational changes as well as ligand–protein interactions were analyzed (Biopolymer and Homology module of Insight 2005; Accelrys, San Diego). In particular, hydrogen atoms were added (pH of 7.2) and the interactions with all protein amino acids and water molecules having at least one atom within a 5 Å radius from any given ligand atom were monitored.

The best solved (i.e., more complete and highest resolution) structure of the LBD of hVDR (in complex with the agonist 2alpha-methyl-AMCR277B; resolution 1.45 Ǻ; PDB ID: 3A40) (Antony et al. 2010) was selected as starting protein conformation in docking studies while two very similar (Cα RMSD = 1.14 Ǻ) and high resolution (≤ 2.20 Ǻ) structures of hRXRα LBD were selected to model the full-length hRXRα LBD (see below).

The presence of RXR-based heterodimer responsive elements on the VDR promoter were predicted using the database of transcription factor binding profiles JASPAR (http://jaspar.genereg.net). In particular, the fragment from -264 to -69 of the VDR promoter was analyzed using the following position frequency matrices (PFMs): MA0074.1 (RXRα-VDR), MA0065.1 (PPARγ- RXRα), MA0115.1 (NR1H2-RXRα) MA0159.1 (RARα-RXRα), MA1146.1 (NR1H4-RXRα), MA1147.1 (NR4A2-RXRα), MA1148.1 (PPARα-RXRα) and MA1149.1 (RARα-RXRG). A relative profile score threshold of 75 was used for the selection.

Modeling of hRXRα ligand-binding domain (LBD)

As above reported, the following template structures were selected to build the molecular model of full length RXRα LBD. The X-ray structure of hRXRα LBD in complex with the agonist 3-(2'-ethoxy)-tetrahydronaphtyl cinnamic acid (Nahoum et al. 2007) (PDB ID 2P1U; resolution: 2.20 Ǻ), lacking residues 243–263, was selected as main template structure while the missing loop was modeled using the X-ray structure of the heterodimer RXRα/PPARγ (Gampe et al. 2000) (PDB ID: 1FM6; resolution 2.10 Ǻ). Ligand molecules were removed and the sequence of 2P1U and 1FM6 were aligned with hRXRα sequence downloaded from the UniProtKB/Swiss-Prot Data Bank (http://www.uniprot.org; entry P19793) by using the Multiple-Alignment algorithm (Homology module, Insight 2005, Accelrys, San Diego). Structurally conserved regions (SCR) were defined as: i) residues 229 − 242 and 264–458 of 2P1U and ii) residues 243–263 of 1FM6. The coordinates of the SCR were transferred to the hRXRα sequence by the SCR-Assign Coords procedure (Homology; Accelrys, San Diego).

The obtained homology model of RXRα LBD was completed inserting the water molecules of RXRα experimentally determined structure (PDB ID: 2P1U) through the UnMerge and Merge commands (Biopolymer module, Insight 2005, Accelrys, San Diego). Atomic potentials and partial charges were assigned using the CVFF force field. The homology model was then subjected to a total energy minimization within Insight 2005 Discover-3 module (Steepest Descent algorithm, maximum RMS derivative = 1 kcal/Å; ε = 1; Cell Multipole method for non-bond interactions (Ding et al. 1992). Only the region aa239 − 271 was left free to move during the minimization, whereas the structurally conserved regions (SCRs) of RXRα LBD were fixed to avoid unrealistic results. The final model was checked by using the Struct_Check command of the ProStat pulldown in the Homology module to verify the correctness of the geometry optimization procedure before moving to the next step. Checks included φ, ψ, χ1, χ2, χ3, and ω dihedral angles, Cα virtual torsions, and Kabsch and Sander main chain H-bond energy evaluation. The RXRα LBD homology model was used for successive dynamic docking studies.

Docking studies on human VDR and RXRα receptors in complex with Aflatoxin B1

The homology model of the full-length hRXRα LBD and the best solved structure of hVDR LBD (PDB ID: 3A40) were employed as starting protein structures in dynamic docking studies. The ligand of 3A40 was removed and atomic potentials and partial charges were assigned using the CVFF force field.

Docking studies were carried out using a Monte Carlo/Simulated Annealing (SA) docking methodology, which considers all the system flexible (Affinity, SA Docking; Insight 2005, Accelrys, San Diego, CA) (Senderowitz et al. 1995) and using the Cell Multipole method for non-bond interactions (Ding et al. 1992). Although all the system (i.e., ligand, protein, and water molecules) is perturbed by Monte Carlo and simulated annealing (SA) calculations in the subsequent dynamic docking protocol, the dynamic docking procedure formally requires a reasonable starting complex structure. To increase the variance of the starting complexes (i.e., starting ligand poses), two AFB1 starting complexes were used for each receptor, for a total of four sets of docking calculations. In particular, AFB1 was positioned: i) in hVDR LBD according to the two superimpositions on 1,25(OH)2D3 (PDB ID: 1DB1) reported in Fig. 2A; ii) in hRXRα LBD according to the two superimpositions on the 9-cis retinoic acid (PDB ID: 1FBY) reported in Fig. 2B. The binding domain area was defined as a flexible subset around the ligand constituted by all residues and water molecules having at least one atom within a 10 Å radius from any given ligand atom. The atoms included in the binding domain area were left free to move during docking calculations. A restrain buffer region was introduced to separate the freely movable atoms and non-movable atoms. If the closest distance of a movable atom to bulk atoms was less than the sum of their van der Waals radii plus the 0.5 Å, that movable atom was restrained to its original position using a harmonic restrain force of 100 kcal mol−1 Å−1.

The docking protocol included a Monte Carlo based conformational search of AFB1 within the defined active site. A Monte Carlo/minimization approach was used for the random generation of a maximum of 20 acceptable complexes. During the first step, starting from the roughly docked structures, the ligand was moved by a random combination of translation, rotation, and torsional changes to sample both the conformational space of the ligand and its orientation to the protein (MxRChange = 3 Å; MxAngChange = 180°). During this step, van der Waals (vdW) and Coulombic terms were scaled to a factor of 0.1 to avoid very severe divergences in the vdW and Coulombic energies. If the energy of a complex structure resulting from the ligand's random moves was higher by the energy tolerance parameter than the energy of the last accepted structure, it was not accepted for minimization. An energy tolerance value of 106 kcal/mol from the previous structure was used to ensure a wide variance of the input structures was successfully minimized. After the energy minimization step (conjugate gradient; 10,000 iterations; ε = 1), the energy test, with an energy range of 50 kcal/mol, and a structure similarity check (rms tolerance = 0.3 kcal/Å) was applied to select the 20 acceptable structures. Each subsequent structure was generated from the last accepted structure. The resulting docked structures were ranked by their conformational energy. Finally, to test the thermodynamic stability of the resulting docked complexes, these latter were subjected to a molecular dynamic simulated annealing protocol using the Cell_Multipole method for non-bond interactions and the dielectric constant of the water (ε = 80*r.) The protocol included 5 ps of a dynamic run divided into 50 stages (100 fs each), during which the system's temperature was linearly decreased from 500 to 300 K (Verlet velocity integrator; time step = 1.0 fs). In simulated annealing, the temperature was altered from an initial temperature to a final temperature in time increments. The temperature was changed by adjusting the kinetic energy of the structure (by rescaling the velocities of the atoms). Molecular dynamics calculations were performed using a constant temperature and constant volume (NVT) statistical ensemble and the direct velocity scaling as temperature control method (temp window = 10 K). In the first stage, initial velocities were randomly generated from the Boltzmann distribution, according to the desired temperature, while during the subsequent stages, initial velocities were generated from dynamics restart data. The temperature of 500 K was applied to surmount torsional barriers, thus allowing an unconstrained rearrangement of the "ligand" and the "protein" active site (initial vdW and Coulombic scale factors = 0.1). Successively temperature was linearly reduced to 300 K in 5 ps, and, concurrently, the vdW and Coulombic scale factors have been similarly increased from their initial values (0.1) to their final values (1.0). A final round of 105 minimization steps (ε = 80*r) followed the last dynamics steps, and the minimized structures were saved in a trajectory file. The complexes obtained by docking studies were ranked by conformational energy values and non-bond interaction energy values (vdW and electrostatic energy contribution; Group-Based method; CUT_OFF = 100; ε = 1; Discover_3 Module of Insight2005). The complex with the best compromise among these two parameters was selected as the structure representing the most probable binding mode.

In order to allow the whole relaxation of the protein, the selected complexes (hVDR and hRXRα) were then subjected to MM energy minimization without restraints (Steepest Descent algorithm; ε = 1) until the maximum RMS derivative was less than 0.1 kcal/Å (Module Discover; Insight 2005). The protein structural quality in the resulting complex was then checked using Procheck (Laskowski et al. 1993).

Ligand-induced protein conformational changes and ligand–protein interactions of the final AFB1/hVDR and AFB1/hRxRα docked complexes were analyzed and compared to those obtained by the analysis of the experimentally determined complexes as reported in the above paragraph.

Cell cultures and treatments

The Saos-2 human osteosarcoma cell line was cultured in DMEM (Gibco-Thermo Fisher Scientific, Inc. Waltham, MA, USA) supplemented with 10% (v/v) fetal bovine serum (Invitrogen-Thermo Fisher Scientific Waltham, Massachusetts, USA) at 37 °C in an atmosphere containing 5% CO2. Cells were passaged according to standard cell culture techniques. Treatments with vitamin D3 and Aflatoxin B1 were performed as follows: Saos-2 cells were plated at a density of 3 × 105cells/well in 6-well plates and were exposed for 24 h with single treatments of AFB1 (0.1 μM) and vitamin D3 (0.1 μM) or different combined treatments with a fixed dose of vitamin D3 (0.1 μM) and increasing amounts AFB1 (0.05, 0.1 and 0.2 μM). A vehicle control (0.05% DMSO) was included in each experiment. Twenty-four hours after treatments, Saos-2 cells were harvested for RNA and protein analysis.

Transient transfections and dual-luciferase reporter assays

Saos-2 cells were seeded at a density of 8 × 104 cells per well onto 12-well culture dishes and transiently transfected using Lipofectamine LTX (Invitrogen, Thermo Scientific) as previously reported (Sodaro et al. 2018a, b). A reporter plasmid containing a 960 bp fragment (-960/ + 1 nt) of the human proximal VDR promoter region was cloned upstream of the luciferase reporter gene (pVDR/Luc). Each well received 490 ng of pVDR/Luc plasmid and 10 ng of a Renilla luciferase construct (pRL-SV40, Promega, Madison, USA) as an internal control. All transfection experiments were conducted in triplicate. Aflatoxin B1 and vitamin D3 treatments were applied 4 h after transfection for each experimental point. After 24 h, cells were lysed and used for the dual-luciferase assays (Dual-Luciferase® Reporter assay system, Promega) as previously described (Sarnelli et al. 2017). All relative luciferase activities were determined by calculating the ratio of the firefly and Renilla luciferase activities, and the results are shown as mean ± SEM (n = 3). For Real-time PCR, RNAs were extracted from Saos-2 cells using Qiazol reagent (Qiagen, GmbH, Hilden, Germany) according to the manufacturer's protocol. One microgram of each RNA was reverse transcribed using QuanTitect Reverse transcription Kit (Qiagen) as reported by manufacturer's protocol and subsequently used for Real-time RT-PCR procedures on a CFX Real-time System (Bio-Rad Laboratories, Hercules, CA, USA). Real-time quantitative analysis of VDR transcripts was performed using primers as previously reported, and β actin mRNA was used as endogenous control (Faniello et al. 2009). Real-time PCR reactions were run in triplicates using the CFX96 Real-Time System (Bio-Rad Laboratories), and CT values were obtained from automated threshold analysis. Data were analyzed with the CFX Manager 3.0 software (Bio-Rad Laboratories) according to the manufacturer's specifications.

Chromatin Immunoprecipitation

Chromatin Immunoprecipitation assays were performed as described (Sodaro et al. 2018a, b). Briefly, Saos-2 cells were chemically cross-linked with 1% formaldehyde, and the reaction was stopped by adding glycine to a final concentration of 125 mM. The fixed cells were washed twice with cold phosphate-buffered saline (PBS 1X) and were lysed using a lysis buffer (5 mM PIPES; 85 mM KCl; 0,5% NP40) supplemented with a protease inhibitor cocktail (Sigma Aldrich). Nuclei were isolated and sonicated in a buffer containing 1% SDS; 10 mM EDTA; 50 mM Tris HCl pH 8.0. The resulting fragments were within the size range of 200–1,000 bp. Samples were then centrifuged at 13,000 × g for 10 min at 4 °C, and the supernatant was pre-cleared with 30 μL protein A/G PLUS-agarose beads for 2 h and incubated with 2 μg of each antibody [RXRα (D-20X) cat. no. sc-553X; RARα(C-20X) cat. no. sc-551X; VDR (C-20X) cat. no. sc-1008X; Santa Cruz Biotechnology, Dallas, TX, USA] overnight at 4 °C. Rabbit IgG antibody (sc-2027X, Santa Cruz Biotechnology) served as a negative control. Following chromatin immune-precipitation, beads were then rinsed five times with buffer A [0,1% SDS; 2 mM EDTA; 20 mM Tris HCl pH 8,0; 1% Triton X-100; 150 mM NaCl], four times with buffer B [0,1% SDS; 2 mM EDTA; 20 mM Tris HCl pH 8,0; 1% Triton X-100; 500 mM NaCl], and once with Tris–EDTA pH buffer. The bound immunocomplexes were eluted by adding 300 μL of fresh elution buffer [10 mM Tris; 1 mM EDTA pH 8.0]. Subsequently, 20 μL of 5 M NaCl was mixed with the eluted product, incubated overnight at 65 °C to reverse the cross-linking. Immunoprecipitated genomic DNA was then purified and dissolved in EB buffer (10 mM Tris; 1 mM EDTA pH 8.0) for ChIP analysis. The immunoaffinity-enriched DNA was subjected to quantitative real-time PCR analysis using SSO Advanced Universal SYBR Green Supermix by CFX96 Detection System (Bio-Rad Laboratories). The primer pairs used in the present study were as follows: VDR promoter For 5′-TCCGCACCTATAATCATCGAC-3′, VDR promoter Rev 5′-GCCACGCTGTAGCCTTAGAT- 3'; VDR enhancer S1 For 5′-CAACTGTCCCAGGCCTGAG-3′, VDR enhancer S1 Rev 5′-GGTGGGGCAACCAAGCTAA-3', HBB LCR region (used as negative control): HS2 For 5’-CCCTGTCGGGGTCAGTGCC-3', HS2 Rev 5’-CACATTCTGTCTCAGGCATCC-3'. The Ct values of specific antibodies and IgG control were normalized to the input values (∆Ct = Ct IpVDR/RXRaα/RARα or IgG-CtInput). The fold enrichment was calculated by the ∆∆Ct cycle threshold method by comparing the ChIP antibody signal to the corresponding IgG negative control (Fold enrichment = ∆∆Ct = 2^-(∆Ct IpVDR/RXRα/RARα-∆Ct IgG). Results are representative of two independent experiments.

Western blot analysis

Saos-2 cells were washed with PBS and lysed in whole-cell extract buffer (50 mM Tris–HCl, pH 8; 10% glycerol; 150 mM NaCl; 1 mM EDTA pH 8; 0.1% Nonidet P-40; and 1 mM NaF) supplemented with a protease inhibitor cocktail (complete cocktail; Sigma Aldrich, St. Louis, USA). According to the manufacturer's protocol, differential nuclear and cytoplasmic protein extracts were carried out using the NE-PER Reagents Kit (Thermo Scientific, Waltham, USA). Western blot analysis was performed as previously described (Di Caprio et al. 2015). Whole-cell extracts (30 μg) and/or differential cytosolic and nuclear extracts (15 μg) were separated by 10% SDS–polyacrylamide gel electrophoresis and electroblotted onto a nitrocellulose membrane. The membranes were then blocked with 5% non-fat milk in Tris-buffered saline for 2 h and hybridized overnight at 4 °C to an anti-VDR rabbit antibody (1:500 dilution; Santa Cruz Biotechnology, Dallas, USA #sc-1008), anti-RXRα rabbit antibody (1:500 dilution; Santa Cruz Biotechnology, #sc-553X), or anti-RARα rabbit antibody (1:500 dilution; Santa Cruz Biotechnology, # sc-551X). Following washing, the membranes were incubated with peroxidase-conjugated mouse anti-rabbit IgG (sc-2357 diluted 1: 5.000; Santa Cruz Biotechnology, Dallas, USA) secondary antibodies for 1 h at room temperature. Anti-GAPDH (1:1000 dilution; Cell Signaling, Danvers, MA, USA #21,118), anti-Vinculin (1:10,000 dilution; Abcam, Cambridge, UK #129,002), and anti-Lamin B1 (1:1000 dilution; Cell Signaling, #13,435) antibodies were used to normalized respectively whole, cytosolic and nuclear extract samples. The blots were developed using the ECL Immobilon Western Chemiluminescent HRP-substrate system (Millipore, Darmstadt, Germany) according to the manufacturer's protocol, and immunoreactive bands were detected by autoradiography according to the manufacturer's instructions or by ChemiDoc XRS Image System (Bio-Rad Laboratories). Quantification of western blots bands was performed using the ImageJ software.

Statistical analysis

All data were assessed as the mean ± standard deviation (SD) of at least three separate experiments performed in triplicate. Graphpad Prism 7 (Graphpad Software, Inc. CA, USA) was used for data analysis. Statistical differences were determined through the One-Way analysis of variance procedure followed by Dunnett's multiple comparison test, comparing results between mock control and treated cells. Differences were considered significant when p < 0.05 and highly significant when p < 0.0001. *p < 0.05, **p < 0.0001 versus mock control; #p < 0.05, ##p < 0.0001 single treatment with vitamin D3 (0.1 μM) versus combined treatment with increasing doses of AFB1.