Historical RNA samples collected from our previous study were used in this study (Marzi et al., 2018). In this previous study, four cohorts of rhesus macaques were infected intramuscularly with 1 × 10³ focus forming units (FFUs) of EBOV-Mayinga (n = 5) or EBOV-Makona isolates: Guinea C07, Mali, or Liberia (n = 3 each). Whole blood (WB) samples were collected at 0, 2, 4, and 6 days post infection (DPI), and RNA was isolated from the WB of EBOV-Mayinga- and EBOV-Makona-infected NHPs using the QIAmp Viral RNA Kit (Qiagen) as previously described (Marzi et al., 2018). This prior work was performed in the maximum containment laboratory at the Rocky Mountain Laboratories (RML), Division of Intramural Research, National Institute of Allergy and Infectious Diseases, National Institutes of Health. RML is an AAALAC accredited institution. All procedures followed standard operating procedures (SOPs) approved by the RML Institutional Biosafety Committee (IBC). Animal work was performed in strict accordance with the recommendations described in the Guide for the Care and Use of Laboratory Animals of the National Institute of Health, the Office of Animal Welfare and the Animal Welfare Act, United States Department of Agriculture. The study was approved by the RML Animal Care and Use Committee (ACUC). Procedures were conducted in animals anesthetized by trained personnel under the supervision of veterinary staff. The humane endpoint criteria for euthanasia were specified and approved by the RML ACUC. All efforts were made to ameliorate animal welfare and minimize animal suffering in accordance with the Weatherall report on the use of NHPs in research¹.

RNA samples from DPI 0, 4, and 6 were used in the current study. Integrity and concentration of historical RNA samples were validated on the Agilent 2100 Bioanalyzer prior to library construction. rRNA was depleted from the samples using the NEBNext rRNA Depletion kit before cDNA library construction with the NEBNext Ultra II Directional RNA Library Prep Kit (Illumina). cDNA library quality and concentration were confirmed on the Agilent 2100 Bioanalyzer before sequencing on the Illumina HiSeq2500 or NovaSeq platforms.

Raw sequences were trimmed to a minimum length of 75 bp and average Phred score of 30 using Trim Galore before alignment to the Macaca mulatta genome “Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa” using tophat. Genes were annotated with the Ensembl annotation files for M. mulatta (Macaca_mulatta.Mmul_8.0.1.97.gtf). Preliminary processing of RNA-Seq data was performed using the systemPipeR package available from Bioconductor (Tw and Girke, 2016).

Three approaches were used to identify differentially expressed genes (DEGs): (1) EdgeR, (2) STEM, and (3) MaSigPro (Ernst and Bar-Joseph, 2006; Robinson et al., 2010; Nueda et al., 2014). EdgeR package, which uses the Trimmed means of M-values (TMM), was used to identify DEGs, and counts were normalized using the reads per kilobase of exon per million (rpkm) (Mortazavi et al., 2008; Menicucci et al., 2017; Versteeg et al., 2019). DEGs were then filtered for only those encoding human protein-coding homologs with an average rpkm ≥ 5, FDR ≤ 0.05, and fold change ≤ −1 or ≥ 1. To identify longitudinal patterns of gene expression changes throughout infection, we used Short Time Series Expression Miner (STEM) software, which clusters genes both by significance and expression patterns (Ernst and Bar-Joseph, 2006). MaSigPro was used for regression analysis with a two-way forward regression strategy comparing the EBOV-Makona Mali, Liberia, and EBOV-Mayinga isolates to EBOV-Makona Guinea C07 (Conesa et al., 2006). Genes that contributed significantly to the most degree of variance in at least 16 different comparisons were fit into the model and were then clustered by temporal expression pattern. Sparse partial least squares discrimination analysis (sPLS-DA) was performed with the mixOmics R package for validation and classification (Rohart et al., 2017). Models were initially built with a minimum of 10 components and validated with five-fold cross-validation and 10 repetitions.

Functional enrichment was carried out using Metascape to identify gene ontology (GO) terms representing biological processes (Zhou et al., 2019). GO network plot was created with Cytoscape visualization software (Otasek et al., 2019; Zhou et al., 2019). All heatmaps, bubbleplots, and Venn diagrams were made and analyzed through R packages (R Core Team, 2015; Wickham, 2009). Digital cell quantification (DCQ) was performed using ImmQuant, a platform that uses a deconvolution algorithm to compare transcriptomic sequencing data to reference datasets in order to predict changes in cell frequencies (Abbas et al., 2005; Novershtern et al., 2011; Frishberg et al., 2016). Both DMAP and IRIS datasets were used.

All statistical analyses were performed in GraphPad Prism (version 7). Genes with similar temporal expression profiles were determined using time course data (analyzed through STEM) and were grouped into clusters (analysis of clusters shown in Figures 2, 4). All clusters for each isolate were regenerated on one plot with multiple biological replicates of the DEGs found through STEM (Figures 2A, 4A). ImmQuant deconvolution data was analyzed using a repeated measure mixed effect model statistical design against 0 DPI, where numbers of samples were not equal between timepoints. Specifically, multiple comparison analyses were performed between subsequent days post infection (4 DPI and 6 DPI) and baseline (0 DPI) and indicated with asterisks (Figures 3, 6 and Supplementary Figure 3). For EBOV-Makona Mali and Liberia datasets, significance was additionally tested using a two-tailed T test for population changes between 0 DPI and both 4 + 6 DPI.

Sequencing data for rhesus macaques is available at BioProject PRJNA718880.

No studies to date have examined the transcriptional response to EBOV-Mayinga, the first EBOV isolate identified in 1976 (CFR of ∼90%). Therefore, we leveraged access to historical WB RNA samples from our previous study to profile host responses at 0, 4, and 6 DPI (Marzi et al., 2018). Significant transcriptional changes were identified with a total of 985 and 1,721 DEGs detected 4 and 6 DPI, respectively, with a substantial overlap between the two time points (Figure 1A). Functional enrichment of these DEGs using Metascape showed that DEGs downregulated at 4 and 6 DPI enriched to GO terms related to adaptive immunity (Figure 1B) such as “thymus development” (e.g., AGER, IL7R, and ZAP70) and “T cell activation” (e.g., AKT1, BTLA, and TCF7) (Figure 1B and Supplementary Figure 1A) (Zhou et al., 2019). DEGs downregulated at 6 DPI enriched to “antigen processing and presentation of exogenous peptide antigen via MHC class II” (e.g., HLA-DMA, -DMB, and -DOB) and “B cell proliferation” (e.g., CD19, CD180, and PLCL2) (Figure 1B and Supplementary Figure 1A). In addition, DEGs downregulated at 6 DPI played a role in “DNA repair” (e.g., ATM, DNA2, and POLE2) and “cell cycle phase transition” (e.g., BUB1B, CDC27, and CKAP5) (Figure 1B).

In contrast to downregulated DEGs, upregulated DEGs detected at 4 and 6 DPI enriched predominantly to GO terms related to innate immunity and inflammation, notably “myeloid cell differentiation” (e.g., BATF2, MAFB, NFE2 < TLR2, TLR3, and RELB) (Figure 1B and Supplementary Figure 1B). DEGs enriching to “type I interferon (IFN) production” consisted primarily of IFN-stimulated genes (ISG; e.g., IFI16, IRF7, ISG15, and STAT1), as well as genes involved in detection and response to pathogenic nucleic acid (e.g., DHX58, DDX58, IRF7, and TBK1) (Supplementary Figure 1C). Finally, some upregulated DEGs at 4 and 6 DPI also enriched to “T cell activation” (e.g., CD274, IL2RA, LFNG, and TCIRG1) (Supplementary Figure 1D).

We next sought to correlate these transcriptional changes with alterations in immune cell frequencies and activation status. Since flow cytometry was not conducted in the earlier study and no PBMC were cryopreserved, we performed DCQ to predict changes in immune cell population frequencies based on our transcriptional data (Figure 1C). Using the IRIS immune cell database, transcriptional findings predicted sharp decreases in the frequencies of T and B cells (particularly antibody-secreting B cells), while the number of myeloid cells and natural killer (NK) cell populations was predicted to increase with disease progression (Figure 1C).

Given that bivariate analysis of the transcriptional changes at each DPI relative to baseline does not consider the dynamic, longitudinal patterns of gene expression, we next used Short Time-series Expression Miner (STEM) to identify clusters of genes, the expression of which changes in a similar manner over time (Figure 2) (Ernst and Bar-Joseph, 2006). We identified three clusters with genes in clusters 0 and 2 significantly upregulated throughout infection, while expression of genes in cluster 1 was slightly decreased 4 DPI before returning to baseline levels of expression at 6 DPI (cluster 1) (Figure 2A). Genes in cluster 0 (n = 2,303) enriched to GO terms related to the innate inflammatory response such as “myeloid leukocyte activation” (e.g., CD14, CD53, and CD55), “cytokine-mediated signaling pathway” (e.g., IFNAR2, IL15RA, JUNB, and LTB), “NIK/NF-kappaB signaling” (e.g., NFKB2, RELA, TRAF2, and TLR9), and “apoptotic signaling pathway” (e.g., CASP1, CASP4, and FAS) (Figure 2B). Genes also enriched to GO term “lymphocyte activation” and were related to both cellular (e.g., CD274, JAK3, TNFRSF14, and TREML2) and humoral (e.g., BCL6, LYN, and PRKCB) immunity (Figures 2B,C). Cluster 2 (n = 495) was comprised of genes that enriched to GO terms indicative of heightened inflammatory state such as “positive regulation of tumor necrosis factor production” (e.g., JAK2 and TNFRSF1A), “myeloid leukocyte activation” (e.g., FCER1G, NLRP3, and TLR4) (Figure 2D), and “defense response to virus” (e.g., BST2, AIM2, and HERC5) (Figure 2E). Lastly, the downregulated genes of cluster 1 (n = 495) were involved in metabolic and cellular processes such as “regulation of cell cycle process” (e.g., ANAPC2, CCND1, and CDC20) and “RNA catabolic process” (e.g., AGO1, LARP1, and PYM1) (Figures 2F,G).

Next, we examined transcriptional changes following challenge with early (Guinea C07) and late (Liberia and Mali) EBOV-Makona isolates (Figure 3). Disease progression and clinical data for these animals were reported in our previous study and indicate uniform lethal infection, albeit slightly delayed following infection with late isolates in rhesus macaques (Marzi et al., 2018). Comparative bivariate transcriptional analysis revealed striking similarities in the overall magnitude and character of the transcriptional response to early and late isolates (Figures 3A–C). The majority of DEGs detected 4 DPI and 6 DPI following challenge with each isolate overlapped, with a greater number of DEGs detected 6 DPI compared to 4 DPI, and most DEGs were upregulated (Figures 3B–D). Liberia exhibited a particularly lower number of DEGs (∼500) at 4 DPI compared to other isolates (>1,000 DEGs) (Figure 3B). These downregulated DEGs enriched to GO terms related to gene expression, including “translation,” “mRNA processing,” and “protein targeting to ER” (Figure 3E). On the other hand, DEGs downregulated by 6 DPI played a role in cell cycle (C07, Mali), adaptive immunity (C07, Liberia, Mali), translation (Liberia), and apoptosis (C07) (Figure 3E). Functional enrichment of upregulated DEGs at 4 and 6 DPI is consistent with a heightened inflammatory response with over-representation of GO terms related to EVD pathology (e.g., “defense response to virus,” “positive regulation of cell death,” and “regulation of cytokine production”) (Figure 3E).

Given that DEGs detected 4 and 6 DPI following infection with all three isolates enriched to similar GO terms, we next examined the overlap between those sets of DEGs (Supplementary Figure 2). DEGs that enriched to GO term “positive regulation of cell death” were upregulated throughout infection and showed significant overlap (Supplementary Figure 2A). These genes played a role in apoptosis (e.g., CASP1, FAS, and TNFRSF1A) and inflammation (e.g., IL6, NLRP3, TLR4, and TNF) (Supplementary Figure 2A). Similarly, many upregulated DEGs enriching to “defense response to virus” were shared between the isolates and encoded nucleic acid sensors (e.g., DDX58, DDX60, and OAS1), components of the IFN signaling pathway (e.g., STAT1, STAT2, and TBK1), and ISGs (e.g., AIM2, IFI16, and ISG15) (Supplementary Figure 2B). Shared DEGS that enriched to “regulation of cytokine production” included a mixture of genes important for inflammation (e.g., IL18R1 and MYD88), chemotaxis (e.g., FLOT1 and ROCK2), and pathogen recognition (e.g., CLEC7A, CLEC4A, and TLR2) (Supplementary Figure 2C). DEGs upregulated at 4–6 DPI in all isolates and enriching to “lymphocyte activation” reflected both B and T cell-mediated immunity (e.g., ARG2, BCL6, LYN, and STAT3) (Supplementary Figure 2D). Similarly, DEGs belonging only to C07 and Liberia for this term also associated with lymphocyte-mediated immunity (e.g., CD81, IL6, and LCP1) as well as T cell regulation (e.g., IDO1, PDCD1LG2, and TCIRG1) (Supplementary Figure 2D).

Complete blood count analysis showed hallmarks of EVD in all four cohorts of animals including lymphopenia, neutrophilia, and thrombocytopenia (Supplementary Figure 3 and Supplementary Table 1). Furthermore, significant declines in platelets were noted as early as 2 DPI but only in Mayinga-infected animals (Supplementary Figure 3 and Supplementary Table 1). Declines in lymphocyte numbers were noted only in Mayinga and C07 infections at 4 DPI (Supplementary Figure 3 and Supplementary Table 1). However, lymphopenia was most pronounced in Mayinga-infected animals (Figure 3B), whereas neutrophilia was not evident in animals infected with the late EBOV-Makona variant Mali (Figure 3C). To gain a deeper understanding of the changes in immune cell frequencies following infection with EBOV-Makona isolates, we performed DCQ. Frequencies of activated NK and DC subsets were predicted to increase over the course of infection for all three isolates (Figures 3F–H). Higher frequencies of monocytes were only predicted following C07 infection (Figure 3F). Significant lymphopenia was predicted to occur following infection with C07 wherein levels of CD4 T cells and several B cell subsets decreased with infection (Figure 3F). In contrast, significant but modest changes in naïve B cells or CD8 T cells were predicted to occur following infection with Mali and Liberia variants, respectively (Figures 3G–H).

We next used STEM to identify groups of genes with similar patterns of longitudinal gene expression (Supplementary Figure 4). In all isolates, we detected a cluster (cluster 0) of genes whose expression progressively increased over the course of infection (Supplementary Figures 4A–C). Genes in cluster 0 from all three isolates enriched to similar GO terms involved in both innate (e.g., “myeloid leukocyte activation” and “regulation of innate immune response”) and adaptive (e.g., “lymphocyte activation” and “T cell differentiation”) immunity (Supplementary Figure 4D). A second cluster (cluster 1) featured genes whose expression levels were slightly downregulated at 4 DPI before returning to near baseline by 6 DPI and enriched to “translation” and “DNA repair” (Supplementary Figure 4E). However, a unique cluster (cluster 2) was identified following infection with the EBOV-Makona Guinea C07 isolate, which consisted of genes robustly upregulated at 4 and 6 DPI. Genes in this cluster played roles in innate immunity (e.g., “myeloid leukocyte activation”), adaptive immunity (e.g., “T cell receptor signaling pathway”), antiviral defense (e.g., “cellular response to type I IFN”), and EVD pathology (e.g., “blood coagulation”) (Supplementary Figure 4F).

We next compared transcriptional changes among animals infected with EBOV-Mayinga or EBOV-Makona isolates (Figure 4). PCA showed a significant overlap between all four infections at each DPI (Figure 4A). We identified a shared core of 704 genes that were involved in antiviral immunity (e.g., IFIH1, IRF7, OAS1, ISG15, and STAT1 – GO terms “defense response to virus,” “cellular response to type I IFN”), inflammation (e.g., NOD2, TICAM1, and TLR4 – GO terms “activation of innate immune response,” “myeloid cell differentiation”), cytokine signaling (e.g., CCL4, INPP5D, and IRAK2 – GO term “cytokine-mediated signaling pathway”), and immune cell activation (e.g., LYN, MAFB, PRKCB, and VAV1 – GO terms “lymphocyte activation”) (Figures 4B,C).

Additionally, there were DEGs unique to each infection. The largest group of unique DEGs was detected following EBOV-Mayinga infection. Those unique DEGs were associated with B cell-mediated immunity (e.g., ATM, CD27, PRKDC, and TCF3), stress responses (e.g., BAX, EYA3, and NUP155), and DNA repair (e.g., COPS2, LIG3, and PYHIN1) (Figure 4D). DEGs detected only following infection with the early EBOV-Makona Guinea C07 isolate played a role in nucleic acid metabolism (e.g., GO term “viral gene expression” and “translation”) and catabolic stress responses (e.g., GO term “positive regulation of apoptotic process”) (Figure 4E). Many DEGs unique to late EBOV-Makona Mali infection enriched to GO terms related to EVD associated processes, notably “blood coagulation,” “regulation of inflammatory response,” and “wound healing” (Figure 4F). The few DEGs unique to late EBOV-Makona Liberia infection played a role in myeloid cell differentiation (e.g., HCLS1 and RARA) and host metabolic processes (e.g., CCS and PN3) (Figure 4G).

To identify clusters of genes with similar patterns of gene expression across all EBOV-Makona and EBOV-Mayinga isolates, we applied a two-way forward regression model using MaSigPro (Figure 5) (Conesa et al., 2006). We retained only those genes that were considered significant in at least 16 comparisons and then clustered them based on temporal expression patterns, which results in four significant clusters (Figures 5A–D). Genes in clusters 1 and 2 exhibited modest upregulation 4 DPI followed by a sharp increase 6 DPI. These 495 genes played a role in the activation and regulation of innate (e.g., “regulation of innate immune response”) and adaptive (e.g., “lymphocyte activation”) immune responses (Figures 5A,E). Notable genes within this cluster are involved in leukocyte–leukocyte interactions (e.g., ICAM1 and ITGB2), myeloid cell signaling (e.g., CD14, RELA, NFKB1, and TRAF3), as well as lymphocyte activation (e.g., BATF, LAT2, and PRKCB) (Supplementary Figure 5A). The 301 genes in cluster 2 mainly played a role in translation (“mRNA processing”), antigen presentation (“autophagy” and “antigen processing and presentation”), and coagulopathy (“blood coagulation”; e.g., A2M, P2RX1, PRKAR2B, and SELP) (Figures 5B,E and Supplementary Figure 5B). Expression of genes in cluster 3 progressively increased from 0 to 6 DPI and enriched to similar GO processes described for clusters 1 and 2 in addition to “type I IFN production”/“response to virus” (e.g., DDX58, IFIH1, IRF1, STAT1/3, and ISG15) (Supplementary Figure 5C). Cluster 4 consisted of 20 genes associated with metabolism that were mostly downregulated over the course of infection (Figure 5D).

To identify gene signatures that can differentiate among the groups, we next applied a sparse partial least-squares discrimination analysis (sPLS-DA) (Figure 6). Contrary to PCA, sPLS-DA permits identification of the minimum number of genes responsible for driving a given component of variation while preserving maximum covariance among defined groups. We identified three components that contributed to the majority of variation among EBOV strains (Figures 6A–C). The six genes that explained component 1 (16% variation) included immunoglobulin light chain IGKV2D-40, intracellular trafficking protein RAB30, and transcription factor TCF4 that regulates lymphoid and plasmacytoid DC development (Figures 6A–D). These six genes were sufficient to separate EBOV-Mayinga and EBOV-Makona isolates, but this separation was heavily influenced by Mali infection (Figures 6D,E). The six genes that explained component 2 (14% variation) separated Liberia from EBOV-Mayinga infections (Figures 6A–D). These six genes played a role in immune signaling (e.g., FCN1 and TYROBP) and ubiquitin-protein ligase complex component (RMND5A) and were upregulated with EBOV-Mayinga and Liberia infection, respectively (Figure 6E). Nearly 300 genes belonging to component 3 (6% variation) distinguished C07 from all other EBOV infections (Figures 6B,C). These genes played a role in host defense as well as cellular homeostasis, as indicated by enrichment to GO terms “myeloid leukocyte activation” (e.g., TLR3, TLR5, CR2, and CCR3), “lymphocyte activation” (e.g., IGHA2, IL7R, ITGA4, and LILRB2), “protein folding” (e.g., UBC and UBE2D1), and “positive regulation of cellular protein localization” (e.g., MAPK1 and TIMP2) (Figures 6F,G). This included a large number of chaperones and ubiquitin-conjugating enzymes (e.g., DNAJA1, HSPA1B, and UBE2D1) and genes modulating lymphocyte-mediated immunity (e.g., IGHA2, IL7R, ITGA4, and LILRB2) that were most prominently expressed in C07 infection (Figure 6G).