Coordinated host-pathogen transcriptional dynamics revealed using sorted subpopulations and single, Candida albicans infected macrophages

The control of fungal infections depends on interactions with innate immune cells including macrophages, which are among the first host cell types to respond to pathogens such as Candida albicans. This fungus is a member of the healthy human microbiome, although also a devastating pathogen in immunocompromised individuals. Consistent with recent findings from studies of other pathogens, we observed that within a population of interacting macrophages and C. albicans, there are distinct host-pathogen subpopulations reflecting cell specific trajectories and infection outcomes. Little is known about the molecular mechanisms that control these different fates. To address this, we developed an experimental system to isolate the major host-fungal pathogen subpopulations observed during ex vivo infection using fluorescent markers. We separated subpopulations of macrophages infected with live C. albicans from uninfected cells and assessed the variability of gene expression in both host and fungal pathogen for each subpopulation across time using RNA-Seq. In infected cells, we observed a coordinated, time-dependent shift in gene expression for both host and fungus. The early response in macrophages was established upon exposure to C. albicans prior to engulfment and involved up-regulation of pathways and regulatory genes required for cell migration, pathogen recognition, activation of engulfment, and phagocytosis; this pro-inflammatory response declined during later time points in parallel with expression changes in C. albicans. After phagocytosis, the initial response of C. albicans was to up-regulate genes related to survival in the nutrient-limited and stressful environment within macrophages; at later time points, gene expression shifted to initiate hyphal growth and escape. To further probe the heterogeneity seen observed in host-pathogen interactions, we performed RNA-Seq of single macrophages infected with C. albicans. We observed that some genes show higher levels of heterogeneity in both host and fungal pathogen cells that we could not detect in subpopulation samples; we observed that the time shift in expression is asynchronous and that expression changes in both the host and pathogen are tightly coupled. This work highlights how analysis of subpopulations and single host-pathogen pairs can resolve population heterogeneity and trace distinct trajectories during host interactions with fungal pathogens.

infected with live C. albicans from uninfected cells and assessed the variability of gene expression in 23 both host and fungal pathogen for each subpopulation across time using RNA-Seq. In infected cells, we 24 observed a coordinated, time-dependent shift in gene expression for both host and fungus. The early 25 response in macrophages was established upon exposure to C. albicans prior to engulfment and 26 involved up-regulation of pathways and regulatory genes required for cell migration, pathogen 27 recognition, activation of engulfment, and phagocytosis; this pro-inflammatory response declined during 28 later time points in parallel with expression changes in C. albicans. After phagocytosis, the initial 29 response of C. albicans was to up-regulate genes related to survival in the nutrient-limited and stressful 30 environment within macrophages; at later time points, gene expression shifted to initiate hyphal growth 31 and escape. To further probe the heterogeneity seen observed in host-pathogen interactions, we 32 performed RNA-Seq of single macrophages infected with C. albicans. We observed that some genes 33 show higher levels of heterogeneity in both host and fungal pathogen cells that we could not detect in 34 subpopulation samples; we observed that the time shift in expression is asynchronous and that 35 expression changes in both the host and pathogen are tightly coupled. This work highlights how 36 macrophages and sorted using fluorescent activated cell sorting (FACS) at time intervals (0,1, 2 and 4 103 hours) to isolate different host-fungal pathogen subpopulations and single-cell pairs ( Figure 1A). These 104 time points were selected to capture the rapid transcriptional changes of C. albicans in response to 105 macrophages 3 . Primary bone derived macrophages were stained with CellMask Deep Red plasma 106 membrane stain prior to sorting. To examine gene expression profiles, both host and fungal RNA were 107 extracted and adapted for Illumina sequencing using Smart-Seq2 (Methods). Both microscopy and 108 FACS analysis revealed that four distinct infection subpopulations can be isolated ( Figure S2A): (i) 109 macrophages infected with live C. albicans (GFP+, mCherry+, Deep red+), (ii) macrophages infected 110 with dead C. albicans (GFP-, mCherry+, Deep red+), (iii) macrophages exposed to C. albicans (GFP-, 111 mCherry-, Deep red+) and (iv) C. albicans exposed to macrophages (GFP+, mCherry+, Deep red-; 112 Figure 1A). 113

114
The number of RNA-Seq reads and transcripts detected for both host and fungal pathogen 115 subpopulations was sufficient for differential expression analysis in most samples and to widely profile 116 parallel transcriptional responses. We aligned reads to a composite reference of both mouse and C. 117 albicans transcriptomes (Methods) and found that the fraction of mapped reads for host and pathogen 118 was correlated with the percent of sorted cells for each subpopulation (Figures 1B, S2B; Table S1). In 119 subpopulations containing both host and pathogen, the fraction of reads averaged 87% macrophage 120 and 13% C. albicans for macrophages infected with live fungal cells, and 95% macrophage and 5% C.  Table S1). Fewer transcripts were detected in 130 subpopulations of macrophages infected with dead C. albicans (an average of 3,214 host transcripts 131 and 983 fungal transcripts Figure S3A) and had modestly correlated biological replicates (e.g. 132 Pearson's r < 0.56). Lower coverage was expected for this subpopulation, as only up to 3% of each 133 sample collected at each time point was comprised of macrophages infected with dead C. albicans 134 ( Figure S2B). Therefore, we primarily focused the differential expression analysis on subpopulations of 135 macrophages or C. albicans exposed, and macrophages infected with live C. albicans, which had high 136 transcriptome coverage, and highly correlated biological replicates (e.g. Pearson's r 0.96 and 0.92 in 137 macrophages and Candida at 4 hours, respectively; Figure S3B). These results indicate that we have 138 established a robust system for measuring host and fungal pathogen transcriptional signal during 139 phagocytosis in sorted infection subpopulations. We next examined how C. albicans gene expression varied across unexposed, exposed and 159 phagocytosed cells over time. Using k-means clustering, we identified sets of genes with similar 160 expression patterns; the major patterns of expression across time were either induced (cluster 1, 2 and 161 3) or repressed (cluster 4 and 5) in the live, phagocytosed subpopulation relative to all other C. albicans 162 infection fates (Figures 2B, S4A). A large number of genes were differentially expressed in 163 phagocytosed C. albicans (732 genes relative to unexposed across all time points) compared with the 164 number of DEGs found in exposed but un-engulfed C. albicans (82 genes relative to unexposed across 165 all time points). Comparing the phagocytosed and un-engulfed C. albicans subpopulations at each time 166 point, the major differential response was found at 1 hour, highlighting a rapid and specific 167 transcriptional response upon macrophage phagocytosis. Many of these genes maintained high 168 expression levels throughout the 4-hour infection time course (Table S2; Figures 2B, 2E). 169 170 Genes highly induced in phagocytosed C. albicans are involved in adaptation to the macrophage 171 environment. This highlights major changes in metabolic pathways; these genes (cluster 1; Figure 2B, 172 2E) are involved in glucose and carbohydrate transport, carboxylic acid and organic acid metabolism, 173 and fatty acid catabolic processes (enriched GO terms corrected-P < 0.05, hypergeometric distribution 174 with Bonferroni correction; Table S3). Prior microarray analysis of C. albicans exposed to macrophages 175 reported that similar changes in metabolism and nutrient uptake allow Candida to utilize the limited 176 spectrum of nutrients available in the phagosome 3 . With RNA-Seq data and sorted infection fates, we 177 observed upregulation of additional genes related to these functions and confirmed which expression 178 changes were specific to the phagocytosed C. albicans subpopulation (Table S2). Genes involved in 179 glyoxylate metabolism, the beta-oxidation cycle and transmembrane transport were significantly 180 induced in phagocytosed C. albicans relative to exposed cells. By contrast, multiple classes of 181 transporters were highly up-regulated in both engulfed and exposed C. albicans subpopulations, 182 including oligopeptide transporters, several high affinity glucose transporters, and amino acid 183 permeases, suggesting these changes are not in response to phagocytosis (Figure 2B; Tables S2 and 184 S3). Cluster 1 includes genes involved in pathogenesis and genes associated with the formation of 185 hyphae, including core filamentous response genes (ALS3, ECE1, HTG2, ORF19.2457; Table S2) 17 . 186 Overall, gene expression levels in cluster 1, including filamentation genes, increased over the time 187 course with the highest expression at 4 hours, regardless of whether C. albicans cells were 188 phagocytosed or remained un-engulfed, suggesting that these genes were induced by the presence of 189 macrophages. While media containing serum can also induce C. albicans filamentation, we found that 190 those genes were more highly induced upon phagocytosis (Table S2; Figure S4B). We also identified 191 genes induced in live, phagocytosed C. albicans and repressed in un-engulfed cells relative to the 192 unexposed subpopulation (cluster 2; Figure 2B). Cluster 2 includes several secreted aspartyl 193 proteases and additional genes involved in transmembrane transport (Opt and Hgt classes; Tables S2 194 and S3). These transporters differ with those found in cluster 1, as the scale of their induction after 195 phagocytosis was more modest. Other sets of genes up-regulated during C. albicans phagocytosis in 196 cluster 3 had lower induction relative to cluster 1 or cluster 2. This set included genes related to 197 oxidation-reduction processes, including dehydrogenases, mitochondrial respiratory response, and 198 transcription factors (enriched GO terms, corrected-P < 0.05, hypergeometric distribution with 199 Bonferroni correction; Figure 2B; Table S3). Cluster 1, 2, and 3 encompassed the major transcriptional 200 modules up-regulated in C. albicans upon phagocytosis. 201

202
Other sets of genes were specifically down-regulated in live, phagocytosed C. albicans at the earliest 203 time point, whereas expression levels of these genes did not change in non-phagocytosed C. albicans 204 (both unexposed and exposed) over the time course (cluster 4 and 5; Figures 2B, 2E). In cluster 4, 205 which displays the strongest signature of repression, genes down-regulated in phagocytosed C. 206 albicans recovered their expression levels at 4 hours. This cluster includes genes related to the 207 translation machinery and peptide biosynthesis, including ribosomal proteins, chaperones and 208 transcription factors that regulate translation (enriched GO terms, corrected-P < 0.05, hypergeometric 209 distribution with Bonferroni correction; Table S3; Figure S4B). Repression of the translation machinery 210 was previously noted as a response of C. albicans to macrophage interaction 3 . Here, we demonstrated 211 that down-regulation of ribosomal proteins, chaperones and translation-regulator transcription factors is 212 specific to phagocytosed C. albicans and that expression of these genes recovered at later time points 213 ( Figures 2E, S4B). Cluster 4 also encompassed highly repressed yeast-phase specific genes, 214 including those involved in ergosterol biosynthesis, cell growth and cell wall synthesis (Figures 2B,  215   S4B; Table S2). Cluster 5, which was not as strongly repressed as cluster 4, includes a set of genes 216 down-regulated in phagocytosed C. albicans that did not recover expression over the length of the time 217 course (Table S2). These genes are largely involved in nucleoside metabolic processes and host 218 adaptation (enriched GO terms, corrected-P < 0.05, hypergeometric distribution with Bonferroni 219 correction; Figure 2B; Table S3). A subset of these down-regulated genes were involved in 220 morphological and cell surface remodeling, including three essential negative regulators of filamentation 221 (SSN6, MFG1, and TUP1) 18-21 . This subset also includes additional repressors of filamentation 22 and 222 transcription factors that regulate the white-opaque phenotypic switch (WOR1 and WOR2) 23,24 and 223 commensalism in the mammalian gut (WOR1) 25 (Figures 2B, S4). These results highlight that a large 224 part of the observed transcriptional variation is likely due to the down-regulation of these transcriptional set 2), which grouped into four clusters with similar expression patterns. PC analysis showed that 243 exposed and infected macrophages clustered closely together, separately from unexposed 244 macrophages along the PC1 and PC2, indicating that exposed and infected macrophages had similar 245 transcriptional responses, which primarily varied across time in the first and second PCs (cluster 1 and 246 2; Figure 2C). This suggests that the early (1 to 4 hours) host transcriptional response to C. albicans is 247 set during exposure and maintained during phagocytosis. For both exposed and infected macrophages, 248 a major difference was found between 2 and 4 hours along PC2, which highlights genes that were 249 differentially induced or repressed at 4 hours in these subpopulations (clusters 3 and 4; Figure 2C,D). 250

251
In all four clusters, infected and exposed macrophages have similar patterns of differential expression 252 relative to unexposed macrophages ( Figure 2D). Overall, these clusters were significantly enriched in 253 activation, migration, phagocytosis, and triggering the innate immune response; this includes the  (Figures 2D, S5). 262 263 In exposed and infected macrophages, many of the genes induced at 1 hour remained relatively stable 264 at 2 and 4 hours (cluster 1; Figure 2D). These genes are related to defense mechanisms such as pro- A second set of genes was up-regulated in both exposed and infected macrophages, with peak 278 expression at 2 hours in exposed macrophages (cluster 2; Figure 2D; Table S4). This set of genes is 279 associated with pathogen recognition, opsonization, and activation of the engulfment (P-value < 0.05 280 right-tailed Fisher's Exact Test; Table S5), including Clec7a (also known as Dectin-1), lectin-like 281 receptors such as galectin 1 (Lgals1) and galectin 3 (Lgals3; Figure 2D; Tables S4 and S5), and other 282 transmembrane receptors (Cd36, Cd74, Fcer1g, Ifnar2, Igsf6, Itgb2, Gm14548, Trem2; Table S4 and 283 S5). In addition, several complement proteins were found more highly induced in exposed 284 macrophages, including complement factor properdin (Cfp), extracellular complement proteins (C1qa, 285 C1qb, and C1qc), and complement receptor proteins (C3ar1 and C5ar1). These results highlight the 286 importance of C. albicans recognition and activation of engulfment in the exposed macrophage 287 subpopulation. Since these genes maintained high expression in infected macrophages, they may also 288 play an important role during phagocytosis or allow for uptake of additional C. albicans cells. 289

290
We also examined subpopulations of macrophages infected with dead C. albicans; this data was 291 analyzed separately, as the total number of cells sorted and therefore the transcriptome coverage was 292 low (i.e. 31% relative to macrophages infected with live C. albicans; see above). While we did not have 293 sufficient data to examine C. albicans transcripts (i.e. 21% transcriptome coverage relative to exposed 294 and live-phagocytosed subpopulations; see above), our analysis of the host transcripts revealed a small 295 set of highly induced genes. This set includes pro-inflammatory cytokines such as Ccl3, Cxcl2, Cxcl14, 296 Il1rn, and Tnf, and transcription regulators such as Cebpb, Irf8, and Nfkbia ( Figure S8). These genes 297 were also induced in macrophages infected with live C. albicans (clusters 1 and 2; Figure 2D), 298 indicating that maintaining expression of these genes may be important for pathogen clearance after 299 phagocytosis. 300 301 Another major shift in gene expression occurred at 4 hours, with sets of genes highly repressed or 302 highly induced at this later time point (cluster 3 and 4, respectively; Figures 2D, 2F) in both the 303 exposed and infected macrophage subpopulations. Highly repressed genes at 4 hours (cluster 3) are 304 enriched in cytokines (Il1a, Irf4) and transmembrane receptors, including intracellular toll-like receptors 305 (Tlr5 and Tlr9), and interleukin receptors (Il1r1, Irf4, Il17ra). This cluster was also enriched in categories 306 associated with proliferation and immune cell differentiation (P-value < 0.05 right-tailed Fisher's Exact 307 Test; Table S5). By contrast, macrophage genes highly induced at 4 hours (cluster 4) were enriched in 308 categories related with phagocytosis and programmed cell death mechanisms, including the chemokine 309 Cxcl10, and transcriptional regulators that play a role in inflammation and programmed cell death 310 (Cebpd, Fos, Irf2bp1, Irf8, Irf9, Nfkbie, Card9, Bcl2; (P-value < 0.05 right-tailed Fisher's Exact Test; 311 Table S5). This suggests that during phagocytosis of C. albicans there was a strong shift toward a 312 weaker pro-inflammatory transcriptional response by 4 hours. Our approach not only identified genes 313 specifically induced in infected macrophages, but also highlighted expression transitions related to time 314 of C. albicans exposure. These patterns were conserved among subpopulations of exposed and 315 successfully isolate single infected macrophages via FACS, we cannot control for the number of C. 329 albicans cells inside of each macrophage with this approach, since macrophages phagocytose variable 330 numbers of C. albicans cells 12 (Methods). We obtained 4.03 million paired-end reads per infected cell 331 on average; a total of 449 single, infected macrophages had more than 1 million paired-end reads, over 332 99% of which passed our sequence quality-control filters ( Figure S9A; Table S1). For macrophages 333 with live C. albicans, we found an average 75% of reads mapped to host transcripts and 11% of reads 334 mapped to C. albicans transcripts (Figure 3A). Although parallel sequencing of host and pathogen 335 decreases the sensitivity to detect both transcriptomes, the number of transcripts recovered (> 1 336 Transcripts Per Million, TPM) for macrophages and C. albicans was sufficient for cell clustering and  We identified differentially expressed genes (corrected-P < 0.05) using a likelihood-ratio test (LRT) for 364 single-cell differential expression 33 (Methods). The transcriptional response among single infected 365 macrophages exhibited two time-dependent states (state 1 M and state 2 M ) associated with expression 366 shifts from 2 to 4 hours (Figure 4, top; Figure S11A). While cells were largely separated into these two 367 states by time, a few macrophages were assigned to the alternate cell state by unsupervised clustering 368 and appeared to be either early or delayed in the initiation of the transcriptional shift in the immune 369 response. Immune cells that display asynchrony in their transcriptional response have previously been 370 reported in LPS-stimulated dendritic cells, where a few "precocious" cells expressed high levels of 371 immune response genes early, leading to cytokine secretion and activation of other cells in the 372 population 30 . However, precocious host cell expression can also be related to pathogen state 14 . In 373 single infected macrophages, we found that 88 and 70 differentially expressed genes (likelihood-ratio 374 test (LRT), P < 0.001) that showed higher expression in the single, infected macrophages assigned to 375 state 1 M , and assigned to state 2 M , respectively (Table S6) (Figure 4, bottom). Notably, independent analysis of the parallel 390 fungal transcriptional response identified two pathogen stages that were also primarily distinguished by 391 time. In live phagocytosed C. albicans, we found a total of 168 differentially expressed genes 392 (likelihood-ratio test (LRT), P < 0.001), 80 and 86 genes showed higher expression in C. albicans 393 phagocytosed by macrophages assigned to state 1 M , and state 2 M , respectively (Table S7). In C. 394 albicans, genes in state 1 C were enriched in organic acid metabolism (P = 6.63e-11; enriched GO term, 395 corrected-P < 0.05, hypergeometric distribution with Bonferroni correction; Table S8), including a strong 396 upregulation of transporters (HGT13) and glyoxylate cycle genes, specifically those from beta-oxidation 397 metabolism (ECI1, FOX3, FOX2, PXP2; Figure 4, bottom). Most macrophages infected by C. albicans 398 in state 1 M induced a strong pro-inflammatory response (co-state 1; Figure 4). At 4 hours, expression 399 of these genes was reduced and we observed an up-regulation of genes enriched in carbon 400 metabolism (P = 2.43e-05; enriched GO term, corrected-P < 0.05, hypergeometric distribution with 401 Bonferroni correction; Table S8), including genes related with a shift to glycolysis and gluconeogenesis 402 (PGK1), fatty acid biosynthesis (FAS1, ACC1), and genes associated with filamentation, such as ECE1, 403 HWP1, OLE1, RBT1; whereas the majority of infected macrophages had down-regulated expression of 404 pro-inflammatory cytokines (co-state 2; Figure 4). In summary, these 2 host-pathogen co-states largely 405  Table S9). 420 Highly expressed (top 5%) unimodal genes with similar expression profiles at 2 and 4 hours 421 encompassed genes involved in the immune response to C. albicans infection, including genes 422 enriched in pathogen recognition, opsonization, and activation of engulfment (cluster 2; Figure 2D), 423 such as complement proteins (C1qb, C1qc) and galectin receptors that recognize beta-mannans 424 (Lgals1, Lgals3; Figure S11). In phagocytosed C. albicans, most genes (an average of 76% of the 425 genes detected; Methods) also had unimodal expression patterns ( Figure S11; Table S10). We also 426 found a set of unimodal genes with similar expression profiles at both time points, which were enriched 427 in the oxidation-reduction process and defense against reactive oxygen species (enriched GO term, 428 corrected-P < 0.05, hypergeometric distribution with Bonferroni correction; Table S10 Table  442 S9). In addition, we further identified differential distributions (e.g. shifts in mean(s) expression, 443 modality, and proportions of cells) across and within 2 and 4 hours as implemented in scDD package 34 . 444 We found a subset of immune genes that showed differential distribution (P < 0.05, Benjamini-445 Hochberg adjusted Fisher's combined test). For example, the Olr1 lectin-like receptor, the tumor 446 necrosis factor receptor superfamily member 12a (Tnfrsf12a), and the transmembrane marker Cd83, 447 had bimodal expression distribution at 2 hours but not at 4 hours; the interleukin 4 receptor, alpha 448 (Il4ra) and the transmembrane marker Cd300a had bimodal distribution at 4 hours but not at 2 hours; 449 and the interleukin 21 receptor (Il21r) had differential distribution and is bimodal at both 2 and 4 hours 450 (Table S9). Prior work demonstrated that macrophages infected with Salmonella and stimulated with 451 LPS also exhibited heterogeneity and bimodal expression patterns in subsets of genes, including some 452 immune response genes 14 . Similar to macrophages infected with Salmonella and stimulated with LPS, 453 Tnf and Il4ra exhibited bimodal expression patterns in single macrophages infected with C. albicans 14 454 ( Figure S13, top). Additionally, we found unique subsets of genes displaying differential distributions in 455 single macrophages infected with C. albicans, but not in response to Salmonella or LPS, including 456 Il17ra and other lectin-like receptors ( Figure S14). These results highlight heterogeneity in gene 457 expression patterns among single macrophages infected with live C. albicans that might indicate 458 distinct cell levels and trajectories in response to fungal infection and suggest that some expression 459 heterogeneity in macrophages is pathogen specific, likely a result of variably expressed pathogen-460 specific receptors. 461 462 We next hypothesized that C. albicans may also demonstrate expression heterogeneity and bimodality 463 that is correlated with expression in the corresponding macrophage cell. We found an average of 23% 464 C. albicans genes that showed patterns of bimodality at 2 and 4 hours (Table S10; Figure S11). We 465 also further examined if genes in phagocytosed C. albicans exhibited differential distributions (e.g. 466 shifts in mean(s) expression, modality, and proportions of cells) across and within 2 and 4 hours using 467 scDD 34 . We found a subset of virulence-associated genes that showed differential distribution (P < 0.05, 468 Benjamini-Hochberg adjusted Fisher's combined test). Genes showing bimodal expression only at 2 469 hours were enriched in regulation of defense response (enriched GO term, corrected-P < 0.05, 470 hypergeometric distribution with Bonferroni correction; Table S10), including genes associated with C. inflammatory cytokine production and phagocytosis (e.g. Tnf, Cxcl10, Irf1). Other sets of genes shifted 496 in gene expression in subpopulations of macrophages exposed to or infected with live C. albicans from 497 2 to 4 hours, including up-regulation of genes required for the activation of the inflammosome (e.g. 498 Nod2, Card9), and the down-regulation of pro-inflammatory cytokines (e.g. Il1a, Il1r1, Il17ra) and 499 intracellular receptors (e.g. Tlr5, Tlr9; Figures 2C, 2D). Live phagocytosed C. albicans within these 500 macrophages showed a rapid shift in gluconeogenic metabolism, amino acid uptake, cell wall 501 remodeling, and initiation of filamentation (Figure 2A, B). Up-regulation of gene clusters increased over 502 the time course with the highest expression at 4 hours, and down-regulation of gene clusters recovered 503 their expression levels at 4 hours, when macrophages up-regulated inflammosome-associated genes 504 and down-regulated of pro-inflammatory cytokines and intracellular receptors (Figure 2). 505 Transcriptional variability in these subpopulations could be further resolved into distinct transcriptional 506 states and cell trajectories at the single-cell level that govern infection fate decisions. We observed that 507 the gene expression in single macrophages infected with live C. albicans at 2 and 4 hours was tightly 508 coupled with that in phagocytosed C. albicans and could be described as two, time dependent host-509 pathogen co-states. At 2 hours, most single infected macrophages induced genes related to a pro-510 inflammatory response and the phagocytosed C. albicans upregulated transporters and glyoxylate cycle 511 genes, specifically those involved in beta-oxidation metabolism (co-state 1; Figure 4). However, a 512 subset of the cells at two hours was more similar to the major co-state found at 4 hours, where 513 expression of these transporters and metabolic C. albicans genes were reduced and we observed an 514 up-regulation of genes related with filamentation, such as ECE1, ALS3, OLE1, RBT1; in parallel, 515 macrophages down-regulated pro-inflammatory genes (co-state 2; Figure 4). Notably, we observed 516 expression heterogeneity of some infection-induced genes that was only detected at the level of single interactions. This also revealed that expression heterogeneity in key genes in both infected 552 macrophages and in phagocytosed C. albicans may contribute to infection outcomes. We identified two, 553 time dependent co-states of host-fungal pathogen interaction. The initial state is characterized by 554 induction of a pro-inflammatory host profile after 2 hours of interaction with C. albicans that then 555 decreased by 4 hours. A previous study indicated that a balance between pro-inflammatory and anti-556 inflammatory responses is necessary for C. albicans to establish infection 35 . A decrease in the pro-557 inflammatory response after C. albicans exposure was also found in human macrophages, where pro-558 inflammatory macrophages that interact with C. albicans for 8 hours or longer skew toward an anti-559 inflammatory proteomic profile 36 . In addition, recent single-cell RNA-Seq analysis of macrophages 560 containing growing Salmonella showed that macrophages shifted to an anti-inflammatory state by 20 561 hours; here, the authors hypothesized that fast-growing intracellular Salmonella overcame host 562 defenses by reprogramming macrophage gene expression 15 . We found that single macrophages 563 infected with C. albicans assigned to co-state 2 decreased expression of pro-inflammatory cytokine 564 genes and upregulated genes involved in the activation of inflammasomes; this response was coupled 565 with activation of filamentation and cell-wall remodeling in phagocytosed C. albicans. Previous work has 566 shown that C. albicans can escape from macrophage phagosomes by lytic or non-lytic mechanisms. 567 eliminated so as to increase fitness of a particular genotype that is fit for that environment. A similar 581 scenario may provide an advantage for phagocytes that have to hedge their bets when they encounter 582 pleomorphic C. albicans. These strategies have been described within microbial populations, where a 583 small fraction of "persister" cells might be transiently capable of surviving exposure to lethal doses of 584 antimicrobial drugs as a bet-hedging strategy 40 . While gene expression patterns are correlated in host 585 and fungal pathogen, it is unclear whether expression heterogeneity among individual infected 586 macrophages results from or results in expression bimodality among phagocytosed C. albicans. 587 588 Our approach represents a novel method to query host-fungal pathogen interaction. The 589 characterization of genes involved in heterogeneous responses is important to consider in the selection 590 of novel antifungal drug targets, as designing therapeutics to target products of genes expressed 591 uniformly among fungal cells in a population may be more effective than targeting the products of 592 genes expressed by only a subset of cells. Additionally, designing combination immunotherapies to 593 affect genes discriminative of distinct infection subpopulations, such as those genes upregulated in 594 macrophages infected with dead C. albicans, could help to increase the proportion of host cells that 595 clear the fungal pathogen. Parallel host-fungal pathogen expression profiling could also allow 596 researchers to not only measure the effects of new drug treatments on the pathogen, but also collect 597 information on how these treatments affect host cells. As single-cell RNA-Seq microfluidic platforms 598 continue to develop and become more cost effective to profile thousands of single cells, it will become 599 tractable to interrogate larger cohorts of host cells infected with fungus. This approach will allow 600 investigation of fungal phenotypic heterogeneity as a driver of different host responses and provide a 601 systems view of these interactions. The reporter construct used in this study was prepared by integrating the GFP and mCherry fluorescent 607 tags driven by the bi-directional ADH1 promoter and a nourseothricin resistance (NAT R ) cassette at the 608 Neut5L locus of Candida albicans strain SC5314 ( Figure S15). Briefly, the pUC57 vector containing 609 mCherry driven by ADH1 promoter (Bio Basic) was digested and this portion of the plasmid was ligated 610 into a pDUP3 vector 41 containing GFP, also driven by ADH1 promoter, a NAT R marker, and homology 611 to the Neut5L locus. The resulting plasmid was linearized and introduced via homologous 612 recombination into a neutral genomic locus, Neut5L using chemical transformation protocol with lithium 613 acetate. Transformation was confirmed via colony PCR and whole-genome sequencing. A whole 614 genome library was created from strain SC5314-Neut5L-NAT1-mCherry-GFP using Nextera-XT 615 (Illumina) and sequenced on Illumina's Miseq (150x150 paired end sequencing). Sequencing reads 616 were de novo assembled using dipSPAdes 42 . Scaffolds were aligned back to the plasmid used to 617 transform wild type SC5314 using BLAST 43 . Sequencing coverage was visualized using IGV 44 . counts paired-reads were aligned to the 'composite reference' as described above. 709

710
For each single macrophage infected with C. albicans, we quantified the number of genes for which at 711 least one read was mapped (TPM > 1). We filtered out low-quality macrophage or C. albicans cells from 712 our data set based on a threshold for the number of genes detected (a minimum of 2,000 unique genes 713 per cell for macrophages, and 600 unique genes per cell for C. albicans, and focused on those single 714 infected macrophages that have good number of transcript detected in both host and pathogen ( Figure  715 S9A). For a given sample, we define the filtered gene set as the genes that have an expression level 716 exceeding 10 TPM in at least 20% of the cells. After cell and gene filtering procedures, the expression 717 matrix included 3,254 transcripts for the macrophages and 915 transcripts for C. albicans. To estimate 718 the number of C. albicans in each macrophage, we measured the correlation of GFP levels from FACs 719 with the total number of transcripts detected in live, phagocytosed C. albicans cells (at least 1 TPM) but 720 found only a modest correlation between these two metrics (R 2 = 0.52). 721

722
To eliminate the non-biological associations of the samples based on plate based processing and 723 amplification, single-cell expression matrices were log-transformed (log(TPM + 1)) for all downstream 724 analyses, most of which were performed using the R software package Seurat 725 (https://github.com/satijalab/seurat). In addition, we do not find substantial differences in the number of 726 sequenced reads and detected genes between samples. We separately analyzed two comparisons of

Detection of variable genes and cell clustering 742
To classify the single cell RNA-Seq from macrophages and C. albicans, the R package Seurat was 743 used 32 . We first selected variable genes by fitting a generalized linear model to the relationship between 744 the squared co-efficient of variation (CV) and the mean expression level in log/log space, and selecting 745 genes that significantly deviated (p-value < 0.05) from the fitted curve, as implemented in Seurat, as 746 previously described 15 . Then highly variable genes (CV > 1.25; p-value < 0.05) were used for principle 747 component analysis (PCA), and statistically significant determined for each PC using JackStraw. 748 Significant PCs (p-value < 0.05) were used for two-dimension t-distributed stochastic neighbor 749 embedding (tSNE) to define subgroups of cells we denominated host-pathogen co-states. We identified 750 differentially expressed genes (corrected-p < 0.05) between co-states using a likelihood-ratio test (LRT) 751 for single-cell differential expression 33 as implemented in Seurat. 752 753

Detection of differential expression distributions and bimodality 754
To detect which genes have different expression distributions single infected macrophage and live C. 755 albicans we compared the distributions of gene expression within single infected macrophage and live 756 C. albicans, across and between 2 and 4 hours, and identified genes showing evidence of differential 757 distribution using a Bayesian modeling framework as implemented in scDD 34 . We used the permutation 758 test of the Bayes Factor for independence of condition membership with clustering (n permutations = 759 1000), and test for test for a difference in the proportion of zeroes (testZeroes=TRUE). A gene was 760 considered differentially distributed using Benjamini-Hochberg adjusted p-values test (p-value < 0.05). The projection score (red: high; blue: low) for each gene (row) onto PC1 and PC2 revealed four clusters 818 in both host and pathogen. For each projection score cluster, immune response genes (macrophages) 819 and functional biological categories (C. albicans; deducted from GO term enrichment analysis) are 820 shown. (b) Transcriptional response in subpopulations of Candida albicans. Heatmap depicts 821 significantly differentially expressed genes for replicates of each C. albicans sorted populations 822 (unexposed, exposed and phagocytosed) at 0, 1, 2 and 4 hours post-infection, grouped by k-means 823 (similar expression patterns) in five clusters. Each cluster includes synthesized functional biological 824 albicans response). Color scheme of gene expression is based on z-score distribution from -1.5 850 (purple) to 1.5 (yellow). Bottom and right margin color bars in each heat-map highlight co-state 1 (dark 851 blue in macrophages, and dark green in C. albicans) and co-state 2 (blue in macrophages, and green in 852 C. albicans), and time post-infection 2 hours (gray) and 4 hours (dark gray). (b) Violin plots at right 853 illustrate expression distribution of a subset six differentially expressed immune response or immune 854 evasion genes for each co-state in macrophages and C. albicans, respectively. 855 856 857