Unravelling Strain-Specific Modifications of Toxoplasma gondii tRNA and sncRNA Using LC-MS/MS

ABSTRACT Many RNA modifications have been detected in rRNA, tRNA and small noncoding RNA (sncRNA) as well as in low-abundance RNA species such mRNA. Although RNA modifications play roles in many cellular and biological processes in various domains of life, knowledge about the diversity and role of RNA modifications in Toxoplasma gondii is limited. In this study, RNA modifications in three T. gondii strains (RH type I, PRU type II, and VEG type III) with distinct virulence abilities were determined by liquid chromatography-tandem mass spectrometry. We compared the levels of modifications of four nucleotides in tRNA and sncRNA, characterized RNA modification patterns of different T. gondii strains, and determined the diversity of RNA modifications. We detected and quantified 22 modified nucleosides in both tRNA and sncRNA. Significant differences in the diversity of the modified nucleosides were found between the three T. gondii strains. RNA modifications were correlated with the expression of many T. gondii virulence proteins. Some of the identified modifications (e.g., 2’-O-methylinosine, pseudouridine) play a role in mediating the host-parasite interaction. These results provide novel insight into the global modifications in tRNA and sncRNA, and the diversity of RNA modifications between T. gondii strains with different virulence backgrounds. IMPORTANCE Although RNA modifications play roles in many cellular and developmental processes in various domains of life, knowledge about the patterns and functions of RNA modifications in T. gondii is limited. Here, a quantitative liquid chromatography-tandem mass spectrometry (LC-MS/MS) approach was used to study global RNA modifications in T. gondii strains of distinct virulence backgrounds. We quantified 22 modified nucleosides in both tRNA and sncRNA. Significant T. gondii strain-specific differences in RNA modifications were detected. More tRNA modifications correlated with T. gondii virulence proteins than sncRNA modifications. RNA modifications were significantly correlated with virulence proteins. Our data provide the first comprehensive profiling of the modifications tRNA and sncRNA in T. gondii, expanding the diversity of RNA modifications in this parasite and suggesting new regulators for modulating its virulence.

T oxoplasma gondii is a widespread protozoan parasite that infects most warmblooded animals and roughly one-third of the world's human population (1). Cats and other felines serve as the definitive host by supporting the formation and excretion of oocysts, which represent the main source of T. gondii infection in humans and animals (2). T. gondii infection can also be acquired via ingestion of undercooked meat containing the parasite cysts (3). Infection in immunocompromised patients may lead to encephalopathy, coma, and death (4) and in pregnant women, can lead to abortion and congenital defects in the newly born infant (5). T. gondii can also cause abortion and infertility in farm animals (6). Given the significant impact of toxoplasmosis on the health of humans and animals, and the challenges to develop a vaccine for humans (7)(8)(9), enhanced understanding of the molecular basis of T. gondii pathogenicity may ultimately lead to the discovery of new treatment options.
RNA modifications and RNA modification enzymes (known as writers) are involved in many biological and pathological processes and play various roles across different domains of life (10)(11)(12). RNA modifications are dynamic and key players in the regulation of gene expression in response to an ever-changing cellular microenvironment. Three classes of molecules are involved in RNA modifications: "writers" (catalyze the deposition of specific modifications), "erasers" (catalyze the removal of specific modifications), and "readers" (for recognition and binding of the modified nucleotides) (13). Understanding RNA modifications in parasites, such as Plasmodium (14), Leishmania donovani (15), Trypanosoma brucei (16), and Echinococcus granulosa (17), has attracted recent attention due to the importance of the regulatory roles of RNA modifications in parasite physiology and development. RNA modifications contribute to stabilizing parasite RNA structures (15) and mediate pathogenicity and host responses to infection (18).
Currently, over 170 different types of posttranscriptional modifications have been identified in different types of RNAs and have energized the growing field of "epitranscriptomics" (19,20). The mRNA includes different epigenetic modifications which can regulate mRNA stability, splicing, and gene expression (19). The methylation of adenosine at position 6, known as N 6 -modification methyladenosine (m 6 A), is the most abundant modification in mRNA and long intergenic ncRNAs (lincRNAs) as well as found in primary miRNA and rRNA (21,22). m 6 A can regulate ncRNAs via processing and maturation of pri-miRNAs (23), stabilizing circRNAs (24), and regulating the interaction between lncRNAs and proteins (25). The role of m 6 A in T. gondii has recently been reported, where deletion of the m 6 A writer components METTL3 and WTAP causes arrest of parasite replication (26).
T. gondii is a genetically diverse organism, encompassing three main genotypes, namely, type I (e.g., RH strain), type II (e.g., PRU strain and ME49 strain), and type III (e.g., VEG strain and CTG strain), whose virulence varies significantly in mice (27). A previous study showed that avirulent isolates of T. gondii exhibit a higher rate of tRNA cleavage than virulent strains (28). Thus, unravelling the differences in RNA modifications between T. gondii strains which differ significantly with respect to their virulence is of particular importance as it may provide a key to hitherto undiscovered RNA modifications that may underpin the regulatory mechanisms of virulence.
Therefore, in the present study we used a liquid chromatography-tandem mass spectrometry (LC-MS/MS) approach to analyze the spectrum of RNA modifications in T. gondii strains RH, PRU, and VEG, belonging to different genotypes I, II, and III, respectively. We also examined the correlation between RNA modifications and T. gondii genotype-specific virulence traits. We tested the hypothesis that the patterns and diversity of RNA modifications would differ as a function of T. gondii virulence. Our comparative and systematic investigation provided evidence that different types of RNA modifications exist among T. gondii strains of distinct virulence backgrounds. The results offered baseline data for further investigations of the role of RNA modifications in T. gondii infection biology.
(17 to 50 nt fragments) between T. gondii type I (RH strain), II (PRU strain), and III (VEG strain) using LC-MS/MS-based approach (29). Our analysis identified 22 types of RNA modifications in both tRNA and sncRNA (Table 1). However, due to the low abundance of four RNA modifications (Um, m 5 Um, ac 4 C, and hm 5 C), only 18 types of tRNA and sncRNA modifications were consistently quantified, and thus, shown in the results.
The modified ratio of adenosine (A), uridine (U), cytidine (C), and guanosine (G) was calculated according to the levels of RNA modification. The result showed that tRNA exhibited more modifications than sncRNA, where 1-3% of each nucleotide was modified in tRNA, while , 1% of each nucleotide was modified in sncRNA, except A nucleotide (.2%) in the three T. gondii strains ( Fig. 1A and B). As shown in Fig. 1C, W was the most marked tRNA modification in the three T. gondii strains, followed by I in RH strain and m 1 A in PRU and VEG strains. In contrast, I was the most marked RNA modification in sncRNA in the three T. gondii strains. Notably, W, I, m 1 A, and m 5 C accounted for over 50% of the modified nucleotides in tRNA, whereas I, W, Cm, and Am represented more than half of the quantified modifications in sncRNA (Fig. 1D). Interestingly, the distribution pattern of the quantified RNA modification percentage was different in tRNA but was similar in sncRNA between T. gondii strains ( Fig. 1C and D), indicating that mechanisms for the control of tRNA modification frequency may be strain (genotype) specific.
Further clustering analysis revealed significant differences in the relative expression levels of the quantified tRNA and sncRNA modifications between T. gondii strains ( Fig. 1E and F). The tRNA modification pattern in the virulent RH strain had the most distinct pattern compared with the other two less virulent strains (Fig. 1E), particularly the cluster, including tRNA modification m 2 2 7 G, Am, and Im whose abundance was proportional to the level of T. gondii virulence. The highest expression level of Am, m 2 2 7 G, and Im was detected in RH strain, followed by the less virulent PRU strain and avirulent VEG strain (Fig. 1G). Some sncRNA modifications, such as m 5 C, m 1 A, and m 2 2 G, were also distinctly expressed between the three strains. However, the expression levels were inversely proportional to the parasite virulence, where the highest expression Modifications of T. gondii tRNA and sncRNA Microbiology Spectrum was detected in VEG strain, followed by PRU strain and the lowest expression was detected in RH strain (Fig. 1H). Differentially abundant RNA modifications were derived by comparing two T. gondii strains of different virulence backgrounds. According to Venn diagram analysis, nine tRNA modifications (m 1 A, Am, m 1 I, m 5 U, Cm, m 7 G, m 2 G, m 2 2 G, and m 2 2 7 G) and two sncRNA modifications (m 1 A and m 5 C) were expressed at significantly different amounts between RH, PRU, and VEG strains ( Fig. 1I and J). Two differentially abundant tRNA modifications (W and m 3 U) were uniquely detected between RH and PRU, whereas two (W and m 3 C), one (Gm), and four (Am, Im, m1G, and m 2 2 7 G) differentially abundant sncRNA modifications were detected between RH versus PRU, RH versus VEG, and PRU versus VEG, respectively. PCoA plots showed a strain-specific RNA modification patterns in tRNA and sncRNA, which clearly separated T. gondii strains into three distinct clusters ( Fig. 1K and L).
Linear correlations of different RNA modifications in T. gondii. The correlation matrices and chord diagrams showed that many RNA modifications are simultaneously expressed in tRNA ( Fig. 2A). and sncRNA (Fig. 2B) in the three T. gondii strains. Among all RNA modifications, linear correlations between some modifications were tRNA-specific, such as m 6 A with m 1 A and m 3 U, and m 2 significant positive linear correlation with m 5 C and m 2 2 G, whereas m 6 A had a significant negative linear correlation with Am and Im (Fig. 2D). These results show that correlations between modifications of the same RNA class (tRNA or sncRNA) do occur.
RNA modifications correlate with the expression of T. gondii virulence proteins. As a strict obligate intracellular parasite, T. gondii possesses an impressive spectrum of effector molecules that enables this parasite to interfere with and hijack the host cell machinery to survive and grow inside the host cell. Some of these effectors are produced by specialized organelles, such as micronemes (MICs) (30), rhoptries (ROPs), and dense granules (GRAs) (31,32). MICs play roles in the motility and adhesion of T. gondii, while ROPs and GRAs enable the parasite to evade and modulate host immune defences. RNA modifications have been shown to play a role in mediating T. gondii interaction with the host cell. For example, deletion of the m 6 A components WTAP and METTL3 causes complete inhibition of parasite replication (26). Additionally, most of W in tRNA and mRNA of T. gondii relies on pseudouridine synthase (TgPUS1), which plays a role in the differentiation of T. gondii from acute to chronic, suggesting its important role in mediating the host-pathogen interaction (33).
To investigate the potential role of RNA modifications in T. gondii virulence, we investigated the correlation between the parasite virulence proteins and RNA modifications. We used LC-MS/MS approach and next-generation sequencing to quantify the expression of RNA modifications (Table S1) and genes encoding virulence-related proteins (Table 2; Fig. 3A). Quantitation from the entire RNAseq experiment is provided in Table S2. We detected marked associations between many tRNA modifications and the expression levels of genes encoding virulence proteins in all three T. gondii strains (Fig. 3B). For example, the tRNA modification m 6 A was negatively correlated with MIC1-MIC6, MIC10-MIC11, GRA1-3, GRA5, SAG2, PKG, and CDPK1, but was positively correlated with GRA9. Additionally, various RNA modifications exhibited various correlations with different T. gondii proteins. For example, m 1 A, m 6 A, m 5 U, W, m 3 U, m 5 C, m 7 G, and m 2 G were negatively correlated with some virulence proteins (e.g., MIC1-6); most of those proteins had positive correlation with I, Am, Im, and m 2 2 7 G. We also analyzed the correlation of sncRNA modifications with T. gondii virulence proteins and found moderate correlations between sncRNA modification patterns and T. gondii virulence factors (Fig. 3C). The only significant negative correlations were detected between W, m 3 U, and m 5 C, and ROP1, ROP4, ROP5, and ROP7 (Fig. 3C).

DISCUSSION
RNA modifications are involved in multiple cellular and biological processes in prokaryotes and eukaryotes (10,34). In the present study, LC-MS/MS-based approach was used for comprehensive profiling of RNA modifications in T. gondii strains RH, PRU, and VEG. We analyzed 22 types of RNA modifications, which are commonly identified in tRNA and sncRNA. Of those, Um, m 5 Um, ac 4 C, and hm 5 C were not quantified in tRNAs or sncRNAs of any of the three T. gondii strains, likely attributed to low expression level.
Over 170 types of RNA modifications have been identified across prokaryotes, protista, and fungi (35). In the present study, only 22 types of modified nucleobase standards were used for detection of RNA modifications as previously described (29), limiting the discovery of more RNA modifications in T. gondii strains.
Compared with other RNA classes, tRNA has more modifications, and are involved in regulating tRNA stability, mRNA translation, and protein synthesis (36,37). The modified nucleotides in tRNA are nearly 10 time more (;17% modified residues) than in the other types of RNAs (1% to 2% modified residues) (38). In agreement with previous studies, the present study revealed that tRNA had more modifications than sncRNA, where 1% to 3% of each nucleotide in tRNA was modified compared with ,1% in sncRNA, except A nucleotide (.2%). Interestingly, pseudouridine (W), which exists in all three domains of life (38), was the most abundant tRNA modification in the three T. gondii strains, suggesting its crucial role in T. gondii biology. This result is consistent with a previous study showing that W is one of the most abundant RNA modifications affecting RNA metabolism in T. gondii and that pseudouridine synthase is essential for parasite differentiation (33).
Our analysis revealed a significant diversity of tRNA and sncRNA modification signatures between T. gondii strains with different genotypes. PCA showed that tRNA and sncRNA modification patterns clearly discriminate between the three T. gondii strains. The distribution patterns of the quantified RNA modifications were more diverse in tRNA compared with sncRNA between T. gondii strains, indicating that mechanisms controlling tRNA modification frequency may be dictated by the parasite genotype. Nine modifications were differentially abundant between all three strains in tRNA compared with only two differentially abundant modifications in sncRNA, suggesting that tRNA is more vulnerable for modifications than sncRNA. Interestingly, RH versus PRU had the highest number of tRNA modifications, but the lowest number of sncRNA modifications. Additionally, the expression of tRNA modifications m 2 2 7 G, Am, and Im in RH virulent strain (genotype I) was significantly higher than that of the less virulent strain PRU (genotype II) and avirulent strain VEG (genotype III).
Altered levels of tRNA modifications in yeasts can promote an efficient stress response (39). A previous study showed that modification of 2'-O-methylinosine (Im) in yeast mRNA is involved in the yeast response to environmental stress (40). Therefore, the high level of Im in RH strain compared to the less virulent PRU strain and the avirulent VEG strain suggests its potential involvement in stress response of T. gondii. Alterations of tRNA modifications can affect codon decoding fidelity (41), tRNA structure stability (42), and the generation of various types of tRNA-derived small noncoding RNAs (tsRNAs) (43), which play a role in cell transcriptional and translational control in response to stress (44).
Some sncRNA modifications such as m 5 C, m 1 A, and m 2 2 G were highly expressed in the avirulent VEG strain, compared to the moderately virulent PRU strain and the virulent RH strain. MicroRNA profiling of different T. gondii genotypes revealed both shared and strain-specific microRNAs (45), which may contribute to strain-specific modifications of sncRNA signature. tsRNA fragment, a novel sncRNA type, has been reported in T. gondii and its expression is related to the parasite growth (28). The level of tsRNA is significantly higher in bradyzoite and sporozoite stages of T. gondii compared to the fast-growing tachyzoites. We speculate that differences in the production of tsRNA between the three T. gondii strains may contribute to the differences in the sncRNA modification signature.
Our results revealed linear correlations between various modifications of the same RNA class (tRNA or sncRNA). A previous study in mice detected a correlation between the level of sperm tsRNA m 5 C and the level of sncRNA m 2 G (46). The interaction between tRNA and mRNA modifications has also been reported (47). Whether there are interactions and cross-talks between the different types of RNA modifications identified in the present study remains to be investigated.
T. gondii is defined by three main clonal lineages which differ significantly in virulence capacity mediated by the expression of effector proteins secreted by the specific organelles MICs, ROPs, and GRAs (48)(49)(50). Revealing the differential expression of RNA modifications across different T. gondii genotypes is therefore important for broad understanding of the role of RNA modifications in host-parasite interaction. We examined the correlation between expression of genes encoding virulence-related proteins and RNA modification levels among strains representatives of the three main T. gondii genotypes with possible implications of RNA modifications in the expression of virulence factors. The results revealed correlations between specific RNA modification patterns and virulence gene expression patterns. Although our data did not demonstrate a causal link between RNA modifications and T. gondii virulence, proteins included in the correlation analysis (Table 2) such as GRA12 and ROP18 are established virulence factors of T. gondii (48)(49)(50).
Given that the expression of virulence factor correlates with parasite virulence, and the virulence gene expression patterns correlate with specific RNA modification patterns, we speculate that RNA modifications may contribute to parasite virulence. Although the expression level of T. gondii virulence proteins is tightly regulated, the molecular mechanisms controlling their regulation remain poorly understood. Various RNA elements operate at different levels of gene expression, ranging from the regulation of transcription and translation, dictating RNA conformation, RNA stability, and modulation of protein binding to RNAs/DNAs as well as activity. Further studies are thus needed to investigate whether RNA modifications can influence any of the aforementioned RNA-mediated processes that regulate the transcription or translation of proteins responsible for the expression of virulence-associated traits in T. gondii.
In conclusion, the present study characterized the modified four nucleotides in T. gondii, identified genotype-specific RNA modification patterns, and detected correlations between RNA modifications of the same RNA class, and between RNA modifications and gene expression of virulence proteins. The discovery of RNA modification patterns distinctive T. gondii strains with different virulence properties should facilitate further investigations to obtain more insight into the molecular mechanisms underlying RNA-based control of virulence gene expression in T. gondii. More research is also needed to investigate whether the three archetypical T. gondii haplotypes encode different copy numbers of specific tRNA genes/species.

MATERIALS AND METHODS
Parasite strains. Toxoplasma strains belonging to three T. gondii genotypes, RH strain (type I), PRU strain (type II), and VEG strain (type III), were provided by Dr. Jin-Lei Wang, Lanzhou Veterinary Research Institute. The growth of the tachyzoite stage of these T. gondii strains in monolayers of human foreskin fibroblasts (HFFs) was performed as previously described (48). Heavily infected HFFs were scraped off and passed through 27-gauge needles. The released tachyzoites were purified using a 5 mm Millipore filter and kept in liquid nitrogen until use. The analysis was performed on three biological replicates per T. gondii strain.
Isolation of~80 nt and 17 to 50 nt RNA fragments. The isolation of ;80 nt and 17 to 50 nt RNA fragments was performed according to the manufacturer's instructions. Briefly, 1 mL RNAiso Plus (TaKaRa, Osaka, Japan) was added to a 1.5 mL microcentrifuge tube containing 1 Â 10 6 T. gondii tachyzoites. The mixture was vortexed vigorously and incubated at room temperature for 5 min. Then, 200 mL of chloroform was added to the sample, followed by vortexing and incubation at room temperature for 10 min. Then, the sample was centrifuged at 12,000 Â g for 15 min at 4°C. The supernatant was transferred to a new microcentrifuge tube, mixed with an equal volume of isopropanol, incubated at room temperature for 10 min, centrifuged at 12,000 Â g for 15 min at 4°C, and finally washed with 75% ethanol to dissolve the RNA pellet in RNase-free water. The isolated RNA was electrophoresed in 15% urea-PAGE at constant pressure (200 v) for 1 h in 1Â TBE buffer (Invitrogen, Waltham, MA, USA). The gel was stained with SYBR GOLD (Invitrogen, Waltham, MA, USA), and based on the position of the RNA maker (NEB #N2102, #N0364 Ipswich, MA, USA), RNA fragments of ;80 nt and 17 to 50 nt were extracted from the gel as previously described (46).
Quantitative analysis of RNA modifications by LC-MS/MS. We used LC-MS/MS approach because mass spectrometry methods have been powerful tools for identification and quantification of RNA modifications (51-53). The isolated RNAs were digested into mononucleotides as previously described (46). Then, 100 ng of RNAs were digested in 30 mL reaction system containing 3 mL 10Â RNA hydrolysis buffer (2,500 mM Tirs-HCl, pH 8.0; 50 mM MgCl 2 and 5 mg/mL BSA), 1 IU benzonase (Sigma-Aldrich, St. Louis, MO, USA), 0.2 IU alkaline phosphatase (Sigma-Aldrich, St. Louis, MO, USA), and 0.05 IU phosphodiesterase I (Thermo Fisher Scientific, Grand Island, NY, USA) at 37°C for 3 h. Then, the enzymes in digestion mixture were removed using a Nanosep 3K spin filter (Pall Corporation, Ann Arbor, MI, USA). Before transferring the digested RNA samples, the spin filter was rinsed twice with 400 mL Watsons distilled water at 14,000 Â g for 15 min at room temperature. The processed RNAs were stored at 220°C for LC-MS/MS analysis. The LC-MS/MS-based RNA modification detection was performed using waters ACQUITY-UPLC I-class carrying Xevo-TQ-S mass spectrometry system (46). The expression level of each RNA modification was quantified according to a standard curve and RNA modification relative quantification methods previously described (46).
RNA-sequencing analysis. We used next-generation sequencing to obtain the transcriptomes of genes associated with virulence properties in T. gondii strains RH, PRU, and VEG. The Nanodrop 2000/ 2000c spectrophotometer was used to determine the concentration and quality of RNA extracted from all samples (Thermo Fisher Scientific, USA). The RNA libraries were prepared from 1 mg total RNA of each T. gondii sample using the NEBNext Ultra RNA Library Preparation Kit (Illumina, NEB, USA) in accordance with the manufacturer's recommendations. Each library was sequenced on the Illumina HiSeq 2500 (Illumina, San Diego, CA, USA). Adapter sequences and low-quality reads were removed using Trimomatic (v.0.39) (54). The filtered reads were aligned onto T. gondii reference genome (ftp://ftp.ensemblgenomes .org/pub/release-31/protists//fasta/toxoplasma_gondii/dna/) using ersion 2.0 with default parameters (55). The HTSeq software was used to count the read numbers mapped into each gene and the fragments per kilobase of transcript per million mapped reads (FPKM) was used to measure the expression level of mRNAs in each sample.
Analysis of RNA modifications and correlations. Statistical analyses were performed using GraphPad Prism 8 (GraphPad Software Inc., CA, USA). Differences in the levels of ribonucleosides in different strains were examined by Student's t test and one-way ANOVA with multiple comparisons. Data are presented as the means 6 standard error of means (SEM) of three independent replicates. P values , 0.05 Modifications of T. gondii tRNA and sncRNA Microbiology Spectrum were considered statistically significant. Linear correlation analysis was performed to study the correlation between RNA modifications or between RNA modifications and RNA modification enzymes. The transcriptomic data of selected genes encoding virulence-related proteins were used to explore the correlation between RNA modifications and virulence-related proteins. Hierarchical clustering heatmaps and correlation matrices were obtained using the R package corrplot and correlations were calculated based on the mean values of the measured modifications using the "Spearman" method. Correlation analysis was used to assess the degree of covariance among the various RNA modifications for each parasite strain, with correlation coefficients calculated using GraphPad Prism 8. Principal coordinates analysis (PCoA) of RNA modification levels was performed using the R package scatterplot3d. The PCoA and chord plots were generated using the R packages scatterplot3d and circlize, respectively. Data availability. RNA-seq reads for the T. gondii have been submitted to NCBI under the accession number PRJNA932999.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.1 MB. SUPPLEMENTAL FILE 2, XLSX file, 0.9 MB. We declare no competing interests.