Single-cell RNA profiling of Plasmodium vivax-infected hepatocytes reveals parasite- and host- specific transcriptomic signatures and therapeutic targets

The resilience of Plasmodium vivax, the most widely-distributed malaria-causing parasite in humans, is attributed to its ability to produce dormant liver forms known as hypnozoites, which can activate weeks, months, or even years after an initial mosquito bite. The factors underlying hypnozoite formation and activation are poorly understood, as is the parasite’s influence on the host hepatocyte. Here, we shed light on transcriptome-wide signatures of both the parasite and the infected host cell by sequencing over 1,000 P. vivax-infected hepatocytes at single-cell resolution. We distinguish between replicating schizonts and hypnozoites at the transcriptional level, identifying key differences in transcripts encoding for RNA-binding proteins associated with cell fate. In infected hepatocytes, we show that genes associated with energy metabolism and antioxidant stress response are upregulated, and those involved in the host immune response downregulated, suggesting both schizonts and hypnozoites alter the host intracellular environment. The transcriptional markers in schizonts, hypnozoites, and infected hepatocytes revealed here pinpoint potential factors underlying dormancy and can inform therapeutic targets against P. vivax liver-stage infection.

Introduction

Malaria is a significant global health burden, with an estimated 421 million cases and 627,000 deaths in 2020 (World Health Organization, 2021). At least five Plasmodium species are known to infect humans, with P. vivax responsible for the majority of cases outside Africa (World Health Organization, 2021). All Plasmodium parasites that infect humans have an obligatory developmental stage in the liver where the parasite undergoes asexual reproduction known as schizogony within a hepatocyte before releasing thousands of merozoites into the blood. In P. vivax, a proportion of parasites known as hypnozoites forgo immediate division in hepatocytes, and persist in the liver for weeks, months, or years before activating to initiate schizogony resulting in a relapse blood infection (Krotoski et al., 1982; Adams and Mueller, 2017). Hypnozoites are thought to cause up to 90% of cases in certain geographical regions (Adekunle et al., 2015; Robinson et al., 2015; Commons et al., 2020), and mathematical models predict that eliminating P. vivax malaria will be difficult without interventions that target hypnozoite reservoirs (White et al., 2018).

The factors influencing the biogenesis, persistence, and activation of hypnozoites are poorly understood. Studies have found that hypnozoites in temperate versus tropical regions have long and short latencies, respectively, suggesting a correlation between geographical region and relapse patterns (Krotoski et al., 1986; Huldén et al., 2008; White, 2011; Battle et al., 2014). These findings suggest the developmental trajectory is genetically encoded and support the hypothesis that it is somehow pre-programmed in the sporozoite (Ungureanu et al., 1976; Lysenko et al., 1977); but other contributing factors from the host could play a role in determining whether the sporozoite becomes a schizont or hypnozoite (Schäfer et al., 2021; Vantaux et al., 2022).

Our understanding of P. vivax liver stage biology at the cellular and molecular level is in its nascency for two overarching reasons. First, protocols to cryopreserve and continuously maintain erythrocytic stage parasites in vitro are lacking, necessitating the collection of gametocyte-infected blood from human volunteers in P. vivax malaria-endemic regions to generate the sporozoites required to initiate an infection. New laboratory models have recently been developed which can partially alleviate the need to obtain gametocytes from human-infected blood (Mikolajczak et al., 2015; Schäfer et al., 2020), yet these sources are expensive, and each have distinct logistical hurdles. In addition to the sporozoite limitation, upon initiating a P. vivax-infected hepatocyte culture, the high ratio of noninfected versus infected hepatocytes has historically made it difficult to isolate and study these forms in sufficient numbers. Robust in vitro pre-erythrocytic platforms capable of long-term culture of sufficient numbers of P. vivax liver stages have only recently been developed (Roth et al., 2018; Maher et al., 2021a). For these reasons, much of our knowledge of hypnozoite biology at the molecular level remains limited to studies using P. cynomolgi—a parasite capable of causing relapsing infections in non-human primates for which molecular tools are available (Cubi et al., 2017; Voorberg-van der Wel et al., 2017; Chua et al., 2019; Voorberg-van der Wel et al., 2020). However, it is unknown whether mechanisms regulating hypnozoite formation and dormancy are consistent across the different relapsing malaria parasites.

Transcriptome-wide profiling of P. vivax liver stages has previously been performed using a bulk RNA sequencing strategy (Gural et al., 2018). While informative, the data provide insights only into average transcript levels derived from either mixed (schizont and hypnozoite) or hypnozoite-only populations, masking any transcriptional variation that might exist between individual parasites. Second, the parasite’s potential impact on the hepatocyte is not considered, barring analysis of host-pathogen interactions. Nevertheless, these shortcomings highlight a need to better understand differences between individual parasites, as well as how the parasite alters the hepatocyte. Advancements in these areas could help elucidate parasite- and host-factors important for the parasite’s liver stage development.

Single-cell RNA-sequencing (scRNA-seq) methods offer new ways of studying cell-to-cell variation (Aldridge and Teichmann, 2020), and have already allowed for novel insights into host-pathogen interactions in various disease contexts (Penaranda and Hung, 2019). In Plasmodium spp., scRNA-seq analyses have so far revealed that individual parasite forms exist on a spectrum of transcriptomic states (Poran et al., 2017; Reid et al., 2018; Howick et al., 2019; Bogale et al., 2021; Real et al., 2021; Ruberto et al., 2021). Their application to study the biology of P. vivax is a burgeoning area of research as only recently have reports described the transcriptomic signatures of sporozoites and blood stages using 10x Genomics’ technology (Sà et al., 2020; Ruberto et al., 2022), and liver forms using Seq-Well S³ technology (Mancio-Silva et al., 2022).

Here, using 10x Genomics’ scRNA-seq platform, we characterize transcriptomic signatures of parasite and host in an in vitro P. vivax liver stage model. We find differences in gene expression between replicating schizonts and hypnozoites, and reveal variation, previously unobservable by bulk sequencing approaches, between individual hypnozoites. In infected hepatocytes, we find that genes associated with energy metabolism and antioxidant stress response are upregulated, and those involved in the host immune response downregulated. Various host genes associated with these processes are found to be exclusively upregulated during hypnozoite infection. Our study elucidates the transcriptional signatures among infected hepatocytes and presents potential mechanisms P. vivax can use to either exploit its host for rapid proliferation or to quietly persist, remaining undetected by the host immune system.

Materials and methods

Blood samples, mosquitos, and infections

The complete protocol used for mosquito rearing and infection has been published (Maher et al., 2021a). In brief, blood samples from symptomatic patients infected with P. vivax were collected at local health facilities in Mondulkiri province (Cambodia) in March 2020 and October 2021. Following a P. vivax gametocyte-containing blood meal, An. dirus mosquitoes were maintained at ambient temperature on a 12:12 light: dark cycle and fed 10% sucrose + 0.5% PABA solution. An. dirus found positive for P. vivax oocysts at six-days post-feeding were transported to the Institut Pasteur Cambodia Insectary Facility in Phnom Penh, Cambodia where they were maintained under the same conditions described above.

Sporozoite isolation and primary human hepatocyte infections

The complete protocol used for seeding primary human hepatocytes, harvesting sporozoites, and infecting cultures has been published (Maher et al., 2021a). In brief, P. vivax sporozoites were isolated from the salivary glands of female An. dirus mosquitoes 18 days after an infectious blood-meal. Primary human hepatocytes (Lot BGW, BioIVT) were seeded 2 days before infection with 30,000 sporozoites (replicate 1), or 15,200 sporozoites (replicate 2) per well. Media was replaced with fresh CP media (BioIVT) containing antibiotics (penicillin, streptomycin, neomycin, and gentamycin) the day after infection and every 2-3 days thereafter. Hepatocytes were treated with 1 µM MMV390048 on days 5, 6, and 7 post-infection to kill replicating schizonts, resulting in cultures enriched with hypnozoites at 9 days post-infection. At days 5 and 9 post-infection, cultures were processed for single-cell partitioning (described in the following section) or immunofluorescence assay (IFA). For the latter, hepatocytes were fixed with 4% paraformaldehyde in 1X PBS. Fixed cultures were then stained overnight at 4°C with recombinant mouse anti-P. vivax UIS4 antibody (Schäfer et al., 2018) diluted 1:25,000 in a stain buffer (0.03% Triton X-100) and 1% (w/v) BSA in 1X PBS. The following day, cultures were washed three times with 1X PBS and then stained overnight at 4°C with rabbit anti-mouse AlexaFluor ™488-conjugated antibody (Thermo Fisher) diluted 1:1,000 in stain buffer. Cultures were washed three times with 1X PBS and counterstained with 1 mg/mL Hoechst 33342 (Thermo Fisher) to detect parasite and host DNA. High content imaging was performed on a Lionheart imaging system (Biotek). Quantification of the number of nuclei and parasite liver forms was performed using Gen5 high content analysis software (Biotek). Confocal images were obtained using an ImageXpress Imaging System (Molecular Devices).

Single-cell partitioning, library preparation, and sequencing

Hepatocyte cultures infected with P. vivax were processed prior to- and post- MMV390048 treatment (5- and 9-days post-infection, respectively). Two replicate experiments, defined as an independent plating of hepatocytes infected with sporozoites from a unique field isolate, were performed several months apart. Additionally, and as a negative control, one replicate of hepatocyte cultures never exposed to P. vivax (naïve) were processed similarly to those that were infected. Naïve cultures processed at the same time point as hepatocytes infected with P. vivax 9 days post-infection also received MMV390048 treatment. At each endpoint for each replicate, approximately 640,000 primary human hepatocytes (16 wells of a 384-well plate each having approximately 10,000 cells) were treated with trypsin (Corning 25-053-Cl) to facilitate their detachment from the plate. Once detached, trypsin was inactivated with 1:1 complete media and the cell suspension was collected into a 5mL protein-low bind centrifuge tube (Eppendorf). To remove residual trypsin, hepatocytes were washed via three sequential rounds of pelleting (50 x g, 3 min, 4°C) and resuspending in 1X PBS containing 0.1% BSA. Cells were then passed through a 35µm strainer (Falcon 352235) to minimize or remove hepatocyte clusters. Prior to partitioning cells on the Chromium controller (10x Genomics), hepatocytes were counted, and the cell concentration was assessed using a hemocytometer. Approximately 10,000 hepatocytes were loaded on a Chromium Chip B (10x Genomics) as per manufacturer’s recommendations. Chips containing hepatocyte suspensions, Gel Beads in Emulsion (GEMs), and reverse transcription reagents were loaded into the Chromium controller (10x Genomics) for single-cell partitioning and cell barcoding. Barcoded cDNA libraries were generated according to the Chromium Single Cell 3’ gene expression protocol (version 3). In total, six cDNA libraries were generated (two replicates at 5 days post-infection, two replicates at 9 days post-infection, one uninfected control for day 5 cultures, and one uninfected control for day 9 cultures). Libraries were loaded on individual flow cell lanes and sequenced at a depth of 400 million reads using a HiSeq X Ten platform (Illumina) at Macrogen (Seoul, Korea).

Computational analysis of single-cell RNA sequencing data

Quality control checks on the raw sequencing data were performed using FASTQC (Andrews, 2010). In a manner similar to other host-pathogen single-cell studies (Ren et al., 2021; Szabo et al., 2021), a multi-species transcriptome index containing the transcriptomes for P. vivax P01 (PlasmoDB.org, version 51) and H. sapiens GRCh8 (Ensembl, version 103) was created using the ‘ref’ function in kb-python (0.26.3)— a python package that wraps the kallisto | bustools single-cell RNA-seq workflow (Bray et al., 2016; Melsted et al., 2021). The P. vivax P01 transcriptome file used for the generation of the multi-species index contained gene regions encoding for UTRs (Siegel et al., 2020). Next, using the same package, reads were pseudoaligned and counted using the ‘count’ function with parameters suitable for sequencing data generated from version 3 of 10x Genomics’ 3’ gene expression technology: kb count -i human_Pv_index.idx -g t2g.txt -x 10xv3 -t8 -o./path/to/read/10xV3sequencesread1.fastq.gz path/to/10xV3sequencesread2.fastq.gz. After performing kb count on the sequencing outputs from the six single-cell libraries, the resulting unfiltered cell-gene count matrices were read into R (v 4.1.1) for further host and parasite-specific processing described in the following sections.

Filtering and normalization of scRNA-seq count matrices

To distinguish between droplets containing cells and those only with ambient RNA, we used DropUtils’ emptyDrops function (Lun et al., 2019) with the lower total UMI count threshold (at or below which all barcodes are assumed to correspond to empty droplets) set to 1,000. The droplets with cells (FDR < 0.001), more specifically, those containing barcoded transcripts derived from either human, parasite, or a combination of the two—were kept for downstream processing.

Assessment of P. vivax liver stage transcriptomes

For the analyses of P. vivax, we filtered out the transcripts derived from human, leaving only transcripts derived from the parasite to be assessed. All samples were processed independently. For each sample, cells with less than 60 genes detected were removed. Next, we filtered out genes with low counts, retaining genes with 10 or greater UMIs across all cells. Post cell and gene filtering, the data were normalized using Seurat (Stuart et al., 2019) ‘LogNormalize’ function with the default parameters selected.

Data integration and clustering

The data were integrated using Seurat’s SCTransform function (Hafemeister and Satija, 2019). This method has been developed to deal with inherent RNA differences in single-cell data and uses a novel method of regularized negative binomial regression that has been demonstrated to prevent sequencing depth and RNA content from biasing downstream analyses. The following parameters indicated: variable.features.n = NULL, variable.features.rv.th = 1.3. Clustering was performed using the Seurat functions FindNeighbours (dims = 1:20) followed by FindClusters (resolution = 0.1, algorithm = 4 for all data or resolution = 0.3, algorithm = 4 for hypnozoite filtered data).

Differential gene expression analysis

To detect differentially expressed genes, the Seurat function FinalAllMarkers and FindMarkers were used. For differential gene expression testing using the FindMarkers function, we used the following parameters: logfc.threshold = 0.0, ident.1 = ‘1’, ident.2 = ‘2’, test.use = “wilcox”, min.pct = 0, only.pos = F, assay = “RNA”. For differential gene expression testing using the FindAllMarkers function, we used the following parameters: test.use = “wilcox”, only.pos = T, assay = “RNA”, log2FC = 0.0. We found negligible differences in outputs between Wilcoxon Rank Sum test and MAST (which uses a hurdle model tailored to scRNA-seq data (Finak et al., 2015)). For the schizont versus hypnozoite comparisons, genes with a Bonferroni adjusted p value < 0.01 & absolute average log2 FC > 0.5 were used RBP and GO analyses. For the hypnozoite versus hypnozoite comparisons, genes with a Bonferroni adjusted p value of < 0.05 and absolute average log2 FC > 0.5 were used for RBP and GO analysis.

RBP analyses

RBPs with cell fating potential reported in humans (Decker and Parker, n.d.; Voronina et al., 2011; Gerstberger et al., 2014; Guallar and Wang, 2014; Youn et al., 2019) were manually curated. The orthologues of these cell fating RBPs in P. vivax were obtained using OrthoMCL (Release-6.8) (Chen et al., 2006). Within the differentially expressed genes, genes encoding for cell fating RBPs were filtered for further assessment.

GO enrichment analyses

GO enrichment analysis was performed on the PlasmoDB (v55) website with the following input parameters selected: P value cutoff = 0.05; Limit to GO Slim terms = No; Evidence: Computed, curated. Enrichment from GO cellular compartment (CC), biological processes (BP), and molecular function (MF) were all assessed. Gene with a Bonferroni adjusted p value < 0.05 were considered significantly enriched.

Assessment of human hepatocyte transcriptomes

For the analyses of hepatocytes, we removed all transcripts derived from P. vivax, leaving only the transcripts derived from hepatocytes. The removal of P. vivax RNA reduced any biases that may arise during the normalization, data reduction, and differential gene expression analyses steps. We retained transcripts detected in at least 30 cells. Next, cells with less than 750 features were removed from further analyses. Post-cell and gene filtering, the data from each replicate were normalized using Seurat’s ‘LogNormalize’ function with the default parameters.

Data integration and low dimensional reduction

Filtered, normalized data matrices were integrated in a manner described in in the Seurat (version 4) vignette, Introduction to scRNA-seq integration, available on the Sajita Lab’s website (https://satijalab.org/seurat/articles/integration_introduction.html). Briefly, integration features were identified in each replicate using the ‘SelectIntegrationFeatures’ function with the nFeatures parameter set to 3000. Next, we used ‘PrepSCTIntegration’ to prepare an object list containing the 3000 features for integration. We then used ‘FindIntegrationAnchors’ to identify a set of anchors between datasets with the following parameters specified: normalization.method = “SCT”, anchor.features = “features”, reference = naïve hepatocytes day 5. Last, using these anchors, the six datasets were integrated using the ‘IntegrateData’ function with the following parameters specified: normalization.method = “SCT”. After integrating the data, we perform dimensionality reduction using the ‘RunPCA’ function, followed by ‘RunTSNE’ with 1:30 dimensions selected.

Differential gene expression analyses

Differential expression between hepatocytes was calculated using the FindMarkers function with the Wilcoxon Rank Sum test implemented in Seurat. Specific parameters implemented for each comparison are indicated in the supplementary tables that the data are reported.

Gene ontology and reactome enrichment analyses

Enriched biological processes and pathways were determined using gProfiler (Raudvere et al., 2019) with the following options selected: ordered query, statistical domain scope—only annotated genes, significance threshold—g:SCS threshold, user threshold—0.05, Gene Ontology—GO biological processes, No electronic GO annotations, biological pathways—Reactome. Transcripts that displayed significant changes in expression in infected hepatocytes versus non-infected hepatocytes were used as input. EnrichmentMap (v 1.1.0) (Merico et al., 2010) for CytoScape (v 3.8.2) (Shannon et al., 2003) was used to create the enrichment networks. Significantly enriched GO biological processes and Reactome pathways (adjusted p value < 0.05) were used as input. Pathways were clustered and annotated using the AutoAnnotate function (Kucera et al., 2016). Node size corresponds to the number of transcripts in the corresponding set, and edge widths correspond to the number of transcripts shared between sets.

Transcription factor enrichment analysis

We downloaded curated sets of known transcription factors from the website (http://amp.pharm.mssm.edu/Enrichr/) (Kuleshov et al., 2016). Transcription factor—target gene sets were obtained by ChIP-X experiments from the ChEA (Lachmann et al., 2010) and ENCODE (Feingold and Pachter, 2004) projects. TF-target relationships using Chip-seq data from ENCODE and CHEA were used to identify overrepresented pathways and TFs. Differentially expressed genes (Bonferroni adjusted p value < 0.05) in each condition versus non-infected cells were used as inputs. Gene sets with an q value < 0.05 were considered significantly enriched.

Carfilzomib dose-response assessment

A carfilzomib dose-response against P. vivax liver forms was obtained from three independent experiments, in which a different P. vivax clinical isolate was allowed to infect BGW human hepatocytes. The 12-day assay was used (Maher et al., 2021a). All independent experiments resulted in sufficient hypnozoite formation within each well for dose-response calculations. Liver stage parasite growth metrics and compound dilutions were loaded into CDD Vault (Collaborative Drug Discovery), and growth data normalized such that zero (negative) inhibition was set as the average of DMSO wells and 100% (positive) inhibition was set to the effective doses of a nigericin control. Normalized values were then combined in Graphpad Prism to fit a dose-response curve and calculated IC₅₀ values from all replicates.

Quantification of AKR1B10 protein expression

To analyze AKR1B10 expression and localization, P. vivax-infected primary human hepatocytes were fixed 12 days post-infection. Five independent hepatocyte infections were used. Infected wells were then permeabilized for 20 minutes with 0.2% TritonX-100, washed thrice with 1X PBS, and blocked in 3% BSA for one hour. Primary antibodies rabbit anti-HsAKR1B10 (PA5-30773, Thermo Fisher) and mouse anti-PvUIS4 were sequentially added, diluted 1:100 and 1:500, respectively. Antibodies were incubated overnight at 4°C and washed three times with 1X PBS before addition of 1:400 Alexa- Fluor 568 donkey anti-rabbit for HsAKR1B10 or 1:400 AlexaFluor 488 goat anti-mouse for PvUIS4. Secondary antibodies were incubated overnight at 4°C. Cells were washed thrice with 1X PBS before adding 0.5 mg/mL DAPI for 10 minutes. Cells were washed three final times with 1X PBS before image collection with an ImageXpress inverted confocal microscope (Molecular Devices). Image analysis was completed using MetaXpress analysis software (Molecular Devices).

Contact for reagent and resource sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Dennis Kyle (email: dennis.kyle@uga.edu).

Results

Strategy used to study P. vivax liver stage biology

Primary human hepatocytes were cultured and infected with P. vivax sporozoites on a 384-well plate platform (Roth et al., 2018; Maher et al., 2021a). Figure 1A outlines the workflow for the experiment. Following infection with sporozoites, cultures were assessed at two endpoints: one containing both schizonts and hypnozoites at day 5 post-infection, and the other containing solely hypnozoites at day 9 post-infection. Schizonts were eliminated from day 9 post-infection cultures by adding the phosphatidylinositol 4-kinase (PI4K) inhibitor MMV390048 (1 µM) to the culture media on days 5, 6 and 7 post-infection (Paquet et al., 2017; Gural et al., 2018). For each collection time point, we simultaneously initiated two sets of hepatocyte cultures: one infected with sporozoites for quantification of immunofluorescent labeled parasites by high-content imaging (HCI), and one infected with sporozoites for transcriptional profiling using a high-throughput, droplet-based scRNA-seq workflow (Zheng et al., 2017). As controls, a third set of naïve cultures were left uninfected and processed for scRNA-seq alongside cultures at days 5 and day 9 post-infection. At each collection point, a target of approximately 10,000 individual hepatocytes were partitioned into droplets prior to the generation of RNA libraries. Droplets contained primers with unique nucleotide sequences (cellular barcodes) that were incorporated into the transcript sequence during the reverse transcription step, allowing for the identification of each transcript’s cellular origin. After sequencing the single-cell libraries, we mapped the resulting reads onto a customized reference transcriptome containing human (Ensembl, v103) and P. vivax PVP01 (PlasmoDB, v51) information. The P. vivax PVP01 transcript information included the nucleotides encoding for the gene-flanking untranslated regions (UTRs) (Siegel et al., 2020), which is important because of the poly-dT capture strategy of the technology used in our workflow. This information partially alleviates the challenges associated with transcription quantification (Sà et al., 2020; Ruberto et al., 2022). Post mapping, we used the transcript information (i.e. human or parasite encoding) with their designated cellular barcodes to assess the number of transcripts in each cell, and to determine which hepatocytes were infected.

FIGURE 1

Figure 1 Design and validation of detection of P. vivax liver stages using high-content imaging and scRNA-seq. (A) Schematic illustrating the sample preparation and data generation pipeline used to assess P. vivax liver stage infection rates and gene expression. (B) Representative high content images of primary hepatocytes infected with P. vivax on day 5 (left) and day 9 (right) post-infection (p.i). Images were captured from one well from a 384-well plate with a 4x objective. Inset: one field of view (white dashed box) from the same well captured with a 20x objective. Cells were stained with DAPI (blue) and PvUIS4 (green). White bar represents 1mm, inset bars represent 10µm. Images are representative of replicate 1. (C) Violin plots displaying the distribution of parasite growth areas fixed at day 5 or day 9 post-infection (left) and net quantity of parasites and hepatocyte nuclei (right). Net infection rate is indicated for each endpoint. Growth and infection metrics are shown for both replicates. ****, p < 1.00E-04; **, 1.00E-03; Wilcoxon test.

As expected, we detected hepatocytes infected with either form of the parasite at 5 days post-infection and enriched for hypnozoites at 9 days post-infection following exposure to MMV390048 (Figure 1B; Figure S1). Visual assessment of infected cultures fixed and stained with an antibody against Upregulated in Infectious Sporozoites 4 (PvUIS4), which localizes to the parasite’s parasitophorous vacuole membrane, confirmed that the parasites we subjected to scRNA-seq had properly invaded hepatocytes and initiated development into hypnozoites and schizonts. Furthermore, hypnozoites displayed typical morphology relative to those found in humanized mice (Mikolajczak et al., 2015) and in vitro (Roth et al., 2018). Figure 1C displays the distribution of individual parasite growth areas, number of parasites detected, and parasitemia levels at each endpoint (Table S1A). Across replicates, the parasitemia ranged from 5.50-6.93% at 5 days post-infection and 1.76-2.77% at 9 days post-infection. Together, the high infection rates coupled with the high-throughput capacity of 10x Genomics’ single-cell isolation and sequencing workflow allowed for sufficient material to be generated for successful profiling of both parasite and host transcriptomes.

P. vivax schizonts and hypnozoites are distinguishable at single-cell resolution

We first analyzed the reads mapping to P. vivax transcripts to identify schizonts and hypnozoites, as well as to assess the extent to which these forms differed at single-cell resolution. In total, we characterized 1,438 parasite transcriptomes across two independent experimental replicates (Figure S2A). The number of parasite transcriptomes profiled corroborated the parasitemia levels quantified using high-content imaging. This correlative evidence suggests that our cell and gene filtering cutoffs effectively eliminated problematic cells (i.e. poorly captured, dying, or dead) (Figure S2B; Table S1B). Low dimensional representation of the data revealed two populations of parasites in samples profiled at 5 days post-infection (Figure 2A, left). In contrast, and as expected, samples profiled at 9 days post-infection contained one population (Figure 2A, right). Given that the population of cells profiled at day 9 post-infection, from which schizonts were chemically removed, overlapped with one of the two cell populations profiled at 5 days post-infection, we deduced that the non-overlapping population of cells encoded for schizonts. Consistent with our reasoning, an unsupervised clustering algorithm (Traag et al., 2019) unbiasedly assigned the cells to one of the two clusters (Figure 2B; Figure S2C). As expected, the number of transcripts detected, the unique molecular identifier (UMI) counts (a readout for absolute transcript abundances (Kivioja et al., 2011)), and the transcript levels of the early liver stage development marker liver stage-specific protein, LISP2 (PVP01_0304700) (Gupta et al., 2019), all displayed greater quantities in the schizont cluster (C2) relative to the hypnozoite cluster (C1) (Figure 2C).

FIGURE 2

Figure 2 Low dimensional representation of P. vivax schizont and hypnozoite transcriptomes. (A) t-SNE plots of integrated P. vivax liver stage scRNA-seq data faceted by parasites isolated on day 5 (left) and day 9 (right) post-infection. (B) t-SNE plot of P. vivax liver stage data colored by parasites encoding for putative hypnozoites (C1) and schizonts (C2). (C) t-SNE plots of P. vivax liver stage scRNA-seq data colored by number of transcripts detected (left), UMI counts (middle), and the expression of the schizont marker, LISP2 (PVP01_0304700) (right). Scales: absolute count (left, middle), normalized expression (right).

Previous scRNA-seq studies in Plasmodium spp. has shown gene expression varies with life-cycle stage, with non-replicating forms displaying less expression per cell compared to replicating forms (Howick et al., 2019). The median number of genes detected in schizonts (2323) was greater than the values observed in other single-cell assessments of replicating stages of the parasite’s life cycle (Figure S3) (Sà et al., 2020; Mancio-Silva et al., 2022). Furthermore, the median number of UMIs per schizont (6905) was also greater in the current dataset versus other reports. The median number of genes detected per hypnozoite (519) and median number of UMIs per hypnozoite (672) were consistent with other single-cell reports assessing non-replicating, uninucleate stages (Bogale et al., 2021; Real et al., 2021; Ruberto et al., 2021), including recent single-cell assessments of P. vivax sporozoites (Ruberto et al., 2022) and hypnozoites (Mancio-Silva et al., 2022) (Figure S3). The intact parasitophorous vacuole membrane of hypnozoites in cultures assessed on days 5 and 9 (Figures 1B, S1), the association between infection rates observed by IFA and the number of transcriptomes assessed, and the per cell metrics typical of Plasmodium spp. life stages with low RNA content provided us with strong assurance that the transcriptomes obtained represented viable hypnozoites. These findings show that the transcriptomes of P. vivax schizonts and hypnozoites can be distinguished using our experimental workflow.

Hypnozoites exhibit distinct transcriptional signatures compared to schizonts

We next performed differential gene expression analysis to identify transcriptome-wide differences distinguishing schizonts and hypnozoites. Using stringent inclusion cutoffs (Bonferroni adjusted p value < 0.01 & absolute (average log₂ fold change [FC]) > 0.5), we found significant differences in the expression of ~18% (836/4722) of genes detected (Figure S4A; Table S2A). Relative to hypnozoites, we found an increase in the expression levels of 368 genes in schizonts. As expected, we observed a significant increase of LISP2 gene expression in schizonts (average log₂ FC: 2.32) compared to hypnozoites. Shortly after invading the hepatocyte, non-relapsing Plasmodium spp. liver forms undergo immediate replication. Studies in rodent malaria models have shown that these forms are metabolically active and that various genes involved in energy metabolism are upregulated (Tarun et al., 2008; Caldelari et al., 2019; Toro-Moreno et al., 2020). Consistent with these prior observations, we found increased transcription associated with glycolysis (PVP01_0816000, PVP01_1229700, PVP01_1244000) and the TCA cycle (PVP01_0106100, PVP01_1332400) in schizonts (Figure S4B; Table S2A). More broadly, we found increased transcription of various metabolic processes, as well as motility, DNA-, RNA-, and protein-processing (Figure S4C; Table S2A). Genes encoding for various proteasome subunits displayed greater transcription in schizonts (Figures S5A, S5B; Table S2A). To further show the importance of this protein complex, we tested a known proteosome inhibitor, carfilzomib (Kortuem and Stewart, 2013), in our standard 12-day P. vivax small molecule screening assay (Maher et al., 2021a), and found it effectively killed schizonts (average growth area IC₅₀ 288nM; quantity IC₅₀ 373 nM) and hypnozoites (IC₅₀ 511nM) with some selectivity over host hepatocytes (IC₅₀ 3.32 µM) (Figures S5C, S5D).

Pseudobulk assessment of the cluster representing hypnozoites revealed robust expression of genes encoding for histones, translation-, and transcription- related proteins (Table S2D). The detection of genes associated with energy metabolism in hypnozoites, albeit at lower levels relative to schizonts, suggested that these forms are also metabolically active (Figure S4B). Relative to schizonts, we found an increase in expression of 468 genes in hypnozoites (Table S2A). A large proportion of these differences (410/468) represented genes with low average log2 fold-changes relative to schizonts (Figures S4A, S4D). As represented in Figure 3A, genes with greater expression in hypnozoites encoded for a diverse set of ontologies—including proteases, membrane proteins, transcription factors, and DNA-/RNA-regulating proteins. Genes encoding for proteases vivapain-1, vivapain-2 (PVP01_1248900, PVP01_0916200), and plasmepsin IV (PVP01_1340900) displayed the greatest differential expression in these forms compared to schizonts (Figures 3A, B; Table S2A). While limited research has reported the functionality of these gene products in P. vivax (Byoung-Kuk et al., 2004), work in other Plasmodium spp. has revealed that they localize to the food vacuole and play an important role in hemoglobin catabolism during the blood stages of the parasite’s life cycle (Dame et al., 2003; Drew et al., 2008).

FIGURE 3

Figure 3 P. vivax hypnozoites have distinct transcriptomic signatures. (A) Dot plots showing transcripts that distinguish all hypnozoites from schizonts. The size of the dot corresponds to the percentage of cells that transcript was detected per condition, colored by average expression. Differentially expressed transcripts were identified using Seurat’s FindMarker function (Wilcoxon rank-sum test, Bonferroni adjusted p values < 0.01). (B) tSNE plots of liver stages colored by transcripts displaying greater levels in hypnozoites. Colored boxes in A and B indicate the family or canonical biological process that the transcript is associated with. Scale: normalized expression; pct. exp: percent of cells expressing transcript. (C) Number of transcripts in schizonts and hypnozoites identified as cellular fating RBPs (left) and GO term enrichment of this subset (right). Adj P: Bonferroni adjusted p value.

We also identified genes with greater expression associated with RNA regulatory mechanisms in hypnozoites. Notably, we found greater expression of CCR4 (PVP01_1429000), a gene encoding for a protein with mRNA deadenylase activity and playing a role in mRNA turnover (Tucker et al., 2001) (Table S2A). In Plasmodium spp., it has been described as important in regulating the expression of genes associated with invasion and egress in blood stages (Balu et al., 2011). RNA-binding proteins (RBPs) also govern gene expression at the post-transcriptional level by transiently storing transcripts until later processing (Hentze et al., 2018). Pumilio domain (PUF) proteins are a major family of RBPs that regulate cellular fating through translational repression in eukaryotes. We found greater expression of the gene encoding for PUF1 (PVP01_1015200) in hypnozoites relative to schizonts (Figures 3A, B; Table S2A). Furthermore, we found greater expression of 7-Helix-1 protein (PVP01_1008400) (Figures 3A, B; Table S2A)—a stress granule component that interacts with the RNA-binding protein PUF2 and serves to be crucial for protein synthesis in P. falciparum (Miao et al., 2010). Numerous studies have shown that Plasmodium spp. use various post-transcriptional regulatory mechanisms to coordinate stage transitions and development (le Roch et al., 2004; Tarun et al., 2008; Foth et al., 2011; Bunnik et al., 2016; Vembar et al., 2016; Baumgarten et al., 2019). Our data corroborate these findings and suggest that similar strategies may be used in P. vivax to regulate liver stage development. Moreover, these data suggest that RBP-mediated cellular fating and translational repression may be involved in regulating hypnozoite quiescence.

Given the potential role in RBP-mediated cellular fating and translational repression governing hypnozoite quiescence, we performed further analyses focused on RBP-mediated regulation. To this end, we scanned the P. vivax genome searching for genes encoding for cellular fating RBPs based on ortholog groups (OGs) of known human cellular fating RBPs (348 RBPs; 227 OGs), identifying 89 OGs (99 RBPs) in P. vivax. We detected 48 genes encoding for cellular fating RBPs with greater expression in schizonts (Figure 3C, left; Table S2B). Notably, we identified eukaryotic translational initiation factors (PVP01_1238100, PVP01_0531200, PVP01_1431600, PVP01_0807600, PVP01_1005600, PVP01_0519700, PVP01_1329200, PVP01_0825000, PVP01_1303500), and splicing factors (PVP01_0716000, PVP01_0521000, PVP01_0804100, PVP01_1111800, PVP01_0607600, PVP01_1012200) suggesting active translation in schizonts. Gene Ontology (GO) term analysis revealed that these transcripts were associated with translation initiation (GO:0005852; GO:0033290; GO:0016282) (adjusted p value < 0.05, Figure 3C, right; Table S2C).

In hypnozoites, we detected 51 genes encoding for cellular fating RBPs with increased expression relative to schizonts (Figure 3C, left; Table S2B). These genes encode for prominent P body markers, namely DOZI/DDX6 (PVP01_0819400), CCR4-NOT complex proteins (PVP01_1429000, PVP01_1014400, PVP01_1331700, PVP01_1453400, PVP01_0929400), mRNA-decapping enzyme subunit 1 (PVP01_0617200), mRNA-decapping enzyme subunit 2 (PVP01_1409900), UPF1 (PVP01_0805200), UPF2 (PVP01_0724300), UPF3 (PVP01_1226700), PUF2 (PVP01_0526500), PABP (PVP01_1442500), Fibrillarin (PVP01_1341600) and trailer hitch homolog (PVP01_1269800) (Table S2B). Genes encoding for cellular fating RBPs in hypnozoite also represented GO terms corresponding to active biological condensate formations, such as the CCR4-NOT complex (GO:0030014; GO:0030015) and the P body (GO:0000932) (adjusted p value < 0.05, Figure 3C, right; Table S2C). Together, these findings suggest that P body-mediated cell fating could be important in regulating hypnozoite persistence.

Transcriptomic heterogeneity exists amongst individual P. vivax hypnozoites

Having described the transcriptomic signatures distinguishing schizonts and hypnozoites, we next sought to assess the extent to which gene expression varied amongst hypnozoites. We thus subsetted the data to include only hypnozoite transcriptomes. In total, 1,147 transcriptomes were assessed (Figure 4A, left). Low dimensional representation of the data revealed slight variation between samples collected at days 5 and 9 post-infection (Figure 4A, right). After regrouping the data using the same unsupervised clustering algorithm described previously, we identified three clusters (Figure 4B, right; Figure S6C). The proportion of hypnozoites in each cluster varied across the day of collection: specifically, the proportion of hypnozoites in cluster 3 (C3) was greater in hypnozoites collected at day 9 compared to day 5 post-infection (Figure 4B, left; Figure S6A). Albeit lower compared to schizonts, we also observed a trend for slightly greater gene expression of LISP2 in C3 relative to the other two clusters (Figure 4C).

FIGURE 4

Figure 4 Transcriptomic signatures of activating and persister hypnozoites. (A) Number of hypnozoites assessed (left) and UMAP of cells identified as hypnozoites colored by day of collection. (B) Proportion of cells in clusters on days 5 and 9 post-infection (left) and UMAP of hypnozoites colored by cluster. (C) Violin plot showing the expression of LISP2 in each cluster. (D) Dot plot showing markers that distinguish hypnozoites in each cluster. The size of the dot corresponds the percentage of cells in the expressing the gene of interest, colored by scaled normalized expression of the transcript in each cluster. Marker transcripts were identified using Seurat’s AllFindMarkers function. Wilcoxon rank-sum test, average log₂FC > 1, Bonferroni adjusted p value < 0.01. (E) Number transcripts in each cluster encoding for cellular fating RBPs sorted by their average log2FC relative to the other clusters (top) and GO Term enrichment of genes in each cluster (bottom). Adj P: Bonferroni adjusted p value.

We next identified markers that define the hypnozoites in each cluster using the FindAllMarkers function in Seurat (Stuart et al., 2019). We defined a marker as a gene displaying greater than 0.5 FC (average log₂) compared to the other clusters. In total, we detected 477 markers (Bonferroni adjusted p value < 0.05) (Figure S6B; Table S3A). Notable markers in C1 included genes associated with drug resistance (PVP01_1447300, PVP01_1259100, and PVP01_1010900), energy metabolism (PVP01_1435000 and PVP01_1122800), protein modification (PVP01_1248900 and PVP01_0916200), and RNA-binding proteins (PVP01_0939900 and PVP01_0819400); in C2, genes associated with translation (PVP01_1255500 and PVP01_0822500); and in C3, genes associated with translation (PVP01_0213700 and PVP01_0203400), fatty acid metabolism (PVP01_1143000 and PVP01_1022800), protein processing (PVP01_1308000 and PVP01_0927700) and ion transport (PVP01_1407500) (Figure 4D; Table S3A).

We next explored the data in the context of cellular fating to elucidate potential patterns in RBP signatures. Among the hypnozoite markers found within our sub-clustering analysis, we identified 80 cellular fating RBPs (Figure 4E, top; Table S3B). Approximately 4% (3/80) of the differentially expressed genes in C3 encoded for genes associated with cellular fating. In this cluster, the increased gene expression of eukaryotic translation initiation factor 4E (PVP01_0203400) further suggests active translation events may be occurring. C2 contained 20% (16/80) of the cellular fating RBPs differentially expressed. In this cluster, we found RBPs which are important determinants of cell proliferation from a state of quiescence such as eukaryotic translation initiation factor 5A (PVP01_1303500), polypyrimidine tract binding protein (PVP01_1142400), GTP-binding nuclear protein RAN/TC4 (PVP01_0918300), and transformer-2 protein homolog beta (PVP01_0802200). C1 contained the majority of the differentially expressed genes encoding for RBPs (~76%, 61/80). Within this cluster, we identified prominent P body markers with cellular fating ability, such as DOZI/DDX6 (PVP01_0819400), CCR4-NOT complex proteins (PVP01_1429000, PVP01_1014400, PVP01_1331700, PVP01_1453400), UPF1 (PVP01_0805200), UPF2 (PVP01_0724300), PABP (PVP01_1442500), Fibrillarin (PVP01_1341600) and trailer hitch homolog (PVP01_1269800). GO analysis revealed enrichment of cellular compartments associated to biological condensate formations, namely CCR4-NOT complex (GO:0030014), P granule (GO:0043186), and germplasm (GO:0060293) (Figure 4E, bottom; Table S3C).

These findings, combined with the shift of hypnozoite proportions and increasing gene expression of LISP2 from C1 to C3, allowed us to infer the phenotype of hypnozoites in these clusters. As these observations likely represent populations of hypnozoites in persisting and activating states, we defined clusters C1, C2, and C3 to represent persisting, early activating, and late activating hypnozoites, respectively. Together, these findings further our understanding of P. vivax hypnozoites through characterization of gene usage and putative regulatory mechanisms in activating and persister forms. Furthermore, they highlight a potential role of P body-mediated cellular fating governing persister hypnozoites.

Strategy used to assess hepatocytes infected with P. vivax parasites

After deciphering the identities of P. vivax transcriptomes representing schizonts and hypnozoites, we next sought to characterize the transcriptional signatures of hepatocytes infected with either of these forms. To this end, we first removed all reads derived from the cell-gene matrices derived from P. vivax to remove their influence on downstream analyses. Next, we integrated the data across all the conditions (Figures S7A, S7B) and examined whether distinct sub-populations of infected hepatocytes existed. As highlighted by the low dimensional visualization of the data, we did not find a prominent pattern of infection (Figure 5). Moreover, we detected parasites in hepatocytes with different metabolic states (MacParland et al., 2018), as represented by their transcriptional signatures associated with different hepatocyte zonation programs (Figures S7C, S7D).

FIGURE 5

Figure 5 Analysis of host-pathogen transcriptional signatures in P. vivax- infected and non-infected hepatocytes. t-SNE plots of hepatocytes are faceted by infection status. Grey, noninfected (naive and exposed) hepatocytes; light blue, hepatocytes infected with schizonts; teal, hepatocytes infected with hypnozoites.

Identification of transcriptional signatures in P. vivax-infected hepatocytes

To identify host transcriptional signatures specific to infection status, we performed differential gene expression analysis, comparing infected to non-infected hepatocytes. Among the hepatocytes infected with schizonts, 106 genes displayed greater expression and 16 displayed decreased expression (Table S4A); while in hepatocytes infected with hypnozoites, 365 genes displayed greater expression and 41 genes displayed decreased expression (Table S4B) (Bonferroni adjusted p value < 0.05). As depicted in the volcano plots in Figure 6A, the magnitude of change in expression in the infected versus non-infected cells (average log₂ fold change) was greater in hepatocytes containing schizonts compared to those containing hypnozoites.

FIGURE 6

Figure 6 Analysis of host transcriptional signatures in P. vivax infected hepatocytes. (A) Volcano plots showing changes in gene expression in infected hepatocytes versus non infected (Wilcoxon rank sum, Bonferroni adjusted p value < 0.05). Positive fold change (FC) values represent genes with greater expression, and negative FC values represent genes with decreased expression. Dashed horizontal red lines: Bonferroni adjusted p value = 0.05; dashed vertical red lines: log₂FC = 0.5. (B) Summary of differential gene expression data shown in the volcano plots and of the enrichment analyses. Hep: hepatocyte. (C) Enrichment map of cellular processes and pathways associated with hepatocyte infection with P. vivax schizonts. (D) Enrichment map of cellular processes and pathways associated with hepatocyte infection with P. vivax hypnozoites. For the maps depicted in panels C and D, node size is proportional to the number of genes identified in each gene set (minimum 3 genes/gene set); each node represents a distinct biological process or pathway derived from gene with decreased expression (red) or increased expression (blue) versus non-infected cells; and edges (blue lines) represent the number of genes overlapping between two processes or pathways. (E) αRepresentative confocal image of P. vivax parasites on day 12 post-infection of hepatocytes (left). Cells were stained with DAPI (blue), PvUIS4 (green), and HsAKR1B10 (red). White bar represents 100µm. Violin plot displaying the distribution of AKR1B10 signal in hepatocyte cultures infected with P. vivax (right). SZ: hepatocyte containing schizont; HZ: hepatocyte containing hypnozoite; EX: hepatocyte not infected but exposed to P. vivax. ****: p <= 0.0001, Wilcoxon test.

To gain greater insights into the biological processes and pathways that characterize a schizont- or hypnozoite-infected hepatocyte, we performed enrichment analysis (Raudvere et al., 2019)—scoring 18,594 gene sets against the genes differentially expressed in each of our transcriptional comparisons. When considering both the schizont- and hypnozoite-infected hepatocyte populations, we detected a total of 254 enriched gene sets (Bonferroni adjusted p value < 0.05) (Figure 6B; Table S4C). Approximately 18% (45/254) of these gene sets overlapped in hepatocytes infected with either schizonts or hypnozoites (Table S4D). These data indicate that while some overlap exists in the gene networks the parasite alters during infection, the majority of processes and pathways altered are specific to schizont versus hypnozoite infected hepatocytes.

Using these gene sets, we then constructed network enrichment maps (Merico et al., 2010) to provide a global overview of the prominent biological themes associated with each infection condition (Figures 6C, D). We detected upregulation of genes associated with pathways related to energy metabolism (oxidative phosphorylation, mitochondrial function, and lipid catabolism) in schizont- and hypnozoite-infected hepatocytes. Analysis of the most differentially expressed genes associated with these pathways included MT-ND2, CHCHD10, and AKR1B10 (Tables S4A, S4B). Corroborating our transcriptomic data, we found an increase in protein levels of AKR1B10 in infected versus non-infected hepatocytes by IFA (Figures 6E and S8). AKR1B10 has been implicated in cell survival through its regulation of lipid metabolism and elimination of carbonyls (Martin and Maser, 2009; Wang et al., 2009; Qu et al., 2021) and has primarily been described as a biomarker of cancer (Fukumoto et al., 2005; Jin et al., 2006; Díez-Dacal et al., 2011). More recently, it has been shown that its expression is positively correlated with non-alcoholic fatty liver disease (Govaere et al., 2020). In schizonts, the upregulation of genes linked to fatty acid metabolism is consistent with observations in hepatocytes harboring developing forms from non-relapsing Plasmodium spp. (Albuquerque et al., 2009; Itoe et al., 2014). We now reveal that infection with hypnozoites also elicits enrichment of genes associated with fatty acid metabolism (Figure 6D), suggesting that these forms may be scavenging host resources to fuel their persistence.

Previous work assessing Plasmodium liver stage biology has revealed that the host is able to detect the parasite (Liehl and Mota, 2012) and that type I interferon plays an integral role in this response (Liehl et al., 2013; Miller et al., 2014). The decrease in expression of genes associated with interferon signaling (REAC:R-HSA-913531, GO:0034340, GO:0034341, GO:0071357, GO:0060337) in hepatocytes infected with either schizonts or hypnozoites suggest that P. vivax manipulates the hepatocytes to evade detection (Figures 6C, D; Tables S4C, S4D). Furthermore, we found that relative to naïve hepatocytes, uninfected hepatocytes in cultures exposed to P. vivax displayed significant downregulation of pathways associated with immune effectors (GO: 0030449, GO:0002697, REAC: R-HSA-166658, R-HSA-168249) (Tables S4E, S4F).

In hepatocytes infected with schizonts, we found a decrease in expression of NFKB-regulated genes (SOD2, TMEM45A, CCL20, and TPM1) (Table S4A), consistent with findings that Plasmodium spp. release signals to block proinflammatory responses (Singh et al., 2007). These findings, in addition to the decrease in the expression of HLA-B—a molecule involved in antigen presentation to T cells that can trigger their activation when ‘self’ peptides are not presented—and the neutrophil chemoattractant CXCL8 suggest that replicating parasites alter hepatocyte signaling pathways to minimize recruitment of innate immune cells and avoid detection by T cells. Similar to hepatocytes infected with schizonts, we found genes associated with interferon signaling downregulated in hepatocytes infected with hypnozoites, namely: UBB, UBC, HSP90AB1, IFITM3, IFITM2, and HLA-A (Table S4B). Furthermore, transcripts associated with the downregulation of class I major histocompatibility complex (CD81 and B2M) were also decreased in these hepatocytes (Table S4B). Together, these findings suggest that P. vivax infection of hepatocytes may suppress the secretion of various chemokines and proinflammatory cytokines to thwart the recruitment of inflammatory cells to the site of infection while simultaneously downregulating MHC- I to reduce the chances they will be detected if these cells are successfully recruited.

Enrichment analysis linking the changes in gene expression in infected hepatocytes with transcriptional factors revealed distinct sets of putative transcriptional regulators. In total, we detected 35 transcription factors with potential roles in governing the gene expression changes in infected hepatocytes (adjusted p value < 0.05, Figure 7; Table S5). In infected hepatocytes, we found enrichment of transcripts governed by NRF2 (Figure 7; Table S5). The regulation of gene expression through NRF2 in infected hepatocytes may serve to minimize reactive oxygen species in the cell (Hayes and Dinkova-Kostova, 2014). Thus, in addition to the decreased expression of genes associated with inflammatory response (Table S4C), NRF2-mediated regulation of antioxidant responsive genes may protect hepatocytes from Plasmodium-induced stress and thus serve as a complementary means of minimizing the host inflammatory response.

FIGURE 7

Figure 7 Transcription factor enrichment in hepatocytes infected with schizonts or hypnozoites. Dot plot of transcription factors enriched in Hepatocyte^+hypnozoite (left) or Hepatocyte^+schizont sets (right). Positive (+) enrichment derived from genes displaying greater expression versus non-infected hepatocytes; negative (-) enrichment derived from genes displaying decreased expression versus noninfected hepatocytes.

Discussion

The infection of the liver by P. vivax constitutes an obligatory step in its life-cycle. In vivo studies have shown that within hepatocytes, P. vivax can exist in one of two phenotypically-distinct forms: schizonts or hypnozoites (Cogswell, 1992; Mikolajczak et al., 2015). They have also been observed ex vivo by infection of hepatocyte monolayers (Hollingdale et al., 1985; Sattabongkot et al., 2006; Gural et al., 2018; Roth et al., 2018). A key problem hindering the study of P. vivax liver stage biology is the difficulty involved in obtaining a sufficient number of infected hepatocytes to make meaningful comparisons between them. We show that by coupling our in vitro liver stage model (Roth et al., 2018; Maher et al., 2021a) with a high-throughput single-cell capture technology (Klein et al., 2015; Macosko et al., 2015; Zheng et al., 2017) that transcriptome-wide characterization of P. vivax-infected hepatocytes is feasible. The findings provide much needed insight into the molecular signatures of P. vivax hypnozoites and the host cells they infect.

Our work corroborates the growing number of Plasmodium spp. scRNA-seq data sets, revealing highly variable gene expression patterns across, and within, the various stages of the parasite’s life cycle. Here, the number of genes detected, as well as the UMIs, per hypnozoite is in line with other single-cell reports assessing the transcriptomic signatures of non-replicating stages of the parasite’s life cycle (Howick et al., 2019; Bogale et al., 2021; Real et al., 2021; Ruberto et al., 2022), including hypnozoites (Mancio-Silva et al., 2022). Our transcriptome-wide assessment supports the notion that hypnozoites are metabolically active (Adams and Mueller, 2017), demonstrated by the transcription of genes associated with various biological functions including glycolysis and fatty acid metabolism. Furthermore, the expression of genes associated with oxidative stress protection, protein export, mitochondrial respiration, and epigenetic mechanisms such as histone methylation and acetylation, are detectable in hypnozoites, suggesting these parasites are both viable and actively transcribing genes to survive.

We also found that hypnozoites exist in various transcriptomic states. We speculate that these differences might represent a spectrum of phenotypes, from persisting to activating. Interestingly, no loss in gene diversity or expression was found in hypnozoites associated with cluster C1 (persisters) compared to C2 and C3 in our hypnozoite-directed analysis, suggesting this subpopulation encodes for viable, distinct forms. In the persister hypnozoite cluster, the differential expression of genes encoding for cellular fating RBPs supports the hypothesis that these forms may use post-transcriptional mechanisms to control gene usage (Mueller et al., 2019). The biological significance of the enrichment of genes associated with translational repression, found exclusively in a subpopulation of hypnozoites (C1), suggests this mechanism may govern quiescence, and warrants further investigation.

Discerning between non-replicating, uninucleate forms and dead parasites poses a challenge. The absence of a validated hypnozoite viability marker further complicates assessments. These limitations notwithstanding, the data reported here seem to suggest that the transcriptomes assessed represent viable parasites. Similar to a recent study using a non-relapsing rodent malaria model (Afriat et al., 2021), we inferred parasite viability based on visual assessment of fixed cell cultures using an immunofluorescence assay (IFA). For each scRNA-seq library processed, we dedicated cultures from the same infection to provide a visual snapshot of the parasite viability at each end point. While Afriat and colleagues (2021) define abortive forms as displaying “vacuole breakdown,” upon inspection of our cultures, we did not find indications of abortive or dead parasites. Moreover, the parasites morphologically identified as hypnozoites were identical in morphology to those described in humanized mice (Mikolajczak et al., 2015) and in vitro (Roth et al., 2018). More specifically, they displayed positive PvUIS4 staining for an intact parasitophorous vacuole membrane and a distinctive prominence. These forms have been shown to activate in vivo (Mikolajczak et al., 2015) as well as in vitro (Maher et al., 2020). Additionally, they have been shown to be susceptible to treatment with primaquine in combination with chloroquine and ionophores, implying they were alive to begin with (Maher et al., 2021b). At the molecular level, we did not find consistencies between our data and the transcriptomic signatures of “abortive” P. berghei liver forms reported by Afriat et al. (2021). For example, orthologous expression of genes encoding for HSP90, HSP70, and HOP—three markers defining abortive parasites in their dataset—displayed greater expression in the schizont cluster relative to hypnozoites. These observations further support the hypothesis that hypnozoite transcriptomes assessed in this study were derived from viable parasites.

Several observations suggest protein degradation may be important for the parasite’s liver stages. Targeting mechanisms related to this process may therefore be an effective means of killing these forms. For instance, the ability of hepatocytes to uptake and break down hemoglobin (Bissell et al., 1972) supports that this protein substrate could be an essential amino acid source for the parasite. Our data suggest that several mechanisms may allow P. vivax to utilize the amino acids derived from this protein: or more broadly, any protein substrate. In hypnozoites, the increased expression of genes encoding for proteases, such as vivapains, supports the idea that these parasites digest host cytosolic proteins. The role of these proteins as hemoglobinases have been defined extensively in P. falciparum blood stages (Salas et al., 1995; Drew et al., 2008; Rosenthal, 2011). A similar role has been described for P. vivax (Byoung-Kuk et al., 2004), but there is limited information about their significance in the liver stages. Interestingly, a recent report using a complementary single-cell capture technology also found that these genes display greater expression in hypnozoites relative to schizonts (Mancio-Silva et al., 2022). Vivapains are highly homologous to cathepsin protease L (TgCpl) in the related apicomplexan Toxoplasma gondii. In these parasites, TgCpl is localized to its plant-like vacuolar compartment—a lysosome-like organelle containing a variety of proteases (Parussini et al., 2010)—and one of its functions is to aid in the digestion of host-derived proteins (Dou et al., 2014). Genetic or chemical disruption of TgCpl negatively impacts the viability of bradyzoites, a persisting form of this parasite (di Cristina et al., 2017). Given these orthologous proteins’ roles in P. falciparum and T. gondii, we speculate that the expression of vivapains may support the persistence of P. vivax hypnozoites.

We highlight distinct transcriptomic signatures of infected versus non-infected hepatocytes, and show that signaling pathways associated with energy metabolism, antioxidant stress, and immune response are altered in infected hepatocytes. The extent to which individual hepatocytes respond to P. vivax schizonts and hypnozoites has previously been unknown. A bulk RNA sequencing approach has been used to identify pathways altered in hepatocytes in response to Plasmodium infection (LaMonte et al., 2019); however, it is limited insofar as it averages gene expression from millions of cells, obscuring any parasite-specific effects on an infected cell. Targeted approaches have revealed P. vivax liver forms recruit host factors to facilitate growth (LaMonte et al., 2019; Posfai et al., 2020), but how certain infected hepatocyte populations contribute to disease transmission, and how hepatocyte subpopulations support development, are intriguing topics for future studies (Vijayan et al., 2021). Our data identify potential host factors essential for P. vivax’s successful liver stage development. We found that some genes of the hepatocyte transcriptome co-vary with the presence of P. vivax mRNA. In particular, the presence of P. vivax in individual hepatocytes is associated with the expression of genes involved in energy metabolism and response to oxidative stress. It will be interesting to determine if these transcriptional signatures are simply consequences of P. vivax infection, or if their higher expression in some hepatocytes is a causative factor that promotes P. vivax infection.

Persistence of hypnozoites is contingent on the longevity of the hepatocyte in which it resides. The increased expression of genes associated with antioxidant stress (Scarpulla, 2002) is consistent with the hypothesis that hypnozoite infection induces changes in the hepatocyte to promote cell survival, in which NRF2 could be playing a regulatory role (Bindschedler et al., 2022). One host gene upregulated in infected hepatocytes validated by IFA, AKR1B10, is regulated by NRF2, which activates protective pathways in response to oxidative stress (Tebay et al., 2015). We also detected genes displaying greater expression in infected hepatocytes regulated by NFYB, a transcription factor important for cell proliferation, mitochondrial integrity, and cellular longevity (Tharyan et al., 2020). At day 5 post-infection, the increased expression of mitochondria-related genes in hepatocytes infected with schizonts may be required to meet the parasite’s energy demands during replication, but the enrichment of mitochondria-related genes in cells infected with hypnozoites could also aid in ensuring longevity of the hepatocyte. In addition to these positive regulators of cell survival, we found a decreased expression of genes associated with cell death such as UBB, UBC, and HSP90AB1, potentially regulated by CEBPD, a transcription factor with a role in cell death (Tsai et al., 2017; Wang et al., 2019).

Persistence of hypnozoites also requires evading host detection. Our analyses indicate the parasite alters the host hepatocyte to evade detection by immune cells. First, there is a decrease in expression of genes encoding for chemoattractants (CCL20, CXCL8, IL32) in infected hepatocytes, and an increase in expression of PTGR1, a negative regulator of the chemotactic factor leukotriene B4 (Yokomizo et al., 1993). These changes would result in decreased secretion of chemokines that recruit innate immune cells to the liver, thereby reducing recruitment of adaptive immune cells like CD8 T-cells that kill infected hepatocytes (Chakravarty et al., 2007). Interestingly, there was also a decrease in expression of MHC-1 molecules (HLA-A and -B). The presentation of parasite peptides to be detected by surveying CD8 T-cells would therefore be reduced (Neefjes et al., 2011), creating an additional layer of protection should innate immune cells be recruited, for example upon lysis of the hepatocyte by schizonts. We propose that P. vivax infected hepatocytes are re-programmed to prevent recruitment of innate immune cells on the one hand, and to reduce the likelihood of their detection by the adaptive immune system on the other, all to enable progress through schizogony and for hypnozoites to persist undetected. Understanding the mechanisms leading up to this warrant further investigation, as their reversal may lead to improved clearance of hypnozoites, thus reducing relapses. Indeed, a vaccine would need to decrease relapses to have an impact on P. vivax transmission and disease (White et al., 2017), but the decrease of MHC-I expression in infected hepatocytes may mean that standard focuses to generate CD8 and CD4 polyfunctional T-cells in the liver may not effectively eliminate P. vivax. Recent attempts to test anti-relapse vaccine candidates with P. cynomolgi did not prevent relapses despite higher inducing IFNγ producing cells, a common metric used as a correlate of protection (Kim et al., 2020), indicating that alternative strategies may be needed.

Recently, another study assessing P. vivax liver stage biology at single-cell resolution was published (Mancio-Silva et al., 2022). On the parasite side, despite the greater expression of nearly all of the genes detected in the current versus the recently published study, our studies are largely complementary at the single common timepoint assessed (day 5 post-infection). The expression of genes encoding for vivapains, PUFs, LISP2, and glycolysis enzymes such as GAPDH and L-lactate dehydrogenase were conserved, adding strong evidence of these genes’ roles in governing liver stage biology. Furthermore, we found differential expression of genes encoding for canonical sexual markers, with greater expression of PVS16 (PVP01_0305600) observed in schizonts relative to hypnozoites, and greater expression of P25 (PVP01_0616100) and G377 (PVP01_1467200) in a subset of hypnozoites. On the host side, however, the transcriptional signatures of hepatocytes in response to P. vivax infection were more divergent. While we detected an upregulation of various antioxidant response genes in P. vivax-infected hepatocytes on day 5 post-infection (i.e. PTGR1, SOD1, GPX2, TXNRD1), and even validate one of these markers (AKR1B10) using IFA, an opposite transcriptional response is observed in the recently published study. More broadly, a comparison of the differentially expressed genes (infected versus exposed but non-infected hepatocytes (i.e. bystander hepatocytes) in each study revealed that only 17 genes overlapped (of 615 genes); and of which, only five (HSP90AB1, SQSTM1, ANXA2, UBC, and QPRT) displayed similar changes in response to P. vivax-infection: downregulation relative to exposed, non-infected hepatocytes. In our data, we also reported the upregulation of genes associated with pathways related to energy metabolism (i.e. oxidative phosphorylation, mitochondrial function). While few genes overlapped between studies, upon reprocessing of the day five post-infection data made available by Mancio-Silva et al. (2022) (see Materials and Methods), we also observed an upregulation of mitochondria-encoding genes (MT-CO2, MT-ND3, MT-RNR1, MT-RNR2, MT-ND2, MT-CYB, MT-CO3, MT-ND4, MT-ND5, MT-CO1, and MT-ATP6) in infected relative to exposed but non-infected hepatocytes, supporting the hypothesis that P. vivax liver forms alter gene networks in hepatocytes associated with energy metabolism.

We attribute the differences in the studies (both parasite- and host- signatures) to their respective culture platforms, single cell capture technologies, and bioinformatic strategies used. In the present study, we used primary human hepatocytes (donor BGW, BioIVT) infected with P. vivax from Cambodian field isolates, while Mancio-Silva et al. (2022) used primary human hepatocytes (donor UBV, BioIVT) co-cultured with murine embryonic fibroblasts and infected with P. vivax from Thailand field isolates. Recent work shows that hepatocyte donor influences P. vivax infection dynamics (Vantaux et al., 2022); therefore, differences in the genetic background of hepatocytes may influence the host’s transcriptional response to infection. Additionally, the infection inoculum has been shown to affect P. vivax infection rates (Vantaux et al., 2022). In our study, inoculums of >15,000 sporozoites per well resulted in P. vivax-infected hepatocytes ranging from 5.50-6.93% at 5 days post-infection and 1.76-2.77% at 9 days post-infection following treatment with MMV390048 to enrich for hypnozoites. The inoculum and infection rates are not reported by Mancio-Silva et al. (2022), but it would nevertheless be interesting to consider the extent to which the differences in host gene signatures could be a result of the studies’ respective strategies.

In addition to the liver stage platforms, there were differences between the respective scRNA capture methods. In the present study we used 10x Genomics’ droplet-based capture method to assess parasite and host gene expression signatures, while Mancio-Silva et al. (2022) used Seq-Well S³ technology followed by customized nucleic acid baits targeting the entire P. vivax genome. Our use of 10x Genomics’ 3’ gene expression technology allowed for a stream-lined workflow with minimal handling: after cells were detached from a 384-well microplate surface they were partitioned, barcoded, and lysed in under 10 minutes. The Seq-Well S³ technology is similarly effective for partitioning thousands of cells for sequencing, but unlike the 10x Genomics’ workflow, its protocol requires greater handling time (minimum of 40 minutes) before the cells are lysed (see Seq-Well S³ protocol version June 2019 from https://shaleklab.com/wp-content/uploads/2017/03/Hughes.et_.al_.Master.Protocol.pdf). The studies’ respective single cell capturing strategies may have further contributed to diverging transcriptional signatures between the two datasets.

Differing bioinformatic strategies could potentially have contributed to divergent gene signature results reported in the two studies. In our workflow, we used a transcriptome alignment protocol (Bray et al., 2016; Melsted et al., 2021), whereas Mancio-Silva et al. (2022) used a genome alignment protocol (Dobin et al., 2013). We tested a genome-alignment approach during the optimization process of our bioinformatic pipeline, finding negligible differences between the two protocols (i.e. the cell clustering and differential gene signatures among the various comparisons highlighted in the present study were conserved).

By examining two early liver stage time points at single-cell resolution, herein we describe a technical advance demonstrating the importance of in vitro models to shed light on the molecular basis of P. vivax liver stage development and infection. We hope this better understanding will facilitate the development of new parasite- and/or host- directed therapies for improved malaria control and elimination. The comprehensive mapping of the transcriptional landscape of both host and P. vivax described here provides an important framework for further investigation.

Data availability statement

All raw sequencing data generated in this study can be found in the Sequence Read Archive (SRA) at the NCBI National Library of Medicine (https://www.ncbi.nlm.nih.gov/sra) under the BioProject code: PRJNA843856. Scripts containing the code used to process the single-cell RNA seq data are available on GitHub at: https://github.com/AnthonyRuberto/Pv_LS_singleCell. Archived scripts (Shell and R) and output files as at time of publication are available at Zenodo (10.5281/zenodo.6463338).

Ethics statement

All research procedures were reviewed and approved by the Cambodian National Ethics Committee for Health Research (approval numbers: #113NECHR, #104NECHR and #089NECHR). The patients/participants provided their written informed consent to participate in this study.

Author contributions

AR, SM, AV, IM, BW, and DK conceived and designed the study; AR, SM, and AV performed the experiments; AR, CB, BB, and CJ curated the data; AR, SM, CB, BB, AJ, and CJ performed the formal analyses; SM, AV, IM, BW, and DK coordinated the research activities; SM, AV, IM, BW, and DK provided resources; SM and DK acquired funding for the work; AR wrote the initial draft of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

High-content imaging data was produced in part by the Biomedical Microscopy Core at the University of Georgia, supported by the Georgia Research Alliance (DK). AJ, BB, and IM acknowledge support from the Victorian State Government Operational Infrastructure Support and Australian Government National Health and Medical Research Council Independent Research Institute Infrastructure Support Scheme. This work was also supported by the Agence Nationale de la Recherche (ANR-17-CE13-0025 to AR and IM); an Australian National Health and Medical Research Council Leadership Fellowship (APP1194330 to AJ); the Bill & Melinda Gates Foundation (OPP1023601 to DK); and Medicines for Malaria Venture (RD/2017/0042 to BW and AV, RD/16/1082 and RD/15/0022 to SM and DK).

Acknowledgments

We thank the patients of Mondulkiri Province, Cambodia, for participating in this study. We thank the Institut Pasteur du Cambodge’s field site manager (Saorin Kim) for logistical assistance, insectary staff (Makara Pring, Koeun Kaing, Nora Sambath) for the An. dirus mosquito colony maintenance, and laboratory staff (Eakpor Piv, Chansophea Chhin, Sreyvouch Phen, Chansovandan Chhun, Sivcheng Phal, Baura Tat) for assistance with the mosquito dissections and the in vitro assays.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2022.986314/full#supplementary-material

Supplementary Table 1 | Summary statistics and cell and gene filtering metrics for each single-cell RNA library.

Supplementary Table 2 | Data associated with schizonts versus hypnozoites analyses. (A) Differentially expressed genes defining schizonts and hypnozoites. Positive Average log 2 FC: Upregulated in hypnozoites relative to schizonts. Rows highlighted in yellow denote inclusion cutoffs for genes considered differentially expressed (Bonferroni adjusted p value < 0.01 & absolute (average log2 fold change [FC]) > 0.5). (B) Differentially expressed genes identified as encoding for cellular fating RNA-binding proteins. (C) GO Term enrichment of Hypnozoite or Schizont transcripts identified as cellular fating RBPs. Rows highlighted in light green are terms with Bonferroni adjusted p value < 0.05. Input parameters: P value cutoff = 0.05; Limit to GO Slim terms = No; Evidence: Computed, curated; PlasmoDB v55. (D) Pseudobulk expression of genes detected in hypnozoites (days 5 and 9 post-infection).

Supplementary Table 3 | Data associated with analyses of hypnozoite-to-hypnozoite heterogeneity. (A) Markers defining each hypnozoite cluster. Marker defined as: Average log 2 FC > 0.5, Adjusted p value < 0.05. Rows colored green are transcripts that meet marker criteria. (B) Differentially expressed genes identified as encoding for cellular fating RNA-binding proteins in hypnozoite clusters. (C) GO Term enrichment of hypnozoite transcripts identified as cellular fating RBPs. Rows highlighted in light green are terms with Bonferroni adjusted p value < 0.05. P value cutoff = 0.05; Limit to GO Slim terms = No; Evidence: Computed, curated; PlasmoDB v55.

Supplementary Table 4 | Data associated with host-specific transcriptional signatures in response to P. vivax infection. (A) Differentially expressed genes in hepatocytes containing schizonts versus non-infected cells. Positive average log2FC: greater expression in infected hepatocytes. (B) Differentially expressed genes in hepatocytes containing hypnozoites versus non-infected cells. Positive average log2FC: greater expression in infected hepatocytes. (C) Enrichment of GO biological processes and Reactome metabolic pathways. (D) Unique and overlapping enriched gene sets in hepatocytes infected with schizonts or hypnozoites. (E) Differentially expressed genes in non-infected hepatocytes exposed to P. vivax (exposed) versus non-infected hepatocytes never exposed to P. vivax (naïve). Positive average log2FC: greater expression in non-infected hepatocytes exposed to P. vivax (exposed). Seurat FindMarkers parameters used to generate gene list: logfc.threshold = 0.25, min.pct = 0.25, only.pos = F, max.cells.per.ident = 20731, assay = “RNA”. (F) Enrichment of GO biological processes and REACTOME metabolic pathways that are downregulated in non-infected hepatocytes exposed to P. vivax (exposed) versus non-infected hepatocytes never exposed to P. vivax (naïve). Genes in Table S4E highlighted in green were used as input to generate enrichment lists.

Supplementary Table 5 | Transcription factor enrichment analysis. Rows highlighted in green are sets that are significantly enriched, adjusted p value < 0.05.

Supplementary Figure 1 | Assessment of P. vivax liver stages using high-content imaging. Representative high content images of primary hepatocytes infected with P. vivax on day 5 (left) and day 9 (right) post-infection. Images were obtained from one well from a 384-well plate with a 4x objective. Inset: one field of view (orange box) from the same well captured with a 20x objective. Cells were stained with DAPI (blue) and PvUIS4 (green). White arrow: liver form assigned as a hypnozoite, yellow arrow: liver form assigned as a schizont. White bar represents 1mm, grey bar represents 200µm. Images are representative of the second biological replicate used in the study.

Supplementary Figure 2 | Metrics associated with P. vivax liver stage scRNA-seq data analyses. (A) Number of P. vivax liver stage transcriptomes assessed post-cell and gene filtering. (B) Scatter plots displaying the number of genes detected versus total number of UMIs for each parasite transcriptome assessed. Red dashed vertical line (60) represents the cut-off used to filter viable from problematic (dead/dying/poorly captured) cells. (C) Number of P. vivax liver stage parasites assigned to cluster 1 (hypnozoites) and cluster 2 (schizonts). SZ: schizont; HZ: hypnozoite; dpi: day post-infection.

Supplementary Figure 3 | Comparison of per cell metrics with other P. vivax scRNA-seq data. Violin plots showing the distribution of genes detected in the current dataset, and other single-cell gene expression assessments of P. vivax liver stages (Mancio-Silva et al., 2022), sporozoites (Ruberto et al., 2022), and blood-stage parasites (Sà et al., 2020). Replicates in sporozoite and blood-stage datasets represent unique 10x Genomics’ scRNA-seq library preparations. Blood-stage scRNA-seq metrics were obtained after realignment of sequencing data to an updated P. vivax transcriptome including UTRs as performed previously (Ruberto et al., 2022). N = number of single-cell transcriptomes assessed for each sample. Light blue dashed horizontal line: median number genes detected or UMIs in schizonts obtained in the current dataset; teal dashed horizontal line: median number genes detected or UMIs in hypnozoites obtained in the current dataset. Day denotes day post-infection.

Supplementary Figure 4 | P. vivax schizonts and hypnozoites have distinct transcriptomic signatures. (A) Scatterplot showing average log₂ fold change versus the difference in the proportion of cells the transcript is detected in clusters encoding for hypnozoites and schizonts. Negative values: greater detection in hypnozoites; positive values: greater detection in schizonts. (B) Violin plots displaying the expression of genes encoding for TCA- and glycolysis-related proteins. (C) Dot plots showing transcripts with decreased expression in hypnozoites relative to schizonts. The size of the dot corresponds to the percentage of cells expressing the gene, colored by average expression. Differentially expressed transcripts were identified using Seurat’s FindMarkers function. Wilcoxon rank-sum test, Bonferroni adjusted p values < 0.01. Scale: normalized expression; pct. exp: percent of cells expressing the gene. ldh: L-lactate dehydrogenase; gapdh: glyceraldehyde 3-phosphate dehydrogenase; atp-synthC: ATP synthase subunit C; gdh1: NADP-specific glutamate dehydrogenase 1.

Supplementary Figure 5 | Carfilzomib inhibits P. vivax schizonts and hypnozoites. (A) Violin plots displaying the distribution of expression of select genes encoding for proteasome subunits in schizonts and hypnozoites. (B) Table of proteasome-associated transcripts differentially expressed between schizonts and hypnozoites. Positive average log2FC: higher expression in schizonts. (C) Structure of carfilzomib. (D) Dose-response curve of carfilzomib-induced inhibition of P. vivax schizont quantity per well (purple, IC₅₀ 373 nM), schizont net growth area per well (yellow, IC50 288nM), hypnozoite quantity per well (blue, IC₅₀ 511nM), and hepatocyte nuclei quantity (orange, IC₅₀ 3.32 µM). Data shown are pooled from three independent experiments, each containing duplicate wells at each concentration. Bars represent ± SEM.

Supplementary Figure 6 | Assessment of hypnozoite transcriptomes. (A) UMAP of hypnozoites colored by cluster and faceted by day post-infection. (B) Violin plot showing the average log₂FC of differentially expressed genes in each cluster. Markers were identified using Seurat’s AllFindMarkers function. Genes with a Bonferroni adjusted p value < 0.05 (Wilcoxon rank-sum) and average log₂FC > 0.5 were deemed markers. Red dashed line: average log₂FC = 0.5. (C) Populations identified in the hypnozoite-directed assessment of the data overlaid on the TSNE plot generated for .

Supplementary Figure 7 | Analysis of host-pathogen transcriptional signatures in P. vivax- infected and non-infected hepatocytes. (A) Number of hepatocytes assessed. Colored by infection status and faceted by sample (left); proportion of infected, exposed, or naive hepatocytes assessed on day 5 or day 9 post-infection. (B) t-SNE plot of hepatocytes colored by day. (C) t-SNE plots of hepatocytes colored by percent expression of zonation markers in human hepatocytes. Markers obtained from Macparland et al. (2018). (D) t-SNE plots of hepatocytes colored by hypnozoite infection status. Grey, uninfected (naive and exposed) hepatocytes; red, hepatocytes infected with C1 hypnozoites; blue, hepatocytes infected with C2; teal, hepatocytes infected with C3 hypnozoites.

Supplementary Figure 8 | Human AKR1B10 is upregulated in hepatocytes infected with P. vivax. (A) Violin plots displaying the distribution of AKR1B10 transcript levels in infected (schizonts and hypnozoites) and non-infected (exposed and naive) hepatocytes. Black dot: mean expression (log₂). Adj p: Bonferroni adjusted p value versus non-infected (naïve and exposed) hepatocytes. (B) Representative confocal image from a second biological replicate of P. vivax parasites on day 12 post-infection of hepatocytes. Cells were stained with DAPI (blue), PvUIS4 (green), and HsAKR1B10 (red). White bar represents 100µm. (D) Associated with (right), violins plots displaying the distribution of AKR1B10 signal across all biological replicates. IFA: immunofluorescence assay; SZ: hepatocyte containing schizont; HZ: hepatocyte containing hypnozoite; EX: hepatocyte not infected but exposed to P. vivax. ****: p <= 0.0001; *: p <= 0.05; n.s: not significant, p > 0.05, Wilcoxon test.

References

Adams, J. H., Mueller, I. (2017). The biology of plasmodium vivax. Cold Spring Harb. Perspect. Med. 7, a025585. doi: 10.1101/cshperspect.a025585

PubMed Abstract | CrossRef Full Text | Google Scholar

Adekunle, A. I., Pinkevych, M., McGready, R., Luxemburger, C., White, L. J., Nosten, F., et al. (2015). Modeling the dynamics of plasmodium vivax infection and hypnozoite reactivation In vivo. PloS Negl. Trop. Dis. 9, e0003595. doi: 10.1371/JOURNAL.PNTD.0003595

PubMed Abstract | CrossRef Full Text | Google Scholar

Afriat, A., Zuzarte-Luís, V., Halpern, K. B., Buchauer, L., Marques, S., Lahree, A., et al. (2021). A spatiotemporally resolved single cell atlas of the plasmodium liver stage. bioRxiv. doi: 10.1101/2021.12.03.471111

CrossRef Full Text | Google Scholar

Albuquerque, S. S., Carret, C., Grosso, A. R., Tarun, A. S., Peng, X., Kappe, S. H. I., et al. (2009). Host cell transcriptional profiling during malaria liver stage infection reveals a coordinated and sequential set of biological events. BMC Genomics 10, 1–13. doi: 10.1186/1471-2164-10-270

PubMed Abstract | CrossRef Full Text | Google Scholar

Aldridge, S., Teichmann, S. A. (2020). Single cell transcriptomics comes of age. Nat. Commun. 11, 1 11:1–1 11:4. doi: 10.1038/s41467-020-18158-5

CrossRef Full Text | Google Scholar

Andrews, S. (2010). FastQC: A quality control tool for high throughput sequence data. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc

Google Scholar

Balu, B., Maher, S. P., Pance, A., Chauhan, C., Naumov, A v., Andrews, R. M., et al. (2011). CCR4-associated factor 1 coordinates the expression of plasmodium falciparum egress and invasion proteins. Eukaryotic Cell 10, 1257. doi: 10.1128/EC.05099-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Battle, K. E., Karhunen, M. S., Bhatt, S., Gething, P. W., Howes, R. E., Golding, N., et al. (2014). Geographical variation in plasmodium vivax relapse. Malaria J. 13, 1–16. doi: 10.1186/1475-2875-13-144

CrossRef Full Text | Google Scholar

Baumgarten, S., Bryant, J. M., Sinha, A., Reyser, T., Preiser, P. R., Dedon, P. C., et al. (2019). Transcriptome-wide dynamics of extensive m6A mRNA methylation during plasmodium falciparum blood-stage development. Nat. Microbiol. 4, 2246. doi: 10.1038/S41564-019-0521-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindschedler, A., Schmuckli-Maurer, J., Wacker, R., Kramer, N., Rehmann, R., Caldelari, R., et al. (2022). Plasmodium berghei-mediated NRF2 activation in infected hepatocytes enhances parasite survival. Cell. Microbiol. 2022, 1–17. doi: 10.1155/2022/7647976

CrossRef Full Text | Google Scholar

Bissell, D. M., Hammaker, L., Schmid, R. (1972). Hemoglobin and erythrocyte catabolism in rat liver: The separate roles of parenchymal and sinusoidal cells. Blood 40, 812–822. doi: 10.1182/BLOOD.V40.6.812.812

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogale, H. N., Pascini, T v., Kanatani, S., JM, Sá, Wellems, T. E., Sinnis, P., et al. (2021). Transcriptional heterogeneity and tightly regulated changes in gene expression during plasmodium berghei sporozoite development. Proc. Natl. Acad. Sci. 118, e2023438118. doi: 10.1073/pnas.2023438118

CrossRef Full Text | Google Scholar

Bray, N. L., Pimentel, H., Melsted, P., Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34, 525–527. doi: 10.1038/NBT.3519

PubMed Abstract | CrossRef Full Text | Google Scholar

Bunnik, E. M., Batugedara, G., Saraf, A., Prudhomme, J., Florens, L., le Roch, K. G. (2016). The mRNA-bound proteome of the human malaria parasite plasmodium falciparum. Genome Biol. 17. doi: 10.1186/S13059-016-1014-0

CrossRef Full Text | Google Scholar

Byoung-Kuk, N. A., Shenai, B. R., Sijwali, P. S., Choe, Y., Pandey, K. C., Singh, A., et al. (2004). Identification and biochemical characterization of vivapains, cysteine proteases of the malaria parasite plasmodium vivax. Biochem. J. 378, 529. doi: 10.1042/BJ20031487

PubMed Abstract | CrossRef Full Text | Google Scholar

Caldelari, R., Dogga, S., Schmid, M. W., Franke-Fayard, B., Janse, C. J., Soldati-Favre, D., et al. (2019). Transcriptome analysis of plasmodium berghei during exo-erythrocytic development. Malaria J. 18, 330. doi: 10.1186/s12936-019-2968-7

CrossRef Full Text | Google Scholar

Chakravarty, S., Cockburn, I. A., Kuk, S., Overstreet, M. G., Sacci, J. B., Zavala, F. (2007). CD8+ T lymphocytes protective against malaria liver stages are primed in skin-draining lymph nodes. Nat. Med. 13, 9 13:1035–1041. doi: 10.1038/nm1628

CrossRef Full Text | Google Scholar

Chen, F., Mackey, A. J., Stoeckert, C. J., Roos, D. S. (2006). OrthoMCL-DB: querying a comprehensive multi-species collection of ortholog groups. Nucleic Acids Res. 34, D363–D368. doi: 10.1093/NAR/GKJ123

PubMed Abstract | CrossRef Full Text | Google Scholar

Chua, A. C. Y., Jie, J., Ong, Y., Malleret, B., Suwanarusk, R., Kosaisavee, V., et al. (2019). Robust continuous in vitro culture of the plasmodium cynomolgi erythrocytic stages. Nat. Commun. 10, 1–13. doi: 10.1038/s41467-019-11332-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cogswell, F. B. (1992). The hypnozoite and relapse in primate malaria. Clin. Microbiol. Rev. 5, 26–35. doi: 10.1128/CMR.5.1.26

PubMed Abstract | CrossRef Full Text | Google Scholar

Commons, R. J., Simpson, J. A., Watson, J., White, N. J., Price, R. N. (2020). Estimating the proportion of plasmodium vivax recurrences caused by relapse: A systematic review and meta-analysis. Am. J. Trop. Med. Hyg. 103, 1094–1099. doi: 10.4269/AJTMH.20-0186

PubMed Abstract | CrossRef Full Text | Google Scholar

Cubi, R., Vembar, S. S., Biton, A., Franetich, J.-F., Bordessoulles, M., Sossau, D., et al. (2017). Laser capture microdissection enables transcriptomic analysis of dividing and quiescent liver stages of plasmodium relapsing species. Cell. Microbiol. 19, e12735. doi: 10.1111/cmi.12735

CrossRef Full Text | Google Scholar

Dame, J. B., Yowell, C. A., Omara-Opyene, L., Carlton, J. M., Cooper, R. A., Li, T. (2003). Plasmepsin 4, the food vacuole aspartic proteinase found in all plasmodium spp. infecting man. Mol. Biochem. Parasitol. 130, 1–12. doi: 10.1016/S0166-6851(03)00137-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Decker, C. J., Parker, R.. (2012). P-bodies and stress granules: Possible roles in the control of translation and mRNA degradation. Cold Spring Harbor Perspectives in Biology 49, a012286. doi: 10.1101/cshperspect.a012286

CrossRef Full Text | Google Scholar

di Cristina, M., Dou, Z., Lunghi, M., Kannan, G., Huynh, M.-H., McGovern, O. L., et al. (2017). Toxoplasma depends on lysosomal consumption of autophagosomes for persistent infection. Nat Microbiol, 1–7. doi: 10.1038/nmicrobiol.2017.96

CrossRef Full Text | Google Scholar

Díez-Dacal, B., Gayarre, J., Gharbi, S., Timms, J. F., Coderch, C., Gago, F., et al. (2011). Identification of aldo-keto reductase AKR1B10 as a selective target for modification and inhibition by prostaglandin A(1): implications for antitumoral activity. Cancer Res. 71, 4161–4171. doi: 10.1158/0008-5472.CAN-10-3816

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Dou, C., Mcgovern, Z. L., Di Cristina, D., Carruthers, M. B. (2014). Toxoplasma gondii ingests and digests host cytosolic proteins. mBio 5, 1188–1202. doi: 10.1128/mBio.01188

CrossRef Full Text | Google Scholar

Drew, M. E., Banerjee, R., Uffman, E. W., Gilbertson, S., Rosenthal, P. J., Goldberg, D. E. (2008). Plasmodium food vacuole plasmepsins are activated by falcipains. J. Biol. Chem. 283, 12870. doi: 10.1074/JBC.M708949200

PubMed Abstract | CrossRef Full Text | Google Scholar

Feingold, E., Pachter, L. (2004). The ENCODE (ENCyclopedia of DNA elements) project. Science 306, 636–640. doi: 10.1126/SCIENCE.1105136/SUPPL_FILE/FEINGOLD.SOM.DC2.PDF

PubMed Abstract | CrossRef Full Text | Google Scholar

Finak, G, McDavid, A, Yajima, M, Deng, J, Gersuk, V, Shalek, AK, et al (2015). MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biology 16, 1–13. doi: 10.1186/s13059-015-0844-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Foth, B. J., Zhang, N., Chaal, B. K., Sze, S. K., Preiser, P. R., Bozdech, Z. (2011). Quantitative time-course profiling of parasite and host cell proteins in the human malaria parasite plasmodium falciparum. Mol. Cell. Proteomics 10, M110.006411. doi: 10.1074/mcp.M110.006411

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukumoto, S. I., Yamauchi, N., Moriguchi, H., Hippo, Y., Watanabe, A., Shibahara, J., et al. (2005). Overexpression of the aldo-keto reductase family protein AKR1B10 is highly correlated with smokers’ non-small cell lung carcinomas. Clin. Cancer Res. 11, 1776–1785. doi: 10.1158/1078-0432.CCR-04-1238

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstberger, S., Hafner, M., Tuschl, T. (2014). A census of human RNA-binding proteins. Nat. Rev. Genet. 15 (12), 15:829–15:845. doi: 10.1038/nrg3813

CrossRef Full Text | Google Scholar

Govaere, O., Cockell, S., Tiniakos, D., Queen, R., Younes, R., Vacca, M., et al. (2020). Transcriptomic profiling across the nonalcoholic fatty liver disease spectrum reveals gene signatures for steatohepatitis and fibrosis. Sci. Transl. Med. 12, eaba4448. doi: 10.1126/SCITRANSLMED.ABA4448

PubMed Abstract | CrossRef Full Text | Google Scholar

Guallar, D., Wang, J. (2014). RNA-Binding proteins in pluripotency, differentiation, and reprogramming. Front. Biol. 9, 389–409. doi: 10.1007/s11515-014-1326-y

CrossRef Full Text | Google Scholar

Gupta, D. K., Dembele, L., Voorberg-Van Der Wel, A., Roma, G., Yip, A., Chuenchob, V., et al. (2019). The plasmodium liver-specific protein 2 (LISP2) is an early marker of liver stage development. Elife 8, :e43362. doi: 10.7554/ELIFE.43362

PubMed Abstract | CrossRef Full Text | Google Scholar

Gural, N., Mancio-Silva, L., Miller, A. B., Galstian, A., Butty, V. L., Levine, S. S., et al. (2018). In vitro culture, drug sensitivity, and transcriptome of plasmodium vivax hypnozoites. Cell Host Microbe 23, 395–406.e4. doi: 10.1016/j.chom.2018.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hafemeister, C., Satija, R. (2019). Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 1–15. doi: 10.1186/S13059-019-1874-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, J. D., Dinkova-Kostova, A. T. (2014). The Nrf2 regulatory network provides an interface between redox and intermediary metabolism. Trends Biochem. Sci. 39, 199–218. doi: 10.1016/J.TIBS.2014.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hentze, M. W., Castello, A., Schwarzl, T., Preiss, T. (2018). A brave new world of RNA-binding proteins. Nat. Rev. Mol. Cell Biol. 19, 5 19:327–341. doi: 10.1038/nrm.2017.130

CrossRef Full Text | Google Scholar

Hollingdale, M. R., Collins, W. E., Campbell, C., Schwartz, A. L. (1985). In vitro culture of two populations (dividing and nondividing) of exoerythrocytic parasites of plasmodium vivax. Am. J. Trop. Med. Hyg. 34, 216–222. doi: 10.4269/AJTMH.1985.34.216

PubMed Abstract | CrossRef Full Text | Google Scholar

Howick, V. M., Russell, A. J. C., Andrews, T., Heaton, H., Reid, A. J., Natarajan, K., et al. (2019). The malaria cell atlas: Single parasite transcriptomes across the complete plasmodium life cycle. Science 365, eaaw2619. doi: 10.1126/science.aaw2619

PubMed Abstract | CrossRef Full Text | Google Scholar

Huldén, L., Huldén, L., Heliövaara, K. (2008). Natural relapses in vivax malaria induced by anopheles mosquitoes. Malaria J. 7, 1–11. doi: 10.1186/1475-2875-7-64

CrossRef Full Text | Google Scholar

Itoe, M. A., Sampaio, J. L., Cabal, G. G., Real, E., Zuzarte-Luis, V., March, S., et al. (2014). Host cell phosphatidylcholine is a key mediator of malaria parasite survival during liver stage infection. Cell Host Microbe 16, 778. doi: 10.1016/J.CHOM.2014.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, J., Krishack, P., Cao, D. (2006). Role of aldo-keto reductases in development of prostate and breast cancer. Front. Biosci. 11, 2767–2773. doi: 10.2741/2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, Y. C., Dema, B., Rodriguez-Garcia, R., López-Camacho, C., Leoratti, F. M. S., Lall, A., et al. (2020). Evaluation of chimpanzee adenovirus and MVA expressing TRAP and CSP from plasmodium cynomolgi to prevent malaria relapse in nonhuman primates. Vaccines (Basel) 8, 1–14. doi: 10.3390/VACCINES8030363

CrossRef Full Text | Google Scholar

Kivioja, T., Vähärautio, A., Karlsson, K., Bonke, M., Enge, M., Linnarsson, S., et al. (2011). Counting absolute numbers of molecules using unique molecular identifiers. Nat. Methods 9, 72–74. doi: 10.1038/nmeth.1778

PubMed Abstract | CrossRef Full Text | Google Scholar

Klein, A. M., Mazutis, L., Akartuna, I., Tallapragada, N., Veres, A., Li, V., et al. (2015). Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell 161, 1187–1201. doi: 10.1016/j.cell.2015.04.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Kortuem, K. M., Stewart, A. K. (2013). Carfilzomib. Blood 121, 893–897. doi: 10.1182/BLOOD-2012-10-459883

PubMed Abstract | CrossRef Full Text | Google Scholar

Krotoski, W. A., Collins, W. E., Bray, R. S., Garnham, P. C., Cogswell, F. B., Gwadz, R. W., et al. (1982). Demonstration of hypnozoites in sporozoite-transmitted plasmodium vivax infection. Am. J. Trop. Med. Hyg. 31, 1291–1293. doi: 10.4269/AJTMH.1982.31.1291

PubMed Abstract | CrossRef Full Text | Google Scholar

Krotoski, W. A., Garnham, P. C. C., Cogswell, F. B., Collins, W. E., Bray, R. S., Gwasz, R. W., et al. (1986). Observations on early and late post-sporozoite tissue stages in primate malaria. IV. pre-erythrocytic schizonts and/or hypnozoites of chesson and north Korean strains of plasmodium vivax in the chimpanzee. Am. J. Trop. Med. Hyg. 35, 263–274. doi: 10.4269/AJTMH.1986.35.263

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucera, M., Isserlin, R., Arkhangorodsky, A., Bader, G. D. (2016). AutoAnnotate: A cytoscape app for summarizing networks with semantic annotations. F1000Res 5, :1717. doi: 10.12688/F1000RESEARCH.9090.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuleshov, M. V., Jones, M. R., Rouillard, A. D., Fernandez, N. F., Duan, Q., Wang, Z., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97. doi: 10.1093/NAR/GKW377

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachmann, A., Xu, H., Krishnan, J., Berger, S. I., Mazloom, A. R., Ma’ayan, A. (2010). ChEA: transcription factor regulation inferred from integrating genome-wide ChIP-X experiments. Bioinformatics 26, 2438–2444. doi: 10.1093/BIOINFORMATICS/BTQ466

PubMed Abstract | CrossRef Full Text | Google Scholar

LaMonte, G. M., Orjuela-Sanchez, P., Calla, J., Wang, L. T., Li, S., Swann, J., et al. (2019). Dual RNA-seq identifies human mucosal immunity protein mucin-13 as a hallmark of plasmodium exoerythrocytic infection. Nat. Commun. 10, 488. doi: 10.1038/s41467-019-08349-0

PubMed Abstract | CrossRef Full Text | Google Scholar

le Roch, K. G., Johnson, J. R., Florens, L., Zhou, Y., Santrosyan, A., Grainger, M., et al. (2004). Global analysis of transcript and protein levels across the plasmodium falciparum life cycle. Genome Res. 14, 2308–2318. doi: 10.1101/GR.2523904

PubMed Abstract | CrossRef Full Text | Google Scholar

Liehl, P., Mota, M. M. (2012). Innate recognition of malarial parasites by mammalian hosts. Int. J. Parasitol. 42, 557–566. doi: 10.1016/J.IJPARA.2012.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Liehl, P., Zuzarte-Luís, V., Chan, J., Zillinger, T., Baptista, F., Carapau, D., et al. (2013). Host-cell sensors for plasmodium activate innate immunity against liver-stage infection. Nat. Med. 20 (1), 20:47–20:53. doi: 10.1038/nm.3424

CrossRef Full Text | Google Scholar

Lun, A. T. L., Riesenfeld, S., Andrews, T., Dao, T. P., Gomes, T., Marioni, J. C. (2019). EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data. Genome Biol. 20, 63. doi: 10.1186/s13059-019-1662-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Lysenko, A. J., Beljaev, A. E., Rybalka, V. M. (1977). Population studies of plasmodium vivax: 1. the theory of polymorphism of sporozoites and epidemiological phenomena of tertian malaria. Bull. World Health Organ 55, 541.

PubMed Abstract | Google Scholar

Macosko, E. Z., Basu, A., Satija, R., Nemesh, J., Shekhar, K., Goldman, M., et al. (2015). Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell 161, 1202–1214. doi: 10.1016/J.CELL.2015.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

MacParland, S. A., Liu, J. C., Ma, X.-Z., Innes, B. T., Bartczak, A. M., Gage, B. K., et al. (2018). Single cell RNA sequencing of human liver reveals distinct intrahepatic macrophage populations. Nat. Commun. 9, 4383. doi: 10.1038/s41467-018-06318-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, S. P., Conway, A. J., Roth, A., Adapa, S. R., Cualing, P., Andolina, C., et al (2020). An adaptable soft-mold embossing process for fabricating optically-accessible, microfeature-based culture systems and application toward liver stage antimalarial compound testing. Lab. Chip 20, 1124–39. doi: 10.1039/C9LC00921C

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, S. P., Vantaux, A., Chaumeau, V., Chua, A. C. Y., Cooper, C. A., Andolina, C., et al. (2021b). Probing the distinct chemosensitivity of plasmodium vivax liver stage parasites and demonstration of 8-aminoquinoline radical cure activity in vitro. Sci. Rep. 11 (1), 1–18. doi: 10.1038/S41598-021-99152-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, S., Vantaux, A., Cooper, C., Chasen, N., Cheng, W., Joyner, C., et al. (2021a). A phenotypic screen for the liver stages of plasmodium vivax. Bio Protoc 11, e4253. doi: 10.21769/BIOPROTOC.4253

PubMed Abstract | CrossRef Full Text | Google Scholar

Mancio-Silva, L., Gural, N., Real, E., Wadsworth, M. H., Butty, V. L., March, S., et al. (2022). A single-cell liver atlas of plasmodium vivax infection. Cell Host Microbe 30, 1048–1060.e5. doi: 10.1016/J.CHOM.2022.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, H. J., Maser, E. (2009). Role of human aldo–keto-reductase AKR1B10 in the protection against toxic aldehydes. Chemico-Biological Interact. 178, 145–150. doi: 10.1016/J.CBI.2008.10.021

CrossRef Full Text | Google Scholar

Melsted, P., Booeshaghi, A. S., Liu, L., Gao, F., Lu, L., Min, K. H. J., et al. (2021). Modular, efficient and constant-memory single-cell RNA-seq preprocessing. Nat. Biotechnol. 39, 7 39:813–818. doi: 10.1038/s41587-021-00870-2

CrossRef Full Text | Google Scholar

Merico, D., Isserlin, R., Stueker, O., Emili, A., Bader, G. D. (2010). Enrichment map: A network-based method for gene-set enrichment visualization and interpretation. PloS One 5, e13984. doi: 10.1371/JOURNAL.PONE.0013984

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, J., Li, J., Fan, Q., Xiaolian, Li, Xinyi, Li, Cui, L. (2010). The puf-family RNA-binding protein PfPuf2 regulates sexual development and sex differentiation in the malaria parasite plasmodium falciparum. J. Cell Sci. 123, 1039. doi: 10.1242/JCS.059824

PubMed Abstract | CrossRef Full Text | Google Scholar

Mikolajczak, S. A., Vaughan, A. M., Kangwanrangsan, N., Roobsoong, W., Fishbaugher, M., Yimamnuaychok, N., et al. (2015). Plasmodium vivax liver stage development and hypnozoite persistence in human liver-chimeric mice. Cell Host Microbe 17, 526–535. doi: 10.1016/J.CHOM.2015.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, J. L., Sack, B. K., Baldwin, M., Vaughan, A. M., Kappe, S. H. I. (2014). Interferon-mediated innate immune responses against malaria parasite liver stages. Cell Rep. 7, 436–447. doi: 10.1016/J.CELREP.2014.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, I., Jex, A. R., Kappe, S. H. I., Mikolajczak, S. A., Sattabongkot, J., Patrapuvich, R., et al. (2019). Transcriptome and histone epigenome of plasmodium vivax salivary-gland sporozoites point to tight regulatory control and mechanisms for liver-stage differentiation in relapsing malaria. Int. J. Parasitol. 49, 501–513. doi: 10.1016/J.IJPARA.2019.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Neefjes, J., Jongsma, M. L. M., Paul, P., Bakke, O. (2011). Towards a systems understanding of MHC class I and MHC class II antigen presentation. Nat. Rev. Immunol. 11 (12), 11:823–11:836. doi: 10.1038/nri3084

CrossRef Full Text | Google Scholar

Paquet, T., le Manach, C., Cabrera, D. G., Younis, Y., Henrich, P. P., Abraham, T. S., et al. (2017). Antimalarial efficacy of MMV390048, an inhibitor of plasmodium phosphatidylinositol 4-kinase. Sci. Trans. Med. 9 (387):eaad9735. doi: 10.1126/SCITRANSLMED.AAD9735

CrossRef Full Text | Google Scholar

Parussini, F., Coppens, I., Shah, P. P., Diamond, S. L., Carruthers, V. B. (2010). Cathepsin l occupies a vacuolar compartment and is a protein maturase within the Endo/Exocytic system of toxoplasma gondii. Mol. Microbiol. 76, 1340. doi: 10.1111/J.1365-2958.2010.07181.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Penaranda, C., Hung, D. T. (2019). Single-cell RNA sequencing to understand host–pathogen interactions. ACS Infect. Dis. 5, 336–344. doi: 10.1021/ACSINFECDIS.8B00369

PubMed Abstract | CrossRef Full Text | Google Scholar

Poran, A., Nötzel, C., Aly, O., Mencia-Trinchant, N., Harris, C. T., Guzman, M. L., et al. (2017). Single-cell RNA sequencing reveals a signature of sexual commitment in malaria parasites. Nature 551, 95–99. doi: 10.1038/nature24280

PubMed Abstract | CrossRef Full Text | Google Scholar

Posfai, D., Maher, S. P., Roesch, C., Vantaux, A., Sylvester, K., Péneau, J., et al. (2020). Plasmodium vivax liver and blood stages recruit the druggable host membrane channel aquaporin-3. Cell Chem. Biol 27:719–727.e5. . doi: 10.1016/j.chembiol.2020.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, J., Li, J., Zhang, Y., He, R., Liu, X., Gong, K., et al. (2021). AKR1B10 promotes breast cancer cell proliferation and migration via the PI3K/AKT/NF-κB signaling pathway. Cell Biosci. 11, 163. doi: 10.1186/S13578-021-00677-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists, (2019 update). Nucleic Acids Res. 47, W191–W198. doi: 10.1093/NAR/GKZ369

PubMed Abstract | CrossRef Full Text | Google Scholar

Real, E., Howick, V. M., Dahalan, F. A., Witmer, K., Cudini, J., Andradi-Brown, C., et al. (2021). A single-cell atlas of plasmodium falciparum transmission through the mosquito. Nat. Commun. 12 (1), 1–13. doi: 10.1038/s41467-021-23434-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Reid, A. J., Talman, A. M., Bennett, H. M., Gomes, A. R., Sanders, M. J., Illingworth, C. J. R., et al. (2018). Single-cell RNA-seq reveals hidden transcriptional variation in malaria parasites. Elife 7, 1–29. doi: 10.7554/eLife.33105

CrossRef Full Text | Google Scholar

Ren, X., Wen, W., Fan, X., Hou, W., Su, B., Cai, P., et al. (2021). COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas. Cell 184, 1895–1913.e19. doi: 10.1016/J.CELL.2021.01.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, L. J., Wampfler, R., Betuela, I., Karl, S., White, M. T., Li Wai Suen, C. S. N., et al. (2015). Strategies for understanding and reducing the plasmodium vivax and plasmodium ovale hypnozoite reservoir in Papua new guinean children: a randomised placebo-controlled trial and mathematical model. PloS Med. 12:e1001891. doi: 10.1371/JOURNAL.PMED.1001891

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenthal, P. J. (2011). Falcipains and other cysteine proteases of malaria parasites. Adv. Exp. Med. Biol. 712, 30–48. doi: 10.1007/978-1-4419-8414-2_3

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, A., Maher, S. P., Conway, A. J., Ubalee, R., Chaumeau, V., Andolina, C., et al. (2018). A comprehensive model for assessment of liver stage therapies targeting plasmodium vivax and plasmodium falciparum. Nat. Commun. 9, 1837. doi: 10.1038/s41467-018-04221-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruberto, A. A., Bourke, C., Merienne, N., Obadia, T., Amino, R., Mueller, I. (2021). Single-cell RNA sequencing reveals developmental heterogeneity among plasmodium berghei sporozoites. Sci. Rep. 11, 4127. doi: 10.1038/s41598-021-82914-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruberto, A. A., Bourke, C., Vantaux, A., Maher, S. P., Jex, A., Witkowski, B., et al. (2022). Single-cell RNA sequencing of plasmodium vivax sporozoites reveals stage- and species-specific transcriptomic signatures. PLOS Neglected Tropical Diseases 16, :e0010633. doi: 10.1371/JOURNAL.PNTD.0010633

PubMed Abstract | CrossRef Full Text | Google Scholar

Sà, J. M., Cannon, M v., Caleon, R. L., Wellems, T. E., Serre, D. (2020). Single-cell transcription analysis of plasmodium vivax blood-stage parasites identifies stage- and species-specific profiles of expression. PloS Biol. 18, e3000711. doi: 10.1371/journal.pbio.3000711

PubMed Abstract | CrossRef Full Text | Google Scholar

Salas, F., Fichmann, J., Lee, G. K., Scott, M. D., Rosenthal, P. J. (1995). Functional expression of falcipain, a plasmodium falciparum cysteine proteinase, supports its role as a malarial hemoglobinase. Infection Immun. 63, 2120. doi: 10.1128/iai.63.6.2120-2125.1995

CrossRef Full Text | Google Scholar

Sattabongkot, J., Yimamnuaychoke, N., Leelaudomlipi, S., Rasameesoraj, M., Jenwithisuk, R., Coleman, R. E., et al. (2006). Establishment of a human hepatocyte line that supports in vitro development of the exo-erythrocytic stages of the malaria parasites plasmodium falciparum and p. vivax. Am. J. Trop. Med. Hyg. 74, 708–715. doi: 10.4269/ajtmh.2006.74.708

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarpulla, R. C. (2002). Transcriptional activators and coactivators in the nuclear control of mitochondrial function in mammalian cells. Gene 286, 81–89. doi: 10.1016/S0378-1119(01)00809-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Schäfer, C., Dambrauskas, N., Steel, R. W., Carbonetti, S., Chuenchob, V., Flannery, E. L., et al. (2018). A recombinant antibody against plasmodium vivax UIS4 for distinguishing replicating from dormant liver stages. Malaria J. 17, 370. doi: 10.1186/S12936-018-2519-7

CrossRef Full Text | Google Scholar

Schäfer, C., Roobsoong, W., Sattabongkot, J., Mikolajczak, S. A., Kappe, S. H. I. (2020). A humanized mouse model for plasmodium vivax to test interventions that block liver stage to blood stage transition and blood stage infection. iScience 23, 101381. doi: 10.1016/j.isci.2020.101381

PubMed Abstract | CrossRef Full Text | Google Scholar

Schäfer, C., Zanghi, G., Vaughan, A. M., Kappe, S. H. I. (2021). Plasmodium vivax latent liver stage infection and relapse: Biological insights and new experimental tools. Annual Review of Microbiology, 87–106. doi: 10.1146/annurev-micro-032421-061155

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, S v., Chappell, L., Hostetler, J. B., Amaratunga, C., Suon, S., Böhme, U., et al. (2020). Analysis of plasmodium vivax schizont transcriptomes from field isolates reveals heterogeneity of expression of genes involved in host-parasite interactions. Sci. Rep. 10, 16667. doi: 10.1038/s41598-020-73562-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, A. P., Buscaglia, C. A., Wang, Q., Levay, A., Nussenzweig, D. R., Walker, J. R., et al. (2007). Plasmodium circumsporozoite protein promotes the development of the liver stages of the parasite. Cell 131, 492–504. doi: 10.1016/J.CELL.2007.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M., et al. (2019). Comprehensive integration of single-cell data. Cell 177, 1888–1902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Szabo, P. A., Dogra, P., Gray, J. I., Wells, S. B., Connors, T. J., Weisberg, S. P., et al. (2021). Longitudinal profiling of respiratory and systemic immune responses reveals myeloid cell-driven lung inflammation in severe COVID-19. Immunity 54, 797–814.e6. doi: 10.1016/J.IMMUNI.2021.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarun, A. S., Peng, X., Dumpit, R. F., Ogata, Y., Silva-Rivera, H., Camargo, N., et al. (2008). A combined transcriptome and proteome survey of malaria parasite liver stages. Proc. Natl. Acad. Sci. 105, 305–310. doi: 10.1073/PNAS.0710780104

CrossRef Full Text | Google Scholar

Tebay, L. E., Robertson, H., Durant, S. T., Vitale, S. R., Penning, T. M., Dinkova-Kostova, A. T., et al. (2015). Mechanisms of activation of the transcription factor Nrf2 by redox stressors, nutrient cues, and energy status and the pathways through which it attenuates degenerative disease. Free Radic. Biol. Med. 88, 108. doi: 10.1016/J.FREERADBIOMED.2015.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Tharyan, R. G., Annibal, A., Schiffer, I., Laboy, R., Atanassov, I., Weber, A. L., et al. (2020). NFYB-1 regulates mitochondrial function and longevity via lysosomal prosaposin. Nat. Metab. 2, 387–396. doi: 10.1038/s42255-020-0200-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Toro-Moreno, M., Sylvester, K., Srivastava, T., Posfai, D., Derbyshire, E. R. (2020). RNA-Seq analysis illuminates the early stages of Plasmodium liver infection. mBio. doi: 10.1128/mBio.03234-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Traag, V. A., Waltman, L., van Eck, N. J. (2019). From louvain to Leiden: guaranteeing well-connected communities. Sci. Rep. 9, 1–12. doi: 10.1038/s41598-019-41695-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, H. H., Lai, H. Y., Chen, Y. C., Li, C. F., Huang, H. S., Liu, H. S., et al. (2017). Metformin promotes apoptosis in hepatocellular carcinoma through the CEBPD-induced autophagy pathway. Oncotarget 8, 13832. doi: 10.18632/ONCOTARGET.14640

PubMed Abstract | CrossRef Full Text | Google Scholar

Tucker, M., Valencia-Sanchez, M. A., Staples, R. R., Chen, J., Denis, C. L., Parker, R. (2001). The transcription factor associated Ccr4 and Caf1 proteins are components of the major cytoplasmic mRNA deadenylase in saccharomyces cerevisiae. Cell 104, 377–386. doi: 10.1016/S0092-8674(01)00225-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ungureanu, E., Killick-Kendrick, R., Garnham, P. C. C., Branzei, P., Romanescu, C., Shute, P. G. (1976). Prepatent periods of a tropical strain of plasmodium vivax after inoculations of tenfold dilutions of sporozoites. Trans. R Soc. Trop. Med. Hyg. 70, 482–483. doi: 10.1016/0035-9203(76)90133-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Vantaux, A., Péneau, J., Cooper, C. A., Kyle, D. E., Witkowski, B., Maher, S. P. (2022). Liver stage fate determination in plasmodium vivax parasites: characterization of schizont growth and hypnozoite fating from patient isolates. Front. Microbiol. 13. doi: 10.3389/fmicb.2022.976606

CrossRef Full Text | Google Scholar

Vembar, S. S., Droll, D., Scherf, A. (2016). Translational regulation in blood stages of the malaria parasite plasmodium spp.: systems-wide studies pave the way. Wiley Interdiscip. Rev. RNA 7, 772. doi: 10.1002/WRNA.1365

PubMed Abstract | CrossRef Full Text | Google Scholar

Vijayan, K., Wei, L., Glennon, E. K. K., Mattocks, C., Bourgeois, N., Staker, B., et al. (2021). Host-targeted interventions as an exciting opportunity to combat malaria. Chem. Rev. 121, 10452–10468. doi: 10.1021/acs.chemrev.1c00062

PubMed Abstract | CrossRef Full Text | Google Scholar

Voorberg-van der Wel, A., Roma, G., Gupta, D. K., Schuierer, S., Nigsch, F., Carbone, W., et al. (2017). A comparative transcriptomic analysis of replicating and dormant liver stages of the relapsing malaria parasite plasmodium cynomolgi. Elife 6, e29605. doi: 10.7554/eLife.29605.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Voorberg-van der Wel, A. M., Zeeman, A. M., Nieuwenhuis, I. G., van der Werff, N. M., Klooster, E. J., Klop, O., et al. (2020). A dual fluorescent plasmodium cynomolgi reporter line reveals in vitro malaria hypnozoite reactivation. Commun. Biol. 3. doi: 10.1038/s42003-019-0737-3

CrossRef Full Text | Google Scholar

Voronina, E., Seydoux, G., Sassone-Corsi, P., Nagamori, I. (2011). RNA Granules in germ cells. Cold Spring Harbor Perspect. Biol. 3, a002774. doi: 10.1101/CSHPERSPECT.A002774

CrossRef Full Text | Google Scholar

Wang, S. M., Lin, H. Y., Chen, Y. L., Hsu, T. I., Chuang, J. Y., Kao, T. J., et al. (2019). CCAAT/enhancer-binding protein delta regulates the stemness of glioma stem-like cells through activating PDGFA expression upon inflammatory stimulation. J. Neuroinflamm. 16, 1–12. doi: 10.1186/S12974-019-1535-Z

CrossRef Full Text | Google Scholar

Wang, C., Yan, R., Luo, D., Watabe, K., Liao, D.-F., Cao, D. (2009). Aldo-keto reductase family 1 member B10 promotes cell survival by regulating lipid synthesis and eliminating carbonyls. J. Biol. Chem. 284, 26742–26748. doi: 10.1074/jbc.M109.022897

PubMed Abstract | CrossRef Full Text | Google Scholar

White, N. J. (2011). Determinants of relapse periodicity in plasmodium vivax malaria. Malaria J. 10, 297. doi: 10.1186/1475-2875-10-297

CrossRef Full Text | Google Scholar

White, M., Amino, R., Mueller, I. (2017). Theoretical implications of a pre-erythrocytic plasmodium vivax vaccine for preventing relapses. Trends Parasitol. 33, 260. doi: 10.1016/J.PT.2016.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

White, M. T., Walker, P., Karl, S., Hetzel, M. W., Freeman, T., Waltmann, A., et al. (2018). Mathematical modelling of the impact of expanding levels of malaria control interventions on plasmodium vivax. Nat. Commun. 9 (1), 1–10. doi: 10.1038/s41467-018-05860-8

PubMed Abstract | CrossRef Full Text | Google Scholar

World Health Organization (2021). “World malaria report 2021,” in World malaria report 2021. Geneva:World Health Organization; 2021. Licence: CC BY-NC-SA 3.0 IGO. Available at: https://www.who.int/publications/i/item/9789240040496.

Google Scholar

Yokomizo, T., Izumi, T., Takahashi, T., Kasamal, T., Kobayashis, Y., Satos, F., et al. (1993). Enzymatic inactivation the porcine kidney O f leukotriene B4 by a novel enzyme found in PURIFICATION AND PROPERTIES OF LEUKOTRIENE bq 12-HYDROXYDEHYDROGENASE*. THE J. OF Biol. Chem. 268, 18128–18135. doi: 10.1016/S0021-9258(17)46820-7

CrossRef Full Text | Google Scholar

Youn, J. Y., Dyakov, B. J. A., Zhang, J., Knight, J. D. R., Vernon, R. M., Forman-Kay, J. D., et al. (2019). Properties of stress granule and p-body proteomes. Mol. Cell 76, 286–294. doi: 10.1016/J.MOLCEL.2019.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, G. X. Y., Terry, J. M., Belgrader, P., Ryvkin, P., Bent, Z. W., Wilson, R., et al. (2017). Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049. doi: 10.1038/ncomms14049

PubMed Abstract | CrossRef Full Text | Google Scholar