Modelling inhibition of avian aromatase by azole pesticides

^†Presented at the 8th International Symposium on Computational Methods in Toxicology and Pharmacology Integrating Internet Resources, CMTPI-2015, 21–25 June 2015, Chios, Greece.

Abstract

The potential effects of pesticides and their metabolites on the endocrine system are of major concern to wildlife and human health. In this context, the azole pesticides have earned special attention due to their cytochrome P450 aromatase inhibition potential. Cytochrome P450 aromatase (CYP19) catalyses the conversion of androstenedione and testosterone into oestrone and oestradiol, respectively. Thus, aromatase modulates the oestrogenic balance essential not only for females, but also for male physiology, including gonadal function. Its inhibition affects reproductive organs, fertility and sexual behaviour in humans and wildlife species. Several studies have shown that azole pesticides are able to inhibit human and fish aromatases but the information on birds is lacking. Consequently, it appeared to be of interest to estimate the aromatase inhibition of azoles in three different avian species, namely Gallus gallus, Coturnix coturnix japonica and Taeniopygia guttata. In the absence of the crystal structure of the aromatase enzyme in these bird species, homology models for the individual avian species were constructed using the crystal structure of human aromatase (hAr) (pdb: 3EQM) that showed high sequence similarity for G. gallus (82.0%), T. guttata (81.9%) and C. japonica (81.2%). A homology model with Oncorhynchus mykiss (81.9%) was also designed for comparison purpose. The homology-modelled aromatase for each avian and fish species and crystal structure of human aromatase were selected for docking 46 structurally diverse azoles and related compounds. We showed that the docking behaviour of the chemicals on the different aromatases was broadly the same. We also demonstrated that there was an acceptable level of correlation between the binding score values and the available aromatase inhibition data. This means that the homology models derived on bird and fish species can be used to approximate the potential inhibitory effects of azoles on their aromatase.

Keywords:

• Homology modelling
• aromatase
• azole pesticides
• human
• bird
• fish
• endocrine disruptors

1. Introduction

Since the 1500s, more than 150 bird species have been lost worldwide. To date, one in eight bird species is susceptible to disappearing, and about 200 species are particularly endangered [1]. Human activities are directly or indirectly responsible for this situation. Forests currently cover about four billion hectares, representing nearly 31% of the Earth’s land surface. While over a period of 5000 years the average net loss of forest was 360,000 hectares per year, the annual destruction rate has now reached about 5.2 million hectares [2]. Numerous bird species can only live in this ecosystem. The situation is particularly critical in the tropical rain forests, in which a square kilometre of forest can contain several hundred species of birds, some of them being highly prone to extinction due to their specialized life history traits [3, 4]. Deforestation impacts the world’s climate. These changes affect the distribution, abundance and behaviour of numerous birds [5,6]. Modern agriculture, which can result from forest destruction, is also one of the greatest threats to bird species due to changes in their habitats and because of the massive use of pesticides that directly affect birds and that can also reduce the availability of their prey. Birds are exposed to pesticides in different ways. Eggs, chicks and adults can be directly contaminated following spray treatments. Indirect contamination can also occur via food absorption, during the preening process and so on. Exposure to pesticides can result in acute poisoning leading to death in a short time-period [7,8]. Pimentel [9] conservatively estimated that if only 10% of the bird population in the US was killed by pesticide treatments, this could lead to the destruction of 72 million birds per year. However, exposure to pesticides does not necessarily lead to acute poisoning. Depending on the chemical nature of the pesticides, their dose and the time of exposure, sublethal and chronic effects can be observed. Among these, those interfering with the endocrine system are of first concern. These chemicals, commonly referred to as endocrine-disrupting chemicals (EDCs), pose a particularly severe threat to bird health. EDCs are believed to exert their effects by (1) mimicking normal hormones such as oestrogens and androgens; (2) antagonizing hormones; (3) altering the pattern of synthesis and metabolism of hormones; (4) modifying hormone receptor levels; or (5) interfering in other signalling systems, which are indirectly in relation with the endocrine system, such as the immune and nervous systems [10,11].

Aromatase is a cytochrome P450 steroidogenic enzyme (CYP19) that catalyses the removal of the 19-methyl group and aromatization of the A-ring of androgens for the synthesis of oestrogens [12]. It is essential for maintaining a physiological balance between androgens and oestrogens in vertebrates and hence the enzyme plays a key role in the physiology of their reproduction [13]. In birds, aromatase is expressed at high levels in different brain regions [14,15]. It intervenes in sexual differentiation in birds, and it has been shown that the use of aromatase inhibitors could induce a permanent female-to-male sex reversal [16–18]. Aromatase has also shown its importance in mediating and activating many male reproductive processes, including courtship display, aggressiveness and copulation behaviour [19–23].

Although steroidal (e.g. exemestane) and non-steroidal (e.g. anastrozole, letrozole) aromatase inhibitors are important in the treatment of oestrogen-dependent postmenopausal breast cancers [24], they are also potential EDCs for human and wildlife. Thus, the non-steroidal aromatase inhibitors are mostly azole-type compounds; this chemical family is also widely used as fungicides in agriculture.

A quarter of the world’s production of fungicides falls to the group of azole compounds developed in the early 1970s. In 2003, one of the leading producers of fungicides in the world had a total revenue of €1.168 billion from this category of pesticides. The second top-selling product in 2003, with a revenue of €315 million, contained tebuconazole [25]. In Switzerland, nearly 40 tonnes of azoles are sold annually, which is approximately 5% of the active ingredient in fungicides [26].

The endocrine disruption potential of azole fungicides, including effects on aromatase, has been the subject of investigations in vitro [27,28] and in vivo [29–33] in different fish species. On the other hand, a very limited number of studies have been made on birds [34–36], although recent studies have shown that wild birds are potentially exposed to azole fungicides [37,38].

In this context, it was of interest to estimate the aromatase inhibition potential of structurally diverse azole pesticides in avian species by using in silico approaches. Because only the crystal structure of human aromatase (hAr) is available, homology modelling was used. Homology modelling approximates 3D structures to gain insights into structure–function relationships and to predict the binding geometries and affinities of molecules [39–41]. It allows a better understanding of the mechanism of action of molecules. Homology models were designed from the aromatase sequences of Gallus gallus, Coturnix coturnix japonica and Taeniopygia guttata. In addition to having sequences available, the choice of these species was governed by the fact that azole drugs used in medicine as aromatase inhibitors have been tested on them to better understand the physiology and ethology of these birds [22,42–47]. This offered invaluable structure–activity relationship (SAR) information as well as the possibility of comparing aromatase inhibition in human versus birds systems. For the same reasons, a homology model was also developed for the fish Oncorhynchus mykiss. The ultimate goal of this paper was to gain insights into the mechanism of action of the azoles and their difference of sensitivity against aromatases, and to see whether their behaviour was different among the studied species in order to derive interspecies structure–activity–activity–relationships.

2. Materials and methods

2.1 Homology modelling

The sequences of the aromatase proteins for the Japanese quail (Coturnix coturnix japonica, Q8JH69), the chicken (Gallus gallus, P19098), the zebra finch (Taeniopygia guttata, Q92112) and the rainbow trout (Oncorhynchus mykiss, Q8JG18) were retrieved from the Universal Protein Resource database. The sequence alignment for each individual species with the human aromatase template (pdb: 3EQM) [48] was performed using the ClustalW server (http://www.ch.embnet.org/software/ClustalW.html). Homology models for the three bird species and fish were developed using Modeller 9v7 [49], and the models were ranked according to the Discrete Optimized Protein Energy (DOPE) scores.

2.2 Docking

Docking of the azole compounds having aromatase inhibitory activity on human aromatase (hAr) was performed on the human aromatase protein and on the homology modelled protein of C. japonica, G. gallus, T. guttata and O. mykiss. To accommodate bulky ligands, open state human aromatase (OhAr) was generated by molecular dynamics simulation, which was further used to model the open state of brain aromatase in O. mykiss (OomAr). The docking experiments [50,51] were performed using the standard precision (SP) module of GLIDE in Schrodinger 9.0 [52]. The active site of the protein for human aromatase for docking the ligands was defined by the amino acids R115, I133, F134, F221, W224, I305, A306, D309, T310, V313, V370, L372, V373, M374, S478 and H480. Due to metal coordination of azole nitrogen, the haeme iron was defined as a metal constraint in the generation of grid. Initially, results obtained with a laboratory test based on human recombinant CYP19 and dibenzylfluorescein as substrate for a series of compounds were used to evaluate the reliability of the docking model. The compounds were also docked in the three avian and fish species to find any similarity and/or dissimilarity in binding of these compounds as compared to human aromatase. For these correlation studies a comparison of the docking scores of each bird or fish species with human docking scores was performed. The validated models were further used to predict the aromatase affinity of structurally diverse azoles and related compounds on human, avian and fish species. The homology-modelled proteins for the three bird and fish species were each subjected to multiple minimization procedures using default settings in the Macromodel suite of Schrodinger 9.0 [53]. The protein preparation wizard and Ligprep were used for protein and ligand preparation.

2.3 Molecular dynamics simulation to model the open state of human aromatase

The crystal structure of 3EQM suggests that the binding pocket is too small in size, because our attempts to dock clotrimazole, prochloraz and ketoconazole resulted in no output. In order to predict an extended binding pocket for the human aromatase that can bind large molecules, we conducted a molecular dynamics simulation (MDS) with ketoconazole by manually placing the molecule at the binding site. The coordinates of the ketoconazole molecule were collected from the crystal structure of lanosterol 14-alpha-demethylase (CYP51) (pdb: 3LD6) after the alignment of both protein structures [54]. Multiple protein minimizations were carried out to remove any steric clash between ketoconazole and the protein structure of hAr. The structure parameterization for the haeme-cysteine complex of hAr protein reported earlier was used for parameterizing same in the present study [55]. The atomic charges and geometric parameters for ketoconazole were optimized using GAMESS [56,57] at HF/6–31G* level. The partial charges were calculated using the RESP [58] program of AMBER11 [59–61]. The MDS was performed using NAMD 2.7 [62] using Amber99SB force field with corrections for Leu, Ile, Asp and Asn to model the protein; TIP3P for water; and GAFF with AM1-BCC for the ligand. The rectangular water box with 14 Å dimensions was embedded around the protein ligand complex with system neutralization by the addition of Na⁺ cations and Cl⁻ anions. The van der Waals cut-off distance was selected at 9.0 Å and the electrostatic interactions were calculated using the particle-mesh Ewald algorithm. In order to eliminate improper contacts, energy minimization was performed for 10,000 steps and then equilibrated for 200 ps at constant temperature (298 K) and pressure (1 bar) conditions via the Langevin dynamics (the collision frequency was 1.0 ps⁻¹), with a coupling constant of 0.2 ps for both parameters. The production simulations were performed for 1.5 ns at 298 K and 1 bar using a timestep of 2 fs. The SHAKE algorithm [63] was applied to all bonds and analysis, and visualization of the trajectories was performed in VMD [64].

2.4 Predicting the open state of O. mykiss brain aromatase

The open state of brain aromatase in O. mykiss (OomAr) was modelled from the open state of hAr by homology modelling using protocols similar to those discussed above for the closed state.

3. Results and discussion

3.1 Homology modelling

The sequence alignment of the three bird species, namely C. japonica, G. gallus and T. guttata, showed high homologous sequence identity and similarity with the human template. This was also the same for O mykiss. In the absence of the crystal structure of aromatase enzyme for these species, homology models for the individual avian and fish species were constructed using the crystal structure of human aromatase (hAr) (pdb: 3EQM). High sequence similarities for G. gallus (82.0%), T. guttata (81.9%), C. japonica (81.2%) and O. mykiss (81.9%) were obtained considering the crystal structure of human aromatase where no structural information for the amino acid residues (1–44) was present (Figure 1). Considering only the common regions of birds with the human aromatase the sequence similarity further increased to 89.8% for both G. gallus and T. guttata and to 88.9% for C. japonica. The 3D structure of the human aromatase enzyme comprises 10 β-strands and 12 alpha helices with the characteristic cytochrome P450 fold. The 12 alpha helices in the human aromatase are defined by A (69–80), B (119–126), C (138–152), D (155–174), E (187–205), F (210–227), G (242–267), H (278–287), I (293–324), J (326–341), K (354–366) and L (440–455) (Figure 1). The β-sheets are defined by sheet2: β4:381–383 and β5:386–388; sheet3: β8:473–475 and β9:479–481; sheet4: β7:458–461 and β10:491–494. The binding region is a putative pocket formed by the helix I305, A306, D309 and T310 from the I-helix, F221 and W224 from the F-helix, I133 and F134 from the B-C loop, V370, L372 and V373 from the K-helix–β3 loop, M374 from β3, and L477 and S478 from the β8–β9 loop. The haeme porphyrin complex present at the active site of the protein forms a chelate complex with the C436 at the L-loop. It has been proposed that the residues A306, D309 and T310 are the key residues involved in C-hydroxylation and A ring aromatization of androstenedione. On the basis of sequence analysis no difference was found regarding the important residues of the I-loop when compared with all the bird species and O. mykiss. Analysis of the other residues at the catalytic cleft also showed no difference regarding the I133 and F134 from the B-C loop, F221 and W224 from the F-helix, V370 from the K-helix—β3 loop and M374 from the β3 loop. However, minor differences were observed regarding L372, which was identical in all bird species except in O. mykiss where it was substituted by F372. Another residue, V373 from the K-helix–β3 loop, was found to be substituted by I373 in three bird species and T373 in O. mykiss. However, the sidechain of V373 faces opposed from the binding site and may have a small contribution regarding selectivity between the avian and fish species (Figures 2(a) and 2(b)).

Figure 1. Sequence alignment of C. japonica, G. gallus, T. guttata, O. mykiss and human aromatases.

Figure 2. (a) The alignment of homology-modelled proteins of C. japonica, G. gallus, T. guttata and O. mykiss with human protein. (b) Comparison of the important residues at the catalytic site of the aromatase enzyme in the studied species. In the three bird species, V373 in humans is replaced by I373, while in O. mykiss it is replaced by T373. L372 in human is replaced by F372 in O. mykiss.

3.2 Docking

3.2.1 Study of fadrozole, vorozole and imazalil

Because it has been demonstrated that fadrozole, vorozole and imazalil presented an inhibitory activity on human, bird and fish aromatases [22,35,42–44,46,47,65–67], these three chemicals were considered for initial docking investigations. These compounds present one asymmetric centre (Figure 3). Consequently, the two enantiomers S and R were docked and the average of the docking scores was used for analysing the activity of the racemic. It was observed that the S-enantiomers showed better binding, as indicated by the scores. The interactions of S-enantiomers are described below.

Figure 3. Structures of 22 azoles and related compounds considered in this study.

The docked conformation of fadrozole (S) was found to bind the human aromatase enzyme with the azole nitrogen chelated with the iron atom, and the phenyl ring was found to have hydrophobic interactions with F134, V370, L372 and L477, while the nitrogen atom of the cyano group showed hydrogen bond interaction with M374 and R115 (Figure 4(a)). Similar hydrogen bond (H-bond) interactions with M374 and R115 were observed for fadrozole in C. japonica (Figure 5(a)). The hydrophobic interactions of fadrozole in C. japonica include interactions with I133, V370, L372 and L477, which closely resemble the hydrophobic interactions observed for human aromatase. In the case of G. gallus, fadrozole formed H-bond contact with M374 (Figure 5(d)) and hydrophobic contacts with F134 and V370. In T. guttata the hydrophobic interactions for the phenyl ring of the fadrozole (S) are V370 and L477, while the nitrogen atom of the cyano group made hydrogen bond interactions with M374 (Figure 5(g)). In O. mykiss the fadrozole (S) phenyl moiety made hydrophobic interactions with F134, V370 and L477 and the nitrogen atom formed a H-bond interaction with M374 (Figure 5(j)).

Figure 4. (a) Docking of fadrozole (S), (b) imazalil (S) and (c) vorozole (S) at the active site of human aromatase.

Figure 5. Docking of fadrozole (S), imazalil (S) and vorozole (S) at the C. japonica aromatase (a–c), at the G. gallus aromatase (d–f), at the T. guttata aromatase (g–i), and at the O. mykiss aromatase (j–l).

Imazalil was found to have numerous hydrophobic interactions at the human aromatase site besides having metal chelation with iron as observed in Figure 4(b). Hydrophobic interactions were observed between the imidazole moiety and I133. The terminal end of the propoxy chain interacted with F221 and W224. The 2,4-dichlorophenyl group was found to be the most important feature contributing to aromatase inhibitory activity of imazalil where the chlorine atom at position four of the phenyl ring had hydrophobic interactions with V370, L372, V373, M374 and L477, while the 2-chloro atom had interaction with F134. The interactions of imazalil with the aromatase of C. japonica were similar to human aromatase, except V373, which was replaced by I373 (Figure 5(b)). In the case of G. gallus aromatase the hydrophobic interactions of the chlorine atom at position four of the phenyl ring in imazalil were with I373, F134, M374 and L477, while the phenyl moiety made hydrophobic contacts with I133, F134, F221 and W224 (Figure 5(e)). The hydrophobic interactions of imazalil in T. guttata include hydrophobic conjunction of the 4-chloro group with F134, L372, I373, M374 and L477 and of the 2-chloro phenyl moiety with I133, F221, W224 and V370 (Figure 5(h)). In O. mykiss, the interaction of the imidazole moiety in imazalil was observed with I133. The terminal end of the propoxy side chain interacted with F221 and I313, while the 4-chloro atom of the phenyl ring showed hydrophobic interactions with I133, F134, W224, V370, M374 and L477 (Figure 5(k)).

In human aromatase the 4-chloro phenyl moiety of vorozole showed hydrophobic interactions with F221, W224, V370 and L477. In addition, the phenyl ring of the benztriazole moiety showed hydrophobic interactions with I133 and F134 and aromatic interactions with the porphyrin ring system (Figure 4(c)). The two unsubstituted nitrogen atoms of the triazole moiety were found to be engaged in H-bond interaction with M374. In C. japonica the 4-chloro phenyl moiety of vorozole made hydrophobic interactions with F221, W224 and V370. The phenyl ring of the benztriazole moiety showed similar hydrophobic and aromatic interactions as observed for human aromatase and aromatic interactions with porphyrin ring system (Figure 5(c)). A single hydrogen bond with M374 was also found for the nitrogen atom situated beside the methyl substituted nitrogen atom. In G. gallus the 4-chlorophenyl moiety showed hydrophobic interactions with F221, W224 and V370, while the phenyl ring in the benztriazole moiety had hydrophobic interactions with I133 and F134 and aromatic interactions with porphyrin ring system. Interestingly, the nitrogen atoms did not form a hydrogen bond with M374 but the 4-chloro phenyl ring made interactions with T310 at I-helix (Figure 5(f)). In case of T. guttata the interactions of vorozole include hydrophobic interactions of the 4-chlorophenyl moiety with F221, V370, L477 and H480. The hydrophobic interactions for the phenyl ring include interactions with I133, W224 and V370 (Figure 5(i)). In O. mykiss the 4-chlorophenyl moiety of vorozole interacted with F221 and I313 while the benztriazole moiety interacted with I133, W224, V370 and L477 (Figure 5(l)). The interaction of fadrozole, vorozole and imazalil with the active site residues of human, bird and fish aromatases described above have been summarized in Table 1. Analysis of this table reveals that they have common interactions pattern with rather insignificant differences in interaction between interspecies aromatases.

Table 1. Interactions of fadrozole, imazalil and vorozole with the active site residues of human, bird, and fish aromatases.

Species	Compound	Site
Human		R115	I133	F134	F221	W224	T310	I313	V370	L372	V373	M374	L477	H480	HEM600
	Fadrozole	HB	–	HY^a	–	–	–	–	HY	HY	–	HB^b	HY	–	Ar^c/M^d
	Imazalil	–	HY	HY	HY	HY	–	–	HY	HY	HY	HY	HY	–	M
	Vorozole	–	HY	HY	HY	HY	–	–	HY	–	–	HB	HY	–	M/Ar
C. japonica		R115	I133	F134	F221	W224	T310	I313	V370	L372	I373	M374	L477	H480	HEM600
	Fadrozole	HB	HY	–	–	–	–	–	HY	HY	–	HB	HY	–	M/Ar
	Imazalil	–	HY	HY	HY	HY	–	–	HY	HY	HY	HY	HY	–	M
	Vorozole	–	HY	HY	HY	HY	–	–	HY	–	–	HB	–	–	M/Ar
G. gallus		R115	I133	F134	F221	W224	T310	I313	V370	L372	I373	M374	L477	H480	HEM600
	Fadrozole	–	–	HY	–	–	–	–	HY	–	–	HB	–	–	M/Ar
	Imazalil	–	HY	HY	HY	HY	–	–	–	–	HY	HY	HY	–	M
	Vorozole	–	HY	HY	HY	HY	HY	–	HY	–	–	–	–	–	M/Ar
T. guttata		R115	I133	F134	F221	W224	T310	I313	V370	L372	I373	M374	L477	H480	HEM600
	Fadrozole	–	–	–	–	–	–	–	HY	–	–	HB	HY	–	M/Ar
	Imazalil	–	HY	HY	HY	HY	–	–	HY	HY	HY	HY	HY	–	M
	Vorozole	–	HY	–	HY	HY	–	–	HY	–	–	–	HY	HY	M/Ar
O. mykiss		R115	I133	F134	F221	W224	T310	I313	V370	F372	T373	M374	L477	H480	HEM600
	Fadrozole	–	–	HY	–	–	–	–	HY	–	–	HB	HY	–	M/Ar
	Imazalil	–	HY	HY	HY	HY	–	HY	HY	–	-	HY	HY	–	M
	Vorozole	–	HY	–	HY	HY	–	HY	HY	–	–	–	HY	–	M

^a hydrophobic interaction;

^b hydrogen bond interaction;

^c aromatic;

^d metal.

3.2.2 Docking of azoles on closed and open states of hAr

In order to consolidate our docking study, an attempt was made to correlate docking scores for 18 azoles (Table 2; Figure 3) in closed state hAR (ChAr) with their human recombinant CYP19 (aromatase) inhibitory activity values [25] (Table 2). A correlation coefficient of −0.92 was obtained between the docking scores of the 15 first racemic azoles in Table 2 and their corresponding pIC₅₀ values. These results clearly suggest that inhibition of human recombinant CYP19 activity of azoles may be estimated from their corresponding docking score values.

Table 2. Docking scores (Kcal/mol) vs. inhibition of human recombinant CYP19 (aromatase) for 18 azoles.

Compound^a	Docking scores (closed)	Docking scores (open)	IC50^b [μM]	pIC50
Flusilazole	−6.12	−4.82	0.055	1.26
Imazalil	−5.37	–	0.072	1.14
Myclobutanil	−5.75	–	0.47	0.33
Penconazole	−6.03	–	0.85	0.07
Epoxiconazole	−5.19	−3.83	1.44	−0.16
Propiconazole	−4.74	–	3.2	−0.51
Tebuconazole	−5.12	–	5.8	−0.76
Cyproconazole	−4.38	–	8.5	−0.93
Triadimenol	−4.32	–	12.6	−1.10
Triadimefon	−4.62	−4.5	17.5	−1.24
Hexaconazole	−4.67	−3.95	35	−1.54
Bifonazole	−7.17	−5.84	0.019	1.72
Miconazole	−6.67	–	0.064	1.19
Fadrozole	−7.78	−5.22	0.0076	2.12
Letrozole	−7.29	−5.17	0.015	1.82
Prochloraz	−6.23^p	−4.81	0.047	1.33
Clotrimazole	−5.8^p	−4.61	0.11	0.96
Ketoconazole	−5.23^p	−4.35	5.6	−0.75

^a See Figure 3 for the structures;

^b Concentration of azole required to reduce product fluorescence by a factor of two [25];

p = predicted docking score.

The same set of 15 azoles was also docked in the C. japonica, G. gallus, T. guttata and O. mykiss aromatase models. The obtained docking scores (Table 3) were then correlated to those obtained for human aromatase leading to pairwise correlation coefficients equal to 0.83, 0.83, 0.82 and 0.88 for C. japonica, G. gallus, T. guttata and O. mykiss, respectively. These results suggest that the 15 studied azoles behave similarly against the studied species.

Table 3. Docking scores (Kcal/mol) of azoles in the aromatases of the different studied species.

Compound	Human	C. japonica	G. gallus	T. guttata	O. mykiss
Flusilazole	−6.12	−5.08	−4.93	−6.24	−6.04
Imazalil	−5.37	−4.67	−4.42	−4.595	−4.78
Myclobutanil	−5.75	−4.52	−4.13	−4.5	−5.21
Penconazole	−6.03	−5.36	−5.6	−5.76	−5
Epoxiconazole	−5.19	−4.77	−4.38	−4.77	−5.48
Propiconazole	−4.74	−4.75	−4.11	−4.83	−4.51
Tebuconazole	−5.12	−5.51	−5.4	−5.12	−5.59
Cyproconazole	−4.38	−4.58	−4.37	−4.25	−4.63
Triadimenol	−4.32	−4.82	−4.65	−5.39	−4.39
Triadimefon	−4.62	−4.96	−5.12	−5.75	−5.81
Hexaconazole	−4.67	−3.5	−4.2	−4.44	−4.15
Bifonazole	−7.17	−7.3	−6.35	−7.45	−6.46
Miconazole	−6.67	−6.3	−6.31	−5.795	−6.5
Fadrozole	−7.78	−6.35	−6.43	−7.04	−7.7
Letrozole	−7.29	−7.04	−7.3	−7.27	−7.43
Prochloraz	−6.23	−5.7	−5.6	−6.0	−6.0
Clotrimazole	−5.8	−5.4	−5.3	−5.6	−5.7
Ketoconazole	−5.23	−5.0	−4.8	−5.2	−5.2

The hAr inhibitory activity of prochloraz, clotrimazole and ketoconazole (Table 2) could not be included in the corresponding correlation exercise because the crystal structure of human aromatase protein used in the study was unable to accommodate them. Consequently, an attempt was made to predict a suitable protein conformation to dock these bulky ligands. It is obvious that such a suitable protein conformation should have an enlarged binding site to accommodate these compounds. To achieve this, we collected the ketoconazole coordinate from lanosterol 14-alpha demethylase (pdb: 3LD6) and manually placed it within the ChAR binding site after protein alignment. A MDS of the ketoconazole bound ChAR for 1.5 ns was performed to attain an enlarged pocket for hAr. The MDS did not indicate any major change in backbone conformation of hAR except a bend of the I-helix to accommodate bulky ligands. The details of the MDS are provided in the supporting information (Figures 1S and 2S in the Supplementary Material). Among the trajectories generated, selection of the enlarged protein pocket was decided by measuring the distance between (1) M374–D309 (2) V370–W224 and (3) M374–G181. After a careful analysis a protein structure was selected with a distance of 20.54 Å between M374–D309, 14.1 Å between V370–W224 and 15.72 Å between M374–G181. The selected protein was considered to be an open state of hAr (OhAr) that can accommodate bulky ligands. By utilizing the OhAr we were able to retrieve the docking conformations of prochloraz, clotrimazole and ketoconazole at the protein active site but with lower docking scores relative to the scores obtained for other compounds at ChAr. We hypothesized that this behaviour may arise due to subtle changes in protein conformation and this decrease in docking scores at OhAr may also happen for other compounds. In order to confirm this we selected seven compounds with high, moderate and low activity and docked them with OhAr. The docking scores of similar azoles are low in the open state of human aromatase (OhAr) as compared to the closed state of human aromatase (ChAr) (Table 2). A significant correlation (r = −0.81) was observed for these seven compounds regarding docking scores and aromatase activity at OhAr.

Further analysis showed a good correlation (r = 0.84) between docking scores at closed and open state for these compounds, which tempted us to predict the closed state docking scores for prochloraz, clotrimazole and ketoconazole. The predicted docking scores at chAr for prochloraz (−6.23), clotrimazole (−5.8) and ketoconazole (−5.23) were next assembled with the 15 compounds to obtain a combined correlation value of −0.91 with pIC₅₀ values for the whole set of 18 compounds (Figure 6). The docking scores for prochloraz, clotrimazole and ketoconazole at ChAr for all of the bird and fish species were also predicted by using individual correlation equations derived from the 15 azoles, and have been depicted in Table 3.

Figure 6. Correlation between the docking scores of 18 azoles on the human aromatase and the corresponding human recombinant CYP19 inhibitory activity values (Table 2).

3.2.3 Docking of azoles on closed and open states of OomAr

Hinfray et al. [28] have designed a rainbow trout microsomal aromatase test to measure the potential adverse effects of xenobiotics on brain and ovarian aromatase activities in vitro. Among the series of 43 structurally diverse chemicals tested, they showed that imidazole fungicides (clotrimazole, imazalil and prochloraz), triazole fungicides (difenoconazole, fenbuconazole, propiconazole and triadimenol), and the pyrimidine fungicide fenarimol were able to inhibit brain and ovarian aromatase activities in a dose-dependent manner. The two steroidal compounds androstatrienedione and 4-hydroxyandrostenedione, which were also tested, showed a ratio of 60 between their brain and ovarian IC₅₀ values (Table 4).

Table 4. Docking scores (Kcal/mol) and inhibition of aromatase activities in brain and ovarian microsomes of O. mykiss [28].

Number	Compound^*	Docking score (close)	Brain IC50 [μM]	pIC50 brain	Ovarian IC50 [μM]	pIC50 ovary
1	Imazalil	−4.78	0.43	0.37	0.32	0.49
2	Propiconazole	−4.51	0.9	0.05	0.9	0.05
3	Triadimenol	−4.39	11	−1.04	26	−1.41
4	Fenbuconazole	−4.5	1.3	−0.11	0.21	0.68
5	Fenarimol	−4.12	6	−0.78	18	−1.26
6	Prochloraz	−6.0^**	1.3	−0.11	1	0
7	Clotrimazole	−5.7^**	0.011	1.96	0.016	1.8
8	Difenoconazole	−5.54	70	−1.85	29	−1.46
9	Androstatrienedione	−8.93	0.015	1.82	0.00025	3.60
10	4-OH androstenedione	−8.58	0.009	2.05	0.00015	3.82

^* Structures displayed in Figure 3

^** predicted docking score (see text).

An attempt was made to derive an equation allowing us to predict the aromatase activity of azoles and two steroidal compounds from their docking scores omAr. Due to the presence of prochloraz and clotrinazole in the data set (Table 4), in a first step, an open state of O. mykiss aromatase (OomAr) was homology modelled using the OhAR. An acceptable correlation was obtained (r = 0.80) between the docking scores at closed and open states of the azole compounds in Table 4, after exclusion of difenoconazole acting as outlier (results not shown). As a result, the closed state docking scores for prochloraz and clotrimazole were estimated to be −6 and −5.7, respectively (Table 4). The values of the docking scores at closed states were then used to compute a model for predicting brain IC₅₀ values. The obtained regression equation presented a modest correlation coefficient (r = −0.69) because clotrimazole (compound #7) and difenoconazole (compound #8) acted as strong outliers (Figure 7). Better results were obtained with the ovarian IC₅₀ values (r = −0.85) even if difenoconazole remained very poorly predicted with a residual value of −1.93. However, it is worth noting that, using a test based on H295R human adrenocortical carcinoma cells, Sanderson et al. [68] found an IC₅₀ value of 4 μM for difenoconazole, which shows a better agreement with the values calculated by both models for this compounds. Obviously, it should be necessary to have more experimental data obtained under the same experimental conditions to derive a more reliable model. Nevertheless, the results obtained in the present study are encouraging and confirm those previously found with the human recombinant CYP19 activity values. They show that docking scores for azoles and related compounds computed from homology models can be used for estimating the aromatase activity of these compounds.

Figure 7. Differences in the residuals obtained for the 10 brain pIC₅₀ values in O. mykiss (Table 4).

3.2.4 Comparative molecular docking analysis

In order to increase the structural diversity of our comparative exercise, 24 additional azoles and related structures (Figure 8) were docked on the aromatases of the different studied species. The obtained scores are listed in Table 5. Again, significant pairwise correlations were obtained between the docking scores for ChAr and those computed on the aromatases of the bird and fish species. The correlation coefficient values ranged from 0.82 for C. japonica to 0.86 for T. guttata (Figure 9).

Figure 8. Structures of 24 additional azoles and related compounds considered in this study.

Table 5. Docking scores (Kcal/mol) of 24 azoles on aromatases of the different studied species.

Compound	Human	C. japonica	G. gallus	T. guttata	O. mykiss
Anastrazole	−6.57	−6.29	−5.9	−6.41	−6.23
Azaconazole	−6.08	−5.47	−5.72	−5.6	−5.01
BitertanolSS	−5.81	−5.73	−5.34	−5.8	−4.94
Bromuconazole	−5.99	−5.82	−5.45	−5.78	−5.49
Climbazole	−5.1	−5.68	−5.91	−5.5	−5.48
Diclobutrazole	−5.99	−5.69	−5.7	−6.13	−5.2
Diniconazole	−5.47	−5.23	−6.16	−6.36	−5.84
Etaconazole	−5.98	−4.87	−5.35	−5.66	−5.56
Fluquinconazole	−4.95	−4.46	−4.84	−4.72	−4.58
Flutriafol	−6.04	−6.86	−6.01	−5.19	−6.32
Imibenconazole	−5.95	−6.12	−5.87	−5.97	−5.75
Metconazole	−5.6	−4.6	−5.94	−6.13	−5.41
Nuarimol	−6.07	−6.09	−6.29	−6.17	−5.96
Oxpoconazole	−3.12	−3.02	−2.93	−3.23	−3.22
Prothioconazole	−5.77	−4.95	−6.23	−4.98	−5.9
Quinconazole	−5.44	−6.13	−5.28	−4.92	−5.88
Simeconazole	−6.1	−5.8	−5.37	−6.05	−5.63
Tetraconazole	−5.34	−5.18	−4.79	−5.57	−4.96
Triflumizole	−4.98	−4.04	−5.5	−4.59	−4.03
Triticonazole	−5.74	−5.74	−5.84	−5.95	−5.84
Uniconazole	−6.1	−5.43	−5.97	−5.89	−6.5
Furconazole	−5.5	−5.49	−5.21	−5.73	−5.38
Ipconazole	−3.61	−3.72	−3.89	−4	−4.03
Pyrifenox	−5.47	−5.12	−5.78	−5.62	−5.21

Figure 9. Correlation matrix of the score values on the aromatases of the different studied species for the 24 compounds in Figure 8.

Finally, when the 46 azoles and related structures (Figures 3 and 8) considered in this study are analysed together (Table 1S in Supplementary Material), the interspecies correlation values remain broadly the same. Indeed, the correlation coefficient values obtained between the docking scores of ChAr and those computed on the aromatases of G. Gallus, T. guttata, C. japonica and O. mykiss equalled 0.82, 0.84, 0.83 and 0.87, respectively (Figure 3S in Supplementary Material). Figure 10 clearly shows that the azoles and related compounds generally behave similarly. Compound #15 in Figure 10 (letrozole) presents very high negative scores for all species. Letrozole is a human aromatase inhibitor commonly used in the treatment of hormonally responsive breast cancers in postmenopausal women. It has also gained recent attention in ovulation induction and infertility treatment [69]. As a result, letrozole can be found in the aquatic environment [70] in which it can induce adverse effects on the biota. Thus, letrozole was shown to impact reproduction, fecundity and endocrine functions in Japanese medaka (Oryzias latipes) [71]. A dose-dependent decrease in fecundity (>25 μg/L) and fertility (>5 μg/L) accompanied by histological changes were observed. At 625 μg/L, the medaka ceased spawning during the third week of exposure. Letrozole (>5 μg/L) reduced plasma vitellogenin levels in females in a dose-dependent manner. Hatchability and time to hatching were detrimentally affected (>5 μg/L), but no morphological deformities were observed. A dose-dependent increase in the proportion of genotypic F1 males was found (>5 μg/L) [71]. In the same way, Liao et al. [72] showed that exposure of medaka fish, at an early stage of sexual development, to continuous chronic environmentally relevant concentrations of letrozole, altered phenotypic sex development and reproduction in adults and skewed the sex ratio in progeny. Sex reversal was also obtained in the European pond turtle (Emys orbicularis) after their thermo-sensitive period by application of letrozole [73]. Sex inversion was also observed in G. gallus by application of letrozole and tamoxifen [74]. Deng and co-workers [75] showed that letrozole also inhibited the development of bone and medullary osteogenesis in prelay pullets by inhibiting the synthesis of oestrogen and its receptor. Impacts on territorial behaviours were also demonstrated in black redstarts (Phoenicurus ochruros) in experiments where letrozole was coupled with the androgen receptor blocker flutamide [76,77]. Much more information was obtained for the adverse effects of fadrozole on birds, an aromatase inhibitor also used against breast cancer and which also presents high negative scores in the studied species (compound #14 in Figure 10). In ovo injection of fadrozole resulted in the development of testes-like gonads in the majority of day-old genetic female chickens and turkey poults [78]. Near 100% of these birds showed masculine-type male genitalia at 1 day of age. Microscopic examination of the gonads of day-old genetic female chicks and female poults hatched from eggs treated with fadrozole showed the presence of atypical seminiferous tubules in about 17% and 73% of the individuals, respectively. Between the day of hatching and 6 weeks (chickens) or 12 weeks (turkeys), gonads in an increasing proportion of females from fadrozole-treated eggs showed normal appearing ovarian follicles [78]. These results were consolidated in other studies [17,44,47,79]. Fadrozole can also induce changes in the behaviour of birds. Male song sparrows (Melospiza melodia morphna) aggressively defend their territories year-round, except briefly during moult. During field experiments, Soma et al. [80] treated wild non–breeding male song sparrows with fadrozole for 10 days. This led to the decrease of several aggressive behaviours. The effects were reversed by 17β–oestradiol replacement. In another study, they showed that an acute exposure for one day decreased nonbreeding aggression [81]. Regarding aquatic ecosystems, a concentration-dependent reduction in fecundity was observed by Ankley et al. [65] in Pimephales promelas exposed for 21 days to 2–50 μg/L of fadrozole. A significant inhibition of brain aromatase activity was observed in both male and female fathead minnows. In females, this inhibition was accompanied by a concentration-dependent decrease in plasma β-oestradiol and vitellogenin concentrations. Histological analysis of ovaries from exposed females showed a decrease in mature oocytes and an increase in preovulatory atretic follicles. Exposed males showed significantly increased plasma concentrations of testosterone and 11-ketotestosterone as well as a significant accumulation of sperm in the testes [65].

Figure 10. Graphical display of the score values obtained for the 46 studied molecules (x-axis) on the aromatases of human, C. japonica, G. gallus, T. guttata and O. mykiss.

Oxpoconazole (compound #33) shows the lowest negative score values for all the species (Figure 10). Unfortunately, they cannot be confronted with in vivo experimental data due to a lack of information on the potential ecotoxicological effects of this fungicide on aquatic and terrestrial biota. However, while the scores for ipconazole (compound # 43 in Figure 10) are broadly the same that those of oxpoconazole, the limited information available on this fungicide tends to prove that it is of limited toxicity for fish and bird species. Ipconazole is practically non-toxic to birds on a dietary basis, with an LC₅₀ of >5620 ppm to both bobwhite quail and mallard duck. Studies on reproduction led to no observed effect concentration (NOEC) values of 50 ppm for the bobwhite quail and 200 ppm (the highest test concentration) for the mallard duck. Early-life stage exposure of ipconazole to newly fertilized fathead minnows showed an LC₅₀ > 2.9 mg/L for length and an LC₅₀ = 1.63 mg/L for weight [82]. Fenarimol (compound #46) shows score values around −4 kcal/M for all the species (Figure 10). Ankley et al. [30] showed that fenarimol reduced fecundity of P. promelas in a 21-day reproduction test at 1.0 mg/L. The fungicide affected several aspects of endocrine function in vivo; however, the suite of observed effects did not reflect either aromatase inhibition or androgen receptor antagonism. Chronic toxicity testing on bobwhite quail and mallard ducks resulted in no observed adverse effect concentrations (NOAECs) at the highest concentrations tested, which were 300 mg a.i./kg diet and 250 mg a.i./kg diet, respectively. No growth and reproductive effects were observed at these concentrations [83].

This comparative exercise shows that the docking behaviour of the azoles and related chemicals on the aromatase of human, G. Gallus, T. guttata, C. japonica and O. mykiss is broadly the same. In addition, a link can be made between the values of the docking scores and the adverse effects observed in these vertebrates.

4. Conclusions

The broad application of azoles in human and veterinary medicine and in agriculture has led to the contamination of the aquatic and terrestrial ecosystems by these chemicals. In silico approaches have proven their interest for assessing the toxicity of these chemicals [84–86], which are potential endocrine disruptors. While there is an increasingly number of valuable in vitro and in vivo studies aimed at estimating the aromatase inhibition potential of azoles to different fish species, the number of experimental studies on terrestrial species, including birds, is very limited. This is particularly annoying because numerous birds, especially grassland species, are in contact with the azole fungicides widely spread on a large number of crops. To overcome this problem, a stepwise strategy, based on homology modelling, was applied.

First, we took advantage of studies on azole drugs acting as strong aromatase inhibitors, which have been used on different species of birds to investigate physiological and ethological functions. This information was completed by the data on the very few azole fungicides for which effects on the endocrine system of birds have been investigated. Lastly, in vitro test results on series of azoles tested for their aromatase activity were also considered. In this context, all of the in vitro systems including fish cell lines were considered.

Because the crystal structure of human aromatase is available, as are the sequences of the aromatase proteins for the Japanese quail (C. coturnix japonica), the chicken (G. gallus), the zebra finch (T. guttata) and the rainbow trout (O. mykiss), homology models were derived for these species. These models were developed to explore the binding sites of aromatase enzyme and the critical interactions of the azoles with them, to estimate whether they behave in a similar way and also to evaluate whether the obtained binding score values can be related to the available aromatase inhibition data. High sequence similarity values were observed between the bird and fish aromatases and the human aromatase. Using a large set of azoles and related compounds we have shown that they behaved broadly similarly. We also demonstrated that there was an acceptable level of correlation between the binding score values and the aromatase inhibition data. This means that the homology models derived for bird species can be used to approximate the potential inhibitory effects of azoles on their aromatase.

Acknowledgements

The authors wish to acknowledge the two anonymous referees for their valuable comments and advice. This research was partly supported by the Pesticide research program (M6P) funded by the French Ministry of Ecology, Sustainable Development and Energy and the National Water and Aquatic Environment Institute (ONEMA).

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

^† Presented at the 8th International Symposium on Computational Methods in Toxicology and Pharmacology Integrating Internet Resources, CMTPI-2015, 21–25 June 2015, Chios, Greece.