iTRAQ-Based Global Phosphoproteomics Reveals Novel Molecular Differences Between Toxoplasma gondii Strains of Different Genotypes

To gain insights into differences in the virulence among T. gondii strains at the post-translational level, we conducted a quantitative analysis of the phosphoproteome profile of T. gondii strains belonging to three different genotypes. Phosphopeptides from three strains, type I (RH strain), type II (PRU strain) and ToxoDB#9 (PYS strain), were enriched by titanium dioxide (TiO2) affinity chromatography and quantified using iTRAQ technology. A total of 1,441 phosphopeptides, 1,250 phosphorylation sites and 759 phosphoproteins were detected. In addition, 392, 298, and 436 differentially expressed phosphoproteins (DEPs) were identified in RH strain when comparing RH/PRU strains, in PRU strain when comparing PRU/PYS strains, and in PYS strain when comparing PYS/RH strains, respectively. Functional characterization of the DEPs using GO, KEGG, and STRING analyses revealed marked differences between the three strains. In silico kinase substrate motif analysis of the DEPs revealed three (RxxS, SxxE, and SxxxE), three (RxxS, SxxE, and SP), and five (SxxE, SP, SxE, LxRxxS, and RxxS) motifs in RH strain when comparing RH/PRU strains, in PRU strain when comparing PRU/PYS, and in PYS strain when comparing PYS/RH strains, respectively. This suggests that multiple overrepresented protein kinases including PKA, PKG, CKII, IKK, and MAPK could be involved in such a difference between T. gondii strains. Kinase associated network analysis showed that ROP5, ROP16, and cell-cycle-associated protein kinase CDK were the most connected kinase peptides. Our data reveal significant changes in the abundance of phosphoproteins between T. gondii genotypes, which explain some of the mechanisms that contribute to the virulence heterogeneity of this parasite.


INTRODUCTION
Toxoplasma gondii is a strictly intracellular protozoan parasite, which can infect almost any warm-blooded animal and over two billion people, worldwide. T. gondii infection is generally mild and asymptomatic in healthy immunocompetent individuals, but can be life-threatening in patients with compromised immune system, such as HIV-infected persons and organ transplant recipients, and in fetuses infected in utero (Elsheikha, 2008). The three clonal lineages of T. gondii commonly found in Eurasia and North America (Type I, II, and II) vary greatly in terms of virulence and phenotypic plasticity (Howe and Sibley, 1995;Sibley and Ajioka, 2008;Sibley et al., 2009a;Khan et al., 2011;Pittman et al., 2014;Shwab et al., 2014). The type I strains are highly virulent for mice, but with limited ability to transform to bradyzoites (Sibley and Boothroyd, 1992;Saeij et al., 2005). Type III strains can transform to bradyzoites and induce chronic infection, but they rarely cause clinical illness in humans. Type II strains are intermediate between type I and type III (Cheng et al., 2015;Ivanova et al., 2016). The Chinese 1 (ToxoDB#9) is the predominant lineage in mainland China and accounts for 78% of total strains isolated from animals and humans (Jiang et al., 2013;Wang et al., 2013Wang et al., , 2015Li et al., 2014).
Microarray and next generation DNA sequencing (NGS) technologies have been used to investigate genetic basis for phenotypic differences between different T. gondii type I strains (Yang et al., 2013). Also, two-dimensional difference gel electrophoresis (2D-DIGE) coupled with MALDI-TOF-MS was used to identify differentially expressed proteins among tachyzoites of different T. gondii genotypes (Zhou et al., 2014). The differential expression of proteins involved in the parasite virulence can even vary between virulent strains, such as Korean Isolate-1 (KI-1) and RH, of the same genotype (Choi et al., 2010). The role of protein phosphorylation in the mediation of virulence of T. gondii has been also reported (Joyce et al., 2011;Lorenzi et al., 2016). Protein phosphorylation is one of the most important post-translational protein modifications (PTMs) and functions in signal transduction pathways by modulating protein activity and protein-protein interactions (Pawson and Scott, 2005;Mithoe and Menke, 2011). Phosphorylation is regulated by protein kinases and protein phosphatases through their dynamically balanced actions. T. gondii genome consists of 8,311 ORFs, including 159 predicted eukaryotic-like protein kinases and harbors 108 active eukaryotic-like protein kinases (Peixoto et al., 2010;Lorenzi et al., 2016). Difference in virulence between T. gondii strains has been attributed to polymorphic variations in rhoptry (ROP) protein kinases and pseudokinases secreted by T. gondii strains to mediate distinct pathogenic processes in the host cell, such as virulence and modulation of host cell signaling (Saeij et al., 2006;Taylor et al., 2006;Khan et al., 2009;Sibley et al., 2009b;Behnke et al., 2011;Reese et al., 2011;Lim et al., 2012). A previous phosphoproteomic study reported that phosphorylation can regulate proteins that play critical role in the invasion and remodeling of the host cell and that T. gondii calcium-dependent protein kinase 1 (CDPK1) contains at least three major phosphosites (serine 25, 61, and 349) (Treeck et al., 2011). Analysis of kinase-deficient mutants T. gondii strains showed that protein phosphorylation plays an important role in the phosphorylation and activation of the host transcription factors STAT3 and STAT6 (Saeij et al., 2007;Ong et al., 2010), TgCDPK3-mediated parasite egress (Treeck et al., 2014), and efficient production of brain cysts (Sugi et al., 2017). Additionally, the rhoptry proteins ROP18 and ROP5 were found to mediate T. gondii evasion of the murine interferongamma response (Niedelman et al., 2012). The virulence of T. gondii Type I strains in mice is attributed to inactivation of host IRG resistance proteins by targeted phosphorylation of threonine residues in the switch I loop by ROP18 kinase (Steinfeldt et al., 2015). Therefore, studying the differences in the abundance of phosphoproteins between major genotypes of T. gondii can complement the picture derived from previous transcriptomic and proteomic profiling studies and provide new insight into the molecular mechanisms controlling parasite virulence.
Recent advances in quantitative mass spectrometrybased proteomics approaches, such as iTRAQ (Isobaric tag for relative and absolute quantitation) combined with phosphopeptide enrichment using titanium dioxide (TiO 2 ) affinity chromatography have enabled the sensitive and accurate identification and characterization of the proteome and phosphoproteome of multiple samples simultaneously and in an unbiased manner (Becker and Bern, 2011;Borchert et al., 2012;Glibert et al., 2015;Solari et al., 2016;Amorim et al., 2017). In the present study, a quantitative phosphoproteomic approach using iTRAQ combined with TiO 2 affinity chromatography was employed to characterize the phosphoproteomic differences between T. gondii tachyzoites of the three major genotypes. Our findings contribute a new resource for understanding strain-specific variations in molecular signaling pathways, which may underpin some of the mechanisms that contribute to the pathogenicity of T. gondii.

Parasite Preparation
Tachyzoites of Type I (RH) and Chinese 1 (PYS) T. gondii were recovered from frozen samples stored in liquid nitrogen as described previously (Zhou et al., 2014). Tachyzoites of each strain were inoculated into female, 6-8 weeks old, specificpathogen-free (SPF) BALB/c mice in order to allow tachyzoites to revive. Mice were provided by Laboratory Animal Center of Lanzhou Veterinary Research Institute. SPF BALB/c mice were infected with tachyzoites (10 mice per strain). Mice were monitored twice daily and euthanized when they exhibited clinical signs indicative of T. gondii infection, such as reduced appetite, ruffled fur and head tilting. The body cavity of mice was flushed with sterile phosphate-buffered saline (PBS) and the peritoneal liquid containing the tachyzoites was pelleted by centrifugation at 1,680 × g for 15 min. The pellet was washed three times with PBS and the supernatant was discarded after the final wash. Finally, the parasite pellet was resuspended in 1 ml PBS and kept in an Eppendorf tube at −80 • C, until use. For collection of tachyzoites of type II PRU T. gondii, SPF BALB/c mice were treated with 0.2 mg of dexamethasone (DSMS) in drinking water on alternate days for three times. Then, 10 mice were infected orally with ∼150 cysts of T. gondii type II PRU strain. Administration of DSMS on alternated days was continued after infection until mice exhibited clinical signs indicative of T. gondii infection at 9 days after infection, then mice were sacrificed and tachyzoites of the PRU strain were harvested from the peritoneal cavity of mice as stated above.

Protein Extraction and Digestion
Protein extraction was obtained from two biological repeats from T. gondii tachyzoites of each of the Type I (RH), Type II (PRU) and Chinese 1 (PYS) strains. Briefly, ∼5 × 10 8 tachyzoites of each strain were lysed in a protein lysis buffer (7 M urea/2 M thiourea, 0.2% SDS, 20 mM Tris-HCl, 10 mM DTT, pH 8.5) containing PMSF (1 mM, Thermo Scientific) and phosphatase inhibitors (PhosSTOP, Roche) by sonication (2 s ON/3 s OFF pulses for 5 min on ice). Lysate debris were eliminated by centrifugation at 25,000 g for 20 min at 4 • C. After centrifugation, samples were deoxidized using 10 mM dithiothreitol (DTT, 60 min, 56 • C) and alkylated by 55 mM IAM (45 min at ambient condition). Subsequently, the concentration of protein was determined using the Bradford assay. Protein of each sample (300 µg) were digested by trypsin (Promega) with an enzyme-to-protein ratio of 1:40 overnight at 37 • C. Samples were treated with 0.5% (v/v) formic acid (FA) to terminate the enzymatic reaction and then were desalted utilizing Strata X solid-phase extraction columns (Phenomenex). Finally, samples were dried using a SpeedVac and used for iTRAQ labeling.

iTRAQ Labeling
Samples were labeled with 8-plexed iTRAQ reagent (Applied Biosystems/MDS Sciex, Foster City, CA) according to the manufacturer's recommendations. The label reagents were blended with 50 µl isopropyl alcohol and peptides were resuspended in 30 µl dissolution buffer. Then, label reagents were added into peptides and incubated for 2 h at ambient condition. Labeled peptides were desalted and concentrated on Sep-Pack C 18 Cartridges (Waters). Proteins of T. gondii strains were labeled with iTRAQ tags as follows: Type I RH strain (113 and 114), Chinese 1 PYS strain (115 and 116), and Type II PRU strain (117 and 118).

Enrichment of Phosphopeptides Using TiO 2 Beads
The labeled phosphopeptides were enriched by Titanium dioxide (TiO 2 , GL Sciences, Tokyo, Japan). The TiO 2 beads were suspended in 1 ml loading buffer (2% glutamic acid/65% ACN/2% TFA/) for 10 min. Then, iTRAQ-labeled peptides were added into the freshly prepared TiO 2 beads at peptides-tobeads ratio of 1:4 (mass/mass). The peptide-beads slurry was shaken on a rotator at 37 • C. After 1 h, the blend was pelleted at 12,000 g for 5 min, and the supernatant was discarded. The pellet was rinsed three times in 2 ml loading buffer by centrifugation at 12,000 g for 5 min. The pellet of the final wash was suspended in 600 µl elution buffer and were shaken for 20 min. Then, the phosphopeptides combined to TiO 2 were eluted with 500 µl elution buffer, vacuum-dried and subjected to LC-MS/MS analysis.

LC-MS/MS
LC-MS/MS analysis was performed with the Q Exactive Mass Spectrometer (Thermal Scientific) combined with HPLC (Phenomenex columns Gemini-NX 3u C18 110A 150 * 2.00 mm). The fractionation of the samples was performed using a High pH Reversed-Phase Peptide Fractionation Kit (Thermo Scientific Pierce) according to the manufacturer's instructions. Briefly, the enriched peptides were redissolved with 300 µl 0.1% TFA and then were loaded into the equilibrated fractionation spin column. Peptides were bound to the resin and desalted by washing the column with water using low-speed centrifugation. The volatile high-pH elution with a step gradient of increasing acetonitrile concentrations were used to elute bound peptides from the spin columns into six different fractions. Each fraction was collected and vacuum-dried by centrifugation (Thermo Scientific). Then, samples were dissolved with 20 µl buffer A (2% ACN, 0.1% FA) and were loaded into a C 18 trap column with a speed of 8 µl/min for 4 min. The linear gradient of buffer B (95% ACN, 0.1% FA) running at 300 nl/min was applied to elute the samples. The elution started from 5% buffer B and increased into 35% buffer B. Then, the concentration of buffer B gradually rose into 60% within 5 min and the linearizing gradient of buffer B ascending into 80% for 2 min was followed. Furthermore, the concentration of buffer B was held at 80% for 2 min and was returned into 5% within 1 min, followed by a maintenance of 5% buffer B for 10 min. Finally, the spectra of first-grade MS were acquired during the scan scope of 350-1,800 m/z in the resolution of 70,000. The AGC (automatic gain control) target and maximum injection time of the first-grade MS were 3e6 and 40 ms, respectively. The spectra of second-grade MS were obtained using the following parameter: resolution = 17,500; maximum injection time = 40 ms; AGC target = 1e5. The mode of second-grade MS spectra was HCD (high-energy collisional dissociation) and the normalized collision energy of which was set as 30 eV.

Bioinformatics
To predict the probable kinase phosphorylation motifs of phosphopeptides identified in our study, the bias of amino acid residues approximating phosphorylation sites was investigated using the online software Motif-X (http://motif-x.med. harvard.edu/). After removal of uncertain phosphorylation sites, phosphopeptides with the length of ±7 aa from the phosphorylation site position were submitted to the Motif-X software. Each motif detected in the alignment was graphically displayed with the logo-like representation. Gene Ontology (GO) analysis (http://www.geneontology.org) was performed to identify functional categories (molecular function, biological process and cell component) of the DEPs. KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis (http://www. genome.jp/kegg/) was carried out to annotate identified DEPs to KEGG pathways. STRING (Search Tool for the Retrieval of Interacting Genes/Protein, http://string-db.org/) database was searched to identify processes that are regulated by phosphorylation and to analyze interaction networks of the DEPs involved. Interactors with a high confidence (combined score ≥ 0.7) were presented in protein-protein interaction (PPI) networks visualized by Cytoscape (version 3.0.2).

Correlation of T. gondii Kinase Peptides With Phosphopeptides
We used HMMER (http://hmmer.org) to search thePfam 31.0 database to identify sequences of the phosphorylated proteins with Pfam matches (Finn et al., 2016). We obtained a list of phosphorylated peptides related to protein kinase and protein tyrosine kinase domains. Then, we performed a Pearson correlational analysis to associate each of the kinase-related phosphorylated peptide with all phosphorylated peptides in the dataset. This analysis was performed using the DEPs ratio values per each T. gondii strain. Only correlations >0.998 were considered. Nodes with log 2 fold change values −0.75 ≥ log 2 RS1 ≤ 1.50, −0.75 ≥ log 2 RS2 ≤ 1.00, and −1.00 ≥ log 2 RS3 ≤ 1.00 were excluded. This network was visualized using Cytoscape 3.6.0 (https://cytoscape.org).

Prediction of Phospho-Enabled and Phospho-Disabled Protein Interactions
To identify potential interactions that are enabled or disabled by phosphorylation of the specific sites we used Mechismo (Betts et al., 2015). Because Mechismo does not support T. gondii proteomes, we used the eggNOG-mapper ([PMC4702882] http://eggnogdb.embl.de/#/app/emapper) to retrieve the human and yeast orthologs of the proteins identified by our phosphoproteomic analysis of T. gondii strains. We first extracted the relevant Orthologous Groups (OG), and then searched for the yeast and human orthologs within these groups using BioPython tools (http://biopython.org/). We then used Clustal Omega (Sievers et al., 2011) to align each orthologous sequence to the corresponding T. gondii protein, and identified the corresponding phosphorylation sites in the human or yeast proteins. These were used as input to Mechismo (http:// mechismo.russelllab.org/) to identify orthologous that enabled or disabled interactions in the human and yeast network. These were then mapped back to T. gondii proteome.

Identification and Quantification of Phosphopeptides
Total proteins obtained from T. gondii tachyzoites of RH, PRU, and PYS strains were digested and phosphopeptides were enriched by TiO 2 (titanium dioxide). The enriched FIGURE 3 | Volcano plots showing the median phosphopeptides log 2 fold change plotted against the -log 10 P-value. The results highlight the phosphopeptides up/down-regulated in RH strain when comparing RH/PRU strains (A), in PRU when comparing PRU/PYS strains (B), and in PYS strain when comparing PYS/RH strains (C), respectively. Red and green circles represent upregulated and downregulated DEPs, respectively. Dotted lines indicate (±) 1.5-fold change.

Prediction of p-Site Motifs Using Motif-x
To evaluate the conservation of sequence surrounding phosphorylation sites, motif analysis was performed between −7 and +7 positions of phosphorylation sites in all identified phosphoproteins and DEPs of RH vs. PRU, PRU vs. PYS, PYS vs. RH to detect the frequencies of amino acid residue at specific positions and identify conserved motifs. This analysis provided insights into classes of kinases differentially activated between tachyzoites of the three T. gondii strains. This analysis discovered 17 motifs in all phosphorylated proteins ( Table 1) and 11 of these represent motifs of the DEPs between parasite strains ( Figure 5). Kinases corresponding to each motif are shown in Table 2. Comparing RH to PRU, three motifs (RxxS, SxxE, and SxxxE) were identified in the up-regulated phosphorylation sites and no motif was detected in the down-regulated phosphorylation sites ( Figure 5A). One RxxS motif was found among the up-regulated phosphorylation sites and two motifs (SxxE, SP) were detected among the down-regulated phosphorylation sites in PRU strain when comparing PRU/PYS strains ( Figure 5B). For PYS vs. RH, three motifs (SxxE, SP, and SxE) had higher expression level in PYS strain and two motifs (LxRxxS and RxxS) had higher expression level in RH strain ( Figure 5C).

Functional Classification of DEPs Using GO Analysis
To gain a better understanding of the functional extent of the differences in protein phosphorylation between strains and  its effects on cellular processes, we performed GO enrichment classification of the differentially expressed phosphorylated proteins (Figures 6-8). Comparing RH to PRU, 86, 185, and 30 GO terms were assigned to molecular function (MF), cellular component (CC), and biological process (BP), respectively. The classification of the MF revealed a significant enrichment of both upregulated and downregulated phosphoproteins involved in binding, heterocyclic compound binding and organic cyclic compound binding. In the CC category, significantly enriched upregulated and downregulated phosphoproteins were involved in membrane, membrane part, intrinsic component of membrane, and integral component of membrane. Interestingly, in the BP category, enriched upregulated phosphoproteins were involved in metabolic process, cellular process and organic cyclic compound metabolic process, whereas enriched downregulated phosphoproteins were involved in cellular process, biological regulation, and regulation of cellular process (Figure 6). In regards to PRU vs. PYS, 69 MF, 130 BP, and 26 CC GO terms were identified. Binding, heterocyclic compound binding and organic cyclic compound binding were the top three MF GO terms, which included most enriched upregulated and downregulated phosphoproteins. In the CC category, the GO term with most enriched downregulated phosphoproteins were cell, cell part, intracellular and intracellular part, whereas the GO term with most enriched upregulated phosphoproteins included membrane, membrane part, intrinsic component of membrane and integral component of membrane. In the BP category, the top three GO terms with most enriched upregulated phosphoproteins were cellular process, metabolic process and organic substance metabolic process. In contract, GO terms cellular process, metabolic process and cellular metabolic process represented the most enriched downregulated phosphoproteins (Figure 7).
Comparing PYS to RH revealed 95 MF, 219 BP, and 46 CC GO terms. The most three abundant GO terms for upregulated or downregulated phosphoproteins under MF category were binding, followed by heterocyclic compound binding and organic cyclic compound binding. For the CC category, membrane, membrane part and intrinsic component of membrane were the mostly prevalent terms for up-and down-regulated phosphoproteins. In the BP category, the most prominent GO terms for upregulated and downregulated phosphoproteins were cellular process and establishment of localization, respectively (Figure 8).
We mapped the extent of association among GO functions by generating co-expression networks of GO terms related to the three functional categories using the R Igraph package. As shown in Figure S1, nine GO terms (heterocyclic compound binding, small molecular binding, organic cyclic compound binding, nucleoside phosphate binding, nucleotide binding, drug binding, ATP binding, adenyl nucleotide binding, and adenyl ribonucleotide binding) under MF category, two terms (membrane and membrane part) under CC category and three terms (localization, establishment of localization and transport) under BP category were significantly enriched in the GO network of upregulated phosphoproteins in RH strain when comparing RH/PRU strains. There were seven, two and three connections among the nine GO terms under MF, two GO terms under CC and three GO terms under BP, respectively. However, only three terms (cell communication, signal transduction and signaling) under BP and no term under MF or CC were significantly enriched in the GO network of downregulated phosphoproteins in RH strain when comparing RH/PRU strains ( Figure S2).
There was no significantly enriched GO term in any category of the GO network of upregulated phosphoproteins in PRU strain when comparing PRU/PYS strains ( Figure S3). Whereas, each category has significantly enriched GO terms in the GO network of downregulated phosphoproteins in PRU strain when comparing PRU/PYS strains ( Figure S4). These included: two GO terms (cation binding and metal ion binding) in the MF category, two terms (nuclear part and nucleoplasm part) in the CC category and 14 enriched terms under the BP category (cellular aromatic compound metabolic process, organic cyclic compound metabolic process, heterocycle metabolic process, nucleic acid metabolic process, nucleobase-containing compound metabolic process, RNA metabolic process, cellular metabolic process, cellular process, cellular macromolecule metabolic process, macromolecule metabolic process, metabolic process, nitrogen compound metabolic process, cellular nitrogen compound metabolic process, and gene expression).
In the GO network of upregulated phosphoproteins in PYS strain when comparing PYS/RH strains, three GO terms FIGURE 5 | Phosphorylation-specific motif analysis using the Motif-X algorithm. The representations of the serine phosphorylation motifs identified from all unambiguous phosphorylation sites are depicted as sequence logos. The height of the residues represents the frequency with which they occur at the respective positions. The color of the residues represents their physicochemical properties. Arrows pointing up or down on the right of each motif denotes that the motif was up-regulated or down-regulated, respectively. (antioxidant activity, peroxidase activity, and oxidoreductase activity, acting on peroxide as acceptor) under MF category were enriched. However, there was no significantly enriched GO term in the CC or the BP category of the GO network of upregulated phosphoproteins ( Figure S5). By contrast, the GO network of downregulated phosphoproteins showed two significantly enriched GO terms (heterocyclic compound binding and organic cyclic compound binding) in the MF category and four enriched GO terms (membrane, membrane part, intrinsic component of membrane and integral component of membrane) in the CC category, and 11 GO terms (localization, macromolecule localization, establishment of localization, protein localization, transport, nitrogen compound transport, establishment of protein localization, organic substance transport, protein transport, amide transport, and peptide transport) in the BP category ( Figure S6).

KEGG Pathway Analysis
A total of 147, 121, and 177 DEPs were annotated in the KEGG database and were mapped to 171, 153, 185 pathways in RH strain when comparing RH/PRU strains, in PRU strain when comparing PRU/PYS strains, and in PYS when comparing PYS/RH strains, respectively. The significantly enriched pathways (p-value ≤ 0.05) between the parasite strains are shown in Figure 9. In RH vs. PRU, 8, 7, and 5 phosphoproteins were significantly enriched in toxoplasmosis, glycerophospholipid metabolism, and glycerolipid metabolism pathways, respectively ( Figure 9A). In PRU vs. PYS, biosynthesis of amino acids, methane  metabolism, taste transduction, microbial metabolism in diverse environments, longevity regulating pathway-multiple species and biosynthesis of antibiotics were the five significantly enriched pathways, and included 6, 3, 3, 8, 5, and 8 phosphoproteins, respectively ( Figure 9B). There were two significantly enriched pathways in PYS vs. RH (viral carcinogenesis and calcium signaling pathway), and included 9 and 5 phosphoproteins, respectively ( Figure 9C).
To examine the association between enriched pathways, the pathway enrichment network was generated for RH vs. PRU, PRU vs. PYS, and PYS vs. RH. In the context of RH vs. PRU, the significantly enriched pathway "RNA transport" was the most independent, because it was only connected to the "mRNA surveillance pathway." Four pathways, including "toxoplasmosis, " "AMPK signaling pathway, " "insulin signaling pathway, " and "cGMP-PKG signaling pathway" were linked ( Figure S7). The significantly enriched pathway "glycerophospholipid metabolism" was connected with "glycerolipid metabolism, " "phosphonate and phosphinate metabolism, " and "ether lipid metabolism" pathways.

Protein-Protein Interactions (PPI)
PPI networks (combined score > 0.7) were constructed to understand the regulatory mechanisms of PTMs and to identify the relevant functional clusters of the DEPs between the parasite strains. The PPI network of up/down-regulated phosphoproteins in RH when comparing RH/PRU included 62 protein nodes and 96 interactor edges. Phosphoproteins involved in RNA splicing or translation, DNA replication or nucleosome remodeling, enzymes, molecular chaperones and RNA synthesis or translation formed major hubs (Figure 10). The PPI network of up/downregulated phosphoproteins in PRU when comparing PRU/PYS involved 59 nodes and 98 edges and contained RNA splicing or translation, DNA replication or nucleosome remodeling, enzymes and molecular chaperones clusters (Figure 11). There were 73 nodes and 133 edges in the PPI network of up/downregulated phosphoproteins in PYS strain when comparing PYS/PRU strains (Figure 12), which included RNA splicing or translation, DNA replication, enzymes, molecular chaperones, and translation functional clusters.

Kinase Associated Network
Based on the assumption that phosphorylation affects the kinase activity of the protein, we performed a correlation analysis to identify potential peptides whose phosphorylation was associated with that of the respective kinase (Table S4). We identified 11 phosphorylated kinases, which after keeping phosphopeptides with > ∼2-fold change abundance, were found to be correlated with 132 phosphopeptides. Positive correlation indicates that phosphorylation of the kinase has a positive effect on its function leading to downstream effects similar to the kinase phosphorylation status, whereas negative correlation implies the opposite effect. While these correlations cannot confirm a regulatory relationship, they may suggest potential kinases and associated substrates that are different among T. gondii strains. Peptides related to T. gondii rhoptry protein 16 A0A125YY08 (connected with 84 phosphorylated peptides), cell-cycle-associated protein kinase CDK S8EVL9 (connected with 84 phosphorylated peptides) and type II rhoptry protein 5 A0A125YQ30 (connected with 89 phosphorylated peptides) were the most connected kinase peptides. Kinases S8G171 and S8GM6 formed a second cluster, and the third cluster of phosphorylated peptides encompassed two proteins: cell-cycle-associated protein kinase CDK A0A125YUD2 and S8FAN7 ( Figure S10). The most connected kinases peptides (ROP5, ROP16 and cell-cycleassociated protein kinase CDK) were enriched in molecular function GO terms related to ATP binding, cyclin-dependent protein serine/threonine kinase activity, MAP kinase activity, and RNA polymerase II carboxy-terminal domain kinase activity.

Potential Phosphor-Dependent Interactions
We then sought to identify the putative role of the different phosphorylation sites in our dataset within the context of the T. gondii protein interaction network. After mapping our T. gondii proteins to their human orthologs and using Mechismo, which predicts whether a mutation or phosphorylation affects the interaction between two proteins by enabling or disabling it, we identified a total of 83 interacting partners of our T. gondii orthologs as enabled or disabled by phosphorylation of our identified sites (Table S5). Among these, 16 could be mapped back to T. gondii proteins and therefore were identified as putatively affected protein interactions in T. gondii (Figure S10). The S700 phosphosite of T. gondii S8EUC5 gene affected 47 interactions and phosphosite T1189 of T. gondii S8EVL9 gene disabled 8 interactions. On the other hand, the S420 phosphosite of T. gondii A0A125YSW4 gene disabled 6 interactions and the S355 phosphosite of S8GUL9 gene enabled 2 interactions and disabled 10 interactions.

DISCUSSION
Pathogenic T. gondii strains seem to rely on the ability to orchestrate the expression of virulence-associated genes (Yang et al., 2013) and proteins (Zhou et al., 2017) for successful infection. Many protein kinases have been also shown to influence the parasite virulence, suggesting a correlation between protein phosphorylation and T. gondii virulence. However, knowledge of the phosphoproteomic profile of strains of different virulence capabilities remains limited. Here, we performed a global phosphoproteomic profiling of T. gondii strains (type I [RH strain], type II [PRU strain], and Chinese ToxoDB#9 [PYS strain]) using iTRAQ-based MS approach combined with TiO 2 affinity chromatography. Our data showed that protein phosphorylation is a key regulatory PTM that plays many roles FIGURE 10 | Interaction network (combined score > 0.7) of the identified phosphoproteins up/down-regulated in RH strain when comparing RH/PRU strains. The nodes represent phosphoproteins and the edges denote interactors between phosphoproteins. The upregulation and downregulation of phosphoproteins were characterized by the color and size of nodes. Red and large size nodes denote up-regulated phosphoproteins in RH, whereas green and small size nodes indicate down-regulated phosphoproteins in RH. The color of edge represents the combined score of interactors.
in the pathobiology of T. gondii, at the individual protein level and at the level of signaling pathways as described below.
We identified 1,441 phosphopeptides, 1,250 phosphorylation sites and 759 phosphoproteins. The phosphorylation sites included 1,131 (90.48%) serine phosphorylation (pSer), 116 (9.28%) threonine phosphorylation (pThr), and 3 (0.24%) tyrosine phosphorylation (pTyr). These results seem consistent with previous phosphoproteomic studies in eukaryotes which identified serine (Ser), threonine (Thr), and tyrosine (Tyr) as the most common sites for phosphorylation (Pearlman et al., 2011) and reported the relative abundance of serine, threonine and tyrosine phosphorylation as 90, 10, and 0.05% (Olsen et al., 2006;Huttlin et al., 2010). The fact that phosphoserine represents the most abundant type of phosphorylation is anticipated because most T. gondii proteins, which play role in the parasite replication, invasion and virulence, contain serine kinase domains. For example, serine rhoptry kinases (ROP5, ROP16 and ROP18) contribute significantly to the differences observed in the virulence between different genotypes Behnke et al., 2011). Also, TgSR3 is a serine-arginine-rich alternativesplicing factor that modulates alternative splicing of over 1,000 T. gondii genes (Yeoh et al., 2015). T. gondii tachyzoites of the highly virulent strain seem to employ a number of proteases to modify the host cells, and to facilitate their dissemination inside the body of the host (Ramírez-Flores et al., 2019). Inhibition of T. gondii serine proteases involved in the processing of MIC proteins and digestion of host proteins suppressed the parasite development and disrupted its secretory pathway (Shaw et al., 2002). Although threonine phosphorylation is less abundant than serine phosphorylation in T. gondii, phosphorylation of both serine and threonine seems to be linked. The activity of type 2C serine-threonine phosphatase can significantly affect the growth of T. gondii in the host cells (Jan et al., 2009). ROP18, an important virulent factor in T. gondii, is also a highly polymorphic serine-threonine kinase (Du et al., 2014). Underrepresentation of tyrosine phosphorylation compared to serine and threonine phosphorylation has been reported earlier in apicomplexan parasites (Treeck et al., 2011). Despite its limited contribution, tyrosine phosphorylation can be involved in important signaling pathways, such as MAPK pathway (Lim and Pawson, 2010;Treeck et al., 2011).

Motifs
Assigning the sequence of phosphosite motifs to special enzyme families can assign phosphoproteins into kinases. Motif analysis identified three, three, and five motifs from the DEPs in RH strain when comparing RH/PRU strains, in PRU strain when comparing PRU/PYS strains, and in PYS strain when comparing PYS/RH strains, respectively (Figures 5A-C). In RH vs. PRU strain, three motifs (RxxS, SxxE, and SxxxE) were upregulated. In PRU vs. PYS strain, the motif (RxxS) was up-regulated and two motifs (SxxE and SP) were downregulated. In PYS vs. RH strain, three motifs (SxxE, SP, and SxE) in PYS strain and two motifs (LxRxxS and RxxS) in RH strain were up-regulated. Each motif is consistent with one or several category enzymes. For example, enzymes consistent with motif RxxS include PKA, PKB/AKT, PKC, PKG, CKI, CKII, IKK, CaM-II, and ATM; enzymes consistent with motif SxxE include PKA, PKG, CKI, CKII, IKK, and ATM; and enzymes consistent with motif SP included PKB, PKG, IKK, ATM, CDC2, CDK, and MAPK (Nakayasu et al., 2013). These data show the differences in the kinases expressed in the tachyzoites of the three T. gondii strains. Motifs are critical amino acid sequences that play roles in the recognition of kinase to its substrate (Schwartz and Gygi, 2005;Pi et al., 2016). Therefore, these strain-specific differences in motifs reflect differences in the substrate-recognition abilities between T. gondii strains of different genotypes.
In T. gondii, cAMP-dependent protein kinase A (PKA) and cGMP-dependent protein kinase G (PKG) play a role in the control of parasite egress from host cells (Jia et al., 2017). CK I and CK II are conserved serine/threonine protein kinases, and are involved in DNA repair and actin recruitment in eukaryotes, respectively. The function of CK I in the lifecycle of T. gondii has not been determined as yet, however, a casein kinase (CK) II-like activity and parasitic type 2C phosphatase were found to modulate the phosphorylation of Toxofilin protein, which plays a role in the actin sequestration and motility of T. gondii (Delorme et al., 2003;Wei et al., 2013). T. gondii IκB kinase (IKK) phosphorylates host IκBα localized at the PVM in infected cells and modulation of the NF-κB pathway in host cell (Molestina and Sinai, 2005). The upregulation of RxxS and SxxE motifs in RH strain, and SxxE and SP motifs in PYS strain, compared to PRU strain, suggest that the functions of PKA, PKG, CKII and IKK kinases representing these motifs (e.g., parasite egress, actin dynamic, and manipulation of NF-κB pathway in host cell) are probably more prominent in the virulent RH and PYS strains than in the avirulent PRU strain.
T. gondii TPK2, a homolog of the CDC2 cyclin-dependent kinase, plays a role in the regulation of the host cell cycle (Khan et al., 2002). In addition, cyclin-dependent kinase 7 (Cdk7) of T. gondii can mediate mRNA synthesis (Deshmukh et al., 2016). The activity of T. gondii cyclin-dependent kinase 9 is regulated by Cdk7 (Deshmukh et al., 2018). T. gondii has two MAPK isomers, TgMAPK1 and TgMAPK2, which play roles in cellular processes, especially in stress response (Wei et al., 2013). The ability of T. gondii to interfere with the cell cycle and adapt to stresses inside the host is crucial for pathogenesis. Given that motif SP was upregulated in PYS strain compared to other strains, it is reasonable to assume that intracellular parasite progression and stress tolerance mediated by CDC2-like TPK2, CDKs, and MAPKs are more prominent in PYS strain; however, this remains to be confirmed.

Functional Characterization of the Phosphorylated Proteins
Here, all identified phosphoproteins were assigned to GO terms in the MF, BP, and CC categories. The top GO terms in the MF were involved in binding and catalytic activities (Figures 6-8), indicating the critical role of phosphorylation in the parasite motility and host cell invasion. Among all phosphorylated proteins of T. gondii strains, the largest percentage of CC terms were involved in the cell, cell part, intracellular, intracellular part, and membrane. Other enriched phosphoproteins were related to intracellular membrane bound organelle, and macromolecular complex and other cellular components, suggesting that T. gondii strains are well-equipped to deal with long-term interaction of the surrogate host cell. Notably, we detected a significant enrichment of phosphoproteins in the membrane even without using any enrichment or isolation procedures for this organelle. The signal transduction originates at the membrane, where phosphorylation of signaling proteins is a key step to mount an effective response to a perceived stress and to transmit a certain message. The highly enriched GO terms in the BP were involved in cellular process, primary metabolic process, organic substance metabolic process, and cellular metabolic process, highlighting the correlation between bioenergetics and parasite virulence. The fact that binding and metabolic processes are the most highly enriched GO terms agrees with the previous finding that T. gondii aldolase (Starnes et al., 2009) and fructose-1,6bisphosphate aldolase-adhesin (Shen and Sibley, 2014) play a role in energy production and host cell invasion.
KEGG pathway analysis showed that DEPs were annotated into pathways with direct relevance to T. gondii high metabolic and bioenergetics demands required to sustain the parasite growth and replication, such as toxoplasmosis, glycerophospholipid metabolism and glycerolipid metabolism, biosynthesis of amino acids, methane metabolism, microbial metabolism in diverse environments, longevity regulating pathway-multiple species and biosynthesis of antibiotics. However, some DEPs (PKA, CAP-Rf, HAUSP, SKIP, and M2-PK) were annotated into pathways that might seem unrelated to T. gondii, such as "viral carcinogenesis" pathway. However, PKA and CAP-Rf was found to participate in the pathway of hepatocellular carcinoma caused by hepatitis C virus, and HAUSP and SKIP were involved in the pathway of lymphoma by Epstein-Barr virus (Saha and Robertson, 2011;Vescovo et al., 2016). Also, M2-PK was involved in pathway of cancer caused by human papillomavirus (He et al., 2003). Given the similarity between T. gondii and cancer (Lun et al., 2015) and the growing evidence for the correlation between T. gondii infection and cancer (Thirugnanam et al., 2013), the role of these DEPs in T. gondii virulence would be interesting to unravel.
Results of PPI network analysis showed that the upregulated phosphoproteins involved in RNA splicing or translation, DNA replication or nucleosome remodeling, molecular chaperones and RNA synthesis or translation in the virulent RH strain, and the upregulated phosphoproteins involved in RNA splicing or translation, molecular chaperones in the virulent PYS strain seem to be overrepresented compared to the avirulent PRU strain. Several epigenetic control mechanisms are responsible for a range of biological processes in the tachyzoite and bradyzoite stages of the life-cycle, from the parasite growth to phenotypic plasticity. For example, TgCARM1, TgGCN5-B, and TgMYST-A and -B are essential histone modifying enzymes that contribute to the propagation of tachyzoites (Dixon et al., 2010). TgSET8 and the phosphorylation of histone H3 on Serine-10 have been implicated in cell cycle regulation of T. gondii. KMTox and TgMYST-B have been linked to the parasite's response to reactive oxygen species and DNA damage, respectively. TgSRCAP and TgGCN5-A have been implicated in the differentiation of tachyzoites into bradyzoites. In a sharp contrast, phosphoproteins involved in only enzymatic processes seem to be more abundant in avirulent PRU strain compared to the virulent PYS strain. Parasite enzymes involved in dopamine production (Prandovszky et al., 2011), metabolism of oxygen radicals, glycolysis and DNA repair are upregulated in bradyzoites to support the parasite dormancy (reviewed in Sullivan and Jeffers, 2012). Also, the eIF2 kinases in T. gondii phosphorylate TgIF2a to elicit translation control in order to endure stresses that the parasite encounters during its life cycle (Joyce et al., 2011). Additionally, T. gondii fructose 1,6-bisphosphatase enzyme is essential in glucose metabolism and plays a role in parasite replication in vitro and virulence in vivo (Blume et al., 2015). These results demonstrate the correlation between the protein phosphorylation and the potential pathogenicity or latency of T. gondii strains.

Kinase-Protein Interactions
Protein kinases phosphorylate specific substrates to control certain cellular processes. We predicted the potential kinaseprotein interactions and identified specific parasite peptides related to the phosphorylated peptides. The parasite kinases (ROP16, ROP5 and cell-cycle-associated protein kinase CDK) had the greatest connected peptides and had regulatory roles in ATP binding, MAP kinase activity, RNA polymerase II carboxyterminal domain kinase activity, and cyclin-dependent protein serine/threonine kinase activity. These findings show the putative role of phosphoproteins in orchestrating T. gondii invasion and metabolic regulation.

CONCLUSIONS
This study profiles the global phosphoproteome of T. gondii strains belonging to three different genotypes. Combined iTRAQ and TiO 2 affinity chromatography analysis revealed marked differences in the phosphorylation patterns of tachyzoites between three T. gondii strains, which could be contributing to the differences in their virulence potential. The role of phosphorylation in the alteration of parasite kinases was evident from the kinase associated network analysis, which showed ROP5, ROP16, and cell-cycle-associated protein kinase CDK to be the most connected kinases peptides. Differentially abundant phosphoproteins (DEPs) were significantly enriched in several pathways related to the pathogenesis of T. gondii. The relevance of the identified DEPs should be characterized in subsequent studies using targeted quantitative phosphoproteomics tools. Because different strains of T. gondii have significantly varied pathogenicity, it will be of interest in the future to interrogate other T. gondii strains and genotypes, which could help to explain inter-strain and inter-genotype differences in virulence.

DATA AVAILABILITY
All the mass spectrometry data have been submitted to the ProteomeXchange Consortium with identifier PXD007777.

ETHICS STATEMENT
All animal experiments were conducted with the approval of the Animal Ethics and Administration Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, and in strict accordance with the Animal Ethics Procedures and Guidelines of the People's Republic of China. Every effort was made to reduce the suffering of the animals used in the experiment.

AUTHOR CONTRIBUTIONS
X-QZ and HE conceived and designed the experiments. Z-XW and C-XZ performed the experiments, analyzed the data, and wrote the paper. GC-M, EP, J-JH, and H-YS contributed the reagents materials analysis tools. HE and X-QZ critically revised the manuscript. All authors read and approved the final version of the manuscript.

ACKNOWLEDGMENTS
We thank BGI-Shenzhen for technical assistance with the LC-MS/MS analysis.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb. 2019.00307/full#supplementary-material Figure S1 | GO function network of upregulated phosphoproteins in RH strain when comparing RH/PRU strains. The clusters coincide with the GO functional categories of the DEPs and are color-coded as indicated. Red, green and blue clusters denote GO terms related to biological process, cellular component and molecular function, respectively. The dark colored nodes represent significantly enriched GO terms (p-value ≤ 0.05). The solid and dashed arrows between nodes represent direct and indirect associations between GO terms, respectively. Figure S2 | GO function network of downregulated phosphoproteins in RH strain when comparing RH/PRU strains. The clusters coincide with the GO functional categories of the DEPs and are color-coded as indicated. Red, green and blue clusters denote GO terms related to biological process, cellular component and molecular function, respectively. The dark colored nodes represent significantly enriched GO terms (p-value ≤ 0.05). The solid and dashed arrows between nodes represent direct and indirect associations between GO terms, respectively. Figure S3 | GO function network of upregulated phosphoproteins in PRU strain when comparing PRU/PYS strains. The clusters coincide with the GO functional categories of the DEPs and are color-coded as indicated. Red, green, and blue clusters denote GO terms related to biological process, cellular component and molecular function, respectively. The solid and dashed arrows between nodes represent direct and indirect associations between GO terms, respectively. Figure S4 | GO function network of downregulated phosphoproteins in PRU strain when comparing PRU/PYS strains. The clusters coincide with the GO functional categories of the DEPs and are color-coded as indicated. Red, green, and blue clusters denote GO terms related to biological process, cellular component and molecular function, respectively. The dark colored nodes represent significantly enriched GO terms (p-value ≤ 0.05). The solid and dashed arrows between nodes represent direct and indirect associations between GO terms, respectively. Figure S5 | GO-network of upregulated phosphoproteins in PYS strain when comparing PYS/RH strains. Red nodes, green nodes, and blue nodes denote GO terms under biological process, GO terms under cellular component and GO terms under molecular function, respectively. Nodes with deeper red, deeper green, and deeper blue represent significantly enriched GO terms (p-value ≤ 0.05) under biological process, cellular component and molecular function, respectively. Edges (solid and dotted arrows between nodes) represent direct (solid lines) and indirect (dashed lines) interactions between molecules. Figure S6 | GO function network of downregulated phosphoproteins in PYS strain when comparing PYS/RH strains. The clusters coincide with the GO functional categories of the DEPs and are color-coded as indicated. Red, green, and blue clusters denote GO terms related to biological process, cellular component and molecular function, respectively. The dark colored nodes represent significantly enriched GO terms (p-value ≤ 0.05). The solid and dashed arrows between nodes represent direct and indirect associations between GO terms, respectively.   Figure S10 | Network of kinase-phosphopeptide associations and potential phosphor-dependent interactions identified by correlation analysis. Log 2 (fold changes) of phosphoproteins up/down-regulated in RH when comparing RH/PRU, PRU when comparing PRU/PYS and PYS when comparing PYS/RH were mapped in the nodes in blue, red, and green, respectively. The gray arrow indicates protein interactions that are predicted to be disabled by the phosphorylation and green arrows indicate protein interactions that are predicted to be enabled by the phosphorylation.
Table S1 | List of phosphoproteins up/down-regulated in RH strain when comparing RH/PRU strains (|log1.5-fold change| > 1 and p-value < 0.05). Table S2 | List of phosphoproteins up/down-regulated in PRU strain when comparing PRU/PYS strains (|log1.5-fold change| >1 and p-value < 0.05). Table S3 | List of phosphoproteins up/down-regulated in PYS strain when comparing PYS/RH strains (|log1.5-fold change| >1 and p-value < 0.05). Table S4 | Correlation values between kinase-related phosphopeptides and associated phosphopeptides. Table S5 | Interactions of human orthologous proteins affected by phosphorylation detected by Mechismo (http://mechismo.russelllab.org/). The interaction effects involved four categories: enabling, enablingWeak, disabling, and disablingWeak and Interaction Confidences involved three categories: High, Medium, and Low. T. gondii proteins and phosphorylation sites measured in our experiment were mapped to their human orthologs and Mechismo was used to predict whether these phosphorylations enable or disable interactions in the human network.