Introduction

The protozoan parasite Trypanosoma cruzi is the etiological agent of Chagas disease, a highly prevalent infectious disease in Central and South America that affects ∼8 million people, with several million more at risk of infection [1].

The life cycle of the parasite involves both an invertebrate vector (triatomine bug) and a vertebrate host. The change in environmental conditions triggers differentiation processes in the parasite developing across four main stages. The epimastigote form replicates within the insect and differentiates into metacyclic trypomastigotes in the insect’s rectal tract. This latter form does not replicate, specializes in host infection, and can invade various cell types. Once inside the cell, the parasite differentiates into an amastigote, the replicative form within the mammalian host. After several rounds of division, amastigotes differentiate into bloodstream trypomastigotes, which, after cell lysis, can either infect new cells or be ingested by the vector [2]. Parasite developmental stages exhibit significant morphological and molecular differences. T. cruzi infection relies on a heterogeneous set of membrane proteins, encoded mainly by large multigene families [3]. Among these are trans-sialidases and trans-sialidase-like proteins (TcS), Mucins, MASPs, GP63, DGFs, and RHS proteins, most of which are involved in infection, tropism, and immune evasion [4].

The trans-sialidase (TcS) superfamily is involved in processes underlying host-parasite interactions [5]. Members that retain enzymatic activity catalyze the transfer of sialyl groups from host glycoconjugates to galactopyranosyl units on the parasite’s surface, an essential activity, as T. cruzi is unable to synthesize de novo sialic acid [6]. This is the largest superfamily, with over 1,400 members, and is subdivided into eight groups based on amino acid sequence [5]. Few members have been functionally characterized, most expressed primarily in mammalian stages [7]. The TcS group I include proteins that retain TS activity, though members from all groups are involved in host-parasite interactions [8].

In recent years, single-cell RNA sequencing (scRNA-seq) has been employed in protozoa, with reports including Plasmodium falciparum, Trypanosoma brucei, and Leishmania major [9-18]. These studies revealed key aspects of the infection process undetectable with conventional methods, highlighting the relevance of this approach for understanding individual variation in gene expression in single-celled organisms [19]. In T. brucei, antigenic variation driven by variant surface glycoproteins (VSGs) has been studied at single-cell resolution to understand the mechanisms that enable subpopulations of this parasite to evade the immune system, as a bet hedging strategy, ensuring parasite survival [14]. In this parasite, scRNA-seq revealed that pre-metacyclic cells express multiple VSG transcripts simultaneously, in contrast to metacyclic forms, which display a protein coat composed of a single VSG type.

While population-level heterogeneity in surface protein expression has been suggested as critical for T. cruzi infection, immune evasion, and persistence [20, 21], this has not been studied at the intra-population level, and the underlying mechanisms remain poorly understood. Here, we present a scRNA-seq analysis of T. cruzi, to understand the heterogeneity in surface protein expression within trypomastigote populations.

Results and discussion

Identification of cell populations

To assess the expression of cell-surface protein-coding genes in Trypanosoma cruzi, we conducted a 10X Chromium Single Cell 3’ assay from a mixed population of amastigotes and cell-derived trypomastigotes, aiming at sequencing the transcriptome of 5000 cells. After low-quality cell filtering and gene expression quantification (see Materials and Methods), we obtained 3192 single-cell transcriptomes with 14321 total genes detected, with a mean of 1088 genes and 2461 UMIs detected per cell. These results were comparable to other scRNA-seq studies done in Trypanosoma brucei and recently reported as a preprint in T. cruzi using 10X Chromium technology [15, 22, 23]. Cell populations (Figure 1a) were defined by identifying cluster-specific gene markers (Figure 1b, Supplementary File 1). Two cell clusters were assigned to trypomastigotes and amastigotes (cluster 0 and 1, respectively). Markers gene expression was consistent with previous bulk RNA-seq data from Dm28c [24] (Figures 1b, 1c and 1d, Supplementary Table 1). In turn, we hypothesize that cluster 2 reflects amastigote-trypomastigote transitioning parasites, as its gene markers are differentially expressed in bulk RNA-seq data, but some are upregulated in amastigotes and others in trypomastigotes.

Figure 1.

Expression pattern of surface protein-coding genes

We analyzed differences in gene expression between single-copy genes and multigene families, between and within the identified cell populations. As expected, single-copy gene expression is high in both amastigotes and trypomastigotes and similar on average between both cell types (Supplementary Figure 1a), while surface protein-coding genes are more expressed in trypomastigotes [25] (Figure 2a).

Figure 2.

Within trypomastigotes, heterogenous expression of surface gene families was observed, with high variation in the number of cells in which each surface gene was detected, as well as the total expression level in each cell (Figure 2b). Even though expression heterogeneity is also observed for single-copy genes (Supplementary Figure 1a and 1b), probably due to the sampling biases that cause gene dropout in 10X Genomics technology, we investigated whether this phenomenon was more pronounced in surface proteins. Therefore, we analyzed differences among single-copy and multigene family genes (together or grouped by multigene family) in trypomastigotes, in terms of the number of cells expressing each of the individual genes of each group (Figure 2c and Supplementary Figure 1c). We additionally calculated their Gini index (Figure 2d and Supplementary Figure 1d) [26] to assess the extent to which the expression levels of a given gene were evenly distributed across individual cells. Ribosomal protein-coding genes were included as a control group for homogeneous expression.

Interestingly, compared to single-copy genes, especially ribosomal genes, multigene family genes showed a greater dispersion regarding the number of cells in which each gene was detected (Figure 2c and Supplementary Figure 1c), and a higher value of its Gini indexes with a lower dispersion (Figure 2d and Supplementary Figure 1d). Both observations indicate a higher expression heterogeneity for surface protein expression in the trypomastigote population. Moreover, when different families were analyzed separately, TcS genes showed a particularly high heterogeneity (Supplementary Figure 1c and 1d).

Expression pattern of the TcS superfamily

When we re-clustered the trypomastigote population based solely on trans-sialidase gene expression, we identified two sub-populations: trypomastigote cluster 0 (“trypo_0” composed of 1186 cells), which over-expressed these genes compared to trypomastigote cluster 1 (“trypo_1” composed of 1015 cells) (Figure 3).

Figure 3.

The two trypomastigote sub-populations segregated by TcS expression, show only slight differences in the expression of other gene categories, such as ribosomal protein-coding genes (FC = 1.03), transporters (FC = 1.09), polymerase-related genes (FC = 1.15), and phosphatases (FC = 1.10) (Supplementary Figure 2). Although there is a tendency for surface protein genes to be more expressed in cluster trypo_0 (FC = 1.38), trans-sialidases displayed the highest fold change between the two trypomastigote subpopulations (FC = 1.53) (Supplementary Figure 2). It is tempting to speculate that this may reflect different infectivity amongst trypomastigote subpopulations, consistent with reports of “broad” and “slender” forms [27]. Recently, while this manuscript was being prepared, a similar approach by Laidlaw et al. described the phenomenon by clustering trypomastigote cells based on the expression of all genes [23]. When applying this strategy to our data, their observation was reproduced, further demonstrating the consistent reproducibility of the scRNA-seq results across studies.

When all cells are considered, most TcS genes are expressed, but in each individual cell only approximately 40 TcS genes were detected (Figure 4d upper panel). Interestingly, we observed that both trypomastigote subpopulations contain a subgroup of TcS genes that are detected in a large portion of cells (>40%) (Figure 4a, Supplementary Table 2), indicating high-level expression at the population level. However, these genes did not account for most of the family’s expression in each individual cell. Even family members detected in only a few cells significantly contribute to the total TcS expression in those cells, being quantified at similar levels as those from the frequently detected group (Figure 4b, 4c and 4d). Gene dropouts might favor random patterns of gene family’s detection in scRNA-seq experiments, particularly affecting genes with low expression. Since genes detected in a high percentage of cells are not more highly expressed than those detected in only a few cells (Figure 4b), we do not believe that this characteristic of scRNA-seq data is responsible for our observation. Actually, at the single-cell level, each detected TcS gene contributes similarly to the overall family expression in those cells, as shown by a small average Gini index for the percentage of the expression of the family in each cell (Figure 4c and 4d). This highlights a key limitation of bulk RNA-seq, as it may wrongly indicate that a set of few genes are highly expressed in each cell, while in fact, TcS expression is evenly distributed among all detected family members, regardless of the number of cells expressing each gene. Indeed, most TcS members detected in more than 40% of cells in this work were found in the top 100 expressed TcS in an independent bulk RNA-seq study (Supplementary Figure 3a).

Figure 4.

When analyzing each trypomastigote subpopulation, no coordinated expression among specific TcS members was observed, as no subclusters of cells were identified based on TcS detection profiles. Even when clustering was restricted to genes detected in more than 40% of cells, no clear subclusters of cells were identified (Supplementary Figure 3b).

The loci for TcS genes and other gene families of surface proteins are mostly grouped in specific genomic compartments [28] that are regulated epigenetically by the activation or silencing of chromatin folding domains [24]. The TcS that are detected in a high percentage of cells are mostly dispersed throughout the genome (Supplementary Table 2). This suggests that their preferential expression is likely not due to colocalization in one or a few ubiquitous activated chromatin-folding domains. These observations raise intriguing questions about the mechanisms governing the selection of TcS members expressed in each cell, warranting further investigation.

Final remarks

The expression of surface protein-coding genes varied between cell stages. Genes that belong to the trans-sialidase-like superfamily, involved in host-parasite interactions, displayed significant heterogeneity in expression levels within trypomastigotes. This variation in TcS expression suggests that not all trypomastigotes express these genes uniformly. The lack of coordinated expression among TcS members, observed in single-cell and bulk RNA-seq studies [3], suggests a complex regulatory mechanism that enables adaptability in different T. cruzi subpopulations. These findings reinforce the importance of TcS variability in driving phenotypic diversity and pathogenicity in T. cruzi. Taken together, these findings not only demonstrate the sensitivity of scRNA-seq in delineating parasite life stages and their associated gene expression profiles, but also reveal a complex regulation of surface proteins and TcS family members, with tentative implications for immune evasion and pathogenicity.

Materials and methods

Trypanosoma cruzi and mammalian cell culture

Epimastigote forms of Trypanosoma cruzi strain Dm28c were derived from axenic cultures cultivated in Brain-Heart Infusion medium (BHI, Oxoid) supplemented with 10% heat-inactivated fetal bovine serum (FBS, Capricorn), penicillin (100 units/mL) and streptomycin (100 µg/mL) as described [29]. Cultures were diluted 1/10 with fresh BHI medium every 3 days and maintained at 28 °C.

Myoblast rat cell line H9c2 (ATCC CRL-1446) was maintained in hgDMEM medium (Gibco) supplemented with 10% heat-inactivated fetal bovine serum, penicillin (100 units/mL) and streptomycin (100 μg/mL) at 37 C in a humidified 5% CO2 incubator. Confluent cells were washed with 1X phosphate-buffered saline (1X PBS), incubated for 5 min with trypsin-EDTA (Gibco), diluted with culture medium and re-plated for maintenance.

Mycoplasma contamination in cell lines was regularly monitored using MycoAlert® Mycoplasma Detection Kit (Lonza), following the manufacturer’s protocol.

Isolation and purification of cellular trypomastigotes and intracellular amastigotes

Late stationary phase epimastigotes of Dm28c strain were used to infect H9c2 cells for a primary infection. Six days post-infection, cellular trypomastigotes were obtained from the supernatant and were used to infect 50% confluent H9c2 cells at a 10:1 rate. Twenty-four hours post-infection the cell culture was washed twice with PBS 1X and maintained with fresh hgDMEM at 37 C in a humidified 5% CO2 incubator.

For amastigote purification, H9c2 cells grown in monolayers and infected for 48 hours were washed with 1X PBS and incubated with trypsin-EDTA for 5 minutes at 37°C. The trypsinization was stopped by adding an equal volume of hgDMEM with 10% FBS. The cell suspension was repeatedly passed through a 27-gauge needle attached to a 30-mL syringe until complete cell disruption was confirmed under the microscope. The supernatant, containing free amastigotes, was collected and centrifuged at 500xg for 10 minutes at 4°C to remove large host-cell debris. The resulting supernatant was then centrifuged at 4,000xg for 10 minutes at 4°C, and the amastigote-containing pellet was washed twice in chilled 1X DPBS (Dulbecco’s Phosphate-Buffered Saline, No Calcium, No Magnesium) and resuspended in 1X DPBS at 200 μL per 1×10⁶ cells, ready for the fixation step.

Cellular trypomastigotes derived from infected H9c2 cells and present in the cell supernatant fraction were collected and centrifuged at 500xg for 10 minutes at 4°C to remove large host-cell debris. The washing and resuspension steps in DPBS were performed for amastigotes as described above.

scRNA-seq library preparation and sequencing

Cell fixation was performed using the Methanol Fixation Protocol for Single-Cell RNA Sequencing [30] as recommended by 10X Genomics technical support. Briefly, chilled 100% Methanol (for HPLC, ≥99.9%, Millipore) was added drop by drop (1×10⁶ cells in 800 μl) and incubated at -20°C for 30 min. For rehydration, fixed cells were first equilibrated at 4°C and then centrifuged at 4,000xg for 5 min at 4°C, the supernatant was discarded and Wash-Resuspension Buffer (3X SSC in Nuclease-free Water, 0.04% UltraPure Bovine Serum Albumin, 1mM DTT, and 0.2 U/ml RNase Inhibitor) was added to the pellet. Cell debris and large clumps were eliminated by passing the sample through a 40 μm Flowmi Cell Strainer.

10X Genomics library preparation was performed at the service provided by the Instituto de Biología y Medicina Experimental (IBYME, Argentina). The library was sequenced by a service provider (Macrogen, Korea) in a HiSeq2500 equipment (two lanes), generating approx. 880 million reads. Raw data is available in the SRA (https://www.ncbi.nlm.nih.gov/sra/) under BioProject PRJNA1200704.

Transcript quantification

T. cruzi 2018 Dm28c genome [28] (release 62, TriTrypDB) and T. cruzi Dm28c maxi circle kDNA sequence [31] were combined and used as reference. To improve the proportion of reads assigned to genes, 11,362 3′UTR regions of the coding sequences (CDS) were annotated using peaks2UTR [32]. Gene expression quantification was performed by pseudoalignment using kallisto bustools [33], with the options --filter to remove potential noise from environmental RNA and --em to apply the Expectation-Maximization (EM) algorithm. This algorithm outperforms other mapping software in handling multimapping reads, enabling more accurate quantification of multigene families [33].

scRNA-seq data processing and analysis

Count matrices from two technical replicates obtained by sequencing the same library were merged: the common barcodes across both datasets were retained, and the count matrices were combined to generate a unified dataset. Subsequently, the following metrics were calculated: nUMI, nGene, and mitoRatio, as well as the log10 of genes per UMI (log10GenesperUMI). A filtering criterion was applied, retaining cells with nUMI > 1200, nGene > 100, log10GenesperUMI > 0.8, and mitoRatio < 0.1. Additionally, ribosomal genes were excluded from subsequent analyses. Data normalization and scaling were performed using the Seurat R version 5 package [34], employing NormalizeData and ScaleData functions.

FindNeighbors function was used to construct a k-nearest neighbors graph of cells using 10 principal components (PCs) and Louvain algorithm was employed for cell clustering using FindClusters function.

To define marker genes, the FindAllMarkers function from Seurat package was used, selecting those with an adjusted p-value < 0.05 and a log2 fold change (log2FC) > 1. To validate these markers, bulk RNA-seq data of Dm28c was incorporated (NCBI BioProject ID PRJNA850400 [24]). Transcripts were quantified using Kallisto [35] (with -b 100 option to perform 100 bootstraps), followed by differential expression analysis conducted with Sleuth [36]. Genes with an adjusted p-value < 0.05 and |log2FC| > 0.25 were filtered for further analysis. Finally, the gene IDs of the Seurat-defined markers were cross-referenced with the IDs of the differentially expressed genes obtained from the bulk RNA-seq analysis to corroborate stage-specific gene expression. For 2D visualization of cell clusters and gene markers expression profiles across cells, UMAP projection was employed [37].

Gene IDs corresponding to multigene families (trans-sialidase and trans-sialidase-like, MASP, Mucins, GP63, RHS and DGF) were obtained by text searches using the current genome annotation on the “description” field. Single and low copy number genes were defined as those that did not belong to the latter gene list. For simplicity “single-copy” will be used to refer to these genes throughout the manuscript.

Data processing was conducted using R version 4.2.0. Statistical analyses were performed using the Wilcoxon rank-sum test and p-values < 0.05 were considered statistically significant.

To study expression inequality, we applied the Gini index, a metric originally developed in the field of economics [26]. The Gini index serves as a measure of expression heterogeneity, and we employed it in two distinct contexts: first, to evaluate the degree of inequality in the distribution of a gene’s expression levels across individual cells (Figure 2d and Supplementary Figure 1d); and second, to assess the extent to which a given cell expresses individual genes at varying rates (Figure 4c).

Supplementary figures

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Additional information

Author Contributions

LI: Data Curation, Formal Analysis, Investigation, Methodology, Visualization, Writing-Original Draft Preparation, Writing-Review & Editing; LB: Methodology, Investigation, Writing-Review & Editing; VAC: Methodology, Investigation, Writing-Review & Editing; JG: Methodology, Formal Analysis, Writing-Review & Editing; GR: Methodology, Writing-Review & Editing; VMH: Methodology, Writing-Review & Editing; MAD: Methodology, Resources, Writing-Review & Editing; JSS: Methodology, Resources, Writing-Review & Editing; JDG: Conceptualization, Methodology, Visualization, Writing-Original Draft Preparation, Writing-Review & Editing; PS: Conceptualization, Funding Acquisition, Methodology, Formal Analysis, Project Administration, Resources, Supervision, Validation, Writing-Original Draft Preparation, Writing-Review & Editing

Financial Support

This project was supported by: CSIC, Universidad de la República, grant number: I+D-2020-505 awarded to PS; LI, LB, JG, MD, JSS and PS received financial support from PEDECIBA. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Additional files

Supplementary File 1.

Supplementary Table 1.

Supplementary Table 2.