Comparative analyses and structural insights of new class glutathione transferases in Cryptosporidium species

Cryptosporidiosis, caused by protozoan parasites of the genus Cryptosporidium, is estimated to rank as a leading cause in the global burden of neglected zoonotic parasitic diseases. This diarrheal disease is the second leading cause of death in children under 5 years of age. Based on the C. parvum transcriptome data, glutathione transferase (GST) has been suggested as a drug target against this pathogen. GSTs are diverse multifunctional proteins involved in cellular defense and detoxification in organisms and help pathogens to alleviate chemical and environmental stress. In this study, we performed genome-wide data mining, identification, classification and in silico structural analysis of GSTs in fifteen Cryptosporidium species. The study revealed the presence three GSTs in each of the Cryptosporidium species analyzed in the study. Based on the percentage identity and comprehensive comparative phylogenetic analysis, we assigned Cryptosporidium species GSTs to three new GST classes, named Vega (ϑ), Gamma (γ) and Psi (ψ). The study also revealed an atypical thioredoxin-like fold in the C. parvum GST1 of the Vega class, whereas C. parvum GST2 of the Gamma class and C. melagridis GST3 of the Psi class has a typical thioredoxin-like fold in the N-terminal region. This study reports the first comparative analysis of GSTs in Cryptosporidium species.


Methods
Species and database. Cryptosporidium species genomes that are available for public use at the Cryptosporidium database or CryptoDB 29 (https ://crypt odb.org/crypt odb/app; release 48 beta, 27 August 2020; accessed on 14 September 2020) and at National Center for Biotechnology information (NCBI) 30 (https ://www.ncbi.nlm. nih.gov/datas ets/genom es/?txid=5806; accessed on 14 September 2020) were used in the study. The Cryptosporidium pathogens examined in this study include ones from both humans and other mammals (Table 1).

Genome data mining, identification and classification of GSTs. Cryptosporidium species genomes
available at CryptoDB 29 were mined for GSTs. Two different methods followed for GST mining. First, the genomes of Cryptosporidium species were mined using the term "glutathione transferase". Second, the species genomes were blasted with GST proteins from Homo sapiens (protein ID: P08263) 53 and C. parvum Iowa II (protein ID: EAK89476.1) 14,38 . The BLASTP mined proteins revealed a range of apicomplexan species which were filtered out to show only Cryptosporidium species. The hit proteins were then collected and subjected to protein family analysis using the Pfam 54 and InterPro 55 programs. The results were analyzed and the hit proteins that were classified as GST by Pfam (PF14497, PF13417 and, PF17172) 54 and InterPro (IPR036282, IPR004045 and IPR010987) 55 were selected.
For the collection of more hits, Cryptosporidium species genomes available at NCBI database 30 was blasted with two GST proteins from C. andersoni 30847 (cand_012830 & cand_023790) and from C. meleagridis UKMEL1 (CmeUKMEL1_05845) that were collected from CryptoDB 29 . The hit proteins were screened for GSTs following the method described above.
A final total count was presented by deleting the duplicated GSTs. The selected GSTs were then grouped into different classes or groups based on their percentage identity, following the conventional criterion of less than 25-30% identity being a new class 20 www.nature.com/scientificreports/ Analysis of homology. The percentage identity between GSTs was deduced using Clustal Omega 56 . The full-length GSTs were subjected to Clustal analysis which produced the percentage identity amongst each of the proteins as matrix identity results. These results were laid out in an Excel spreadsheet where the results were analyzed to identify the percentage identity between GSTs.
Collection of different GST classes' protein sequences. For comparative analysis, GST protein sequences belonging to different GST classes were collected using multiple methods to build a library for phylogenetic analysis. On the European Molecular Biology Laboratory (EMBL) site 57 , GSTs sequences that are placed under the GST superfamily (IPR040079) were retrieved. The GST classes namely CLIC (IPR002946), Alpha (IPR003080), Mu class (IPR003081), Pi (IPR003082), Omega (IPR005442), Zeta (IPR005955) and Sigma (IPR003083) were collected under EMBL. More sequences were obtained through text search using the UniProt protein knowledge base 58 . A specific GST class was searched on the site and the hits obtained were further verified using Pfam 54 and InterPro 55 to ensure uniformity with the GSTs collected from the EMBL site 57 . The remaining GSTs that were not in the databases were retrieved from published articles. The Cryptosporidium species GST sequences along with protein sequences of different GST classes used in the phylogenetic analysis are presented in Supplementary Dataset 1.
Phylogenetic analysis. The GST sequences in supplementary dataset 1 were used to make a phylogenetic tree for inferring their evolutionary relationship. First, all the GST protein sequences were aligned by MAFFT v6.864 embedded on the Trex-online server 59 . Then, the alignment was automatically submitted to the server for inferring the tree with different models and the optimized tree was selected. Finally, the tree was submitted to iTOL for viewing and annotation 60 . Thioredoxin from Oryctolagus cuniculus (protein ID: P08628) was used as an outgroup.
For the construction of the phylogenetic tree of the Cryptosporidium GST proteins, the protein sequences were aligned using MUSCLE software 61 embedded in MEGA7 62 . The evolutionary history was inferred by using the maximum likelihood method with 100 bootstrap replication based on the JTT matrix-based model 63 . Evolutionary analyses were conducted in MEGA7. Cellular localization and transmembrane helices prediction. Cellular localization of GSTs was predicted using the Bologna Unified Subcellular Component Annotator (BUSCA) 64 . BUSCA is the latest, accurate program available for the prediction of proteins' subcellular localization; it integrates different computational tools such as identifying signal and transit peptides (DeepSig and TP-pred3), GPI-anchors (PredGPI) and transmembrane domains (ENSEMBLE3.0 and BetAware) with tools for discriminating subcellular localization of both globular and membrane proteins (BaCelLo, MemLociand SChloro) 64 . The outcomes of these different programs were processed and integrated to predict subcellular localization of both eukaryotic and bacterial proteins 64 . Prediction of transmembrane helices in GSTs was done using TMHMM Server v. 2.0 65 . This program is well known for its high degree of accuracy in the prediction of transmembrane helices and discrimination between soluble and membrane proteins.

Template identification.
To construct 3D models of proteins, reference protein structures previously solved by crystallization or Nuclear Magnetic Resonance are needed. These would serve to simulate not only the fold of a protein but also a full atom model to build. These proteins are referred to as templates. Either single or multiple templates can be used in constructing the 3D model of a protein 66 . In this study, three different web servers, namely NCBI BLAST (v2.10.1) 67 , i-TASSER (v5.1) 68 and PHYRE (v2.0) 69 , were consulted to identify the Table 1. Cryptosporidium species used in the study and their major host specificity.
Protein sequence alignment for modeling. T-COFFEE webserver 70 was used for aligning the GST proteins and the template sequences. The aligned files were downloaded in FASTA format and modified to generate files to be used for protein modelling 71 .
Protein modeling, optimization and validation. The MODELLER v9.21 program 71 was used to build GST models. Multiple structures were produced by Modeller 9.21. The model with the best DOPE assessment was selected as the output structure to be used. The structures modeled were viewed using PyMOL 72 . The model for each GST was then subjected to evaluation for stereochemistry and energetic quality at the Structural Analysis and Verification Server (http://servi cesn.mbi.ucla.edu/SAVES /) and ProsaII (https ://prosa .servi ces.came. sbg.ac.at/) 73 . Based on the validation results, the protein models were then refined on the GalaxyWeb Refiner server 74 . After refinement, the models were again subjected to evaluation and validation using programs such as ERRAT 75 , Verify3D 76 , PROCHECK 77,78 , and RAMPAGE 79 and ProsaII 73 .

Results and discussion
Two different sizes of GSTs present in Cryptosporidium species. Genome data mining of 15 Cryptosporidium species revealed the presence of 3 GST genes in each of the species genomes ( Table 2). The presence of more than one GST gene is common in eukaryotic species 23 . Among 45 GSTs, 30 were found to have the characteristic GST motifs 20,27 , such as the N-terminal domain, which houses the G site, and C terminal domain, which determines the substrate specificity (H-site) ( Table 2 and Fig. S1). The remaining 15 GSTs have one of the characteristics GST motifs indicating either these sequences are diverse or fragmented or not properly annotated ( Table 2). These GSTs were considered incomplete and were not included for further analysis unless indicated. Future genome editing and better gene prediction programs will help in getting the complete sequences for these GSTs and possibly predicting characteristic N-and C-terminal motifs. In total, 30 GSTs were taken for further analysis. Analysis of GST protein sizes revealed the presence of two different lengths of GSTs in Cryptosporidium species (Table 2). One type of GST protein is shorter in size with amino acids ranging between 157 and 268, and another type of GST protein is longer in size, with amino acids ranging between 373 and 466 ( Table 2). GSTs from Cryptosporidium species seem to be the longest in amino acid length, as most of the GSTs reported in other organisms to date are 200-250 amino acids in length 23 . Furthermore, it can be noted that the addition in length is found only on the outer N-and C-terminal regions, with the center of the protein containing the GST-superfamily domains ( Table 2). In order to assess whether Cryptosporidium species GST proteins are indeed properly annotated gene products, we further analyzed the gene structure. Interestingly, all the longer GSTs had a single exon, thus no introns, but shorter GSTs were the products of 1-4 exons ( Table 2). This could be indicative of shorter GSTs being prone to having multiple isoforms owing to gene shuffling. Due to presence of these multiple introns, the production of more diverse short GSTs can be expected compared to longer GSTs 80 .
Cryptosporidium species GSTs are cytosolic in nature. Most of the GSTs identified in organisms are cytosolic in nature, with the exception of GSTs belonging to the classes MAPEG and Kappa (mitochondrial) ( Table S1). In order to identify the cellular localization, we subjected Cryptosporidium species GST protein sequences to the TMHMM Server v. 2.0 for the prediction of transmembrane helices in their structure 65 and the BUSCA server 64 for identifying possible localization in a cell. TMHMM prediction revealed that none of the Cryptosporidium species GSTs had transmembrane helices, indicating they were soluble and thus possibly cytosolic (Table S2). To authenticate our results, we also subjected 395 GSTs belonging to 17 different classes to TMHMM prediction (Table S3). The TMHMM predicted the presence of no transmembrane helices in previously designated cytosolic GSTs, whereas transmembrane helices were predicted for previously designated microsomal GSTs (Table S3). This indicated that the TMHMM results on the prediction of no transmembrane helices in Cryptosporidium species GSTs were in agreement with previous annotations. Furthermore, BUSCA indicated that all 30 Cryptosporidium species GSTs were cytosolic (Table S4). Based on these in silico results, we concluded that the 30 Cryptosporidium species GSTs were cytosolic in nature.

Cryptosporidium species GSTs belongs to new classes. Phylogenetic analysis of Cryptosporidium
species GSTs revealed that the 30 GSTs could be grouped into three different groups (Fig. 1). The shorter GSTs were grouped together (Group 1) and so were the longer GSTs (group 2). Interestingly, despite the short amino acid length, four GSTs diverged from these two groups (Group 3) (Fig. 1). Analysis of the amino acid percentage identity among Cryptosporidium species GSTs further confirmed that they indeed belonged to three different groups. Group 1 GSTs shared an amino acid percentage identity of 54-100%, whereas groups 2 and 3 shared identities of 48-100% and 42-71%, respectively. Group 3 GSTs had 13-21% identity with Group 2 GSTs and 14-22% identity to Group 1 GSTs. The percentage identity between Groups 1 and 2 was 17-25%. This indicates that all three groups of Cryptosporidium species GSTs indeed belonged to three different classes as the percentage identity between these groups was below 25-30%, qualifying them to be their own class 20,23,26-28 .
Although the above results clearly indicated that Cryptosporidium species GSTs belong to three different groups, it was still not clear whether they fell under one of the GST classes described in the literature (Table S1). Thus, the comprehensive phylogenetic analysis of proteins belonging to 17 known GST classes and Cryptosporidium species GSTs was carried out (Fig. 2). Phylogenetic analysis revealed that Cryptosporidium species GSTs did not align with any of the 17 pre-existing GST classes and formed three new groups (Fig. 2) www.nature.com/scientificreports/ indicates that Cryptosporidium species GSTs belong to three different new GST classes. Thus, we named groups 1, 2 and 3 of Cryptosporidium GSTs Vega (ϑ), Gamma (γ) and Psi (ψ), respectively. A point to be noted is that all the GST proteins aligned together as per their GST class on the phylogenetic tree, indicating our phylogenetic analysis is correct and thus we conclude that Cryptosporidium species GSTs indeed belong to new GST classes.

Cryptosporidium parvum GST1 of Vega class has atypical thioredoxin-like fold. Identification
of three new GST classes in Cryptosporidium species in this study necessitated examination of the structural aspects of these new classes to see if any deviations or novel folds might be present, compared to the canonical structure of GSTs 20,27 . Analysis of the primary structure revealed that all Cryptosporidium species GSTs have Nand C-terminal regions characteristic of GSTs that usually contain a G-site and H-site 20,27 , respectively (Table 2 and Fig. S1). All GSTs have the highly conserved proline amino acid residue (Fig. S1) that is part of the cis-Pro loop responsible for connecting the N-and C-terminal regions in order to maintain the GST structural integrity 81 . It was observed from Fig. S1 that Psi class GSTs have a Tyr residue in the N-terminal domain in close proximity to the expected active site Tyr. The same was observed with the Vega class GSTs with the expectation of C. muris and C. baileyi. Vega and Psi GSTs have a few tyrosine residues in the N-terminal region, but they are not at a position that is considered part of an active site 20,27 (Fig. S1). Similarly, the majority of the Gamma class GSTs consist of an active site Tyr residue with the exception C. andersoni, C. baileyi, C. ryanae and C. bovis species. In these species, Phe replaces the active site Tyr residue. Mutagenesis studies have shown that the presence of Phe at the supposed position of the active site Tyr significantly reduces the catalytic activity. This highlights the critical role played by the active site Tyr in the catalytic activity of GST 82,83 . The effect of these mutations in the context of Cryptosporidium GSTs is yet to be studied. Multiple sequence alignments of Vega and Gamma GSTs revealed that amino acids in the N-and C-terminal regions of these GSTs are highly conserved (Fig. S1). For this reason, we selected C. parvum GSTs 1 and 2 www.nature.com/scientificreports/ (CpGST1 and CpGST2) as representative of the Vega and Gamma GST classes for structural analysis along with C. meleagridis UKMEL1 GST3 (CmGST3) for the Psi class. Structural analysis of the three GSTs was carried out using in silico homology modeling. The structural analysis was aimed at assessing only the secondary structural elements that are characteristic of GST proteins 20,27 . These GST models are not aimed to assess the binding affinities or the residues involved in binding to different ligands. In order to build 3D models we performed a template search at three different webpages, namely NCBI 67 , PHYRE 69 and I-TASSER 68 . The templates found were of low sequence identity but had relatively good coverage (Table S5). This was expected, since these GSTs are new. We then proceeded to build 3D models using a multiple template method, as this approach is known to improve the quality of homology models 84 . We built 3D models for all three GSTs, attempting single and multiple templates, while also using different combinations of the available templates listed in Table S5. The best 3D models with good quality closest to the templates were chosen for the structural analysis.
For each GST, 20 models were built using the MODELLER v9.21 program 71 . The best model evaluated by DOPE score was selected and subjected to structural quality analysis. The selected model for each GST was then refined on the GalaxyWeb Refiner server 74 and further subjected to structural quality evaluation using different programs such as ERRAT 75 , Verify3D 76 , PROCHECK 77,78 , RAMPAGE 79 and ProsaII 73 . The overall quality of the models was assessed by the combination of these programs' values and by comparing these with the templates' structural evaluation scores (Tables S6 and S7). The models generated for CpGST1 and CpGST2 were found to be of good quality, as different structural validation programs indicated that the quality of the model structures was close to the quality of the template structures (Tables S6 and S7). The model generated for CmGST3 had www.nature.com/scientificreports/ all parameters in acceptable range including Z-score of − 3.68 indicating the model is of good quality with the exception of Verify3D where 26% residues had an average 3D-1D score > 0.2 (Tables S6 and S7). The three GST models generated in the study, along with their corresponding sequence alignments with their templates, are presented in Figs. 3, 4 and 5. Structural analysis revealed the presence of 2β-sheets and 3α-helices in the N-terminal region and 6α-helices in the C-terminal region of CpGST1 (Fig. 3). The overall structure of CpGST1 at the N-terminal domain seems completely different compared to the canonical GST N-terminal domain 20,27 . The N-terminal region of CpGST1 did not have the typical thioredoxin-like fold, nor did it follow the βαβ-α-ββα motif; it was rather composed of two antiparallel β-sheets and 3α-helices (Fig. 3). It is rare to find GSTs that do not possess the conventional thioredoxin βαβ-α-ββα motif. Kappa class GSTs, which are mitochondrial GSTs, are the closest GSTs that do not follow the traditional thioredoxin fold but have still been found to carry out a similar molecular function  www.nature.com/scientificreports/ as conical GSTs [89][90][91] . This is also common for MAPEG GST and the mPGES-1 (microsomal ProstaGlandin E-Synthase type 1) subfamily of proteins, as they too are a group of structurally unrelated proteins with GSH transferase activities 23,91 . Because the GST superfamily shares such vast variations in terms of their structural conformation, this ααββα conformation of CpGST1 can be considered a unique Vega class feature.
In contrast to the CpGST1 model, the CpGST2 and CmGST3 models N-terminal domain follows the thioredoxin-like fold, which is characteristic of cytosolic enzymes in the GST superfamily 20,22,27 . The N-terminal domain was complete with 4β-sheets and 3α-helices following a βαβ and ββα arrangement with the two motifs linked by an α2 (Figs. 4 and 5). The C-terminal domain contains helices with each model CpGST2 and CmGST3 having a varying number of helices (Figs. 4 and 5). It has been suggested that an increase in the number of helices in the C-terminal domain, may allow for a broader substrate range and/or offer a deeper catalytic pocket that facilitates the conjugation of larger substrates 92,93 .

Conclusions
In this genomic era, in silico based comparative studies at genome level or at protein family level have become an important tool to uncover novel aspects in organisms. This study is such an example, where genomes of Cryptosporidium species were mined for glutathione transferases (GSTs), enzymes playing a key role in cellular defense and detoxification that are also a potential drug target against pathogens and metabolic disorders. Analysis revealed an interesting feature, namely the presence of two different sizes of GSTs (short and long) in these species. The longer GST proteins were found to be longer than the GSTs found in other organisms, with the size attributed to C-and N-terminal extensions. One of the major findings of the study is the identification of GSTs belonging to three new GST classes in Cryptosporidium species. In addition, Cryptosporidium parvum GST1 had an atypical thioredoxin fold in the N-terminal region with an αα-ββ-α motif rather than the typical thioredoxin-like fold with a βαβ-α-ββα motif. Future study includes functional and structural (X-ray or NMR) characterization of Cryptosporidium species GSTs. The study results serve as reference for future mining and annotation of GSTs Cryptosporidium species.
Received: 25 June 2020; Accepted: 28 October 2020 Figure 5. In silico structural analysis of Psi class representative Cryptosporidium meleagridis strain UKMEL1 GST3 glutathione transferase 3 (CpGST3). 3D model of CmGST3 (A) and its amino acid sequence alignment with templates (B). Secondary structural annotations were done as per modeled structure where α-helices and corresponding amino acids are colored in red while the β-sheets and their corresponding amino acids are colored in yellow. The active-site tyrosine and the cis-proline residues are boxed in purple and grey respectively. The template Protein Data Bank codes, 1K3Y-B and 1ZL9-A, represents GST protein crystal structures from Human (Alpha class GST) and Caenorhabditis elegans (Pi class GST). www.nature.com/scientificreports/