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

Anthony Ruberto, Steven Maher, Amélie Vantaux, Chester Joyner, Caitlin Bourke, Balu Balan, Aaron Jex, Ivo Mueller, Benoit Witkowski, Dennis Kyle

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anthony Ruberto, Steven Maher, Amélie Vantaux, Chester Joyner, Caitlin Bourke, et al. Single-cell RNA profiling of Plasmodium vivax-infected hepatocytes reveals parasite- and host-specific transcriptomic signatures and therapeutic targets. Frontiers in Cellular and Infection Microbiology, 2022, 12, pp.986314. 10.3389/fcimb.2022.986314. [hal-04697382](https://hal.science/hal-04697382v1)

HAL Id: hal-04697382

https://hal.science/hal-04697382v1

Submitted on 24 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Xinzhuan Su,
National Institutes of Health (NIH),
United States

REVIEWED BY
Ian Cheeseman,
Texas Biomedical Research Institute,
United States
Gigliola Zanghi,
Seattle Children's Research Institute,
United States

*CORRESPONDENCE
Dennis E. Kyle
dennis.kyle@uga.edu

SPECIALTY SECTION
This article was submitted to
Parasite and Host,
a section of the journal
Frontiers in Cellular and
Infection Microbiology

RECEIVED 04 July 2022
ACCEPTED 08 August 2022
PUBLISHED 25 August 2022

CITATION
Ruberto AA, Maher SP, Vantaux A, Joyner CJ, Bourke C, Balan B, Jex A, Mueller I, Witkowski B and Kyle DE (2022) Single-cell RNA profiling of *Plasmodium vivax*-infected hepatocytes reveals parasite- and host- specific transcriptomic signatures and therapeutic targets. *Front. Cell. Infect. Microbiol.* 12:986314. doi: 10.3389/fcimb.2022.986314

COPYRIGHT
© 2022 Ruberto, Maher, Vantaux, Joyner, Bourke, Balan, Jex, Mueller, Witkowski and Kyle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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

Anthony A. Ruberto¹, Steven P. Maher¹, Amélie Vantaux², Chester J. Joyner¹,³, Caitlin Bourke⁴,⁵, Balu Balan⁴,⁶, Aaron Jex⁴,⁵,⁶, Ivo Mueller⁴,⁵, Benoit Witkowski² and Dennis E. Kyle¹*

¹Center for Tropical and Emerging Global Diseases, University of Georgia, Athens, GA, United States, ²Malaria Molecular Epidemiology Unit, Institut Pasteur du Cambodge, Phnom Penh, Cambodia, ³Department of Infectious Diseases, College of Veterinary Medicine, University of Georgia, Athens, GA, United States, ⁴Population Health & Immunity Division, Walter and Eliza Hall Institute, Parkville, VIC, Australia, ⁵Department of Medical Biology, The University of Melbourne, Parkville, VIC, Australia, ⁶Faculty of Veterinary and Agricultural Sciences, University of Melbourne, Parkville, VIC, Australia

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.

KEYWORDS

*Plasmodium vivax*, liver stage malaria, host-parasite interactions, hypnozoite, single-cell RNA sequencing (scRNAseq), 10X Genomics, primary human hepatocyte (PHH)
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* U1S4 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 рах -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 'FindIntegration Anchors' 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-*Pv*UIS4 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 *Hs*AKR1B10 or 1:400 AlexaFluor 488 goat
anti-mouse for *Pv*UIS4. 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.

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><img src="image_4.png" /><figcaption>FIGURE 1</figcaption></figure>

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.

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

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
<figure><img src="image_5.png" /></figure>

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

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

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
A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

C

<figure><img src="image_12.png" /></figure>

<table><caption>FIGURE 3</caption><thead><tr><th>Cluster</th><th>GO Term</th><th>Adj P</th></tr></thead><tbody><tr><td rowspan="4">Hypnozoites</td><td>CCR4-NOT complex</td><td>2.20e-9</td></tr><tr><td>Nucleus</td><td>1.22e-2</td></tr><tr><td>P body</td><td>2.28e-2</td></tr><tr><td>Regulation of gene expression</td><td>1.97e-6</td></tr><tr><td rowspan="4">Schizonts</td><td>Eukaryotic translation initiation factor 3 complex</td><td>5.22e-6</td></tr><tr><td>RNP complex</td><td>2.35e-5</td></tr><tr><td>Translation preinitiation complex</td><td>1.46e-4</td></tr><tr><td>Translation</td><td>7.83e-7</td></tr></tbody></table>

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) ISNE 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.

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

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
A

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

B

<figure><img src="image_15.png" /></figure>

Cluster

<figure><img src="image_16.png" /></figure>

C

PVP01_0304700::liver specific protein 2

<figure><img src="image_17.png" /></figure>

D

<figure><img src="image_18.png" /></figure>

E

<figure><img src="image_19.png" /></figure>

<table><caption>FIGURE 4</caption><thead><tr><th>Cluster</th><th>GO Term</th><th>Adj P</th></tr></thead><tbody><tr><td rowspan="3">C1</td><td>CCR4-NOT complex</td><td>1.15e-04</td></tr><tr><td>Germ plasm</td><td>2.62e-03</td></tr><tr><td>P granule</td><td>2.62e-03</td></tr><tr><td rowspan="3">C2</td><td>Regulation of RNA splicing</td><td>5.05e-05</td></tr><tr><td>Ribonucleoprotein complex</td><td>3.54e-04</td></tr><tr><td>Translation</td><td>4.83e-02</td></tr><tr><td rowspan="2">C3</td><td>Translational initiation</td><td>5.36e-02</td></tr><tr><td>RNA binding</td><td>3.01e-03</td></tr><tr><td></td><td>Nucleic acid binding</td><td>6.71e-03</td></tr></tbody></table>

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 to 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 log2FC > 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.

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

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

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 (Mérico 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
<figure><img src="image_20.png" /></figure>

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.

described as a biomarker of cancer (Fukumoto et al., 2005; Jin
et al., 2006; Diez-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
A

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

B

Differential
Expression vs

non-infected hepatocytes

Hep+schizont

Increase: 106

Decrease: 16

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

Hep+hypnozoite

Increase: 365

Decrease: 41

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

\# of genes

\# of biological processes &
pathways enriched

C

Hepatocyte+schizont

vs non-infected hepatocytes

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

D

Hepatocyte+hypnozoite

vs non-infected hepatocytes

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

E

<figure><img src="image_31.png" /></figure>

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<sub>2</sub>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.
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.

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
<figure><img src="image_32.png" /><figcaption>FIGURE 7</figcaption></figure>

Transcription factor enrichment in hepatocytes infected with schizonts or hypnozoites. Dot plot of transcription factors enriched in Hepatocyte<sup>+hypnozoite</sup> (left) or Hepatocyte<sup>+schizont</sup> 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.

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 Protocols.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 threhold = 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* (exposed).
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 PyUIS4 (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 UMs 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 UMs in schizonts obtained in the current dataset; teal dashed horizontal line: median number genes detected or UMs 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. Idh: 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, IC₅₀ 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 naïve 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 (naïve 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 naïve) 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), PyUIS4 (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., and Mueller, I. (2017). The biology of plasmodium vivax. *Cold Spring Harb. Perspect. Med.* 7, a025585. doi: 10.1101/cshperspect. a025585

Adekunle, A. I., Pinkeyvych, 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

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

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
Aldridge, S., and 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

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

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

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

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

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

Bissell, D. M., Hammaker, L., and 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

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

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

Bunnik, E. M., Batugeda, G., Saraf, A., Prudhomme, J., Florens, L., and 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

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

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

Chakravarty, S., Cockburn, I. A., Kuk, S., Overstreet, M. G., Sacci, J. B., and 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

Chen, F., Mackey, A. J., Stoeckert, C. J., and 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

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 cynomolgii erythrocytic stages. *Nat. Commun.* 10, 1–13. doi: 10.1038/s41467-019-11332-4

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

Commons, R. J., Simpson, J. A., Watson, J., White, N. J., and 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

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

Dame, J. B., Yowell, C. A., Omara-Opyene, L., Carlton, J. M., Cooper, R. A., and 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

Decker, C. J., and 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

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

Diez-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-3186

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

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

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

Feingold, E., and 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

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

Foth, B. J., Zhang, N., Chaal, B. K., Sze, S. K., Preiser, P. R., and 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

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

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

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

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

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

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

Hafemeister, C., and 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

Hayes, J. D., and 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

Hentze, M. W., Castello, A., Schwarzl, T., and 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

Hollingdale, M. R., Collins, W. E., Campbell, C., and 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

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

Huldén, L., Huldén, L., and 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

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

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

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 cynomolgii to prevent malaria relapse in nonhuman primates. *Vaccines (Basel)* 8, 1–14. doi: 10.3390/VACCINES8030363

Kiviova, 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

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
Kortuem, K. M., and Stewart, A. K. (2013). *Carfilzomib*. *Blood* 121, 893–897. doi: 10.1182/BLOOD-2012-10-459883

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

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

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

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

Lachmann, A., Xu, H., Krishnan, J., Berger, S. I., Mazloom, A. R., and 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

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

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

Liehl, P., and 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

Liehl, P., Zuzarte-Luis, 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

Lun, A. T. L., Riesenfeld, S., Andrews, T., Dao, T. P., Gomes, T., and 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

Lysenko, A. J., Beljaev, A. E., and 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.

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

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

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

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

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/BIOPROTOCOL.4253

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

Martin, H. J., and 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

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

Merico, D., Isserlin, R., Stueker, O., Emili, A., and 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

Miao, J., Li, J., Fan, Q., Xiaolian, Li, Xinyi, Li, and 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

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

Miller, J. L., Sack, B. K., Baldwin, M., Vaughan, A. M., and 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

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

Neefjes, J., Jongsma, M. L. M., Paul, P., and 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

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

Parussini, F., Coppens, I., Shah, P. P., Diamond, S. L., and Carruthers, V. B. (2010). Cathepsin J1 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

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

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

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

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

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

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

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

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

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

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

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

Ruberto, A. A., Bourke, C., Merienne, N., Obadia, T., Amino, R., and 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

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

Sà, J. M., Cannon, M. v., Caleon, R. L., Wellems, T. E., and 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
Salas, F., Fichmann, J., Lee, G. K., Scott, M. D., and 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

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

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

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

Schäfer, C., Roobsoong, W., Sattabongkot, J., Mikolajczak, S. A., and 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

Schäfer, C., Zanghi, G., Vaughan, A. M., and 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

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

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

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

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

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

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

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

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

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

Traag, V. A., Waltman, L., and 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

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 CEPBD-induced autophagy pathway. *Oncotarget* 8, 13832. doi: 10.18632/ONCOTARGET.14640

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

Ungureanu, E., Killick-Kendrick, R., Garnham, P. C. C., Branzei, P., Romanescu, C., and 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

Vantaux, A., Pèneau, J., Cooper, C. A., Kyle, D. E., Witkowski, B., and 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

Vembar, S. S., Droll, D., and 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

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

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 cynomolgii. *Elife* 6, e29605. doi: 10.7554/eLife.29605.001

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 cynomolgii reporter line reveals *in vitro* malaria hypnozoite reactivation. *Commun. Biol.* 3. doi: 10.1038/s42003-019-0737-3

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

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

Wang, C., Yan, R., Luo, D., Watabe, K., Liao, D.-F., and 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

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

White, M., Amino, R., and 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

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

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.

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

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

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