# A single-cell liver atlas of Plasmodium vivax infection

Liliana Mancio-Silva, Nil Gural, Eliana Real, Marc Wadsworth, Vincent Butty, Sandra March, Niketa Nerurkar, Travis Hughes, Wanlapa Roobsoong, Heather Fleming, et al.

► To cite this version:

Liliana Mancio-Silva, Nil Gural, Eliana Real, Marc Wadsworth, Vincent Butty, et al. A single-cell liver atlas of Plasmodium vivax infection. Cell Host & Microbe, 2022, 30 (7), pp.1048-1060.e5. 10.1016/j.chom.2022.03.034. hal-03700204

HAL Id: hal-03700204

https://hal.science/hal-03700204

Submitted on 1 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cell Host & Microbe

## A single-cell liver atlas of *Plasmodium vivax* infection

### Graphical abstract

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

### Authors

Liliana Mancio-Silva, Nil Gural,
Eliana Real, ...,
Jetsumon Sattabongkot,
Alex K. Shalek, Sangeeta N. Bhatia

### Correspondence

liliana.mancio@inserm.fr (L.M.-S.),
shalek@mit.edu (A.K.S.),
sbhatia@mit.edu (S.N.B.)

### In brief

*Plasmodium vivax* parasites cause relapsing human malaria, which is considered a barrier to its eradication. Mancio-Silva et al. perform time-resolved, dual single-cell transcriptional profiling of this elusive organism to garner insights into hypnozoite dormancy and activation, hepatic sexual commitment, and host immune response to infection.

### Highlights

* MPCC and Seq-Well platforms enable single-cell survey of *P. vivax* liver infection

* Dual transcriptome analysis links dormancy and sexual conversion to host state

* Single-cell profiling reveals dynamic host responses in infected and bystander cells

* Dataset is validated via *in situ* transcript detection and functional assays

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

# A single-cell liver atlas of *Plasmodium vivax* infection

Liliana Mancio-Silva,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="11" />,* Nil Gural,<fnref n="1" />,<fnref n="2" />,<fnref n="11" /> Eliana Real,<fnref n="3" /> Marc H. Wadsworth II,<fnref n="1" />,<fnref n="4" /> Vincent L. Butty,<fnref n="2" />,<fnref n="5" /> Sandra March,<fnref n="1" />,<fnref n="2" />,<fnref n="8" />,<fnref n="9" />

Niketa Nerurkar,<fnref n="1" />,<fnref n="2" /> Travis K. Hughes,<fnref n="1" />,<fnref n="4" /> Wanlapa Roobsoong,<fnref n="6" /> Heather E. Fleming,<fnref n="1" />,<fnref n="2" /> Charlie A. Whittaker,<fnref n="2" />,<fnref n="5" />

Stuart S. Levine,<fnref n="2" />,<fnref n="5" /> Jetsumon Sattabongkot,<fnref n="6" /> Alex K. Shalek,<fnref n="1" />,<fnref n="2" />,<fnref n="4" />,<fnref n="7" />,<fnref n="8" />,* and Sangeeta N. Bhatia<fnref n="1" />,<fnref n="2" />,<fnref n="8" />,<fnref n="9" />,<fnref n="10" />,<fnref n="12" />,*

## SUMMARY

Malaria-causing *Plasmodium vivax* parasites can linger in the human liver for weeks to years and reactivate to cause recurrent blood-stage infection. Although they are an important target for malaria eradication, little is known about the molecular features of replicative and non-replicative intracellular liver-stage parasites and their host cell dependence. Here, we leverage a bioengineered human microliver platform to culture patient-derived *P. vivax* parasites for transcriptional profiling. Coupling enrichment strategies with bulk and single-cell analyses, we capture both parasite and host transcripts in individual hepatocytes throughout the course of infection. We define host- and state-dependent transcriptional signatures and identify unappreciated populations of replicative and non-replicative parasites that share features with sexual transmissive forms. We find that infection suppresses the transcription of key hepatocyte function genes and elicits an anti-parasite innate immune response. Our work provides a foundation for understanding host-parasite interactions and reveals insights into the biology of *P. vivax* dormancy and transmission.

## INTRODUCTION

*Plasmodium* parasite, the causative pathogen of malaria, has a complex life cycle that spans multiple hosts. Disease transmission is initiated upon the bite of an infected *Anopheles* mosquito, which deposits infective parasites, called sporozoites, into the human host. Sporozoites travel to the liver, invade hepatocytes, and replicate, forming thousands of new parasites called merozoites, which eventually break out into the bloodstream, cyclically invading erythrocytes and initiating clinical symptoms. Parasite transmission back to the mosquito is ensured by the development of sexual gametocyte forms during the asexual erythrocytic cycle that are taken up upon bite to restart the life cycle in the insect host. In the case of *P. vivax*, a subset of liver-stage parasites develops into dormant forms called hypnozoites, which can reactivate weeks to years after the initial infection to cause relapsing disease. Thus, the liver stage, which is obligate yet clinically silent and includes relapse-causing hypno-

zoites, presents a unique opportunity for malaria intervention before the onset of symptoms. However, our knowledge of liver-stage malaria and the response of its hepatocyte host to infection is sparse due to difficult access to the parasite and lack of suitable human liver models. To date, much of our historical knowledge has been based on liver biopsies of infected patients, making it challenging to perform mechanistic studies on liver-stage forms, especially the quiescent hypnozoites. Transcriptomic studies hold promise for unveiling mechanistic insights into the liver-stage biology of relapsing *P. vivax*, but the low infection rate and reduced quantity of parasite transcripts in a transcriptionally active host cell environment has made it difficult to perform these studies.

We recently provided some of the first insights into the transcriptional features of liver stages of *P. vivax* by leveraging an *in vitro* primary human liver platform (micropatterned co-cultures [MPCCs]) that recapitulates key aspects of *P. vivax* liver-stage biology, including the establishment of persistent dormant

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

<footnote n="1">¹Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA</footnote>

<footnote n="2">²David H. Koch Institute for Integrative Cancer Research, MIT, Cambridge, MA 02139, USA</footnote>

<footnote n="3">³Institut Pasteur, Université Paris Cité, Inserm U1201, CNRS EMR9195, Unité de Biologie des Interactions Hôte-Parasite, 75015 Paris, France</footnote>

<footnote n="4">⁴Department of Chemistry, MIT, Cambridge, MA 02139, USA</footnote>

<footnote n="5">⁵BioMicro Center, MIT, Cambridge, MA 02139, USA</footnote>

<footnote n="6">⁶Mahidol Vivax Research Unit, Faculty of Tropical Medicine Mahidol University, Bangkok 10400, Thailand</footnote>

<footnote n="7">⁷Ragon Institute of Massachusetts General Hospital, MIT, and Harvard, Cambridge, MA 02139, USA</footnote>

<footnote n="8">⁸Broad Institute of MIT and Harvard, Cambridge, MA 02139, USA</footnote>

<footnote n="9">⁹Howard Hughes Medical Institute, Chevy Chase, MD 20815, USA</footnote>

<footnote n="10">¹⁰The Wyss Institute for Biologically Inspired Engineering Harvard University Boston, MA 02215, USA</footnote>

<footnote n="11">¹¹These authors contributed equally</footnote>

<footnote n="12">¹²Lead contact</footnote>

<footnote>*Correspondence: liliana.mancio@inserm.fr (L.M.-S.), shalek@mit.edu (A.K.S.), sbhatia@mit.edu (S.N.B.)
https://doi.org/10.1016/j.chom.2022.03.034</footnote>
<figure><img src="image_5.png" /></figure>

forms, growing schizonts, merozoite release, and subsequent infection of overlaid erythrocytes. Our previous work revealed reduced transcriptional activity in hypnozoite-enriched samples. Specifically, transcripts encoding functions related to cell division and invasion machinery were suppressed, consistent with a quiescent state (Gural et al., 2018). However, the single-time-point bulk sequencing used in the previous study prevented us from capturing the inherent heterogeneity of the distinct liver-stage parasite forms as well as host-target and bystander cell responses. Achieving a deeper understanding of the pathogen-host interactions has the potential to provide insights into mechanisms that could be leveraged to treat or prevent infection, as suggested by innate interferon (IFN) responses to infection by rodent malaria parasites (Liehl et al., 2014) and a number of viruses (Schneider et al., 2014). However, a closer look into the infection-specific host responses and potential protective responses in uninfected hepatocytes requires single-cell resolution. Recently, diverse single-cell technologies have revealed stage-specific transcriptional signatures in mosquito (Real et al., 2021), blood and gametocyte (Poran et al., 2017) stages from non-relapsing human parasites, and the entire life cycle of rodent parasites (Hewick et al., 2019), all of which can be easily cultured and propagated in laboratories. For *P. vivax*, single-cell transcriptomic studies have been conducted in blood stages collected from infected monkeys (Sà et al., 2020), but the analysis of liver infection has not been performed to date.

Here, we present a comprehensive view of the liver-stage transcriptomes of a human-infecting malaria parasite and its host cells at single-cell resolution. To achieve this goal, we coupled MPCC with Seq-Well, a recently developed low-cost and portable single-cell platform that does not require fluorescent labeling and is compatible with the use of samples collected in endemic settings (Gierahn et al., 2017). With this combined platform, we defined distinct signatures between early-, dormant-, mid-, and late-stage parasites and identified sub-populations of sexually committed forms in the liver, which were previously thought to appear only during erythrocytic infection. We interrogated the pathogen-host interactions and described dynamic innate immune responses in infected and uninfected bystander cells, suggesting a likely mechanism for endogenous protection from secondary infections. We report cytokine- and stage-dependent anti-parasite activity by the induction of IFN signaling and validated host factors involved in cell autonomous innate defense. Together, the data presented here unveil the *P. vivax* liver-stage development with unprecedented resolution and offer insights into the host pathways that might impact their unique biology.

## RESULTS

### Profiling liver-stage *P. vivax* infection by pairing targeted sequencing with single-cell analysis

To comprehensively profile *P. vivax* liver-stage infection, we collected cells from MPCCs of primary human hepatocytes at multiple time points following infection (Figures 1A-1C). Sampling spanned the full liver-stage developmental period (days 1-11), comprising a mix of both replicating schizonts and non-replicative hypnozoites. To obtain hypnozoite-enriched samples, MPCCs were treated from days 5 to 8 with a phospha-

tidylinositol 4-kinase (PI4K) inhibitor compound, a dosing regimen that eliminates the replicative parasites while preserving the dormant parasites (Gural et al., 2018). For hypnozoite-enrichment without drug treatment, we collected cultures on day 14. To obtain a baseline reading of the host transcriptome, naive and mock-exposed MPCCs were prepared and collected in parallel. Mock samples were obtained by exposing MPCCs to non-infected mosquito material. Samples from two independent infections performed with patient-derived *P. vivax* sporozoites were processed for high-throughput single-cell RNA sequencing (scRNA-seq) using Seq-Well S³ (Hughes et al., 2020). Samples were also collected in bulk for RNA sequencing (RNA-seq) analysis.

Because whole transcriptome scRNA-seq resulted in a poor representation of *P. vivax* genes and transcripts, we incorporated an additional step whereby barcoded parasite transcripts were enriched using previously validated nucleic acid baits targeting the entire *P. vivax* genome (Gural et al., 2018). Capture and re-sequencing of samples remarkably increased the efficiency of *P. vivax* gene and transcript detection (Figure S1A). After quality control filtering (Figures S1B and S1C; Tables S1 and S2), we obtained single-cell transcriptomes from 1,494 individual parasites. The daily tally fairly reflects the number of hepatic infections, except for day 11 samples from which a higher number of parasites was recovered (likely representing free merozoites released from mature schizonts during sample processing).

### Single-cell profiling *P. vivax* liver infection defines stage-specific gene signatures

Clustering of parasite transcriptomes with Seurat yielded 8 *P. vivax* clusters, corresponding to distinct developmental liver stages that were visualized by uniform manifold approximation and projection (UMAP; Figure 1D; Table S2). Early-stage parasites from day 1 collection assembled in a single cluster Pv_C1, whereas parasites collected on days 4-8 spread into separate clusters Pv_C2-4. Notably, Pv_C3 and Pv_C4 clusters comprise most of the PI4K-enriched samples, as well as day 14 cultures, and likely represent the hypnozoite parasite population. Late-stage liver schizonts (day 11) scattered across clusters Pv_C5-C8 (Figure 1E). The different samples and developmental stages are concordant with bulk RNA-seq profiling (Figures 1F and S1D; Table S3) and the rodent single-cell malaria atlas, respectively (Figure S1E).

Based on the gene patterns that define each cluster, early liver-stage parasites (Pv_C1) are largely characterized by the expression of conserved genes with unknown function as well as likely residual expression of sporozoite-specific genes (gene ontology [GO] terms related gliding activity and cytoskeleton organization process) (Figures 1G and 1H; Table S4). The Pv_C2 cluster appears to represent a core mid-stage liver program, comprising well-characterized genes, such as those encoding liver-specific protein 1 (LISP1) and 2 (LISP2), and genes implicated in housekeeping functions important for parasite growth, such as translation (EIF3C) and metabolism (ENO and LDH). From Pv_C2 cluster, two branches appear to emerge, one comprised of Pv_C3 and Pv_C4 and the other Pv_C5 and Pv_C6 clusters. PUF1, a Pumilio RNA-binding protein known to be involved in translational repression, and genes encoding proteins with peptidase
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_13.png" /><figcaption>**Figure 1.** Single-cell profiling of *P. vivax* liver-stage infection</figcaption></figure>

(A) Timeline of infection, drug treatment, and sample collection for scRNA-seq (days 0–14). PI4K inhibitor treatment (days 5–8) depletes replicative stages. (B and C) Parasite numbers (B) and size distribution (C). Each dot represents an individual parasite (n = 3–8 wells pooled from 2 independent infections, n = 1 well shown for day 14 PI4K). The traced line indicates the 10-µm cutoff used for hypnozoite staging. The percentage of hypnozoites is as follows: day 5, 62%; day 8, 46%; day 8 PI4K, 83%; day 11, 36%; day 11 PI4K, 90%; day 14, 78%; day 14 PI4K, 67%.

(D and E) UMAP of 1,494 individual *P. vivax* parasites colored by parasite cluster type (D) and sample identity (E) using Seurat. Each dot represents an individual parasite.

(F) Predicted developmental stage based on Spearman correlations with samples collected in bulk and analyzed by RNA-seq.

(G) Heatmap showing top 5 marker genes for each cluster.

(H) UMAP with the expression of PVP01_1262900 (conversed Plasmodium protein, unknown function), LISP2, MSP1, and TRAG31 highlighted.
See also Figure S1 and Tables S1, S2, and S4.

activity (*vivapain-1* and -2) appear in the top 5 marker genes for clusters Pv_C3 and Pv_C4 (Figure 1G). The Pv_C5 and Pv_C6 clusters, on the other hand, are represented by GO

terms related to mitotic division and adhesion of symbiont to host. Gene cluster identifiers in late-stage Pv_C5 and Pv_C6 include several members of merozoite surface proteins
<figure><img src="image_14.png" /><figcaption>**Figure 2. Characterization of P. vivax signature of hepatic sexual commitment**</figcaption></figure>

(A) Average expression within each cell cluster of *P. vivax* orthologs of 591 *P. falciparum* genes previously associated with sexual development, from sexual ring stages to mature gametocytes (Pelle et al., 2015). Parasite clusters indicated on top are temporarily ordered as in Figure 1G.

(B and C) UMAP and relative expression of G377, PVS16, and GEXP5 by quantitative RT-PCR (mean ± SEM; n = 2 independent infections). The 3 genes have been shown to be direct AP2-G targets.

(D) RNA *in situ* hybridization of *P. vivax* parasites at day 10. GEXP5 transcripts shown in magenta. Scale bars, 10 µm. Parasite numbers and size distribution of GEXP5-positive and -negative schizonts (mean ± SEM; n = 4 wells; t test: ****p < 0.0001).

(E) Pie charts showing the proportion of AP2-G positive parasites by cluster. The values are as follows: Pv_C1, 15.4%; Pv_C2, 1.4%; Pv_C3, 16.5%; Pv_C5, 25.5%; Pv_C6, 4.2%; Pv_C8, 0.7%; and Pv_C4 and Pv_C7, 0%.

(F) Average expression of differentially expressed genes between AP2-G positive and AP2-G negative parasites by infection time point. The genes shown are associated with the top GO terms.

(MSP), serine repeat antigen (SERA), and rhoetry-associated membrane antigen (RAMA) and rhoetry neck (RON) multigene families that are necessary for red blood cell invasion (Figure 1G). Finally, among late-stage parasites, Pv_C7 and Pv_C8 represent a previously unidentified population of hepatic parasites that co-express merozoite- and gametocyte-specific genes. Marker genes for these clusters include multiple copies of the tryptophan-rich protein family (TRAG, which are expressed by merozoites of multiple *Plasmodium* species and have been implicated in erythrocyte invasion), the gametocyte antigen G27/25, the gamete release protein (GAMER), and the homolog of gametocyte exported protein 5 (GEXP5, aka PHISTc). Interestingly, several members of the apicomplexan AP2 family of transcription factors (AP2 TFs) were found to be highly expressed all throughout the liver-stage developmental phase. Multiple members appeared specif-

ically assigned to individual clusters, suggesting distinct roles that could contribute to the observed cascade of differential gene expression in the liver (Figures 1G and S1F).

Transcriptional profiling of *P. vivax* liver stages reveals early sexual commitment

To gain further insight into the subpopulation of parasites defined by clusters Pv_C7 and PV_C8, we searched the dataset for additional gametocyte-specific genes obtained from various blood-stage datasets (Obaldia et al., 2018; Pelle et al., 2015; Sà et al., 2020). The expression of canonical sexual markers (PVS16, P28, and P25), including female-specific (RUBV1 and G377), mature and immature gametocyte-specific genes, were detected throughout the developmental liver stages in multiple clusters (Figure 2A). Differential expression patterns were confirmed by quantitative RT-PCR analysis for a small subset
<figure><img src="image_15.png" /></figure>

of gametocyte genes (Figures 2B and 2C). For further validation, we performed in situ hybridization for GEXP5, which appeared highly expressed in day 11 samples (Figures 2B and 2C) and is known as an early sexual stage marker in blood stages of the human malaria parasite *P. falciparum* (Tibúrcio et al., 2015). We found GEXP5 transcripts in 5%–15% of schizonts, with positive signal detected in all individual nuclei of the segmented schizont (Figure 2D).

Moreover, we found that up to 25% of parasites across different clusters express the gene encoding AP2-G, the transcriptional master regulator of sexual development in *Plasmodium* blood stages (Kafsack et al., 2014; Sinha et al., 2014) (Figure 2E). Exploring the dynamics of this TF, we found that the proportion of cells expressing AP2-G peaked on days 1, 5, 11, and 14 (Figure 2F), a wave-like pattern reminiscent of that observed in blood-stage parasites undergoing sexual commitment (Poran et al., 2017). A set of genes encoding other AP2 TFs and the chromatin remodeler (imitation switch [ISWI]) appeared to follow a similar pattern of transcription, suggesting that AP2-G might coordinate with other transcriptional regulators to accomplish sexual commitment in the liver (Figure 2F). On day 11, AP2-G induction is also linked to the genes implicated in egress and invasion, as seen previously for blood stages (Josling et al., 2020).

Our findings corroborate that commitment to gametocytogenesis might occur early during liver-stage development in a subset of parasites, likely leading to the formation of sexually committed merozoites. These results are consistent with previous studies that detected PVS16 protein in a fraction of *P. vivax* hepatic schizonts (Roth et al., 2018; Schafer et al., 2020) and others that found gametocyte-specific transcripts in the blood of volunteers infected with *P. vivax* by mosquito bites (Adapa et al., 2019). In the controlled human infection, gametocyte transcripts were identified in the top-ranking genes as early as day 9 after infection, when the first merozoites emerge from the liver.

## Non-replicative *P. vivax* liver-stages depend on proteolytic activity and are sexually committed

To inform the identification of hypnozoite-specific gene signature(s), we first compared the distribution of non-treated versus the PI4K-treated hypnozoite-containing samples across the UMAP, which revealed that the only cluster significantly enriched for PI4K-treated cells was Pv_C3. This enabled us to calculate a hypnozoite score based on the top markers of that cluster and precisely define replicative and non-replicative hypnozoite states across the UMAP (Figures 3A–3C and S2A). Hypnozoites were significantly enriched in the genes encoding proteins with proteolytic activity (26 genes) and nucleic-acid-binding activity (141 genes), including TFs, DNA and RNA binding, translation, and helicase activity (Figure 3D; Table S4).

Among the TFs, we found the expression of gametocytogenesis regulators AP2-G, AP2-FG (Yuda, 2020), and AP2-O3 (Li et al., 2021), suggesting that a subpopulation of *P. vivax* hypnozoites could be pre-committed to become sexually transmissive forms. The presence of gametocyte-specific transcripts in the hypnozoite-enriched samples was confirmed by quantitative RT-PCR analysis (Figure S2B). By combining the hypnozoite score with AP2-G expression as proxy for sexual commitment, we further identified four distinct parasite populations and states

occurring at varying proportions during liver-stage development (replicative AP2-G positive, hypnozoite AP2-G positive, replicative AP2-G negative, and hypnozoite AP2-G negative, which exhibit varying proportions during liver-stage development; Figure 3E).

To further explore the transcriptional differences among individual hypnozoites, we next re-clustered Pv_C3 and PV_C4, which revealed four diverse subclusters with a distinctive pattern of AP2 TFs expression (Figures 3F and S2C). Interestingly, the Pv_H_C4 subcluster appeared to be depleted of the markers of translational repression PUF1 and DOZI but significantly enriched in a gene encoding a 7-helix protein. This protein has been postulated to operate in the re-initiation of the translation of repressed transcripts during *P. falciparum* gametogenesis (Bennink et al., 2018). Other marker genes for this subcluster are ATG7, which encodes a protein necessary for parasite growth in *P. falciparum* blood stages (Walker et al., 2013), and PHD1, a member of the GCN5 histone acetyltransferase complex that controls the expression of erythrocyte invasion genes (Miao et al., 2021). The Pv_H_C4 subcluster could, therefore, represent a subpopulation of *P. vivax* hypnozoites undergoing reactivation, i.e., shifting from dormancy to schizogony.

Protein degradation is important for the maintenance of a quiescent state induced by nutrient limitation during the *P. falciparum* erythrocytic cycle (Babbitt et al., 2012). To functionally validate our findings and determine whether proteolytic activity (Figure 3G) could be required to maintain the viability of quiescent hypnozoites, we exposed *P. vivax*-infected MPCCs to E64, a membrane-permeable protease inhibitor. No effect was found on schizont forms; however, the number of hypnozoites was significantly reduced in the E64-treated cultures (Figures 3H and S2D). Altogether, these results suggest that *P. vivax* hypnozoites in the liver could represent a reservoir of sexually committed parasites and that dormancy might depend on transcription and translational repression mechanisms and proteolysis for long-lasting viability.

## Dormancy and hepatic sexual differentiation are associated with distinct host metabolic states

To characterize host cell signatures linked to infection and host-dependent differences in parasite transcriptomes, we profiled the samples of unexposed, mock-, and *P. vivax*-exposed MPCCs across all time points (Figure 1A). From the 33,379 high-quality human transcriptomes obtained after quality control filtering (Figure S3A; Table S2), we created a subset encompassing the *P. vivax*-exposed samples across all time points, corresponding to 12,841 cells. Clustering using Seurat revealed 7 human clusters (H_C1-7) with differential enrichment in *P. vivax*-associated hepatocytes (Figures 4A and 4B). Cluster H_C2 harbored only a small number of infected hepatocytes, whereas H_C6 appeared to be highly enriched in infected cells.

Comparing *P. vivax*-positive and -negative cells demonstrated that infected cells are transcriptionally distinct in a time-dependent manner (Figure 4C; Table S5). GO term enrichment analysis revealed strong activation of innate immune defense genes during early infection and suppression of liver-specific functions in *P. vivax*-infected cells in later infection. Examples of differentially expressed genes on day 1 include the markers of inflammation and stress response, such as *LRG1*, *TNFSF10*, and *NFKB1*, and
<figure><img src="image_16.png" /></figure>

A

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

B

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

C

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

D

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

E

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

F

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

G

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

H

<figure><img src="image_24.png" /><figcaption>Figure 3. Characterization of *P. vivax* hypnozoite-stage transcriptomes</figcaption></figure>

(A) Average expression of the 10 genes in the hypnozoite score.

(B and C) UMAP with cells colored by their hypnozoite score (B) and replicative state (C). Cells with a score >0.75 were labeled as hypnozoite.

(D) GO term analysis of differentially expressed genes between hypnozoite and replicative.

(E) UMAP with hypnozoite/replicative states and AP2-G expression highlighted.

(F) UMAP of hypnozoite clusters (re-clustering Pv_C3 and Pv_C4). Gene markers for each hypnozoite subcluster are shown on the bottom.

(G) Violin plot showing the expression of 4 cysteine peptidases upregulated in hypnozoite compared with replicative parasites.

(H) Treatment (days 0–8) of *P. vivax*-infected cultures with protease inhibitor E64 (1 µM; mean ± SEM; n = 3–5 wells; t test: *p* < 0.05).
See also Figure S2.

the family of acute-phase serum amyloid A proteins (SAA-1, -2, and -4) with increased expression in the *P. vivax*-infected cells. LRG1 is a recently identified biomarker that is able to distinguish between malaria and dengue infection (Kumar et al., 2020). A significant number of genes encoding transporters, complement proteins, and genes related to lipid metabolic processes appeared to be induced in the infected population at this time point. Notably, the cholesterol transporter *SRB1* was found expressed in higher levels in 63% of the infected cells on day 1 (Figure S3B), consistent with a recent report identifying it as host factor for *P. vivax* parasites (Manzoni et al., 2017). Other host factors previously identified in other species do not show significant differential expression (Figure S3B; Table S5). Later in infection, however, IFN response and lipid metabolic processes become suppressed. Other significant GO terms appearing dysregulated include apoptosis, translation, and viral replication. These results are

consistent with a model suggesting the shutdown of hepatic function during infection. Albumin (ALB), arginase 1 (ARG1), and fatty-acid-binding protein 1 (FABP1) appeared in the top 10 of most downregulated genes. IFN-related genes with significantly low expression in later infection, such as *ISG15* and *IFI6*, might be involved in creating a permissive environment for late-stage parasite development and replication.

Given that sexual differentiation and gametocytogenesis in malaria parasites have been linked to host-derived physiological signals (Brancucci et al., 2017), we next compared the transcriptomes of the infected hepatocytes based on the sexual and replicative status. To identify hepatocytes hosting replicative and non-replicative parasites, we applied the hypnozoite score defined in Figure 3. For the sexual state, we compared AP2-G-positive and AP2-G-negative *P. vivax* parasites. Interestingly, we found that the hepatocytes bearing non-replicative or
<figure><img src="image_25.png" /><figcaption>B</figcaption></figure>

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

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

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

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

<figure><img src="image_30.png" /><figcaption>Figure 4. Dual scRNA-seq analysis of *P. vivax* liver-stage infection</figcaption></figure>

(A) UMAP of 33,379 individual cells colored by cell-type cluster using Seurat.

(B) UMAP of 12,841 *P. vivax*-exposed hepatocytes showing the distribution of the different samples by day (left) and infection status (right). *P. vivax*-positive hepatocytes colored in dark red.

(C) Volcano plots showing differentially expressed genes in *P. vivax*-positive and -negative cells. Genes colored in green or yellow show a fold change of at least 1.5 (FDR < 0.05). No significantly different genes detected on day 8.

(D and E) Average expression of the top 20 differentially expressed genes (FC > 4 and FDR < 0.05) between infected host cells carrying replicative versus non-replicative parasites (D) and *AP2-G*-positive versus *AP2-G*-negative parasites (E) shown at the indicated infection time points. UMAPs of the different cell types shown on the right.

See also Figure S3 and Table S5.
<figure><img src="image_31.png" /></figure>

A

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

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

B

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

C

<figure><img src="image_35.png" /><figcaption>**Figure 5. IFN responses in P. vivax-infected and bystander hepatocytes**</figcaption></figure>

(A) Dot plots showing the expression of IFN-related genes throughout infection time course in naive, mock-exposed, and P. vivax-exposed samples (top) and in the P. vivax-exposed samples by infection status (bottom).

(B) UMAP highlighting ISG5 and IFI6 expression levels.

(C) Expression of IFITM3 protein (green) in P. vivax-infected cultures. Representative images of parasites on day 8. Scale bars, 50 µm. Bar plot shows frequency distribution of IFITM3-positive hepatocyte layer relative to the infected cell by infection time point for a representative experiment.
See also Figure S4.

AP2-G-positive parasites exhibit significantly lower expression of the genes encoding the proteins implicated in iron transport and storage (TF, FTL, and FTH1) and lipid metabolism (APOA1, APOC1, APOB, FABP1, and CYP3A4). The genes encoding ALB, haptoglobin (HP), and alpha-1-antitrypsin (SERPINA1) proteins are also among the top 20 of differentially expressed (Figures 4D and 4E; Table S5). These observations point to common mechanisms, suggesting that deficiency in specific intracellular metabolites (possibly iron or lipids) could serve as triggers of both quiescence and sexual commitment.

**Interferon responses in uninfected bystander hepatocytes**

Taking advantage of our single-cell approach, we also investigated the impact of P. vivax infection on the bystander hepatocytes.

Comparing the transcriptomes of the cells exposed to P. vivax with those of the mock-exposed or naive hepatocytes, we detected a widespread induction of the alpha/beta IFN response, mainly in the P. vivax-exposed uninfected cells (Figures 5A, 5B, and S4). Day 1 samples stood out given that the induction of multiple IFN-responsive genes and gene families, including IRF7, IFI6, IFI27, IFIT1-5, IFI44, and IFIH1, was observed in more than 40% of the hepatocytes. The transcripts of effector IFN-induced transmembrane proteins IFITM2 and IFITM3 were increased throughout the infection time course (Figures 5A and 5B). Although the strong induction observed at early time points could be linked to the sporozoite cell traversal and wounding of cells prior infection (Mota et al., 2001), the temporal order and gene composition of the IFN response to P. vivax are remarkably similar to those of the IFN response to hepatitis C
<figure><img src="image_36.png" /><figcaption>**Figure 6. IFN responses control *P. vivax* hepatic infection**</figcaption></figure>

(A) IFN treatment (days 5–8) of *P. vivax*-infected MPCCs.

(B) siRNA treatment 3 days prior to infection. SRB1 and CD81 are known host factors for *P. vivax* and *P. falciparum*, respectively, and were used as controls. Bar plots show the quantification of parasite numbers at day 8 (mean ± SEM; n = 5–6 wells pooled from 2 independent infections; one-way ANOVA test: \*p < 0.05; \*\*p < 0.01; \*\*\*p < 0.001; \*\*\*\*p < 0.0001).

See also Figure S5.

virus (HCV) infection (Sheahan et al., 2014), rather suggesting a common hepatocyte defense mechanism in response to human hepatotropic pathogens.

To validate our findings with an independent approach, we performed immunofluorescence analysis of IFITM3 in *P. vivax*-infected MPCCs. We found that IFITM3 protein localizes in the bile canaliculi and tight junctions and confirmed the upregulation of IFITM3 in the uninfected bystander cells, with higher expression detected in the vicinity of the infected hepatocytes (Figure 5C). IFITMs restrict infections by many viruses, employing diverse mechanisms (Diamond and Farzan, 2013). Based on its localization in the host membrane and lack of association with the parasite membrane, IFITM3 is likely to be involved in gap junction communication and propagation of host responses from infected to uninfected cells (Luther et al., 2020; Patel et al., 2009).

## Host interferon responses control *P. vivax* infection

Given the innate immune responses observed during early *P. vivax* infection (Figures 4 and 5), we hypothesized that the activation of this pathway could be associated with the significant drop in the parasite numbers observed from day 1 to 11 (Figure 1B). To assess the impact of activated IFN signaling

on *P. vivax* infection and provide a mechanistic insight into IFN-mediated responses, we employed two approaches: (1) treatment with IFN cytokines and (2) silencing of downstream effector proteins. For the first approach, *P. vivax*-infected MPCCs were treated with IFN alpha (IFNa) and beta (IFNb) cytokines. As control, we used IFN gamma (IFNg), which is known to reduce *P. vivax* infection by activating autophagy (Boonhok et al., 2016; Ferreira et al., 1986). Although no effect was detected in terms of the kinetics of parasite development (Figure S5A), the number of both *P. vivax* schizonts and hypnozoites were significantly reduced upon treatment with IFNa. IFNb appeared to have a stronger effect on schizont stages. By contrast, IFNg treatment revealed a potent anti-hypnozoite activity (Figure 6A).

For the gene silencing approach, we utilized siRNA oligonucleotides to specifically reduce the expression of IFI6 and IFI27 prior to infection. IFI6 and IFI27 belong to the same family of IFNa-induced effector proteins exerting opposing activities on apoptosis. They can exhibit antiviral activity depending on the virus (Qi et al., 2015; Ullah et al., 2021). Both *IFI6* and *IFI27* transcripts were increased upon infection with *P. vivax*; however, only 23% of the *P. vivax*-positive cells co-expressed the 2 genes (Figures 5A, S5B, and S5C). Consistent with an
antagonistic function, the proteins they encode for appeared to have opposing effects on infection. The reduction of IFI6 in MPCCs led to a modest increase in the numbers of *P. vivax* schizonts, whereas the reduction of IFI27 showed a significant decrease in the numbers of schizonts. No effect was found on parasite size or hypnozoite numbers (Figures 6B and S5A). These results support an antiparasitic role for IFI6, likely promoting parasite elimination by the activation of apoptosis in the host cell. IFI27, on the other hand, appears to enhance *P. vivax* infection in human hepatocytes. Our siRNA experiment also validated the postulated role of SRB1 as a species-specific entry host factor for *P. vivax* (Manzoni et al., 2017) (Figure 6B). Collectively, our results suggest that successful *P. vivax* infection likely involves robust mechanisms for subverting host IFNa/b responses, which are distinct in the hepatocytes harboring dormant hypnozoites.

DISCUSSION

Here, we comprehensively surveyed the molecular composition of *P. vivax* parasites in distinct liver stages, as well as host or bystander hepatocytes, throughout the course of infection. The presented data enabled us to assemble a single-cell liver atlas of relapsing human malaria. Leveraging the portability of two established platforms (MPCC and Seq-Well) to culture and collect individual human hepatocytes infected in an endemic location, we performed dual single-cell sequencing analysis and further developed a method to selectively enhance the capture of parasite transcripts. We demonstrated the robustness and utility of this approach by sequencing 1,494 *P. vivax* parasites and 33,383 hepatocyte transcriptomes and provided validations and mechanistic insights for the identified gene signatures and parasite-host interactions.

Dual transcriptional profiling of *P. vivax* infection revealed host- and stage-dependent gene expression patterns in both parasites and hepatocytes. Our data are consistent with the following model: upon invasion, a subset of *P. vivax* sporozoites in cells with reduced hepatic metabolism activate AP2-G. Of those, a fraction will remain dormant, whereas the remaining portion will activate the schizogony program. Thus, at the point of egress, two populations of hepatic *P. vivax* merozoites will emerge from the liver: an asexual population developing into blood stages and a sexually committed population developing into gametocytes for mosquito transmission (Figure S6). This model, although distinct from the existing understanding of the *Plasmodium* life cycle, is in agreement with the historical observations of rapid *P. vivax* transmission that occur before the onset of symptoms (Baker, 2010) and is consistent with the high abundance of gametocyte transcripts detected in the blood collected from malaria patients (Adapa et al., 2019; Kim et al., 2019). In our model, parasites that do not activate schizogony, hypnozoites, remain in a dormant state via transcriptional/translational repression mechanisms and rely on proteolytic activity to sustain viability. Bypassing the need for an asexual replication phase might represent an evolutionary imperative to preserve genome integrity, given the prolonged period between infection and transmission. From a clinical perspective, our results advocate for the use of gametocyte-targeting molecules as a way to reduce the hypnozoite reservoir in the liver. Supporting this strategy, a recent *in vitro*

drug screen reported anti-hypnozoite activity for several previously identified *P. falciparum* gametocidal compounds (Maher et al., 2021). An alternative would be the development of a “commit and kill” clinical approach. Instead of employing epigenetic modulators to first awaken the dormant parasites and then kill the replicating forms with a conventional drug as previously proposed (Dembélé et al., 2014), we would leverage drugs that enhance gametocyte commitment and then use a gametocidal drug to kill. Parasite-specific protease inhibitors could also be screened and investigated for anti-relapsing activity.

Additionally, our work reveals the dysregulation of IFN and inflammatory signaling pathways in *P. vivax*-infected hepatocytes. After a strong initial activation, the reduced expression of IFN-responsive genes in infected cells suggests a protection mechanism for the parasite from host-cell-mediated killing. The downregulation of IFN signaling may also contribute to the under-activation of adaptive immune responses by antigen-presenting cells, leading to parasite persistence in the organism. The mechanisms by which IFN suppression is achieved in infected cells or the impact of host genetics and host-dependent IFN responses, which are key determinants of the clinical outcome in HCV infections (Sheahan et al., 2014), remain unknown. Nevertheless, induction treatments with IFN and gene silencing experiments in our immune-cell-free microliver system began unveiling restricting and enhancing factors of *P. vivax* infection. Furthermore, we described the upregulation of innate immune response pathways in uninfected bystander cells, similar to those observed in viral infections that help control the spread of the virus to neighboring cells (Kotliar et al., 2020; Sheahan et al., 2014). In the context of malaria, this phenomenon might prime or protect cells from a secondary infection, as suggested by the inhibition of malaria re-infection by IFN in a rodent model of *Plasmodium* (Liehl et al., 2014, 2015). Future studies employing spatial single-cell transcriptomics could reveal *P. vivax*-specific gene signatures and spatially heterogeneous responses in bystander cells.

In summary, this study presents a transcriptional description of individual, patient-derived *P. vivax* liver-stages, their host cells, and uninfected bystander cells over the course of infection in human hepatocytes. Leveraging multiple single-cell technologies (host-pathogen sequencing, *in situ* hybridization, and immunofluorescence), we reveal earlier than expected sexual commitment during *P. vivax* liver-stage development in both replicative and non-replicative parasites and a dominant innate immune response that exhibits distinct signatures in infected and uninfected bystander hepatocytes. Taken together, these clinically relevant insights provide a framework for characterizing host-parasite interactions in *P. vivax* infections. We expect the methods described here to be applicable for profiling not only other *Plasmodium* parasites at various stages of the life cycle but also other intracellular pathogens where the low abundance of transcripts or host contamination make it difficult to perform single-cell studies.

STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

● KEY RESOURCES TABLE
## RESOURCE AVAILABILITY

* Lead contact

* Materials availability

* Data and code availability

## EXPERIMENTAL MODEL AND SUBJECT DETAILS

* Fibroblasts and primary cell cultures

* *P. vivax* parasites

## METHODS DETAILS

* MPCCs and *P. vivax* infection

* siRNA and drug treatments

* Seq-Well and *P. vivax* capture

* Tri-genome mapping target generation

* scRNA-seq processing

* scRNA-seq analysis of *P. vivax* transcriptomes

* Bulk RNA-seq analysis

* Comparison between *P. vivax* single-cell and bulk transcriptomes

* Comparison between *P. vivax* and *P. berghei* single-cell transcriptomes

* scRNA-seq analysis of human transcriptomes

* Quantitative RT-PCR

* Immunofluorescence analysis

* Fluorescence in situ hybridization

## QUANTIFICATION AND STATISTICAL ANALYSIS

## SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.chom.2022.03.034.

## ACKNOWLEDGMENTS

We are grateful to Sabrina Hawthorne and Alexei Stortchevoi for help with the targeted scRNA-seq; Sabine Thiberge and Patty Chen for antibody production; and Dan Neafsey, Bronwyn MacInnis, Artur Scherf, and Sandra Pellegrini for insightful discussions. This work was supported by the Bill & Melinda Gates Foundation (OPP1023607, OPP1202327, and INV-027498), the US National Institute of Health (5U19AI110818-07, 5U24AI118672, and P30-CA14051), a Broad Institute BN10 grant, a Sloan Fellowship in Chemistry (to A.K.S.), the French Institut National de la Santé et de la Recherche Médicale (to L.M.-S.), and Agence Nationale de la Recherche (ANR-20-CE15-0031). V.L.B. and C.A.W. were partially supported by Cancer Center Support grant P30-CA14051 from the NCI to the Barbara K. Ostrom (1978) Bioinformatics and Computing Core Facility of the Swanson Biotechnology Center. S.N.B. is a HMMI investigator. We also acknowledge the UtechS Photonic Biolmaging, C2RT, Institut Pasteur, supported by France Biolmaging, ANR-10-INBS-04, and Investments for the Future. The following reagent was obtained through BEI Resources, NIAID, NIH: hybridoma 2F2 anti-*Plasmodium vivax* circumsporozoite protein (CSP), MRA-184, contributed by Elizabeth Nardin.

## AUTHOR CONTRIBUTIONS

Conceptualization, L.M.-S., N.G., M.H.W., T.K.H., S.M., S.S.L., A.K.S., and S.N.B.; investigation, L.M.-S., N.G., M.H.W., T.K.H., and N.N.; resources, W.R. and J.S.; software and formal analysis, E.R., M.H.W., V.L.B., and C.A.W.; writing – original draft, L.M.-S. and N.G.; writing – review & editing, L.M.-S., N.G., E.R., H.E.F., A.K.S., and S.N.B.; supervision and funding acquisition, L.M.-S., A.K.S., and S.N.B.

## DECLARATION OF INTERESTS

A.K.S. reports compensation for consulting and/or SAB membership from Merck, Honeycomb Biotechnologies, Cellarity, Repertoire Immune Medicines, Hovione, Ochre Bio, Third Rock Ventures, Relation Therapeutics, Empress

Therapeutics, FL82, and Dahlia Biosciences. S.N.B. reports compensation for cofounding, consulting, and/or board membership in Glympse Bio, Satellite Bio, CEND Therapeutics, Xilio Therapeutics, Catalio Capital, Intergalactic Therapeutics, Port Therapeutics, Vertex Pharmaceuticals, Danaher, and Moderna and receives sponsored research funding from Johnson & Johnson, Revitope, and Owstone.

Received: July 30, 2021

Revised: January 31, 2022

Accepted: March 25, 2022

Published: April 19, 2022

## REFERENCES

Adapa, S.R., Taylor, R.A., Wang, C., Thomson-Luque, R., Johnson, L.R., and Jiang, R.H.Y. (2019). Plasmodium vivax readiness to transmit: implication for malaria eradication. BMC Syst. Biol. 13, 5.

Ashburner, M., Ball, C.A., Blake, J.A., Botstein, D., Butler, H., Cherry, J.M., Davis, A.P., Dolinski, K., et al. (2011). Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29.

Babbitt, S.E., Altenhofen, L., Cobbold, S.A., Istvan, E.S., Fennell, C., Doerig, C., Llinás, M., and Goldberg, D.E. (2012). *Plasmodium falciparum* responds to amino acid starvation by entering into a hibernatory state. Proc. Natl. Acad. Sci. USA 109, E3278–E3287.

Baker, D.A. (2010). Malaria gametocytogenesis. Molecular and Biochemical Parasitology 172, 57–65.

Bennink, S., Bohl, A. Von, Ngwa, C.J., Henschel, L., Kuehn, A., Pilch, N., Weißbach, T., Rosinski, A.N., Scheuermayer, M., Repnik, U., et al. (2018). A seven-helix protein constitutes stress granules crucial for regulating translation during human-to-mosquito transmission of *Plasmodium falciparum*. PLoS Pathog. 14, e1007249.

Boonhok, R., Rachaphaew, N., Duangmanee, A., Chobson, P., Pattaradilokrat, S., Utaisincharoen, P., Sattabongkot, J., and Ponpuak, M. (2016). LAP-like process as an immune mechanism downstream of IFN-γ in control of the human malaria *Plasmodium vivax* liver stage. Proc. Natl. Acad. Sci. USA 113, E3519–E3528.

Brancucci, N.M.B., Gerdt, J.P., Wang, C., De Niz, M., Philip, N., Adapa, S.R., Zhang, M., Hitz, E., Niederwieser, I., Boltryk, S.D., et al. (2017). Lysophosphatidylcholine regulates sexual stage differentiation in the human malaria parasite *Plasmodium falciparum*. Cell 171, 1532–1544.e15.

Chen, A.X., Chhabra, A., Song, H.G., Fleming, H.E., Chen, C.S., and Bhatia, S.N. (2020). Controlled apoptosis of stromal cells to engineer human micro-livers. Adv. Funct. Mater. 30, 1–10.

Dembélé, L., Franetich, J.-F., Lorthiois, A., Gego, A., Zeeman, A.-M., Kocken, C.H.M., Le Grand, R., Dereuddre-Bosquet, N., van Gemert, G.-J., Sauerwein, R., et al. (2014). Persistence and activation of malaria hypnozoites in long-term primary hepatocyte cultures. Nat. Med. 20, 307–312.

Diamond, M.S., and Farzan, M. (2013). The broad-spectrum antiviral functions of IFIT and IFITM proteins. Nat. Rev. Immunol. 13, 46–57.

Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., and Gingeras, T.R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21.

Ferreira, A., Schofield, L., Enea, V., Schellekens, H., Meide, P. van der, Collins, W.E., Nussenzweig, R.S., and Nussenzweig, V. (1986). Inhibition of development of exoerythrocytic forms of malaria parasites by gamma-interferon. Science 232, 881–884.

Gene Ontology Consortium (2021). The Gene Ontology resource: enriching a GOLD mine. Nucleic Acids Res. 49, D325–D334.

Gierahn, T.M., Wadsworth, M.H., Hughes, T.K., Bryson, B.D., Butler, A., Satija, R., Fortune, S., Love, J.C., and Shalek, A.K. (2017). Seq-Well: portable, low-cost RNA sequencing of single cells at high throughput. Nat. Methods 14, 395–398.

Gural, N., Mancio-Silva, L., Miller, A.B., Galstian, A., Butty, V.L., Levine, S.S., Patrapuvich, R., Desai, S.P., Mikolajczak, S.A., Kappe, S.H.I., et al. (2018).
<figure><img src="image_37.png" /></figure>

In vitro culture, drug sensitivity, and transcriptome of Plasmodium vivax hypnozoites. Cell Host Microbe 23, 395-406.e4.

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W.M., Zheng, S., Butler, A., Lee, M.J., Wilk, A.J., Darby, C., Zager, M., et al. (2021). Integrated analysis of multimodal single-cell data. Cell 184, 3573-3587.e29.

Howick, V.M., Russell, A.J.C., Andrews, T., Heaton, H., Reid, A.J., Natarajan, K., Butungi, H., Metcalf, T., Verzier, L.H., Rayner, J.C., et al. (2019). The malaria cell atlas: single parasite transcriptomes across the complete Plasmodium life cycle. Science 365, eaaw2619.

Hughes, T.K., Wadsworth, M.H., Gierahn, T.M., Do, T., Weiss, D., Andrade, P.R., Ma, F., de Andrade Silva, B.J., Shao, S., Tsoi, L.C., et al. (2020). Second-strand synthesis-based massively parallel scRNA-seq reveals cellular states and molecular features of human inflammatory skin pathologies. Immunity 53, 878-894.e7.

Josling, G.A., Russell, T.J., Venezia, J., Orchard, L., van Biljon, R., Painter, H.J., and Llinás, M. (2020). Dissecting the role of PfAP2-G in malaria gametocytogenesis. Nat. Commun. 11, 1503.

Kafsack, B.F.C., Rovira-Graells, N., Clark, T.G., Bancells, C., Crowley, V.M., Campino, S.G., Williams, A.E., Drought, L.G., Kwiatkowski, D.P., Baker, D.A., et al. (2014). A transcriptional switch underlies commitment to sexual development in malaria parasites. Nature 507, 248-252.

Kim, A., Popovici, J., Menard, D., and Serre, D. (2019). Plasmodium vivax transcriptomes reveal stage-specific chloroquine response and differential regulation of male and female gametocytes. Nat. Commun. 10, 25-28.

Kiselev, V.Y., Yiu, A., and Hemberg, M. (2018). Smap: projection of single-cell RNA-seq data across data sets. Nat. Methods 15, 359-362.

Kotliar, D., Lin, A., Logue, J., Hughes, T., Khoury, N., Raju, S., Wadsworth, M., Chen, H., Kurtz, J., Dighero-Kemp, B., et al. (2020). Single-cell profiling of Ebola virus infection in vivo reveals viral and host transcriptional dynamics. Cell 183, 1383-1401.e19.

Kumar, V., Ray, S., Aggarwal, S., Biswas, D., Jadhav, M., Yadav, R., Sabnis, S.V., Banerjee, S., Talukdar, A., Kochar, S.K., et al. (2020). Multiplexed quantitative proteomics provides. Commun. Biol. 3, 683.

Li, B., and Dewey, C.N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323.

Li, Z., Cui, H., Guan, J., Liu, C., Yang, Z., and Yuan, J. (2021). Plasmodium transcription repressor AP2-O3 regulates sex-specific identity of gene expression in female gametocytes. EMBO Rep. 22, e51660.

Liehl, P., Meireles, P., Albuquerque, I.S., Pinkevych, M., Baptista, F., Mota, M.M., Davenport, M.P., and Prudêncio, M. (2015). Innate immunity induced by Plasmodium liver infection inhibits malaria reinfections. Infect. Immun. 83, 1172-1180.

Liehl, P., Zuzarte-Luís, V., Chan, J., Zillinger, T., Baptista, F., Carapau, D., Konert, M., Hanson, K.K., Carret, C., Lassnig, C., et al. (2014). Host-cell sensors for Plasmodium activate innate immunity against liver-stage infection. Nat. Med. 20, 47-53.

Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550.

Luther, J., Khan, S., Gala, M.K., Kedrin, D., Sridharan, G., Goodman, R.P., Garber, J.J., Masia, R., Diagacomo, E., Adams, D., et al. (2020). Correction for Luther et al., Hepatic gap junctions amplify alcohol liver injury by propagating cGAS-mediated IRF3 activation. Proc. Natl. Acad. Sci. USA 117, 16704.

Macosko, E.Z., Basu, A., Satija, R., Nemesh, J., Shekhar, K., Goldman, M., Tirosh, I., Bialas, A.R., Kamitaki, N., Martersteck, E.M., et al. (2015). Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell 161, 1202-1214.

Maher, S.P., Vantaux, A., Chaumeau, V., Chua, A.C.Y., Cooper, C.A., Andolina, C., Péneau, J., Rouillier, M., Rizopoulos, Z., Phal, S., et al. (2021). Probing the distinct chemosensitivity of Plasmodium vivax liver stage parasites and demonstration of 8-aminoquinoline radical cure activity in vitro. Sci. Rep. 11, 19905.

Mancio-Silva, L., Fleming, H.E., Miller, A.B., Milstein, S., Liebow, A., Haslett, P., Sepp-Lorenzino, L., and Bhatia, S.N. (2019). Improving drug discovery by

nucleic acid delivery in engineered human microlivers. Cell Metab. 29, 727-735.e3.

Manzoni, G., Marinach, C., Topçu, S., Briquet, S., Grand, M., Tolle, M., Gransagne, M., Lescar, J., Andolina, C., Franetich, J.F., et al. (2017). Plasmodium P36 determines host cell receptor usage during sporozoite invasion. eLife 6, e25903.

March, S., Ramanan, V., Trehan, K., Ng, S., Galstian, A., Gural, N., Scull, M.A., Shlomai, A., Mota, M.M., Fleming, H.E., et al. (2015). Micropatterned coculture of primary human hepatocytes and supportive cells for the study of hepatotropic pathogens. Nat. Protoc. 10, 2027-2053.

McCarthy, D.J., Campbell, K.R., Lun, A.T.L., and Wills, Q.F. (2017). Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics 33, 1179-1186.

Miao, J., Wang, C., Lucky, A.B., Liang, X., Min, H., Adapa, S.R., Jiang, R., Kim, K., and Cui, L. (2021). A unique GCN5 histone acetyltransferase complex controls erythrocyte invasion and virulence in the malaria parasite Plasmodium falciparum. PLoS Pathog. 17, e1009351.

Mota, M.M., Pradel, G., Vanderberg, J.P., Hafalla, J.C.R., Frevert, U., Nussenzweig, R.S., Nussenzweig, V., and Rodríguez, A. (2001). Migration of Plasmodium sporozoites through cells before infection. Science 291, 141-144.

Obaldia, N., Meibalan, E., Sa, J.M., Ma, S., Clark, M.A., Mejia, P., Moraes Barros, R.R., Otero, W., Ferreira, M.U., Mitchell, J.R., et al. (2018). Bone marrow is a major parasite reservoir in Plasmodium vivax infection. mBio 9, 1-16.

Patel, S.J., King, K.R., Casali, M., and Yarmush, M.L. (2009). DNA-triggered innate immune responses are propagated by gap junction communication. Proc. Natl. Acad. Sci. USA 106, 12867-12872.

Pelle, K.G., Oh, K., Buchholz, K., Narasimhan, V., Joice, R., Milner, D.A., Brancucci, N.M.B., Ma, S., Voss, T.S., Ketman, K., et al. (2015). Transcriptional profiling defines dynamics of parasite tissue sequestration during malaria infection. Genome Med. 7, 19.

Poran, A., Nötzel, C., Aly, O., Mencia-Trinchant, N., Harris, C.T., Guzman, M.L., Hassane, D.C., Elemento, O., and Kafsack, B.F.C. (2017). Single-cell RNA sequencing reveals a signature of sexual commitment in malaria parasites. Nature 551, 95-99.

Qi, Y., Li, Y., Zhang, Y., Zhang, L., Wang, Z., Zhang, X., Gui, L., and Huang, J. (2015). IFI6 inhibits apoptosis via mitochondrial-dependent pathway in dengue virus 2 infected vascular endothelial cells. PLoS One 10, e0132743.

Real, E., Howick, V.M., Dahalan, F.A., Witmer, K., Cudini, J., Andradi-Brown, C., Blight, J., Davidson, M.S., Dogga, S.K., Reid, A.J., et al. (2021). A single-cell atlas of Plasmodium falciparum transmission through the mosquito. Nat. Commun. 12, 3196.

Roth, A., Maher, S.P., Conway, A.J., Ubalee, R., Chaumeau, V., Andolina, C., Kaba, S.A., Vantaux, A., Bakowski, M.A., Thomson-Luque, R.T., et al. (2018). A comprehensive model for assessment of liver stage therapies targeting Plasmodium vivax and Plasmodium falciparum. Nat. Commun. 9, 1837.

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.

Schäfer, C., Kappe, H.I., Roobsoong, W., Kangwanrangsan, N., Bardelli, M., Rawlinson, T.A., Dambrauskas, N., Trakhimets, O., Parthiban, C., Goswami, D., et al. (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.

Schneider, C.A., Rasband, W.S., and Eliceiri, K.W. (2012). NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671-675.

Schneider, W.M., Chevillotte, M.D., and Rice, C.M. (2014). Interferon-stimulated genes: a complex web of host defenses. Annu. Rev. Immunol. 32, 513-545.

Sheahan, T., Imanaka, N., Marukian, S., Dorner, M., Liu, P., Ploss, A., and Rice, C.M. (2014). Interferon lambda alleles predict innate antiviral immune responses and hepatitis C virus permissiveness. Cell Host Microbe 15, 190-202.
Sinha, A., Hughes, K.R., Modrzynska, K.K., Otto, T.D., Pfander, C., Dickens, N.J., Religa, A.A., Bushell, E., Graham, A.L., Cameron, R., et al. (2014). A cascade of DNA-binding proteins for sexual commitment and development in Plasmodium. Nature 507, 253–257.

Tibúrcio, M., Dixon, M.W.A., Looker, O., Younis, S.Y., Tilley, L., and Alano, P. (2015). Specific expression and export of the Plasmodium falciparum Gametocyte Exported Protein-5 marks the gametocyte ring stage. Malar. J. 14, 334.

Ullah, H., Sajid, M., Yan, K., Feng, J., He, M., Shereen, M.A., Li, Q., Xu, T., Hao, R., Guo, D., et al. (2021). Antiviral activity of interferon alpha-inducible protein

27 against hepatitis B virus gene expression and replication. Front. Microbiol. 12, 656353.

Walker, D.M., Mahfooz, N., Kemme, K.A., Patel, V.C., Spangler, M., and Drew, M.E. (2013). Plasmodium falciparum erythrocytic stage parasites require the putative autophagy protein PfAtg7 for normal growth. PLoS One 8, e67047.

Yuda, M., Kaneko, I., Iwanaga, S., Murata, Y., and Kato, T. (2020). Female-specific gene regulation in malaria parasites by an AP2-family transcription factor. Mol. Microbiol. 113, 40–51.
## STAR★METHODS

### KEY RESOURCES TABLE

<table><thead><tr><th>REAGENT or RESOURCE</th><th>SOURCE</th><th>IDENTIFIER</th></tr></thead><tbody><tr><td colspan="3"><strong>Antibodies</strong></td></tr><tr><td>Alexa 635-conjugated anti-mouse</td><td>Thermo Fisher Scientific</td><td>Cat# A-31575, RRID:AB_2536185</td></tr><tr><td>Alexa 546-conjugated anti-mouse</td><td>Thermo Fisher Scientific</td><td>Cat# A-11030, RRID:AB_144695</td></tr><tr><td>Alexa 647-conjugated anti-rabbit</td><td>Thermo Fisher Scientific</td><td>Cat# A-21246, RRID:AB_1500778</td></tr><tr><td>Alexa 488-conjugated anti-rabbit</td><td>Thermo Fisher Scientific</td><td>Cat# A-11008, RRID:AB_143051</td></tr><tr><td>IFITM3 rabbit monoclonal</td><td>Cell Signaling Technology</td><td>Cat# 59212 RRID:AB_2799561</td></tr><tr><td>PvCSP mouse monoclonal (clone 2F2)</td><td>BEI Resources</td><td>MRA-184</td></tr><tr><td>PvBip rabbit polyclonal</td><td>This manuscript</td><td>N/A</td></tr><tr><td>PvUIS4 mouse polyclonal</td><td>Gift of Stefan Kappe (Seattle Children's Research Institute)</td><td>N/A</td></tr><tr><td colspan="3"><strong>Biological samples</strong></td></tr><tr><td>Primary human hepatocytes</td><td>BioIVT</td><td>Cat# M00995-P</td></tr><tr><td colspan="3"><strong>Chemicals, peptides, and recombinant proteins</strong></td></tr><tr><td>MMV390048</td><td>Malaria Medicines Venture</td><td>N/A</td></tr><tr><td>IFN alpha 2a</td><td>R&amp;D Systems</td><td>11100-1</td></tr><tr><td>IFN beta</td><td>R&amp;D Systems</td><td>8499-IF-010</td></tr><tr><td>IFN gamma</td><td>R&amp;D Systems</td><td>285-IF-100</td></tr><tr><td>E-64</td><td>Sigma</td><td>E3132</td></tr><tr><td>B/B Homodimerizer</td><td>Clontech</td><td>635059</td></tr><tr><td colspan="3"><strong>Deposited data</strong></td></tr><tr><td>Raw sequencing data</td><td>This manuscript</td><td>GSE197459</td></tr><tr><td>Count matrices of the quality-controlled data</td><td>This manuscript</td><td>zenodo.6280956</td></tr><tr><td colspan="3"><strong>Experimental models: Cell lines</strong></td></tr><tr><td>J2-3T3 murine embryonic fibroblasts</td><td>gift of Howard Green (Harvard Medical School)</td><td>N/A</td></tr><tr><td colspan="3"><strong>Experimental models: Organisms/strains</strong></td></tr><tr><td>Plasmodium vivax</td><td>Mahidol Vivax Research Unit (Bangkok, Thailand)</td><td>N/A</td></tr><tr><td colspan="3"><strong>Oligonucleotides</strong></td></tr><tr><td>ON-TARGETplus CD81 siRNA</td><td>Horizon</td><td>Cat# L-017257-00-0005</td></tr><tr><td>ON-TARGETplus SRB1 siRNA</td><td>Horizon</td><td>Cat# L-010592-00-0005</td></tr><tr><td>ON-TARGETplus IFI6 siRNA</td><td>Horizon</td><td>Cat# L-003672-00-0005</td></tr><tr><td>ON-TARGETplus IFI27 siRNA</td><td>Horizon</td><td>Cat# L-006465-00-0005</td></tr><tr><td>ON-TARGETplus Non-Targeting Control Pool</td><td>Horizon</td><td>Cat# D-001810-10-05</td></tr><tr><td>Primer sequences in Table S6</td><td>Integrated DNA Technologies</td><td>N/A</td></tr><tr><td colspan="3"><strong>Software and algorithms</strong></td></tr><tr><td>Image Processing</td><td>https://imagej.net/software/fiji/ (Schneider et al., 2012)</td><td>N/A</td></tr><tr><td>NIS-Elements Advanced Research</td><td>www.microscope.healthcare.nikon.com/ products/software</td><td>N/A</td></tr><tr><td>Zeiss Zen</td><td>https://www.zeiss.com/microscopy/int/ products/microscope-software/zen.html</td><td>N/A</td></tr><tr><td>GraphPad Prism 9</td><td>www.graphpad.com/scientific-software/ prism</td><td>N/A</td></tr></tbody></table>

(Continued on next page)
Continued

<table><thead><tr><td>REAGENT or RESOURCE</td><td>SOURCE</td><td>IDENTIFIER</td></tr></thead><tbody><tr><td>RStudio (v 1.4.1717)</td><td>Rstudio<br>https://www.rstudio.com</td><td>N/A</td></tr><tr><td>DropSeq (v 1.12)</td><td>(Gierahn et al., 2017)<br>https://github.com/broadinstitute/Drop-seq</td><td>N/A</td></tr><tr><td>STAR (v 2.5.3a)</td><td>(Dobin et al., 2013)<br>https://github.com/alexdobin/STAR</td><td>N/A</td></tr><tr><td>Picard 533 v. 2.9.0-1-gf5b9f50-SNAPSHOT</td><td>Broad Institute<br>https://github.com/broadinstitute/picard</td><td>N/A</td></tr><tr><td>scater (v 1.20.1)</td><td>(McCarthy et al., 2017)<br>https://github.com/jimhester/scater</td><td>N/A</td></tr><tr><td>Seurat (v 4.0.5)</td><td>(Hao et al., 2021)<br>https://github.com/satijalab/seurat</td><td>N/A</td></tr><tr><td>Scmap (v 1.14.0)</td><td>(Kiselev et al., 2018)<br>https://github.com/hemberg-lab/scmap</td><td>N/A</td></tr><tr><td>Hmisc (v 4.6-0).</td><td>https://cran.r-project.org/web/packages/Hmisc/index.html</td><td>N/A</td></tr><tr><td>pheatmap (v 1.0.12)</td><td>https://cran.r-project.org/web/packages/pheatmap/pheatmap.pdf</td><td>N/A</td></tr><tr><td>RSEM (v 1.3.0)</td><td>(Li and Dewey, 2011)<br>https://github.com/deweylab/RSEM/releases/tag/v1.3.0</td><td>N/A</td></tr><tr><td>DESeq2</td><td>(Love et al., 2014)<br>https://bioconductor.org/packages/release/bioc/html/DESeq2.html</td><td>N/A</td></tr></tbody></table>

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Sangeeta N. Bhatia (sbhatia@mit.edu).

Materials availability

Antibodies generated in this study will be shared by the lead contact upon request.

Data and code availability

Bulk RNA-seq and scRNA-seq data have been deposited at GEO (GSE197459). Gene-by-cell expression matrices have been deposited at Zenodo (doi:zenodo.6280956). Any additional information required to reanalyze the data reported in this paper will be available from the lead contact upon request. Microscopy data reported in this paper will be shared by the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Fibroblasts and primary cell cultures

J2-3T3 male murine embryonic fibroblasts (gift of Howard Green, Harvard Medical School) were cultured at < 20 passages in medium comprising of Dulbecco's Modified Eagle Medium (DMEM, Corning), 10% (v/v) bovine serum (Thermo Fisher Scientific), and 100 mg/ mL penicillin-streptomycin (Corning) and were kept at 37°C in a 5% CO2 environment.

Cryopreserved primary human hepatocytes were purchased from BioVT, a vendor permitted to sell products derived from human organs procured in the United States of America by federally designated Organ Procurement Organizations. Human hepatocytes (male donor, age 57) were maintained in DMEM with 10% fetal bovine serum (FBS, GIBCO), 1% ITS (insulin/transferrin/selenous acid and linoleic acid, BD Biosciences), 7 ng/mL glucagon (Sigma-Aldrich), 40 ng/mL dexamethasone (Sigma-Aldrich), 15 mM HEPES (GIBCO), and 100 mg/mL penicillin-streptomycin (Corning). Hepatocyte cultures were kept at 37°C in a 5% CO2 environment.

P. vivax parasites

Anopheles dirus mosquitoes were fed on blood collected from symptomatic patients attending malaria clinics in Tak, Songkla, and Ubon-Ratchathani Provinces in Thailand, confirmed positive for only *P. vivax* via microscopy and RT-PCR. The *P. vivax* infected blood
samples were collected from patients under the approved protocol by the Ethics Committee of the Faculty of Tropical Medicine, Mahidol University (MUTM 2018-016). The written informed consent was obtained prior to blood collection.

Briefly, *P. vivax* infected blood was drawn into heparinized tubes and kept at 37°C until processing. Infected blood was washed once with RPMI 1640 incomplete medium. Packed infected blood was resuspended in warm non-heat inactivated naïve human AB serum for a final hematocrit of 50%. Resuspended blood was fed to laboratory reared female *Anopheles dirus* mosquitoes for 30 minutes via an artificial membrane attached to a water-jacketed glass feeder kept at 37°C. Engorged mosquitoes were kept on 10% sugar at 26°C under 80% humidity at the designated insectary at the Mahidol Vivax Research Unit. Sporozoites were dissected from the salivary glands of infected mosquitoes 14-21 days after blood feeding and pooled in DMEM supplemented with 200 mg/mL penicillin-streptomycin.

## METHODS DETAILS

### MPCCs and *P. vivax* infection

Primary human hepatocytes were seeded on collagen-micropatterned 96-well plates and surrounded with murine embryonic fibroblasts 3T3-J2s as detailed previously (March et al., 2015). For the scRNA-seq analysis, MPCCs were established using 3T3-J2s expressing an inducible apoptosis switch (Chen et al., 2020). MPCCs were infected with fresh sporozoites obtained through hand dissection of *P. vivax*-infected mosquitoes. For mock samples, MPCCs were exposed to material from non-infected mosquito salivary glands (matched number of dissected mosquitoes and processed similarly as the infected counterparts). Naïve samples received the culture media used for the infections.

### siRNA and drug treatments

To obtain hypnozoite-enriched samples, *P. vivax*-infected MPCCs were dosed with a schizont-specific drug Pi4K inhibitor (MMV390048, 1 µM) for 3 days starting on day 5 after infection. Similar dosing schedule was used for treatments with IFN alpha 2a, beta and gamma. E64 treatment started at 4h post-infection after media washing the cultures.

siRNA oligonucleotides were added to patterned hepatocytes as described previously (Mancio-Silva et al., 2019). ON-TARGET-plus SMARTpool siRNA oligonucleotides (Dharmacon) were delivered to hepatocytes by using RNAiMAX Transfection Reagent (Thermo Fisher Scientific) per manufacture's protocols at final concentration 50 nM. *P. vivax* sporozoites were added at day 3 post siRNA treatment.

### Seq-Well and *P. vivax* capture

For the scRNA-seq analysis, an updated Seq-Well protocol including a second-strand synthesis step was employed (Hughes et al., 2020). Briefly, MPCCs were first treated with B/B Homodimerizer (0.5 µM) for 30 minutes at 37°C to partially deplete the fibroblast cells via apoptosis (Chen et al., 2020). After washing, the cultures were dissociated by Trypsin (0.25%) 5-minute treatment at 37°C. A suspension with 10,000 - 15,000 cells (pooled from duplicate MPCC wells) was then loaded onto a functionalized-polydimethylsiloxane array preloaded with uniquely barcoded mRNA capture beads. For each time-point in the two independent infections, two separate arrays were loaded with *P. vivax*-infected samples, and single arrays for mock and naïve samples. After cells had settled into wells, the array was sealed with a hydroxylated polycarbonate membrane with a pore size of 10 nm, facilitating buffer exchange while permitting cell lysis, mRNA transcript hybridization to beads, and reverse transcription. The obtained bead-bound cDNA product then underwent Exonuclease I treatment to remove excess primer before proceeding with second-strand synthesis and PCR amplification. Illumina sequencing (NOVAseq6000) was performed for one of the two infection batches.

To capture parasite reads from Seq-Well, full length cDNAs were amplified with an additional 5 cycles using Kapa HiFi polymerase (3-minute extension time was used to increase concentration for capture). 200 - 300 ng of cDNA was concentrated using a speedvac, reconstituted in 3.4 µL water and hybridized for 32 hours as previously described (Gural et al., 2018). 10 µL of capture material was amplified for 15 cycles following the standard protocol. Captured cDNA was then prepared into Illumina libraries using NexteraXT (Illumina). Final libraries were quality controlled using Fragment Analyzer (Agilent) and qPCR prior to Illumina sequencing (NextSeq500). Sequencing of captured parasite transcripts was performed for the two infection batches.

### Tri-genome mapping target generation

Chromosomes and contigs from human (hg19), murine (mm10) [same ENSEMBL releases as in (Macosko et al., 2015)] and *P. vivax* genome (PvP01 v1 release) were renamed with a species-specific prefix and concatenated into a fasta file. Corresponding gtf files were adapted to match prefixed chromosome names, and genes names were adjusted as an ENSEMBLID_SPECIES_GeneSymbol concatenated string.

### scRNA-seq processing

For each sample, fastq files originating from multiple sequencing runs were concatenated and processed using an analytical pipeline derived from the DropSeq pipeline v. 1.12, as described in (Gierahn et al., 2017) (https://github.com/broadinstitute/Drop-seq). Briefly, reads were converted to a bam file using picard v. 2.9.0-1-gf5b9f50-SNAPSHOT, tagged with cell and transcript barcodes, and subsequently sequencing adapters and polyadenosine tracts were trimmed. Upon regenerating fastq files, reads were aligned with STAR v. 2.5.3a (Dobin et al., 2013) against the aforementioned tri-genome reference. Genomic features of the aligned reads were
annotated using the combined species-gene symbol nomenclature including gene and exon of origin when relevant. Bead synthesis errors were assessed and when possible, altered unique molecular identifiers (UMIs) were repaired. Cell barcode abundance was tallied, and gene expression was called for the top 10,000 cell barcodes. Count matrices of genes x cells were imported in the R v. 3.6 statistical environment and Seurat (v 4.0.5) was used as the primary analytical package (Hao et al., 2021). Count matrices were merged into a single Seurat object, which was split by the sample of origin. Human, murine and *P. vivax* genes and transcripts counts were tallied for each cell barcode.

### scRNA-seq analysis of *P. vivax* transcriptomes

The Seq-Well analytical pipeline was run on the samples as described above. For samples for which pre and post-capture data were available (batch 2), matching cell barcodes were identified in the corresponding genes x cell matrices and UMI counts were summed over pre and post-capture libraries. The resulting matrices were merged into a Seurat (v 4.0.5) object (2,171 cells, 85,276 genes). Only *P. vivax* genes were retained for subsequent analysis. Cell outliers were initially removed from the dataset by running the runColData-taPCA() function in scater (v 1.20.1) (McCarthy et al., 2017). Samples with > 5% of mitochondrial genes and > 30% of total reads mapping to rRNA genes were also excluded from the analysis. Finally, the Seurat object was split by batch and each individual object was further filtered by excluding cells with fewer than 20 genes and 50 reads for batch 1 and with fewer than 100 genes and 200 reads for batch 2. The two objects were then normalized with SCTransform and integrated using the following functions in Seurat: Select-tlntegrationFeatures(), PrepSCTintegration(), FindIntegrationAnchors(), RunPCA(), IntegrateData(). Transcriptomes were visualized using UMAP dimensionality reduction in Seurat with n_neighbors = 20, min.dist = 0.3, repulsion.strength = 0.02, and localconnec-tivity = 3. For the identification of cluster-specific gene markers and differential gene expression between AP2G-positive and AP2G-negative cells, the FindAllMarkers() and FindMarkers() functions in Seurat were used with a MAST test.

### Bulk RNA-seq analysis

fastq files were mapped using STAR v. 2.5.3a (Dobin et al., 2013) against PVP01 *P. vivax* v1 genome assembly and annotation, and quantitated by RSEM v. 1.3.0 (Li and Dewey, 2011). Differential expression analysis was performed in DESeq2 on protein-coding genes (Love et al., 2014).

### Comparison between *P. vivax* single-cell and bulk transcriptomes

*P. vivax* single cells were assigned one of several possible developmental stages based on gene expression correlations between scRNA-seq and bulk transcriptomes. Pearson and Spearman correlations were computed using the Hmisc package in R (v 4.6-0).

### Comparison between *P. vivax* and *P. berghei* single-cell transcriptomes

*P. vivax* cells were compared with a subset of the mouse Malaria Cell Atlas dataset (liver and blood asexual and sexual stages; https://zenodo.org/record/2843883#.Ygz7wC-I1bU) using the sRAINmap() function in sRAIN (v 1.14.0) (Kiselev et al., 2018) with the MCA data as the reference index and the *P. vivax* data as the query dataset. Each cell was given a stage assignment based on the top matched cell from the reference index.

### scRNA-seq analysis of human transcriptomes

Genes x cell data matrices were called for 10000 cells in each sample and merged in the R computational environment (R v. 3.6.0) using the Seurat v. 3.2.3 computational package. Cell barcodes had to meet the following conditions to be retained: displaying a minimum of 300 transcripts and 200 genes (across the three genomes), with less than 5 percent murine mitochondrial transcripts and less than 30 percent of human mitochondrial transcripts (to account for the mitochondrial richness of hepatocytes). Additionally, log₁₀ ratios of human and murine-mapped transcripts per cell barcodes were calculated, and only cell barcodes with a ratio greater than 0 were retained. The resulting merged Seurat object was log-normalized to 10,000 transcripts (NormalizeData function), and variable features were selected using Seurat's FindVariableFeatures() function, picking 2,000 genes based on the variance-stabilizing transform. The merged object was scaled and centered using a linear model, with scale.max=10, block.size =1,000 and min cells.to.block=3,000. Principal component analysis was performed using RunPCA, retaining the top 12 components. The number of principal components retained was picked based on the inspection of the “elbow plot” and JackStraw procedure (n=100 replicates, 50 PCs) in Seurat. A UMAP embedding was calculated on the top 12 principal components using the RunUMAP() function with n_neighbors=50, min.dist = 0.3, spread=1, metric = “correlation”. Cluster identification was based on the Louvain algorithm with the top 12 principal components and 20 nearest neighbors at a resolution of 0.3. Gene and transcript coverage plots per cell barcodes were generated on the merged dataset prior to the final filtering step (for barcodes with a 10-fold enrichment for human transcript). Cluster markers were identified using the FindAllMarkers() function with logfc thresholds set to 0.58, min.diff.pct = 0.20, and test.use = “MAST”. For differential gene expression between infected and non-infected cells or between cells carrying AP2-G positive and AP2-G negative parasites the Find-Markers() function set to logfc thresholds = 0.25 and test.use = “MAST” was used. Functional class enrichment on differentially expressed gene sets was carried out using the Gene Ontology (GO) tool (Gene Ontology Consortium, 2021; Ashburner et al., 2011).

### Quantitative RT-PCR

Total RNA from pooled triplicate wells of *P. vivax*-infected MPCCs was extracted with TRIzol (Thermo Fisher), DNAse treated and purified using the RNeasy MinElute Cleanup Kit (Qiagen). cDNA synthesis was performed using SuperScript II (Thermo Fisher)
and RT-PCR was carried out using PowerUp SYBR Green Master Mix (Applied Biosystems) in a Roche Light Cycler 480 Real-Time PCR Detection System according to the manufacturer's instructions. The primers used are listed in Table S6. Relative gene expression was calculated with the delta-delta Ct method, using PVP01_1213400 as housekeeping gene.

## Immunofluorescence analysis

P. vivax-infected MPCCs were fixed in ice-cold methanol or 4% paraformaldehyde (PFA), washed in phosphate-buffered saline (PBS) and stored at 4°C. Parasites were detected using P. vivax-specific antibodies (PvUIS4 and PvBip) on methanol fixed cells as described in (Gural et al., 2018). For IFITM3 staining, PFA-fixed cells were permeabilized with 0.2% TritonX100 for 10 minutes at room temperature, washed in PBS and blocked with 2% bovine serum albumin (BSA) in PBS for 30 minutes at room temperature. IFITM3 rabbit monoclonal antibody was incubated overnight at 4°C (1:100). For parasite detection, cells were co-incubated with anti-PvCSP monoclonal antibodies. Alexa-conjugated 488 anti-rabbit (1:1000) and alexa-conjugated 635 anti-mouse (1:1000) antibodies were incubated for 1 hour at room temperature, followed by nuclear staining with Hoechst. Images were captured on a Nikon Eclipse Ti or Zen-ApoTome inverted wide-field microscopes using 20x objectives. To quantify parasite size, the area of parasite defined by the PvUIS4 staining was measured using NIS-Elements Microscope Imaging Software and automatically converted to equivalent diameter.

## Fluorescence in situ hybridization

P. vivax-infected MPCCs were fixed in 3.7% PFA for 10 minutes at room temperature, washed in PBS, immersed in 70% ethanol and stored at 4°C. Custom labelled probes set specific to *Pv18S rRNA* (FAM dye) and *PvGEXP5* (Quasar 670 dye) purchased from Stellaris were hybridized overnight in the dark at 37°C following the manufacturer's instructions. After nuclear staining and washing, cells were imaged in a Nikon Eclipse Ti fluorescence microscope by using 60× objectives. To quantify parasite size, the area of parasite defined by the *Pv18S rRNA* staining was measured using Fiji Software (Schneider et al., 2012) and converted to diameter.

## QUANTIFICATION AND STATISTICAL ANALYSIS

n represents the number wells from each plate as described in the figure legends. Exception for Figure 2C where n represents 2 independent infections. Methods used for computing statistical significance are indicated in figure legends. Statistical significance was considered for p values below 0.05. Data was analyzed using GraphPad Prism Software.