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

A novel FACS-based single-cell RNA sequencing pipeline captures transcriptomic signatures of individual effector compartments.According to the current model of parasite effector secretion (Fig. 1A), the host response to infection with Toxoplasma tachyzoites can be attributed to 5 parasite-dependent stimuli: (i) the secretion of paracrine effectors into the extracellular milieu, (ii) injection of ROPs, (iii) activity of MIGs, (iv) secretion of MDGs, and (v) parasite invasion and establishment of the PV. To resolve the individual impacts of each of these classes of stimulus, we transcriptomically profiled infected, uninfected, and U-I host cells from monolayers exposed to parasites. To this end, we devised a novel pipeline in which (i) infections were designed to generate a heterogeneous pool of host cells, each impinged upon by one of several combinations of parasite-dependent stimuli; (ii) a fluorescence-activated cell sorting (FACS)-based protocol was employed to separate host cells of interest from infected monolayers based on the parasite-dependent stimuli by which they were affected; and (iii) the host cell transcriptomes were profiled by full-length scRNA-seq. Given the novelty and complexity of this pipeline, before we delineate our biological findings, here we present an account of our strategy and technical findings during pipeline development.

(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₁, S, or G₂/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₂ 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–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₁ 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 parasite-dependent 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₁-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.

• Open in new tab
• Download powerpoint

FIG 3

Infection response recapitulates previously identified signatures of early infection. (A) Gene set enrichment analysis (GSEA) of the ranked list of differentially expressed genes between Wt I-I and mock cells at 2 hpi (false discovery rate [FDR], <0.1; fold change [FC], >1.5) using the Hallmark gene sets from the Molecular Signatures Database. (B) GSEA as in panel A of the ranked differentially expressed genes between Wt I-I and mock cells at 3 hpi. Gray cells indicate pathways common to those identified with an FDR of <10⁻⁵ by Naor et al. (31). NES, normalized enrichment score; TNFA, tumor necrosis factor alpha; IL6, interleukin-6; TGF, transforming growth factor; UP, upregulated by the stimulus indicated in the gene set name; DN, downregulated by the stimulus indicated in the gene set name.

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.

• Open in new tab
• Download powerpoint

FIG 4

Host response to injected parasite effectors. (A) GSEA of the ranked list of differentially expressed genes (FDR, <0.1; fold change, >1.5) between Wt U-I and Wt U-U cells at 3 hpi using the Hallmark gene sets. Green rows indicate gene sets enriched from genes expressed at higher levels in Wt U-I than in Wt U-U cells, and red rows indicate gene sets enriched in genes expressed at lower levels in Wt U-I than in Wt U-U cells. Gene sets are designated as being part of the infection response if they were also found to be enriched from DEGs between Wt I-I and mock cells at 2 hpi or 3 hpi. (B) GSEA as in panel A between various U-I cells (i.e., under Wt, CytD, and Dmyr1 infection conditions) and their cognate U-U cells at 3 hpi. Gene sets in bold gray text were also found to be enriched in the infection response. Positive normalized enrichment scores represent enrichment in genes expressed at higher levels in U-I cells than in U-U cells.

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 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 parasite-dependent 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).

• Open in new tab
• Download powerpoint

FIG 5

Comparison of host responses to injection versus infection (i.e., signatures in Wt U-I versus Wt I-I cells). (A) Scatterplot comparing injection-associated genes in Wt I-I cells and Wt U-I cells, as determined by the Spearman correlation of gene expression to the CCF2-AM ratio in Wt I-I and Wt U-I cells. Contours reflect the density of points, with central contours being in the area of highest density. Labeled points on the plot are the top 20 genes with the lowest probability of belonging to a Gaussian model fit with parameters that best describe the observed Spearman correlation data. Yellow and blue shading marks regions of the scatterplots where gene expression trends in Wt U-I cells are either enhanced or dampened upon parasite invasion and the release of dense granule effectors into the host cell (as seen in Wt I-I cells). (B) GSEA of the ranked list of differentially expressed genes (false discovery rate [FDR], <0.1; fold change [FC], >1.5) between Wt U-I and Wt I-I cells at 3 hpi, where green rows denote gene sets enriched from genes expressed at higher levels in Wt U-I cells than in Wt I-I cells and red rows represent genes expressed at higher levels in Wt I-I than in Wt U-I cells. NES, normalized enrichment score; PI3K, phosphatidylinositol 3-kinase.

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.

At 3 hpi, 122 DEGs and 423 CCF2-AM correlation deviants were identified for Wt I-I versus Wt U-I cells, pointing to a divergence in the injection versus infection responses at this later time point. Of the DEGs, 100 (∼82.0%) exhibited opposing trajectories in Wt I-I versus Wt U-I cells (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 injection-associated 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 ROP-induced 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 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.

• Open in new tab
• Download powerpoint

FIG 6

GSEA of the ranked list of differentially expressed genes between Dmyr1 U-I and Dmyr1 I-I cells at 3 hpi (false discovery rate [FDR], <0.1; fold change [FC], >1.5), which illustrates the impact of MYR-independent dense granule proteins (MIGs) and parasite invasion on the host response to infection with Toxoplasma. Green rows denote gene sets enriched from genes expressed at higher levels in Dmyr1 U-I cells than in Dmyr1 I-I cells, and red rows represent genes expressed at higher levels in Dmyr1 I-I cells than in Dmyr1 U-I cells. NES, normalized enrichment score.

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

• Open in new tab
• Download powerpoint

FIG 7

Impact of MYR-dependent dense granule proteins (MDGs) and MYR-dependent paracrine effectors, as illustrated by comparing Dmyr1 I-I and Wt I-I cells. (A) GSEA of the ranked list of differentially expressed genes between Dmyr1 I-I and Wt I-I cells at 3 hpi (false discovery rate [FDR], <0.1; fold change [FC], >1.5). Green rows denote enrichment from genes expressed at higher levels in Dmyr1 I-I cells than in Wt I-I cells, and red rows represent enrichment from genes expressed at higher levels in Wt I-I than in Wt U-I cells. NES, normalized enrichment score. (B) Scatterplot comparing injection-associated genes in Dmyr1 I-I cells and Wt I-I cells, as determined by the Spearman correlation of gene expression to the CCF2-AM ratio in Dmyr1 I-I and Wt I-I cells. Contours reflect the density of points, with central contours being in the area of highest density. Labeled points on the plot are the top 20 genes with the lowest probability of belonging to a Gaussian model fit with parameters that best describe the observed Spearman correlation data. Yellow and blue shading mark regions of the scatterplots where gene expression trends in Dmyr1 I-I cells are either enhanced or dampened by MDGs or paracrine effects (as seen in Wt I-I cells).

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 (which captures all paracrine effects) (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.

• Open in new tab
• Download powerpoint

FIG 8

Comparison of Dmyr1 U-U versus mock cell types reveals the impact of MYR-dependent paracrine effects on host cell transcription during Toxoplasma infection. (A) GSEA of the ranked list of differentially expressed genes between Dmyr1 U-U and mock cells at 3 hpi (false discovery rate [FDR], <0.1; fold change [FC], >1.5). sc indicates the current single-cell RNA-seq data set. NES, normalized enrichment score. (B) Trends in the expression of Dmyr1 U-U versus mock DEGs are conserved in U-I and I-I cells from Wt and Dmyr1 infections, as illustrated by strip plots depicting the expression of DEGs between Dmyr1 U-U and mock cells across all cell types. Each point represents a single cell, and its y axis displacement reflects the average of the normalized expression scores for all DEGs between Dmyr1 U-U and mock cells. The normalized expression score of a given gene is calculated by scaling the log₂ counts per median (cpm) for that gene such that the cell with the lowest cpm receives a score of 0 and the cell with the highest cpm receives a score of 1.

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.

• Open in new tab
• Download powerpoint

FIG 9

Summary of the host transcriptional response to parasite-dependent stimuli. (A) Bubble plot of parasite-dependent host transcriptional signatures at 3 hpi that depicts the behavior of 11 representative Hallmark gene sets from the Molecular Signatures Database (columns) across 8 key comparisons between experimental conditions (rows) (indicated on the left of each plot), each of which encapsulates the impact of certain parasite-dependent stimuli on the 11 gene sets (indicated on the right of each plot). Colors indicate normalized enrichment scores from gene set enrichment analysis of differentially expressed genes between each comparison. An FDR of <0.25 is considered significant (as is standard for GSEA). Parasite-dependent stimuli include ROPs (rhoptry proteins), MDGs (MYR-dependent dense granule proteins), MIGs (MYR-independent dense granule proteins), P(MD) (MYR-dependent paracrine effects), P(MI) (MYR-independent paracrine effects), and all stimuli [ROPs, MDGs, MIGs, invasion, P(MD), and P(MI)]. (B) Bubble plot as in panel A where colors indicate log₂-normalized fold changes in the expression of DEGs that fall under each comparison and each gene set. (C) Model of individual impacts of parasite-dependent stimuli on host transcription, which include immune pathways (i.e., TNFA signaling via NF-κB, IL6 JAK STAT3 signaling, the inflammatory response, IL2 STAT5 signaling, the interferon gamma response, and allograft rejection), cellular stress (i.e., hypoxia and apoptosis), and cell proliferation-associated pathways (i.e., KRAS signaling and the p53 pathway). Curved lines indicate modulation of host transcription, and jagged lines indicate modulation without effector translocation to the nucleus.

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 MYR-dependent paracrine factors, both appear to rein in transcriptional signatures induced by the other stimuli [Fig. 7B; Fig. 9A and B, MDGs + P(MD); and Fig. 9C, MYR-dependent stimuli). Together, all five parasite classes produce transcriptional signatures that veer toward the induction of genes pertaining to inflammation and cellular stress but that are less pronounced than the response to the injection of ROPs.