Single-Cell Profiling of Ebola Virus Disease In Vivo Reveals Viral and Host Dynamics

1Ebola virus (EBOV) is among the world’s most lethal pathogens, with estimated case fatality rates of 40%–66% in recent epidemics. EBOV infection in humans causes Ebola virus disease (EVD), characterized by fever, malaise, muscle aches, and gastrointestinal distress, rapidly progressing to coagulopathy, shock, and multi-organ failure. While recently developed vaccines and monoclonal antibody therapeutics have shown great promise, case fatality rates in treated patients still exceed 30%, highlighting the need for further research into EVD.

2Studies of EVD pathogenesis face numerous logistical challenges that have limited their scope. Experiments with live EBOV require maximum containment (biosafety level 4 [BSL-4]), restricting them to a few specialized research facilities. In vivo studies are especially challenging: human EVD is difficult to study during deadly outbreaks in resource-limited settings, necessitating animal models. Commonly used laboratory mice do not display key features of human EVD when exposed to naturally occurring EBOV. While other small animal models such as ferrets mostly recapitulate human EVD, they are genetically distant from primates and lack the primate-specific genotype for NPC1, the key cellular receptor for EBOV entry. EVD in nonhuman primates (NHPs) most closely resembles human EVD, but NHP studies are often limited to small sample sizes.

3The two main approaches to studying EVD—analyzing infected cells in culture and infected animals in vivo—have revealed important, if somewhat contradictory, aspects of EBOV’s impact on the immune system. In culture, EBOV infects myeloid cells, potently inhibiting both production of type 1 interferon and signal transduction downstream of interferon receptors. Under-activation of this key innate antiviral response hinders the activation of the adaptive immune system by antigen-presenting cells, a key determinant of fatal outcomes. EVD in vivo, by contrast, is characterized by high fever and dramatic upregulation of hundreds of interferon-stimulated genes (ISGs), in response to dozens of inflammatory cytokines, suggesting that an aberrant over-activation of innate and adaptive immunity underlies much of EVD pathology.

4High-throughput single-cell technologies, such as single-cell RNA-sequencing (scRNA-seq) and protein quantification by CyTOF, have enabled the analysis of viral infection at unprecedented resolution. These methods quantify the cell-type composition and expression programs of individual cells—signals that are obscured in bulk measurements. By quantifying viral RNA within cells, scRNA-seq allows comparison of gene expression between infected and uninfected (bystander) cells in a diseased host, giving a far more nuanced view of host and viral gene expression within infected cells. This approach can also disentangle direct effects of infection within a cell from the effects of the inflammatory cytokine milieu, which bulk approaches cannot easily do. However, to date, technical constraints have prevented use of these approaches in a BSL-4 setting: many scRNA-seq technologies require droplet generators and inactivation protocols not well suited to that environment, while CyTOF instruments involve high-volume exhaust and superheated components incompatible with BSL-4 installation. Thus, new protocols compatible with sample inactivation are needed.

5Here, we describe the first BSL-4 investigation of a pathogen with high-dimensional single-cell technologies. We applied CyTOF and Seq-Well—a portable single-cell RNA-seq platform—to a total of 90 peripheral blood mononuclear cell (PBMC) samples from 21 rhesus monkeys prior to, or during, lethal EBOV challenge in vivo. We further inoculated PBMCs with EBOV ex vivo and profiled their gene expression with Seq-Well. These data allow us to dissect host-virus interactions and comprehensively catalog changes in cell-type abundance and cell state over the course of EVD. Moreover, as EBOV has an RNA genome and transcribes polyadenylated mRNAs, we detected viral RNA within individual cells, allowing us to define EBOV tropism with high resolution and identify EBOV-associated transcriptional changes in putative pro- and antiviral genes.

6To comprehensively profile EBOV-induced immune dysfunction in vivo, we collected peripheral immune cells prior to, and at multiple days post-infection (DPI), corresponding to different stages of acute EVD (Figure 1). Cohorts of ≥3 NHPs were sacrificed as baseline uninfected controls, at pre-defined DPI, or upon reaching humane euthanasia criteria. These cohorts were recently characterized for viral load, clinical score, blood chemistry, and liver pathology. Viral load first became detectable in all NHPs on DPI 3, preceding detectable clinical signs (e.g., fever) by 1–2 days (Figure 2A). Clinical signs of EVD progressed until euthanasia criteria were reached at DPI 6–8. Cells collected at all time points were used for CyTOF, while cells collected at baseline and at sacrifice were also used for Seq-Well (Figure S1A).

7After standard quality control filters (STAR Methods), we obtained single-cell transcriptomes from ∼58,000 PBMCs and 42-protein CyTOF profiles from ∼15 million (M) PBMCs. We visualized these data with uniform manifold approximation and projection (UMAP) non-linear dimensionality reduction (Figures 2B–2H). Unsupervised clustering of either the transcriptomes or a down-sampled set of 1.1 M protein profiles (STAR Methods) yielded clusters that could be readily identified from well-known RNA and protein markers as the major circulating immune cell types (Figures 2B, 2C, S1B, and S1C) and agreed well with manual gating (Data S1). After batch correction of the CyTOF data and integration of the transcriptomes (STAR Methods), samples were well distributed across cell-type clusters (Figures S1D–S1G) but separated by DPI (Figures 2E and 2F), indicating dynamic cell states over the disease course. By sub-clustering within cell types, we identified subtypes based on identifying marker genes (Figure S2).

8In addition to the major PBMC cell types, a cluster of immature neutrophils emerged during EVD, marked by high gene expression of CD177 and SOD2, and protein expression of CD66 and CD11b. Based on scRNA-seq, neutrophils increased from 0.2% of cells at baseline to 65.1% in late EVD (by CyTOF, from 9.3% to 49.8%; Figures S3A and S3B). Though neutrophils are typically removed during density-based PBMC isolations, immature neutrophils (band cells) are less dense than mature neutrophils and can be released from the bone marrow and co-isolate with PBMCs in infections, including during EVD.

9We quantified absolute abundance of each cell type during EVD by combining CyTOF or scRNA-seq data with complete blood counts (STAR Methods); both estimates were in general agreement (Figure S3A).

10Neutrophil abundance increased from baseline by 5-fold by DPI 4, before returning to baseline levels in late EVD (p < 0.05 for DPI 3–4, p = 0.059 for DPI 5, Wilcoxon signed-rank test, Figures 2D and S3B). Lymphocyte abundance decreased, with natural killer (NK) cells declining 1 day before the other cell types (p < 0.05 on DPI 3–6 for B, NK, CD8+ T, and CD4+ T; except for CD4+ T on DPI 5); both observations are in agreement with previous NHP studies. All lymphocyte populations slowly recovered after DPI 4 (Figure S3B). Monocyte abundance initially increased >2-fold before declining precipitously between DPI 4 and 5.

11Changes in circulating cell-type abundance could reflect cell proliferation, death, and movement into and out of bone marrow, lymph, and tissues. While we could not directly quantify death or movement between compartments, we estimated the fraction of dividing cells using the proliferation marker Ki67 in both the CyTOF and scRNA-seq data (Figure S3C). The fraction of Ki67+ monocytes increased from 17% at baseline to 56% at DPI 5 and remained > 40% thereafter (p = 1.1 × 10–5, rank-sum test of DPI 5–8 versus baseline), suggesting increased proliferation (Figures 2G and S3D). The maximum on DPI 5 coincided with a decline in absolute abundance (Figures 2D and S3B), suggesting that losses due to extravasation or death dominate increases from proliferation in late EVD. By contrast, neutrophil proliferation remained roughly constant despite the dramatic changes in abundance, further evidence that immature neutrophils were released from the bone marrow during disease. Intriguingly, the fraction of dividing T and NK cells stayed relatively constant until increasing dramatically on DPI 8 for both NHPs that survived until then (RNA: Figure S3C; protein: Figure S3D, p = 0.022 rank-sum test of DPI 8 versus baseline for NK, CD8+ T, and CD4+ T). Proliferation is a core component of T cell mediated viral clearance but requires time for activated T cells to accumulate; the observation that significant T cell proliferation only occurred in the 2 NHPs at the latest DPI suggests that they may have begun to mount a T cell response.

12Not all cell types support EBOV entry and replication. We used scRNA-seq to identify infected cells in vivo based on the presence of EBOV’s RNA genome and poly-adenylated mRNA (Figure 2H). We developed a statistical model to identify infected cells containing more EBOV reads than expected due to ambient RNA contamination (STAR Methods) and estimated a mean sensitivity of 51% when ≥0.1% of cellular transcripts derived from EBOV (false positive rate = 1%, Figure S1H). We also spiked uninfected Madin-Darby canine kidney (MDCK) cells into a subset of PBMC samples to serve as a negative control (Table S1).

13Monocytes comprised the main infected cell population in vivo, first detectable at DPI 4 and with an increasing frequency thereafter (Figure 2I). Consistent with previous studies, T cells, B cells, and neutrophils were not identified as infected more often than expected by chance or than MDCK control cells. We did not observe any infected plasmacytoid (pDCs) or conventional dendritic cells (cDCs) in circulation, though infected DCs have been observed in culture and in lymph nodes.

14Monocytes are avid phagocytes and their intracellular viral RNA could reflect internalization of infected cell debris rather than infection. However, neutrophils also phagocytose debris from infected cells but do not support productive EBOV replication. Indeed, we do not detect neutrophils to be infected more often than expected by chance (Figure 2I), consistent with rapid RNA degradation inside of phagosomes. Moreover, infected monocyte transcriptomes predominantly contained viral mRNA rather than genomic RNA from extracellular virions (see later Results). These observations suggest that most cells with viral RNA reflect viral transcription rather than uptake of infected cell debris or viral particles, though we cannot exclude the latter possibilities in every cell.

15Having examined the varying abundance of each immune cell type, we next cataloged changes in their gene-expression profiles throughout EVD. We grouped cells into EVD stages: “incubation,” preceding detectable viral load or clinical signs (DPI 1 and 2; CyTOF only); “early,” detectable viral load but no clinical signs (DPI 3); “middle,” both viral load and clinical signs (DPI 4 and 5); and “late,” when NHPs reached euthanasia criteria (DPI 6–8) (Figure 1).

16We compared transcriptomes from each EVD stage to baseline for each cell type (STAR Methods) and identified 1,437 genes that were differentially expressed in at least one cell type and stage with a false discovery rate (FDR) corrected q value <0.05 and a fold change of >30% (Table S2). To identify patterns of gene expression associated with cell type and time, we performed unsupervised clustering of the differential expression signatures and identified 11 gene modules (Figures 3A and 3B; Table S3). We excluded neutrophils, pDCs, cDCs, and plasmablasts because of small sample sizes. Three modules that we term “global” were broadly up or downregulated compared to baseline across cell types, and the remaining modules were cell-type specific.

17The “global up” module contained 136 genes, consisting mostly of regulators and targets of the interferon (IFN) alpha (α) and gamma (γ) signal transduction cascade such as STAT1 and MX1. Gene sets labeled “response to IFN-α,” and “response to IFN-γ” were significantly enriched in this module (IFN-α: odds ratio [OR] = 69.5, q = 8 × 10−39; IFN-γ: OR = 45.9, q = 1 × 10−39; Fisher’s exact test; Table S3). This IFN response was consistent with a >10-fold increase of IFN-γ mRNA in CD8+ T cells from 0.4 transcripts per ten thousand (TP10K) at baseline to 5.0 at mid EVD (Figure S4A). Type 1 IFN (α/β) mRNAs also increased, but the change was not statistically significant, as IFN-α/β mRNAs are expressed transiently. To further characterize the IFN response, we scored each cell by the expression of literature-annotated ISGs. The median ISG score increased for all cell types and periods (p < 1 × 10−5, rank-sum test, Figure 3C).

18To determine whether IFN-α, IFN-γ, or both were driving this response, we identified genes that were annotated as regulated by only one of the two. Both uniquely IFN-α- and uniquely IFN-γ-regulated genes were significantly enriched in the “global up” module (q < 0.01, IFN-α: OR = 20.9; IFN-γ: OR = 16.6), a pattern that held true for each cell type and EVD stage (Figure S4B). These results suggest that both IFN-α and IFN-γ substantially and independently influenced the gene-expression profiles of circulating cells during EVD.

19The “global late down” module contained 144 genes that were predominantly downregulated during late EVD. It contained numerous regulators of translation (REACTOME_TRANSLATION gene set enrichment q = 5.2 × 10−7, Table S3), which is consistent with a core antiviral function of IFN being to downregulate translation.

20After elucidating the global effects of EBOV infection on immune cells, we next investigated the transcriptional responses specific to each cell type.

21The 2 modules “B/T early up” and “lymph late up” reflect changing gene expression of lymphocytes at different stages of EVD. “B/T early up” is strongly associated with the gene set HALLMARK_TNFA_SIGNALING_VIA_NFKB (q = 1.3 × 10−9), including the canonical activation marker CD69, CD48, and the transcription factor FOS. This module is unlikely to represent antigen-dependent activation as it occurs in most lymphocytes, and several of the top upregulated genes (e.g., GADD45B, ZFP36L2) are associated with growth arrest rather than proliferation. In addition, the 5th most enriched gene set in the “B down” module reflects BCR signaling (q = 0.00017), i.e., reduced BCR activation during EVD. Thus, the “B/T early up” module likely represents a cytokine-mediated, non-antigen-dependent activation of lymphocytes.

22The “lymph late up” module is upregulated in late EVD across all lymphocytes. The top associated gene sets implicate DNA repair (q = 0.00031) and apoptosis via TRAIL (q = 0.00032). This latter gene set is consistent with previous reports of T cell apoptosis in EVD and with the decline in lymphocyte abundance in our dataset (Figure 2D).

23Monocytes were of particular interest as they are well known to be central to EVD pathogenesis. Consistent with previous observations, we found that monocytes had far more significant gene-expression changes than the other cell types: 1,020 genes (11.6% of genes tested) in at least one EVD stage versus baseline, compared to 505 genes (6.6%) for B cells, which ranked second. We therefore focused our attention on characterizing monocytes.

24One prominent feature of the monocyte differential expression profile was the downregulation of several major histocompatibility complex (MHC) class II genes by mid and late EVD (Figure 3D). Monocytes and professional antigen-presenting cells display viral antigens on MHC class II proteins at the cell surface to stimulate the adaptive immune response. While IFN-γ typically upregulates MHC class II gene expression, we observed decreased MHC class II on monocytes despite elevated IFN-γ mRNA levels in T cells (Figure S4A) and widespread IFN-γ transcriptional response in monocytes (Figure S4B). Previous reports have described loss of HLA-DR, one of the 4 MHC class II proteins, during EBOV infection of monocytes ex vivo, and NHPs and humans in vivo. However, the specific MHC genes affected, cell-type specificity, temporal dynamics, and relationship with EBOV infection status, have not been previously described.

25We observed widespread changes in MHC expression throughout EVD (Figure 4A). The most striking decreases occurred in MHC class II genes of monocytes (>5-fold for DPA, DPB, and DRA by late EVD, all q < 1 × 10−21), with smaller changes in MHC class I genes (<1.7-fold increase for A, A3, and B, q < 5 × 10−4). B cells displayed modest reductions in MHC class II genes as well (>1.9-fold for DPA, DPB, DRA, and DQA1, q < 1 × 10−22). pDCs and cDCs showed no statistically significant reduction of any MHC class II gene (q > 0.05), but our dataset contained few DCs (Table S1) so we had less power to detect these effects. CyTOF revealed similar patterns in protein levels (p < 1 × 10−61 for monocytes at all stages, rank-sum test; p = 0.0012 for B cells at late EVD, Figures 4B and S4C). This phenomenon held true for each individual NHP: even as monocytes became activated, demonstrated by upregulation of the canonical activation marker CD38 (Figure 4C; Data S1), they showed dramatic downregulation of HLA-DR protein expression at DPI 5–8 (p = 9.5 × 10–7, Figure 4D).

26Reduced MHC class II expression in monocytes was not a direct consequence of EBOV infection of those cells. Only a small (∼5%) percentage of monocytes were infected at mid EVD (Figure 2I), too small to explain the size of the reduction. In fact, expression of MHC class II genes was comparable or higher in infected cells than in uninfected (bystander) cells (Figure 4E).

27We identified genes correlated with MHC class II in monocytes (STAR Methods), as they may be part of a co-regulated expression program. Many of the most correlated genes are functionally involved in antigen presentation, including CD74, LGMN, and B2M (Figure 4A; Spearman ρ = 0.42, 0.41, 0.33, respectively, all p < 10−170). One of the most associated genes was ZFP36 (Spearman ρ = 0.43, p < 10−296), which directly regulates turnover of MHC class II and other immune-related mRNAs. Thus, MHC class II and other genes involved in antigen presentation may be part of a single transcriptional module, co-regulated by ZFP36 and/or other genes.

28Next, we characterized genes that were differentially expressed between infected and bystander monocytes, as these could represent host entry factors, restriction factors, or genes that are regulated by infection within a cell. For this and subsequent analyses, EBOV transcripts were excluded from library size normalization to avoid bias in estimating host gene expression in infected cells (STAR Methods). 505 genes were differentially expressed between infected and bystander monocytes (q < 0.05) of which 276 changed by >30% (Figure 5A; Table S4).

29We observed 3 broad categories of differentially expressed genes: those associated with monocyte subtypes, monocyte-to-macrophage differentiation, and ISGs.

30CD14 and FCGR3 (which codes for CD16) were overexpressed in infected monocytes relative to bystanders. These two genes define classical (CD14+) and non-classical (CD16+) monocytes, which are the dominant blood monocyte subsets in healthy individuals. Classical monocytes are highly phagocytic scavenger cells, while non-classical monocytes are involved in complement and antibody-mediated phagocytosis. At baseline, these subsets were easily detectable as 2 distinct clusters, marked by high expression of CD14 (Figure 5B, bottom half of the blue lobe of DPI panel) or CD16 (top half of the blue lobe).

31However, monocyte subsets changed dramatically during EVD. Both single-positive CD14+ and CD16+ monocytes declined, while 2 unusual populations emerged: a large population of CD14– CD16– cells (double negatives [DNs]) and a smaller population of CD14+ CD16+ cells (double positives [DPs]) (Figure 5C). While 87.6% of cells fell into single CD14+ or CD16+ bins at baseline, this dropped to 33.8%, 35.6%, and 15.2% in early, middle, and late EVD, respectively. We confirmed this pattern at the protein level by CyTOF (Figure 5D).

32At late EVD, the most frequent population was DN cells, which rose to make up 70.4% (by scRNA-seq) or 56.7% (by CyTOF) of monocytes. As they expressed neither CD14 nor CD16, we confirmed that their gene-expression profiles were most correlated with conventional monocytes and not other cell types (Figure S5A). DNs first emerged on DPI 3, coinciding with the 2-fold increase in monocytes (Figure 2D), and with the increased monocyte proliferation observed in mid EVD (Figure 2G). The DN population is in fact highly proliferative: while 0% of monocytes at baseline expressed elevated levels of Ki67 (smoothed log TP10K >1), 37% of DN monocytes exceeded this threshold by late EVD (Figures S5B and S5C). Therefore, DNs underlie the increased monocyte proliferation in mid EVD.

33To determine whether DNs arise in human EVD, we re-analyzed published CyTOF data from 4 patients during the 2013–2015 epidemic (STAR Methods). All of the human cases showed a loss of conventional CD14+ and CD16+ single-positive monocytes and an emergence of proliferative DN monocytes during acute disease (Figures 5E and S5D). Then, the DNs disappeared and were replaced by conventional CD14+ and CD16+ monocytes as the patients recovered. Thus, DNs are a feature of human sub-lethal cases and our lethal NHP model.

34Mature monocytes in circulation are non-dividing, but infection and cancer can induce the release of proliferating immature myeloid cells from the bone marrow, a process known as emergency myelopoiesis. We therefore hypothesized that the DN population may derive from emergency myelopoiesis.

35If DNs result from emergency myelopoiesis, their gene expression may be more similar to bone marrow monocyte precursors (BM-MPs) than circulating monocytes. To test this, we analyzed a reference scRNA-seq dataset of BM-MPs from healthy human bone marrow and mature monocytes from human PBMCs (Figures S5E–S5H; STAR Methods). Compared to mature monocytes, BM-MPs showed lower expression of CD14 and FCGR3A (the human CD16 gene) and higher expression of MPO, AZU1, S100A8, and S100A9 (Figures S5I and S5J). In our data, the expression of those genes in DNs at late EVD relative to baseline monocytes mirrored BM-MPs, suggesting their similarity (Figure S5K). Moreover, in a nearest-neighbor test (STAR Methods), DNs were more likely, than conventional subtypes, to be matched with BM-MPs (Figure 5F). Thus, DNs may represent immature monocytes released from the bone marrow in response to the EVD cytokine milieu.

36In addition to DNs, we also observed CD14+ CD16+ DP cells, which rose to make up 20.4% of the monocytes by mid EVD (Figure 5C). A similar increase in DP monocytes has been observed in other viral infections. We found that the DP population harbored a disproportionately high percentage of EBOV-infected cells (Figure 5G), consistent with the fact that CD14 and CD16 were each higher in EBOV-infected monocytes than in bystanders (Figure 5A). At late EVD, 22.1% of DPs were infected compared to only 1.7% of DNs. Thus, the differential expression of CD14 and CD16 in infected cells results from increased infection of the DP cells, rather than increased expression of CD14 on classical and CD16 on non-classical monocytes.

37Differentially expressed genes between infected and bystander monocytes (Figure 5A; Table S4) and between DP and DN monocytes (Figure S6A; Table S5) were enriched for monocyte-to-macrophage differentiation associated genes, including known EBOV entry factors. Freshly isolated monocytes are largely refractory to EBOV infection in cell culture, but EBOV entry factors are upregulated during in vitro macrophage differentiation, allowing increased infection. In vivo, we observed higher levels of macrophage differentiation markers in infected cells than in bystanders, including the known EBOV entry factors cathepsin L (CTSL) and B (CTSB), and GNPTAB (q = 6.7 × 10−9, 3.8 × 10−7, and 2.1 × 10−3, respectively). By contrast, the cellular receptor NPC1 was not significantly differentially expressed, suggesting that natural variability in NPC1 mRNA does not influence EBOV infectivity within circulating monocytes in vivo.

38To determine whether upregulation of these factors was part of a general macrophage differentiation program, we used gene sets derived from published bulk RNA-seq data of primary blood monocytes before and after differentiation into macrophages in vitro. Genes that are upregulated during in vitro differentiation were significantly enriched in infected cells (OR = 3.5, p = 3.1 × 10−11, Fisher’s exact test), while downregulated genes were significantly enriched in bystanders (OR = 3.7, p = 4.2 × 10−8, Fisher’s exact test; combined chi-square goodness of fit test p = 2.2 × 10−30). Additional differentiation studies yielded similar findings (chi-square goodness of fit test p = 2.6 × 10−9, Fisher’s exact test p = 5.7 × 10−12). Genes associated with differentiation into M2-polarized macrophages were more enriched among infected cells than those of M1-polarized macrophages (OR = 7.8, p = 1.3 × 10−10 versus OR = 3.3, p = 1.6 × 10−3, Table S4).

39The proportion of infected cells increased along with a score reflecting the expression of genes associated with macrophage differentiation (STAR Methods). Over the range of scores, we observed that the percentage of infected cells rose ∼5-fold from 3% to 15% (Figure 5H).

40Among the monocyte subsets, DPs generally had the highest macrophage scores (Figures S6A and S6B). Thus, some of the enrichment of infected cells among the DP subset could be attributed to their more macrophage-like gene expression.

41However, there was substantial heterogeneity in macrophage scores within monocyte subsets, and macrophage score and CD14/CD16 subset were independently predictive of infectivity. Stratifying cells in each subset by macrophage score, less macrophage-like DPs were still more likely to be infected than more macrophage-like DNs (p = 1.9 × 10−25, Fisher’s exact test, Figure S6C). Yet, within the DPs, more macrophage-like cells were more likely to be infected than less macrophage-like cells (p = 0.0003, Fisher’s exact test, Figure S6C). Similarly, when we fit a logistic regression predicting the infection status of each cell using macrophage score, smoothed CD14 and CD16 expression values, and a CD14xCD16 interaction term (STAR Methods), both the CD14xCD16 interaction term (which is highest in DPs) and macrophage score were positively associated with infection, while the CD14 and CD16 terms were negatively associated (p < 0.01 for all coefficients).

42Finally, we noticed that several key ISGs such as MX1 were expressed at lower levels in infected cells than in bystanders (MAST q = 7.7 × 10−14, rank-sum test p = 2.0 × 10−5, Figures 5A and 5I). To determine whether infection suppressed overall ISG expression, we assessed the magnitude of the IFN response between infected and bystander cells at late EVD. While both bystander and infected monocytes at late EVD had higher ISG scores than monocytes at baseline, ISG scores were lower in infected cells than bystanders (not statistically significant by rank-sum test, Figure S6D). More strikingly, there was a significant negative correlation between ISG score and the percentage of cellular transcripts derived from EBOV (i.e., the intracellular viral load) (Spearman ρ = −0.62, p = 1.1 × 10−11, Figure 5J). This suggests that ISGs are downregulated during viral replication within infected cells (see later Results).

43In order to more thoroughly probe viral and host gene expression during the viral life cycle, we isolated PBMCs from 2 healthy rhesus monkeys (NHP1 and NHP2) and inoculated them ex vivo with live EBOV, EBOV rendered replication-incompetent by gamma irradiation, or only media as a control (Figure 6A). We performed scRNA-seq using Seq-Well at 4 h or 24 h post-infection (HPI), corresponding to early (start of viral transcription) and middle-to-late stages (viral genome replication, virion assembly) of the viral life cycle. Unlike alternative systems such as virion-like particles, inoculation with gamma-irradiated EBOV allowed us to profile cells containing viral RNA and characterize the host response in the absence of effective viral transcription and translation.

44We obtained single-cell transcriptomes from 50,646 PBMCs ex vivo and observed similar cell-type representation, clustering by treatment condition, and distribution of EBOV-infected cells as with the in vivo collections (Figures S7A–S7C). ISG scores were higher in NHP1 than NHP2 (Figures S7D–S7G), so we analyzed cells from each NHP separately and jointly to avoid potential artifacts. While we did not observe infected DCs in vivo, we found that 16% of DCs (16.0% in NHP1, 15.7% in NHP2) inoculated with live virus ex vivo were infected by 24 HPI (Figure S7H).

45Ex vivo, monocytes were the predominant infected cell type, with >65% infected by 24 HPI (Figures S7H and S7I). 11.8% of the monocytes treated with irradiated virus also contained a statistically significant number of viral reads, even though gamma irradiation damages the viral genome and eliminates productive replication. As expected, these cells had a significantly lower fraction of EBOV reads per cell than those treated with live virus (Figure S8A). Moreover, viral RNAs from the cells treated with irradiated virus were substantially less likely to be coding-sense mRNA transcripts than those from cells with live virus (44% versus 92%, Figure S8B). Thus, our method can detect fragments of viral genomic RNA from irradiated virus in cells, but these are not productive infections.

46The intracellular viral load varies over several orders of magnitude (Figures 6B and 6C). While most cells harbored viral loads <0.1%, a substantial minority had loads >10%, with a maximum of 57.5% in vivo and 52.3% ex vivo. The observed heterogeneity was not due to differential transcriptome coverage, because cells with low and high viral load had a similar range of transcripts detected (Figures S8C and S8D).

47The ex vivo data gave us an opportunity to test predictions based on established models of EBOV transcription. Transcription of EBOV’s 7 genes by the viral RNA-directed RNA polymerase L follows the canonical stop-start mechanism of filoviruses and other non-segmented negative-strand (NNS) RNA viruses. L initiates transcription at the 3′ end of the genome, and processes from 5′ to 3′; at each gene’s transcription termination signal, L either falls off the genomic RNA template or reinitiates transcription of the subsequent gene (Figure 6D; ). NP is transcribed first and at the highest level, proceeding down the genome to L last and at the lowest level.

48When we quantified the relative expression of EBOV genes as a function of viral load, we observed an unexpected accumulation of GP mRNA (Figures 6E and 6F) in both NHP1 and NHP2 (Figures S8E and S8F). At low viral loads, both in vivo and ex vivo, the gene-expression distribution roughly matched the expected pattern, with most transcripts derived from NP, and the fewest from the 5′ genes VP30, VP24, and L. In agreement with this pattern, cells inoculated with irradiated virus were highly enriched in NP mRNA (Figure S8B), suggestive of RNA fragment transcription. However, as viral load increased in cells infected with live EBOV, GP unexpectedly became the most abundant viral transcript. While higher levels of GP than NP RNA have been observed in previous bulk RNA-seq datasets, our data indicate that this phenomenon occurs within individual cells, as a life-cycle-dependent regulation of viral gene expression.

49Next, we exploited natural variability in viral load to identify host gene-expression changes correlated with viral replication, which may represent pathways directly regulated by infection. Instead of testing for differential expression between infected and bystander cells as we did previously, we looked for continuous association between viral abundance and host transcript levels in infected monocytes (STAR Methods). This identified 264 genes that were negatively correlated and 211 genes that were positively correlated with viral load ex vivo (q < 0.05), of which 34 changed by >30% per 10-fold increase in viral load (Figure 7A; Table S6).

50Consistent with our previous observation that ISG score decreased with intracellular viral load, individual ISGs decreased dramatically with EBOV levels, both ex vivo and in vivo (e.g., MX1, q = 1.5 × 10−24, 7.4 × 10–9, respectively, Figures 7A and 7B). Ex vivo, the most negatively associated gene was STAT1, the master transcriptional regulator of the IFN response (q = 4.9 × 10−41 ex vivo, p = 0.0072 in vivo but not significant by FDR). Previous experiments have shown that the EBOV protein VP24 inhibits STAT1 activity by blocking its translocation to the nucleus. However, this is the first observation that STAT1 mRNA levels decrease with viral replication within infected cells in vivo.

51The expression of STAT1 and other negatively regulated ISGs remained relatively constant as EBOV levels rose to 1% of cellular transcripts but declined precipitously with higher EBOV levels (Figure 7C). This suggests that there is a delay before EBOV can downregulate host antiviral genes since it must transcribe and translate VP24 and other immunomodulatory proteins before they can act. The trajectories of these host antiviral genes were consistent between donor NHPs ex vivo (Figure S8G), even though more cells from NHP1 mounted an IFN response than NHP2.

52Only a handful of host genes increased in expression alongside viral load, but their trajectories were consistent between the two NHPs (Figure S8H). The most dramatically upregulated gene was DYNLL1 both ex vivo and in vivo (q = 2.5 × 10−27, 1.5 × 10–5, respectively, Figures 7A and 7B). DYNLL1 is a multi-functional protein involved in intracellular transport. Intriguingly, DYNLL1 was previously shown to increase EBOV replication in a minigenome reporter assay. Our data show that DYNLL1 mRNA is upregulated, starting when EBOV RNA constitutes 0.1%–1% of transcripts, before ISGs decrease (Figure 7D).

53Several other genes that we identified as upregulated alongside viral replication have known or speculated pro-viral functions in protein folding and synthesis. For example, HSPA5 (q = 4.5 × 10−22 ex vivo, p = 0.04 in vivo but not significant by FDR) encodes a chaperone protein that is an essential host factor for EBOV, but upregulation of its mRNA has not been previously observed. Other hits such as DDIT3 and NFE2L1 (q = 1.3 × 10–10, 1.8 ×10−8 ex vivo, respectively) are sensors of ER and oxidative stress and have been implicated in cell lines infected with Marburg virus and in monocytes in EBOV-infected NHPs. We observe a corresponding enrichment of 3 gene sets associated with ER stress response, most notably BUYTAERT_PHOTODYNAMIC_THERAPY_STRESS_UP (q = 4.7 × 10−14, Fisher’s exact test, Table S6). Increased translation due to viral replication can drive cell stress and also deplete tRNAs and free amino acids. Upregulation of IARS—isoleucine tRNA synthetase (q = 8.4 × 10−9 ex vivo)—may reflect this cellular burden. 3 gene sets associated with depletion of amino acids were significantly upregulated, most notably KRIGE_RESPONSE_TO_TOSEDOSTAT (aminopeptidase inhibitor)_24HR_UP (q = 2.3 × 10−6, Fisher’s exact test, Table S6), suggesting that viral replication exhausts cellular amino acid stores. This hypothesis is consistent with prior observations of depleted amino acids in the plasma of fatal human EVD cases.

54Despite recurrent outbreaks, the molecular basis of EVD pathogenesis remains understudied due to the biosafety and logistical challenges of such research. By adapting CyTOF and scRNA-seq approaches for use in BSL-4 containment, we comprehensively surveyed the molecular correlates of disease progression and viral replication in circulating immune cells in a NHP model of EVD. This study, which is the first high-parameter, single-cell investigation under BSL-4 containment, shed new light on changes in peripheral cell-type and -state abundance throughout lethal EVD, characterized the EBOV-infected populations of circulating immune cells, and identified genes regulated by the cytokine milieu or by direct EBOV infection.

55We characterized transcriptional- and protein-level changes in monocytes during EVD in NHPs, some of which reflect disruption of their physiological antiviral function. In agreement with the well-established importance of monocytes in EVD, we observed that they had over twice as many differentially expressed genes as other cell types, including genes involved in IFN response, cytokine production, myeloid differentiation, and antigen presentation. Monocytes became activated by IFN during EVD, which normally upregulates MHC class II genes. Surprisingly, almost all MHC class II and related antigen-presentation genes were strikingly downregulated. Moreover, MHC class II expression decreased in both infected and uninfected monocytes, suggesting that the decrease was caused by the cytokine milieu. Reduced monocytes antigen presentation might explain why a failed or delayed adaptive immune response is a hallmark of fatal EVD in humans.

56As EVD progressed, conventional CD14+ and CD16+ monocyte subsets largely disappeared and were replaced by two populations: CD14+ CD16+ (DP) monocytes, which increase in other infections, and an unexpected CD14– CD16– (DN) population, that, to our knowledge, has not been previously described in viral infection. The DN monocytes were highly proliferative, with transcriptomes more similar to bone marrow monocyte precursors than circulating monocytes. This suggests that they are the product of emergency myelopoiesis, in which cytokines stimulate the bone marrow to release immature myeloid-lineage cells. DNs had high expression of neutrophil granule genes such as MPO and AZU1, suggesting that they represent myeloid progenitors prior to the branching of neutrophil and monocyte lineages. We also identified the emergence of an analogous DN population in human acute EVD cases. This finding highlights the power of high-parameter methods such as scRNA-seq and CyTOF: despite little-to-no detection of CD14 or CD16, there were enough other RNA and protein markers to reliably identify these DNs as monocyte precursors.

57Our data refine the picture of EBOV’s tropism in NHPs, demonstrating that the predominant infected population in circulation are DP monocytes expressing markers of macrophage differentiation. Myeloid cells are known to be major targets of EBOV, including DCs ex vivo and in lymph nodes in vivo. While we too observe infected DCs ex vivo, in vivo we only detected infected monocytes among circulating immune cells. This might reflect biological phenomena required for DC infection, such as cell density, cell-to-cell contact, or viral dose. Viral tropism heavily depends on viral entry, during which the cellular receptor NPC1 and other entry factors are required. While we did not observe differences in NPC1 mRNA expression between infected and uninfected monocytes, expression of known EBOV entry factors like cathepsin B and L were associated with infection. Cultured monocytes are only susceptible to EBOV infection upon differentiation, and our data further show that infectivity in vivo strongly correlates with physiological variability of monocyte differentiation state. Furthermore, our data show that the relative abundance of DP monocytes, the preferred circulating cell targets of EBOV, increases over the course of infection, perhaps driven by the cytokine milieu of EVD.

58Among infected monocytes, we observed substantial heterogeneity in intracellular viral load (up to 57.5% of transcripts), which we exploited as a proxy for staging EBOV’s progression through its life cycle. This enabled us to better understand the over-representation of GP mRNA that had been observed in bulk RNA-seq datasets. In our data, cells with low viral load (reflecting early infection) had EBOV transcript abundances that mirrored the genome organization as expected—i.e., expression was highest for NP and lowest for L. But cells with high viral load had higher abundance of GP than NP mRNA, suggesting that alternate transcription or post-transcriptional regulation, such as increased mRNA stability, may account for accumulation of GP late in the viral life cycle. Translation of GP, a structural protein, at the end of the viral life cycle coincides with virion assembly and thus may benefit EBOV fitness, as has been observed for other viruses.

59Many host ISGs negatively correlated with intracellular viral load, and our data suggest that viral infection downregulates ISG expression in vivo, rather than preferentially infecting cells with low ISG levels. First, there are multiple established mechanisms by which EBOV downregulates ISGs, such as by preventing STAT1 from translocating to the nucleus. Second, monocytes mount a strong ISG response by DPI 3, yet the percentage of infected cells rises gradually from DPI 4 onward, implying that EBOV is able to replicate despite the inhibitory activities of ISGs. Third, we do not observe any cells with low (0.01%–0.1%) viral load and low ISG levels. Since a newly infected cell begins with low viral load, this observation implies that EBOV is infecting monocytes that are mounting a full IFN response. Because EBOV potently suppresses IFN production in infected cells, it remains unclear how the IFN response is triggered. Possibilities include incomplete IFN antagonism in infected cells, defective particles, or genomes that fail to suppress IFN, pathogen- and danger-associated molecular patterns released from dying infected cells, and other classes of cytokines produced from cells that do not support EBOV replication. Better understanding this antiviral response could be the key to novel therapeutics or vaccines.

60Several putative pro-viral host genes were positively correlated with intracellular viral load, suggesting that they are directly responding to, or are regulated by, the presence of virus in infected cells (Figures 7A, 7B, and 7D). DYNLL1 was the gene most associated with viral load. Previous work identified an interaction between DYNLL1 and EBOV VP35 that increased EBOV RNA synthesis in a minigenome assay. DYNLL1 increases replication of rabies virus an NNS virus similar to EBOV. Here, we showed that DYNLL1 expression is upregulated within infected cells in vivo, suggesting that EBOV manipulates cellular pathways to encourage a pro-viral cellular environment. Nuclear DYNLL1 typically represses its own transcription factor ATMIN; we hypothesize that EBOV VP35 sequesters DYNLL1 protein in the cytoplasm, relieving repression of ATMIN, thus upregulating DYNLL1 mRNA. Additional genes that were upregulated with intracellular viral load were associated with translation and cell stress (e.g., HSPA5, IARS), reflecting the strain placed on cells by viral replication. By computationally staging cells within the viral life cycle, we nominated several putative pro-viral genes for further study, highlighting the utility of single-cell profiling to study host-virus interactions.

61Our data provide insights into circulating immune cells in a lethal model of EVD in rhesus monkeys; however, there are several key aspects of pathogenesis that are not reflected in this study. Much of the clinical syndrome of EVD is due to viral replication and pathology in tissues, including the liver, vasculature, lung, kidney, lymph nodes, and bone marrow. Studies of these compartments would be complementary to ours as immune cells pass through the circulation and extravasate into tissues to fight infection. An additional limitation of our study is that human EVD is not always lethal and has slower clinical kinetics than the uniformly lethal model used here. Studying EVD in survivors, either in human clinical cases or in non-uniformly lethal animal models, could substantially broaden our understanding of pathogenesis and successful adaptive immunity.

62The accumulation of additional host-pathogen single-cell datasets promises to greatly enhance our understanding of infection by allowing us to determine which features of pathogenesis are shared between, or specific to, individual pathogens. For example, our scRNA-seq and CyTOF data identified several molecular commonalities between EVD and immunosuppressive septic shock, which is also characterized by loss of MHC class II expression in monocytes, increased DP monocytes, and emergency myelopoiesis. Soluble mediators, including cytokines and glucocorticoids, could be key drivers of both EVD and sepsis pathophysiology. Tumor necrosis factor (TNF)-α signaling has been extensively implicated as a driver of systemic loss of vascular resistance and shock during EVD. Glucocorticoids have been less well studied in EVD but decrease MHC class II during sepsis and reduce CD14 expression while increasing the abundance of DP monocytes. Indeed, the connection between EVD and sepsis may be direct: studies have found bacterial invasion during EVD in NHPs and in humans, with immune signatures that resemble sepsis.

63In summary, this work expands our understanding of EVD and provides a general paradigm for exploring molecular features of host-pathogen interactions, such as tropism and dysregulation of cell circuitry in infected cells, which can be applied to other emerging pathogens.

64Further information and requests for resources and reagents should be directed to Aaron Lin (alin@broadinstitute.org).

65This study did not generate new unique reagents.

66The analysis scripts used in this study are available at https://github.com/dylkot/SC-Ebola.

67The accession number for single-cell RNA-Seq datasets reported in this paper is GEO: GSE158390, with raw sequence data available on SRA.

68The accession number for CyTOF.fcs files reported in this paper is FlowRepository: FR-FCM-Z2LX.

69This study included a subset (21 of 27) outbred rhesus monkeys (Macaca mulatta) of Chinese origin described recently. These 27 nonhuman primates (NHPs) were randomized into cohorts (Figure S1A,), balancing age, weight, and sex across 7 groups. All work was approved and performed in accordance with the Guide for the Care and Use of Laboratory Animals of the National Institute of Health, the Office of Animal Welfare, and the US Department of Agriculture.

70This study utilized the Ebola virus/H. sapiens-tc/COD/1995/Kikwit-9510621 (EBOV/Kikwit; GenBank accession MG572235.1; Filoviridae: Zaire ebolavirus) isolate for the in vivo and ex vivo challenges, obtained from the Biodefense and Emerging Infections Research Resources Repository (BEI Resources, Manassas, VA, USA). It is the standard challenge stock defined by the filovirus animal non-clinical group (FANG) for testing product efficacy for FDA approval and is well characterized.

71For all 21 outbred rhesus monkeys, two baseline blood samples were collected between 0–14 and 14–30 days prior to infection (Figure S1A). 18 NHPs were exposed to the EBOV/Kikwit isolate diluted to a target concentration of 1,000 plaque forming units (PFU) in a volume of 1 mL/dose. All NHPs were inoculated within a 5 month period. This same cohort has already been described recently.

72Beginning on day post-infection (DPI) 0, NHPs were observed 1–3 times daily and given a clinical score based on five criteria: overall clinical appearance and signs of hemorrhage; respiratory rate, mucous membrane color, and dyspnea; recumbency; non-responsiveness; and core temperature. Each criterium was assigned a score between 1 and 10, and all scores were added together. Once an NHP reached a combined score of > 10, the animal was humanely euthanized.

73Blood was drawn from anesthetized animals into BD vacutainer plastic serum separator tubes (SST) for serum viral load quantification, or in BD vacutainer plastic blood collection tubes with K3EDTA for hematology and peripheral blood mononuclear cell (PBMC) purification (Becton Dickinson, Franklin Lakes, NJ, USA). SST tubes were centrifuged at room temperature for 10 minutes (min) at 1800 x g to isolate serum. K3EDTA tubes were mixed by gentle inversion prior to hematology and PBMC purification.

74250 μL of each whole blood sample was analyzed on a Sysmex 2000i XT (Sysmex Corporation, Kobe, Hyogo Prefecture, Japan). Parameters analyzed by this instrument were: counts of basophils, eosinophils, lymphocytes, monocytes, neutrophils, white blood cell count; percentages of each cell type; and mean platelet volume.

75To estimate the abundance of lymphocyte cell types, we multiplied the CBC lymphocyte count by the proportion of lymphocytes of each cell type (CD8 T cells, CD4 T cells, NK cells, and B cells) which was calculated from the unsupervised clustering of the CyTOF data (see ‘CyTOF’ below).

7670 μL of sample inactivated by TRIzol LS was added to 280 μL of Buffer AVL (QIAGEN, Hilden, Germany) with carrier RNA. Samples were then extracted using the QIAamp Viral RNA Mini Kit (QIAGEN) in accordance with the manufacturer’s instructions, eluted in 70 μL of Buffer AVE, aliquoted, and frozen. Viral load was determined using the BEI Resources Critical Reagents Program experimental EZ1 reverse transcription qPCR kit assay in accordance with the manufacturer’s instructions.

77We centrifuged whole blood in K3EDTA tubes at 1800 x g for 10 min at room temperature, removed EDTA plasma, and added phosphate buffered saline (PBS, Thermo Fisher Scientific, Waltham, MA, USA) to the pelleted cells to double the original whole blood volume. We gently poured the PBS-blood cell mixture into an Accuspin tube containing Histopaque (Sigma-Aldrich, St. Louis, MO, USA) and centrifuged at 1000 x g for 10 min at room temperature with the brake set to 1. Following centrifugation, we removed the top, clear supernatant layer to within 0.5 cm of the cloudy white layer containing PBMCs.

78We transferred the cloudy PBMC layer to a clean 15 mL conical tube and increased the volume to 10 mL using PBS supplemented with 2% heat-inactivated fetal bovine serum (PBS/2%HI-FBS) and mixed by inversion. We then centrifuged at 300 x g for 10 min at 4°C with the brake set to 1. Following centrifugation, we removed the supernatant, resuspended the cell pellet with PBS/2%HI-FBS to a final volume of 10 mL, and mixed using gentle raking to wash the cells. We repeated the wash step 2 more times with the centrifuge set to 200 x g for 10 min at 4°C with the brake set to 1. We then resuspended the cell pellet in 9.5 mL PBS/2%HI-FBS for counting using the Countess Cell Counting system (Thermo Fisher Scientific). We aliquoted 0.5 mL for Seq-Well, and used the remaining 9 mL volume for CyTOF.

79We performed Seq-Well as described previously, with the S3 protocol and some controls and modifications to adhere to the BSL-4 environment.

80As an experimental negative control to test our statistical model, we spiked Madin-Darby canine kidney (MDCK) cells, constituting ∼5% of the total sample, into a subset of PBMC samples (Table S1) from EVD NHPs immediately before scRNA-Seq. As MDCKs were not exposed to EBOV, viral reads in these transcriptomes should be due to ambient RNA contamination.

81After loading and sealing beads and cells in Seq-Well arrays, we placed them in a −80°C freezer until further processing – this step was required due to time constraints in the BSL-4. Later, we removed sealed Seq-Well arrays from the −80°C freezer, placed them in 4-well dishes, and allowed them to equilibrate to room temperature for at least 30 min. We then covered arrays in 5 mL Seq-Well Lysis Buffer per protocol.

82We performed RNA hybridization and RT as specified in the protocol. After RT, we collected beads by centrifugation at 1000 x g for 1 min at room temperature. We resuspended beads with GeneXpert Lysis Buffer (Cepheid, Sunnyvale, CA, USA) for inactivation, which was required prior to removal from the BSL-4 laboratory according to standard operating procedures. After removal, we washed beads thrice with TE buffer containing 0.01% Tween 20 and shipped at 4°C for further library construction and sequencing, which was performed with the S3 protocol. We sequenced all libraries on either NextSeq 550 High Output or NovaSeq 6000 S2 flowcells (Illumina, San Diego, CA, USA), with 20 cycles for Read 1 (cell barcode and unique molecular index [UMI]) and 88 cycles for Read 2 (cDNA of interest). In some cases, we merged fastq files from multiple sequencing runs for increased coverage.

83We observed long concatemers of the common scRNA-Seq adaptor sequence (SeqB, 5′-AAGCAGTGGTATCAACGCAGAGTAC-3′) at high frequency in some Seq-Well libraries, likely owing to low RNA input and the challenging environment of processing samples in the BSL-4 suite. These concatemers disrupted Illumina sequencing runs because the Read 1 sequencing primer annealed to multiple SeqB sequences on a single template, allowing multiple sequencing-by-synthesis reactions simultaneously. We therefore devised a strategy to remove SeqB concatemers using depletion of abundant sequences by hybridization (DASH), a CRISPR-based method to degrade target DNA sequences prior to sequencing. SeqB lacked a ‘NGG’ protospacer adjacent motif (PAM) for the common S. pyogenes Cas9 (SpyCas9) for which DASH was originally described; therefore, we modified DASH to use S. aureus Cas9 (SauCas9). Moreover, in contrast to SpyCas9, SauCas9 is a multi-turnover enzyme, suggesting that SauCas9 would have higher cleavage efficiency, which was important since SeqB was present in multiple copies within a concatemer.

84First, we designed a guide RNA (gRNA) to target SeqB. Based on the position of the SauCas9 PAM and the length of SeqB, only 17 nucleotides of the gRNA protospacer could anneal to SeqB. Because gRNA length is critical to SauCas9 cleavage efficiency, we prepended 4 random bases as a 5′ overhang (Key Resources Table). We in vitro transcribed this gRNA using the MEGAshortscript T7 Transcription Kit (Thermo Fisher Scientific), purified it using the RNA Clean & Concentrate Kit (Zymo Research, Irvine, CA, USA), and verified the correct RNA length on a 15% TBE-urea gel (Bio-Rad Laboratories, Hercules, CA, USA).

85We performed SauCas9 DASH according to reaction conditions laid out for in vitro SauCas9 cleavage assays. We incubated 10 pmol gRNA with 5 pmol SauCas9 (New England Biolabs [NEB], Ipswich, MA, USA) at 25°C for 10 min, and then added up to 5 fmol DNA (2000:1000:1 RNA:protein:DNA ratio) and NEBuffer 3.1 (NEB) to 1X. We incubated this reaction at 37°C for 2 hours (h), and quenched by adding EDTA to 50 μM, SDS to 1%, and 4 total U of Proteinase K (NEB) at room temp for 10 min. We removed degraded concatemers with two consecutive 0.8X SPRI purifications using Ampure XP beads (Beckman Coulter, Brea, CA, USA), eluted, and performed 6–9 cycles of PCR with the NEBNext Ultra II Q5 Master Mix (NEB) using Illumina P7 and the Seq-Well P5-TSO hybrid primer.

86We added 1 mL of 16% paraformaldehyde (PFA, Electron Microscopy Sciences, Hatfield, PA, USA) to 9 mL PBMCs to fix the cells. We incubated samples at room temperature for 10 min followed by a final centrifugation at 600 x g for 5 min at 4°C with the brake set to 9. We removed the supernatant, added 1 mL of PBS/5%HI-FBS for every 3 × 106 cells (e.g., 2 mL for a sample containing 6 × 106 cells), aliquoted into 1 mL aliquots in cryovials, and stored in a −80°C freezer.

87We equilibrated fixed PBMCs to come to room temperature before transferring approximately 2 × 106 cells per sample into 1.2 mL cluster tubes in a 96-tube rack. We barcoded samples and multiplexed them into 6 batches of 16 samples using a previously described 16-plex palladium-based mass-tag cell barcoding scheme. We aspirated pelleted barcoded cells to a volume of 50 μL and incubated them with 15 μL of Human TruStain FcX (Biolegend, San Diego, CA) for 10 minutes. We stained cells for 30 min with 175 μL of a reconstituted lyophilized cocktail of metal-tagged cell-surface antibodies described previously supplemented with two additional antibody channels (Data S1). We washed and permeabilized surface-stained cells with methanol before aspirating down to 50 μL and staining for 60 min with 190 μL of reconstituted lyophilized intracellular staining antibody cocktail (Data S1). We washed and resuspended fully-stained cells in a volume of 750 μL.

88Following staining, we inactivated samples by adding 250 μL of 16% PFA to 750 μL of each sample, for a final concentration of 4% PFA, and incubated at 4°C overnight. The next day, we centrifuged samples 600 x g for 5 min at 4°C and aspirated down to 100 μL. We resuspended samples in 1 mL 4% PFA in PBS and transferred to a clean 2 mL cryovial. We then removed samples from the BSL-4 using a dunk tank and froze at −80°C within 30 min of PFA addition.

89Following inactivation, we thawed and processed inactivated samples within a BSL-2 lab for iridium intercalation, then mixed with 1xEQ beads (Fluidigm, South San Francisco, CA, USA) and run on a CyTOF Helios (Fluidigm) instrument using a Super-Sampler introduction system (Victorian Airship & Scientific Apparatus LLC, Alamo, CA, USA).

90FCS data files were normalized across all runs using the data normalization software and debarcoded using the single-cell debarcoder tool as previously described. Data were uploaded and analyzed using CellEngine software (https://cellengine.com, Primity Bio, Fremont, CA, US), and a gating strategy was applied to identify cell populations using canonical markers. Frequencies for each population were determined as a function of total CD66-CD45+ cell events and reported marker intensities are expressed as medians. Also see ‘CyTOF data preprocessing, clustering, and dimensionality reduction’ section below for details of unsupervised analysis of CyTOF data.

91PBMCs were isolated from healthy NHPs as previously described. We diluted cells to 3.3 × 106 cells/mL RPMI/10%HI-FBS and transferred 900 μL each to 2 mL external thread cryogenic vials (Corning, Corning, NY) for each experimental condition. We inoculated cells with 100 μL live virus (EBOV/Kikwit, the same stock used for in vivo NHP inoculation) for a final MOI of 0.1 PFU/cell, an equivalent dose of irradiated virus (EBOV/Kikwit treated with 5 mRads gamma irradiation), or media only, and incubated for 4 or 24 hours with slow rocking prior to Seq-Well processing.

92Raw sequencing files were demultiplexed and converted to fastq using bcl2fastq version 2.20. Reads were then trimmed, aligned to a reference transcriptome, and parsed into a digital gene expression matrix using the previously published Dropseq-tools pipeline version 2.0. We used a Snakemake wrapper around Dropseq-tools that is available in the open source Github repository https://github.com/Hoohm/dropSeqPipe. In brief, we trimmed adaptor sequences using Cutadapt version 1.16, performed spliced alignment of trimmed reads using STAR aligner version 2.6.1b, identified core barcodes using the whitelist function in umi_tools version 0.5.5, and used Dropseq-tools to correct barcodes and extract digital count matrices. A frozen version of the pipeline used to process the Seq-Well data is available at the Github repository https://github.com/dylkot/dropSeqPipe-dak.

93Sequencing reads were aligned to a hybrid genome/genebuild of Macaca mulatta (genome assembly Mmul_8.01, Ensembl gene build 92) and EBOV/Kikwit (GenBank accession KU182905.1).

94The scRNA-Seq data was preprocessed, clustered, and visualized using Scanpy. We removed cells with < 300 genes detected, > 10% of their UMIs derived from mitochondrial genes, or > 95% of UMIs mapped to non-genic regions. We excluded ribosomal genes, genes correlated with the percentage of UMIs assigned to mitochondrial genes (Pearson R > 0.1), and HBB as these were largely driven by the technical covariate of whether cells had loaded into Seq-Well arrays fresh or had undergone a freeze-thaw cycle with cryoprotectant. We also excluded EBOV genes and cell-cycle genes (defined by correlation with TOP2A, Pearson R > 0.1) prior to clustering so that these signals would not influence identification of cell types.

95We performed multiple iterations of clustering to detect and exclude doublets and to identify distinct cell populations at multiple levels of granularity. In each iteration, we filtered genes detected in < 10 of the cells being clustered. We then transformed raw UMI counts by normalizing the sum of counts of each cell to 10,000 (TP10K), adding 1 to each expression value, and taking the natural logarithm. In each clustering iteration, we identified and subsetted the data to highly variable genes using the highly_variable_genes function in Scanpy with the default parameters. We Z-normalized each gene and set transformed values exceeding 10 to 10. Z-normalized data was used as input to principal component analysis. We determined the number of principal components to use for downstream analysis by identifying an elbow on the Skree plot of the eigenvalues associated with each principal component. For the in vivo Seq-Well data, we used the Harmony algorithm to remove variation due to whether a PBMC sample had been processed fresh, or following a freeze-thaw. No Harmony adjustment was used for the ex vivo EBOV dataset. The Harmony-adjusted or raw principal components were then used to construct a nearest neighbor graph with the number of neighbors set to the maximum of 30 or 0.001 x the number of cells. Lastly, we clustered cells using the Leiden community detection algorithm.

96We annotated broad PBMC clusters in the in vivo and ex vivo EBOV datasets based on the following marker genes: CD8+ T cells (CD3D, GZMB, GNLY), CD4+ T cells (CD3D, IL7R), B cells (MS4A1, IGHM), Monocytes (CFD, LYZ), cDCs (FLT3, IRF8), pDCs (IRF8, GZMB), Neutrophils (CD177, LCN2), Platelets (PF4, CAVIN2), Plasmablasts (MZB1, IGHM, IGHA), and spike-in control MDCK cells (COL5A2, SLC20A1).

97For the in vivo EBOV dataset, the first clustering iteration was used to identify and filter a cluster of multiplets (expressing high levels of B cell, T cell, Neutrophil, and Monocyte genes) and MDCK control cells. A second clustering iteration was run on the filtered data to identify broad cell type clusters of T/NK-cells, B cells, and myeloid cells (Monocytes, cDCs, pDCs, and neutrophils). Sub-clustering was performed on each broad cell type in two iterations, the first to identify remaining doublets to exclude, and the second to cluster cells into final cell-types and sub-types based on annotation of marker genes (Figure S2).

98An analogous sequence of clustering iterations was used for the ex vivo EBOV dataset. However, as there were no MDCK cells spiked in, we proceeded straight to sub-clustering of the T/NK, B, Monocyte/DC, and multiplet populations from the initial clustering iteration.

99Doublets and multiplets identified during any clustering iteration were excluded. Then visualization in 2 dimensions was accomplished by computing the nearest neighbor graph with 0.001 x the number of cells as nearest neighbors, followed by Uniform Manifold Approximation and Projection.

100Gene expression values were smoothed to facilitate direct visualization of CD14, FCGR3 (which codes for CD16), and MKI67 (which codes for Ki67) by running MAGIC on log TP10K expression values with the ‘cosine’ distance metric and 3 diffusion steps.

101We performed all differential expression tests using MAST on log(TP10K + 1) normalized data. For all differential expression tests, we included (1) the percentage of mitochondrial reads and (2) the number of genes detected in a cell as covariates. For tests in the in vivo dataset, we additionally included a binary indicator covariate of whether the cell was derived from a sample that had been processed fresh, or had undergone freeze-thaw. For tests in the ex vivo dataset, we additionally included a binary indicator covariate of which NHP donor the sample was derived from.

102For viral load comparisons, we used the log10 viral load as a continuous exogenous variable and only considered cells with ≥ 1 viral read. For all other comparisons, we used a binary exogenous variable indicating the reference and the query group. “Viral load” and “bystander vs. EBOV infected cells” in the in vivo dataset were conducted only considering cells from DPI 5–8. “Viral load” comparisons in the ex vivo dataset only considered cells from the 24 HPI time point.

103Differential expression p-values were corrected for multiple hypothesis testing using the method of Benjamini and Hochberg.

104We clustered the log fold-change profiles of 1,437 differentially expressed genes (rows of Figure 3A) considering B, CD4+ T, CD8+ T, NK, and monocyte populations at the three EVD stages. Prior to clustering, insignificant values (p > 0.2) were first set to 0 and genes were normalized to unit L2 norm. We performed k-means clustering with K = 11 and the default parameters in Scikit-learn version 0.22.2.post1. Varying K above and below 11 led to highly consistent results with splitting or merging individual clusters at the margin, and we picked K = 11 as the lowest value that yielded the “B down” module.

105Our infection detection model assumes that EBOV reads assigned to a cell are either due to true viral RNAs inside that cell, or due to ‘ambient’ extracellular EBOV RNA in the sample that, by chance, are captured in a well of the Seq-Well array along with the cell. Our null model for how many EBOV reads would be expected in a cell by chance therefore depends on two main parameters which we estimate from the data: (1) what proportion of ambient RNAs in a sample are due to EBOV, and (2) what proportion of a cell’s expression profile is due to ambient RNA. Our method therefore precedes through the following steps which we describe below:Estimate an ambient RNA profile for each Seq-Well arrayEstimate the proportion of each cell’s transcripts that are due to ambient RNADetermine if there are more EBOV transcripts in each cell than would be expected based on the proportion of EBOV in the ambient RNA and the cell’s estimated level of ambient RNA contamination.

106Estimate an ambient RNA profile for each Seq-Well array

107Estimate the proportion of each cell’s transcripts that are due to ambient RNA

108Determine if there are more EBOV transcripts in each cell than would be expected based on the proportion of EBOV in the ambient RNA and the cell’s estimated level of ambient RNA contamination.

109We assume that cell barcodes with few transcripts detected correspond to wells in the Seq-Well array that lack a cell, and therefore, that any transcripts assigned to those cell barcodes derive from ambient RNA. We consider all cell barcodes with fewer than 50 transcripts detected to be empty wells and compute an ambient RNA profile of the proportion of transcripts from these barcodes assigned to each gene. This is analogous to the approach used in several ambient RNA correction approaches such as.

110We denote the number of genes in the dataset as G and define the ambient RNA profile for a given Seq-Well array, a, as a G-dimensional vector. denotes the estimated proportion of ambient RNAs that are assigned to gene g based on transcripts from all cell barcodes with fewer than 50 UMIs. Since is a proportion, the following constraints hold:

111We adapt the previously published Consensus Non-negative Matrix Factorization (cNMF) method to estimate the ambient RNA contamination level of each cell. In brief, cNMF learns a user-specified number of gene expression programs (GEPs), each a non-negative G-dimensional vector representing the average expression profile of an individual cell type or cellular activity (e.g., cell-cycle or interferon response) that are present in the data. In addition, it learns a “Usage” matrix reflecting the % contribution of each GEP in each cell (i.e., the % of each cell’s transcripts derived from each GEP). For the following notation, matrices are indicated in bold. Denoting the number of cells in the dataset as C and the user-specified number of GEPs as K, cNMF is given an input CxG matrix of transcript counts and returns a non-negative KxG matrix of GEPs reflecting the relative contribution of gene g in GEP k, as well as a non-negative CxK matrix reflecting the percentage of transcripts in cell c that are due to GEP k.

112We adapt this approach to return an updated, Cx(K+1) dimensional usage matrix that includes an additional column reflecting the usage of the ambient RNA profile. We run cNMF as published to obtain the GEP matrix. We then append the ambient RNA profile for a given array as an additional row to and normalize each row to sum to 1 like so:where denotes the L1 norm, a tilde is used to denote an L1 normalized vector, and denotes the ith row of. We then run a final iteration of NMF with the GEP matrix fixed to. This jointly estimates a usage of the ambient RNA profile and the other GEPs using non-negative least-squares. For a given cell c from an array a, this amounts to solving the following non-negative least-squares optimization:We combine these coefficients for all cells and programs into a single matrix and L1 normalize the usages to sum to 1:where denotes the number of cells derived from array a. The “K+1”th column of reflects the estimated proportion of transcripts due to ambient RNA in all of the cells from array a. We repeat this calculation for each array separately and denote the estimated contribution of ambient RNA for a cell c as.

113For each cell, we determine a threshold,, for the number of EBOV transcripts required to call that cell infected while keeping the false positive rate below f, a user specified threshold (f = 0.01 in all of our analyses). We first calculate, the proportion of ambient RNA transcripts in array a that are due to EBOV as follows:Then, for a cell with transcripts total; a given proportion,, of its reads derived from ambient RNA;, the proportion of ambient RNA contaminating reads expected to map to EBOV, we compute using binomial statistics as follows:Where F is the inverse survival function of the Binomial distribution with event probability p and N trials. We identify cells with reads as infected.

114We estimate sensitivity to call an infected cell with either 1% or 0.1% of its reads due to EBOV across a range of 60 false positive rate thresholds (f). We first randomly sampled 2,000 cells from the live EBOV treatment samples in the ex vivo data, or the non-baseline samples from the in vivo data, to serve as an empirical distribution for and. For each cell and specificity threshold f, we then calculate the probability of correctly calling a true positive cell as positive. We model the distribution of the number of EBOV reads in a cell as the convolution of 2 binomial distributions: is the binomially distributed number of true cell-derived reads mapping to EBOV and is the analogous distribution for ambient RNA-derived reads mapping to EBOV. is 0.01 or 0.001 by assumption, and is estimated empirically for each cell as. We calculate the expected values for and and round up to the nearest integer, as follows:We then directly calculate the convolution of the two binomial distributions for each cell:We then calculate the sensitivity for each cell as the probability that the convolution distribution is greater than the empirical threshold:We plot the average sensitivity across the 2000 randomly sampled cells as a function of f for equal to either 0.01 or 0.001.

115We downloaded gene sets from the Molecular Signatures Database version 6.2 for gene set enrichment testing. We considered all Hallmark or C2 gene sets containing greater than 10 genes that were present in our expression data. We tested expression modules for enrichment using Fisher’s exact test and corrected for multiple hypothesis testing using the method of Benjamini and Hochberg.

116To test continuous expression profiles for gene set enrichment (Figure S4B), we used the rank-sum test comparing genes in the gene set to all genes not in the set.

117We identified 58 genes in the “Global up” module that were also included in one or more of the following gene sets from the molecular signatures database: HECKER_IFNB1_TARGETS, BROWNE_INTERFERON_RESPONSIVE_GENES, MOSERLE_IFNA_RESPONSE, HALLMARK_INTERFERON_ALPHA_RESPONSE, HALLMARK_INTERFERON_GAMMA_RESPONSE (Table S3). We then scored cells for the average expression of these genes using the score_genes function in Scanpy with 58 control genes, as this was the number of genes in the ISG set, and otherwise default parameters.

118We computed a macrophage score based on the set of 618 genes annotated as significantly up or downregulated during in vitro monocyte-to-macrophage differentiation (Table S4). We computed each cell’s macrophage score as the dot-product of its expression profile for the 618 genes (in log TP10K) with the log fold-change reported for each gene in. This effectively weights genes by both the direction and magnitude of their change during in vitro macrophage differentiation.

119We obtained all of the human PBMC datasets produced using v3 or v3.1 chemistry from the 10X website (Key Resources Table), aggregated them together, and processed the resulting dataset using the same pipeline as the NHP Seq-Well data. Briefly, we first filtered out genes detected in fewer than 10 cells before converting to log TP10K and performed PCA as described above. Then, we used Harmony to integrate out variation due to the different samples of origin and used 30 nearest neighbors for Leiden community detection and UMAP dimensionality reduction (Figure S5E). We did not perform any sub-clustering on this dataset.

120We obtained Human Cell Atlas bone marrow data from the Human Cell Atlas data portal and processed it according to the same pipeline as the NHP data with a few modifications. We filtered doublets prior to clustering by running Scrublet separately within each of 8 donor batches with an expected doublet rate parameter of 6%. We identified and excluded cell-cycle associated genes, as those with a Pearson correlation > 0.3 with TOP2A. We integrated data from the different donor batches using Harmony and used 30 nearest neighbors for Leiden community detection and UMAP dimensionality reduction. We performed 3 rounds of sub-clustering: First we clustered all of the cells to identify monocyte and dendritic lineage cells (Figure S5F). Second, we clustered just hematopoietic stem cells (HSCs) and monocyte/dendritic progenitor cells to identify doublets (as those falling into a cluster characterized by T cell marker genes such as CD3D and CD3E). Finally, we re-clustered this set with the doublets excluded to identify monocyte lineage cells, plasmacytoid dendritic cells, and conventional dendritic cells (Figure S5G).

121We confirmed our marker gene-based annotations of the myeloid cell populations by comparing these cells to the circulating human PBMC dataset. We identified the nearest neighbor of each bone marrow myeloid progenitor cell in the PBMC dataset based on Euclidean distance of TP10K-normalized cells, considering overdispersed genes identified in the PBMC dataset based on the V-score (baseline-corrected Fano factor). We then visualized the nearest-PBMC assignment of the bone marrow myeloid cells on a UMAP embedding (Figure S5H).

122Finally, we combined the monocytes and monocyte precursor cells from the human PBMC and bone marrow datasets into a single reference. We again normalized the data to log TP10K and computed UMAP embeddings following the same procedure as for the individual datasets (Figure S5I), using Harmony to remove variation due to donor sample. We then down-sampled this data so that there would be equivalent numbers of cells of each of the bone marrow and PBMC clusters (i.e., 982 cells per cluster as that was the number of cells in the smallest cluster). We identified the nearest neighbor of each NHP monocyte in the down-sampled reference dataset, as described above and computed the percentage of NHP monocytes assigned to CD14+ or CD16+ clusters from either human bone marrow or PBMC (Figure 5F).

123We used Statsmodels version 0.11.1 to fit a logistic regression predicting EBOV infection status among all monocytes from late EVD, based on the following features: macrophage score, MAGIC smoothed values of CD14, and CD16 (FCGR3), and an interaction term for the product of the MAGIC smoothed values for CD14 and CD16.

124Our clustering and dimensionality analysis of the NHP CyTOF data was analogous to the Seq-Well pipeline with a few adaptations. We down-sampled a total of 1.1 million cells, consisting of 300,000 baseline cells and 100,000 cells from each DPI, selecting a uniform number of cells from each sample at a given DPI. We used the Arcsinh transformation of CyTOF raw intensity values divided by 5, which is standard in the field. We set a ceiling of transformed intensity values for each gene at the 99.999th percentile to reduce the effect of very rare outliers. We mean-centered the data but did not variance normalize prior to PCA. We then performed multiple clustering iterations of the data using the Leiden algorithm with the number of nearest neighbors set to the maximum of 30 or 0.01% of the number of cells in the dataset.

125We annotated broad PBMC clusters in the CyTOF datasets based on the following marker proteins: CD8+ T cells (CD3, CD8), CD4+ T cells (CD3, CD4), NK cells (CD16, CD161), B cells (CD19, IgM), Monocytes (CD11b, BDCA3, CD14, CD16, HLA-DR), cDCs (HLA-DR, CD11c, CD1c), pDCs (HLA-DR, CD4, CD123), Neutrophils (CD11b, CD66), Platelets (CD61, BDCA3), Plasmablasts (IgM high, but little or no CD19), Basophils (CD11b, CD123, low HLA-DR). There was also a cluster of cells characterized by high HLA-DR, Ki67, and CD38 which we annotate as “Unassigned APC.”

126In the first clustering iteration, we identified and excluded clusters of doublets as those expressing markers of two or more broad cell-types. We also excluded a cluster characterized by high expression of CD3 but low expression of both CD4 and CD8 that we speculated was a technical artifact as it was nearly exclusive to a single CyTOF batch.

127Next, to address batch effects between the different CyTOF runs, we grouped cells into the broad categories of Monocytes and DCs, CD4+ T, CD8+ T, NK, Neutrophils, B, Plasmablasts, and Platelets, and used COMBAT to adjust for batch effect separately within each broad category, using the default COMBAT parameters in the Scanpy implementation.

128We then ran cNMF to identify and regress out artifact signals in the data. We ran cNMF using all 42 markers as inputs and adapted the method to not perform any variance scaling prior to NMF; this is because CyTOF intensities of different markers are already expressed on a comparable scale, unlike genes in single-cell RNA-Seq data which can vary over different orders of magnitude. We also set a floor on the input data to be non-negative (because after COMBAT, a small percentage of the values were slightly below 0). We selected K = 7 as the cNMF dimensionality as the solution stability fell off dramatically at higher values. One gene expression program (GEP) was characterized by high levels of the platelet markers CD61 and BDCA3. It was mostly enriched in the platelet cluster but was also elevated in a subset of cells of all major cell types, which suggests that it reflects an artifact of platelets sticking to cells. A second GEP was characterized by high levels of all of the intracellular markers and it was also distributed throughout cells of multiple clusters. We interpreted this GEP as a cell permeabilization artifact reflecting the relative accessibility of a cell’s intracellular proteins to CyTOF antibody staining. We regressed these 2 GEPs out of the data by subtracting the matrix (outer) product of the Usage and GEP matrices. We denote the number of cells as C, the number of genes as G (42 in our data), and the number of programs selected as K (7 in our data). The CxG input data matrix is denoted as X, the KxG GEP matrix returned by cNMF as and the CxK usage matrix as. Then the correction is as follows:Where and are the Cx1 dimensional matrices representing the usage of the platelet and permeabilization GEPs and and are the 1xG dimensional matrices representing the spectra of the platelet and permeabilization GEPs, and multiplication is the outer product. We also filtered cells assigned to the platelet cluster in the initial clustering, since their predominant signal had been regressed out.

129We then repeated clustering of the corrected data to identify broad cell types (Neutrophils, Monocytes, and DCs), followed by sub-clustering within each broad cluster to generate the sub-clusterings in Figure S2. The data was visualized using the UMAP algorithm as described for the scRNA-Seq data.

130The human PBMC CyTOF data was processed with the same pipeline as the NHP data with a few modifications. The data were down-sampled to 280,000 cells total (20,000 per sample), batch correction was performed using Harmony, and 0.001 x the number of cells was used for K nearest neighbor graph construction. No cNMF or COMBAT adjustment steps were performed.

131Details of statistical testing, sample size, center, and dispersion can be found in the figure legends, the main text, and STAR Methods.

132This study did not generate new additional resources (website, forum, clinical trial).