Author response: A single-parasite transcriptional atlas of Toxoplasma Gondii reveals novel control of antigen expression (2020)

Article Figures and data Abstract eLife digest Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Toxoplasma gondii, a protozoan parasite, undergoes a complex and poorly understood developmental process that is critical for establishing a chronic infection in its intermediate hosts. Here, we applied single-cell RNA-sequencing (scRNA-seq) on >5,400 Toxoplasma in both tachyzoite and bradyzoite stages using three widely studied strains to construct a comprehensive atlas of cell-cycle and asexual development, revealing hidden states and transcriptional factors associated with each developmental stage. Analysis of SAG1-related sequence (SRS) antigenic repertoire reveals a highly heterogeneous, sporadic expression pattern unexplained by measurement noise, cell cycle, or asexual development. Furthermore, we identified AP2IX-1 as a transcription factor that controls the switching from the ubiquitous SAG1 to rare surface antigens not previously observed in tachyzoites. In addition, comparative analysis between Toxoplasma and Plasmodium scRNA-seq results reveals concerted expression of gene sets, despite fundamental differences in cell division. Lastly, we built an interactive data-browser for visualization of our atlas resource. eLife digest Toxoplasma gondii is a single-celled parasite that can infect most warm-blooded animals, but only reproduces sexually in domestic and wild cats. Distantly related to the malaria agent, it currently infects over a quarter of the world’s human population. Although it is benign in most cases, the condition can still be dangerous for foetuses and people whose immune system is compromised. In the human body, Toxoplasma cells infiltrate muscle and nerve cells; there it undergoes a complex transformation that helps the parasites to stop dividing quickly and instead hide from the immune system in a dormant state. It is still unclear how this transition unfolds, and in particular which genes are switched on and off at any given time. To understand this transformation, scientists often measure which genes are active across a group of parasites. However, this approach gives only an ‘average’ picture and does not allow each parasite to be profiled, missing out on the diversity that may exist between individuals. One area of particular interest, for example, is a set of genes called SAG1-related sequences. They code for the ‘molecular overcoat’ of the parasite, an array of proteins that sit on the surface of Toxoplasma cells. More than 120 SAG1-related genes exist in the genome of each Toxoplasma parasite, creating a whole wardrobe of proteins that potentially hide the parasites from the immune system. Here, Xue et al. harnessed a technique called single-cell RNA sequencing, which allowed them to screen which genes were active in 5,400 individual Toxoplasma parasites from different strains. The analysis included both the rapidly dividing form of the parasite (present in the initial stage of an infection), and the slowly dividing form found in people who carry Toxoplasma without any symptoms. The resulting ‘atlas’ contains previously hidden information about the genes used at each stage of parasite development: this included unexpected similarities between Toxoplasma and the malaria agent, as well as subtle differences between two of the Toxoplasma strains. The atlas also sheds light on how individual parasites turns on SAG1-related sequences. It reveals a surprising diversity in the composition of the protein coats sported by Toxoplasma cells at the same developmental stage, a strategy that may help to thwart the immune system. One individual parasite in particular had an unusual combination of coat and other proteins found in both the fast and slow-dividing human forms. This parasite had been grown in human cells, yet a closer analysis revealed that it had activated several genes (including ones encoding the protein coat) that are normally only ‘on’ in the parasites going through sexual reproduction in domestic and wild cats. This new data atlas helps to understand how Toxoplasma are transmitted to and grow within humans, which could aid the development of treatments. Ultimately, a better knowledge of these parasites could also bring new information about the agent that causes malaria. Introduction Toxoplasma gondii is an intracellular protozoan parasite that is thought to infect over a quarter of the world’s population (Pappas et al., 2009). Like some of its Apicomplexan cousins, Toxoplasma undergoes a complex developmental transition inside the host. In intermediate hosts, including humans and virtually all other non-feline, warm-blooded animals, Toxoplasma parasites remain haploid and transition from a replicative, virulent tachyzoite to an encysted, quasi-dormant bradyzoite. This asexual developmental transition is tightly coupled to the clinical progression of Toxoplasma infection. Although acute infection with tachyzoites produces few if any symptoms in healthy human children and adults, infected individuals, if left untreated, progress to a chronic stage wherein tachyzoites transition to bradyzoites that can persist for life in neurons and muscle cells. When infected individuals become immunocompromised, such as in chemotherapy, HIV infection, or organ transplantation (Rabaud et al., 1994; Robert-Gangneux et al., 2015), bradyzoites can reactivate to become tachyzoites, causing severe neurological damage and even death. While no causal link has been established, a population-wide study has uncovered significant association of Toxoplasma infection with schizophrenia in chronically infected humans (Sutterland et al., 2015). Chronic infection in mice has been observed to induce behavioral changes such as loss of aversion to cat urine, which is hypothesized to increase the transmission rate of Toxoplasma to its definitive feline host where sexual reproduction occurs (Vyas et al., 2007). As there are no therapeutic interventions to prevent or clear cysts in infected individuals, understanding how Toxoplasma transitions through its life stages remains of critical importance. The development of in vitro methods to induce Toxoplasma differentiation have facilitated investigation of several aspects of chronic infection, including transition of tachyzoites to bradyzoites (Soête et al., 1994; Jeffers et al., 2018). Bulk transcriptomic analyses of Toxoplasma gondii at distinct asexual stages reveal genetic modules that are expressed in each stage (Buchholz et al., 2011; Manger et al., 1998a; Pittman et al., 2014; Yip, 2007; Cleary et al., 2002; Radke et al., 2005; Chen et al., 2018; Fouts and Boothroyd, 2007), including AP2 transcription factors that are thought to play a role in differentiation (Hong et al., 2017; White et al., 2014); however, transitioning parasites convert to the bradyzoite stage asynchronously and display a high degree of heterogeneity along the developmental pathway and in gene expression (Soete et al., 1993; Watts et al., 2015). Furthermore, parasites within the same tissue cysts have been shown to display heterogeneity in the expression of bradyzoite marker proteins (Ferguson et al., 1994). The transition of tachyzoites to the bradyzoite stage results in an overwhelming majority of mature bradyzoites in the G1 phase of the cell cycle that divide slowly, if at all (Radke et al., 2003; Sinai et al., 2016). Furthermore, tachyzoites exhibit slower growth kinetics immediately prior to the bradyzoite transition (Radke et al., 2003; Jerome et al., 1998). This suggests that parasites exit the cell cycle to differentiate into bradyzoites, a pattern consistent with developmental processes in several other eukaryotic organisms (Ali et al., 2011; Kim et al., 2010). Dissecting these cell cycle aspects of stage conversion requires a more detailed analysis than has been possible with bulk measurement of tachyzoite or bradyzoite populations, or with the use of genetically modified parasites coupled with chemical synchronization of cell cycle progression (Radke et al., 2005; Radke and White, 1998; Conde de Felipe et al., 2008; Behnke et al., 2010). This is because the latter approaches require large quantities of synchronized parasites and can potentially introduce artificial perturbations. Furthermore, bulk measurement fails to distinguish parasite-to-parasite variation that is independent of cell cycle or known developmental processes, potentially missing the phenotypic diversity intrinsic to a population of cells. Single-cell RNA sequencing (scRNA-seq) offers a powerful and unbiased approach to reveal the underlying heterogeneity in an asynchronous population of cells. Droplet and FACS-based approaches have already been applied towards multicellular parasites such as Schistosoma to reveal developmental changes within different hosts (Wang et al., 2018). Recently, scRNA-seq has revealed a surprising degree of heterogeneity in another apicomplexan parasite, Plasmodium (Reid et al., 2018; Ngara et al., 2018; Poran et al., 2017; Howick et al., 2019). Analyses derived from these single-parasite measurements uncovered rare and critical transition events in parasite development that were undetectable in bulk measurements. Combined with novel analytical tools and increase in measurement throughput, scRNA-seq can help facilitate the discovery of regulatory factors that mediate these developmental transitions in a system-wide fashion. Here, we performed scRNA-seq to reconstruct transcriptional dynamics of asynchronous Toxoplasma parasites in the course of cell cycle and asexual development in vitro. Our analysis reveals the existence of hidden cell states and rare parasites that show highly unusual patterns of gene expression associated with specific transcription factors. We also discover that individual parasites vary substantially in the expression of surface antigen genes, suggesting the possibility of a novel form of antigenic variation that may play a crucial role in host immunity evasion. Importantly, we identified a single parasite from our scRNA-seq dataset that displayed an unusual expression pattern of surface antigens, leading us to identify and validate the regulatory role of a previously uncharacterized AP2 transcription factor typically associated with parasites in sexual development. Lastly, we show that despite fundamental differences in their modes of cell division, there are conserved transcriptional programs between the asexual life cycles of Toxoplasma gondii and Plasmodium berghei. Our results combined provide the first comprehensive single-cell atlas of Toxoplasma in the course of asexual development and help reveal that the antigenic repertoire of this parasite is much more heterogeneous than previously appreciated. Results Technical validation of single-parasite sorting and sequencing There are more than a dozen approaches available for single-cell isolation and transcriptome amplification. Based on benchmark comparisons, Smart-seq2 generally has higher sensitivity than competing droplet-based approaches (Svensson et al., 2017; Ziegenhain et al., 2017). We reasoned that sensitive measurement is crucial in our study, given that single Toxoplasma gondii parasites are at least 50-fold smaller in volume than a typical mammalian cell, and thus the average parasite gene is likely expressed with much lower copy number per cell than a typical mammalian gene. For our initial studies, we used the common Type I lab strain of Toxoplasma, RH, grown in vitro in human foreskin fibroblasts (HFFs). Following such growth, individual tachyzoites were released by passage through a narrow-gauge needle and then purified by fluorescence activated cell sorting (FACS) into 384-well or 96-well plates. We then synthesized, amplified, and barcoded cDNA using Smart-seq2 and Illumina Nextera protocols. We reduced the reagent cost in 384-well plates effectively by ten-fold compared to the 96-well format. The sequenced reads were bioinformatically deconvolved and grouped into individual parasites for analysis using modified bcl2fastq and custom python scripts (Materials and methods). A schematic to illustrate our experimental workflow is shown in Figure 1. Figure 1 with 2 supplements see all Download asset Open asset Schematic of single-cell RNA-sequencing (scRNA-seq) based on a modified Smart-seq2 protocol for 384-well or 96-well plate. A table of strain types with the number of sequenced cells and cells that passed quality checking (QC) is provided. To ensure that our workflow efficiently captures single Toxoplasma parasites, we mixed equal numbers of two transgenic lines of RH, one expressing GFP and the other expressing mCherry, and sorted individual parasites into a 384 well plate based on the presence of either green or red fluorescence. After Smart-seq2 amplification, we quantified the expression of GFP and mCherry mRNAs using quantitative polymerase chain reaction (qPCR). Across all 301 wells that we measured, we observed the presence of both GFP and mCherry mRNA in only one well, indicating that the rate of doublet events is below 1% (Figure 1—figure supplement 1a). To address the possibility that the reduced reagent volume in the 384-well format could potentially saturate the reaction chemistry and thus limit quantification range, we sorted varying numbers of RH and quantified with qPCR the mRNA of a gene encoding the abundantly expressed surface protein, SRS29B (SAG1) (Figure 1—figure supplement 1b). The detected amount of SAG1 mRNA present in wells containing single, eight, or fifty RH fell into the expected distributions based on the number of parasites sorted, without signs of saturation, indicating that the assay is capable of quantitative measurement at the single Toxoplasma level. We then proceeded to sort parasites into 384-well based on live/dead staining and sequence 729 RH (612 passed quality control) strain single Toxoplasma parasites from asynchronous populations grown under tachyzoite conditions. We also sorted and sequenced 287 RH parasites (237 passed quality control) into 96-well plate for comparison, which we will discuss further in another section. For Pru and ME49 strains, we collected parasites at several time points post alkaline treatment which induces differentiation from tachyzoites to bradyzoites to follow changes in their expression profiles during in vitro development (Materials and methods), yielding 2655 Pru (2198 passed quality control) and 1828 ME49 (1552 passed quality control) single parasites. RH reads were aligned to the GT1 strain genome, which is the most complete reference for Type I parasites, while Pru and ME49 were aligned to the ME49 strain Type II genome reference. Because many genes encoding Toxoplasma secretion factors and surface proteins are evolutionary products of gene duplication events (Reid, 2015), we expected high sequence similarity amongst a substantial portion of the parasite genes. Thus, we modified our gene counting pipeline to account for duplicated genes by distributing reads across all regions with equal alignment score that passed thresholds (Materials and methods). The reason why we adopted a correction scheme for multiply-mapped reads is because the analysis of co-occurrence, or lack thereof, of pathogenic factors (e.g. surface antigens) hinges on sensitive detection of their expression. We thus faced a choice of increasing the false positive in gene alignment by correcting for multiply mapped reads or increasing its false negative by counting only uniquely aligned reads. We chose to account for multiply mapped genes, which would otherwise be obscured, as we hypothesized that surface antigen expression may vary between individual parasites, which we further discuss in another section. A comparison of counting methods does not reveal significant differences in the observed counts (Figure 1—figure supplement 1c). Further analysis reveals that our modified pipeline recovered the detection of more parasite genes than default parameters (Figure 1—figure supplement 1d). To ensure that poorly amplified or sequenced parasites did not confound our downstream analysis, we filtered samples based on several quality metrics including percent reads mapping to ERCC spike-in sequences, number of genes detected, and sequencing depth (Materials and methods; Figure 1—figure supplement 2a). On average, each sequenced parasite contains 30–50% reads that mapped to Toxoplasma genes encoding proteins (top panel in Figure 1—figure supplement 2b). Most of the unmapped reads are from Toxoplasma’s 28 s ribosomal RNA. The relatively high rate of rRNA contamination was also observed in single-parasite RNA sequencing of Plasmodium (Reid et al., 2018). We suspect this occurred due to non-specific priming as protozoan cells have low RNA input. We normalized for sequencing depth across cells by dividing each read count by the read sum of each corresponding cell and then multiplied by the median of read sum within each dataset to yield ‘count per median’ (CPM). After filtering ERCC spike-in and rRNA genes, we detected on average 862, 1290, and 970 genes per parasite with greater than or equal to two read counts (Materials and methods) in the RH, Pru, and ME49 datasets, respectively (bottom panel in Figure 1—figure supplement 2b). Characterization of our measurement sensitivity based on logistic regression modeling of ERCC spike-in standards (Materials and methods) (Lönnberg et al., 2017) reveals a 50% detection rate of 17, 17, and 26 molecules for RH, Pru, and ME49 datasets, respectively (top panels in Figure 1—figure supplement 2c). The sensitivity of our 384-well Smart-seq2 measurement is comparable to the previously reported range for the 96-well format (Ziegenhain et al., 2017). As expected from our qPCR titration experiment, scRNA-seq measurement of gene expression is quantitative at single parasite resolution based on ERCC standards. We determined that the linear dynamic range of our scRNA-seq measurement spans over three orders of magnitude (bottom panels in Figure 1—figure supplement 2c). Taken together, we demonstrate a scalable and cost-effective approach to measure the transcriptomic changes of individual parasites with high sensitivity and accuracy. Cell cycle landscape of asynchronous Toxoplasma Previous work posited a potential link between bradyzoite development and cell cycle, which poses a significant challenge to the bioinformatic analysis of either process (Radke et al., 2003). To characterize cell cycling changes without confounding contributions from developmental processes, we first analyzed an asynchronous population of Type I RH strain parasites grown under tachyzoite conditions; this extensively passaged lab strain is known to have little propensity to switch to bradyzoites under such conditions (Soête et al., 1994) (Materials and methods). After filtering out genes whose expression levels did not vary significantly between individual parasites, we projected the data with principal component analysis (PCA) (Materials and methods). Interestingly, the first two principal components (PCs) reveal a circular trajectory that coincides with relative DNA content, determined using a cell permeable DNA content stain (top panel in Figure 2a). Unsupervised neighborhood clustering identified five distinct clusters of parasites based on their transcriptional profiles (middle panel in Figure 2a) (Materials and methods). We computed RNA velocity to infer transcriptional dynamics (La Manno et al., 2018; Wolf et al., 2018) and the velocity vector field indicates a net ‘counter-clockwise’ flow of transcriptional changes (bottom panel in Figure 2a) (Materials and methods). We assigned cell cycle phase to the clusters based primarily on change in DNA content (Figure 2—figure supplement 1a) but also considering previous bulk transcriptomic characterization (Behnke et al., 2010; Figure 2—figure supplement 1b). Unsupervised clustering identified two distinct clusters in G1 state, which we have designated as G1a and G1b. We found a list of differentially expressed genes between the two G1 clusters. The G1a cluster is highly enriched for the expression of metabolic genes such as phenylalanine hydroxylase (TGGT1_411100) and pyrroline-5-carboxylate reductase (TGGT1_236070), as well as invasion-related secreted factors such as MIC2 (TGGT1_201780), and On the other cluster is enriched for the expression of reductase and and as well as uncharacterized proteins The relative of and states were determined to be and chemical the between the scRNA-seq data of asynchronous parasites and previously bulk transcriptomic measurement suggests in cluster and cell cycle (Figure 2—figure supplement 1b). This a of as it of cell cycle of a parasite without on chemical which may to Figure 2 with 1 supplement see all Download asset Open asset scRNA-seq dynamics of Toxoplasma cell cycle in population. of the first two principal components in RH data RH cells are by fluorescence measurement of a cell permeable DNA content cells are by cluster and by the state. RNA velocity vector field is on of the with in the of net transcriptional of the most gene expression are displayed along the by in along the of the left to the of and the relative fluorescence of DNA content using the same scheme as in To the of gene expression through the cycle, we a of RH using 2 (Materials and methods). The results show a clear expression pattern for the expressed genes along the (Figure 2b). To further characterize cell cycle expression we genes based on and clustering (Materials and methods). of the in tachyzoites are at different in the cell cycle (Behnke et al., 2010). To and this we the expression for each set of genes based on their in et al., 2008; This the with (bottom panel in Figure 2—figure supplement which also with the of DNA and mRNA content (top panels in Figure 2—figure supplement 1c). On the other we also observed where a given expression was to the of its set in Figure 2—figure supplement 1d). For example, of genes as or are assigned cluster while the genes are assigned clusters 1 or 2 (Figure 2—figure supplement genes as and based on their to genes encoding known are assigned to cluster 2 instead of cluster panel in Figure 2—figure supplement have determined these two to be This is consistent with our of between their and known expression profiles along the et al., 2019). analysis of we also identified genes not as within the cluster such as and that have been shown to and designated as and respectively et al., 2014; et al., panel in Figure 2—figure supplement As another example, in expression in while the majority of complex genes are expressed towards the phase of the cell cycle panel in Figure 2—figure supplement A has of as a protein based on et al., which is consistent with our group of as in cluster one where genes A list of RH genes with their corresponding clustering is We high of expression for several genes in each sets, suggesting that the Toxoplasma may significant Our scRNA-seq data provide an to help identify genes and infer of uncharacterized heterogeneity in Toxoplasma Toxoplasma has one of the most developmental programs of any single-celled however, it is how synchronized the transition is between developmental To address we the heterogeneity within Pru, a II strain that is capable of tissue cysts with that from in growth in in vitro alkaline conditions (Soete et al., 1993; et al., 2017). We applied scRNA-seq to measure and Pru parasites grown in as tachyzoites and the switch to bradyzoites by growth in alkaline for and of the first two of Pru tachyzoites reveals the expected circular (Figure supplement cell cycle progression as for the RH tachyzoites, To validate we a based on our cell cycle in RH (Materials and methods). to we observed in RH, cell cycle reveals that the population of Pru is of and parasites in and with previous et al., our data show most Pru parasites are in the G1 with a of (Figure supplement 1b). To identify transcriptomic changes associated with the we projected data from both and Pru parasites two using a (Materials and methods) et al., 2018). Unsupervised clustering revealed distinct clusters of parasites, which we (Figure with treatment time points and cell cycle states (Figure Figure supplement suggesting that the asexual differentiation with cell cycle in Toxoplasma, as previously (Radke et al., 2003). We the by post alkaline and observed expression of previously bradyzoite marker genes in parasites, including and with a in expression of SRS29B a surface marker gene (Figure supplement The of parasites in the population suggests two possible the of SAG1 mRNA is relatively and we are SAG1 the parasites were still tachyzoites, or the asexual transition by alkaline treatment is highly Interestingly, RNA velocity analysis suggests that may be a as the trajectory into either or as by the net transcriptional flow Figure to panel in Figure Figure with supplements see all Download asset Open asset scRNA-seq dynamics of asexual development in vitro. of and Pru parasites with indicating cluster panel the number of parasites in each of Pru parasites or by post cell cycle and RNA velocity from left to of differentially expressed genes across clusters of cells by clustering The five most enriched genes from each cluster are of Pru by the expression of bradyzoite (top and tachyzoite (bottom red marker genes. of average expression of AP2 transcription factors for each normalized by the cluster expression within each and green enriched in clusters and To the gene modules specific to a given we gene expression for each cluster (Figure and