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

Published:November 06, 2020 Cell

Highlights

• Interferon response is suppressed in infected cells but activated in bystander cells

• EBOV represses antiviral genes and upregulates pro-viral genes in infected cells

• Proliferative CD14– CD16– monocyte precursors expand in circulation during EVD

• Identification of expression markers of EBOV tropism for circulating cells in vivo

Summary

Ebola virus (EBOV) causes epidemics with high mortality yet remains  understudied due to the challenge of experimentation in high-containment and outbreak settings. Here, we used single-cell transcriptomics and  CyTOF-based single-cell protein quantification to characterize  peripheral immune cells during EBOV infection in rhesus monkeys. We  obtained 100,000 transcriptomes and 15,000,000 protein profiles, finding that immature, proliferative monocyte-lineage cells with reduced  antigen-presentation capacity replace conventional monocyte subsets发现抗原呈递能力降低的未成熟增殖单核细胞谱系细胞取代了常规单核细胞亚群,  while lymphocytes upregulate apoptosis genes and decline in abundance 而淋巴细胞则上调了凋亡基因并降低了丰度.  By quantifying intracellular viral RNA, we identify molecular  determinants of tropism among circulating immune cells and examine  temporal dynamics in viral and host gene expression 研究了病毒和宿主基因表达的时间动态. Within infected  cells, EBOV downregulates STAT1 mRNA and interferon signaling, and it upregulates putative pro-viral genes (e.g., DYNLL1 and HSPA5), nominating pathways the virus manipulates for its replication. This  study sheds light on EBOV tropism, replication dynamics, and elicited  immune response and provides a framework for characterizing host-virus  interactions under maximum containment 为表征宿主-病毒在最大遏制下的相互作用提供了框架.

Results

Figure 1Study Design

Under BSL-4 containment, we collected blood samples from a total of 21 rhesus monkeys at multiple days post-EBOV inoculation, extracted peripheral blood mononuclear cells (PBMCs), and profiled single-cell transcriptomes and 42 protein markers using Seq-Well and CyTOF. Seq-Well quantifies both host (black) and viral (red) RNA expression, allowing comparisons between infected and bystander cells. Daily clinical parameters (body temperature, clinical signs, and body weight) were also collected for each animal, and complete blood counts were obtained for each blood draw. See also Figure S1A and Table S1.

Figure 2 Changing Cell-Type Abundance, Proliferation Rate, and Infection Status during EVD

(A) Time course of viral load (red, left y axis, log10 scale) and clinical score (blue, right y axis). Markers: mean; error bars: minimum and maximum; LOD, limit of detection by reverse transciption quantitative PCR.

(B and C) Uniform Manifold Approximation and Projection (UMAP) embedding of Seq-Well (B) and CyTOF (C) data, colored by annotated cluster assignment. See also Figures S1 and S2, and Data S1.

(D) Fold change (log2 scale) in the absolute abundance (cells/μL of whole blood) of each cell type relative to baseline based on CyTOF clusters. Error bars: mean ± 1 SE. See also Figures S3A and S3B.

(E and F) UMAP embedding of Seq-Well (E) and CyTOF (F) data, colored by the day post-infection (DPI) on which each cell was sampled.

(G) Percentage of Ki67-positive cells (CyTOF intensity >1.8) of each cell type. Error bars: mean ± 1 SE. See also Figures S3C and S3D.

(H) UMAP embedding of Seq-Well data, colored by the percentage of cellular transcripts mapping to EBOV.

(I) Percentage of infected cells by cell type based on Seq-Well. Dashed line: 1% false positive rate threshold for calling infected cells. Error bars: 95% CI on the mean based on 1,000 bootstraps. See also Figures S1H.

Figure S1 Cell-Type Markers for Seq-Well and CyTOF Clusters, Related to Figure 2

(A) Overview of study cohorts and blood draw timelines. Animals were grouped into cohorts with pre-scheduled necropsy times (at baseline, or day post infection [DPI] 3, 4, 5, 6 - n = 3 each), or allowed to progress until clinical score exceeded 10 (terminal), predetermined euthanasia criteria. Dots: scheduled blood draws for each cohort; red: intermediate (non-necropsy) draw; gray: draw that coincided with euthanasia and necropsy. Necropsy and baseline normal draws were used for Seq-Well and CyTOF, while intermediate post-infection draws were available only for CyTOF.

(B) Expression profiles of cell-type marker genes (columns) for cell-type clusters (rows) based on the in vivo Seq-Well data. Circle area represents the percentage of cells in each group in which the gene was detected, and color denotes the average expression level (loge TP10K).

(C) Average expression (Z-normalized CyTOF intensity) profiles of cell-type marker genes (columns), for cell-type clusters (rows), based on the CyTOF data.

(D) Uniform Manifold Approximation and Projection (UMAP) embedding of post-integration Seq-Well data, colored by the sample source (NHP, DPI, and whether the sample was loaded for Seq-Well without any freezing [.fresh] or was frozen with cryoprotectant and thawed prior to Seq-Well [.FRZ]). A maximum of 500 cells per sample is plotted to increase representation across samples.

(E) UMAP embedding of Seq-Well data, colored by whether cells were processed fresh (orange) or after freeze/thaw (blue) prior to Seq-Well.

(F) UMAP embedding of Seq-Well data, colored by depletion of abundant sequences by hybridization (DASH) treatment. We developed a DASH-based method to remove a PCR adaptor artifact from some Seq-Well sequencing libraries (STAR Methods), and performed this 0 times (No DASH, blue), 1 time (DASH, orange), or 2 times sequentially (DASHx2, red). For a few samples, we sequenced ‘No DASH’ and ‘DASH’ libraries and merged the reads (mixed, green).

(G) UMAP embedding of batch-corrected CyTOF data, colored by the multiplex batch in which it was pooled and analyzed by CyTOF.

(H) Receiver operating characteristic curves for identifying EBOV-infected cells. Estimates of sensitivity to detect an infected cell at various false positive rate thresholds in vivo (left) and ex vivo (right). Curves are estimated separately for a hypothetical viral load of 0.1% (blue line) and 1% (orange line).

Figure 3 Patterns of Differential Expression across EVD Stages and Cell Types

(A) Fold changes (loge scale) of 1,430 differentially expressed genes (rows) in each cell type at early (E), middle (M), and late (L) EVD (columns), relative to baseline, with insignificant values (p > 0.2) set to 0. Genes were grouped into modules through unsupervised k-means clustering. See also Tables S2 and S3.

(B) Same as (A) but displaying the average loge fold change of each module.

(C) Distribution of interferon-stimulated gene (ISG) scores for each cell type. White markers: median; bars: interquartile range. See also Figures S4A and S4B.

(D) Differential expression of monocytes in late EVD compared to baseline.

Figure S4 Quantification of Cytokine Expression and Enrichment of Response Signatures, Related to Figures 3 and 4

(A) Average expression values (loge TP10K) of literature-annotated cytokines (columns) across cell types and stages of acute EVD (rows). Values are plotted as a ratio relative to the maximum across cell types and stages. Values that are statistically different from baseline (p < 0.05) are indicated with a blue star.

(B) Heatmap of rank-sum test statistics for comparison of differential expression log fold-changes of genes in a gene set (rows) compared to genes not in the set. The log fold-changes were defined from differential expression profiles of each cell type at each EVD stage (columns) relative to baseline. Five gene sets were tested — three from the Hallmark database (IFN ALPHA, IFN GAMMA, and TNF ALPHA VIA NFKB) (

Liberzon et al., 2015) and 2 constructed from the hallmark sets, as uniquely IFNα-regulated genes in “IFN ALPHA” but not “IFN GAMMA” (“IFN ALPHA - GAMMA”), and vice versa for uniquely IFNγ-regulated (“IFN GAMMA - ALPHA”). See also Table S3.

(C) Fold change (log2 scale) in average HLA-DR CyTOF intensity on B cells at each DPI relative to baseline for each PBMC sample. Colored lines connect serial samples from the same NHP.

Figure 4 Monocytes Dramatically Reduce Expression of MHC Class II Proteins Independent of Infection Status

(A) Expression of major histocompatibility (MHC) or MHC-associated genes (rows) in key cell types at baseline (B), early (E), middle (M), or late (L) EVD (columns). Circle size: percentage of cells in that group in which the gene was detected; color: mean expression in Z score normalized, loge transcripts per 10,000 (TP10K). The “MAMU-” prefix, which designates MHC genes in rhesus monkeys, was removed; the “HLA-” prefix is indicated by “(H).”

(B) CyTOF intensity of HLA-DR protein in antigen-presenting cells. Boxes: median and interquartile range; whiskers: 2.5th and 97.5th percentiles. Colored stars indicate significant decreases from baseline (rank-sum test p < 0.05) with color corresponding to stage.

(C and D) Fold change (log2 scale) in average CD38 (C) and HLA-DR (D) CyTOF intensity on monocytes at each DPI relative to baseline, connected by colored lines for each NHP. See also Figure S4C and Data S1.

(E) Average gene expression (loge TP10K) for four MHC class II genes in monocytes, stratified by cell-infection status. Error bars: 95% CI on the mean based on 200 bootstraps.

Figure 5 ISG Suppression, Co-expression of CD14 and CD16, and Expression of Macrophage Genes Are Associated with Monocyte Infectivity

(A) Differential expression between infected and bystander monocytes from DPI 5–8. Genes are colored by membership in sets of genes (Mac. Up/Down = up- or downregulated during in vitro differentiation of monocytes into macrophages). See also Table S4.

(B) UMAP embedding of monocyte gene expression data, colored by (left-to-right) DPI, CD16 expression (loge TP10K), CD14 expression (loge TP10K), and percentage of cellular transcripts mapping to EBOV.

(C) Smoothed expression (loge TP10K) of CD14 and CD16 for monocytes during EVD. Boxes: CD14+, CD16+, DN, and DP subsets described in the text; numbers: percentage of cells in each subset at that EVD stage. See also Figures S5A and S5B.

(D) CD14 and CD16 protein expression (CyTOF intensity) on monocytes at each DPI. Bivariate kernel density plot with 200 randomly sampled cells is overlaid as a scatterplot. See also Figure S5C.

(E) CD14 and CD16 protein expression (CyTOF intensity) on monocytes in a case of human EVD, colored by Ki67 protein expression for multiple days after symptom onset. See also Figure S5D.

(F) Percentage of assignment of NHP CD14/CD16 subsets at each EVD stage to human myeloid reference populations (BM-MP: bone marrow monocyte progenitors, PBMC-CD16+: circulating CD16+ monocytes, PBMC-CD14+: circulating CD14+ monocytes). See also Figures S5E–S5K.

(G) Percentage of infected monocytes in each CD14/CD16 subset in late EVD. Error bars: 95% CI on the mean based on 1,000 bootstraps.

(H) Association between macrophage score (x axis) and percentage of infected cells (left y axis, red) and expression of the differentiation marker NR1H3 (right y axis, blue, loge TP10K). We ordered monocytes from late EVD by macrophage score, and averaged percentage of infected cells and NR1H3 expression within 400-cell sliding windows. See also Figures S6A–S6C.

(I) MX1 expression (loge TP10K) in monocytes at baseline, and uninfected bystanders or infected cells in late infection. Boxes: median and interquartile range; whiskers: 2.5th and 97.5th percentiles. Statistical significance was assessed by rank-sum test. See also Figure S6D.

(J) Scatterplot of ISG score (y axis) versus percentage of cellular transcripts mapping to EBOV (x axis) for infected monocytes in late EVD (DPI 6–8). Statistical significance was assessed by Spearman ρ.

Figure S5 Extended Characterization of Interferon and Double-Negative CD14– CD16– Monocytes, Related to Figure 5

(A) Clustermap of pairwise Pearson correlations between cell type clusters at baseline and late EVD. Correlations are computed on average loge TP10K expression values of overdispersed genes. DN and DP monocytes at late EVD are more similar to monocytes (including baseline CD14+s) than other cell types.

(B) Scatterplot of MAGIC-smoothed expression values (loge TP10K) of CD14 and CD16 for monocytes in baseline, early, mid, and late disease stages. Cells are colored by smoothed expression levels of MKI67 (the gene coding for Ki67 protein). Boxes: CD14+, CD16+, DN, and DP subsets described in the text; numbers: percentage of cells falling into each subset.

(C) Scatterplot of protein expression (CyTOF intensity) of CD14 and CD16 for 1,000 randomly sampled monocytes at each DPI. Cells are colored by Ki67 expression. Boxes: CD14+, CD16+, DN, and DP subsets described in the text; numbers: percentage of cells falling into each subset.

(D) Scatterplot of protein expression (CyTOF intensity) of CD14 and CD16 for monocytes during human EVD. Left: monocytes from healthy human controls. Right: monocytes from 3 EVD cases (S1, S2, and S3) at various days post symptom onset. Cells are colored by Ki67 marker intensity. Boxes: CD14+, CD16+, DN, and DP subsets described in the text; numbers: percentage of cells falling into each subset.

(E) UMAP embedding of healthy human PBMCs dataset, colored by annotated cluster assignment, based on known marker genes. (Plasma.: Plasmablast).

(F) UMAP embedding of healthy bone marrow cells, colored by cluster assignment, based on marker genes. (HSC: hematopoietic stem cell, Plasma.: Plasmablast, Megakar.: Megakaryocyte, Mono/DC: monocyte and dendritic cell, BM-Macro: bone marrow macrophage).

(G) UMAP embedding of sub-clustered HSC and monocyte/dendritic lineage cells. (BM: bone marrow, MP: monocyte progenitor)

(H) Same UMAP embedding as Figure S5G, but colored by the cluster identity of their nearest neighbor in the human PBMC dataset (Figure S5E).

(I) UMAP embedding of the merged reference dataset of healthy bone marrow HSCs and monocyte lineage cells and PBMCs. Left sub-panel is colored by cluster assignment. Right sub-panels are colored by marker gene expression (loge TP10K).

(J) Expression profiles of selected genes for human bone marrow monocyte progenitors (BM-MPs) and human circulating monocytes (PBMC-Monos). Circle area: percentage of cells in which the gene was detected; color: average expression (Z-normalized loge TP10K).

(K) Expression profiles of selected genes for NHP monocyte subsets at baseline or late EVD for orthologs of the genes in (J). Circle area: percentage of cells in which the gene was detected; color: average expression level (Z-normalized loge TP10K). CD34 is grayed out because it is detected in <10 cells.

Figure 6 Viral Transcriptional Dynamics of Infected Monocytes In Vivo and Ex Vivo

(A) Schematic of EBOV challenge of PBMCs ex vivo. See also Figure S7.

(B and C) Percentage of cellular transcripts derived from EBOV  (intracellular viral load) in monocytes from PBMCs inoculated with live  virus ex vivo (B) or from PBMCs of NHPs infected in vivo (C). See also Figures S8A–S8D.

(D) Schematic of EBOV transcription. The viral RNA-directed RNA-polymerase  transcribes each gene sequentially but occasionally releases the genomic RNA template, ending transcription. As a result, transcription  frequency decreases from NP to L.

(E and F) Proportion of each EBOV gene versus viral load (log10 scale), ex vivo (E) or in vivo (F). We ordered infected monocytes by viral load and averaged the  percentage of each viral gene over 50-cell sliding windows. Bands: mean ± 1 SD. See also Figures S8E and S8F.

Methods

Scoring cells for interferon response and macrophage differentiation

We 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 (

Satija et al., 2015) with 58 control genes, as this was the number of genes in the ISG set, and otherwise default parameters.

We computed a macrophage score based on the set of 618 genes annotated as significantly up or downregulated during in vitro monocyte-to-macrophage differentiation (Dong et al., 2013) (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 (Dong et al., 2013). This effectively weights genes by both the direction and magnitude of their change during in vitro macrophage differentiation.

Comparison of EVD monocyte subsets with human bone marrow and PBMC data

We 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 (

Korsunsky et al., 2019) to integrate out variation due to the different samples of origin and used 30 nearest neighbors for Leiden community detection (

Traag et al., 2019

) and UMAP dimensionality reduction (

Becht et al., 2018

) (Figure S5E). We did not perform any sub-clustering on this dataset.

We obtained Human Cell Atlas bone marrow data from the Human Cell Atlas data portal (Hay et al., 2018) 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 (Wolock et al., 2019) 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).

We 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) (Klein et al., 2015). We then visualized the nearest-PBMC assignment of the bone marrow myeloid cells on a UMAP embedding (Figure S5H).

Finally, 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).