A single-cell liver atlas of Plasmodium vivax infection

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 proﬁling. 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 deﬁne host- and state-dependent transcriptional signatures and identify unappreciated populations of replicative and non-replicative parasites that share features with sexual transmissive forms. We ﬁnd 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. With this combined platform, we deﬁned distinct signatures between early-, dormant-, mid-, and late-stage parasites and identiﬁed sub-populations of sexually committed forms in the liver, which were previously thought to appear only during erythrocytic infection. We interro-gated 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. 3 days starting on day 5 after infection. dosing schedule was treatments with IFN alpha 2a, beta and gamma. E64 treatment at 4h post-infection media 30 percent of human mitochondrial transcripts (to account for the mitochondrial richness of hepatocytes). Additionally, log 10 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 identiﬁcation 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 ﬁnal ﬁltering step (for barcodes with a 10-fold enrichment for human transcript). Cluster markers were identiﬁed using the FindAllMarkers() function with logfc.threshold 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.threshold = 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 2021;


In brief
Plasmodium vivax parasites cause relapsing human malaria, which is considered a barrier to its eradication. Mancio-Silva et al. perform timeresolved, 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.

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 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-timepoint bulk sequencing used in the previous study prevented us from capturing the inherent heterogeneity of the distinct liverstage parasite forms as well as host-target and bystander cell responses. Achieving a deeper understanding of the pathogenhost 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 infectionspecific 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 (Howick 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 stagedependent 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 liverstage 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 nonreplicative 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 3 . 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 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 S1 and Tables S1, S2, and S4.
(MSP), serine repeat antigen (SERA), and rhoptry-associated membrane antigen (RAMA) and rhoptry 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 bloodstage 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 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 mm. Parasite numbers and size distribution of GEXP5-positive and -negative schizonts (mean ± SEM; n = 4 wells; t test: ****p < 0.0001). 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 , 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 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 fattyacid-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-Gpositive and AP2-G-negative P. vivax parasites. Interestingly, we found that the hepatocytes bearing non-replicative or  Figure S3 and Table S5.  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 hepato-cytes. 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

OPEN ACCESS
Resource 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 A B Figure 6. IFN responses control P. vivax hepatic infection (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.

Resource
Cell Host & Microbe 30, 1-13, July 13, 2022 9 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 hostand 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 IFNresponsive 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 adaptative 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(Liehl et al., , 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:  (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 naive 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 . 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 mM) 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-TARGETplus 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 . Briefly, MPCCs were first treated with B/B Homodimerizer (0.5 mM) for 30 minutes at 37 C to partially deplete the fibroblast cells via apoptosis . 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 mL water and hybridized for 32 hours as previously described (Gural et al., 2018). 10 mL 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