Elucidating structure and dynamics of glutathione S-transferase from Rhipicephalus (Boophilus) microplus

Abstract Rhipicephalus (Boophilus) microplus is tick parasite that affects the cattle industry worldwide. In R. (B.) microplus, acaricide resistance develops rapidly against many commercial acaricides. One of main resistance strategies is to enhance the metabolic detoxification mediated by R. (B.) microplus glutathione-S-transferase (RmGST). RmGST detoxifies acaricides by catalyzing the conjugation of glutathione to acaricides. Although structural and dynamic details of RmGST are expected to elucidate the biologic activity of this molecule, these data have not been available to date. Thus, Molecular Dynamics simulations were employed to study ligand-free RmGST at an atomic level. Like other m-class GSTs, the flexible m loop (m1) of RmGST was observed. M1 seems to shield the active sites from the bulk. A RmGST dimer is stabilized by the lock-and-key motif (F57 as “key”) and hydrogen bonds of R82-E91 and R82-D98 at the dimer interface. Without substrates, conserved catalytic Y116 and N209 can interact with V112, G210 (for Y116) and F215 (for N209). Overall, most residues involving in RmGST function and stability are similar to other m-class GSTs. This implies similar structural stability and catalytic activity of RmGST to other GSTs. An insight obtained here will be useful for management of acaricide resistance and tick control. Communicated by Ramaswamy H. Sarma


Introduction
Ticks are blood-sucking ectoparasites that can transmit etiologic agents of human and animal diseases, and that impact livestock production worldwide (Narasimhan et al., 2021). Acaracides are a frontline tool for tick control (Walker, 2014), however the continuous use has led to tick resistance to acaricides, and importantly, the accumulation of poison residues in cattle products and the environment (Graf et al., 2004). Tick resistance to acaricides has become a major problem throughout the world. There are numerous reports of tick resistance to commonly used pesticides such as organochlorines (Miller et al., 2001), organophosphates (Klafke et al., 2017), amidines (Miller et al., 2001;Rodriguez-Hidalgo et al., 2017), synthetic pyrethroids (Klafke et al., 2017;Miller et al., 2001;Rodriguez-Hidalgo et al., 2017), macrocyclic lactones (Klafke et al., 2017;Rodriguez-Hidalgo et al., 2017), and phenylpyrazoles (Klafke et al., 2017). Such acaricide failures cause economic losses due to losses of livestock production. Vaccination is a promising and sustainable approach to counteract acaricide resistance in ticks, and the search for novel anti-tick vaccine candidates is ongoing. However, vaccination has not yet reached optimum performance, because such products are not yet sufficient for simultaneous protection against multiple tick species in diverse geographical regions (Schetters et al., 2016).
Rhipicephalus (Boophilus) microplus, the southern cattle tick, is one of the most widely distributed ticks worldwide, including Asia. This tick species is a competent vector of pathogens that cause babesiosis (Babesia bovis and Babesia bigemina) and anaplasmosis (Anaplasma marginale), which are considered the two most impactful tick-borne diseases of cattle worldwide (Shahein et al., 2013). Chemical pesticides serve as the primary method for controlling R. (B.) microplus, which results in drastic increases in the prevalence of acaricide resistance among these tick populations (Dzemo et al., 2022). R. (B.) microplus was reported to be resistant to pyrethroid, organophosphate, amitraz, fipronil, and macrocyclic lactone acaricides (Dzemo et al., 2022). Resistance in ticks is an adaptation which results from intensive exposure to acaricides (Aguilar et al., 2018). The long-term use of acaricides increases the frequency of naturally present resistance genes among tick populations, leading to the reduction of acaricide efficacy (Kumar, 2019). Increased detoxification and target site insensitivity have been reported to contribute to an increase in acaricide resistance (Shahein et al., 2013). The primary mechanism of R. (B.) microplus resistance to acaricides is metabolic detoxification mediated by cytochrome P450, esterase, and glutathione-S-transferase (GST) (Chevillon et al., 2007;Dzemo et al., 2022;Le Gall et al., 2018;Li et al., 2004;Miller et al., 1999). Such enzymes hydrolyze, oxidize, and conjugate toxicants. Conversely, some xenobiotics are transported to the exterior of cells by ATP-binding cassette transporter (ABC transporter) (Le Gall et al., 2018;Pohl et al., 2012). GSTs of R. (B.) microplus have attracted attention because of their involvement in resistance mechanisms of many common acaricides such as permethrin, amitraz, and ivermectin (Daborn et al., 2002;Li et al., 2004;Pasay et al., 2009). Notably, GST was also found to be immunostimulatory to cattle and was tested as an anti-tick vaccine candidate (Parizi et al., 2011;Shahein et al., 2013;Veerapathran et al., 2009). Although the inhibition of GST function was suggested as one of novel tick control strategies (Hernandez et al., 2018), the information on GST structure and function remains limited. Specifically, to the best of our knowledge, the precise nature of GST structure, dynamics at the molecular level, and their role in detoxification have been unavailable.
GSTs are multifunctional enzymes that protect cells against chemical toxicity and oxidative stress Hernandez et al., 2020), thus contributing to detoxification of acaricides (Hernandez et al., 2018). In mammals, GSTs are categorized into five families, namely alpha (a), mu (l), pi (p), theta (h), and sigma (r), based on their sequence similarity and cross-immunoreactivity (Le Gall et al., 2018). These different GST classes reflect the chromosomal location of their genes (Shahein et al., 2013)  GSTs catalyze the nucleophilic conjugation of glutathione (GSH) to a wide range of electrophilic compounds such as endo-and xenobiotics (Gulc¸in et al., 2018;Muller et al., 2008;Temel et al., 2018). A thioether linkage is formed between the sulfur of GSH and xenobiotics, which makes a resulting conjugate more water-soluble and easy to be excreted from cells (Le Gall et al., 2018). In insects, the GST system is a major mechanism that detoxifies insecticides and reduces oxidative stress (Hernandez et al., 2021;Ku et al., 1994). In mosquitos, the depletion of GST increases superoxide production and apoptosis (Chen et al., 2012;Laborde, 2010). The functional forms of GSTs are dimers comprising identical or structurally related subunits. Each subunit (chains A and B) contains two domains: domains I (N-terminal) and II (C-terminal) ( Figure 1A). The N-terminal domain (I) is rather conserved, while domain II is variable (Deponte, 2013). Each chain contains two active sites that are a highly conserved G-site (GSH binding site) in domain I and a diverse H-site (hydrophobic co-substrate binding site) in domain II ( Figure  1A). The different H-sites across GST families reflect the great variability of co-substrates.
In R. (B.) microplus, the 223 amino acid-long RmGST is related to the mammalian m class GST (He et al., 1999;Shahein et al., 2013). The N-terminus shows higher sequence similarity than the C-terminus, which is common in GST families (Sinning et al., 1993). Although several tick GST studies have been reported (Guerrero et al., 2012;Reddy et al., 2011;Wei et al., 2001;Yang et al., 2017), there is a lack of structural and dynamic properties of RmGST. Such information is necessary for understanding the nature of GSTs that affect acaricide resistance in R. (B.) microplus. Thus, in this work, the structure and function of RmGST were investigated. To obtain molecular insights, Molecular Dynamics (MD) simulations were performed. MD simulations have been successfully used to reveal the behaviors of other tick proteins (Pongprayoon et al., 2020;. This insight is expected to provide better understanding of an important acaricideresistance mechanism, which may in turn facilitate development of novel approaches to tick control.

Materials and methods
The three dimensional structure of RmGST (Accession No. AAD15991) was constructed using MODELLER software (Eswar et al., 2006). The crystal structure of chain B of human m class GST (hGSTM3) was used as a template (PDB code: 3GTU) (www.rcsb.org) because this template has the highest sequence identity (68.66%) to RmGST. The quality of RmGST was checked using RAMPAGE (http://mordred.bioc.cam.ac.uk/ $rapper/rampage.php) (Biasini et al., 2014). The models show good proportions of residues in favoured (98.48%) and outlier (1.52%) regions ( Figure S1 in supplementary information). A protein was placed in a cubic box (with a dimension of 5 Â 5 Â 5 nm 3 ) and solvated with TIP3P water molecules (22,500 molecules). The protonation states of charged amino acids were set at physiologic pH. Counter ions (4 Na þ ions) were added to neutralize the system. The energy minimization of 1000 steps was performed to remove close contacts, using steepest descent algorithm followed by 100-ps NVT and NPT runs. Two copies with different random seeds of 500-ns production runs were performed (RmGST1 and RmGST2). The endings of "1" and "2" refer the simulations 1 and 2. The data shown are the average between the two simulations. Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) were calculated using an initial structure at 0 ns from each production run as a reference.
The GROMACS 5 package (www.gromacs.org) (Lindahl et al., 2001) was employed with Amber99SB force fields. The particle mesh Ewald (PME) techniques (Darden et al., 1993) with a Fourier spacing of 0.12 nm and a short range cut-off of 1 nm were used for electrostatic treatment. The simulations were conducted in the constant number of particles, pressure, and temperature (NPT) ensemble. The Berendsen algorithm at 1 bar with a coupling constant s p ¼ 1 ps was used for pressure coupling. The temperature of the protein, ligand, and solution were coupled separately using the vrescale thermostat (Berendsen et al., 1984) at 300 K with a coupling constant s t ¼ 0.1 ps. The time step of 2 fs was used for integration. The coordinates were recorded every 2 ps.
All data were analyzed by GROMACS tools and in-house codes. Gmx hbond with default parameters. The hydrogendonor-acceptor cutoff angle was set to 30 and the cutoff radius (X-acceptor) was 0.35 nm. VMD was used for visualization and graphic images (Humphrey et al., 1996).

Results and discussion
At the beginning, the structural flexibility is determined by C-alpha RMSDs and RMSFs. C-alpha RMSDs of whole RmGST1 and RmGST2 are in a range of $0.3-0.4 nm (Figure 2A). Chains A and B are expected to have lower RMSDs than the whole protein ($0.2-0.3 nm in Figure 2A). In addition, RMSFs indicate clearly that residues 35-45 (a turn connecting B1 and A2; m1) are the main contributor for the protein dynamics in both chains ( Figure 1A and 2B); however, residues 110-130 (part of A4 and its connected loop) in RmGST2 were also found to be mobile. Both regions were reported to participate in the active site. The high fluctuation of m1 and residues 110-130 here may be because of the absence of substrate. Principal Component Analysis (PCA) was also performed to further explore protein dynamics. PCA can confirm the high mobility of residues 35-45 (m1) on A2 region in both cases ( Figure 3A). Only the motion obtained from the first principal component 1 (PC1) is considered here, because the first principal component 1 (eigenvector) accounts for major motions ( Figure S2 in supplementary information). The m1 mobility was also reported to be the most prominent motion in ligand-free vertebrate GSTs (Ricci et al., 1996;Stella et al., , 1999. Such m1 is the "m loop" which is the same as the distinctive and conserved feature of the mammalian m class (Dourado et al., 2010;Hearne & Colman, 2006). The m loop, together with the C-terminal tail and loop, connected A4 and A5 from the active site cleft (their location can be seen in Figure S3), where m1 acts as a lid that shields the active site from the bulk (Perbandt et al., 2004). No hydrogen bonding of m1 with the C-terminus and A4 was captured. This can be one of reasons for high m1 flexibility. The dynamics of the m1 region was evaluated by the distance between the tip of m1 (P38) and P118 on A4 ( Figure  3B). The fluctuated distances observed confirmed the high flexibility of the m1 zone. Comparing between the two chains, chain A showed a more mobile m1 than that of chain B ( Figure 3B). The mobility of the m1 loop may influence the ligand-binding affinity, warranting further study.
Interactions at the dimer interface were also investigated. The hydrophobic lock-and-key feature is characteristic for GST dimerization, which is common in a-, m -, and p-classes (Sinning et al., 1993). Phenylalanine or tyrosine (the "key") in one subunit is wedged into a hydrophobic pocket of the other unit (the "lock") (Alves et al., 2006;Hegazy et al., 2004;Sayed et al., 2000). This lock-and-key motif was also reported to be important for catalytic activity (Codreanu et al., 2005;Hornby et al., 2002). For RmGST, F57 acts as the "key" which fits into the "lock" provided by F104, F138, and Y141 on A4 and A5 ( Figure 1A and 4A). The distances between the "key" (F57) with the hydrophobic pocket-lining residues are measured in Figure 4A. The shorter distances of F57-F138 and F57-Y141 in chain B demonstrate the closer lock-and-key packing of chain B than chain A ( Figure 4A). Furthermore, in most chains, the constant distances between phenyl carbons (CZ) of F104 and F138 illustrate the F104-F138 stacking pose ( Figure 4B, position 2, and 4 C). Unlike others, F104 of chain B in RmGST1 seems to move away from the "lock" site ( Figure 4B, position 1), whereas the rest shift their phenyl ring toward the "lock" pocket ( Figure 4B, position 1!2). The drift of F104 in RmGST1 may reduce the binding affinity of F57 to the hydrophobic "lock" pocket. Further study of the ligand-GST interaction is required to confirm the F104 reorientation.
Apart from the inter-subunit hydrophobic lock-and-key interaction network, the inter-hydrogen bonds at the dimer interface are also crucial to stabilize a protein dimer. Approximately 5-9 hydrogen bonds are employed for GST dimerization ( Figure 5A). At the dimer interface, A3 (residues 72-85) in one subunit and A4 (residues 90-115) from the other seem to be the main contributors for stabilizing a GST dimer. For each chain, A3 and A4 can form 2-3 inter-hydrogen bonds with their adjacent chain ( Figure S4). This interaction network is formed by R82 on A3 interacting with E91 and D98 on A4 from the neighboring subunit ( Figure 5B). It can be seen in Figure 5C that E91 and D98 play a role in trapping R82 from the adjacent chain. This cluster produces 4-8 inter-ionic hydrogen bonds in total, which seems to be the major polar forces accounting for GST dimer formation. The similar electrostatic network between chains was also identified in chicken m GST (Sun et al., 1998).
In the case of binding sites, RmGST and hGSTM3 share highly conserved ligand-binding residues (Patskovsky et al., 1999). Compared to hGSTM3, key polar residues for GSH binding (G-site) in RmGST are conserved (Y7, R43, W46, K50, Q72, S73, Y116, and N209) ( Figure 6A). In contrast, more sequence variance is found in the xenobiotic site (H-site). It was reported that I13, L16, L114, Y119, I211, and W218 are the H-site-lining residues in hGSTM3 (Patskovsky et al., 1999), while most of these residues are preserved in RmGST (corresponding to I10, L13, Y116, and I208). However, positions 111, 112, and 215 in RmGST are replaced by W111, V112, and F215 ( Figure 6B). The replacement by tryptophan at position 111 can enhance the aromaticity of the H-site pocket, which may facilitate the binding of a wider range of aromatic compounds, but further study is required. Furthermore, a pocket-lining tyrosine (Y119 in hGSTM3) was reported to be the active residue in a, m, p, and s classes (Deponte, 2013). This tyrosine is also found in RmGST (Y116). In hGSTM3, the template in this work, not only Y119 but also N212 were reported to play an important role in the catalytic activity (Patskovsky et al., 1999). N212 is located at the C-terminus near a hydrophobic pocket, while Y119 is at the base of A4. These two residues are conserved in RmGST (N209 and Y116) ( Figure 6D). In an absence of substrate, Y116 employs its backbone to hydrogen bond with the V112 backbone and leaves its hydroxyl group to interact with G210 ( Figure 6C). N209 is trapped by the F215 backbone ( Figure 6C). Y10, reported to stabilize a thiol group of GSH in GST, is also seen here as Y7 ( Figure 6A) (Angelucci et al., 2005;Armstrong, 1991;Parsons & Armstrong, 1996). Y10 was suggested to either point towards the catalytic site or shift away and form the cation-p interaction with R21 (R18 in RmGST) (Angelucci et al., 2005;Johnson et al., 2003). Herein, both residues were conserved (Y7 and R18), but no Y7-R18 interaction was identified because Y7 was directed to the active site throughout a course of simulations ( Figure 6A). Further experimental study is required to explore the conformational change of Y7 in RmGST.
Seemingly, most residues involving in the binding of GSH, and xenobiotic and dimerization in RmGST, are similar to other m-class GSTs (Armstrong, 1991;Codreanu et al., 2005;Hornby et al., 2002). This implies similar structural stability and catalytic activity of RmGST to other GSTs. However, some variations of residues inside the active site are captured, which should influence the catalysis. The understanding of RmGST function will be important for tick control. Future work should include the study of ligand-bound RmGST.

Conclusions
In this work, the structural and dynamic properties of RmGST were investigated for the first time. Seemingly, RmGST shows similar features to m-class GSTs. The RmGST m1 was the most mobile, similar to many m-class GSTs . The high mobility of m1 herein supports the ligand-induced conformational change of RmGST, as reported with other GSTs (Stella et al., 1999). To stabilize its dimer, RmGST employs the lock-and-key motif and electrostatic interactions at the dimer interface. In the case of the catalytic sites, like other m GSTs, polar residues in the G-site were conserved, while the hydrophobic H-site in RmGST seems to show more aromaticity due to a presence of tryptophan. A presence of tryptophan may accommodate a wider range of aromatic insecticides in R. (B.) microplus. Overall, conserved key active site-lining residues in RmGST suggests a similar catalytic mechanism to other m-class GSTs. Since RmGST is part of acaricide resistance mechanism, therefore the findings obtained here can serve as a base for developing and designing potent pesticides. Nonetheless, an additional study is needed, especially with respect to the binding mechanism of RmGST to glutathione and toxicants. RmGST was not only suggested as an anti-tick vaccine candidate, due to its ability to elicit a bovine immune response, but this molecule also has a role in the metabolics of acaricide resistance. The acaricide binding to RmGST will be our further work. An insight into structural and dynamic properties of RmGST thus can be useful for both anti-tick vaccine design and management of acaricide resistance.