Unraveling the Plasmodium vivax sporozoite transcriptional journey from mosquito vector to human host

Malaria parasites transmitted by mosquito bite are remarkably efficient in establishing human infections. The infection process requires roughly 30 minutes and is highly complex as quiescent sporozoites injected with mosquito saliva must be rapidly activated in the skin, migrate through the body, and infect the liver. This process is poorly understood for Plasmodium vivax due to low infectivity in the in vitro models. To study this skin-to-liver-stage of malaria, we used quantitative bioassays coupled with transcriptomics to evaluate parasite changes linked with mammalian microenvironmental factors. Our in vitro phenotyping and RNA-seq analyses revealed key microenvironmental relationships with distinct biological functions. Most notable, preservation of sporozoite quiescence by exposure to insect-like factors coupled with strategic activation limits untimely activation of invasion-associated genes to dramatically increase hepatocyte invasion rates. We also report the first transcriptomic analysis of the P. vivax sporozoite interaction in salivary glands identifying 118 infection-related differentially-regulated Anopheles dirus genes. These results provide important new insights in malaria parasite biology and identify priority targets for antimalarial therapeutic interventions to block P. vivax infection.


Selection of media and study design to emulate insect and human host-like microenvironment.
The microenvironment composition was selected based on media components previously determined to resemble mosquito and human host-like environments that salivary gland sporozoites encounter upon transmission 34,35 . Schneider's insect media (pH 7) was used to mimic the mosquito microenvironment whereas incomplete RPMI 1640 media with the addition of 3% (w/v) of bovine serum albumin (BSA) represented the human microenvironment (Supplementary Table S1). Our study design focused on the serum protein albumin as it is the most abundant (>50%) of total serum proteins and has important roles in transport as well as modulating osmotic pressure 36 . Additionally, previous studies have shown that addition of albumin significantly contributes to sporozoite activation 8,35 . Further, we used BSA rather than human serum albumin as they are chemically similar, and to avoid serum-donor bias and non-specific protein binding. In total, four different buffered media compositions were tested including RPMI (B1), RPMI + 3% BSA (B2), Schneider's (B3), and Schneider's + 3% BSA (B4) (Fig. 1a).
In order to recapitulate the sporozoite's transitional journey from vector to host, our study design considered initial time zero (0 h) immediately after salivary gland isolation (manual dissection) with subsequent 4 h post-dissection incubations at two different temperatures, room temperature (RT) or 37 °C (Fig. 1a). We hypothesized that the prolonged 4 h post-dissection sampling would capture sporozoite gene expression specific to the microenvironment in which the sporozoite was maintained. For instance, transformation into early liver-stage parasites when held in a human-like microenvironment, as previously shown in Kaiser et al. 17 , or preservation of a quiescent transcriptional state when maintained in an insect-like microenvironment. Additionally, the 4 h post-dissection samples at RT or 37 °C were used to investigate the effect on the sporozoite when exposed to temporal temperature shift. Further, our salivary gland isolations typically consisted of 10-20 mosquitoes per microenvironment condition and required no additional sporozoite purification, reducing sample manipulation. Lastly, P. vivax salivary gland sporozoites (PvSGSs) in the various microenvironment conditions were phenotypically assessed at each time point using in vitro assays measuring percent gliding motility and invasion rate of PHHs (Fig. 1a).
All P. vivax patient isolate-derived samples (48 total) were analyzed for expression of genes showing low experimental variation indicated between replicates ( Supplementary Fig. S1a). These samples collected for RNA-seq included 5 biological replicates (patient isolates) for 0 h post-dissection, 3 biological replicates at 4 h post-dissection kept at RT, and 4 biological replicates at 4 h post-dissection kept at 37 °C. The quantified transcript abundance of PvSGS-infected samples showed a mean of 14% ± 9% (± standard deviation.) mapped to P. vivax reference genome P01 (PlasmoDB, version 34), a mean of 50% ± 12% mapped to the WRAIR2 An. dirus reference genome (Vectorbase, VB-2018-02, AdirW1.7), and a mean of 36% ± 5% of unmapped reads (mean transcript count ≥20 FPKM) ( Supplementary Fig. S1b). We achieved very high coverage of the P. vivax genome, where approximately 10-30% of the total reads (on average half a million malarial parasite reads per sample) are mapped back to the P. vivax reference sequences covering over 4000 genes with high confidence, confirming the high-quality preparation of our samples. Furthermore, an average of 70% of RNA transcripts detected from PvSGS samples and uninfected salivary gland controls mapped to An. dirus reference genome with, >93% of transcripts mapped to salivary gland-associated genes previously identified from an An. gambiae proteomic study 37 (Supplementary Fig. S1c).
Sporozoite transcriptome signatures linked to host microenvironments. Initially, sporozoite transcripts were identified through a comprehensive comparison of our sporozoites transcriptome to datasets previously published for P. vivax blood-stages 38 . We identified 3,510 up-regulated genes and 1,923 down-regulated genes in the sporozoite stage compared to blood-stages, of which 903 genes were identified as specific to sporozoites (Supplementary Fig. S2 and Supplementary Data S2). As a comparison of sporozoite to blood-stage gene expression across species, we additionally examined previously generated sporozoites transcriptome data from P. falciparum 39,40 , P. vivax 41 , and P. yoelii 42 (Fig. 1b). The datasets were re-normalized to a scale from 1 to 100, and similarity was calculated using a non-parametric method (Kendall's Tau). There are overall similarities between sporozoites data across different species with a conserved set of highly expressed genes, such as Cell Traversal protein for Ookinetes and Sporozoites (CelTOS) and Circumsporozoite Protein (CSP) (Fig. 1c). The overall correlation of transcriptome expression is stronger between the same stage of different species (i.e., sporozoite vs. sporozoite) than the correlation between different stages of the same species (i.e. oocyst or hypnozoite) (Fig. 1b).
Coupling our in vitro sporozoite bioassays allowed us to categorize sporozoite genes into three phenotype conditions: (i) quiescent; (ii) invasion-activated; and (iii) early liver-stage ( Fig. 2a and Supplementary Data S3). Genes preferentially expressed in salivary glands at RT in insect media without mammalian serum components (BSA) defined the quiescent phenotype. These quiescent-state genes were identified as genes most highly expressed in unperturbed PvSGSs by multiple pairwise comparisons of different microenvironment activators between the key insect-like and mammalian-like variables, including temperature (37 °C vs RT °C, log2FC < −1, P (p value) ≤ 0.05) (Supplementary Data S4) and insect or mammalian media (with or without BSA, log2FC < −2, P ≤ 0.05) (Supplementary Data S5). High transcript abundance implicated gliding motility proteins such as TRAP-like protein (TLP) and others (e.g., ROM1) as representative of the quiescent stage. Invasion-activated genes, which were genes preferentially expressed in the early phase (4 h vs. 0 h, log2FC < −2, P ≤ 0.05, Supplementary Data S6) and induced by BSA, included > 100 invasion-and motility-related genes (RPMI +3% BSA, 0 h vs. RPMI +3% BSA, 4 h, 37 °C, log2FC < −2, P ≤ 0.05, Supplementary Data S7). Many leading and attempted vaccine candidates are included in the invasion group, such as CelTOS, Apical Membrane Antigen 1 (AMA1) and Thrombospondin-related Adhesive Protein (TRAP). Finally, the early liver-stage is represented by genes with enriched expression at core human body temperature 37 °C (4 h, 37 °C vs. 4 h, RT, log2FC > 2, P ≤ 0.05, Supplementary Data S8). The early liver-stage is also the largest group, representing a major transcriptome shift when the motile parasite transitions to the intracellular microenvironment and initiates liver-stage parasite development. As expected for intracellular developmental requirements, several hundred up-regulated genes are related to protein biosynthesis, fatty acid synthesis, and mitochondrial function. To analyze these important parameters in our experimental setup, i.e. elapsed time (0 h and 4 h), host temperature shifts (RT and 37 °C), and host chemical environment change (with and without BSA), we calculated the overall transcriptome shift as indicated by Pearson's r between different experimental conditions. We used the transcriptomes of freshly-dissected sporozoites (0 h post-dissection) as a baseline of comparison and calculated the similarities between 4 h post-dissection with varying BSA and 37 °C conditions. We found a large overall change of transcriptome induced by temperature, with an r decrease to 0.2, when the parasite is maintained at 37 °C. The most significant transcriptome change occurred when both BSA and 37 °C are combined, with an r decrease of 0.3 to 0.5 (Fig. 2b).

An insect-like microenvironment improves sporozoite vitality and hepatocyte invasion.
To define the most suitable conditions to sustain viability and infectivity of freshly-dissected PvSGSs, pH and time were experimentally assessed using dissection into Schneider's insect media in comparison to RPMI. For all experiments, sporozoite viability and infectivity were measured using in vitro assays of sporozoite invasion of cultured PHHs and quantification of development of liver-stage parasites after 6-8 days post-infection. Collection of freshly dissected PvSGSs into Schneider's, pH 4 adversely affected sporozoite viability, resulting in no hepatocyte invasion and liver-stage parasite development, whereas Schneider's, pH 7 resulted in the highest rates of PHH infection and liver-stage parasite development (3-fold higher than RPMI) ( Supplementary Fig. S3a). When collected PvSGSs were maintained for 24 hours (at 4 °C or RT), only PvSGSs in Schneider's, pH 7 or Schneider's +3% BSA, pH 7 retained the ability to invade hepatocytes showing a ~2-fold increase in liver-stage parasites when the PvSGSs were maintained at RT (Supplementary Fig. 3b). These studies demonstrated that sporozoites-infected salivary glands collected into Schneider's, pH 7 are most efficient for hepatocyte infection, and these conditions are capable of maintaining sporozoite viability for at least 24 hours. The media Schneider's, pH7 was thus selected as the insect-like environment for our RNA-seq study and is referred to as "Schneider's" hereafter.
The PvSGSs samples collected for RNA-seq were also simultaneously subjected to bioassays measuring sporozoite motility and invasion of PHHs. We found PvSGSs exposed to Schneider's insect media showed a prolonged gliding motility phenotype across all time points (>60%) in comparison to RPMI which had a large reduction to <10% gliding after time 0 (Fig. 3a-e and Supplementary Fig. S4). Also, significant differences in infectivity and liver-stage parasite development resulted from exposure to the different microenvironments ( Fig. 4a-d). Of highest importance, PvSGSs dissected into Schneider's with subsequent addition of serum-containing hepatocyte A standard sandwich immunofluorescence assay (IFA) based on the monoclonal anti-circumsporozoite protein (CSP) antibody was used to visualize sporozoite motility where manual quantification (from 10 fields of view or 1,000 total PvSGSs) was used to measure sporozoite metrics (percent gliding, motility path). All images were captured at 20x magnification, 0.74 NA. Scale bars (white) represent 10 µM. Graph bars represent mean with s.d for experimental replicates (n = 10) and biological replicate (n = 3) from 10 fields of view or 1,000 PvSGSs. Statistical significance was calculated using a two-way ANOVA with Tukey's multiple comparisons test to all means where statistical significance values are represented as P < 0.05 (*), P < 0.001 (***), and P < 0.0001 (****). (Blue = RPMI, Purple = RPMI + 3% BSA, Green = Schneider's, Orange = Schneider's + 3% BSA).  115 ± 30 liver-stage parasites (Fig. 4b). Similarly, at 4 h post-dissection, 37 °C only PvSGSs in Schneider's and Schneider's + BSA mediums successfully infected PHHs (Fig. 4d). Immunofluorescence staining of the parasitophorous vacuole membrane (PVM) biomarker Upregulated in Infective Sporozoites 4 (UIS4) and active glycolysis biomarker Glyceraldehyde 3-Phosphate Dehydrogenase (GAPDH) were used to identify liver-stage parasites and distinguish between developing schizonts and hypnozoites ( Fig. 4e-f) 43 . All liver-stage parasites were imaged, classified, and quantified on a high-content imaging system (Operetta, Perkin Elmer) following an established methodology previously described 44 .
Identification of genes involved in sporozoite quiescent state. To discover the molecular metabolic activities of the sporozoite quiescent state, differentially expressed genes were identified between insect-like sporozoites at 0 h, (Schneider's, RT) with environmentally activated sporozoites (4 h post-dissection, RPMI + 3% BSA, 37 °C). Transcriptome profiling revealed 197 genes down-regulated upon entry into the human-like environment with the highest annotated differentially expressed genes summarized in Table 1 . Notably, our differential-expression analysis revealed genes involved in intracellular calcium signaling cascades, including Calcium-Dependent Protein Kinase (CDPK) 4, centrin (CEN) 2, and CDPK6, indicating proteins are maintained by sporozoites prior to stimulation from the mammalian host ( Fig. 2c) (FPKM ≥ 20, log2FC ≥ 1, P ≤ 0.05). Further, our results implicate genes with unexplored but hypothesized roles in sporozoites such as the recently characterized proteolytic cascade in P. falciparum blood-stage involving Plasmepsin X (PMX), subtilisin 2 (SUB2), and AMA1 used for invasion (FPKM ≥ 20, log2FC ≥ 1, P ≤ 0.05) 45 . In total, we identified 74 genes with unknown functions that we predict may have a role in the sporozoite pre-invasive, quiescent state. Among the top most highly expressed were PVP01_0943600, PVP01_0210600, and PVP01_0319000 (FPKM ≥ 20, log2FC ≥ 2, P ≤ 0.05) (Supplementary Data S3).
The function of the Plasmodium apicoplast is indispensable for the parasites survival, but the organelle's function is not fully defined. The Suf family is one of the four known pathways for Fe-S cluster biogenesis and is believed to provide the Fe-S clusters to many proteins functional in the apicoplast, including proteins involved in the isoprenoids biosynthesis pathway. Our analysis indicated an increased abundance of apicoplast genes, suggesting the organelle may have a role in both quiescent sporozoites (SufE) and early liver-stage parasite formation. To validate RNA-seq data, four genes of the Suf family (SufE, SufA, SufS, SufC) and the most conserved protein of the isoprenoid pathway (IspD) 46 were analyzed by qPCR (n = 2) which confirmed SufE and SufS are actively expressed ( Supplementary Fig. S5a). Further, high-resolution fluorescent microscopy of P. vivax sporozoites using mouse immune sera raised against proteins SufS and SufE stained positive for SufS in a single sporozoite organelle, similarly to an Acyl Carrier Protein (ACP) control with known localization to the apicoplast. Alternatively, SufE had a diffused staining pattern with localization in multiple sporozoite organelles ( Supplementary Fig. S5b). These data indicate the potential importance of this pathway during sporozoite and liver-stages, warranting further investigation.

Specific human factors induce sporozoite transition into an early liver-stage parasite.
After a sporozoite invades the host hepatocyte, formation of a PVM ensues for transition into a highly metabolically-active intracellular liver-stage parasite preparing for rapid replication. To explore the early liver-stage parasite transcriptomic profile, an analysis of differentially expressed genes was performed using our PvSGSs exposed to the mammalian-like environment (4 h post-dissection, 37 °C). We identified 518 differentially expressed genes noting 30 highly abundant annotated differentially expressed genes in Table 3 (Supplementary Data S3). Interestingly, the ookinete surface protein p25 (Pvs25), potentially a false-positive, followed by the Parasitophorous Vacuole protein 1 (PV1) showed the highest differential expression with temperature change and BSA induction. Another intriguing find was the expression of Sortilin, recently characterized in P. falciparum asexual blood-stage as a transmembrane protein involved in protein trafficking to rhoptries of the apical complex, suggesting a role in the sporozoites secretory system during and after host hepatocyte invasion 49,50 . Additional genes implicated in initial liver-stage parasite development include Liver Stage Associated Protein 1 (LSAP1), Sporozoite Invasion-Associated Protein 2 (SIAP2) and four PVM-associated genes (PVP01_0504800, PVP01_1271000, PVP01_073480, PVP01_0602100) termed Early-Transcribed Membrane Proteins (ETRAMPs), and UIS4. Activation of sporozoite intracellular calcium signaling cascades in response to the mammalian microenvironment was evident with increased abundance of calmodulin, CDPK7, protein kinase 4 (PK4), and Receptor of Activated Protein C Kinase 1 (RACK1). Our results also support recent reports that heat shock proteins (HSP) are critical for liver-stage parasite development as HSP70, HSP70-2, HSP90, and PVP01_1405400 were amongst the highest differentially expressed genes 51 . In order to rapidly transition between pivotal life-cycle stages, Plasmodium uses translational repression to regulate pre-transcribed mRNA stockpiles that encode required proteins. Our analysis of differentially expressed genes revealed high expression of known translational regulators ALBA-1 52 , ALBA-2, and ALBA-4 53 indicating these genes as potential global regulators of developing liver-stage parasites (FPKM ≥ 20, log2FC ≥ 2, P ≤ 0.05).

Characterization of the Anopheles dirus salivary gland transcriptome. Experimental PvSGS-
infected An. dirus samples and uninfected An. dirus salivary gland controls were analyzed for expression of genes on day 14 (post-blood meal). Transcriptome analysis showed similar expression levels for salivary gland proteins in all samples including salivary gland protein 6 (SGP6), D7 long-form salivary protein, salivary gland protein 1-like, and TRIO salivary gland protein (Supplementary Data S9). The top expressed known gene encoded was anopheline antiplatelet protein (AAPP), which is secreted by blood-feeding females to block platelet adhesion to collagen and prevent aggregation 54 . To explore the transcriptome of PvSGS-infected An. dirus salivary glands, differentially expressed genes were compared between infected and uninfected salivary glands. The results showed 71 genes up-regulated and 47 genes down-regulated differentially expressed genes, but only 43% of these differentially expressed genes have a known function (uninfected vs infected, log2FC ≥ 2 or ≤−2, P ≤ 0.05) (Fig. 5, Table 4, Supplementary Data S10). The top three genes with significantly excess abundance were uncharacterized (ADIR008862, ADIR014872, and ADIR010223) while the most highly-expressed known genes (ADIR007940 and ADIR007941) are members of the HSP20/alpha-crystallin family which have not been previously reported. As expected, results showed high induction of immune response-related genes, such as the anti-microbial proteins (AMPs): gambicin, defensin, cecropin, attacin, and leucine-rich immune protein (LRIM1) (Fig. 5) 55 . Further assignment of biological relevance using Gene Ontology (GO) analysis revealed differences in gene expression of biological processes where infected salivary glands showed an increase of protein catabolism, apoptosis inhibition, structural composition, and olfaction. On the other hand, down-regulation was observed in genes relating to proton and amino acid transport, protein synthesis, and DNA repair (Fig. 5, Table 4) 56 . Interestingly, enrichment in gene expression for uninfected mosquitoes showed increased expression of three genes (ADIR014650, ADIR003224, and ADIR004887) annotated as salivary gland proteins with significantly lower levels of expression in PvSGS-infected salivary gland samples (Supplementary Data S10).

Discussion
Transmission of Plasmodium sporozoites from mosquito to human involves a transition between two highly distinctive environments, insect to mammal, requiring rapid adjustment for parasite survival. Inside the mosquito's salivary gland, the sporozoite exists in a quiescent state primed for introduction to the mammalian host but not yet activated. Upon inoculation, sporozoites must immediately react to changes in their microenvironment, especially temperature and serum components, to initiate gliding motility and migration to the liver. In conjunction with this activation, changes in gene expression reflect future developmental requirements for sporozoite cell traversal and then parasite development within hepatocytes. In this study, we used microenvironments with extrinsic stimuli to mimic in vitro infection process to define key relationships with distinct biological functions. Most notable, post-dissection preservation of quiescence by placing sporozoites in an insect-like microenvironment and delaying strategic activation until inoculating hepatocyte cultures prevented untimely activation of invasion-associated genes to dramatically increase invasion rates into primary human hepatocytes. Additional studies will be needed to determine the specific components that modulate sporozoite activation. Our experimental method allowed us to examine factors of the different biotic microenvironments in multiple pair-wise comparisons (insect vs mammal, BSA vs no BSA), time (0 h vs 4 h), and temperature (RT vs 37 °C) where differential expression revealed 3 phenotypic states (i) quiescent; (ii) invasion-activated; and (iii) early liver-stage. We found that across different Plasmodium species the sporozoite transcriptome have strong similarities with a conserved set of highly expressed genes such as CelTOS, CSP, and TRAP (Fig. 1b,c). However, this study is the first to differentiate between activated vs non-activated sporozoites in relation to key host physiological triggers accompanied with complete transcriptomic analysis (Fig. 2) and quantitative phenotypic assessment (Figs 3, 4). Previous studies in rodent malaria's and P. falciparum have under-explored the host physiological effects on sporozoites, required laborious purification, and are unable to link expression with in vitro phenotypes 39,40,42 . Within the last decade, several blood-stage P. vivax sequencing projects in genomics 57,58 , transcriptomics 38,59,60 , and proteomics 61,62 have advanced our insight of parasite biology 63 , however, P. vivax sporozoite biology is underrepresented with only one recent publication on the P. vivax sporozoite proteome 63,64 . Our analysis provides more understanding of transitional host effects as the first published P. vivax sporozoite RNA-seq study and addresses how to improve in vitro infectivity and intracellular development of primary human hepatocytes.
The RNA-seq analysis revealed 197 genes expressed by quiescent PvSGSs maintained in an insect-like environment including genes associated in motility (TLP, TRAP, profilin, CLAMP), cell traversal (GEST, CelTOS, SPECT1), and host invasion (CSP, SIAP1, SPATR). Further, our analysis suggests genes with roles in blood-stage motility and invasion such as Plasmepsin X45 and merozoite TRAP-like protein (MTRAP) 65,66 to potentially have uncharacterized functions in sporozoites or early liver-stage parasites. In agreement with the RNA-seq results, our bioassays revealed the most notable finding that PvSGSs collected in Schneider's insect media remain quiescent until activation by exposure to hepatocyte culture media (Fig. 3a-b). Strategic activation by serum-containing   67,68 , or temperature 9 to cause an intracellular calcium flux followed by exocytosis 69 . In agreement, our results demonstrate increased expression of 16 genes including calmodulin, CDPK1, SSP3, SPELD, and regulatory proteins when sporozoites are activated by BSA. CDPK1 orchestrates many processes required for invasion of host cells, including microneme release and activation of actin-myosin motor in PvSGSs 70 , while SSP3 and SPELD are sporozoite surface proteins with functions in motility 71 and liver-stage parasite maturation 72 . Based on these results it is evident that premature activation of PvSGSs in mammalian-based mediums mistimes the just-in-time cascade of events leading to reduced sporozoite viability and infectivity. The transcriptome shift from PvSGSs into early liver-stage parasites, mimicked experimentally by exposure to BSA and human body temperature, identified 517 genes with different levels of expression in the transition. Our GO and KEGG enrichment analysis associated the genes with translation and DNA replication with implications in glycolysis and pentose phosphate pathways. As expected, the differentially expressed genes included many liver-stage-associated genes previously determined to be essential for establishment and development (PV1, ETRAMPs, LSAP1) 73 . Much speculation exists of what genes are present in early liver-stage parasite development where our results more conclusively reveal through the up-regulation of genes involved in the secretory pathway, such as Sortilin, and regulatory processes of gene expression, such as translational repressors ALBA 1, 2, and 4. Furthermore, we attempted to discern the role of the Plasmodium apicoplast in sporozoites as the organelle is vital for the parasites survival due to the prokaryotic type metabolic pathways functional within. The importance of Fe-S cluster biogenesis has been deemed essential for the apicoplast maintenance in the erythrocytic stages of P. falciparum 74,75 and in oocyst development of P. berghei 76 but yet to be studied in the pre-erythrocytic stages. Our RNA-seq analysis and qRT-PCR determines that initial components of the Suf pathway (Suf S and E) are actively up-regulated in infectious sporozoites and early liver-stage parasites where diffused immunofluorescence staining of SufE alludes to the possibility of multiple roles for the protein as previously reported in Arabidopsis thaliana 77 . Lastly, our analysis helps prioritize development of new vaccine candidates, such as SIAP2, by its link to hepatocyte invasion and subsequent development 9,78,79 , and new drug targets in signaling cascade, such as calmodulin, CPDK7, PK, and RACK1.
We offer the first RNA-seq analysis of uninfected and Plasmodium-infected An. dirus, a major SE Asia vector for malaria. While the mosquito salivary gland main purpose is to produce saliva and digestive enzymes to breakdown ingested nutrients, it must tolerate the presence of PvSGS for extended periods in order for transmission to be successful. Our study indicates the presence of the parasite has a profound effect on the An. dirus with induction of innate immune responses genes to P. vivax infections similar to other studies 55,[80][81][82][83][84][85] ; however, we identified many genes not yet associated with Plasmodium infections, including the odorant binding protein (OBP) ADIR010042, the serine protease (SP24D) ADIR007695, and the small heat shock associated proteins  (alpha-crystallin B chain) ADIR007940 and ADIR007941. Heat shock proteins are essential to mosquito survival and induced by biotic factors such as temperature or colonized parasite [86][87][88] . Presumably, higher expression in P. vivax-infected mosquitoes may contribute to the sporozoites need to mediate stress and thermotolerance when transitioning between hosts during multiple mosquito blood meals. Further, HSP20 has specifically been implemented in regulation of sporozoite adhesion and locomotion 89 . Additionally, we found down-regulation of three salivary gland proteins (ADIR014650, ADIR003224, and ADIR004887) in P. vivax-infected salivary glands previously unreported. The parasite manipulation of salivary gland proteins may be a strategy used to bypass the activation of the host innate immune response and avoid detection. Overall, our identification of novel mosquito biomarkers can increase understanding of transmission, improve metrics for measuring transmission, and act as targets for transmission blocking. Research of P. vivax usually lags significantly behind that of P. falciparum primarily due to the inability to continuously culture P. vivax in the laboratory. However, our research in this study with P. vivax sporozoites provide innovative changes to broadly enhance malaria sporozoite research. Our experimental analyses of P. vivax sporozoites revealed that the microenvironmental changes associated with transition from mosquito into human host trigger distinct processes for initial activation, hepatocyte infection, and liver-stage parasite development. Importantly, distinct calcium-dependent protein kinase signaling pathways associated with these phases of infection represent attractive new drug targets that might be capable of halting migrating sporozoites before they can infect and hide within the liver. Additionally, we offer the first study on mosquito salivary gland gene and proteins involved in P. vivax infection and transmission. Altogether our study represents an important advance in understanding mosquito transmission and the key processes regulating sporozoite infectivity of P. vivax.

Mosquito colonization, rearing and P. vivax infection. For colonization, 3-5-day old adult An. dirus
was provided blood meals preserved using citrate phosphate dextrose adenine (CPDA-1) via an artificial membrane feeder containing human blood (Interstate Blood Bank, Inc., Memphis, TN), as previously described 90 . After blood-feeding, females were manually mated, eggs obtained, and larvae and adults maintained as described earlier 91 . The blood-fed adults were maintained in the AFRIMS insectary at 25 °C ± 1 °C and 80% relative humidity until transport to the Mae Sod Malaria Clinic, Mae Sod District (Thailand Ministry of Public Health. A 5% solution (v/v) of commercially produced multivitamin syrup (MULTILIM, Atlantic Pharmaceutical Co., Bangkok, Thailand) that included 5% sucrose (LIN, TTR Group, Col, Uthai Thani, Thailand) was as a food source and changed daily 92 .
The P. vivax -infected blood was obtained from Thai patients who reported to medical clinics along the Thai-Myanmar border with signs and symptoms of malaria infection and confirmation of vivax malaria parasites in the blood. Thick and thin blood smears were prepared and then examined microscopically for P. vivax parasites and gametocyte densities calculated per 200 white blood cells. In total, 20 ml of venous blood was collected in a heparinized tube and maintained at a constant temperature of 37 °C using a water bath circulation system. After starvation for 8 h, 5-7-day old An. dirus females were allowed to blood-feed for 30 min using an artificial membrane feeder. Unfed and partially fed mosquitoes were removed, counted, and discarded in accordance with standard protocols. Engorged mosquitoes were then secured in a container and transported to the AFRIMS insectary where they remained until use at AFRIMS or shipment to USF (Tampa, FL). The adults were maintained at temperature and humidity as described earlier and provided with 10% multivitamin sugar solution that was changed daily to reduce fungal and bacteria growth. P. vivax infected An. dirus mosquitoes were shipped to USF for experiments RNA-seq experiments (n = 1) and preliminary experimental optimization ( Supplementary Fig. S3).  Supplementary Table S1. The microenvironments were termed as the following: B1 = RPMI, B2 = RPMI + 3% BSA, B3 = Schneider's, B4 = Schneider's + 3% BSA. Time points were termed as the following: 0 h refers to time after dissection was finished and sporozoites were taken off ice (4 °C) where 4 h refers to time after dissection was finished and sporozoites are incubated at respective temperatures. Temperature differences at 4 h post-dissection occurred at either room temperature (RT) or 37 °C, 5% CO 2 . The 4 h RT with 37 °C induction refers to a 30 min incubation (the amount of time to run the gliding assay) of 4 h RT samples with 37 °C, 5% CO 2 incubation.

Dissection of P. vivax salivary gland sporozoites.
For all experiments, salivary glands were dissected from An. dirus mosquitoes 14-16 days post-infectious blood meal. Mosquitoes were immobilized by cold incubation (−20 °C, 5 min), transferred to ice (4 °C), and individually washed in a 3-step wash: 70% ethanol, 0.5 mg ml −1 penicillin-streptomycin and 1 mg ml −1 of neomycin (PSN), and 2.5 µg ml −1 Fungizone TM diluted in 1× Dulbecco's Phosphate Buffered Solution (DPBS). Intact, clean salivary glands were collected into 100 μl of collection media with dissection time under an hour. Salivary glands were centrifuged (1600 G, 3 min) with gland disruption using plastic pestle and manually pipetting then counted using a hemocytometer. A total of 200,000-2,000,000 sporozoites were collected per condition for RNA-seq.
Sporozoite gliding motility assay. The gliding motility assay used was described in a previous study 34 .
Briefly, 20,000-40,000 sporozoites were allowed to glide for 30 min in respective media at either room temperature (23 °C) or 37 °C on a glass coverslip coated with 10 µg ml −1 of P. vivax anti-CSP monoclonal antibodies (mAB) 2F2 (subtype VK210) or 2E10.E11 (subtype VK247) (BEI Resources, NIAID, NIH) 93 . Subsequently, sporozoites were fixed with 4% paraformaldehyde (PFA) for 10 minutes and blocked with 1% bovine serum albumin (BSA) in 1× phosphate buffered saline (PBS) overnight at 4 °C. Sporozoites were stained with 10 µg ml −1 of mAB as primary antibody for 1 hour at room temperature followed by staining with 10 µg ml −1 secondary goat anti-mouse . Cryopreserved PHH were thawed following manufacturer's instructions with post-thaw viability measured by trypan blue exclusion with > 5 million viable cells per vial and diluted in HCM supplemented with a final concentration of 50 µg ml −1 penicillin-streptomycin, 100 µg ml −1 of neomycin, and 10 µg ml −1 of gentamicin (HCM + ) to a concentration of 900 cells per μl. The cell mixture was seeded as a confluent monolayer with a total of 18,000 cells per well in a 384-well plate then incubated at 37 °C, 5% CO2. Following 2 days post-seed, each well received a complete media change (40-50 µl) with HCM + then ready for sporozoite infection.
Microscopy and image analysis of P. vivax liver-stage parasites. In brief, a standard immunofluorescence assay (IFA) protocol was used to stain and quantify P. vivax liver-stage parasites. Wells were incubated in blocking buffer (0.03% TritonX-100, 1% (w/v) BSA in 1× DPBS) with rabbit anti-UIS-4 polyclonal antibody or mouse anti-ACP antibody at 1:1,000-fold dilution (kindly provided by Center for Infectious Disease Research, WA, U.S.A) and mouse anti-GAPDH monoclonal antibody at 1:50,000-fold dilution (The European Malaria Reagent Repository, Cat. No. 7.2) overnight at 4 °C 43,94 . After, wells were washed thrice with 1× DPBS and incubated for 1 hour at room temperature with secondary goat anti-mouse or anti-rabbit Alexa Fluor ® 488 or 568 conjugate at 1:1,000-fold dilution (2 µg ml −1 ; Cat. No. A11001, A11034, A11004, A11011, Molecular Probes, ThermoFisher, Waltham, MA, U.S.A.) antibody and Hoechst at 1:1,000-fold dilution (10 µg ml −1 ) then washed thrice and filled with 1× DPBS until imaging and for long-term storage. Liver-stage parasite imaging, identification, and quantification was performed on the Operetta Imaging System with Harmony software 4.1 (20x, 0.4 NA; Perkin Elmer, Waltham, MA, U.S.A.) following previously described protocol 44 . Each individual well of a 384-well plate had 35 fields of view collected, capturing anti-UIS4 signal using FITC channel and hepatocyte nuclei using DAPI channel. Predetermined image analysis parameters were set to collect the total liver-stage parasite population. Large, developing liver-stage schizonts and small, dormant hypnozoites were distinguished based upon fluorescent intensity, area, and roundness. Statistical analysis of PHH invasion assays were performed using Prism 7 (Graphpad, La Jolla, CA, U.S.A.).
High-resolution IFA images of the P. vivax liver-stage parasites were captured with a 100x oil objective, 1.4 NA on a DeltaVision Core system (GE Healthcare Life Sciences, Piscataway Township, NJ, U.S.A.). The same IFA protocol was followed as described above. Mouse immune sera samples collected from mice raised against protein to SufE and SufS were diluted in blocking buffer at 1:50 or 1:200 (kindly provided by V.S. and G.S.S., Z.P. and S.G). All high-resolution images were deconvoluted using the softWoRx ® image analysis package (GE Healthcare Life Sciences, Piscataway Township, NJ, U.S.A.) and merged z-stacks with added scale bars were processed in (ImageJ) 95 . Images captured with secondary antibody Alexa Fluor ® 568 were re-colored to magenta using FIJI in order to meet publication image standards. were then measured on a Nanodrop (ThermoFisher, Waltham, MA, U.S.A.). 0.5 μg-1.0 μg total RNA of each sample was used for library preparation. Libraries were prepared by using TruSeq Stranded mRNA Preparation Kit v2 (Illumina, San Diego, CA, U.S.A.) following manufacturer's recommendations. Library quantification was conducted by qPCR and TapeStation (Agilent Technologies, Santa Clara, CA, U.S.A.) measurements. Libraries were subsequently pooled and sequenced across 4 or 12 lanes of a MiSeq at a read length of 100 bp paired end using 300 cycle V2 MiSeq reagent kit (Illumina, San Diego, CA, U.S.A.).
Mapping, data processing, and transcriptome profiles. The sequencing raw reads from each sporozoite sample were aligned to the P. vivax reference genome P01 58 of PlasmoDB version 34 96 . The sequencing raw reads from each infected and uninfected An. dirus salivary gland samples were aligned to the An. dirus WRAIR2 reference genome of VectorBase release VB-2018-02 (AdirW1.7) 56 . A maximum of one mismatch per read was allowed. The mapped reads from TopHat 97 were used to assemble known transcripts from the reference and their abundances were estimated using Cufflinks 98 . The expression level of each gene was normalized as Fragments Per Kilobase of exon per Million (FPKM) mapped reads for each condition. All the RNA sequencing raw reads have been deposited into NCBI's Gene Expression Omnibus which are accessible through GEO Series accession number GSE117538.
Differential gene expression analysis. The difference of gene expression levels across the different time points and buffered media conditions were determined from fragment counts using EdgeR 99 in R software 100 . The aligned read pairs from sorted bam files were counted against P. vivax P01 (PlasmoDB v34) and An. dirus WRAIR2 (AdirW1.7) gene annotations using featureCounts application of the Subread suite 101 , defining features at the gene-level. TMM (trimmed mean of M-values) normalization 100 of the counts was performed across the five runs and three biological replicates to eliminate compositional biases between the samples. For each pairwise comparison (Supplementary Dataset S4-S9) performed, a gene was inferred differentially expressed if the FPKM ≥ 20, log2FC of the mean FPKM across biological replicates is greater 2 or less than −2, P ≤ 0.05 and FDR (False Detection Rate) < 0.1.

Gene Ontology (GO) term enrichment analysis.
Predetermined gene lists from each pairwise differential gene expression comparison were used to test the enrichment of specific GO term in each group. The GO enrichment analysis of the molecular function, cellular component, and biological processes were conducted using the Gene Ontology (GO) tool on PlasmoDB 96 (www.plasmodb.org) and analyzed using a P of ≤0.05 cutoff. REVIGO program 102 was used to reduce the redundancy in GO terms. GO terms with a frequency ≥20% were reported and fold enrichment was used to generate figures.
Gene expression comparison of human and mosquito stages of P. vivax. Gene expression comparison between recently published P. vivax blood-stage 38 and our current sporozoite stage transcriptomic data, to understand potential underlying variations and identify sporozoite specific genes (Supplementary Dataset S3, Supplementary Fig. S2). Raw reads from sporozoite samples are re-aligned 98 to P. vivax reference genome Sal1 of PlasmoDB version 12 to be consistent with blood-stage RNA-seq data. Peak gene expression across the sporozoite (buffered media and time point) and blood (time point) samples are respectively profiled along with peak expression point, to classify the genes and rank genes based on sporozoite stage expression.
Quantitative real time RT-PCR (qRT-PCR). Total RNA was extracted from samples preserved in Trizol as described above. After giving DNase I (Invitrogen, Carlsbad, CA, U.S.A.) treatment, the total RNA (about 50-100 ng) was subjected to first-strand cDNA synthesis using the Superscript II Reverse Transcriptase kit (Invitrogen, Carlsbad, CA, U.S.A.) and gene-specific primers following manufacturer's protocol. The resulting cDNA was diluted 1:5 with nuclease-free water. For real time analysis, the primers were designed based on mRNA of desired genes available from PlasmoDB (Supplementary Table S2) using freely available software GeneRunner and cross verified using Primer Blast (NCBI). Real-time PCR analysis was performed on Agilent Technologies Stratagene Mx3005P System using the Brilliant II SYBR Green QPCR Master Mix with low ROX (Agilent Technologies, Santa Clara, CA, U.S.A.). Briefly, the PCR reaction consisted of 12.5 µl of Brilliant II SYBR Green QPCR Master Mix, 10 µmole of forward and reverse primers and 2 µl of diluted cDNA in a total volume of 25 µl. PCR cycling conditions were performed using the default conditions of the Mx3005P Software: 95 °C, 10 minutes followed by 40 cycles of 94 °C, 30 seconds; 55 °C, 1 minute; and 72 °C, 30 seconds. The housekeeping P. vivax gene, seryl-tRNA-synthetase, known to be transcribed stably throughout different erythrocytic stages was used as an endogenous control (normalizer) for every run. Two buffered media compositions (RPMI + 3% BSA, Schneider's) were used from 2 biological samples and run in triplicates for the following timepoints, 0 h post-dissection and 4 h post-dissection, 37 °C. For data analysis, the normalized quantity of each target gene was expressed as the ratio of the relative amount of target gene over the quantity of the housekeeping gene. Then, log2FC was used calculated from qPCR data (Schneider's, 0 h post-dissection vs. RPMI + 3% BSA, 4 h post-dissection, 37 °C) with results compared to the RNA-seq generated FPKM (Supplementary Fig. S5a).
Statistical analysis. All graph bars are represented by means with standard deviation (s.d.). Significance of mean was assessed by two-way ANOVA with Tukey's multiple comparisons (Figs 3 and 4) When standard deviation is reported, n represents biological replicates and experimental replicates per biological replicate. When p value is reported, P is used for representation.
Disclaimer. Material has been reviewed by the Walter Reed Army Institute of Research. There is no objection to its presentation and/or publication. The opinions or assertions contained herein are the private views of the SCIENTIFIC REPORTS | (2018) 8:12183 | DOI:10.1038/s41598-018-30713-1 authors, and are not to be construed as official, or as reflecting true views of the Department of the Army or the Department of Defense. The investigators have adhered to the policies for protection of human subjects as prescribed in AR 70-25.