Differential Impacts on Host Transcription by ROP and GRA Effectors from the Intracellular Parasite Toxoplasma gondii

This work performs transcriptomic analysis of U-I cells, captures the earliest stage of a host cell’s interaction with Toxoplasma gondii, and dissects the effects of individual classes of parasite effectors on host cell biology.


(i) FACS-based isolation of host cells for single-cell RNA sequencing purifies uninfected-injected cells.
To generate host cells impacted by various combinations of parasite-dependent stimuli for downstream RNA-seq, we subjected the host cells to infection with either wild-type RH (type I) strain parasites, mutant Δmyr1 parasites lacking the MYR1 protein (constructed from an RH Δmyr1 mCherry parental strain [25]) (see Fig. S1 in the supplemental material), or a parasite-free cell lysate (i.e., mock infection) from feeder human foreskin fibroblast (HFF) cells used to maintain both parasite lines. From here on, we designate host cells that arose from the wild-type, Δmyr1, and mock infections Wt, Dmyr1, and mock cells, respectively. All parasites were engineered to express the ROP fusion protein toxofilin-beta-lactamase (Tfn-BLA) as well as a constitutive red fluorescent protein. In addition, we chose 10 T1/2 mouse fibroblasts as the preferred host cell type over conventionally used HFFs, which we found unsuitable for single-cell sorts due to their propensity to remain clumped in FACS buffer.
As each true infection results in a heterogeneous monolayer containing U-I cells, infected cells (designated I-I for being infected and injected), and uninfected cells (denoted U-U for being uninfected and uninjected) (Fig. 1B), the infection conditions yielded seven key species of host cell: Wt U-I, Wt I-I, Wt U-U, Dmyr1 U-I, Dmyr1 I-I, Dmyr1 U-U, and mock. Each of these species was expected to be influenced by a specific combination of parasite-dependent stimuli (Fig. 1C). By comparing the transcriptomes of key pairs of these seven cell species (see "Data availability" in Materials and Methods, below), we expected to uncover the impact of previously unexamined parasite-dependent stimuli, e.g., ROPs (summarized in Fig. 1D). Of note, these experimental conditions enabled transcriptomic assessment of paracrine effects during Toxoplasma infection (by comparison of U-U versus mock-infected cells) and generated pure I-I and U-I cell populations at extremely short infection durations, an improvement over traditional methods that rely on high multiplicities of infection (MOIs) and long infection times to distinguish infected versus uninfected transcriptomic signatures.
To capture individual cells from each cell species, we devised a FACS-based protocol that purified each species from the host cell monolayers ( Fig. 2A). The pipeline employed the Tfn-BLA assay (8), in which the (mock-)infected host cell monolayers were stained for 30 min with CCF2-AM, a reporter dye that is taken up by live host cells and shifts from green to blue fluorescence when cleaved by beta-lactamase. CCF2-AM staining occurred at 30, 90, or 150 min postinfection, to yield cells at 1, 2, or 3 h postinfection (hpi). In addition, CCF2-AM staining was performed at room temperature (23°C) to preserve the integrity of the CCF2-AM dye, which decomposes at 37°C. U-I cells were then sorted by FACS based on their blue (cleaved) CCF2-AM fluorescence and their absent red fluorescence (since they lacked internalized parasites). In addition, I-I cells were sorted for their double-positive blue (from cleaved CCF2-AM) and red (from internalized parasites) fluorescence, and U-U and mock cells were sorted for their green (cleaved) CCF2-AM fluorescence. Single cells and bulk populations of 50 to 100 cells of each cell species were collected.
To ensure confidence in the identity of the sorted cells and to limit the number of parasites per I-I cell to one, we gated stringently during FACS, particularly for the "U-I cell" gate and the "parasite-associated" gate (which included I-I cells) (Fig. 2B). We also subjected all cells to extensive washes before FACS, as debris from the human feeder cells used to culture the parasites contaminates the U-I gate (Fig. S2), likely due to the retention of both CCF2-AM dye and the parasite fusion protein Tfn-BLA. Of note, the debris posed significant challenges to downstream bioinformatic analyses of bulk U-I cells, but they registered as single cells in the scRNA-seq pipeline and were automatically discarded during quality control due to their high percentage of human reads. At a multiplicity of infection of ϳ6, U-I cells constituted ϳ8.8% of the injected (i.e., U-I plus I-I) cell population at 1 hpi, ϳ2.0% at 2 hpi, and ϳ3.4% at 3 hpi (Fig. 2C). Abundances of sorted U-I, I-I, and U-U cells relative to all sorted cells are depicted in Fig. 2D.
To validate U-I cells as appropriate models for host responses to injected parasite effectors, we also limited the possibility of U-I cells arising by mechanisms other than aborted invasion events. One such mechanism is host immune clearance of internalized parasites. However, this route is an unlikely source of U-I cells in our pipeline for two reasons: (i) immune clearance of intracellular pathogens in mouse cells is activated by exogenous IFN-gamma, which we did not add to the host cells, and (ii) all known strains of Toxoplasma escape IFN-gamma-mediated pathogen clearance by suppressing the IFN-gamma signaling pathway if they infect the host cell before it is exposed to extracellular IFN-gamma (33,34). In another potential mechanism, an infected host cell may divide and produce an uninfected daughter (a U-Id cell, for U-I by division). To limit the occurrence of U-Id cells, we serum starved all host cells for the 24 h preceding infection to induce cell cycle arrest (Fig. S3A) and limited infection durations to 3 hpi or less, as live-video microscopy of 200 parasite-infected, serum-starved 10 T1/2 host cells for 16 h revealed that no infected cell divided before 3.67 hpi (Fig. S3B). Finally, a newly internalized parasite might spontaneously exit early from the host cell, a rare and poorly characterized process that is therefore difficult to control. Although we cannot exclude the possibility that a very small number of our U-I cells arose from this mechanism of premature parasite egress, and although U-Id cell production cannot be completely ruled out, an advantage of scRNA-seq is its ability to distinguish differences between cells arising via different mechanisms (assuming that they differ transcriptomically).
(ii) Technical validation of single-cell sequencing for cells exposed to or parasitized by Toxoplasma gondii. Given the limited duration of infection (Յ3 h), we expected relatively subtle transcriptional changes in U-I and I-I cells. To maximize the sensitivity of scRNA-seq analyses to detect such signatures, we employed Smart-seq2 library preparation and sequenced to a depth of 1 million reads per cell. All reads were aligned to a concatenated mouse-Toxoplasma genome using the genomic sequence for the GT1 parasite strain, a clonal relative of RH (35,36). In addition, we identified and discarded 43 mouse genes to which reads from extracellular RH parasites aligned (mostly representing evolutionarily conserved genes) (Table S1).
To ensure that poorly amplified or poorly sequenced host cells did not confound downstream analysis, we filtered samples based on several quality metrics (see Materials and Methods), yielding 453, 2,026, and 2,875 cells at 1, 2, and 3 hpi, respectively, for downstream analysis (Fig. S4A). The mapping efficiency for the analyzed cells was Ͼ85% in all experimental trials (Fig. S4B). Characterization of measurement sensitivity based on logistic regression modeling of ERCC spike-in standards revealed a 50% rate of detection of 31, 11, and 34 RNA molecules per cell at 1 hpi, 2 hpi, and 3 hpi, respectively (Fig. S4C), amounting to a level of sensitivity comparable to that previously reported (37). After filtering out genes with expression levels above the detection limit in Ͻ6 cells, we normalized for sequencing depth across cells by dividing each read count by the median read sum to yield counts in units of counts per median (cpm) (see Materials and Methods). Gene expression in the scRNA-seq data set exhibited a strong positive correlation to expression in bulk RNA-seq samples processed in the same way, specifically for differentially expressed genes (DEGs) between experimental conditions (i.e., U-U, I-I, and U-I cells) (Fig. S4D). Furthermore, each bulk experiment's expression data exhibited the best correlation with their cognate single-cell data (Fig. S4E). Taken together, these data demonstrate that scRNA-seq captures transcriptomic signatures in host cells similar to those detected in bulk RNA-seq experiments, validating the scRNA-seq platform as an approach to further characterize host cell transcriptomic responses during Toxoplasma infection.
(iii) Single-cell resolution reveals cell-to-cell heterogeneity inaccessible to the bulk RNA sequencing platform. A key advantage of single-cell resolution is that it enables the identification and separation of parasite-independent host cell heterogeneity. Accordingly, single-cell resolution facilitates an extra checkpoint to validate the identities of U-I, I-I, and U-U cells by quantifying the percentage of Toxoplasma-derived reads in each cell. As expected, the Toxoplasma read content across all cells exhibited a bimodal distribution, where most uninfected cells contained Ͻ0.01% Toxoplasma reads, while most infected cells possessed higher percentages, i.e., 0.5 to 4%. However, small proportions of U-I cells and I-I cells exhibited unexpectedly high (for U-I) or low (for I-I) Toxoplasma read counts and were considered to be misclassified (Fig. S5A). This may have resulted from a low rate of TdTomato or mCherry loss in some of the parasites (for cells misclassified as U-I) or from attached but not fully invaded parasites being dislodged from the host cell at some point between fluorescence detection and the deposition of the cell into lysis buffer (for cells misclassified as I-I). Such misclassified cells may have contributed significant and potentially misleading signatures to their respective samples. To preclude this possibility in our single-cell data set, we excluded all U-I cells with Ͼ0.04% Toxoplasma reads and all I-I cells with Ͻ0.32% Toxoplasma reads from downstream analysis.
Next, because the experimental pipeline examined cells during their earliest interaction with Toxoplasma, we expected the subtle, parasite-dependent transcriptional signatures to be potentially eclipsed by intrinsic host processes such as the cell cycle, which still progressed in some cells even with serum starvation (Fig. S3A). To separate the host cell cycle from other biological processes that the parasites could potentially modulate, we determined the phases of all single cells (i.e., G 1 , S, or G 2 /M) based on the expression of 175 curated cell cycle marker genes ( Fig. S5B and Table S2) (see Materials and Methods). A breakdown of the cell cycle phase composition under each experimental condition across all time points revealed a consistent pattern at 2 hpi in which the proportion of cells in G 2 or M phase increased in a manner that appeared dependent on injected effectors, i.e., from 13.3% in mock cells to 33.9% in WT U-I, 41.4% in Dmyr1 U-I, and 22.5% in cytochalasin D-treated U-I (CytD U-I) cells (Fig. S5C). These findings are consistent with the potential induction of cell cycle arrest in U-I cells by injected effectors, as previously reported (38)(39)(40)(41).
To identify other sources of heterogeneity in the data set, we used the uniform manifold approximation and projection (UMAP) algorithm for dimensionality reduction and for visualization of relationships between all cells based on the most dispersed (i.e., variable) genes. Leiden clustering revealed three distinct clusters, here designated populations 1, 2, and 3 (Fig. S5D). Curiously, cells from different time points exhibited differences in the proportion of cells belonging to each population, where cells at 2 hpi fell entirely in population 1 and cells at 3 hpi fell in all three populations. In conventional single-cell analysis, dimensionality reduction and cell clustering are used to identify novel cell types. Although the 10 T1/2 host cell line used in our experiments is clonal, it is derived from a pluripotent stem cell population that has a propensity to differentiate such that the three identified populations could conceivably represent distinct differentiation states (42).
To limit cell-to-cell heterogeneity from overshadowing potentially subtle transcriptomic signatures induced by parasite effector secretion and parasite invasion, we limited our remaining analyses to host cells in G 1 phase (to limit cell cycle signatures) and from UMAP population 1 (to factor out potential cell type signatures).
Injected parasite effectors drive inflammatory transcriptional signatures associated with parasite infection. Having established a robust pipeline to isolate, RNA sequence, and bioinformatically analyze the transcriptomes of individual cells from each of seven relevant cell species (i.e., Wt U-I, Wt I-I, Wt U-U, Dmyr1 U-I, Dmyr1 I-I, Dmyr1 U-U, and mock cells) (Fig. 1C), we next sought to determine host responses to each of five individual classes of parasite-dependent stimuli, namely, rhoptry proteins (ROPs), MYR-dependent dense granule proteins (MDGs), MYR-independent dense granule proteins (MIGs), parasite invasion, and paracrine effects. To isolate each parasitedependent stimulus, we used the model-based analysis of single-cell transcriptomics (MAST) algorithm to perform differential gene expression analysis on key pairs of the seven relevant cell species from the pool of G 1 -phase, UMAP population 1 cells. Across all pairwise comparisons between conditions within each time point, 39, 2,252, and 10,995 DEGs were detected at 1 hpi, 2 hpi, and 3 hpi, respectively (Table S3A). These results show that parasite effectors and invasion result in almost no detectable transcriptional response at 1 hpi, a modest response involving the regulation of a core group of genes at 2 hpi, and ramping up of this response at 3 hpi, in which transcription in a much larger set of host cells is modulated. Because of the negligible response at 1 hpi, we excluded data from this time point from the remaining analyses.
For additional validation, we assessed the 2-hpi and 3-hpi data sets for their agreement with data from two previous bulk RNA-seq studies that captured the host response to infection at 6 hpi in nominally the same parasite strains and under the same culture conditions albeit in HFFs instead of mouse 10 T1/2 fibroblasts (10,31). Because these studies measured the infection response by comparing infected versus mock-infected cells, we examined DEGs between Wt I-I and mock-infected (i.e., mock) cells in our own data sets. Gene set enrichment analysis (GSEA) of the resulting DEGs at 2 hpi and 3 hpi (Table S3B) using the Molecular Signatures Database's Hallmark gene sets (43) revealed that nearly all significantly enriched gene sets at 2 hpi and 3 hpi were previously identified in the reference data set of Naor et al., and many pertained to immune processes (Fig. 3) (31). Of note, gene sets were considered to be significantly enriched if their false discovery rates (FDRs) were Ͻ0.25, a standard cutoff for GSEA given the lack of coherence in most transcriptional data sets and the relatively low number of gene sets being analyzed. Nearly all gene sets not preserved in the data set of Naor et al. were enriched from genes downregulated upon infection. The lack of agreement between our downregulated gene sets and those obtained previously likely reflects the infection response's tendency toward gene upregulation and stochasticity in expression levels for the substantially fewer downregulated genes detected, an interpretation corroborated by the downregulated genes' higher P values and lower fold changes. Overall, these results indicate that the 2-hpi and 3-hpi data sets capture well-known host responses to infection with Toxoplasma gondii.
Next, to determine the host response to injection of parasite effectors, we compared U-I cells from wild-type parasite infection (i.e., Wt U-I cells) to uninfected cells from the  same monolayer (i.e., Wt U-U cells). At 2 hpi, only 10 DEGs were identified (Table S3C), and GSEA of these DEGs revealed no significant enrichment of the Hallmark gene sets. Therefore, injection of parasite effectors without invasion appears to elicit only a trace response at 2 hpi. At 3 hpi, 156 DEGs were detected between Wt U-I and Wt U-U cells (Table S3C), for which GSEA revealed a total of 17 gene sets (Fig. 4A); 14 corresponded to genes expressed at higher levels in Wt U-I cells than in Wt U-U cells, and several were associated with inflammation. In addition, 10 gene sets were common to the infection response, i.e., also enriched in Wt I-I versus mock DEGs (Fig. 3), which suggests that much of the early response to parasite infection is driven by the injection of parasite effectors (likely ROPs) into host cells prior to invasion.
To confirm that the transcriptional signatures in Wt U-I cells originated from ROP injection, we made two additional comparisons. For the first, we compared U-I cells from Δmyr1 parasite infection (i.e., Dmyr1 U-I cells) to uninfected cells from the same monolayer (i.e., Dmyr1 U-U cells). As MDGs fail to traverse the PVM during Δmyr1 parasite infections, the absence of MYR1 should effectively limit the parasites' effector repertoire to ROPs and MIGs; consequently, transcriptional signatures from Dmyr1 U-I cells and Wt U-I cells should largely resemble one another. In the second comparison, to account for the possibility that Wt U-I signatures originated from a preexisting difference in those host cells from their neighbors, rather than from injected effectors, we collected U-U, U-I, and I-I cells from host monolayers exposed to wild-type parasites pretreated with the invasion inhibitor cytochalasin D, which allows parasite attachment but blocks subsequent invasion. As expected, cytochalasin D treatment increased the proportion of U-I cells in infected monolayers by ϳ6-fold at 2 hpi and by ϳ9.3-fold at 3 hpi (Fig. S6), so at least 85 to 90% of the artificially induced U-I cells (i.e., CytD U-I cells) presumably arose from drug-induced abortion of parasite invasion events rather than from parasite-independent host cell differences. We predicted that transcriptional Impact of Toxoplasma ROPs/GRAs on Host Transcription ® signatures detected in CytD U-I cells would be driven almost entirely by injected parasite effectors and would mirror those identified in both Wt U-I and Dmyr1 U-I cells.
As expected, at 2 hpi, very few DEGs (29 and 51, respectively) and no enriched gene sets were identified for Dmyr1 and CytD U-I cells versus their corresponding U-U cells, while at 3 hpi, substantially more DEGs (103 and 174 for Dmyr1 and CytD U-I versus U-U cells, respectively) and gene sets were identified (Tables S3D and E). The 3-hpi gene sets exhibited a strong degree of overlap with the original Wt U-I-versus-Wt U-U comparison, such that 14 of the 17 Wt U-I-versus-Wt U-U gene sets were also identified in one or both of the corresponding Dmyr1 and CytD comparisons (Fig. 4B). These data are consistent with an inflammatory response to ROP injection that drives much of the host cell's total response to parasite infection.
In light of the many genes and gene sets shared between the infection and injection responses, we next sought to define the distinction between these responses by comparing Wt U-I (injected) to Wt I-I (infected) cells using two complementary approaches. In the first approach, we used MAST to identify DEGs between Wt U-I and Wt I-I cells and subjected the DEGs to GSEA. To better compare the trajectories of Wt I-I and Wt U-I gene expression, we also computed fold changes in DEG expression between each of these conditions and Wt U-U cells (Table S3F). In the second approach, we identified genes with significant differences between Wt U-I and Wt I-I cells in their correlation to a quantity called the CCF2-AM ratio, i.e., the ratio of blue (cleaved) to green (uncleaved) CCF2-AM indicator dye detected during FACS.
In this second approach, the CCF2-AM ratio was used as a quantitative proxy for the influence of parasite-dependent effectors on individual host cells. Briefly, in the strictest sense, the CCF2-AM ratio is a biological readout for the penetration of a given host cell by the injected ROP fusion protein Tfn-BLA, as the intracellular CCF2-AM dye in our pipeline is cleaved by the beta-lactamase in Tfn-BLA. Since the host cells should exhibit more or less equal loading of the CCF2-AM substrate, we presumed that the extent of conversion of the CCF2-AM signal from green (uncleaved) to blue (cleaved) reflected both the concentration of Tfn-BLA protein introduced into the cell and the amount of time that it had spent within the cell. In Wt U-I cells, Tfn-BLA penetration occurs concomitantly with the injection of the other ROP effectors; therefore, the CCF2-AM ratio can be interpreted as a quantitative measure for ROP penetration in Wt U-I cells.
Wt I-I cells, however, are presumably penetrated not only by ROPs but also by MIGs, MDGs, and the parasites themselves. Accordingly, since the amount of ROPs injected into cells and the time that ROPs spend inside cells likely track with the same quantities for the remaining parasite-dependent stimuli, we used the CCF2-AM ratio in Wt I-I cells as a proxy for the presence of all four parasite stimuli inside each cell. Next, we computed the Spearman correlation between each gene's expression and the CCF2-AM ratio separately in both Wt U-I cells and Wt I-I cells. Of note, because we calculated the correlation by incorporating the CCF2-AM ratios from Wt U-U cells (which should be devoid of ROPs, MIGs, MDGs, and parasites) as negative controls for both the Wt U-I and Wt I-I analyses, the correlation data do not encapsulate the influence of parasitedependent stimuli that Wt U-U cells have in common with Wt U-I and Wt I-I cells, i.e., paracrine factors. Finally, we identified genes that exhibited significant shifts between Wt U-I and Wt I-I cells in their relationship to the CCF2-AM ratio (Table S4A) by first modeling a Gaussian distribution of their Spearman correlations to the CCF2-AM ratio (Fig. 5A). Genes that significantly deviated from the Gaussian distribution and that also exhibited a sufficient difference in their CCF2-AM correlation scores in Wt U-I versus Wt I-I cells were interpreted to be associated with either ROP injection alone (in Wt U-I cells) or with all secreted effectors plus parasites (in Wt I-I cells, in which ROP penetration should be accompanied by all other parasite-induced insults).
Analysis of Wt I-I versus Wt U-I cells at 2 hpi revealed no significant DEGs, pointing to a profound similarity in host transcription in Wt I-I and Wt U-I cells at 2 hpi. However, the corresponding CCF2-AM correlation data exhibited higher sensitivity for detecting differences between Wt I-I and Wt U-I cells and yielded 81 significant deviants from the Gaussian distribution (Fig. 5A, left). Of the deviants, 73 (ϳ90.1%) fell in regions of the scatterplot where genes either (i) were negatively correlated with the CCF2-AM ratio under one condition and positively correlated under the other or (ii) had stronger positive or negative correlations in Wt U-I cells than in Wt I-I cells as determined by a Gaussian model (Fig. 5A, left, blue regions; Table S4A). This provides evidence for downstream events in infection (i.e., the release of dense granule proteins or parasite invasion) suppressing the effects of genes induced by ROP injection by as early as 2 hpi.

Impact of ROPs (in Wt U-I cells) enhanced by MIGs + MDGs + Parasites (in Wt I-I cells)
Wt U-I Correlation  (Table S3F); i.e., they were either (i) upregulated in one cell type and downregulated in the other compared to a common Wt U-U standard or (ii) more upregulated or downregulated in Wt U-I cells compared to the standard than in Wt I-I cells. In addition, of the CCF2-AM correlation deviants, 333 (ϳ79%) exhibited at least a moderate negative correlation with the CCF2-AM ratio in Wt I-I cells, an observation consistent with these genes being downregulated in response to the combination of all parasite-derived insults and effectors during infection (Fig. 5A, right; Table S4A). Furthermore, 282 (ϳ67%) of the deviants fell in regions of the scatterplot in which the impact of ROPs is dampened by dense granule proteins (GRAs) and parasite invasion (Fig. 5A, right, blue regions; Table S4A), which is consistent with host responses in Wt U-I cells being counteracted by parasite effectors that are operating in Wt I-I cells only. This lends further support to the notion that effectors injected circa invasion (i.e., ROPs) are counterbalanced by subsequently introduced effectors (i.e., GRAs). Of note, the remaining 141 CCF2-AM deviant genes at 3 hpi fell within the sections of the scatterplot in which the effects of ROPs are enhanced by those of the GRAs and parasite invasion (Fig. 5A, right, yellow regions; Table S4A). This suggests that while GRA release and parasite penetration may counterbalance some ROP-induced host responses, these events may also enhance other genes induced by ROP injection.
To discern the biological significance of the differences between the infection and injection responses, we performed GSEA on Wt U-I versus Wt I-I DEGs. The resulting enriched gene sets (Fig. 5B) included many that were previously identified in the host infection response, with gene sets enriched from genes expressed at higher levels in Wt I-I cells corresponding to more inflammatory processes. Because the vast majority of genes with higher expression levels in Wt U-I cells than in Wt I-I cells exhibited evidence of ROP effectors being counteracted by subsequently secreted effectors such as MDGs and MIGs (84 out of 93 genes) (Table S3F), these findings are consistent with injectionassociated inflammatory host processes ramping up upon parasite penetration (and GRA release), while other injection-associated processes are dampened by these later events.
Parasite effectors counteract one another to yield a modest host response to early parasite infection. Next, we sought to determine the contribution of individual parasite effector compartments to the difference between the host infection (in Wt I-I cells) and ROP injection (in Wt U-I cells) responses. Because Wt U-I and Wt I-I cells are both injected with ROPs and originate from the same monolayer, ROPs and paracrine effectors likely do not explain the differences between these two cell types. Instead, MDGs, MIGs, and the act of parasite penetration itself are, a priori, most likely to explain these differences. Our data set, which includes infected (I-I), bystander uninfected (U-U), and U-I cells originating from monolayers infected with wild-type and Δmyr1 parasites, presents a unique opportunity to examine the impacts of these compartments in isolation.
To determine how individual effector compartments contribute to the transcriptomic differences between Wt I-I and Wt U-I cells, we made two key comparisons. In the first comparison, we examined Dmyr1 U-I cells, which presumably respond to ROP injection, and Dmyr1 I-I cells, which are thought to respond to ROP injection plus MIGs and parasite penetration ( Fig. 1C and D). We predicted that MIGs and parasite penetration would account for the differences between Dmyr1 I-I and Dmyr1 U-I cells and for a subset of the differences between Wt U-I and Wt I-I cells. In the second comparison, we examined Wt I-I versus Dmyr1 I-I cells, which originate from monolayers infected with wild-type and Δmyr1 parasites, respectively. Signatures detected in Dmyr1 I-I cells likely reflect the host response to the combination of ROP injection, MIG activity, parasite penetration of host cells, and paracrine effectors secreted into the extracellular milieu, while those detected in Wt I-I cells likely reflect the response to these elements plus MDG secretion (Fig. 1C). Therefore, we predicted that comparing Dmyr1 I-I versus Wt I-I cells would illustrate the impact of MDGs as well as any paracrine effects dependent on the presence of MYR1 (Fig. 1D) and that this would account for a second subset of the differences between the infection and injection responses showcased in the Wt I-I-versus-Wt U-I comparison.
(i) MYR-independent GRAs and parasite penetration collectively enhance ROPinduced inflammatory responses. To examine the host response to the combination of MIGs and parasite penetration, we compared Dmyr1 U-I cells, which should be penetrated by ROPs, and Dmyr1 I-I cells, into which parasites should secrete ROPs and MIGs but not MDGs ( Fig. 1C; Table S3G). At 2 hpi, only 5 DEGs and no GSEA-enriched gene sets were identified, implying a profound similarity between Dmyr1 U-I and Dmyr1 I-I cells at this time point (as was also the case for Wt U-I versus Wt I-I cells). At 3 hpi, 80 DEGs were identified, pointing to a slight divergence between these two cell types at this time point. Forty-eight (60%) of the DEGs exhibited evidence of being influenced by effectors that counteract one another's effects (i.e., their expression exhibited either [i] opposing trends in Dmyr1 U-I versus Dmyr1 I-I cells compared to an uninfected Dmyr1 U-U standard or [ii] stronger induction or suppression in Dmyr1 U-I cells than in Dmyr1 I-I cells), which suggests that MIGs may play a role in neutralizing the effects of the ROPs preceding them. Note, however, that we cannot exclude the possibility that this effect is attributable to stimuli related to physical penetration by the parasites.
GSEA of the 3-hpi DEGs identified 15 enriched gene sets (Fig. 6), 9 of which were also found to be enriched in the Wt U-I-versus-Wt I-I comparison (Fig. 5B). For gene sets enriched from genes expressed at higher levels in Dmyr1 I-I cells than in Dmyr1 U-I cells, nearly all corresponded to gene sets already identified as part of the host response to Impact of Toxoplasma ROPs/GRAs on Host Transcription ® ROP injection and/or were associated with inflammatory processes. In contrast, none of the gene sets enriched from genes expressed at higher levels in Dmyr1 U-I cells were identified as part of the injection response and instead corresponded to other processes, i.e., the complement cascade, coagulation, MTORC1 signaling, and myogenesis. Taken together, these results suggest that MIGs and parasite penetration indeed account for some of the difference between the infection response (in Wt I-I cells) and the injection response (in Wt U-I cells). More specifically, genes corresponding to inflammatory processes that are already induced upon ROP injection (in Wt U-I and Dmyr1 U-I cells) appear to be further induced, i.e., enhanced, by the combination of MIGs and parasite penetration, while genes for which the effect of MIGs plus parasites dampens the influence of the ROPs appear to correspond to a different set of cellular processes. Although the latter set of genes accounts for the majority (60%) of Dmyr1 U-I versus Dmyr1 I-I DEGs, the majority of the total enriched gene sets correspond to cellular processes enhanced rather than dampened by MIGs and parasites, and these enhanced gene sets collectively exhibit much higher statistical significance. This raises the possibility that at least some processes dampened by MIGs and parasites may not be adequately captured by the Hallmark gene sets.
(ii) MYR-dependent GRAs counterbalance parasite effectors released earlier in the lytic cycle. A previously reported bulk RNA-seq experiment comparing host transcription in monolayers infected with wild-type versus Δmyr1 parasites at 6 hpi (31) revealed a set of genes such that in the absence of MYR1, expression changes were unmasked, while in the presence of MYR1 (i.e., during wild-type infection), there was no net change. This work implies that, collectively, MDGs and associated paracrine effectors secreted during infection with wild-type but not Δmyr1 parasites participate in counterbalancing transcriptional signatures induced by prior parasite-dependent stimuli.
To determine whether MDGs and associated paracrine effects play a similar role in the present single-cell data sets and to ascertain their contribution to the difference between the infection and injection responses, we compared Dmyr1 I-I and Wt I-I cells.
Analysis of this comparison using the MAST algorithm identified 46 and 367 DEGs at 2 hpi and 3 hpi, respectively (Table S3H). While GSEA of the 2-hpi DEGs returned no significantly enriched gene sets, GSEA at 3 hpi revealed 16 significantly enriched gene sets, all but one of which were enriched from genes expressed at higher levels in Dmyr1 I-I cells than in Wt I-I cells (Fig. 7A). Furthermore, 8 gene sets were identified in the previously reported bulk RNA-seq experiment comparing infections with wild-type versus Δmyr1 parasites (31), while 6 corresponded to infection-associated gene sets (i.e., they were enriched from Wt I-I versus mock DEGs), and 8 were identified as part of the ROP injection response (i.e., they were enriched from Wt U-I versus Wt U-U DEGs). Taken together, these data are consistent with MDGs and/or their associated paracrine effects selectively impinging on Wt I-I cells and dampening the effects of parasite effectors introduced into host cells earlier in the lytic cycle. The overlap between these gene sets and those of the injection response shows that some of the effectors whose responses were dampened by MDGs and associated paracrine effects include injected effectors (i.e., ROPs), which was previously suspected (31) but never before explicitly demonstrated.
Although comparing Dmyr1 I-I versus Wt I-I cells illustrates the collective impact of MDGs and MYR-dependent paracrine effects on host transcription (Fig. 1D), this comparison cannot distinguish the individual impacts of these two stimuli. To determine the effect of MDGs alone, we performed CCF2-AM ratio correlation analysis on Wt I-I versus Dmyr1 I-I cells (Table S4B), as was performed for Wt I-I versus Wt U-I cells. As was the case for the computation of both Wt I-I and Wt U-I CCF2-AM correlations, this type of analysis excludes the impact of paracrine effects because each set of CCF2-AM correlation data is computed using the CCF2-AM ratios of both I-I cells and U-U cells from the same infected monolayer. Totals of 224 and 343 CCF2-AM ratio correlation deviants were identified between Wt I-I and Dmyr1 I-I cells at 2 hpi and 3 hpi, respectively. Totals of 218 (ϳ97.3%) and 297 (ϳ86.6%) deviants at 2 hpi and 3 hpi, respectively, fell in the regions of the CCF2-AM correlation scatterplot in which the impacts of effectors found in Dmyr1 I-I cells (i.e., ROPs plus MIGs) are dampened by effectors found in Wt I-I cells (i.e., MDGs) (Fig. 7B, blue regions; Table S4B), which implies that signatures induced by effectors released into Dmyr1 I-I cells are counteracted specifically by MYR-dependent effectors released into Wt I-I cells.
(iii) MYR-dependent paracrine effects also counterbalance parasite effectors released earlier in the lytic cycle. To determine the impact of specifically paracrine factors on the host response, we examined two key comparisons: Wt U-U versus mock  (Table S3I) and Dmyr1 U-U versus mock (which captures MYR-independent paracrine effects) ( Table S3J). As illustrated in Fig. 1C and D, differences between these two comparisons should be attributable to paracrine factors released from host cells in a MYR-dependent fashion. At both 2 hpi and 3 hpi, paracrine effects in Δmyr1 parasite infections exhibited robust differences from paracrine effects in wild-type infections. At 2 hpi, ϳ196 and 54 DEGs were identified for the Dmyr1 and Wt cell comparisons, respectively; at 3 hpi, the difference was even more pronounced, with 1,864 DEGs identified for the Dmyr1 comparison, versus only 20 DEGs for the Wt. In addition, GSEA of Dmyr1 U-U versus mock DEGs exposed an abundance of gene sets, many of which corresponded to inflammatory processes and other pathways found to be part of the infection response, whereas GSEA of Wt U-U versus mock DEGs enriched for few, if any, gene sets (3 and 0 at 2 hpi and 3 hpi, respectively) (Fig. 8A). Furthermore, the responses of these DEGs were reproduced consistently not only between U-U cells and mock-infected cells but also between U-I or I-I cells and mock-infected cells (Fig. 8B), which establishes that these trends affect all cell types within a given infected monolayer. Taken together, these results suggest that gene expression trends induced during Toxoplasma infection (i.e., those encapsulated by the DEGs that arise from comparing Dmyr1 U-U and mock cells), including those corresponding to pathways that respond to injected ROPs, are suppressed via a MYR-dependent paracrine mechanism. Accordingly, supernatants taken from host cell cultures infected with wild-type and Δmyr1 parasites could, in theory, exert a transcriptional influence on fresh host cell monolayers, an interesting avenue for future investigation.
Model of host responses to individual Toxoplasma parasite-dependent stimuli and effector compartments. The above-described analyses have accounted for host responses to 5 parasite-dependent stimuli: (i) ROP injection, (ii) MIG secretion, (iii) MDG secretion, (iv) paracrine effects (which can be further subdivided into MYR-independent versus MYR-dependent paracrine effects), and (v) parasite invasion. To succinctly represent the interplay between these stimuli, we curated a list of gene sets that best represented the expression trends captured in our analyses. Gene sets were included if they were significantly enriched (false discovery rate of Ͻ0.25) from DEGs in a majority of the following 8 key comparisons: (i) Wt U-I versus Wt U-U (ROP injection), (ii) Wt I-I versus Dmyr1 I-I (MDGs plus MYR-dependent paracrine effects), (iii) Wt U-U versus Dmyr1 U-U (MYR-dependent paracrine effects), (iv) Dmyr1 I-I versus Dmyr1 U-I (MIGs plus parasite invasion), (v) Dmyr1 U-U versus mock (MYR-independent paracrine effects), (vi) Wt I-I versus Wt U-I (MIGs , MDGs, and parasite invasion), (vii) Wt I-I versus Wt U-U (all stimuli except paracrine effects), and (viii) Wt I-I versus mock (all stimuli). The 12 gene sets selected included those pertaining to immune responses, cell proliferation, cellular stress, and the complement pathway. Expression trends for the DEGs within these gene sets are summarized in Fig. 9.
Based on our analyses, we propose the following model of parasite-driven effects on host cell transcription. First, because Wt U-I cells appeared to induce DEGs in nearly all the showcased gene sets compared to Wt U-U cells, ROP injection likely induces many of the pathways that these gene sets capture, particularly immune-related and cellular stress pathways ( Fig. 9A and B, ROPs, and Fig. 9C, ROP injection). These signatures may arise in response to the ROPs themselves or may be due to cell trauma secondary to perforation of the host cell membrane by the parasite during ROP injection. Next, parasites penetrate the host cell and secrete MIGs into the PVM. The host response to these two stimuli together is relatively mild, as evidenced by the high false discovery rates and small number of genes per gene set, even at 3 hpi. Although MIGs and parasites appeared to counteract the effects of ROP injection at the gene level, their most significant impact at the gene set level appears to be on the enhancement of inflammatory transcriptional signatures induced by ROP injection (Fig. 9A and B, MIGs plus invasion, and Fig. 9C, MIG secretion and parasite penetration). This does not exclude the possibility that other host processes not covered by the Hallmark gene sets are suppressed by MIGs and parasites. Next, in cells infected with Δmyr1 parasites, MYR-independent paracrine factors secreted by neighboring infected cells appear to enhance the effects of ROP injection even more so than MIGs plus parasites and may do so not only for immune-related genes but also for genes pertaining to cellular stress, complement, and cellular proliferation [ Fig. 9A and B, P(MI)]. In contrast, during wild-type infections, two additional parasite-dependent stimuli, MDGs and MYRdependent paracrine factors, both appear to rein in transcriptional signatures induced by the other stimuli [ Fig. 7B; Fig. 9A   Impact of Toxoplasma ROPs/GRAs on Host Transcription ® toward the induction of genes pertaining to inflammation and cellular stress but that are less pronounced than the response to the injection of ROPs.

DISCUSSION
In this study, we examined host responses to infection with the parasite Toxoplasma gondii using scRNA-seq on in vitro-infected, uninfected, and uninfected-injected (U-I) host cells, the latter of which arise from aborted invasion events and have until very recently (44) been characterized primarily morphologically. Key fixtures of our experimental pipeline included (i) early time points, to limit isolating false-positive U-I cells arising from mechanisms besides aborted invasion; (ii) FACS, which purified rare U-I cells and relatively rare infected cells at early time points; and (iii) single-cell resolution, which enabled bioinformatic validation of all cells' infection status. The level of confidence lent by these measures to the validity of the captured U-I cells enabled interrogation of host responses specifically to ROP injection, an aspect of parasite infection previously inaccessible for study due to the rapid kinetics of effector secretion at the time of invasion. Note that while others have also used scRNA-seq to measure host responses to Toxoplasma infection at a single-cell level, those studies analyzed cells from animals at many days postinfection and therefore did not assess the earliest impacts of infection or particularly ROP injection (45).
Because the experimental pipeline also leveraged infections with Δmyr1 parasites, our data set is a comprehensive resource for the individual impacts of not only ROPs but also MIGs, MDGs, and paracrine stimuli on host transcription. Our analyses revealed an early response to Toxoplasma infection with subtle yet clear signatures overlapping inflammatory and cellular stress signaling axes. The induction of these axes appears to be (i) driven primarily by ROP injection; (ii) enhanced somewhat by the combination of MIGs, parasites, and MYR-independent paracrine factors; and (iii) counterbalanced by MDGs and downstream paracrine effects, i.e., factors secreted during wild-type but not Δmyr1 infection. These findings substantiate previous evidence that at least some MDGs suppress host responses induced by other parasite-driven transcriptomic perturbations (31). They may also explain the recent finding that the avirulent phenotype of Δmyr1 parasites during in vivo mouse infections is rescuable by coinfecting animals with both wild-type and MYR1-deficient parasites (46); parasites expressing MYR1 may induce host cells to secrete paracrine factors that suppress the transcription of inflammatory gene products that would otherwise limit Δmyr1 parasite infections. Note that while this article was in preparation, a transcriptomic study of U-I macrophages was reported by Chen and colleagues (44). While they used different strains (Pru and CEP), looked at later time points (20 to 24 hpi), and did not look at MYR1-dependent effects, their primary conclusions that U-I cells experience a major impact of rhoptry effectors (in their case, specifically ROP16) and that paracrine effects are also in play are similar to the conclusions reached here.
Of note, the host response to MIGs plus parasites was especially subtle, given that parasite invasion involves dramatic mechanical perturbations to host cells that might be expected to trigger transcriptional responses. As we did not use host cells containing MIGs and no parasites or vice versa, we could not discern the impact of the MIGs in isolation. Nonetheless, it is tempting to speculate that the response to MIGs plus parasites may reflect MIGs counteracting the effects associated with parasite penetration and PVM formation, which may explain the subtle net inflammatory response.
In addition, the presence of ROPs and MDGs that respectively activate and suppress certain host processes might be interpreted as energetically wasteful. Why might Toxoplasma, an obligate intracellular organism likely under selective pressure for transcriptional efficiency, expend extra resources on effectors that negate one another's effects? One possibility is that ROPs may have undergone selection optimizing for functions required to establish and protect the parasite's intracellular niche (4)(5)(6)(7)(8)(9)(47)(48)(49)(50)(51)(52) but that may also trigger unavoidable host responses detrimental to the parasite. In this scenario, MDGs and paracrine effects could ameliorate such ROP-triggered side effects while theoretically leaving processes beneficial to the parasite intact. Another possibility is that ROPs and GRAs could grant Toxoplasma the ability to fine-tune host responses in terms of timing and/or magnitude, likely an advantage to an organism that must be equipped to encounter a diversity of intracellular environments due to its extraordinarily broad host range (1,53,54).
Our analysis also reveals a striking effect on the host cell cycle, in which U-I cells exhibit an enrichment of G 2 /M-phase cells. Curiously, this enrichment appeared not to be preserved at 3 hpi, which suggests that the responsible parasite factors exert only a transient influence on host cell cycle-related genes. Of the parasite effectors currently known to modulate the host cell cycle (19, 22, 23, 40, 55), ROP16 is most likely to explain these results: ROP16 phosphorylates ubiquitin-like containing PHD and RING finger domain 1 (UHRF1) in a manner that peaks at 3 hpi, which leads to epigenetic silencing of cyclin B1 (40), a component of the cyclin B1/Cdk1 complex required for the G 2 /M transition. However, ROP16 is likely not the only ROP to impinge on the host cell cycle, as a comparison of G 1 -phase U-I versus U-U cells also revealed enrichment in the "p53 Pathway" and "KRAS Signaling Up" gene sets ( Fig. 4 and Fig. 9A and B, ROPs), whose corresponding pathways promote cell cycle progression. Further analysis of U-I cells from yet more time points and particularly those not in G 1 may shed more light on this interplay.
Finally, our analyses thus far reflect a fraction of the possible uses of our data set, which includes variables of time, parasite strain, infection status, cell cycle phases, and UMAP populations. For example, our data set includes reads from not only host cells but also the Toxoplasma parasites, rendering the data a cotranscriptomic resource that will likely illuminate novel host-parasite interactions. Furthermore, because of limiting numbers, the analyses described here dealt primarily with G 1 -phase cells in UMAP population 1 and largely excluded cells in the remaining cell cycle phases and UMAP populations. Future analyses of such cells may reveal roles for the cell cycle phase, host cell type, or other host-dependent processes.

MATERIALS AND METHODS
Cell and parasite culture. All Toxoplasma gondii strains were maintained by serial passage in human foreskin fibroblasts (HFFs) cultured at 37°C in 5% carbon dioxide (CO 2 ) in complete Dulbecco's modified Eagle medium (cDMEM) supplemented with 10% heat-inactivated fetal bovine serum (FBS), 2 mM L-glutamine, 100 U/ml penicillin, and 100 g/ml streptomycin.
A CRISPR-Cas9 strategy was used to construct RH Δmyr1 mCherry Tfn-BLA parasites from an RH Δmyr1 mCherry parental strain (25). The parental strain was transfected with the plasmid pSAG1::Cas9-U6::sgUPRT and a linear construct that contained toxofilin-hemagglutinin-beta-lactamase (Tfn-HA-BLA) expressed under the control of the toxofilin gene's endogenous promoter. The Tfn-HA-BLA construct was PCR amplified from the plasmid SP3 (8) such that the final amplicon was flanked by 20-nucleotide (nt) homology arms to the UPRT gene that are identical to those used to target constructs to the UPRT locus previously (56). Fifteen micrograms of pSAG1::Cas9-U6::sgUPRT and 3 g of the Tfn-HA-BLA linear amplicon were transfected into ϳ10 7 parental strain parasites using the Amaxa Nucleofector 4D system as described above, except for the setting on the Nucleofector (T-cell human unstimulated HE setting). Selection for parasites with 5 M 5-fluorodeoxyuridine (FUDR) began at 1 day posttransfection and proceeded for 3 lytic cycles in monolayers of human HFFs. The transfected, selected parasite populations were subjected to two rounds of single cloning, one to generate populations of parasites enriched for the presence of the Tfn-HA-BLA construct and one to purify for individual parasite clones containing the construct, where the readout for the presence of the construct was cleavage of the BLA-cleavable fluorescence resonance energy transfer (FRET)-based dye CCF2-AM, which results in a fluorescence color shift.
An IFA of monolayers infected with newly constructed RH Δmyr1 mCherry Tfn-BLA parasites was used to validate the correct localization of the Tfn-HA-BLA construct to the parasite rhoptry organelles and the absence of expression of host factors known to be induced in the host nucleus as a result of MYR-dependent GRAs. Coverslips seeded with confluent HFF monolayers were infected with putative RH Δmyr1 mCherry Tfn-BLA clones for ϳ24 h, fixed in 4% formaldehyde for 20 min, permeabilized in 0.2% Triton X-100 for 20 min, blocked in 1ϫ phosphate-buffered saline (PBS) containing 3% bovine serum albumin (BSA) (3% BSA solution) for 1 h at room temperature, stained with primary and secondary antibodies in the 3% BSA solution, mounted onto glass slides using DAPI (4=,6-diamidino-2phenylindole)-containing VectaShield (catalog number H-1200; Vector Laboratories), and sealed with colorless nail polish. For Tfn-HA-BLA localization, primary antibodies were mouse anti-ROP2/3/4 (1:250 dilution) (clone 4A7) (58) and rat anti-HA (1:500) (clone 3F10, catalog number 11867431001; Sigma-Aldrich), and secondary antibodies were goat anti-mouse IgG-Alexa Fluor 647 (1:2,000 dilution) (catalog number A-21235; Thermo Fisher) and goat anti-rat IgG-Alexa Fluor 488 (1:2,000 dilution) (catalog number A-11006; Thermo Fisher). To verify the absence of the MYR1 protein, the primary antibody rabbit anti-c-myc (1:600 dilution) (catalog number M5546; Sigma) and the secondary antibody goat anti-rabbit IgG-Alexa Fluor 594 (1:2,000 dilution) (catalog number A-11012; Thermo Fisher) were used, and coverslips were seeded with either Ͼ3-week-old confluent HFF monolayers or younger monolayers that had been serum starved for at least 24 h (using 0.5% FBS instead of 10% FBS in the culture medium) to ensure that there was no spurious induction of c-myc in uninfected host cells. Only clones where the anti-HA and anti-ROP2/3/4 signal colocalized and where no expression of c-myc was detected in the infected host nucleus were selected. RH Tfn-HA-BLA TdTomato parasites (32) and the parental RH Δmyr1 mCherry parasites were also subjected to both protocols as positive and negative controls, respectively. Coverslips were imaged on either the Stanford Neuroscience Imaging Service core's Zeiss LSM 710 confocal microscope or an Olympus BX60 upright fluorescence microscope.
Western blotting was used to verify the correct size of the Tfn-HA-BLA fusion protein expressed in RH Δmyr1 mCherry Tfn-BLA parasite clones. Lysates were generated by the treatment of parasite pellets with SDS-PAGE loading dye containing 10% beta-mercaptoethanol. The lysates were separated by SDS-PAGE and transferred to a polyvinylidene difluoride (PVDF) membrane, and the membrane was blocked with TBST (Tris-buffered saline, 0.05% Tween 20) containing 5% milk for 1.25 h. The membrane was incubated with horseradish peroxidase (HRP)-conjugated rat anti-HA (clone 3F10) monoclonal antibodies (Roche, Indianapolis, IN) at a dilution of 1:5,000 for 2 h and with a 1:10,000 dilution of rabbit anti-SAG1 for 30 min followed by a 1:20,000 dilution of goat anti-rabbit-HRP for 30 min and developed using the ECL Prime Western blot system (catalog number RPN2232; Sigma-Aldrich).
Determination of time to division for infected 10 T1/2 cells. To determine the time postinfection at which infected 10 T1/2 cells divided, 10 T1/2 host cell monolayers were infected with RH Tfn-BLA TdTomato parasites, and the infected monolayers were imaged by time-lapse microscopy.
To prepare the parasites for infection, the parasites were released from heavily infected monolayers of HFFs by mechanical disruption of the monolayers using sterile, disposable scrapers and passage at least 6 times through a 25-gauge syringe. The parasites were washed by pelleting out HFF debris (133.5 to 208.5 ϫ g for 5 min) and resuspending the parasite pellet generated by spinning the remaining supernatant at 469.2 ϫ g in phenol red-negative low-serum DMEM.
Time-lapse, epifluorescence images of the infected monolayers were acquired over 16 h in a controlled (37°C and 5% CO 2 ) environment using a Nikon Eclipse inverted microscope (Julie Theriot laboratory). Images were acquired every 20 min using 100 ms of exposure at 25% power for the phase, mCherry (to visualize parasites), and green fluorescent protein (GFP) (to visualize host cytoplasmic CTG-CMFDA) channels. Cells that were uninfected at the start of the time-lapse were monitored to determine whether they were infected by the end of the time-lapse, as determined by (i) the clearing of cytoplasmic CTG-CMFDA in the exact position of the parasite, (ii) the disappearance of birefringence in the phase channel upon parasite invasion, and (iii) the parasite tracking with the cell at all time points following presumed infection. For each cell for which the precise moment of infection was captured in the live-video footage, the time to division was determined using the time when the cell was infected as the start time and the time when the cell divided (if applicable) as the end time.

FACS of single uninfected-injected and control cells for single-cell RNA sequencing. (i) Preparation of single cells for FACS.
To generate uninfected-injected (U-I) cells for FACS, 10 T1/2 host cells approximately 1 week from their date of thaw were seeded into 6-well tissue culture plates at a density of approximately 2.6 ϫ 10 5 cells per well, incubated for Ͼ2 h in cDMEM, and serum starved by incubation in low-serum DMEM (i.e., cDMEM containing 0.5% FBS instead of 10% FBS) for 24 h. 10 T1/2 host cells were chosen as they are monolayer-forming, contact-inhibited fibroblasts that are also suitable for single-cell sorting due to adequate dissociation into individual cells by a combination of mechanical and chemical disruption. Next, either RH Tfn-BLA TdTomato or RH Δmyr1 mCherry Tfn-BLA parasites were released from heavily infected monolayers of HFFs by mechanical disruption of the monolayers using sterile, disposable scrapers and passage at least 6 times through a 25-gauge syringe. A parasite-free lysate was similarly generated by mechanical disruption of uninfected HFFs. Parasites and the parasitefree lysate were washed by pelleting out HFF debris (133.5 to 208.5 ϫ g for 5 min) and resuspending the parasite pellet generated by spinning the remaining supernatant at 469.2 ϫ g in low-serum DMEM containing either 1 M dimethyl sulfoxide (DMSO) (for wild-type and Δmyr1 conditions) or 1 M the invasion inhibitor cytochalasin D (CytD) (for CytD-treated wild-type conditions), incubated at room temperature for 10 min, and applied to the serum-starved 10 T1/2 monolayers at an MOI of 6, which maximized the abundance of U-I cells in tissue culture. All 10 T1/2 monolayers were spun at 469.2 ϫ g for 5 min to synchronize parasite contact with the monolayer. Infections were allowed to proceed for 30 min, 1.5 h, or 2.5 h at 37°C. Because CytD is a reversible inhibitor, extra CytD-containing low-serum DMEM was added to each of the host monolayers infected with drug-treated parasites so as to maintain the concentration of CytD at 0.5 to 1 M for the entire infection duration.
To identify host cells injected by parasite proteins, a 6ϫ stock solution of the beta-lactamase substrate CCF2-AM (catalog number K1032; Thermo Fisher) was added to the medium over 10 T1/2 cells so that the final concentration of CCF2-AM was 1ϫ. CCF2-AM-treated monolayers were incubated under foil (to protect them from light) for 30 min at room temperature (to prevent the breakdown of CCF2-AM, which degrades at 37°C), bringing up the total infection durations to 1 h, 2 h, and 3 h.
To harvest the 10 T1/2 monolayers for subsequent FACS analysis, the monolayers were (i) washed 3 times in 1ϫ PBS to remove any HFF debris adhering to the monolayer, (ii) incubated in trypsin (prepared in plastic vessels only) at room temperature for 6 to 10 min, (iii) quenched in an equal volume of FACS buffer (1ϫ PBS, 2% FBS, 50 mM MgCl 2 ·6H 2 O, 50 g/ml DNase I), (iv) passed 3 times through an 18-gauge syringe to break any residual cell clumps, and (v) washed in FACS buffer to remove excess trypsin. Of note, DNase I and MgCl 2 were included in the FACS buffer to prevent clumping of cells from cell death. The cells were then stained with a viability dye and an extracellular parasite stain by (i) resuspending them in 500 l of chilled 4°C 1ϫ PBS containing 3% BSA and a 1:500 dilution of rabbit anti-SAG2A primary antibody (gift of C. Lekutis) for 30 min on wet ice, (ii) washing them in 5 ml of ice-cold 1ϫ PBS and spinning them at 133.5 ϫ g (lowest setting) at 4°C for 5 min, and (iii) resuspending the pellets in chilled 1ϫ PBS containing 3% BSA, a 1:1,000 dilution of goat anti-rabbit IgG-Alexa Fluor 647 (Thermo Fisher), and 3 l/ml near-infrared live/dead fixable viability dye (catalog number L94375; Thermo Fisher) for 30 min on wet ice. Samples were then washed as described above, the pellets were resuspended in 1 ml of chilled FACS buffer, and the cell suspension was transferred through a nylon filter cap (catalog number 4620F40; Thomas Scientific) into polypropylene FACS tubes stored on wet ice and protected from light until FACS sorting.
(ii) FACS of single cells. To prepare the multiwell lysis plates into which cells were deposited during FACS, lysis buffer was dispensed either by the Mantis liquid-handling robot (Formulatrix) at 0.4 l per well into 384-well hard-shell low-profile PCR plates (Bio-Rad) for single-cell RNA sequencing (scRNA-seq) or by hand at 5 l per well into 96-well hard-shell low-profile PCR plates (Bio-Rad) for bulk RNA sequencing. Lysis buffer was prepared in batches of 8 ml by mixing 5.888 ml of water, 160 l of a recombinant RNase inhibitor (TaKaRa Clontech), 1.6 ml of 10 mM deoxynucleoside triphosphate (dNTP) (Thermo Fisher), 160 l of 100 M oligo(dT) (iDT), 1:600,000-diluted ERCC spike-in RNA molecules (Thermo Fisher), and 32 l of 10% Triton X-100. All reagents were declared RNase free. Lysis plates were prepared the night before each FACS sort, stored overnight at Ϫ80°C, and kept on dry ice during the FACS sort.
All host cell samples were sorted at the SSFF by the BD Influx special-order sorter using the following channels: forward scatter (FSC) (488-nm blue laser, side-scatter [SSC] detector), side scatter (488-nm blue laser, FSC detector), Brilliant Violet 421 (BV421) (405-nm violet laser, V460 detector, which detected cleaved CCF2-AM), Brilliant Violet 510 (BV510) (405-nm violet laser, V520 detector, which detected uncleaved CCF2-AM), mCherry (561-nm yellow laser, Y610 detector, which detected parasite-associated cells), allophycocyanin (APC) (640-nm red laser, R670 detector, which detected the extracellular parasite stain versus anti-SAG2A), and APC-Cy7 (640-nm red laser, R750 detector, which detected dead cells). The gating strategy used to obtain U-I, I-I, and U-U cells is indicated in Fig. 2B. More specifically, cells without red fluorescence (i.e., parasite free) but with an enhanced signal from cleaved CCF2-AM (i.e., injected) were sorted as U-I cells, while those with red fluorescence (i.e., parasite associated), an enhanced signal from CCF2-AM (i.e., injected), and low extracellular parasite staining were sorted as I-I cells. Of note, the parasite-associated gate, from which I-I cells were obtained, was intentionally kept narrow to ensure that host cells were each infected with approximately 1 parasite apiece, which limited confounding downstream analysis with the penetration of Ͼ1 parasite at two different time points. Single-color and colorless controls were used for compensation and adjustment of channel voltages. Fluorescence data were collected with FACSDiva software and analyzed with FlowJo software. Cells were index sorted such that each cell's fluorescence data were recorded for subsequent analysis. For single-cell experiments, cells were sorted into 384-well lysis plates at 1 cell per well. For bulk experiments, cells were sorted into 96-well lysis plates at 50 to 100 cells per well. Plates were sealed with foil plate sealers and immediately placed on dry ice until the completion of the sort. Plates were then stored at Ϫ80°C until library preparation.
cDNA synthesis from single-cell RNA, library preparation, and sequencing. To convert the RNA obtained from single cells and bulk samples to cDNA, we employed the Smart-seq2 protocol (59). For single-cell library preparation, the liquid-handling robots Mantis (Formulatrix) and Mosquito (TTP Labtech) were employed to transfer and dispense small volumes of reagents, and the final reaction mixture volume was 2 l per well. For bulk sample library preparation, liquid handling was performed with standard multichannel pipettes, and the final reaction volume was 25 l per sample. cDNA was subjected to 19 rounds of preamplification, quantified using EvaGreen, and diluted in EB buffer to obtain a final concentration of 0.4 to 0.8 ng/l per sample. Library preparation continued using in-house Tn5 tagmentation. For single-cell libraries, we used custom-barcoded indices for each cell, and for bulk libraries, we used Nextera XT indices. Libraries were submitted to the Chan Zuckerberg Biohub Genomics Core for sequencing. Single-cell libraries were sequenced on the NovaSeq 6000 system by 2-by 150-bp paired-end sequencing aiming at ϳ1 million reads per cell. Bulk libraries were sequenced on the NextSeq system by 2-by 150-bp paired-end sequencing at ϳ10 million reads per sample.
Sequence alignment. Read outputs from sequencing were aligned to a concatenated genome composed of the mouse genome (GRCm build 38) and the Toxoplasma gondii GT1 genome (ToxoDB version 36), which is the most complete reference for type I parasite strains such as the RH strains used in this work. Alignment was performed using STAR, and transcript counting was performed using HTSeq-count, with standard parameters used for both packages. A custom Python script was used to sum transcript counts to yield a final gene count matrix consisting of all sequenced cells and the number of reads detected for each gene.
Data preprocessing. To filter out cells of poor quality from the analysis, we excluded cells based on the following metrics: gene count, total read sum, percentage of reads that mapped to the mouse-Toxoplasma concatenated genome, percentage of reads derived from spiked-in ERCC standards, and percentage of reads derived from rRNA.
The gene count matrices were then normalized as counts per median (cpm). Briefly, we first calculated the sum of reads for all cells. We then divided the read counts by the corresponding sum of reads in each cell and multiplied the fractional count by the median of the sum of reads as a scaling factor. Normalized data were transformed to log 2 space after adding a pseudocount of 1 for each gene of each cell.
To determine the detection limit of each experimental trial (e.g., the 50% detection rate), we computed a logistic regression model from a plot of the detection probability for spiked-in ERCC standards. We then excluded genes where Ͻ5 cells in the experimental trial expressed that gene at a level above the detection limit.
To identify host genes associated with infection, we first excluded mouse genes to which Toxoplasma sequences erroneously map (see Table S1 in the supplemental material) by aligning RNA sequences obtained from single-cell extracellular RH parasites to the concatenated mouse-Toxoplasma genome and eliminating all mouse genes with an average log 2 (cpm ϩ 1) expression value of 0.2 or higher.
Single-cell versus bulk sample correlation analysis. To validate the single-cell expression data, we plotted the log 2 mean expression value of each differentially expressed gene, as identified by the MAST algorithm using the single-cell data, in single cells versus bulk samples. The Sklearn package was applied to compute linear regression and the corresponding R 2 (coefficient of determination) values.
Cell cycle analysis and annotation. To predict the cell cycle phase of individual single cells, we curated a list of 175 cell cycle marker genes from the literature (60-67) and from the CycleBase 3.0 database (68). We computed the first two principal components of the gene count matrix using principal-component analysis (PCA) and projected the cells. We partitioned the cells using K-means clustering and assigned the clusters with their predicted cell cycle phases (G 1 , S, and G 2 /M) based on the expression of 175 cell cycle marker genes curated from the literature (Table S2).
Dimensionality reduction. To visually represent relationships between single cells across all experimental trials based on their transcriptional variation regardless of the experimental conditions, we identified and filtered for the top 1,000 genes with the highest dispersion (i.e., the genes with the most variable expression for their bin groups with similar expression levels) across all data sets, applied mutual nearest-neighborhood batch correction (MNNPY) to correct for batch effects, and projected the data onto two-dimensional space with the uniform manifold approximation and projection (UMAP) algorithm using default parameters in Scanpy (69). Leiden clustering using the top 1,000 dispersed genes enabled partitioning of cells into three populations (1, 2, and 3), which we separated into individual data sets for downstream analysis.
Infection status classification. We determined the host infection load by quantifying the percentage of reads that mapped to Toxoplasma in a given sample. We filtered samples that were shown by FACS to exhibit one presumed infection status (based on red fluorescence from internalized parasites) but were determined to exhibit the opposite or an ambiguous infection status otherwise (based on the percentage of reads derived from Toxoplasma).
Differential expression analysis. To obtain differentially expressed genes (DEGs) between all pairs of conditions, we used the model-based analysis of single-cell transcriptomics (MAST) algorithm (70) to compute the results on all G 1 , correctly classified, and UMAP population 1 cells. We used the default settings, except with an adaptive conditional mean of expression based on 20 bins, with at least 30 genes in each bin, and we did not filter out any gene with a nonzero expression frequency in the samples.
Gene set enrichment analysis. We performed gene set enrichment analysis (GSEA) on the lists of differentially expressed genes between all pairs of experimental conditions, ranked by their relative expression under each of the two conditions, for each experimental trial using the fast preranked gene set enrichment analysis (fgsea) package (71). Genes were compared to the Molecular Signature Database's Hallmark gene sets. Pathways with an adjusted P value, i.e., a false discovery rate (FDR), of Ͻ0.25 were considered to be significantly enriched at the top or the bottom of the ranked list of differentially expressed genes.
Identification of differentially regulated genes between conditions using CCF2-AM ratio correlation analysis. To identify host genes associated with injection, we first computed the CCF2-AM ratio for each cell, a metric that serves as a readout for parasite effector injection and that was calculated by dividing the log-transformed cleaved CCF2-AM fluorescence by the log-transformed uncleaved CCF2-AM fluorescence. Next, we excluded mouse genes below the detection limit and mouse genes to which Toxoplasma sequences erroneously map. Finally, for each infection or injection condition (i.e., Wt U-I, Wt I-I, Dmyr1 U-I, Dmyr1 I-I, CytD U-I, and CytD I-I), we computed the Spearman correlation of each of the remaining genes to the CCF2-AM ratio using cells from the condition of interest and cells from the cognate U-U condition (e.g., to calculate Spearman correlations for Wt U-I cells, we correlated gene expression to the CCF2-AM ratio in Wt U-I and Wt U-U cells).
To identify genes differentially regulated between pairs of conditions using the CCF2-AM correlation data, we generated scatterplots where each data point represented a gene and its displacement on each of the x and y axes represented the Spearman correlation under each of the conditions. We modeled a Gaussian distribution using the Sklearn package with default parameters with settings "n_compo-nentsϭ1" and "covariance_typeϭ'full'." Genes with a difference in the CCF2-AM correlation score of greater than 0.2 or less than Ϫ0.2 and whose probability densities were more than 3 standard deviations from the mean probability density were considered to exhibit significant differential expression between the pair of conditions under scrutiny.
Generation of strip plots. The strip plots in Fig. 8B were generated using seaborn's catplot and boxplot packages to plot the normalized expression scores for each cell across 11 experimental conditions for DEGs between Dmyr1 U-U and mock samples. The normalized expression score for each cell was calculated by (i) subtracting the minimum log 2 cpm for that gene across all cells in the experiment, (ii) dividing the difference by the maximum log 2 cpm across all cells such that the cell with the lowest count received a score of 0 and the cell with the highest count received a score of 1, and (iii) computing the average of the normalized cpm for each cell.
Generation of bubble plots. Bubble plots in Fig. 9A and B were generated using a custom Python script in which columns indicated the gene set and rows indicated the comparison (taken to signify the impact of one to a few parasite-dependent stimuli on host transcription) from which the gene sets were enriched. The bubbles in each plot were color-coded and sized based on one of two schemes. In the first scheme, bubble color indicated the normalized enrichment score, where scores were positive if the genes corresponding to a given gene set were expressed at a higher level in the first member of the pair in the given comparison than in the second member of the pair, and bubble size indicated the significance (i.e., FDR) of the enrichment, where the absolute size of each bubble corresponded to the reciprocal of the FDR. In the second scheme, bubble color indicated the log 2 -normalized fold change in the expression of DEGs from a given comparison that also corresponded to a given Hallmark pathway, where positive fold changes indicated that genes were expressed at higher levels in the first member of the pair in the comparison than in the second member, and bubble size indicated the number of DEGs used to calculate the fold change. To calculate the fold change for each given comparison's gene set, DEGs from the comparison that fell under the gene set were identified, the expression of these DEGs in cpm was normalized across all cells to a maximum value of 1 and averaged across all the cells under each condition, and the average normalized expression under the first experimental condition was divided by the average normalized expression under the second experimental condition.
Data availability. The RNA sequencing data set produced in this study has been uploaded in its entirety to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database under accession number GSE145800. The data set includes the raw fastq files, processed gene count files (in counts per median), an anndata file containing the processed gene count files and other metadata such as the cell cycle phase and percentage of Toxoplasma-derived reads for each cell, and results of differential gene expression analysis and gene set enrichment analysis for up to 11 distinct species of host cell (i.e., Wt U-I, Wt I-I, Wt U-U, Dmyr1 U-I, Dmyr1 I-I, Dmyr1 U-U, CytD U-I, CytD I-I, CytD U-U, mock, and CytD mock) at each of three time points (i.e., 1, 2, and 3 h) postinfection. Here, Wt, Dmyr1, and CytD designate host cells arising from monolayers infected with wild-type (RH Tfn-BLA TdTomato), Δmyr1 (RH Δmyr1 mCherry Tfn-BLA), and cytochalasin D-treated wild-type parasites, respectively, while mock and CytD mock refer to cells arising from monolayers mock infected with a parasite-free lysate (where the lysate was pretreated with cytochalasin D under the CytD mock condition).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.