key: cord-0930551-5oossfn9
authors: Steuerman, Yael; Cohen, Merav; Peshes-Yaloz, Naama; Valadarsky, Liran; Cohn, Ofir; David, Eyal; Frishberg, Amit; Mayo, Lior; Bacharach, Eran; Amit, Ido; Gat-Viks, Irit
title: Dissection of Influenza Infection In Vivo by Single-Cell RNA Sequencing
date: 2018-06-27
journal: Cell Systems
DOI: 10.1016/j.cels.2018.05.008
sha: 7a4dfe6bb349a27e0cf9dbb7f98daa944b3f1650
doc_id: 930551
cord_uid: 5oossfn9

Summary The influenza virus is a major cause of morbidity and mortality worldwide. Yet, both the impact of intracellular viral replication and the variation in host response across different cell types remain uncharacterized. Here we used single-cell RNA sequencing to investigate the heterogeneity in the response of lung tissue cells to in vivo influenza infection. Analysis of viral and host transcriptomes in the same single cell enabled us to resolve the cellular heterogeneity of bystander (exposed but uninfected) as compared with infected cells. We reveal that all major immune and non-immune cell types manifest substantial fractions of infected cells, albeit at low viral transcriptome loads relative to epithelial cells. We show that all cell types respond primarily with a robust generic transcriptional response, and we demonstrate novel markers specific for influenza-infected as opposed to bystander cells. These findings open new avenues for targeted therapy aimed exclusively at infected cells.

Correspondence ido.amit@weizmann.ac.il (I.A.), iritgv@post.tau.ac.il (I.G.-V.) In Brief Simultaneous mapping of host and viral transcriptomes at the same single cell provides a framework to study host-viral interactions. Analysis of cells derived from lungs of in vivo influenza infection reveals both generic and cell-typespecific infection response. By analyzing the cellular heterogeneity of both bystander (exposed but uninfected) and infected cells, we characterize novel markers specific for influenza-infected as opposed to bystanders. These findings open new avenues for targeted therapy aimed exclusively at infected cells.

Influenza infection is a major cause of severe morbidity and mortality. Understanding the interplay between the influenza virus and its host cells is critical for the development of successful therapeutic approaches. In numerous studies, the host response to influenza infection has been characterized by measuring bulk cell populations and their further experimental validation. Collectively, these have provided rough models of the host response (e.g., Altboum et al., 2014; Shapira et al., 2009 ), but it is clear that the in situ infected lung is far more complex than the current models. For example, while epithelial cells are known to be the main targets of the influenza virus, several studies have documented that other cell types, such as endothelial cells, natural killer (NK) cells, macrophages, and dendritic cells (DCs), are also susceptible to the influenza virus, with potential implications of intracellular infection for their functionality (Manicassamy et al., 2010; McFadden et al., 2009) . Complexity may also be related to a wide range of viral transcriptional states within the infected cells (Russell et al., 2018; Xin et al., 2018; Zanini et al., 2018) , as well as to the heterogeneity of host-response states (Avraham et al., 2015) . Another complication may be attributable to the dual role of the metabolic machinery in supporting the host while also limiting the energetic demands of the viral life cycle (Jovanovic et al., 2015; Kissig et al., 2017) . Yet another source of heterogeneity may derive from the possibility that only a subset of the cells are actually infected by influenza, while most of the cells in the population are bystanders (exposed but uninfected) and typically respond to defensive host signals such as type I interferons (IFNs) (McFadden et al., 2009) .

Multiple key questions related to the complexity of the in vivo influenza infection have yet to be answered. In particular, the extent and nature of intracellular infection in different cell types has not been systematically elucidated. Furthermore, systematic characterization of host-response heterogeneity in bystander and infected cells across various cell types remains uncharacterized. Advances in the fluorescent labeling of influenza viruses have opened the way to an unbiased separation of infected and bystander cells (De Baets et al., 2015; Manicassamy et al., 2010) . A potential source of limitation with this technique is that reporter proteins cannot be expressed in the influenza virus without changing its infectivity, probably because of the significant burden of the reporter on viral fitness (Breen et al., 2016) . As a result, it has been a challenging task to characterize in vivo the full repertoire of viral-host interactions using these strains. An opportunity to tackle this challenge was provided by recent advances in single-cell genomic technologies (Jaitin et al., 2014; Paul et al., 2015) , allowing for simultaneous mapping of both the host and the viral transcriptome in the same single cell.

Here, we report the first comprehensive single-cell mRNA profiling of influenza-treated lungs. We characterized the host and the viral transcriptomes for 7,325 single cells from nine immune and non-immune cell types in both wild-type and Irf7knockout (KO) mice. We found that substantial proportions (22%-62%) of all cell types are infected by the influenza virus, and showed that infected cells predominantly have low intracellular viral transcriptome loads. Exceptional in this respect are epithelial cells, which are infected with massive viral loads (reaching up to 20% of all transcripts). Using Irf7-KO mice, which are incapable of launching an effective type I IFN response, we found that infectivity is a basic cellular property that is independent of the type I IFN host response. Moreover, our analysis highlights a core module of type I IFN genes (''IFN module'') that is shared by nearly all cell types. By comparing infected and bystander cells in both wild-type and Irf7-KO mice, we showed that this module is predominantly affected by Irf7dependent bystander exposure, consistently with the established IFN feedback loop. Another notable finding was that influenza infection also leads to a generic repression of mitochondrially encoded genes and that this response is largely uncoupled from induction of the IFN module: unlike the IFN module, the mitochondrial module responds to both Irf7-dependent and Irf7-independent signals and manifests its peak of repression only after intracellular viral invasion.

Our results open the way to interventions targeted specifically on infected cells, suggest principles of viral-host interactions that are applicable in influenza and possibly also in other pathogens, and highlight the power of simultaneous single-cell measurements of both hosts and viral transcriptomes in delineating a comprehensive map of in vivo infection.

Dissecting In Vivo Influenza Infection Using Combined Single-Cell Mapping of Host and Viral Transcriptome To simultaneously study, in an unbiased way, both host and viral transcriptional states after influenza infection, we used singlecell RNA sequencing (scRNA-seq) of cells from mouse lungs (female C57BL/6J mice aged 5.5-6.5 weeks) obtained 2 days after influenza infection. We first used fluorescence-activated cell sorting (FACS) to isolate immune (CD45 + ) and non-immune (CD45 À ) cells derived from the lungs of influenza-treated and PBS (''control'')-treated mice. Next, we used a massively parallel single-cell RNA-seq approach (MARS-seq; Jaitin et al., 2014) that profiles all polyadenylated mRNAs (STAR Methods, Figure 1A) . A total of 2,034 and 2,146 cells were analyzed from influenza-treated and control mice, respectively. To evaluate the quality of the data, we identified genes that are differentially expressed in cells isolated from the influenza-treated mice compared with the control (binomial test, STAR Methods). As expected, these genes (P binomial <10 À5 ) were overrepresented in antiviral-related signaling pathways, such as IFN signaling (P hg <10 À8 ) and Irf7 signaling (P hg <10 À7 ).

Clustering of single cells enabled us to distinguish between different cell types while taking into account the possible confounding effect of antiviral-related programs. We did this by first excluding genes that correlated strongly with the top two prin-cipal components (top 1.5% of principal-component-correlated genes), which were strongly enriched for antiviral, anti-proliferative, and inflammatory processes ( Figure S1A ). Next, we applied the expectation maximization clustering algorithm as previously described (Paul et al., 2015) , and manually annotated each cluster based on the expression of canonical markers. The analysis produced nine robust clusters with a total of 4,064 cells, which comprised four clusters of non-immune cell types: epithelial cells (EP), lymphatic and blood endothelial cells (LEC and BEC, respectively), and fibroblasts (FIB), and five specialized immune cell types: granulocytes (GN), T cells, B cells, NK cells, and a joint cluster of dendritic cells, monocytes, and macrophages collectively termed the mononuclear phagocyte system (MPS) (Figure 1B,  Table S1 , STAR Methods). To determine whether cells were indeed clustered by cell type rather than by treatment or by intracellular viral gene transcription, each individual cell was annotated according to its best match with the global expression of previously sorted bulk cell subsets. As shown in Figures S1B and S1C, the annotation of individual cells based on the bulk transcriptome data was in agreement with the marker-based annotation of these clusters. We further used t-distributed stochastic neighbor embedding (t-SNE) to project the different cell subpopulations onto a 2-dimensional space. As expected, cells annotated with the same cell type formed a separate t-SNE cluster independently of treatment (influenza-treated versus control) ( Figures 1C and 1D ).

Since the viral mRNA (vmRNA) is polyadenylated, MARS-seq captures both viral and host mRNAs within each individual cell ( Figure 1A ). We inferred the viral load within a cell from the amount of unique molecular identifiers (UMIs) that are aligned with viral segments. Specifically, the viral load was defined as the proportion of viral UMIs (that is, UMIs aligned with one of the viral segments) in the total UMI content of a given cell. Accordingly, we defined three groups of cells: (1) ''infected cells,'' cells whose viral load is higher than a certain selected viral-load cutoff; (2) ''bystanders,'' cells that are derived from the influenzatreated mouse but do not detectably express any vmRNA (viral load is zero); and (3) ''unexposed cells,'' cells derived from a control mouse that was never infected ( Figure 1A ). Utilizing a model that takes cross-contamination and various technical artifacts into account, we showed that infected cells can be recovered with high precision over a range of viral-load cutoffs (Supplemental Information and Figures S2A-S2C ). For example, among the inferred infected granulocytes we correctly predicted 99% by using a strict viral-load cutoff of 1%, and 88% by using a more permissive (0.01%) cutoff. Utilizing a reporter model of in vitro infection (Shnayder et al., 2018) , we were able to assess how . In each single cell, the host and the vmRNA were simultaneously measured, allowing identification of infected as opposed to bystander cells as well as retrospective annotation of cell types based on transcriptional identities. (B) Transcriptional signatures of single cells associated with nine major cell types. Gene expression of annotated cell-type-specific genes (rows) is shown across 4,064 single cells (columns). The expression matrix displays the clustering of cells into nine cell-type groups (first row), where each cell type consists of cells derived from different treatments (light/dark gray respectively denote control/influenza-treated animals; second row) and shows heterogeneity of cell infection states (white/yellow/brown respectively denote unexposed/bystander/infected cells, third row). (CÀE) Visualization of single cells across different cell types, treatments, and intracellular infection states. t-SNE plots show the cell-type annotation of single cells (C); the mouse to which the cell belongs (D), and the intracellular infection state (E). Color-coding is as in (B). See also Figure S1 and Table S1 . well vmRNA detection recapitulates the actual infected cells (see details in Supplemental Information). We observed that our cell grouping correctly predicts a high fraction of bona fide infected cells ( Figure S2D ): using viral-load cutoffs of 0.05%, we correctly recovered 79% of the infected cells. Infected cells that were not identified were those with low viral loads that lead to dropout events; it should, however, be possible to identify many of those cells by increasing the RNA-seq coverage.

We first explored the prevalence of infected cells within each of the nine major cell types. We focused on infected cells carrying a viral load of at least 0.05%, as this load faithfully represents a reasonable precision-recall trade-off (Supplemental Experimental Procedures, Figure S2E , and Table S2 ), high purity of influenza-infected cells (Russell et al., 2018, Supplemental Experimental Procedures and Figure S2F ), and relatively high viral expression level ( Figures S2G and S2H ). Based on this cutoff, we calculated for each cell type the percentage of infected cells among the cells derived from the influenza-treated mouse. Interestingly, relatively large percentages of infected cells were observed in each of the nine cell types, ranging between 22% and 31% in adaptive immune cells, between 25% and 36% in innate immune cells, and between 32% and 62% in non-immune cells, with the highest percentage of virus-infected cells (62%) observed in epithelial cells (Figures 2A and 1E ). Even after subtracting the expected percentages of erroneous predictions, intracellular viral infection was still prevalent in all cell types (Table S2) .

Several lines of evidence suggested that the high predicted prevalence of infection was not due to technical artifacts. First, we confirmed the high prevalence of infected cells using flow cytometry of immune cell subsets that were stained for the intracellular viral nucleoprotein (NP) (STAR Methods). In agreement with our predictions, we found a high prevalence of NP-positive cells within the infected lung tissue, with observed percentages of 48%, 43%, 55%, 61%, and 43% NP-positive cells in Cd45 + , Cd45 + Ly6g + , Cd45 + Cd11b + , Cd45 + Ly6g À Cd11b À , and Cd45 À cells, respectively (using false detection rate <5%; Figures (B) Extent to which the abundance of infected cells is an intrinsic property of the cell type. Shown is the percentage of infected cells in each cell type from an influenza-treated wild-type mouse (x axis) compared to its biological replicate (y axis). (C) Analysis of potential implications of the type I IFN pathway in cell-type-specific infection. Shown are the percentages of infected cells in an Irf7 wild-type (x axis) versus an Irf7-KO mouse (y axis). (D) Single-cell heterogeneity of intracellular viral load within the influenza-treated host. Shown are the percentages of low (yellow), medium (light brown), and high (dark brown) viral-load states (y axis) within the population of infected cells, as identified for each of the nine cell types (x axis; total numbers of infected cells are indicated). See also Figures S1-S4 , Tables S2 and S3. and S3B). In comparison, only 2% of NP-positive cells were observed when cells were stained for extracellular NP (Figure S3B, STAR Methods) , indicating that the observed high prevalence of infected (NP-positive) cells could not be attributed to free NP proteins or NP-RNA complexes. Second, we observe a low amount of background noise: only 25 of 2,075 unexposed cells (1.2%) were erroneously classified as infected cells (Figure 2A and Table S1 ). Third, an authentic influenza replication involves the generation of spliced isoforms for vmRNAs derived from specific segments. Indeed, we could detect the typical spliced form of the vmRNA derived from Segment 7 in cDNA libraries of Cd45 + cells ( Figure S3C , STAR Methods), strengthening the notion that viral replication is not limited to epithelial cells but can also occur in immune cells. Fourth, UMIs were specifically aligned to the 3 0 end of the polyadenylated viral RNA ( Figure S1D ) and cells were found to simultaneously express different viral genes ( Figure S1E ), confirming that the results were not attributable to non-specific predictions. Fifth, detection of vmRNA could, in principle, be confounded by the quality and complexity of cells. In that case, the infected and bystander cell populations should differ in their cell qualities. Such a difference is implausible, however, because the quality metrics of infected cells resemble those seen in bystanders ( Figures S3D and S3E) , with the single notable exception of cell size, which, in agreement with earlier studies (Xin et al., 2018) , is higher in the infected cells (Supplemental Information and Figures S3D-S3G ).

The widespread presence of infected cells found in all major cell types raises a fundamental question: do the differences in fractions of infected cells in different cell types ( Figure 2A ) reflect a true biological signal (rather than experimental noise)? To address this question we performed a second biological replicate of an influenza-treated animal (using the same pipeline as in Figure 1A ). This biological replicate consisted of the same cell-type populations, with totals of 833 bystanders and 402 infected cells (Table S1 ). Importantly, cell-type-specific percentages of infected cells were found to be highly reproducible in the two biological replicates (Spearman correlation = 0.93; Figure 2B), although there were differences in absolute fractions between experiments. In particular, T cells were associated with the lowest infection, while the highest percentage of infected cells was observed in epithelial cells. Among the rest, non-immune cells were generally associated with higher infection rates than those of immune cells. Notably, these results were reproducible in replicate mice examined 72 hr after infection (Spearman correlation = 0.83; Figure S4A ) and were consistent with previous reports ( Figure S4B ). Since all of the single-cell experimental and computational procedures were applied on each mouse independently, our results suggested that the observed differences in fractions of infected cells between cell types reflect a true cell-type-specific property.

One mechanism that might explain this cell-type-specific property is type I IFN signaling driven by the master transcription regulator protein Irf7 (Ciancanelli et al., 2015) . To test this possibility we generated single-cell RNA-seq data of lung tissue from an influenza-treated Irf7-KO mouse at 48 hr post infection (Table  S1 ). Intriguingly, despite the major difference in IFN signaling cir-cuits observed between Irf7-KO and WT mice (Figure S1A), we found good correlation between the cell-type fractions of infected cells in the two strains (Spearman correlation = 0.84; Figure 2C ). Together, these data suggested that the observed differences in infection percentages result from one or more cell-type-intrinsic properties, independently of the type I IFNmediated host response.

We next focused on the composition of viral loads in the subpopulation of infected cells. A recent report described a wide range of viral loads during influenza infection of epithelial cells in vitro (Russell et al., 2018) . We wanted to find out whether this epithelial cell property also applies during in vivo infection, and if so, we further aimed to assess the relevance of this finding to additional cell types. To address this, we divided the group of infected cells into three categories of viral-load states: low (<0.5%), medium (between 0.5% and 5%), and high (>5%).

We found that epithelial cells manifest a wide range of viralload states spanning two orders of magnitude, with 33%, 23%, and 44% of infected epithelial cells characterized respectively by low, medium, and high viral-load states ( Figures 2D and  S4C ). Our data therefore confirm the previously reported epithelial heterogeneity in vitro (Russell et al., 2018) and, in addition, extends this result to include the case of in vivo infection. Interestingly, we found a clear distinction between epithelial cells and the rest of the cell types, with the observation that infected cells, in all of the non-epithelial cell types, mainly displayed low heterogeneity of viral-load states: the vast majority (85%-100% of the infected cells in the different cell types) maintained low viral-load states, whereas high viral-load states were maintained by only up to 1.1% of those cells ( Figure 2D and Table S3 ). Similar results were obtained in additional infected mice ( Figure S4D ), further suggesting that viral-load heterogeneity is higher in epithelial cells than in other cell types, at least during early stages of infection.

Characterizing a Core Signature of the Antiviral Response Shared Across All Cell Types Evidence from diverse studies points to a powerful cellular response to influenza infection, involving the induction of multiple IFN-responding genes (Ivashkiv and Donlin, 2013) . However, current studies have typically used bulk-tissue data or focused on one or a few cell types (Altboum et al., 2014; Mostafavi et al., 2016) . Those studies were therefore limited in their ability to identify the ''generic'' influenza response that is presumably essential in all immune-and non-immune cell types. To characterize this response systematically, we calculated the differential expression between cell populations of the influenza-treated and control tissues for each of the nine cell types separately. We found a total of 450 genes to be differentially expressed (P binomial <10 À6 ) in at least one cell type ( Figure 3A , Table S4 ; STAR Methods). Clustering of these genes revealed a large generic module of 101 genes that showed consistent upregulation during influenza infection in nearly all immune-and non-immune cell types ( Figure 3A , e.g., Irf7 and B2m in Figure 3B ). As shown in Figure S5A , comparison of expression levels in the influenza-treated and control mice revealed that the module (legend continued on next page) was indeed expressed at higher levels during infection in all genes and cell types. As an additional support, the module was validated in an independent biological replicate experiment ( Figure S5B) . Notably, the majority of genes in this module represented known IFN-stimulated genes (e.g., P hg <10 À39 and P hg <10 À51 for the ''IFNb1'' and ''IFNR'' signaling categories, respectively) and targets of various antiviral transcription factors such as Irf7 and Stat1 (P hg <10 À41 and P hg <10 À46 , respectively), indicating that the module is part of a general type I IFN response that is independent of the cell type. For brevity, we refer to this generic type I IFN module as the ''IFN module.'' Having defined the IFN module, we were able to characterize the antiviral response in each single cell. We did this by generating an antiviral progression trajectory for each cell type by arranging the order of the individual cells according to their ''antiviral score,'' calculated as the average expression of a single cell across the IFN module. The trajectory revealed strong co-regulation in all cell types: genes of the IFN module were either co-expressed or anti-correlated through the trajectory ( Figure 3C ), as if the trajectory manifests a gradual progression of cells from a homeostatic to an antiviral state. Consistently with this notion, the influenza-treated cells were more strongly biased than the unexposed cells toward the high antiviral range (p < 10 À13 in all cell types; Figure S5C ). To validate the antiviral trajectory, we used the single-cell RNA-seq data from influenza-treated Irf7-KO mice. We hypothesized that in the absence of Irf7, the cells would exhibit an impaired trajectory of the antiviral response. We found that cells from the Irf7-KO sample indeed showed a clear bias toward the high antiviral range, but did not develop to their full antiviral potential ( Figure S5C ). This showed that the trajectory indeed represents a snapshot of the cells' progression through the antiviral response. We note that the plasticity of this antiviral trajectory probably reflects reversible regulatory alterations (unlike irreversible developmental processes), as previously observed in immune cells (Lavin et al., 2014; Okabe and Medzhitov, 2014) .

Recent studies demonstrated substantial suppression of mitochondrial-encoded genes in response to IFN and lipopolysaccharide treatments, which was attributed to either a metabolic shift in oxidative phosphorylation or removal of the entire mitochondrial fraction through mitophagy (Jovanovic et al., 2015; Kissig et al., 2017) . However, the generality of these findings in different cell types and their relevance for influenza stimulation remained unknown. Interestingly, analysis across the cell types revealed seven mitochondrially encoded genes that were consistently repressed in nearly all cell types ( Figures  3A and S5A ; the ''mitochondrial module''). Two lines of evidence further supported this module. First, we analyzed bulk data from influenza-treated lungs and found substantial repression of the mitochondrial module ( Figure S5D ). Secondly, the data clearly demonstrated anti-correlation of the IFN and mito-chondrial modules through the antiviral cell trajectory (Figure 3D ). Taken together, our data suggested that the generic IFN and the mitochondrial modules constitute two arms of the host response that are independent of the particular cell type.

Further inspection of three selected cell types (granulocytes, T cells, and fibroblasts, which had the largest amount of profiles among the innate, adaptive, and non-immune cell categories, respectively) indicated that about 70% of the expressed genes were dynamically up-or downregulated through the antiviral response trajectory (Figures S6A and S6B) . Interestingly, trajectory-regulated genes of different cell types were enriched with the same underlying regulatory programs ( Figures S6C and  S6D ), although the specific target genes of each regulatory program differed drastically between cell types (e.g., Daxx was regulated in granulocytes and T cells but not in fibroblasts, whereas Gbp6 was regulated in T cells and fibroblasts but not in granulocytes; Figures S6E and S6F ). Overall our data provide comprehensive information about which gene is regulated over the antiviral trajectory in particular cell types (Table S5 ). Since the cell-type specificity of the antiviral dynamics is currently known for only a limited number of genes, this mapping is critical for understanding the mechanisms of antiviral response progression in different cell types.

A high-resolution map of expression changes in infected and bystander cells can convey valuable information about the host-response circuits. In particular, differences in expression levels can be affected by extracellular exposure to the infectious lung environment, or may be associated with intracellular viral invasion. Previous studies have characterized many of the signaling cascades that operate in influenza infection, such as IFN and NF-kB signaling (Shapira et al., 2009 ), but the contribution of each stimulus to the host response was not well understood. We hypothesized that differential expression between the populations of bystanders and unexposed cells would identify transcriptional regulation associated with exposure to extracellular environmental signals, since both cell subsets lack intracellular host-pathogen encounters. Similarly, differential expression between infected and bystander cell populations could identify transcriptional regulation associated with the intracellular virus infection, since both subsets are essentially exposed to the same extracellular signals. For simplicity, we termed bystander/unexposed differential expression as ''bystander response'' and infected/bystander differential expression as ''infection response'' (Figure 4A, top; in both cases, a signed log binomial test p value was used). While it would have been interesting to test the effect across the entire spectrum of viral-load state ( Figure 2D ), here we focused on (B) Shown are the distributions of single-cell expression levels (y axis) of selected generic-response genes (subpanels) in each cell type (x axis) for cells derived from control (light gray) and influenza-treated mice (dark gray). Asterisks denote the significance of differential expression. (C) Antiviral cell trajectory. Shown is the expression of the IFN module of genes (rows, Z-normalized for each gene) across single cells ranked by their position along the antiviral trajectory (columns). Cells of different types appear in separate panels. (D) Opposite dynamics of the generic IFN and mitochondrial modules. Expression of the modules (y axis) in each cell type (color-coded) across the antiviral trajectory (x axis) is shown. See also Figures S5 and S6 , and Tables S4 and S5. the low/medium range that was represented in our data across all cell types (Tables S2 and S3 ).

The generic IFN module showed significant bystander response in all cell types, indicating a generic extracellular expo-sure effect. The infection response of this module, in contrast, was substantially smaller and was not consistently observed or reproduced across cell types ( Figures 4A-4C, top, and S7A) . It is well recognized that environmental signals during in vivo A C D B

Bystander and infection types of responses are illustrated (top), and are shown as signed log 10 p value (positive/negative, induced/repressed): the ''bystander response'' refers to the bystander versus unexposed differential expression and ''infection response'' refers to the infected versus bystander differential expression.

Shown is the average expression across the IFN module (top) and the mitochondrial module (bottom) for representative cell types. Asterisks denote significant bystander response or infection response (*p < 10 À3 , **p < 10 À13 ). infection act through the IFN feedback loop (Iwasaki and Pillai, 2014) . In agreement, our analysis of the infected Irf7-KO mouse revealed a significant reduction in the bystander response of the IFN module ( Figures 4C and 4D, left, and S7A) . The IFN module, therefore, is primarily regulated by extracellular exposure in a manner dependent on Irf7 activity.

The mitochondrial module demonstrated the effects of both extracellular exposure and intracellular infection by exhibiting a significant infection response in addition to a significant bystander response (Figures 4A and 4B ). In all cell types except for epithelial cells, the mitochondrial module was typically downregulated in bystander cells relative to unexposed cells, and the repression was further enhanced in the infected cell population (e.g., Figures 4B and S7BÀS7D) . The knockout effect of Irf7 ablated the bystander response of the mitochondrial module (as in the case of the IFN module), but had no substantial effect on its infection response ( Figures 4C and 4D , middle). Conceptually, this suggested that although the virus is maintained at a relatively low viral load, it is still sensed through Irf7-independent intracellular viral sensors. Furthermore, on the basis of asynchronous regulation, our observations clearly pointed to uncoupling between the mitochondrial and the IFN programs: whereas the peak of the IFN response was typically induced in the absence of intracellular invasion, the mitochondrial repression peak appeared in the presence of an intracellular influenza virus.

While the bystanders of most cell types exhibited considerable downregulation of mitochondrially encoded genes, the bystanders of epithelial cells displayed the opposite trend ( Figures  4D right, and S7C) . Notably, intracellular infection reversed the epithelial-cell-specific bystander response (with consistent results along the entire range of low/medium viral loads; Figure S7E) , and thereby all cell types converged to similar expression levels ( Figures 4D and S7F ).

Whereas our observations pointed to two distinct bystander mechanisms acting on the mitochondrial module (induction in epithelial cells and repression in other cell types), comparison with the Irf7-KO findings provided indications that different cell types differ only in their balance between the two underlying programs: (1) In epithelial cells, overexpression of the mitochondrial module in wild-type bystander cells was maintained (and even slightly elevated) in Irf7-KO bystander cells (Figures 4D, right, and S7B) , suggesting that the regulation in epithelial cells is achieved by strong Irf7-independent activation and weak Irf7-dependent inhibition. (2) In the remaining cell types, the bystander response in the wild-type was associated with repression while the bystander response in Irf7-KO showed the opposite effect (Figures 4C, left, and 4D , middle), consistently with combined control of a strong Irf7-mediated inhibition and a weak Irf7-independent activation. Similar results were obtained with an independent biological replicate ( Figures S8A and S8B) .

In summary, our data point to several generic mechanisms that act in all cell types. The IFN module responds predominantly to extracellular exposure, mediated mainly by Irf7 in all cell types, in agreement with previous findings (Iwasaki and Pillai, 2014) . The mitochondrial module, in contrast, is regulated by three distinct generic mechanisms: (1) bystander exposure leading to an Irf7-mediated inhibition (strong inhibition in all cell types except for epithelial cells); (2) bystander exposure that leads to an Irf7-independent induction of genes (strongest in epithelial cells); and (3) a generic repression, observed in all cell types in response to intracellular viral invasion.

In view of the generic response associated with intracellular infection, we next examined whether some programs display cell-type-specific effects. As statistical power is reduced upon testing a single gene in a single cell type, we analyzed three cell types with the largest amount of cell profiles-granulocytes, T cells, and fibroblasts-as representatives of innate, adaptive, and non-immune cells. A total of 234 genes attained significant differential expression between infected cells and bystanders in one of these cell types (false discovery rate [FDR] < 0.01, Table  S6 ), suggesting a transcriptional response associated with intracellular viral infection. To validate the identified genes, we carried out an experiment using an independent biological replicate. The results showed that the candidate genes significantly overlap with those identified by the replicate mouse. In granulocytes, for instance, of 47 candidates, 17 (36%) attained significant scores in another mouse (P hg <10 À7 , Figure S8C ). Differential expression scores were also strongly correlated between the two replicates ( Figure S8D ). We therefore focused on the set of 34 validated genes ( Figure 5A ). Examination of the expression pattern along the antiviral trajectory showed that the identified targets consisted primarily of responding genes, with stronger response amplitudes in infected cells than in bystanders ( Figures  5A and 5B ). In particular, in 24 of the 34 genes (70%), the bystander and infection responses followed the same direction (either enhanced repression or enhanced induction). Thus, intracellular infection was generally associated with an elevated peak of response of cell-type-specific genes, irrespective of the direction of response. This is not unexpected, given that the intracellular viral sensing through Toll-like receptor 3 and RIG-I probably resulted in an enhanced antiviral response.

Overall, the 34 targets showed enrichment in antiviral and inflammatory response pathways (p < 10 À20 ) and in functional categories related to cell movement, migration, and homing genes (p < 10 À11 ). In addition, we found that the gene targets were enriched in extracellular proteins (either secreted or attached to the external side of the plasma membrane; p < 0.03). In conclusion, our results point to cell-type-specific functional influences of intracellular viral infection on the composition and functionality of both intracellular and extracellular environments.

In this study, we used single-cell RNA sequencing of polyadenylated mRNA to simultaneously study the viral and host transcriptomes in the same individual cell. The biology of the influenza virus offers valuable insights into the interpretation of these data. Assembled influenza virions possess eight segments of a non-polyadenylated negative-stranded RNA genome; inside the infected cell, the viral genome is replicated and also serves as a template for transcription of positive-sense polyadenylated vmRNAs (York and Fodor, 2013) . Taken together, the polyadenylation signature provides specific detection of vmRNAs that were produced during intracellular viral transcription. Importantly, our strategy is capable of identifying intracellular viral products even in the absence of productive infection (Brooke et al., 2013; Heldt et al., 2015) . This property is particularly important to reconstruct the host response to any type of intracellular viral constituents.

Our dual viral-host scRNA-seq mapping should be broadly applicable to a wide range of viruses. In particular, our mapping is tailored for viruses that generate polyadenylated mRNAs during their intracellular replication. scRNA-seq is especially effective for viruses in which the RNA genome is embedded within free viral particles that are not polydenylated. Our approach should be therefore applicable to a variety of medically important negative-sense single-stranded RNA (ssRNA) viruses (such as RSV, influenza, Ebola, and measles) and double-stranded DNA viruses (such as herpesviruses, adenoviruses, and poxviruses). Our dual-transcriptome mapping is less effective in reverse transcribing viruses and doublestranded RNA viruses, which have two sources of positive-stranded RNAs (both genome-and transcriptome-derived RNAs); further optimization will be needed to apply our methodology in these cases. Of note, in the special case of positivesense ssRNA viruses, our approach is not applicable due to the lack of polyadenylated transcriptome (e.g., Zika, dengue, and yellow fever viruses) or the polyadenylated genome (e.g., foot-and-mouth disease virus, rubella, and severe acute respiratory syndrome). Tailored technologies have been developed for this specific case (e.g., Zanini et al., 2018) . Our dual-transcriptome approach shares with these technologies the challenge of distinguishing intracellular viral transcription from exogenous acquisition (by phagocytosis) of infected cells. Combining single-cell transcriptomics data with flow cytometry or mass cytometry (CyTOF) measurements could help determine the cell's infection state.

Infectivity of the influenza virus has traditionally been ascribed to epithelial cells of the airways, although low amounts of A B Figure 5 . Cell-Type-Specific Transcriptional Changes after Intracellular Infection (A) List of the top 34 genes manifesting significant infection response (false discovery rate, p < 0.01) in a particular cell type. The bystander response (bystanders versus unexposed) and infection response (infected versus bystanders) of each gene in its relevant cell type are shown on a red/blue color scale (signed log 10 p value ranging from À50 to +50; positive/negative, induced/repressed). The three bottom rows indicate membership in functional categories. Directionality of effects; enhanced induction, enhanced repression, and other dynamics are marked. (B) Representative genes displaying significant infection response in a particular cell type. Shown are the moving averages of expression levels in bystander (yellow) and infected (brown) cells (y axis) along the antiviral trajectory (x axis) for representative gene-cell type pairs from A. See also Figure S8 and Table S6. infection have also been reported in other cell types (De Baets et al., 2015; Manicassamy et al., 2010; McFadden et al., 2009 ). Our findings suggested that viral infection in vivo is highly prevalent in many additional cell types (Figure 2A ). In agreement with the transcriptome analyses, our flow cytometry assays confirmed that immune cells indeed harbor a large percentage of cells expressing the viral NP during in vivo infection ( Figures  S3A and S3B ). If the high prevalence of infected cells is not due to technical noise, this model would predict that the percentage of infected cells would be a cell-type-specific property. Quantitatively similar fractions of infected cells were indeed observed across the different biological replicates examined ( Figures 2B, 2C, S4A, and S4B ). The prediction that infected cells would be prevalent in all cell types was consistent with the observation that much of the wide spectrum of antiviral-related cell states was generic and appeared in nearly all cell types (Figures 3 and 4 and Supplemental Experimental Procedures). Thus, instead of a different antiviral response being tailored to ''susceptible'' and ''resistant'' cell types, the host response confers the first line of defense in all key cell types.

In agreement with our findings, infectivity levels based on a GFP-expressing virus have been previously reported across all cell types (De Baets et al., 2015; Manicassamy et al., 2010; McFadden et al., 2009) . Notably, those studies could identify only a low fraction of infected cells in each cell type, probably owing to a combination of several differences: GFP detection is likely to be less sensitive than detection of amplified vmRNA, and GFP-expressing viruses rapidly lose their reported gene (see, e.g., Manicassamy et al., 2010) ; Furthermore, the expression levels of different viral genome segments may vary significantly within cells (Heldt et al., 2015) , implying that some cells may display low GFP levels despite their substantial expression of other viral genes. It should be borne in mind that our study was limited to a particular experimental setting, and further experimentation will therefore be needed to investigate the effects of various host genetic backgrounds and viral strains.

Recent studies have demonstrated a wide range of intracellular virus-related states, including heterogeneity in viral transcription and virus titers. For example, in vitro infection of epithelial cells results in wide heterogeneity of transcriptional viral loads (Russell et al., 2018) , and during in vitro infection, low viral loads are prevalent in Madin-Darby canine kidney epithelial (MDCK) cells (Brooke et al., 2013) . Here, focusing on in vivo infection, we found that epithelial cells indeed display high viral-load heterogeneity, whereas the remaining cell types maintain low viral-load states ( Figure 2D) . Notably, high vmRNA expression does not necessarily imply high viral progeny, probably owing to failure to express one or more viral segments (Heldt et al., 2015) ; in fact, up to 90% of cells infected at low multiplicity of infection fail to release infectious progeny (Brooke et al., 2013) . Changing intracellular host response and viral loads could affect viral titers and should be further explored (Cristinelli and Ciuffi, 2018) . Novel experimental techniques will be needed to reveal the underlying mechanisms of virus production levels through joint analysis of host transcriptional responses ( Figure 1A ) and viral titer assays (Heldt et al., 2015) at single-cell resolution. Revealing such mechanisms will be of great interest, since their modification may lead to more effective therapeutic interventions.

Epithelial cells are considered to be a major source of influenza virion progeny (Manicassamy et al., 2010; McFadden et al., 2009) , but the contribution of intracellular viral invasion relative to that of the intracellular viral-load state has not been investigated in detail. Our data suggested that the intracellular viral state is a predominant player in determining the higher viral progeny of epithelial cells compared with any other cell types: the fraction of high viral-load cells was found here to be 30-fold larger in epithelium ( Figure 2D ), while the fraction of infected epithelial cells was only $1.9-fold larger (Figure 2A) .

Studies of several cell types have demonstrated that the type I IFN response is accompanied by an antagonistic suppression of mitochondrially encoded genes (Jovanovic et al., 2015; Kissig et al., 2017) . However, those studies relied on data from a single cell type, ignoring the possible generality of the findings. By analyzing influenza infection over all cell types from the entire lung, we were able to establish the general principles of antiviral induction and mitochondrial suppression across all of the abundant cell types of influenza-treated lungs. The existence of two antagonistic generic programs may be an efficient antiviral strategy: whereas the activation of the type I IFN response provides antiviral defense mechanisms such as biosynthesis and secretion of cytokines, suppression of mitochondrial activity typically limits energetic resources, which, in turn attenuates the efficiency of viral replication and manipulation (Goswami et al., 2013) .

Suppression of mitochondrial energetic resources during in vivo influenza infection may be both beneficial (limiting intracellular viral replication) and detrimental (limiting rapid host response). It is therefore likely that the different cell types maintain an appropriate balance between antiviral response and mitochondrial suppression. Our data suggested (1) that the heterogeneity of the generic IFN module typically reaches its main potential in bystander cells (Figure 4) , and upon intracellular viral invasion the antiviral program is fine-tuned in a cell-typespecific manner ( Figure 5 ); (2) that the peak of mitochondrial generic suppression typically arises in infected cells, and that in the absence of intracellular viral invasion, bystander cells demonstrate cell-type-specific control of mitochondrially encoded genes (upregulation in epithelial cells and downregulation in the rest of the cell types; Figure 4 ). Given these observations, it is tempting to speculate that the mechanism of mitochondrial energy balance is addressed by means of temporal uncoupling of the IFN module from the mitochondrial module: bystanders exploit the ''window of opportunity'' prior to intracellular infection to achieve their shift toward high IFNresponse states, while supporting the energetic demands with an adequate amount of energetic resources. Mitochondrial suppression reaches its full capacity only in infected cells, when viral restriction is crucial in a ''just-in-time'' mechanism (Zaslaver et al., 2004) . Accordingly, after intracellular viral invasion, when energetic resources are limited, substantial adjustment of the host gene expression is limited to specific genes and cell types.

Although epithelial cells have been studied in depth, it is still unclear how their response program differs from that of other cell types. Such comparison has been difficult because of the need to isolate multiple cell types simultaneously from the same infected organ. This question is particularly relevant because exposure of the epithelial barrier to the invading pathogen probably precedes the exposure of other cell types; thus, the ability of epithelial cells to maintain substantial secretion of antiviral and pro-inflammatory cytokines may be crucial for an efficient and timely host response. The present study shows that epithelial cells are phenotypically distinct: bystander epithelial cells display substantial induction of mitochondrially encoded genes, rather than the observed characteristic partial suppression by the rest of the cell types studied here ( Figures 4C and 4D) . Thus, our data support a model in which bystander epithelial cells display tailored changes in energetic resources to support the feasibility of their excessive energetic demands. While these observations are consistent with the experiments described here, additional experimentation will be needed to determine whether the mitochondrial transcriptome is decoupled from the cell size during viral infection in vivo.

An important therapeutic challenge is to identify particular cell functionalities that can distinguish between infected and bystander cells. Our data indicated that mitochondrial functionality is an attractive target, given the observation that in all cell types the mitochondrial module was more strongly suppressed in infected cells than in bystanders ( Figure 4A ). Another promising observation was that mobility-and extracellular matrixrelated genes are differentially expressed in infected and bystander cells ( Figure 5) . Furthermore, infected cells tended to be larger than bystanders ( Figure S3G ). Determining whether such functional differences exist and whether there might be interrelationships among the various functionalities will require further experimentation. In such experiments, both viral load and functional characterizations will need to be measured simultaneously at the resolution of single cells. The ability to relate cell functionalities with viral-load states should prove useful in the future development of infected-cell-specific therapeutic strategies.

Detailed methods are provided in the online version of this paper and include the following: sample was then linearly amplified by T7 in-vitro transcription, and the resulting RNA was fragmented and converted into a sequenceready library by tagging the samples with pool barcodes and illumina adapter sequences during ligation, RT, and PCR. Each pool of cells was tested for library quality, and concentration was assessed as described (Jaitin et al., 2014; Paul et al., 2015) .

mRNA-Seq Pre-processing The analysis starts with a joint alignment of MARS-seq reads for both the host and viral genomes, followed by the joint quantification of host and viral transcripts. We performed this joint analysis as previously described (Shnayder et al., 2018) , with a few modifications. We used HISAT (Kim et al., 2015) to map the single-cell reads to a joint reference of both the mouse genome build mm9 and the influenza genome (NCBI records: NC_002016, NC_002017, NC_002018, NC_002019, NC_002020, NC_002021, NC_002022, NC_002023). To identify authentic intracellular viral replication (rather than free virions), specificity for the positive-sense viral mRNA is imperative; in accordance, we mapped reads to the positive-strand-specific isoforms of the viral genes. Reads that were mapped at multiple (host or viral) positions were omitted. We recorded reads that were mapped to either influenza segments or to exonic regions of mouse genes based on the UCSC transcript annotation. The expression level for each single cell was calculated on the basis of its number of UMIs. To account for spurious mapping, we filtered reads with cell-barcoding errors and UMIs that contained more than one error that might be due to sequencing or amplification errors. Quantification of technical noise associated with the assessment of viral transcripts is described in Supplemental Experimental Procedures. Poor-quality genes (with <5 recorded UMIs across all cells) were excluded from further analysis. We also filtered poor-quality cells (with < 200 recorded UMIs across all genes; Table S1 ). Additional cells were filtered at the cell-type annotation stage (see below).

The 'viral load' of a cell was defined as the number of UMIs that map to one of the eight viral segments, expressed as a percentage of all mapped UMIs. A cell was considered 'infected' if its viral load was higher than a certain cutoff. 'Bystanders' were defined as cells derived from an influenza-treated mouse in which no viral segments were detected. These definitions were previously used in a combined single-cell mapping of host and viral transcriptome (Russell et al., 2018) . To determine the quality of this categorization we used precision and recall (sensitivity) measures, as well as the percentage of false negatives (Supplemental Experimental Procedures).

Cell types were annotated in four steps. First, we identified and filtered host genes that were highly variable across cells and thus might blur the cells' identity. We did this by performing principal component analysis (PCA) on all cells, and omitting from further analysis genes that correlated with the top two PCs (top 1.5% of PC-correlated genes). In the second step, cells were clustered using the multinomial mixture model clustering algorithm (Paul et al., 2015) . The first two steps were performed separately on CD45 + and CD45 À cells. In the third step, each of the resulting cell groups was annotated on the basis of two criteria: (1) by the presence of a set of well-known cell surface markers, as demonstrated in Figure 1B , and (2) by resemblance to profiled bulk populations of cell types from the ImmGen compendium (Heng et al., 2008) . For each group, the resemblance was estimated using deconvolution on the group's average profile as the dependent variable and bulk profiles as explanatory variables. In the fourth step, groups were merged into 12 clusters based on their cell-type annotation. The nine largest clusters that appeared in all mice were used for further analysis, and the remaining three small clusters were discarded (Table S1) . As an example, we analyzed a total of 2075 cells (from control mice) and 1989 cells (from infected mice, 48 h, replicate I) within these nine clusters (Table S1 ). The quality of the clustering solution was evaluated according to the homogeneity of each cluster at single-cell resolution. In particular, we tested the similarity between each single cell and a bulk cell-type profile (using Pearson's correlation; Figure S1B ). For visualization purposes, the data were reduced to two dimensions using the t-SNE algorithm (Van Der Maaten and Hinton, 2008) , based on the known cell-type markers from Figure 1B .

Detection of the Influenza NP Lung samples from C57BL/6J-infected or PBS-treated mice were prepared for flow-cytometric analysis as described above. Cd45 + cells were enriched by the use of Cd45 microbeads (130-052-301 MACS, Miltenyi Biotec). Prior to fixation, cells were stained with the following antibodies (all from Biolegend): anti-mouse I-A/I-E (BLG-107621, 1:50), anti-mouse Ly-6g (BLG-127607, 1:100), antimouse/human Cd11b (BLG-101215, 1:100) and anti-mouse Cd64 (BLG-139305, 1:50). For intracellular NP staining, cells were permeabilized using an intranuclear staining kit (00-5523, Invitrogen) and were then incubated with NP for 1 h. No permeabilization was performed in the case of extracellular NP staining. Gates were set at standard strict FSC and SSC levels. After exclusion of doublets, NP levels were monitored in populations of Cd45 + , Cd45 + Ly6g + , Cd45 + Ly6g -Cd11b + , Cd45 + Ly6g -Cd11b -, and Cd45cells.

Viral spliced forms could not be identified based on our 3'-end single-cell mRNA libraries (Jaitin et al., 2014; Paul et al., 2015) since they share their 3' sequence. To identify viral spliced forms within immune cells, Cd45 + cells were taken from lung samples of infected and control mice as described above. Total RNA was then extracted from these cells and reverse-transcribed using qPCRBIO cDNA Synthesis Kit (PCR Biosystems Cat# PB30.11-10). The generated cDNA was used as a template in PCR reaction, using GoTaq Flexi DNA Polymerase (M8295, Promega) with the following combinations of primers that are specific for the spliced or unspliced forms of

Digital cell quantification identifies global immune cell dynamics during influenza infection

Pathogen cell-to-cell variability drives heterogeneity in host immune responses

A systems analysis identifies a feedforward inflammatory circuit leading to lethal influenza infection

Replication-competent influenza A viruses expressing reporter genes

Most influenza A virions fail to express at least one essential viral protein

Infectious disease. Life-threatening influenza and impaired interferon amplification in human IRF7 deficiency

The use of single-cell RNA-Seq to understand virus-host interactions

A GFP expressing influenza A virus to report in vivo tropism and protection by a matrix protein 2 ectodomain-specific monoclonal antibody

Single-cell analysis and stochastic modelling unveil large cell-to-cell variability in influenza A virus infection

The Immunological Genome Project: networks of gene expression in immune cells

Classification of low quality cells from singlecell RNA-seq data

Regulation of type I interferon responses

Innate immunity to influenza virus infection

Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types

Dynamic profiling of the protein life cycle in response to pathogens

HISAT: a fast spliced aligner with low memory requirements

PRDM16 represses the type I interferon response in adipocytes to promote mitochondrial and thermogenic programing

Tissue-resident macrophage enhancer landscapes are shaped by the local microenvironment

Analysis of in vivo dynamics of influenza virus infection in mice using a GFP reporter virus

Cytokine determinants of viral tropism

Parsing the interferon transcriptional network and its disease associations

Tissue-specific signals control reversible program of localization and functional polarization of macrophages

Transcriptional heterogeneity and lineage commitment in myeloid progenitors

Extreme heterogeneity of influenza virus infection in single cells

A physical and regulatory map of host-influenza interactions reveals pathways in H1N1 infection

Defining the transcriptional landscape during cytomegalovirus latency with single-cell RNA sequencing

Reconstructing lineage hierarchies of the distal lung epithelium using single-cell RNA-seq

Visualizing data using t-SNE

Single-cell analysis of the impact of host-cell heterogeneity on infection with foot and mouth disease virus

Biogenesis, assembly, and export of viral messenger ribonucleoproteins in the influenza A virus infected cell

Singlecell transcriptional dynamics of flavivirus infection

Just-in-time transcription program in metabolic pathways

The authors declare no competing interests. 

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact Irit Gat-Viks (iritgv@post.tau.ac.il).

Adult female C57BL/6J mice aged 5.5À6.5 weeks were used in the study. Wild-type animals were obtained from Harlan and Irf7-KO mice were purchased from RIKEN BioResource Research Center. All mice were housed on hardwood chip bedding under specificpathogen-free conditions at the Animal Breeding Center of the Weizmann Institute of Science (Rehovot, Israel). They were held in individually ventilated cages and lights were on a 12-h light/dark cycle. Rodent chow and water were given ad libitum throughout the experiment. All animals were handled according to the regulations formulated and approved by the Institutional Animal Care and Use Committee (IACUC).

Mouse-adapted PR8 H1N1 influenza virus (A/Puerto Rico/8/34) was grown in hen egg amnion. Anesthetized mice were inoculated intranasally with 40 mL of diluted virus (6310 3 PFU/mouse). Control mice received 40 mL of sterile PBS intranasally. Body weight was measured daily. All of the infected mice showed 8% to 11% body weight loss 48 hr after infection. Mice were killed 48 hr or 72 hr after treatment (see Table S1 for a detailed list of experiments).

Infected or PBS-treated mice were killed by anesthesia overdose. Their lungs were then perfused with PBS via the right ventricle, and then divided into two groups that were subjected to different dissociation techniques. Half of the entire pool of lungs was dissected and dissociated into single-cell suspensions, using the Miltenyi Lung Dissociation Kit (Miltenyi Biotec) in combination with a gentleMACS dissociator (Miltenyi Biotec) and enzymatic dissociation. The other half of the lung pool was minced, suspended in 5mL of digestion buffer consisting of elastase (3 U/mL; Worthington Biochemical Corporation) and DNase I (0.33 U/mL; Sigma-Aldrich) in DMEM/F12 medium, incubated with frequent agitation at 37 C for 20 min, and briefly triturated (Treutlein et al., 2014) . Whereas the Milteny kit is optimized for most immune cell types, the elastase-based protocol allows for a higher representation of epithelial cells. Next, an equal volume mixture of DMEM/F12 supplemented with 10% fetal bovine serum and penicillin-streptomycin (1 U/mL; Biological Industries) was added to the single-cell suspensions in both groups. Following their enzymatic incubations, cells derived from the same lungs (after the two dissociation techniques) were merged and forced through a 100-mm mesh and then through a 70-mm mesh into ice-cold FACS buffer (2 mM EDTA pH 8.0, 0.5% bovine serum albumin in PBS). Prior to sorting, all samples were stained for DAPI, filtered through a 40-mm mesh, and stained with F450-conjugated TER-119 and APC-conjugated CD45 antibodies (BioLegend). Cell populations were sorted using SORP-aria (BD Biosciences) and gates were set at standard FSC and SSC levels. Following erythrocyte and doublet exclusion (by application of an FSC-H/SSC-W singlets gate), isolated DAPI -TER-119 -CD45 + and DAPI -TER-119 -CD45cells were single-cell-sorted into 384-well cell-capture plates containing 2 mL of lysis solution and barcoded poly(T) reverse transcription (RT) primers for single-cell RNA-seq (Jaitin et al., 2014; Paul et al., 2015) . Four empty wells served as controls in each 384-well plate. Immediately after sorting, each plate was spun down to ensure immersion of cells into the lysis solution, snap-frozen on dry ice, and stored at À80 C until processed.Massively Parallel Single-Cell mRNA-Seq Library Preparation (MARS-Seq) Single-cell mRNA libraries were prepared as previously described (Paul et al., 2015; Shnayder et al., 2018) . Briefly, mRNA from singlesorted cells was barcoded using poly(T) RT primers, converted into cDNA, and pooled using an automated pipeline. 

Differential expression between two cell populations, A and B (e.g., influenza-treated mouse versus control; unexposed cells versus bystander cells) was calculated for each gene independently. The differential expression test makes the null assumption that the two cell populations are derived from the same binomial distribution of gene expression: x i $ B(n i , p AB ), where x i is the number of UMIs of a certain gene in the i-th cell, n i is the total UMI count in the i-th cell, and p AB is the probability of sampling a UMI associated with that gene in both populations A and B. The alternative hypothesis is that the expression in each cell population is derived from a different binomial distribution: x i $ B(n i , p A )and x i $ B(n i , p B ), respectively, in cells from populations A and B, where p A and p B are the probabilities of sampling a UMI associated with the relevant gene in the respective cell populations. Accordingly, here we calculated the maximum likelihood of the null and the alternative hypotheses by fitting the binomial parameters b p A , b p B and b p AB to each cell population separately. Finally, a likelihood ratio test (LRT) was applied and a P-value was approximated in accordance with Wilks' theorem using a chi-squared distribution with one degree of freedom. We report the 'differential expression' (DE) either as a P-value or as a signed log 10 (LRT P-value), where the sign corresponds to the difference between the maximum likelihood parameters (namely, the success probabilities b p A and b p B ) of the two cell populations. Note that calculation of the differential expression takes into consideration differences in total UMI counts (n i values) among the two populations of cells. Here no adjustment for inflation of observed P-values is needed ( Figure S7G ).In all cases, we applied differential expression tests only on genes that were associated with reported mapped UMIs in at least 20% of the cells derived from one of the two compared cell populations. In particular, the 'bystander response' is defined as the differential expression between bystander and unexposed cells. Similarly, the 'infection response' is defined as the differential expression between infected and bystander cells. As the power of infection response testing largely depends on the purity of the population of infected cells, all infection response tests relied on improved precision that is associated with compromised recall: the cells that are considered 'infected' exclude 1-viral-UMI cells and include cells carrying viral load > 0.03%. This rule indeed attained substantially higher precision ( Figure S2I , top) while maintaining a reasonable number of inferred infected cells ( Figure S2I, bottom) .

To identify the generic response, we compared cells derived from the influenza-treated mouse and the control mouse. We identified 450 nuclear genes that satisfied the following four criteria: (i) genes that were differentially expressed (P binomial <10 -6 ) in at least one cell type; (ii) genes in which the fold change in mean expression was greater than 2 in at least one cell type; (iii) genes that were expressed in more than 20% of either the control or the influenza-treated mice; (iv) genes in which the mean difference in expression between the two cell populations was higher than 0.0004 in at least one cell type. The resulting 450 responding genes ( Figure 3A and Table S4) were then subjected to hierarchical clustering, which highlighted a group of 101 genes that were consistently upregulated in all cell types. This gene group is referred to as the 'generic IFN module'.Given this module, the 'antiviral state' of a cell was defined as the average expression of the generic IFN module in each single cell. Ordering of cells by their state is referred to as the 'antiviral trajectory'. This trajectory was validated in cells of each cell type separately by demonstrating: (i) enrichment of antiviral-related functionalities in trajectory-correlated genes ( Figure S6D ); (ii) distinct distributions of Irf7-KO-derived and wild type-derived single cells through the trajectory ( Figure S5C ); and (iii) co-regulation of genes along the trajectory (Figures 3C and 3D) . Transcriptional dynamics is represented as a moving average of expression over the antiviral trajectory, sliding one cell at a time with a window size of one-sixth of the trajectory length ( Figures 3D, 5B , S6F, and S7D).

We limited our analysis to 3 key cell types: granulocytes, T cells and fibroblasts as the largest representatives of innate, adaptive and non-immune cells, respectively (Table S1 ). To characterize the dynamics of genes along the antiviral trajectory, we first selected genes that were expressed in at least 20% of the cells and then clustered the genes based on their expression in the infected tissue along the antiviral trajectory (independently in each cell type). To restrict our analysis to a robust signature of transcriptional dynamics, we assembled only the genes that reside in large clusters (>460 genes, Figure S6 and Table S5 ).

Enrichment of regulators and of pathways was calculated using a hyper-geometric test using the Ingenuity Knowledge Base (QIAGEN) gene annotation. All reported P-values were Bonferroni corrected.

The single cell data reported in this study have been deposited in the Gene Expression Omnibus (GEO) database, under GEO accession number GSE107947.