Hormone Activity of Hydroxylated Polybrominated Diphenyl Ethers on Human Thyroid Receptor-β: In Vitro and In Silico Investigations

Background Hydroxylated polybrominated diphenyl ethers (HO-PBDEs) may disrupt thyroid hormone status because of their structural similarity to thyroid hormone. However, the molecular mechanisms of interactions with thyroid hormone receptors (TRs) are not fully understood. Objectives We investigated the interactions between HO-PBDEs and TRβ to identify critical structural features and physicochemical properties of HO-PBDEs related to their hormone activity, and to develop quantitative structure–activity relationship (QSAR) models for the thyroid hormone activity of HO-PBDEs. Methods We used the recombinant two-hybrid yeast assay to determine the hormone activities to TRβ and molecular docking to model the ligand–receptor interaction in the binding site. Based on the mechanism of action, molecular structural descriptors were computed, selected, and employed to characterize the interactions, and finally a QSAR model was constructed. The applicability domain (AD) of the model was assessed by Williams plot. Results The 18 HO-PBDEs tested exhibited significantly higher thyroid hormone activities than did PBDEs (p < 0.05). Hydrogen bonding was the characteristic interaction between HO-PBDE molecules and TRβ, and aromaticity had a negative effect on the thyroid hormone activity of HO-PBDEs. The developed QSAR model had good robustness, predictive ability, and mechanism interpretability. Conclusions Hydrogen bonding and electrostatic interactions between HO-PBDEs and TRβ are important factors governing thyroid hormone activities. The HO-PBDEs with higher ability to accept electrons tend to have weak hydrogen bonding with TRβ and lower thyroid hormone activities.


Research
Polybrominated diphenyl ethers (PBDEs) have become ubiquitous environmental pollutants because of their historical and widespread use as flame retardants, and they have received great attention from ecological health and environmental perspectives (Athanasiadou et al. 2008;He et al. 2008;Lagalante et al. 2009). Most recently, hydroxylated PBDEs (HO-PBDEs) have caused increasing concern because of reports of their natural production and metabolism (Hakk and Letcher 2003;Malmberg et al. 2005;Malmvarn et al. 2005;Mercado-Feliciano and Bigsby 2008;Wan et al. 2009). HO-PBDEs have been detected in the blood of fish, birds, and mammalian species and even in the abiotic environ ment (Athanasiadou et al. 2008;Cantón et al. 2008;Lacorte and Ikonomou 2009). Recent studies (Dingemans et al. 2008;Qiu et al. 2009) suggest that some of the toxic effects of PBDEs might be due to their HO metabolites.
In vitro tests have shown that certain PBDEs can affect thyroid hormone homeostasis by acting as potent competitors of thyroid hormones for binding to human transthyretin (TTR) and thyroid hormone receptors (TRs) (Cantón et al. 2006;Darnerud et al. 2007;Zhou et al. 2002). For instance, levels of serum thyroxine (T 4 ) were significantly decreased when rats were exposed to PBDEs (Zhou et al. 2002). The effects of PBDEs on T 4 levels may require metabolic activation because HO-PBDEs, but not the PBDE congeners themselves, behave as ligands for human TTR in vitro (Meerts et al. 2000;Qiu et al. 2007). Kitamura et al. (2008) also reported that 4-OH-BDE-90 and 3-OH-BDE-47 markedly inhibited the binding of triiodo thyronine (T 3 ) to TRα and acted as thyroid hormone-like agents. Recently, Kojima et al. (2009) reported that 4-HO-BDE-90 significantly inhibited TRα-and TRβ-mediated transcriptional activity induced by T 3 . Consequently, HO-PBDEs have attracted great attention (Routti et al. 2009). However, there is lack of systematic investigation into the mechanisms by which HO-PBDEs interfere with hormonal actions (Zhao et al. 2008).
There are mainly two subtypes of TRs, TRα and TRβ, expressed from two different genes. TRα mediates the effects of thyroid hormones on the heart, in particular on heart rate and rhythm, whereas most actions of the hormones on the liver and other tissues are mediated more through TRβ (Forrest and Vennstrom 2000). The initial step for chemical mode of action is binding to an intra cellular receptor (Kavlock et al. 1996). Given the large number of compounds that may bind to the receptors, there is increasing interest in developing computational methods (in silico) to predict affinity of compounds with the receptors, including quantitative structure-activity relationships (QSARs) (Du et al. 2008;Valadares et al. 2007). Furthermore, molecular docking and virtual screening have become an integral part of many modern structure-based computational simulations of chemicals (Martinez et al. 2008). Docking methodologies use the knowledge of three-dimensional structure of a receptor in an attempt to optimize the bound ligand or a series of molecules into the active site. Combinational use of docking with QSAR can provide more information on the interaction between the ligand and the receptor (Sippl 2002;Soderholm et al. 2005).
In this study, we applied an integrated in vitro and in silico approach to evaluate the thyroid hormone-disrupting potency of HO-PBDEs and to identify critical structural elements and physico chemical properties of HO-PBDEs related to their hormone activity. The hormone activities of 18 HO-PBDEs to human TRβ were determined using the recombinant two-hybrid yeast assay. Molecular docking was performed to find the significant ligand-receptor interactions in the binding site of TRβ, which provided a better understanding of inter actions between the HO-PBDEs and TRβ. Molecular structural descriptors were computed, selected, and employed to characterize the interactions, and finally we constructed QSAR models. We also assessed the applicability domain (AD) of the QSAR model.
Recombinant two-hybrid yeast assay and statistical analysis. The recombinant twohybrid yeast system employed a yeast cell transformed with the human TRβ plasmid, coactive plasmid, and the reporter gene expressing β-galactosidase (Li et al. 2008). We examined the specificity of the yeast two-hybrid assay for TRβ ligand using DMSO (control), T 3 , and other steroid hormones. T 3 induced β-galactosidase activity, whereas 17β-estradiol, dihydro testosterone, and progesterone did not. Thus, the recombinant two-hybrid yeast assay was highly specific for TRβ ligand without cross-talk to other receptor agonists.
We performed the recombinant twohybrid yeast assay as described previously by Li et al. (2008). Briefly, yeast transformants were grown overnight at 30°C, with vigorous orbital shaking (130 rpm). For the assay, exponentially growing overnight cultures were diluted with synthetic dextrose/leucine/ tryptamine medium to an optical density at 600 nm (OD 600 ) of 0.75. All the samples were determined at least in triplicate. Each triplicate included a positive control (T 3 ) and a negative control (DMSO). Each tested chemical was serially diluted in DMSO for a total of 7-11 concentrations. Serial dilutions (5-µL steps) were combined with 995 µL medium containing 5 × 10 3 yeast cells/mL, resulting in a test culture in which the volume of DMSO did not exceed 0.5% of the total volume. For each test culture, 200 µL was transferred into a well of the 96-well plate and incubated at 30°C with vigorous orbital shaking (800 rpm) on a titer plate shaker (TITRAMAX 1000; Heidolph Instruments GmbH, Hamburg, Germany) for 2 hr, and the cell density of the culture was measured at OD 600 (GENios A-5002; Tecan Austria GmbH, Salzburg, Austria). Then, 50 µL test culture was transferred to a new 96-well plate, and after addition of 120 µL Z-buffer (21.51 g/L Na 2 HPO 4 ·12H 2 O; 6.22 g/L NaH 2 PO 4 ·2H 2 O; 0.75 g/L KCl; 0.25 g/L MgSO 4 ·7H 2 O) and 20 µL chloroform, the cultures were carefully mixed and pre incubated for 10 min at 30°C and 13,000 rpm. The enzyme reaction was started by adding 40 µL o-NPG (13.3 mM, dissolved in yeast-based buffer). The assay culture was further incubated at 30°C for 1 hr. Finally, the reactions were terminated by the addition of 100 µL sodium carbonate (1 M). The resulting absorption was meas ured at 420 nm. The β-galactosidase activity (U) was calculated according to the following equation: where U is the activity of β-galactosidase, t is the incubation duration of the enzyme reaction, V is the volume of the test culture, D is the diluting factor (6.6), OD 600 is the cell density measured at 600 nm, and OD 420 and OD´4 20 are the cell density of the enzymic reaction supernatant and the blank, respectively, measured at 420 nm. The dose-response curves for U of the tested compounds were fitted by iterative fourparameter curve fit method using SigmaPlot, version 10.0 (Systat Software Inc., Chicago, IL, USA). The concentration inducing 20% of the maximum effect (REC 20 ) value was calculated from the fitted dose-response curves. We evaluated the statistical significance of differences by analysis of variance (we considered p < 0.05 significant). Molecular docking. We adapted the CDOCKER algorithm to find the binding mode for HO-PBDEs to TRβ. CDOCKER is an implementation of a CHARMM (Chemistry at HARvard Macromolecular Mechanics)-based docking tool that has been shown to be viable (Wu et al. 2003). It has been incorporated into Discovery Studio 2.1 (Accelrys Software Inc., San Diego, CA) through the Dock Ligands (CDOCKER) protocol. We extracted the crystal structure of TRβ (Thyroid receptor beta1 in complex with a beta-selective ligand; PDB ID 1NAX) from the RCSB Protein Data Bank [RCSB (Research Collaboratory for Structural Bioinformatics) PDB; http://www.rcsb.org/ pdb]. In CDOCKER, random ligand conformations were generated from the initial ligand structure through high-temperature molecular dynamics followed by random rotations. Then, the random conformations were refined by grid-based simulated annealing, which makes the results more accurate. The CDOCKER interaction energy between the ligand and TRβ (E binding ) was then computed. The docking analysis provided insights into the inter actions between the ligands and the receptor, which facilitated the selection of appropriate molecular parameters to charac terize the interactions in the QSAR studies.
Molecular structural descriptor generation and QSAR development. We hypothe sized that the thyroid hormone activities of HO-PBDEs were dependent on a) the partition of the compounds between water and the biophase, and b) the interaction between the ligands and the receptor TRβ. Thus, 12 theoretical parameters were computed and selected to characterize the processes: the logarithm of octanol/ water partition coefficient (logK ow ), average molecular polarizability (α), molecular volume (V), dipole moment (µ), energy of the highest occupied molecular orbital (E HOMO ), energy of the lowest unoccupied molecular orbital (E LUMO ), formal charge on hydroxyl hydrogen atoms (q OH -), formal charge on hydroxyl oxygen atoms (q O -H ), formal charge on ether oxygen atoms (q O ), electrophilicity index (ω), harmonic oscillator model of aromaticity index (I A ), and the number of bromine atoms (n Br ) [see also Supplemental Material, Table S1 (doi:10.1289/ehp.0901457)]. We purposely selected logK ow to describe the partition process. The parameters V, n Br , α, and µ also partly describe partition because many of these parameters correlate with logK ow (Nguyen et al. 2005). The parameters E HOMO , E LUMO , q OH -, q O -H , q O , ω, and I A were purposely selected to describe the inter molecular electro static inter actions between the ligands and TRβ. The quantum chemical parameters E HOMO , E LUMO , q OH -, q O -H , and q O proved successful in many QSAR studies for characterizing inter molecular electro static inter actions (Colosi et al. 2006). ω measures the ability of a compound to soak up electrons. The relative binding affinity of some estrogen derivatives correlated strongly with ω (Chattaraj 2006). The aromaticity of compounds (I A ) may influence their non covalent inter actions with the receptor, and I A has been used to charac terize halogenated biphenyls (Alonso et al. 2008).
We computed logK ow values using the EPI Suite, version 4.0 (U.S. Environmental Protection Agency 2009). V (defined as the volume inside a contour of 0.001 electrons/ bohr 3 density) and the quantum chemical parameters were computed by the Gaussian 03 programs (Frisch et al. 2004). Initial geometries of HO-PBDEs were pre optimized by semi empirical PM3 Hamiltonian and then optimized by density functional theory at the B3LYP/6-311+G(d,p) level. Frequency analysis was performed on the optimized geometries to ensure the systems had no imaginary vibration frequencies. I A was calculated by DRAGON software (Todeschini and Consonni 2000).
The 18 HO-PBDEs were randomly divided into a training set (80%) and a validation set (20%), as listed in Table 1. Partial least squares (PLS) regression was performed in volume 118 | number 5 | May 2010 • Environmental Health Perspectives developing the model because PLS can analyze data with strongly collinear, noisy, and numerous predictor variables (Wold et al. 2001). We used Simca-S (version 6.0; Umetri AB, Umea, Sweden) for the PLS analysis. Simca-S uses leave-many-out cross-validation to determine the number of PLS components (A). Crossvalidation simulates how well a model predicts new data and gives a statistical fraction of the total variation of the dependent variables that can be predicted by all the extracted components (Q 2 CUM ) for the final model. The PLS analysis was performed repeatedly to eliminate redundant molecular structural parameters, as done in our previous studies (Chen et al. 2004). Model predictability was evaluated by external validation, which was characterized by the determination coefficient (R 2 ), root mean square error (RMSE), and external explained variance (Q 2 EXT ), which are defined as follows (Schüürmann et al. 2008): where y i fit is the fitted -logREC 20 value of the ith compound; ȳ is the average response value in the training set; y i and ŷ i are the observed and predicted values for the ith compound, respectively; ȳ EXT is the average response value of the validation set; n is the number of compounds in the training set; and n EXT is the number of compounds in the validation set.
We assessed the AD of the developed QSAR model using the Williams plot, that is, the plot of standardized residuals (σ) versus leverage (hat diagonal) values (h i ) (Eriksson et al. 2003). We calculated σ as follows: where y i and ˆy i are the observed and predicted values for the ith compound, respectively, and n is the number of compounds in the training set.
The h i value of a chemical in the original variable space and the warning leverage value (h*) are defined as follows: where x i is the descriptor vector of the considered compound, X is the model matrix derived from the training set descriptor values, and p is the number of predictor variables.

Results and Discussion
Thyroid hormone activity determined by recombined yeast. Based on a plot of U versus log T 3 concentrations, the maximal induction of T 3 was achieved at 1.0 × 10 -6 M [see Supplemental Material, Figure S2 (doi:10.1289/ ehp.0901457)]. From the dose-response curve, the median effective concentration value of T 3 was 1.4 × 10 -7 M, which was similar to that reported by Li et al. (2008). The 18 tested HO-PBDEs induced β-galactosidase activity in a concentration-dependent manner in the concentration range from 10 -11 to 10 -6 M (see Supplemental Material, Figure S3). Supplemental Material, Table S2 lists the determined REC 20 values for the 18 HO-PBDEs. With 2 × 10 -7 M of the tested compounds, the PBDEs showed no significant β-galactosidase activity compared with DMSO, and HO-PBDEs exhibited significant activity [p < 0.05; see Supplemental Material, Figure S4 (doi:10.1289/ehp.0901457)]. Previous studies (Hamers et al. 2008;Meerts et al. 2000) indicated that HO metabolites of PBDEs could compete with T 4 for binding to TTR and exert thyroid hormone activity. Kitamura et al. (2008) reported that 4-OH-BDE-90 and 3-OH-BDE-47 markedly inhibited the binding of T 3 to TRα and acted as thyroid hormone-like agents. Kojima et al. (2009) reported that 4-HO-BDE-90 significantly inhibited TRα-and TRβ-mediated transcriptional activity induced by T 3 . Thus, with respect to the observed thyroid hormone activity of HO-PBDEs, our results are consistent with previous findings.
Docking analysis. Figure 1 shows the docking view of T 3 and representative HO-PBDEs and PBDEs (4-OH-BDE-42, 4´-OH-  in the binding site of TRβ. At the deep end of the pocket, Arg282 and Ile275 serve as anchoring points for the ligands. The ligands also interact with the second polar region within the binding pocket, Leu341. We observed hydrogen bonding to be a characteristic inter action. As shown in Figure 1B and C, there are mainly two types of hydrogen bonds: those formed between the hydroxyl oxygen of HO-PBDEs and the hydrogen of Arg282 and Ile276, and those between the hydroxyl hydrogen of HO-PBDEs and the carbonyl oxygen of Leu341. However, for BDE-116, we could find no hydrogen bonds with the amino acid residues in TRβ. We also observed π-π inter actions between the phenyl of HO-PBDEs and Phe272, Phe442, and Phe455. The ligand-receptor binding energy (E binding ) of the 18 HO-PBDEs is listed in Table 1. As shown in Figure 2, we obtained a simple linear free energy relationship between -logREC 20 and E binding , further proving that binding to TRβ is the key step for the HO-PBDEs to exert their thyroid hormone activity. However E binding itself was not a good predictor for -logREC 20 , as indicated by the big prediction residuals for some HO-PBDEs (Figure 2). This was expected because E binding values calculated by CDOCKER differed from the real binding energy because the environmental factors (e.g., solvents, pH, and ions) were not considered in the modeling, and the ability of β-galactosidase expression for different HO-PBDEs may differ. Thus, it is necessary to develop multi parameter QSAR models for -logREC 20 prediction. Development, validation, and AD of the QSAR. PLS analysis with -logREC 20 as the dependent variable and the molecular structural parameters as predictor variables resulted in the following optimal QSAR model: -logREC 20 = 5.73 × 10 + 8.01 × 10 -1 n Br + 9.62 × 10 -1 logK ow -4.95 × 10 I A + 2.84 E LUMO -1.66 ω + 3.26 × 10 -2 µ 2 , where n = 14, A = 3, R 2 = 0.913, Q 2 CUM = 0.873, RMSE (training set) = 0.418, n EXT = 4, Q 2 EXT = 0.500, and RMSE (validation set) = 0.731 (p < 0.0001).
Q 2 CUM of the QSAR was as high as 0.873, implying good robustness of the model. The differences between R 2 and Q 2 CUM did not exceed 0.3, indicating no over fitting in the model (Golbraikh and Tropsha 2002). Figure 3 shows that the predicted -logREC 20 values are consistent with the observed values for both the validation and training sets. The model has good predictive abilities, as indicated by Q 2 EXT = 0.500 and RMSE = 0.731. The AD is shown by the Williams plot ( Figure 4); h i values of all the compounds in the training and validation sets were lower than the warning value (h* = 1.500). Thus, none of the compounds are particularly influential in the model space, and the training set has great representativeness. The standardized residuals of all the compounds in the training and validation sets are < 3, so there are no out liers in the developed QSAR model. Considering the mechanism of action, we can infer that the developed QSAR model can be used to predict thyroid hormone activity on TRβ of other HO-PBDEs similar to those used in the present study. To discriminate between active and inactive compounds (e.g., PBDEs and HO-PBDEs), discriminant models should be developed in further studies.
Mechanistic implication of the developed model. The developed PLS model extracted on three PLS components loaded primarily on six predictor variables. Values of the variable importance in the projection (VIP) and PLS weights (W*) are listed in Table 2. From the W* values, one can see how the predictor variables and the response variable combine in the projections (PLS components) and how they relate to each other.
The first PLS component (W*c[1]) is loaded primarily on three descriptors, n Br , log-K ow , and I A (Table 2). These three descriptors remarkably govern -logREC 20 values, as indicated by their large VIP values among the predictor variables. The coefficients in the QSAR model indicate that -logREC 20 values of HO-PBDEs increase with n Br and logK ow values but decrease with increasing I A values. The observation is reasonable because log-K ow correlates with n Br positively (r = 0.999, Green dashed lines indicate hydrogen bonds between HO-PBDEs and amino acid residues, gray is carbon, red is oxygen, blue is nitrogen, purple is iodine, and white is hydrogen.   TRβ. The deviation from planarity of an aromatic ring is a structural measurement of aromaticity: The higher the planarity, the higher the π electron delocalization and the higher the aromaticity (Alonso et al. 2008). Even for a single ring of polychlorinated biphenyl, aromaticity varied with the number of chlorine atoms (Alonso et al. 2008). Thus, the involvement of I A in the QSAR model may imply the effects of planarity (or non planarity) of HO-PBDE molecules on thyroid hormone activity. According to the QSAR model, I A plays a negative effect on thyroid hormone activity (n = 18; r = 0.66; p < 0.005). E binding also correlates positively with I A values (n = 18; r = 0.54; p < 0.05). The second PLS component (W*c[2]) is loaded primarily on I A , E LUMO , ω, and µ 2 , and the third PLS component (W*c[3]) is loaded primarily on n Br , logK ow , E LUMO , and ω. E LUMO itself has a negative value and measures the ability of a molecule to accept electrons. The docking analysis showed that hydrogen bonds were the charac teristic inter actions between the hydroxyl oxygen of HO-PBDEs and the hydrogen of Arg282 and Ile276. Because HO-PBDE molecules with lower E LUMO values tend to accept electrons easily, accept protons with difficulty, and accordingly form the hydrogen bonds with difficulty, the developed PLS model shows that -logREC 20 increases with E LUMO values. Likewise, because ω meas ures the ability of a molecule to soak up electrons (Chattaraj 2006), in the developed PLS model -logREC 20 increases with decreasing ω values. Finally, µ measures the dipole-dipole and dipole-induced interactions between inter acting molecules (Hurth et al. 2004), so HO-PBDEs with large µ values may exhibit strong dipole-dipole interactions with the receptor (TRβ), leading to large -logREC 20 values.

Conclusion
We determined thyroid hormone activities of selected HO-PBDEs by the recombinant two-hybrid yeast assay. Docking analysis indicated that hydrogen bonding and electrostatic interactions are the key steps for HO-PBDEs to exert thyroid hormone activities. We developed a QSAR to characterize the inter actions and model the thyroid hormone activities. The HO-PBDEs with higher ability to accept electrons (as indicated by E LUMO and ω) tend to have weak hydrogen bonding with the receptor and lower thyroid hormone activities. The developed QSAR model has good robustness, predictive ability, and mechanism interpretability.