Homology modeling to screen for potential binding of contaminants to thyroid hormone receptor and transthyretin in glaucous gull (Larus hyperboreus) and herring gull (Larus argentatus)

Thyroid hormone disrupting chemicals (THDCs) are of major concern in ecotoxicology. With the increased number of emerging chemicals on the market there is a need to screen for potential THDCs in a cost-efficient way, and in silico modeling is an alternative to address this issue. In this study homology modeling and docking was used to screen a list of 626 compounds for potential thyroid hormone disrupting properties in two gull species. The tested compounds were known contaminants or emerging contaminants predicted to have the potential to reach the Arctic. Models of transthyretin (TTR) and thyroid hormone receptor α and β (TRα and TRβ) from the Arctic top predator glaucous gull (Larus hyperboreus) and temperate predator herring gull (Larus argentatus) were constructed and used to predict the binding affinity of the compounds to the thyroid hormone (TH) binding sites. The modeling predicted that 28, 4 and 330 of the contaminants would bind to TRα, TRβ and TTR respectively. These compounds were in general halogenated, aromatic and had polar functional groups, like that of THs. However, the predicted binders did not necessarily have all these properties, such as the perand polyfluoroalkyl substances that are not aromatic and still bind to the proteins.


Introduction
Certain environmental pollutants have chemical structures resembling thyroid hormones (THs). These pollutants can act as TH analogs and compete with the THs for the binding to the serum TH transport protein transthyretin (TTR) or thyroid hormone receptors (TRs). Hydroxyl-polybrominated diphenyl ethers (OH-PBDEs), hydroxylpolychlorinated biphenyls (OH-PCBs), per-and polyfluoroalkyl substances (PFASs), bisphenol A (BPA), tetrabromo-BPA (TBBPA) and halogenated phenols are among the compounds reported to bind to the TH binding site in human TTR (hTTR) and/or TR (hTR) (see Fig. 1 for illustrations of molecular structure for these groups of chemicals) [1][2][3][4]. These compounds may reduce serum TH concentrations by displacing THs from the transport proteins, hence increasing hepatic excretion. They may alternatively bind to TRs and directly disrupt normal TH signaling [5,6].
THs are essential for several physiological processes such as reproduction, metabolism and development. They are found circulating in the plasma of all groups of vertebrates and are under control of the hypothalamus-pituitary-thyroid (HPT) axis. TTR and albumin are the major serum transport proteins of THs in birds [7]. TTR is synthesized in the liver and the brain choroid plexus, and transports TH in the blood and cerebrospinal fluid. TRs belong to the large family of nuclear receptors and regulate gene expression in response to THs. There are several TR isoforms and the major subtypes are thyroid hormone receptor α and β (TRα and TRβ). The isoforms differ in ligand affinity and specificity and are differentially expressed in different organs [8].
The glaucous gull (Larus hyperboreus) is a top predator in the Arctic marine food web and has therefore been used as a bioindicator species for long-range transported contaminants [9]. The herring gull (Larus argentatus) is in the same genus as the glaucous gull. Herring gulls are predators as well and live in temperate regions closer to big cities and point sources of pollution. Studies of glaucous gull have revealed relationships between exposure to contaminants and altered circulatory levels of THs [10,11]. Herring gull from highly contaminated sites in the Great Lakes had depleted thyroid gland hormone stores, cases of hypothyroidism and enlarged thyroid glands with depressed hormone stores [12]. Ucan-Marin et al. [13] found that OH-PBDEs and to a smaller extent MeO-PBDEs and OH-PCBs bind to TTR and albumin in glaucous gull and herring gull and are TH competitors. Furthermore Ouyang et al. [14] tested herring gull egg extracts in an hTTR in vitro assay and found that they interfered with T4-TTR binding.
Most often substances are identified as endocrine disruptors based on in vitro and in vivo studies [15]. These methods are time consuming and costly, limiting the number of chemicals that can be tested for being potentially toxic. At the same time new chemicals are continuously produced and for the design of safer chemicals it is crucial that potential harmful compounds are identified early during development of new chemicals. Rapid and cost-effective predictions by in silico methods can be useful for predicting putative harmful effects and prioritizing chemical entities prior to in vitro and in vivo testing.
There are several in silico methods for toxicity prediction, and ligand based quantitative structure-activity relationship (QSAR) and docking are two traditional methods. QSAR models correlates physiochemical properties of chemicals with biological activity but are limited by being derived from relatively small ligand datasets or are focused on specific chemo-types of compounds. Docking based affinity predictions are often used in drug development and require structural knowledge about the targeted protein. If the three-dimensional (3D) structure of the target or a related homologues protein is known, docking based methods can be used to screen huge libraries of compounds for putative binding. This will provide information on the TH-disruptive potential of a range of pollutants. In the present study, legacy and emerging compounds were docked into homology models of the glaucous gull and herring gull TRs and TTR to predict their binding affinities and potential thyroid disruptive properties.

Materials and methods
Since the 3D structures of the gull TRα, TRβ and TTR are not known, the homology modeling approach was used to construct 3D models. Homology modeling as well as docking was performed using the internal coordinate mechanics (ICM) software version 3.7 (http:// www.molsoft.com).

Construction of homology models
The homology modeling approach consists of four steps: 1) template identification, 2) amino acid sequence alignment, 3) model construction 4) refinement and evaluation of the model. To identify X-ray structures to use as templates the Protein Data Bank (PDB) was searched.
The amino acid sequences of glaucous gull TRα and TRβ were not available, but the complete sequences from chicken (Gallus gallus) and fragments of the sequences from herring gull were accessible at UniProt Knowledgebase (UniProtKB, http://www.uniprot.org/). In chicken and herring gull respectively, the UniProt ids of TRα are P04625 and Q5D226, and of TRβ P68306 and Q5D225. The glaucous gull and herring gull TTR amino acid sequence is identical [13] and accessible at UniProt id B0FWC5. This sequence is almost complete except for a small truncation in the N-terminal. However, this part was available for chicken with id P27731. From these sequences, hybrid sequences were constructed and used for homology modeling. Since the sequences of TTR are identical, it was reasonable to assume that the sequences of TRα and TRβ should be very similar between these species. The selected templates for TRα and TRβ where human, with PDB id 2 h79 and 2j4a respectively. Three TTR models were constructed from the templates 1kgi and 1kgj from rat (Rattus norvegicus) and 1sn0 sea bream (Sparus aurata).
The ICM software was used to construct amino acid sequence alignments and the 3D-homology models. ICM used a rigid body homology modeling method where the target model was constructed by transferring the backbone conformation of the core regions from the template to the target. The non-conserved loop regions were constructed through PDB loop searching by matching the loop regions with respect to sequence similarity and steric interactions with the surroundings of the model. The side chains of identical amino acids were transferred directly from the template, while the side chains of the nonconserved amino acids were either modeled or added to the target without reference to the template, using the most probable rotamer of side chains [16]. The ICM refineModel macro was used to energy optimize the constructed models and find the most likely conformation.

Evaluation of the models
To assure the quality and prediction capability of the models, separate compound test sets for TRα, TRβ and TTR were built, including potent binders and poor binders or decoys. Decoys are theoretical compounds that may function as negative controls resembling the active ligands in physiochemical properties but are topologically dissimilar.
From the ChEMBL database, strong binders and poor binders of the hTRs were found. The X-ray structure of hTRα was bound to the agonist T3, and hTRβ was bound to an antagonist (3,5-dibromo-4-(3-isopropyl-phenoxy)benzoic acid). Therefore, binding affinity values of the half maximal effective concentration EC 50 < 5 nm and the half maximal inhibitory concentration, IC 50 < 2 nm, Ki < 2 nm, log IC 50 > 8.7 were used for identifying TRα and TRβ binders, respectively. The ten most structurally different binders of each receptor were chosen for the test by using the Cluster Set function in ICM. To find poor TR-binders affinity values IC 50 > 2000 nm, Ki > 2000 nm, log IC 50 < -5.7 were selected, and 90 of the structurally most divergent compounds were chosen for each of the test sets.
For the TTR test set, TTR ligands were selected from the PDB database and clustered to select 20 structurally different compounds. To find decoys, a list of SMILES (simplified molecular-input line-entry system) codes of the ligands was sent to the database DUD-E (http:// dude.docking.org/). For each ligand 50 decoys were constructed, with 82 structurally different decoys selected for the test set. All the homology models were evaluated on the ability to separate strong binders from decoys/poor-binders using Receiver Operator Characteristics (ROC) curves calculated by the ICM nosauc macro.

Construction of thyroid hormones and contaminant-set
The chemical structures of the contaminants were constructed in MarvinSketch. Formal charge of the compounds was set to correspond to pH 7.4. The THs have a hydroxyl group (-OH) and a carboxyl group (-COOH) that both can be deprotonated. For this reason, three forms of both hormones where docked into the models. A dataset of 626 contaminants was constructed based on Howard and Muir [17] and on Vorkamp and Riget [18] studies. The former focused upon potentially persistent and bioaccumulative organic chemicals not considered in Great Lakes, North American, and Arctic contaminant measurement programs, and the latter upon potential Arctic contaminants. The present dataset focused on compounds not under regulations, but also included known suspected TTR-binding compounds.

Docking and scoring
The icmPocketFinder macro and the ICM Receptor Setup function were used to detect and define possible binding pockets in the models. The docking produces a target-ligand complex that should resemble the "native" complex in the biological system. Scoring predicts the binding energy between the protein and the ligand. For predicting the binding pose a Monte Carlo global optimization procedure was used. The protein was included as set of rigid pre-calculated grid potential maps representing van der Waal, hydrogen bonding, electrostatic and hydrophobic ligand-receptor interaction terms. During the global optimization procedure, the low energy conformations were saved and ranked based on scoring. Each compound was docked in three parallels in each model, and the conformations with the best scoring value was selected. The score was calculated with the ICM docking energy function.

Model evaluation
The amino acid sequence of the binding site in both the TRs and the TTR is highly conserved through evolution [13]. In the TRs there was only one amino acid that differed between the chicken and gull sequences: an aspartate in chicken was changed to glutamate in gull TRα outside of the ligand binding pocket (BP). The TR models had the same number of helixes and inner ligand binding domain (LBD) covered by the helixes of the C-terminus as the templates. The TRα and the TRβ BP consists of 26 amino acids and 23 amino acids respectively. Only one amino acid differed in the LBD of the isoforms where TRα had a serine and TRβ had an asparagine.
The constructed gull TTR (gTTR) models consisted of a homotetramer with a dimer-dimer interface forming a central binding channel with two LBDs. For easier comparison between templates and targets, the numbering of hTTR was used. The BP consisted of thirteen, ten and nine amino acids in each subunit of the protein Models 1, 2 and 3 respectively (Table S2). The conformation of the binding site differed slightly between TTR models since the co-crystallized ligands are causing differences in side chain conformations between the templates and thereby between our models. These differences may account for conformational flexibility of the binding site. There were differences between the modeled BPs and the BP of hTTR. The amino acids E54 and M13 were in the BP of hTTR and Model 1 but not in the BP of Models 2 and 3. Models 1 and 3 included V16 and Models 1 and 2 included T118 that were not a part of the hTTR BP. The docking evaluation showed that the ROC area under the curve (AUC) was 0.74 for TRα, and 0.70 for TRβ. For the homology Models 1, 2 and 3 of gTTR the AUC was 0.96, 0.95 and 0.88 respectively. These values qualify that the built homology models were fair in predicting which ligands are binders and not.
The score is not a quantitative measure of ligand binding and it was therefore necessary to estimate a scoring threshold based on the score of confirmed binders. In gTTR Models 1, 2 and 3 the threshold was defined as the highest score (where higher scores equal lower affinities) of the TTR binders in the test set; −18, −17 and −16 respectively. In the TRβ model two of the binders displayed scoring values of 7.7 and −4.5 and hence were predicted not to bind. They were subsequently excluded when determining the thresholds. The mean scoring value of the binders was set as threshold for TRα and TRβ, −29.4 and −35.1 respectively.

Docking of the contaminants and the thyroid hormones
The contaminants with a score lower than the thresholds were considered as binders. Of the 626 contaminants, 28 (4.47%) and four (0.64%) bound to gull TRα and TRβ, respectively, while the corresponding numbers for the TTR models 1, 2 and 3 were 144 (23.00%), 146 (23.32%), and 230 (36.74%), respectively (Table S4, S5 and S6). In total 330 compounds, 52.72% of all the tested compounds, were predicted to bind to one or more of the TTR models. The ten best scoring compounds in each model is shown in Table 1. There was no clear relationship between protonation of the THs and the scoring values for the gTTR models, but the THs were in general bound the strongest in Model 1 (Table S3).

TRα model
The TRα model predicted that several PFASs could bind to the LBD, including two perfluoroalkyl carboxylic acids (PFCAs), four perfluoroalkane sulfonamides (FASAs), one perfluoroalkyl alcohol/ketone, two fluorotelomer acrylates and two fluorotelomer methacrylates (scores −30.03 to −38.28). Long-chained PFCAs were overall the strongest binders -although in general PFASs registered higher scores (and hence lower apparent affinities) than THs. There are few studies on PFASs and TH disruption, but one study by Ren et al. found that PFASs bound to TRs and affect TR-signaling in humans [19]. The present gTRα model suggests that TR binding is a possible mode of action (MOA) for PFASs. However, observed effects of PFASs on thyroid homeostasis may also be caused by PFAS activation of peroxisome proliferator-activated receptors which heterodimerize with retinoid X receptor causing reduced activity of TRs [20]. Binding of PFASs to TR could explain the observed change in TH-responsive gene expression in primary cultures of herring gull neuronal cells as reported in a study by Vongphachan et al. [21], although these changes was observed when exposed to the short-chained and not long-chained PFASs.
The PBDEs, PCBs and their metabolites are known as TH mimics. Modeling studies, binding assays and signaling transduction assay of hTRα show that OH-PBDEs and OH-PCBs could affect TH signaling through binding to TH BP [22,23]. However, there is evidence that the observed response in signal transduction in humans is caused by PCBs and OH-PCBs suppressing transcription through partial dissociation of the receptor complex from the T3-response element in DNA [23]. To the author's knowledge experimental avian studies of binding of these compounds to TRα have not been published. However, because the nuclear receptor is highly conserved through evolution binding of the compounds to the receptor is likely. This is shown in the gTRα model (Table S5), four PBDEs (scores −30.73 to −31.90), six OH-PBDEs (scores −29.81 to −32,98) and two OH-PCBs (score −31.50 to −31.49) were predicted to bind with scores similar the PFASs. Of the OH-PBDEs and PBDEs, mainly tetra-and penta-brominated compounds bound to the model. Of the OH-PCBs, penta-and hexa-brominated congeners were the best binders. The results were therefore in line with human studies [22,23].
It therefore seems that the optimal ligand structure should mimic THs. In contrast, PFASs that are not aromatic and have a very different structure from the THs were also predicted to be equally good binders. However, their structures are similar to fatty acids and acyl-CoA esters which also are TR ligands [24,25]. The TRα model predicted that the contaminants formed hydrogen bonds between the polar functional groups and R228 just like the -COOH of THs (Fig. 2). The contaminants also had hydrogen bonds with other amino acids (N179, A180, R262, R266, T275, S277, G278, L287 and G290), where S277 was the most important forming bonds to most of the compounds. In agreement with our results Ren et al. [19] found that the polar end of PFASs had hydrogen bonds with R228, and that some PFASs also formed hydrogen bonds with other residues such as R262 and S277. Longer carbon chains enhance the hydrophobicity of PFASs. Long-chained PFCAs therefore interact and bind more strongly than their shorter-chained counterparts to the highly hydrophobic LBD of TRα.
PFASs are used in many different classes of man-made substances  1 Compound was in top ten for different targets. 2 Compound was in top ten in several of the TTR models.
such as water-, soil-and stain-resistant coatings for different materials, aviation hydraulic fluids, fire-fighting foams, adhesives, waxes, polishes, paints and many other products. In industry PFASs are used as surfactants, emulsifiers, additives, coatings and wetting agents. The strong and stable carbon-flourine bonds make the compounds very resistant to environmental degradation. There is large variation among the different groups of PFASs in chain-length, degree and pattern of fluorination, molecular weight, presence of polar functional groups and so on that makes it hard to determine their general physio-chemical properties and environmental fate [30].

TRβ model
In the TRβ model only PFASs were predicted to bind, three longchained PFCAs and a single FASA having score beneath the threshold ( Table 1). None of the other compounds were predicted to bind, although several of the PBDEs and OH-PBDEs that were docked had scores just above the threshold. This indicates that there was a difference in affinity to compounds between the different gTR isoforms. Studies on binding of contaminants to the TRβ have reported that an array of halogenated or aromatic compounds could bind to TRβ in humans and fish, and in vitro studies using reporter gene assays or GH3 cell proliferation assay revealed that several compounds can affect hTRβ-dependent gene expression [31][32][33][34][35]. Many of these studies used docking studies to support the claim that the effects are caused by direct binding to the TRβ [36,37]. This diversity in binding is not reflected in the present data from docking of contaminants to the gTRβ model as less than 1% of the compounds bound to the receptor. However, Kollitz et al. [32] showed that the affinity of zebrafish TRβ differed from hTRβ. Species-differences could possibly explain the lower affinity shown by the gTRβ model compared to hTRβ.
Docking in the gTRβ model showed that the amino acids R282 and R320 formed hydrogen bonds with the -COOH of THs and PFCAs (Fig. 3). The -OH of THs formed hydrogen bonds to H435. R320, A279 and R316 formed a hydrogen bond with the amine group in the FASA diammonium N-ethylheptadecafluoro-N-[2-(phosphonatooxy)ethyl]octanesulfonamidate. Alanine and arginine form a small region within the BP with strong hydrogen bond donating potential. The rest of the LBD consists of hydrophobic amino acids forming a highly hydrophobic pocket which is important for binding of hydrophobic compounds [36]. Other studies have identified R282, R320, N331 G332, T329 and H435 as important for binding to TRβ, which is in agreement with the gTRβ models identifying R282, R320 and H435 as important for forming hydrogen bonds [36,38].

TTR models
The legacy pollutants PBDEs, PCBs and their metabolites have previously experimentally been proven to bind to TTR both in humans as well as glaucous gull and herring gull [5,13]. These compounds were therefore included in the docking to validate the prediction of the TTR models. The predicted binding affinities from best to poorest were OH-PBDEs (fourteen congeners, top score −25.72) > PBDEs (twelve congeners, top score −26.25) > OH-PCBs (twelve congeners, top score −25.63) > methoxyl (MeO)-PBDEs (five congeners, top score −25.16) > methyl sulfone (MESO 2 )-PCBs (22 congeners, top score −21.61) > PCBs (five congeners, top score −21.05). Of these compounds 4′-OH-BDE-49 had the lowest score on average in the models and the strongest TTR-binder. Hydroxylated metabolites with OH-in the para position and around five halogen atoms were the strongest predicted binders. The MeSO 2 -PCBs were predicted to be TTR-binders primarily in Model 3, only six of these were binders in Model 1 and none of them bound in Model 2. The highest scoring OH-PBDEs, PBDEs and OH-PCBs had scores equal to T4 and T3 in Models 2 and 3, but much poorer in Model 1. These results are in accordance with Ucan-Marin et al. [5,13] results from in vitro competitive binding of PBDEs, PCBs and their metabolites to recombinant gTTR, and validated the model predictions.
For the PFAS groups with compounds of different length there was a trend of decreasing scores (indicating higher affinity) with increasing chain length up to a certain length. The longer-chained PFASs are more hydrophobic and interact with the hydrophobic part of the BPs, however with longer chains the compounds start to be too large for the BP. This is consistent with studies on hTTR [4,39,40]. Binding of the longchained PFCAs to TTR are a concern since they are some of the PFASs found in high concentrations in wildlife across the globe, and they biomagnify through the food chains [41][42][43]. Long-chained PFCAs such as PFTriA and PFTeA have been found to positively correlated with the  TH levels in glaucous gull from Svalbard [44]. Binding of these compounds to TTR as well as to TRs, potentially causing increased excretion of T4 and T3 can be part of the explanation for the detected positive correlation. Because of the considerable evidence showing that PFSAs and long-chained PFCAs can bioaccumulate and has toxic effects, the production and use of these compounds have been phased out or regulated [30,45]. However, new PFASs have been taken into use as substitutes. Their persistency in the environment, potential for longrange transport, potential to bioaccumulate and toxicity is partly unknown. Some studies indicate that these substitutes are broken down to PFSAs and PFCAs in the environment [30]. The potential risk and toxic effects on wildlife is uncertain. The TTR and TR models identify several of these new PFASs as potential THDCs that can disrupt TH homeostasis by interfering with their circulatory transport of THs and TH-signaling indicating that they can be just as toxic to wildlife.
There are large chemical variations among the PFASs, especially in the molecular size and partial charge characteristics such as their relative polar and hydrophobic surface area [4]. In contrast to the present study, human studies have shown that fluorotelomer alcohols and/or Nsubstituted perfluoroalkyl sulfonamides did not bind to hTTR. For hTTR only the PFASs with acidic functional groups and not those with nonacidic/neutral functional groups could bind [4,39,40]. For gTTR both group of PFASs were predicted to bind. These species differences indicate the need for studies on specific keystone wildlife species for assessing the potential effects of THDCs on biodiversity and ecosystem functioning.
Binding between PFASs and gTTR showed that the functional group of many neutral PFASs and some longer-chained acidic PFASs were oriented to the inner part binding to S117 (Fig. 4) and not towards the outer part and K15. The outer part of the binding pocket is wider than the inner part, giving more space for the long hydrocarbon chains. Previous studies on the binding of PFASs to hTTR found that acidic functional groups formed strong hydrogen bonds with K15 [39,40,46]. The hydrocarbon tail then faces the interior of the pocket and the compounds with eight carbons fit the best within this hydrophobic region of hTTR, therefore having the highest affinity. The PFASs with longer chains were too large and required that the fluorinated carbon tail bend for the compounds to fit inside the pocket. The neutral PFASs will locate their functional group towards the inner part of the hTTR pocket forming hydrogen bonds with S117, however the binding energy were lower [39,40,46]. Yang et al. [46] suggested that this is because PFASs are aliphatic and not aromatic, and hence are incapable of forming cation-π interactions with S117. In the TTR models the PFASs interacting with S117 didn't have a much lower predicted binding energy than then PFASs facing K15 even do they are not aromatic.
In addition to the contaminants discussed above, 168 other contaminants were predicted to bind to gull TTR (Table S6), mainly in Model 3. The best binders of these compounds in each TTR model are listed in Table 1, and 4-Chloro-3-nitrobenzoic acid had the lowest score (indicating higher affinity), −31.36 in Model 2. However, most of the scores were just below the threshold, which indicate that they potentially are only very weak binders. Structural similarities of these compounds with THs are apparent: almost all are aromatic, with the majority likewise containing one or two benzene rings. Most are halogenated (Cl, Br, I and F), but 61 of the compounds do not contain halogen atoms such as the phosphorus flame retardant triphenyl phosphate (TPhP) and the phathalate butyl benzyl phthalate (BBP) ( Table S6). The compounds consist of functional groups such as -OH, amino groups, -COOH and nitro groups, and some have ether, ester or thiol groups, which can be polar and form hydrogen bonds to amino acids within the BP. Other in silico and in vitro studies that have showed in agreement with the present study that many structurally different compounds have affinity for TTR [35,47,48]. S117 formed hydrogen bonds with the amino and carboxylic acid units of the THs, with the TH hydroxyl further interacting with K15. The opposite orientation was also observed with the -COOH facing towards the outer part of the binding pocket interacting with K15. In the mid-region of the binding pocket, the pocket is highly hydrophobic and interacts with the halogen atoms in the THs. This is consistent with previous studies on the binding of THs to TTR [2,49,50]. The TTR models showed that the amino acids S117 and K15 were very important for binding of the ligands. Both binding modes are observed for the structurally diverse group of compounds predicted to bind to the gull TTR. This is in accordance with the study of Xhang et al. [47] who also found that THDCs interact like THs with hydrogen bonds to S117 in hTTR.
Many of the contaminants listed within Table S6 are in use and considered as emerging, their environmental distribution and toxicity is largely unknown. The compounds have been detected in the temperate latitudes of the Northern Hemisphere where the majority of global industry is located. It is also predicted for some of the emerging compounds that they can be persistent and be subjected to long-range transport to polar regions [18]. Thus, these potentially emerging THDCs may pose risk to wildlife in temperate and Arctic ecosystems such as different gull species. Modeling studies, like the present study, can be applied to identify possible MOAs of emerging compounds and indicate the potential risk of compounds being TH disrupting.
Of the top ten strongest binders identified through the TTR models, nine are compounds that are not PFASs, PCBs, PBDEs or their metabolites. Amongst these is bis(2,4-dichlorophenyl)peroxyanhydride associated with plastic packing [51]. 4,4′-dibromobenzil is used as colorant and heat stabilizer [52]. 3-(2-chloro-4-trifluoromethylphenoxy) benzoic acid is a pesticide intermediate which also appeared amongst the top ten compounds displaying greatest predicted binding affinity within the gTRα model. Benzoic acid, 2,3,4,5-tetrachloro-6-cyano-, methyl ester is used in coating products, inks, toners and polymers and is likely to be released to environment from outdoor use of materials like treated wood products, treated textile, vehicles or other products based on metal, wood, paper and plastic [53]. Phenyl 1-hydroxy-4nitro-2-naphthoate is likely used as dye [54]. BAPP is used as a heatresistant solidifying agent in polyester-type materials. 4-Chloro-3-nitrobenzoic acid is an intermediate in production of dyes. Tetrachloro-obenzoquinone is used as an oxidant in biochemical and chemical synthesis, and N-(4-bromo-2,6-dichloro-3-methylphenyl)acetamide is likely a herbicide intermediate [17].
Most investigations of binding to TTR focus on humans. The present study shows that compounds with high binding affinity to hTTR in some cases had high scores, and thus low binding affinity in the gTTR models. Examples of this is TBBPA and chlorophenols which were shown to bind relatively strongly to hTTR [55], but not to the gTTR models. Ucan-Marin et al. [5] showed species-specific differences between hTTR and gTTR for competitive binding of PBDEs, PCBs and their metabolites. Computer-based models can be helpful in investigating when the species-specific differences are high and for extrapolation of results between species.

Conclusion
Overall the results indicated that a diverse group of compounds can bind to gull TTR, TRα and TRβ, and that in silico modeling is a good tool for rapidly and cost efficiently identifying these compounds in large databases of chemicals. The models identify binding of compounds to TTR and TRs as a possible MOAs that could lead to TH disruptions, whilst detecting structural properties which are of importance when interacting with these proteins. This information is valuable for planning further studies on THDCs.

Author statement
All authors have contributed to the manuscript.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.