Unravelling the Role Played by Non-Covalent Interactions in the Action Mechanism of PCDDs within Cells

Unravelling the Role Played by Non-Covalent Interactions in the Action Mechanism of PCDDs within Cells Lorena Ruano,†,‡Álvaro Pérez-Barcia,‡Vito F. Palmisano,†,¶Juan J. Nogueira,∗,†,§Marcos Mandado,∗,‡and Nicolás Ramos-Berdullas∗,‡ †Department of Chemistry, Universidad Autónoma de Madrid, 28049, Madrid, Spain. ‡Department of Physical Chemistry, University of Vigo, Lagoas-Marcosende s\n, ES-36310-Vigo, Galicia, Spain ¶Theoretical Chemistry Group, Zernike Institute for Advanced Materials, University of Groningen, Groningen, The Netherlands. §Institute for Advanced Research in Chemistry (IAdChem), Universidad Autónoma de Madrid, 28049 Madrid, Spain E-mail: [email protected]; [email protected]; [email protected] Abstract The aryl hydrocarbon receptor (AhR) is a ligand-activated transcription factor that mediates biological signals and regulates diverse cellular functions. Of particular concern are the effects triggered by dioxins and dioxin-like compounds (DLCs), whose toxicological outcomes arise through both canonical and non-canonical pathways, leading to the designation of AhR as the "dioxin receptor". However, conventional risk assessment approaches based on toxic equivalency factors (TEFs), which primarily reflect the capacity of these compounds to bind and activate AhR, do not fully account 1 for critical aspects such as environmental concentration and bioavailability, potentially underestimating their true impact. In this work, we present a comparative analysis of polychlorinated dibenzo-p-dioxins (PCDDs) with varying degrees of chlorination, focusing on their interactions with the AhR at the ligand-binding domain and on their permeation abilities across a model lipid membrane. To this end, we combine classical molecular dynamics (CMD) simulations with a hybrid quantum mechanics/molecular mechanics energy decomposition analysis (QM/MM-EDA) framework. This integrated approach enables a molecular-level characterization of receptor binding affinities and membrane permeation efficiencies. Our findings provide novel insights into the mechanisms underlying the relative toxicity of DLCs and highlight the need for integrative assessment strategies that encompass both receptor-ligand interactions and physicochemical behavior in biological environments. It is noteworthy that the toxicity of these compounds, as quantified by the pEC50 index, correlates with the membrane permeation barrier rather than with AhR binding affinity, identifying permeation as the key mechanistic step in the toxicological process of these compounds. Introduction The aryl hydrocarbon receptor (AhR) is a cytosolic transcription factor characterized by ligand-dependent basic helix-loop-helix (bHLH) and Per-Arnt-Sim (PAS) domains.1The AhR exhibits affinity for a broad spectrum of structurally diverse exogenous and endogenous ligands, including both agonists and antagonists.2,3 Upon activation, the AhR regulates the expression of numerous genes, including those involved in xenobiotic metabolism and others that influence key physiological processes4,5 and pathological conditions such as cancer.6,7 Specifically, the binding of dioxin-like compounds (DLCs), including polychlorinated dibenzo-p-dioxins (PCDDs), polychlorinated dibenzofurans (PCDFs), and polychlorinated biphenyls (PCBs),8,9 to the PAS-B domain of AhR initiates the dissociation of chaperone proteins and the formation of a heterodimer with the AhR nuclear translocator (ARNT). This 2 complex translocates to the nucleus, where it binds to xenobiotic response elements (XREs) in deoxyribonucleic acid (DNA), promoting the transcription of metabolizing enzymes. This signaling cascade is widely recognized as the canonical AhR activation pathway.10 Despite this relatively detailed mechanistic understanding and nearly fifty years of research following the discovery of the high-affinity binding between the most toxic dioxin,11 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), and AhR, which is often referred to as "the dioxin receptor", the scientific community continues to investigate the toxicological mechanisms and long-term health effects associated with chronic exposure to DLCs.12–15 Indeed, emerging evidence suggests that AhR signaling is not restricted to this canonical mechanism. Non-canonical pathways have been reported, varying by ligand structure, cellular context and environmental conditions.16,17 Moreover, several ligands are capable of eliciting toxic responses via AhR-independent mechanisms.18 These findings challenge the prevailing assumption that AhR activation alone fully accounts for the toxicity of DLCs and underscore the need for a broader understanding of the mode of action of the pollutants.19 In addition to receptor-mediated effects, the toxicological impact of dioxins is significantly influenced by their physicochemical properties as persistent organic pollutants (POPs).9,20 Routine toxicity assessment of DLCs is based on the toxic equivalency factor (TEF), which reflects their relative potency compared to TCDD and accounts only for AhR-mediated mechanisms. However, many DLCs with lower TEF values than TCDD are released into the environment in significantly larger quantities. Their high hydrophobicity, chemical stability and resistance to degradation lead to bioaccumulation, particularly in adipose tissues, facilitate biomagnification through the food chain and contribute to long-term toxic burdens in exposed organisms.21,22 Consequently, relying solely on TEF values may underestimate the actual risk posed by these compounds, as TEFs do not capture the effects of concentration and distribution.8Therefore, a more comprehensive assessment of dioxin toxicity should integrate both receptor binding and physicochemical behaviors, including absorption, membrane transport and context-specific potency estimates beyond the TEF framework.23 3 The initial events involving both membrane transport and receptor engagement that underlie the mechanism of action of DLCs have rarely been investigated in detail at the molecular level. Early studies using molecular docking techniques provided valuable structural insights into the ligand-binding domain (LBD) of AhR and identified key amino acid residues involved in TCDD binding.24–26 Advanced computational tools to capture receptor conformational flexibility, such as classical molecular dynamics (CMD) combined with enhanced sampling or metadynamics techniques have been used almost exclusively on TCDD in order to calculate binding energies and identify potential ligand-binding pathways to the human AhR.27–29 Similarly, studies of membrane absorption and diffusion processes have primarily used CMD simulations.30–32 In both receptor and membrane contexts, these classical methods previously employed lack the precision necessary to characterize the interactions with high accuracy, for which a quantum mechanical description would be beneficial. To address these limitations, hybrid simulation approaches that integrate quantum mechanics (QM) for the interaction region with molecular mechanics (MM) for the surrounding environment have proven effective in accurately describing non-covalent interactions in complex biological systems. In particular, hybrid energy decomposition analysis (EDA) frameworks employing QM/MM calculations allow interaction energies to be partitioned into fundamental contributions: electrostatics, Pauli repulsion, induction and dispersion. Recent studies by the authors have successfully used this QM/MM-EDA approach to investigate the nature of intermolecular interactions in biologically relevant systems.33,34 Notably, one such study focused on the absorption and diffusion of two highly toxic compounds, TCDD and 2,3,7,8-tetrachlorodibenzofuran (TCDF), through lipid membranes.35 In the present study, a comparative analysis of PCDDs (see Figure 1) with different levels of chlorination was conducted and their interactions with AhR (see Figure 2) and their permeation through a model lipid membrane were examined. To this end, a combination of CMD simulations and QM/MM-EDA was used to characterize these processes at the molecular level, complemented by the classical molecular mechanics generalized Born surface 4 area (MM/GBSA) approach used for comparative analysis. The goal was to provide, for the first time, a deeper understanding, supported by accurate quantum chemical results, of the influence of PCDDs physicochemical properties on both receptor binding and membrane permeation, thereby offering new insights into the factors underlying their relative toxicity. 0 O OCl Cl Cl Cl 0 O OCl Cl Cl Cl Cl 0 O OCl Cl Cl Cl Cl Cl m3 m2 m1 0 O OCl Cl Cl Cl Cl Cl m4 0 O OCl Cl Cl Cl Cl m6 Cl 0 O OCl Cl Cl Cl Cl Cl m5 Cl 0 O OCl Cl Cl Cl Cl Cl m7 Cl Cl Figure 1: Polychlorinated dibenzo-p-dioxins (PCDDs) referred to as m1 through m7. Figure 2: AhR complex in blue with the chaperone Hsp90 and the co-chaperone XAP2 in green. PDB ID: 7ZUB. 5 Methodology and Computational Details Prior to the CMD simulations, the structures of the PCDDs (see Figure 1), hereafter referred to as m1 through m7, were optimized using density functional theory (DFT) with the B3LYP36–39 functional and the cc-pVDZ40 basis set. Additionally, Merz–Singh–Kollman charges were computed at the M06-2X41/cc-pVDZ40 level for subsequent use in the electrostatic component of the force field. CMD simulations were then performed with Amber2042 with the ligand integrated into the protein and the lipid membrane (see details below). After that, the interaction energy between the ligands and the protein and the membrane were characterized by two different methodologies: by the classical approach MM/GBSA and by a QM/MM-EDA scheme. Classical Molecular Dynamics Simulations Lipid Membrane A dioleoylphosphatidylcholine (DOPC) lipid bilayer composed of 64 DOPC molecules per layer was constructed using the CHARMM-GUI43 platform. A 22.5 Å-thick water layer was added above and below the bilayer to solvate the system. To mimic physiological conditions, a 0.15 M KCl concentration was included. The membrane was described using the Lipid17 force field, an updated version of the earlier Lipid1144 and Lipid1445 force fields. The TIP3P model46 was used for water molecules, while K+and Cl−ions were described using appropriate Amber parameters.47 The system underwent a multi-step relaxation protocol. Initially, 10000 steps of energy minimization were performed, transitioning from steepest descent to conjugate gradient halfway the simulation. This was followed by a two-step heating phase: first, a 5 ps NVT simulation increased the temperature to 100 K using a Langevin thermostat with a collision frequency of 1 ps−1; second, a 100 ps NPT simulation further heated the system to 303 K, applying a Berendsen barostat (2 ps pressure relaxation time) to maintain the pressure at 1 bar alongside the Langevin thermostat. During heating, 6 a positional restraint of 10 kcal/molÅ2was applied to the lipid molecules. Subsequently, a 5 ns NPT equilibration run were performed to stabilize the periodic box dimensions. A 125 ns production simulation followed under the same NPT conditions as in the second heating step. Finally, a randomly selected snapshot from the production trajectory was extracted for use as the final system configuration. After the equilibration of the lipid membrane, each PCDD molecule was positioned 32 Å above the center of the relaxed DOPC membrane. Force field parameters for the PCDDs were derived from the general Amber force field48 (GAFF2), with Merz–Singh–Kollman charges computed at the M06-2X41/cc-pVDZ40 level of theory with Gaussian1649 software. The system then underwent a 20000-step energy minimization, transitioning from a steepest descent algorithm to conjugate gradient halfway the simulation. Subsequently, the system was heated to 303 K under NVT conditions using a Langevin thermostat with a collision frequency of 1 ps−1. Positional restraints of 10 kcal/molÅ2and 5 kcal/molÅ2were applied to the lipid bilayer and the PCDD molecule, respectively, during this 1 ns heating phase. Following heating, four consecutive 5 ns equilibration runs were performed in the NPT ensemble to prepare the system for umbrella sampling simulations. These simulations employed a Berendsen barostat with a pressure relaxation time of 1 ps and the Langevin thermostat with a collision frequency of 1 ps−1to maintain a pressure of 1 bar and a temperature of 303 K, respectively. During the equilibration process, positional restraints of 10.0, 5.0, 2.5 and 1.0 kcal/molÅ2on lipids and 5.0, 2.5, 1.0 and 1.0 kcal/molÅ2on PCDDs were applied respectively. After equilibration, each ligand was steered along the z-axis, defined as the reaction coordinate, toward the center of the bilayer using umbrella sampling simulations lasting 32 ns. A pulling force constant of 1.1 kcal/molÅ2was applied, while a restraint force constant of 2.5 kcal/molÅ2maintained each PCDD within its respective sampling window. The resulting trajectory was divided into 32 windows along the reaction coordinate. For each window, a 30 ns CMD simulation was performed under the same conditions as the pulling step. All 7 simulations employed a 2 fs time step, with bonds involving hydrogen atoms constrained using the SHAKE50 algorithm. Long-range electrostatic interactions were treated with the particle-mesh Ewald51 (PME) method with a cutoff of 10 Å. At the bilayer center of each system, the total binding free energy (∆Gtot) was computed using the MM/GBSA method. The decomposition of ∆Gtot into van der Waals (vdW), electrostatic, polar solvation, and non-polar solvation components was also evaluated. The first 5 ns of each trajectory were discarded for analysis. A salt concentration of 0.15 M was used and solvation radii were assigned using the mbondi2 radii set (igb=5), as developed by Onufriev et al.52 Aryl Hydrocarbon Receptor The crystal structure of the protein with PDB ID: 7ZUB was obtained from the Protein Data Bank. The AhR corresponds to chain D (residues 286–399) of this structure and was used for all subsequent simulations (see Figure 2). Rigid molecular docking was performed to determine the optimal binding pose of the ligands within the receptor’s binding site.10,53 A grid size of 16 ×16 ×16 Å3with a grid spacing of 1 Å was used. The pose with the most favorable free energy was then solvated in a truncated octahedral water box, ensuring a minimum distance of 14 Å between any protein atom and the box edge and Na+and Cl− ions were added using the tleap module of AmberTools20.42 The system components were described as follows: the protein with the ff19SB force field,54 the ligand with GAFF2,48 the water molecules using the TIP3P model,46 and ions with parameters developed by Joung and Cheatham.47 Then, in order to equilibrate the system structure and density, CMD simulations were carried out. Energy minimization was first performed using 5000 steps conducted by the steepest descent approach, followed by 5000 steps using the conjugate gradient method. The system was then gradually heated from 0 K to 303.15 K over 40 ps in the NVT ensemble, using a Langevin thermostat with a collision frequency of 1 ps−1. Finally, a 500 ns production simulation was run in the NPT ensemble, applying a Berendsen barostat to maintain the 8 pressure at 1 bar. A 2 fs time step was used throughout all stages. Electrostatic interactions were treated with the PME51 method, using a 12 Å cutoff. All bonds involving hydrogen atoms were constrained using the SHAKE algorithm.50 The total binding free energy of each of the ligands with the protein was calculated using the MM/GBSA approach. After an initial 20 ns period in which the protein structure stabilized, as assessed by root mean square deviation (RMSD) analysis, 200 evenly spaced snapshots were extracted from the remaining 480 ns of the production trajectory. These geometries were used to compute ∆Gtot and its decomposition into vdW, electrostatic, polar solvation, and non-polar solvation components. A salt concentration of 0.15 M was used, and mbondi2 solvation radii were applied (igb = 5) following the model developed by Onufriev et al.52 Additionally, a pairwise ligand–residue energy decomposition analysis was performed to identify individual contributions from protein residues within 5 Å of each ligand’s center of mass. This analysis provided detailed insight into the energetic contributions of specific amino acids to ligand binding. Quantum Mechanics/Molecular Mechanics Energy Decomposition Analysis In order to properly account for quantum effects, such as Pauli repulsion and dispersion forces, QM/MM calculations were performed with Gaussian1649 on a selected set of geometries, equally spaced along the CMD simulations. Specific details regarding the geometries selected for both membrane permeation and AhR binding are provided below. Subsequently, an EDA of the intermolecular interaction energies was carried out on these geometries, along with a statistical analysis of each energy component. The EDA scheme employed in this study follows the methodology developed by Mandado et al.,55–57 which is based on the fragmentation of the complex’s deformation density into Pauli and polarisation contributions. Specifically, its extension to the QM/MM level, recently implemented in the EDANCI program57,58 was used. In this methodology, the different energy components include 9 over a set of representative equally spaced geometries sampled at the center of the membrane. Figure S1 in the supplementary material shows the convergence analysis of the interaction energy and its components across the set of geometries and residues. The mean value of the total interaction energy in m1 is -31.93 kcal/mol, progressively decreasing (i.e., becoming more stabilizing) with increasing chlorine substitution, reaching -37.80 kcal/mol for m7. These interaction energies correlate quite well with the classical MM/GBSA free energies, as shown in Figure 6A. 47 44 41 38 G total ( MM / GBSA ) (kcal/mol) 37 35 33 31 E int (QM/MM) r 2 =0.941 A 44 42 40 38 E int ( MM / GBSA ) (kcal/mol) 37 35 33 31 E int (QM/MM) r 2 =0.951 B Figure 6: QM/MM total interaction energies vs A) MM/GBSA ligand/protein binding free energy and B) MM/GBSA total interaction energy. However, the explicit inclusion of QM effects, which are absent from classical force fields, reveals important discrepancies in how specific non-covalent interactions are represented. Notably, while polarisation dominates the interaction energy, the electrostatic component emerges as a non-negligible contributor to the stabilization of the pollutants in the membrane. In order to compare the electrostatic energy contribution against the non-electrostatic one (vdW), as it was done in the analysis of the classical energies, the vdW QM energy is computed as the sum of the Pauli and polarisation terms, which is analogous to the vdW energy provided by the classical Lennard-Jones potential from the force field. When compared with the QM vdW energy, the electrostatic QM contribution shows a comparable or even 16 slightly greater magnitude, particularly in the highly chlorinated species m6 and m7. This stands in contrast to the classical MM/GBSA calculations, which largely underestimate the electrostatic contribution and, in the case of m7, even predict a slightly destabilizing effect. Additionally, dispersion energy constitutes the most stabilizing contribution, accounting for approximately 69-70% of the total stabilization energy, followed by electrostatic (2122%) and induction contributions (8-9%). These relative proportions remain quite invariant across the PCDDs series. As anticipated, all energetic components increase in absolute magnitude with progressive chlorination. Consequently, the enhanced energetic stability observed with increasing chlorine substitution is primarily due to a smaller relative increase in Pauli repulsion compared to the overall rise in attractive interactions, rather than to the amplification of any specific attractive component. AhR Binding In the previous section, it was shown that the lower the free energy barrier for membrane permeation, the higher the toxicity of a given PCDD. Both classical and quantum results initially suggest that membrane translocation is the most critical step in determining toxicity within the studied series. These findings also indicate that translocation across the membrane is a thermodynamically unfavorable process, meaning that only a small fraction of PCDD molecules reach the intracellular space to interact with the AhR. The next question, therefore, concerns how differences in AhR binding affinity among the studied PCCDs influence their toxicological activity. This factor is anticipated to be critical, given that AhR activation requires ligand binding to the PCDD molecule. As a first step, an MM/GBSA analysis was conducted to see if the binding energies correlate with the toxicity and to identify the key residues of AhR contributing to the binding free energy. As shown in Figure 7, phenylalanine at position 295 (PHE295) emerged as the most important interacting residue across all PCDDs, followed by another phenylalanine, PHE351. Eight additional residues exhibiting notable interactions with the pollutants are 17 also presented in the figure. Notably, the MM/GBSA results align well with the AhR receptor cavity previously described in the literature for planar ligands,53 with all but PRO297 appearing among the top ten interacting residues. These ten residues have been selected for inclusion in the QM region of the QM/MM-EDA calculations as described above. m1 m2 m3 m4 m5 m6 m7 Ligands PHE295 PHE351 ILE325 HIE291 LEU353 CYS333 TYR322 HIE337 VAL381 ILE349 Residues 10 8 6 4 2 0 G total (kcal/mol) Figure 7: Pairwise residue decomposition of the ligand/protein binding free energy. The free energy values, along with their decomposition into electrostatic, vdW, and solvation components using the MM/GBSA method, are presented in Figure 8A and Table S3. As observed in the membrane calculations, the electrostatic contribution appears to be significantly underestimated, contributing only marginally to the total pollutant/protein interaction energy. In contrast to the membrane results, however, the solvation component here is repulsive, resulting in total free energy values that are very similar to those obtained in the membrane. Specifically, the free energy in the membrane ranges from -37.0 kcal/mol in m1 to -48.1 kcal/mol in m7, while in the protein, it ranges from -37.7 kcal/mol in m1 to -49.3 kcal/mol in m7. In both environments, the free energy becomes progressively more favorable with an increasing number of chlorine substituents. 18 A B Figure 8: Decomposition of the ligand/protein (A) MM/GBSA binding free energy and (B) QM/MM-EDA total interaction energy. These strongly stabilizing free energy values indicate that the binding of PCDDs to the AhR is highly thermodynamically favorable in all cases. This suggests that even the small fraction of molecules that permeate the cell membrane will readily bind to the receptor –or to other biological targets with higher affinity – potentially masking the specific contribution of AhR binding to the relative toxicity of different PCDDs. In addition, the finding that increased chlorination enhances binding affinity to the AhR, as indicated by more negative binding free energies, leads to an inverse correlation between binding free energies and pEC50 values. This initially seems contradictory, as receptor activation is generally assumed to require ligand binding. However, this apparent paradox may be explained by recognizing that while receptor binding is a prerequisite for activation, it is not the sole determinant in the action mechanism. Highly chlorinated ligands, despite their strong affinity, may impair the receptor’s conformational flexibility or interfere with subsequent steps in the activation cascade, such as nuclear translocation, dimerization with ARNT, or transcriptional complex assembly. These mechanistic constraints suggest that excessive chlorination may decouple binding from functional activation, thereby weakening the predictive power of affinity alone in toxicodynamic assessments.61 19 To further explore the nature of the binding between PCDDs and the AhR, QM/MMEDA calculations were performed on a series of equally spaced geometries sampled from the production phase of the CMD simulations. The mean values of the ligand/protein total interaction energy, ∆Eint, and its individual components are summarized in Figure 8B and Table S4. On the other hand, Figure S2 in the supplementary material presents the convergence analysis of the interaction energy and its components across the sampled geometries. The QM region employed for the calculations encompassed the ten amino acids that interact most strongly with the PCDD ligands, based on the MM/GBSA analysis presented in Figure 7. According to the QM/MM-EDA calculations, the total interaction energy generally increases with the number of chlorine atoms, similar to the trend observed in the membrane, though the difference between ligands is about half of that seen in the membrane. The interaction energy ranges from -21.2 kcal/mol in m1 to -24.4 kcal/mol in m6. However, m7 deviates from this trend, exhibiting an interaction energy (in absolute value) lower than that of m4,m5, and m6. This deviation may reflect an insufficiently sized QM region for this molecule, suggesting the need for a larger number of amino acids to more accurately represent the interaction at the quantum level. In fact, m7 is the pollutant that shows the strongest classical interactions with the ten amino acids selected for the QM region, as indicated by the color scale in Figure 7. This suggests that some important amino acids may be missing from the QM region, leading to discrepancies in the predicted interaction energy compared to the correct value. Consequently, the expected energy order may be altered, as the differences between pollutants are very small in this case. The interaction energy components presented in Figure 8B and Table S4 account for the differences observed in the total energies with respect to the membrane. While the electrostatic energies are similar to those calculated in the membrane, especially for the least chlorinated compounds, the polarisation energies in AhR are much less negative, and this reduction is not compensated by a corresponding decrease in Pauli energy. Dispersion 20 energy remains the most stabilizing component, but its values in AhR are notably lower than those observed in the membrane because the middle of the membrane is much more non-polar than the binding pocket of the AhR protein. The induction energy continues to be the least significant stabilizing term, and, like the electrostatic energy, the differences with the membrane are relatively small. Since the toxicity correlates well with the energy needed to cross the membrane but does not correlate with the binding energy to the AhR protein, it seems that membrane permeation is the relevant step, although binding to the protein is obviously necessary. However, this binding is very favorable for all ligands investigated here. Concluding Remarks and Future Perspectives This work provides compelling evidence that membrane translocation constitutes the central determinant of the relative toxicity of PCDDs. The progressive decrease in free energy with increasing chlorination, together with the linear correlation between permeation energies and experimental pEC50 values, supports the notion that intracellular accumulation critically regulates toxic potency. Regarding AhR binding, our results confirm that ligandreceptor association is thermodynamically highly favorable across all congeners and becomes progressively stronger with chlorination. Nevertheless, the observed inverse relationship between binding free energy (and interaction energy) and pEC50 values suggests that toxicity cannot be explained solely by receptor recognition. In this sense, highly chlorinated congeners may exhibit strong binding but limited activation, reducing the predictive power of AhR binding affinity alone. Although a reasonably good correspondence is observed between the quantum and classical total interaction energies, pronounced discrepancies are observed in the contribution of electrostatic interactions to the stabilization of PCDDs. Specifically, within both the membrane environment and the AhR binding pocket, the classical model markedly under21 estimates the electrostatic component. These discrepancies underscore the need for hybrid approaches that incorporate quantum effects when characterizing noncovalent interactions relevant to the toxicity of these compounds. Despite the computational complexity of the biological systems analyzed here, future work should aim to incorporate entropic contributions and dynamical factors into affinity assessments at the hybrid QM/MM level, as well as to extend the analysis to additional nuclear receptors, thereby providing a more comprehensive framework for predicting the toxicity of persistent halogenated pollutants. Acknowledgement This work was partially supported by the Spanish Ministry of Science and Innovation MCIN/AEI/10.13039/501100011033 through the projects PID2020-117806GA-I00, PID2022138023NB-I00 and CNS2022-135720, by the Xunta de Galica through the project GRC2024/27, and by the Comunidad de Madrid through the Attraction of Talent Program (Grant reference 2022-5A/BMD-24244). N.R.B. thanks the University of Vigo for a postdoctoral fellowship under the "Retención de Talento Investigador da Universidade de Vigo 2023" program. References (1) Schmidt, J. V.; Bradfield, C. A. Ah receptor signaling pathways. Annual Review of Cell and Developmental Biology 1996,12, 55–89. (2) Nguyen, L. P.; Bradfield, C. A. The Search for Endogenous Activators of the Aryl Hydrocarbon Receptor. Chemical Research in Toxicology 2008,21, 102–116, PMID: 18076143. (3) Denison, M. S.; Nagy, S. R. Activation of the Aryl Hydrocarbon Receptor by Struc22 turally Diverse Exogenous and Endogenous Chemicals. Annual Review of Pharmacology and Toxicology 2003,43, 309–334. (4) Larigot, L.; Juricek, L.; Dairou, J.; Coumoul, X. AhR signaling pathways and regulatory functions. Biochimie Open 2018,7, 1–9. (5) Larigot, L.; Benoit, L.; Koual, M.; Tomkiewicz, C.; Barouki, R.; Coumoul, X. Aryl Hydrocarbon Receptor and Its Diverse Ligands and Functions: An Exposome Receptor. Annual Review of Pharmacology and Toxicology 2022,62, 383–404. (6) Baker, J. R.; Sakoff, J. A.; McCluskey, A. The aryl hydrocarbon receptor (AhR) as a breast cancer drug target. Medicinal Research Reviews 2020,40, 972–1001. (7) Leclerc, D.; Staats Pires, A. C.; Guillemin, G. J.; Gilot, D. Detrimental activation of AhR pathway in cancer: an overview of therapeutic strategies. Current Opinion in Immunology 2021,70, 15–26. (8) Haws, L. C.; Su, S. H.; Harris, M.; DeVito, M. J.; Walker, N. J.; Farland, W. H.; Finley, B.; Birnbaum, L. S. Development of a Refined Database of Mammalian Relative Potency Estimates for Dioxin-like Compounds. Toxicological Sciences 2005,89, 4–30. (9) Van den Berg, M. et al. The 2005 World Health Organization Reevaluation of Human and Mammalian Toxic Equivalency Factors for Dioxins and Dioxin-Like Compounds. Toxicological Sciences 2006,93, 223–241. (10) Pandini, A.; Soshilov, A. A.; Song, Y.; Zhao, J.; Bonati, L.; Denison, M. S. Detection of the TCDD Binding-Fingerprint within the Ah Receptor Ligand Binding Domain by Structurally Driven Mutagenesis and Functional Analysis. Biochemistry 2009,48, 5972–5983, PMID: 19456125. (11) Poland, A.; Glover, E.; Kende, A. S. Stereospecific, high affinity binding of 2,3,7,8tetrachlorodibenzo-p-dioxin by hepatic cytosol. Evidence that the binding species is 23 receptor for induction of aryl hydrocarbon hydroxylase. Journal of Biological Chemistry 1976,251, 4936–4946. (12) Sorg, O.; Saurat, J.-H. Development of skin diseases following systemic exposure: example of dioxins. Frontiers in Toxicology 2023,Volume 5 - 2023. (13) Gao, J.; Xu, Y.; Zhong, T.; Yu, X.; Wang, L.; Xiao, Y.; Peng, Y.; Sun, Q. A review of food contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin and its toxicity associated with metabolic disorders. Current Research in Food Science 2023,7, 100617. (14) Hoyeck, M. P.; Matteo, G.; MacFarlane, E. M.; Perera, I.; Bruin, J. E. Persistent organic pollutants and -cell toxicity: a comprehensive review. American Journal of Physiology-Endocrinology and Metabolism 2022,322, E383–E413. (15) Torti, M. F.; Giovannoni, F.; Quintana, F. J.; García, C. C. The Aryl Hydrocarbon Receptor as a Modulator of Anti-viral Immunity. Frontiers in Immunology 2021,Volume 12 - 2021. (16) Wright, E. J.; Pereira De Castro, K.; Joshi, A. D.; Elferink, C. J. Canonical and noncanonical aryl hydrocarbon receptor signaling pathways. Current Opinion in Toxicology 2017,2, 87–92. (17) Sondermann, N. C.; Faßbender, S.; Hartung, F.; Hätälä, A. M.; Rolfes, K. M.; Vogel, C. F.; Haarmann-Stemmann, T. Functions of the aryl hydrocarbon receptor (AHR) beyond the canonical AHR/ARNT signaling pathway. Biochemical Pharmacology 2023, 208, 115371. (18) Santini, G.; Bonati, L.; Motta, S. Computational discovery of novel aryl hydrocarbon receptor modulators for psoriasis therapy. Scientific Reports 2025,15, 19963. (19) Sorg, O. AhR signalling and dioxin toxicity. Toxicology Letters 2014,230, 225–233. 24 (20) Safe, S. Polychlorinated Biphenyls (PCBs), Dibenzo-p-Dioxins (PCDDs), Dibenzofurans (PCDFs), and Related Compounds: Environmental and Mechanistic Considerations Which Support the Development of Toxic Equivalency Factors (TEFs). Critical Reviews in Toxicology 1990,21, 51–88. (21) Merrill, M. L.; Emond, C.; Kim, M. J.; Antignac, J.-P.; Bizec, B. L.; Clément, K.; Birnbaum, L. S.; Barouki, R. Toxicological Function of Adipose Tissue: Focus on Persistent Organic Pollutants. Environmental Health Perspectives 2013,121, 162–169. (22) Regnier, S. M.; Sargis, R. M. Adipocytes under assault: Environmental disruption of adipose physiology. Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease 2014,1842, 520–533. (23) Zhang, W.; Xie, H. Q.; Li, Y.; Zhou, M.; Zhou, Z.; Wang, R.; Hahn, M. E.; Zhao, B. The aryl hydrocarbon receptor: A predominant mediator for the toxicity of emerging dioxin-like compounds. Journal of Hazardous Materials 2022,426, 128084. (24) Cozzini, P.; Kellogg, G. E.; Spyrakis, F.; Abraham, D. J.; Costantino, G.; Emerson, A.; Fanelli, F.; Gohlke, H.; Kuhn, L. A.; Morris, G. M.; Orozco, M.; Pertinhez, T. A.; Rizzi, M.; Sotriffer, C. A. Target Flexibility: An Emerging Consideration in Drug Discovery and Design. Journal of Medicinal Chemistry 2008,51, 6237–6255. (25) Motto, I.; Bordogna, A.; Soshilov, A. A.; Denison, M. S.; Bonati, L. New Aryl Hydrocarbon Receptor Homology Model Targeted To Improve Docking Reliability. Journal of Chemical Information and Modeling 2011,51, 2868–2881. (26) Bordogna, A.; Pandini, A.; Bonati, L. Predicting the accuracy of protein–ligand docking on homology models. Journal of Computational Chemistry 2011,32, 81–98. (27) Casalegno, M.; Raos, G.; Sello, G. Identification of viable TCDD access pathways to human AhR PAS-B ligand binding domain. Journal of Molecular Graphics and Modelling 2021,105, 107886. 25