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₂ 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.

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.

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.

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₄CO₃ 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.

The tryptic peptides were fractionated by high pH reverse-phase 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 (Liu et al., 2016).

The tryptic peptides dissolved in NETN buffer (100 mM NaCl, 1 mM EDTA, 50 mM Tris–HCl, 0.5% NP-40, pH 8.0) were incubated with pre-washed anti-malonyl-K antibody beads (PTM Biolabs, Hangzhou, China) at 4°C overnight with gentle shaking (Zhang et al., 2016). Then, the beads were washed four times with NETN buffer and twice with H₂O. The bound peptides were eluted from the beads with 0.1% trifluoroacetic acid (TFA; Sigma-Aldrich, Saint Louis, United States) and were cleaned using C18 ZipTips (Merck Millipore, Billerica, United States).

The tryptic peptides were dissolved in 0.1% formic acid (solvent A), directly loaded onto a reverse-phase analytical column (15-cm 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 ultra-performance 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 data-dependent procedure was performed, selecting the 20 most intense ions for MS/MS with 30 s dynamic exclusion at a resolution of 17,500.

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.

Subcellular localization of all malonylated proteins identified in three separate experiments was predicted using WOLF PSORT¹. 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². 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)³, 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⁴ 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 Neighboring Genes (STRING) database⁵ 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).

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 experiments were performed, which identified 294 Kmal sites on 203 proteins in Exp.1 (Supplementary Table S1), 345 Kmal sites on 236 proteins in Exp.2 (Supplementary Table S2), and 352 Kmal sites on 230 proteins in Exp.3 (Supplementary 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.

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.

One motif, C (X₃) Kmal (Kmal indicates the lysine malonylation site, and X represents a random amino acid residue), was obtained (Figure 1D). Except C (X₃) Kmal, K (X₆) Kmal, Kmal (X₆) K, and K (X₅) Kmal were also detected in S. erythraea (Xu et al., 2016). However, the number of Kmal motifs was much lower than lysine acetylation motifs (Li et al., 2016; Lv et al., 2016; Zhang et al., 2016) and succinylation motifs (Pan et al., 2015; Jin and Wu, 2016). The sequence logos show a strong bias for cysteine (C), which was observed at the −3 and −4 positions (Figure 1C), similar to the Kmal bias for arginine (R) detected in the common wheat (Liu et al., 2018).

To understand the primary function of malonylated proteins, enriched GO terms at 2nd level in the three GO categories: biological process (BP), cellular component (CC) and molecular function (MF) were identified (Supplementary 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).

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-3-phosphate 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.

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.