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

1Polybrominated 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. Most recently, hydroxylated PBDEs (HO-PBDEs) have caused increasing concern because of reports of their natural production and metabolism. HO-PBDEs have been detected in the blood of fish, birds, and mammalian species and even in the abiotic environment. Recent studies suggest that some of the toxic effects of PBDEs might be due to their HO metabolites.

2In 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). For instance, levels of serum thyroxine (T4) were significantly decreased when rats were exposed to PBDEs. The effects of PBDEs on T4 levels may require metabolic activation because HO-PBDEs, but not the PBDE congeners themselves, behave as ligands for human TTR in vitro. also reported that 4-OH-BDE-90 and 3-OH-BDE-47 markedly inhibited the binding of triiodothyronine (T3) to TRα and acted as thyroid hormone-like agents. Recently, reported that 4-HO-BDE-90 significantly inhibited TRα- and TRβ-mediated transcriptional activity induced by T3. Consequently, HO-PBDEs have attracted great attention. However, there is lack of systematic investigation into the mechanisms by which HO-PBDEs interfere with hormonal actions.

3There 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β. The initial step for chemical mode of action is binding to an intracellular receptor. 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). Furthermore, molecular docking and virtual screening have become an integral part of many modern structure-based computational simulations of chemicals. 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.

4In 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 physicochemical 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 interactions 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.

5We selected 18 HO-PBDEs for the present study [see Supplemental Material, Figure S1 (doi:10.1289/ehp.0901457)]; most have been detected in environmental/biological samples. The HO-PBDEs congeners (> 97% purity) were purchased from AccuStandard (New Haven, CT, USA). 3,3′,5-Triiodothyronine (T3; 95% purity), dimethyl sulfoxide (DMSO; GC, 99.5% purity), o-nitrophenyl β-d-galactopyranoside (o-NPG; ≥ 98% purity), sodium dodecyl sulfate (99% purity), leucine (99% purity), tryptophan (99% purity), yeast-based nitrogen (99% purity), and β-mercaptoethanol (99% purity) were purchased from Sigma Chemical Company (St. Louis, MO, USA). Stock solutions of HO-PBDEs were prepared in DMSO.

6The recombinant two-hybrid yeast system employed a yeast cell transformed with the human TRβ plasmid, coactive plasmid, and the reporter gene expressing β-galactosidase. We examined the specificity of the yeast two-hybrid assay for TRβ ligand using DMSO (control), T3, and other steroid hormones. T3 induced β-galactosidase activity, whereas 17β-estradiol, dihydrotestosterone, and progesterone did not. Thus, the recombinant two-hybrid yeast assay was highly specific for TRβ ligand without cross-talk to other receptor agonists.

7We performed the recombinant two-hybrid yeast assay as described previously by. 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 (OD600) of 0.75. All the samples were determined at least in triplicate. Each triplicate included a positive control (T3) 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 × 103 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 OD600 (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 Na2HPO4·12H2O; 6.22 g/L NaH2PO4·2H2O; 0.75 g/L KCl; 0.25 g/L MgSO4·7H2O) and 20 μL chloroform, the cultures were carefully mixed and preincubated 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 measured at 420 nm. The β-galactosidase activity (U) was calculated according to the following equation:

8where 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), OD600 is the cell density measured at 600 nm, and OD420 and OD′420 are the cell density of the enzymic reaction supernatant and the blank, respectively, measured at 420 nm.

9The dose–response curves for U of the tested compounds were fitted by iterative four-parameter curve fit method using SigmaPlot, version 10.0 (Systat Software Inc., Chicago, IL, USA). The concentration inducing 20% of the maximum effect (REC20) 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).

10We 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. 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β (Ebinding) was then computed. The docking analysis provided insights into the interactions between the ligands and the receptor, which facilitated the selection of appropriate molecular parameters to characterize the interactions in the QSAR studies.

11We hypothesized 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 (logKow), average molecular polarizability (α), molecular volume (V), dipole moment (μ), energy of the highest occupied molecular orbital (EHOMO), energy of the lowest unoccupied molecular orbital (ELUMO), formal charge on hydroxyl hydrogen atoms (qOH), formal charge on hydroxyl oxygen atoms (qOH), formal charge on ether oxygen atoms (qO), electrophilicity index (ω), harmonic oscillator model of aromaticity index (IA), and the number of bromine atoms (nBr) [see also Supplemental Material, Table S1 (doi:10.1289/ehp.0901457)]. We purposely selected logKow to describe the partition process. The parameters V, nBr, α, and μ also partly describe partition because many of these parameters correlate with logKow. The parameters EHOMO, ELUMO, qOH, qOH, qO, ω, and IA were purposely selected to describe the intermolecular electrostatic interactions between the ligands and TRβ. The quantum chemical parameters EHOMO, ELUMO, qOH, qOH, and qO proved successful in many QSAR studies for characterizing intermolecular electrostatic interactions. ω measures the ability of a compound to soak up electrons. The relative binding affinity of some estrogen derivatives correlated strongly with ω. The aromaticity of compounds (IA) may influence their noncovalent interactions with the receptor, and IA has been used to characterize halogenated biphenyls.

12We computed logKow values using the EPI Suite, version 4.0. V (defined as the volume inside a contour of 0.001 electrons/bohr3 density) and the quantum chemical parameters were computed by the Gaussian 03 programs. Initial geometries of HO-PBDEs were preoptimized by semiempirical 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. IA was calculated by DRAGON software.

13The 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 developing the model because PLS can analyze data with strongly collinear, noisy, and numerous predictor variables. 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). Cross-validation 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 (Q2CUM) for the final model. The PLS analysis was performed repeatedly to eliminate redundant molecular structural parameters, as done in our previous studies. Model predictability was evaluated by external validation, which was characterized by the determination coefficient (R2), root mean square error (RMSE), and external explained variance (Q2EXT), which are defined as follows:

14where yifit is the fitted − logREC20 value of the ith compound; ȳ is the average response value in the training set; yi 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 nEXT is the number of compounds in the validation set.

15We assessed the AD of the developed QSAR model using the Williams plot, that is, the plot of standardized residuals (σ) versus leverage (hat diagonal) values (hi). We calculated σ as follows:

16where yi and ŷi are the observed and predicted values for the ith compound, respectively, and n is the number of compounds in the training set.

17The hi value of a chemical in the original variable space and the warning leverage value (h*) are defined as follows:

18where xi 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.

19Based on a plot of U versus log T3 concentrations, the maximal induction of T3 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 T3 was 1.4 × 10−7 M, which was similar to that reported by. 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 REC20 values for the 18 HO-PBDEs.

20With 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 indicated that HO metabolites of PBDEs could compete with T4 for binding to TTR and exert thyroid hormone activity. reported that 4-OH-BDE-90 and 3-OH-BDE-47 markedly inhibited the binding of T3 to TRα and acted as thyroid hormone-like agents. reported that 4-HO-BDE-90 significantly inhibited TRα- and TRβ-mediated transcriptional activity induced by T3. Thus, with respect to the observed thyroid hormone activity of HO-PBDEs, our results are consistent with previous findings.

21Figure 1 shows the docking view of T3 and representative HO-PBDEs and PBDEs (4-OH-BDE-42, 4′-OH-BDE-17, and BDE-116) 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 interaction. 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 π–π interactions between the phenyl of HO-PBDEs and Phe272, Phe442, and Phe455.

22The ligand–receptor binding energy (Ebinding) of the 18 HO-PBDEs is listed in Table 1. As shown in Figure 2, we obtained a simple linear free energy relationship between –logREC20 and Ebinding, further proving that binding to TRβ is the key step for the HO-PBDEs to exert their thyroid hormone activity. However Ebinding itself was not a good predictor for − logREC20, as indicated by the big prediction residuals for some HO-PBDEs (Figure 2). This was expected because Ebinding 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 multiparameter QSAR models for − logREC20 prediction.

23PLS analysis with − logREC20 as the dependent variable and the molecular structural parameters as predictor variables resulted in the following optimal QSAR model:

24where n = 14, A = 3, R2 = 0.913, Q2CUM = 0.873, RMSE (training set) = 0.418, nEXT = 4, Q2EXT = 0.500, and RMSE (validation set) = 0.731 (p < 0.0001).

25Q2CUM of the QSAR was as high as 0.873, implying good robustness of the model. The differences between R2 and Q2CUM did not exceed 0.3, indicating no overfitting in the model. Figure 3 shows that the predicted − logREC20 values are consistent with the observed values for both the validation and training sets. The model has good predictive abilities, as indicated by Q2EXT = 0.500 and RMSE = 0.731.

26The AD is shown by the Williams plot (Figure 4); hi 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 outliers 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.

27The 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.

28The first PLS component (W*c[1]) is loaded primarily on three descriptors, nBr, logKow, and IA (Table 2). These three descriptors remarkably govern − logREC20 values, as indicated by their large VIP values among the predictor variables. The coefficients in the QSAR model indicate that − logREC20 values of HO-PBDEs increase with nBr and logKow values but decrease with increasing IA values. The observation is reasonable because logKow correlates with nBr positively (r = 0.999, p < 0.001) and because HO-PBDEs with large logKow and nBr values tend to partition into the biophase (yeast cells). The harmonic oscillator model of aromaticity index (IA) may characterize the noncovalent interactions with 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. Even for a single ring of polychlorinated biphenyl, aromaticity varied with the number of chlorine atoms. Thus, the involvement of IA in the QSAR model may imply the effects of planarity (or nonplanarity) of HO-PBDE molecules on thyroid hormone activity. According to the QSAR model, IA plays a negative effect on thyroid hormone activity (n = 18; r = 0.66; p < 0.005). Ebinding also correlates positively with IA values (n = 18; r = 0.54; p < 0.05).

29The second PLS component (W*c[2]) is loaded primarily on IA, ELUMO, ω, and μ2, and the third PLS component (W*c[3]) is loaded primarily on nBr, logKow, ELUMO, and ω. ELUMO itself has a negative value and measures the ability of a molecule to accept electrons. The docking analysis showed that hydrogen bonds were the characteristic interactions between the hydroxyl oxygen of HO-PBDEs and the hydrogen of Arg282 and Ile276. Because HO-PBDE molecules with lower ELUMO values tend to accept electrons easily, accept protons with difficulty, and accordingly form the hydrogen bonds with difficulty, the developed PLS model shows that − logREC20 increases with ELUMO values. Likewise, because ω measures the ability of a molecule to soak up electrons, in the developed PLS model − logREC20 increases with decreasing ω values. Finally, μ measures the dipole–dipole and dipole-induced interactions between interacting molecules, so HO-PBDEs with large μ values may exhibit strong dipole–dipole interactions with the receptor (TRβ), leading to large − logREC20 values.

30We 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 interactions and model the thyroid hormone activities. The HO-PBDEs with higher ability to accept electrons (as indicated by ELUMO 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.