Global Proteomic Analysis of Lysine Malonylation in Toxoplasma gondii

Lysine malonylation (Kmal) is a new post-translational modification (PTM), which has been reported in several prokaryotic and eukaryotic species. Although Kmal can regulate many and diverse biological processes in various organisms, knowledge about this important PTM in the apicomplexan parasite Toxoplasma gondii is limited. In this study, we performed the first global profiling of malonylated proteins in T. gondii tachyzoites using affinity enrichment and Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis. Three experiments performed in tandem revealed 294, 345, 352 Kmal sites on 203, 236, 230 malonylated proteins, respectively. Computational analysis showed the identified malonylated proteins to be localized in various subcellular compartments and involved in many cellular functions, particularly mitochondrial function. Additionally, one conserved Kmal motif with a strong bias for cysteine was detected. Taken together, these findings provide the first report of Kmal profile in T. gondii and should be an important resource for studying the physiological roles of Kmal in this parasite.


INTRODUCTION
Toxoplasmosis, caused by the protozoan parasite Toxoplasma gondii, is estimated to affect approximately one-third of the world population (Montoya and Liesenfeld, 2004). This parasite has the ability to infect almost all mammalian and avian species (Webster, 2010;Robert-Gangneux and Darde, 2012;Liu et al., 2015). The life cycle of T. gondii involves definitive host (members of the cat family; Felidae) and intermediate (mammals, including humans) host. During its development, the parasite progresses through three main morphological stages, including one replicative stage (tachyzoite), which is associated with the acute phase of infection, the dormant stage (bradyzoitescontaining tissue cyst), which is associated with latent form of infection, and the environmentally resistant oocyst stage. The life cycle of T. gondii includes asexual reproduction, which involves the formation of tachyzoites and bradyzoites-containing cysts in the intermediate host and sexual reproduction which involves the formation of oocysts in the feline intestinal epithelium. In order to adapt to different environments and survive inside various tissues within different hosts, the parasite tightly regulates its metabolic and protein functions at the post-translational level (Xiao et al., 2010;Dubey et al., 2017).
Lysine post-translational modifications (PTMs), such as acetylation (Choudhary et al., 2009), methylation (Peng et al., 2011), succinylation (Lin et al., 2012;Hirschey and Zhao, 2015), and ubiquitination (Hershko and Ciechanover, 1998), play key roles in broadening the functional diversity of proteins and impact significantly on the regulation of protein functions in prokaryotic and eukaryotic organisms (Lin et al., 2012;Hirschey and Zhao, 2015). Lysine succinylation has been investigated in T. gondii, where it was found to be involved in a broad range of cellular functions (Li et al., 2014). Also, lysine acetylation has been studied in three T. gondii strains belonging to three different genotypes and the level of acetylation was found to correlate with the parasite strain virulence and has been found widespread on proteins of diverse functions in T. gondii (Jeffers and Sullivan, 2012;Wang et al., 2019).
Post-translational modification of proteins via lysine malonylation (Kmal) has been reported across many metabolic pathways, such as fatty acid synthesis and oxidation (Hirschey and Zhao, 2015), mitochondrial respiration (Xie et al., 2012), glycolysis (Peng et al., 2011;Lin et al., 2012), and modification of histones (Hershko and Ciechanover, 1998). Lysine malonylation was firstly observed in mammalian cells and bacterial cells (Lin et al., 2012). Since then, there have been growing interests in exploring the regulatory roles of Kmal in various microbial species, such as Escherichia coli (Qian et al., 2016), Cyanobacteria (Ma et al., 2017), and Saccharopolyspora erythraea (Xu et al., 2016). Although Kmal can regulate many crucial and diverse cellular processes , its existence and function in T. gondii remain unknown.
In the present study, we characterized malonylated proteins in T. gondii tachyzoites using Liquid chromatography-tandem mass spectrometry (LC-MS/MS) coupled with sensitive immuneaffinity purification. Three parallel experiments were performed, which identified 294, 345, 352 Kmal sites on 203, 236, 230 proteins, respectively. Functional analyses showed predominant presence of malonylated proteins in metabolic processes, such as glycolysis/gluconeogenesis, aminoacyl-tRNA biosynthesis, pentose phosphate pathway, and fatty acid biosynthesis. To our knowledge, this study is the first to characterize protein malonylation in T. gondii. Our data lay the foundation for future investigations into the biological functions of malonylated proteins in T. gondii and in the context of hostparasite interaction.

Parasite Culture
Tachyzoites of T. gondii RH strain were maintained by serial passage in human foreskin fibroblast (HFF) monolayers, which were grown in Dulbecco's modified Eagle's medium (DMEM, Gibco, United States) supplemented with 10% fetal calf serum (FBS, Gibco, United States), 100 U/ml antibiotics (penicillin-streptomycin solution). T. gondii tachyzoites and HFF monolayers were cultivated in a 5% CO 2 humidified incubator at 37 • C. The tachyzoites and cell debris were harvested when the infected HFF monolayer was lysed. The mixture was washed several times with phosphate buffered saline (PBS) and passed through a 25 gauge needle. Then, the parasites were filtrated using a 3 µm membrane filters (Millipore) in order to remove the cell debris, and stored at −80 • C until use.

Protein Extraction
The parasite pellets were sonicated with 12 short bursts of 10 s, followed by intervals of 10 s on ice for three times in lysis buffer [8 M urea, 1% Protease inhibitor cocktail and for PTM experiments, deacetylase inhibitors were also added to the lysis buffer, e.g., 3 µM trichostatin A (TSA) and 50 mM nicotinamide (NAM)], and the remaining debris was removed by centrifugation at 12,000 g for 10 min at 4 • C. Finally, the pellet was discarded, and the protein concentration was examined with BCA kit.

Western Blotting Analysis
Proteins (20 µg) of tachyzoites were electrophoresed on 12% SDS-PAGE and transferred to PVDF membrane. TBST buffer (25 mM Tris-HCl, 150 mM NaCl, 0.1% Tween20, pH 8.0) with 5% BSA was used to block the membrane for 60 min. Then, the membrane was incubated with anti-malonyllysine antibody (1:500, catalog no. PTM-901; PTM Biolabs, Hangzhou, China) in TBST with 2.5% BSA overnight at 4 • C. The membrane was washed three times with TBST, followed by incubation with horseradish-peroxidase-conjugated Goat anti-Rabbit IgG (1:5000; Thermo) for 60 min at room temperature. After washing the membrane three times, an ECL kit was used for protein visualization.

Trypsin Digestion
The protein solution was reduced with 5 mM dithiothreitol (DTT) for 30 min at 56 • C and then alkylated with 11 mM iodoacetamide (IAA) for 15 min at 25 • C away from light. The proteins were diluted by adding 100 mM NH 4 CO 3 until the urea concentration became <2 M. Finally, trypsin was added at 1:50 trypsin-to-protein mass ratio for the first digestion overnight and 1:100 trypsin-to-protein mass ratio for a second 4 h-digestion.

HPLC Fractionation
The tryptic peptides were fractionated by high pH reversephase HPLC using Thermo Betasil C18 column (5 µm particles, 10 mm ID, 250 mm length). The peptides were separated using a gradient of 8-32% acetonitrile (pH 9.0) over 60 min into 60 fractions. Thereafter, the separated peptides were combined into four fractions as described previously .

Liquid Chromatography-Tandem Mass Spectrometry Analysis
The tryptic peptides were dissolved in 0.1% formic acid (solvent A), directly loaded onto a reverse-phase analytical column (15cm length, 75 µm i.d.) at a constant flow rate of 700 nL/min on an EASY-nLC 1000 UPLC system. The gradient included an increase from 6-23% solvent B (0.1% formic acid in 98% acetonitrile) over 26 min, followed by 8 min from 23-35%, and reaching to 80% in 3 min and holding at 80% for the last 3 min. The peptides were subjected to a nanospray ionization (NSI) source followed by tandem mass spectrometry (MS/MS) in Q Exactive TM Plus (Thermo) coupled online to an ultraperformance liquid chromatography (UPLC). The m/z scan range was 350-1800 for full scan, and intact peptides were detected in the Orbitrap at a resolution of 70,000. A datadependent procedure was performed, selecting the 20 most intense ions for MS/MS with 30 s dynamic exclusion at a resolution of 17,500.

Database Search
The MS/MS data were processed using MaxQuant search engine (v.1.5.2.8). Tandem mass spectra were searched against Uniprot T. gondii database (ToxoDB 46, 8322 sequences downloaded on March 18, 2020) concatenated with reverse decoy database. Trypsin/P was specified as cleavage enzyme allowing up to 4 missing cleavages. The mass tolerance for fragment ions was set as 0.02 Da. Cysteine carbamidomethylation was specified as fixed modification. Methionine oxidation, and Kmal were specified as variable modifications. False discovery rate (FDR) was adjusted to <1% and minimum score for modified peptides was set >40. All other parameters in MaxQuant were set to default values. The site localization probability was set as >0.75. Label-free quantification was performed using the MaxQuant label free quantification (LFQ) algorithm (Cox et al., 2014), by comparing the abundance of the same peptides across runs, with both ion intensities and spectral counts used for this purpose.

Bioinformatic Analysis
Subcellular localization of all malonylated proteins identified in three separate experiments was predicted using WOLF PSORT 1 . For Motif analysis, the sequence models contained amino acids at specific positions of the malonyl-21-mers (10 amino acids upstream and downstream of the Kmal sites) in all protein sequences were predicted by MoMo software (motif-x algorithm). The T. gondii protein sequences from the UniProt database was used as background database parameter, and other parameters were set as default.
GO annotation was performed against UniProt-GOA database 2 . Firstly, the identified protein identifiers (IDs) were converted to UniProt IDs and then mapped to GO IDs. If the identified malonylated proteins were not annotated by UniProt-GOA database, the InterProScan software were used to annotate protein's GO function based on protein sequence alignment method. Also, an online service tool, KEGG Automatic Annotation Server (KAAS) 3 , was used to annotate protein's KEGG database description. Then, the annotation result was mapped to the KEGG pathway database using KEGG mapper. The protein domains were annotated by using InterProScan software and the InterPro domain database 4 based on protein sequence alignment. A two-tailed Fisher's exact test was employed for the GO/KEGG/Domain enrichment analysis of the differentially malonylated proteins against all identified proteins. The terms with P values <0.05 were considered significant.
All differentially abundant malonylated proteins were searched against the Search Tool for Recurring Instances of 2 http://www.ebi.ac.uk/GOA/ 3 http://www.genome.jp/kegg/kaas/ 4 http://www.ebi.ac.uk/interpro/ Neighboring Genes (STRING) database 5 for Protein-Protein interaction (PPI). The interaction score was set with high confidence ( 0.7), and all other parameters were set as default. The PPI network was visualized by Cytoscape (version 3.5.0).

RESULTS AND DISCUSSION
Our study provides the first characterization of Kmal, a newly discovered type of lysine acylation, in T. gondii. We have identified multiple malonylated proteins with a wide range of functions, such as metabolism and immune response, as discussed below.

Global Detection of Kmal Sites on T. gondii Proteins
To detect the abundance of malonylated proteins in T. gondii tachyzoite extracts, western blotting was performed prior to the proteomic experiments ( Figure 1A). The global malonylome analysis of T. gondii was performed using affinity enrichment followed by high-resolution LC-MS/MS. Three parallel 5 https://string-db.org/   Table S3). A total of 160 proteins were overlapped between Exp.1 and Exp.2, 151 proteins were overlapped between Exp.1 and Exp.3, 170 proteins were overlapped between Exp.2 and Exp.3, and 138 proteins were overlapped among these three experiments. The low number of overlapped proteins may be attributed to the low incidence of this PTM and the resolution of the methods used, resulting in less repetitions. Overall, 506 Kmal sites across 326 proteins were detected. The number of Kmal proteins in T. gondii was higher than that reported in common wheat, S. erythraea (Xu et al., 2016;Liu et al., 2018), but was lower than that detected in E. coli (Qian et al., 2016), rhizobacterium Bacillus amyloliquefaciens FZB42 (Fan et al., 2017) and mammals (Colak et al., 2015), however was close to the number of malonylated proteins identified in Cyanobacteria (Ma et al., 2017). As shown in Figure 1B, a total of 326 proteins were modified/malonylated at more than two lysine residues. For example about 73% (239 out of 326) of the identified malonylated proteins were modified at one lysine residue, 14% (46 out of 326) were modified at two lysine residues, and 13% (41 out of 326) were modified at ≥ two lysine residues. The finding that 27% of the malonylated proteins were modified at multiple lysine residues/sites is in agreement with the results detected in wheat and S. erythraea (Xu et al., 2016;Liu et al., 2018), indicating some conservation in the level of protein malonylation between T. gondii and other organisms, such as bacteria and plants, however the level and site of malonylation may vary among different organisms.
Interestingly, T. gondii Elongation Factor 1-Alpha (TgEF-1α) was the largest malonylated protein with 10 malonylated sites, which seems to have functions related to cell motility, protein turnover, cell growth, and signal transduction (Ridgley et al., 1996), DNA replication and repair protein networks (Toueille et al., 2007), suggesting an important, conserved role for Kmal in the activities of this protein. The second highly malonylated protein was the eukaryotic mitochondrial porins which plays a role in transporting tRNAs and apoptosis (Flinner et al., 2013). T. gondii 14-3-3 protein had six Kmal sites, this protein is found in the parasitophorous vacuole and the excreted/secreted parasite antigens. T. gondii 14-3-3 protein has been suggested to play a role in the modulation of the migratory properties of host cells (Weidner et al., 2016) and host immunity (Assossou et al., 2004). In the present study, five Kmal sites were identified in heat shock protein 70 (Hsp70). A previous study showed that the heat shock protein was the most widely succinylated protein, which can be succinylated on up to 17 independent lysine residues (Li et al., 2014), and the heat shock protein can also be ubiquitinated (Salmon de Monerri et al., 2015). The heat shock protein was shown to play a critical role in modulating host immune response during Plasmodium falciparum infection (Pooe et al., 2017) and was also associated with protective immunity against T. gondii (Kikumura et al., 2010). These results suggest that Kmal participates in the regulation of various biological processes of different cellular components and influences several biological functions in T. gondii.

Motifs of Malonylated Peptides
To further examine the properties of Kmal sites in T. gondii, the flanking sequences from ten amino acids upstream to ten amino acids downstream of the modified site were studied using the MoMo software and hierarchical cluster analysis. As shown in Figure 1C, the frequency of cysteine (C) residue at the position −4 to −3 was highest, the positively charged lysine (K) residue was enriched at the positions −10 to −6 and +10 to +5, whereas the negatively charged residues glutamate (E) and phenylalanine (F) were not observed, this result agree with Kmal reported in common wheat (Liu et al., 2018). Previous study showed that G residues and E residues are highly enriched surrounding the acetylated lysine, while G is overrepresented in the positions on the N-terminal side (Jeffers and Sullivan, 2012). Therefore, the additional characteristics of the surrounding amino acids are not obvious from the motifs, except in some specific positions.

Primary GO Enrichment Analysis of the Malonylated Proteins
To understand the primary function of malonylated proteins, enriched GO terms at 2nd level in the three GO categories:  Table S4). In terms of BP, 57, 56, 56 and 56 malonylated proteins were enriched in cellular metabolic process, organic substance metabolic process, nitrogen compound metabolic process and primary metabolic process, accounting for 11% of all identified malonylated proteins in T. gondii, respectively (Figure 2A). Regarding the MF category, 16, 12, and 12% of the malonylated proteins were related to protein binding, organic cyclic compound binding and heterocyclic compound binding, respectively ( Figure 2B). Also, we used WoLF PSORT to further predict the subcellular locations of the identified proteins (Horton et al., 2007). As shown in Figure 2C, most of the malonylated proteins in T. gondii were located in the cytoplasm (30%) and nucleus (28%), the rest of malonylated proteins were predicted to be located in plasma membrane (15%) extracellular (11%) and mitochondria (10%) (Figure 2C).

Advanced Functional Enrichment of Malonylated Proteins
The GO enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed in T. gondii, as shown in Supplementary  Table S4. Based on the GO enrichment results of the CC category, malonylated proteins were significantly located in the cell body ( Figure 3A). In terms of MF, the most significantly enriched malonylated proteins were related to ligase activity and small molecule binding ( Figure 3B). Regarding BP, malonylated proteins were mainly implicated in purine metabolism and monocarboxylic acid biosynthetic  processes processes (Figure 3C). Interestingly, seven malonylated proteins in T. gondii belonged to the pentose phosphate pathway (Figure 4A). A highly proliferative organism such as T. gondii is an avid auxotrophic for many metabolites and must consume these, and carbon sources (e.g., glucose) from the host cell to sustain its growth and proliferation. The pentose phosphate pathway along with glycolysis and mitochondrial respiration are used by the intracellular tachyzoites of T. gondii to generate energy and essential anabolic precursors from glucose (MacRae et al., 2012). In addition, eight malonylated proteins in T. gondii belonged to the fatty acid biosynthesis (Figure 4B). These are involved in the growth, development, and reproduction of T. gondii (Fu et al., 2019). Malonyl-CoA can serve as a substrate of fatty acid synthesis and as an inhibitor of fatty acid oxidation (Foster, 2012). Thus, malonylation might indirectly regulate these important processes in T. gondii. KEGG pathway analysis showed that most malonylated proteins are enriched in glycolysis/gluconeogenesis, aminoacyl-tRNA biosynthesis, fatty acid biosynthesis ( Figure 5A). Furthermore, enriched protein domain analysis showed that malonylated proteins are significantly enriched in ribosomal protein L7Ae/L30e/S12e/Gadd45 family, glutamate/Leucine/ Phenylalanine/Valine dehydrogenase, phosphoglycerate kinase and pyridine nucleotide-disulfide oxidoreductase were significantly enriched in the malonylated proteins ( Figure 5B).
The malonylated protein fructose-bisphosphate aldolase, which plays a central role in glycolysis and gluconeogenesis pathways, has been considered as a potential target for drug development against pathogenic bacteria (Ziveri et al., 2017). Another malonylated protein (leucyl aminopeptidase; LAP) is a functional aminopeptidase in the cytoplasm of T. gondii (Jia et al., 2010). Lactate dehydrogenase (LDH1) was also identified as a malonylated protein, and has been previously shown to catalyze the conversion of pyruvate and lactate in anaerobic growth conditions and used for energy supply (Abdelbaset et al., 2017). T. gondii expresses two different lactate dehydrogenase enzyme genes, LDH1 is only expressed in the tachyzoite stage whereas LDH2 is preferentially expressed in the bradyzoite stage. Given the important role of this enzyme in the parasite virulence, phenotypic differentiation, and latent infection, the potential of using T. gondii LDH mutant strains as vaccine candidates has been suggested (Abdelbaset et al., 2017). Five Kmal sites were identified in glyceraldehyde-3phosphate dehydrogenase (GAPDH) in the present study, this protein is also succinylated and ubiquitinated in T. gondii, and the glyceraldehyde-3-phosphate dehydrogenase I (GAPDH I) was found in the parasite cytoplasm, succinylated on two lysines (Li et al., 2014;Salmon de Monerri et al., 2015). Previous reports showed that phosphorylation and glycation of GAPDH are involved in the regulation of (PPI) and intracellular localization of the enzyme, and that the aggregation of GAPDH can be affected by all types of PTMs (Muronetz et al., 2016). Lysine 213 of GAPDH can be malonylated in macrophages following lipopolysaccharide (LPS) stimulation, which causes its dissociation from several inflammatory mRNAs, promoting translation (Galván-Peña et al., 2019), implying that Kmal has important roles in these processes. A total of 83 malonylated proteins identified in the present study had also been identified as acetylated proteins in a previous study (Jeffers and Sullivan, 2012). KEGG pathway analysis showed that these 83 proteins were mainly enriched in glycolysis/gluconeogenesis, carbon fixation in photosynthetic organisms, aminoacyl-tRNA biosynthesis and biosynthesis of antibiotics (Figure 6), which indicate that malonylation and acetylation may have similarities or synergistic effects on some functions and mechanisms. Further studies are warranted to elucidate the relationship between malonylation and acetylation in T. gondii.

PPI Network of the Malonylated Proteins in T. gondii
The PPI analysis was performed by searching the STRING database and PPI networks were visualized using Cytoscape software in order to further identify the major biological processed affected by Kmal in T. gondii (Shannon et al., 2003; Figure 7, Supplementary Table S5). Using MCODE (Minimal Common Oncology Data Elements), a number of highly associated subnetworks of Kmal proteins were identified, including glycolysis/gluconeogenesis, ribosome, proteasome and aminoacyl-tRNA biosynthesis processes, which are in agreement with the KEGG pathway enrichment analysis, suggesting that these processes play important roles in shaping the proteomic landscape of T. gondii.

CONCLUSION
In this study, we described the first malonylome profile in T. gondii using mass spectrometry in combination with immune-affinity purification. Over 300 proteins were found to be malonylated and were localized in different cellular compartments. These malonylated proteins were associated with diverse biological and metabolic processes, suggesting that Kmal is involved in the regulation of T. gondii physiology. We also found that some important proteins related to T. gondii invasion, emergence and virulence were malonylated, which suggest that Kmal plays an important role in the regulation of functions of these critical proteins of T. gondii. Our findings will serve as a useful resource for further investigation of the functions of Kmal, especially in the context of parasite interaction with the host cells. Further studies should focus on the experimental verification of the functions of these malonylated proteins to reveal the exact roles of malonylation in T. gondii. Analysis of the extent of differences and similarities in Kmal between different developmental stages (bradyzoites, tachyzoites and oocysts) and between distinct genotypes of T. gondii is warranted, which in turn will facilitate the elucidation of mechanisms underlying the phenotypic transformation and varying virulence of T. gondii.

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

AUTHOR CONTRIBUTIONS
RD, F-CL, and X-QZ conceived and designed the research. L-BN performed the research, analyzed the data and drafted the manuscript. Q-LL, HE, and N-JH contributed materials, reagents, and analysis tools. HE, F-CL, RD, and X-QZ critically revised the manuscript. All authors reviewed and approved the final version of the manuscript.