Research ArticleEMERGING INFECTIONS

A conserved transcriptional response to intranasal Ebola virus exposure in nonhuman primates prior to onset of fever

See allHide authors and affiliations

Science Translational Medicine  28 Mar 2018:
Vol. 10, Issue 434, eaaq1016
DOI: 10.1126/scitranslmed.aaq1016

An enlightening Ebola model

Using large doses of highly pathogenic agents in animal models can ensure consistent and fully penetrant infection but does not recapitulate human exposure and possibly immune responses. Speranza and colleagues intranasally infected cynomolgus macaques with a small dose of Ebola virus, which led to various disease presentation and allowed for detection of early immunity. Regardless of the timing of symptoms, analysis of the host response revealed a conserved interferon signature that preceded fever by several days. These results shed light on the host response to Ebola virus and established a model that could be used for future studies of pathogenesis, treatment, or prevention.

Abstract

Ebola virus disease (EVD), caused by Ebola virus (EBOV), is a severe illness characterized by case fatality rates of up to 90%. The sporadic nature of outbreaks in resource-limited areas has hindered the ability to characterize the pathogenesis of EVD at all stages of infection but particularly early host responses. Pathogenesis is often studied in nonhuman primate (NHP) models of disease that replicate major aspects of human EVD. Typically, NHP models use a large infectious dose, are carried out through intramuscular or aerosol exposure, and have a fairly uniform disease course. By contrast, we report our analysis of the host response to EBOV after intranasal exposure. Twelve cynomolgus macaques were infected with 100 plaque-forming units of EBOV/Makona through intranasal exposure and presented with varying times to onset of EVD. We used RNA sequencing and a newly developed NanoString CodeSet to monitor the host response via changes in RNA transcripts over time. When individual animal gene expression data were phased based on the onset of sustained fever, the first clinical sign of severe disease, mathematical models indicated that interferon-stimulated genes appeared as early as 4 days before fever onset. This demonstrates that lethal EVD has a uniform and predictable response to infection regardless of time to onset. Furthermore, expression of a subset of genes could predict disease development before other host-based indications of infection such as fever.

INTRODUCTION

Ebola virus disease (EVD) is a serious illness that can cause up to 90% case fatalities in humans and nonhuman primates (NHPs) (1). The causative agent, Ebola virus (EBOV), spills over into human and NHP populations on a sporadic basis and typically causes limited outbreaks. A notable exception was the 2013–2016 outbreak in West Africa, which caused about 28,000 cases (2). This event highlighted the need to better understand the pathogenesis and host response to EVD as a means to create better therapeutics, vaccines, and diagnostics.

The study of EVD is difficult because human clinical samples available for research purposes are rare, and priority is not usually placed on preserving samples for subsequent experimentation during an outbreak. To fill this gap, animal models have been developed to study EVD in vivo. The current model most widely used for understanding pathogenesis and testing of therapeutics and vaccines is NHPs, specifically cynomolgus and rhesus macaques (3). Upon infection with EBOV, NHPs display many of the same symptoms as humans including coagulopathy, lymphocytosis, and rash (4). The NHP EVD model typically uses an infectious dose of 1000 plaque-forming units (PFU) via an intramuscular (IM) injection or aerosol exposure resulting in disease symptoms within 3 to 5 days and 100% lethality within 6 to 10 days after exposure (5, 6). Intranasal exposures at 250 or 500 PFU also result in a uniformly lethal model of EVD but with a delay in time to death compared to 100 PFU IM exposures (7). By contrast, human EVD is thought to result from exposure to lower infectious doses via parenteral or mucosal routes (8, 9) with the incubation period ranging widely from 2 to 21 days after exposure to EBOV (10). These differences highlight the need for experimental NHP models that more closely recapitulate natural human exposure and ensuing disease.

Important data regarding the host response to infection have been gained from investigating circulating immune responses in 100% lethal NHP models (11). For example, transcriptomics has been used to explore early host responses to infection in time-resolved NHP data sets (1215) and to define important aspects of the host response that contribute to survival in mice (16), NHPs (17, 18), and humans (1921). These analyses have pointed to a robust and early immune response after EBOV exposure largely characterized by the up-regulation of interferon-stimulated genes (ISGs).

Left unaddressed in these analyses has been whether the host response to infection in a uniform disease onset model is the same as that seen in a more variable disease onset that is representative of human cases of EVD. We were able to approach this question through a recently developed model of EVD. This model used a mucosal exposure route that was less than 100% fatal and showed variable onset of disease. We investigated the host response to infection in this model using telemetry devices, blood chemistry analysis, and transcriptomic approaches of RNA sequencing (RNA-seq) and NanoString. This allowed us to map the NHP response to infection at both symptomatic and asymptomatic phases and to relate these findings to human EVD.

RESULTS

In the context of ongoing efforts to develop NHP models that more accurately replicate EVD observed in humans, 12 cynomolgus macaques were exposed to EBOV/Makona (Mak) via the intranasal route. Virus (target dose of 100 PFU) was administered either via a pipette or mucosal atomization device (MAD) to each NHP. After exposure, disease onset was not uniform, with animals showing variable onset of disease symptoms and time to death. Two animals were asymptomatic for disease during the course of the experiment.

Description of four distinct response groups after exposure

The NHPs clustered into four distinct groups based on the appearance of disease symptoms, onset of viremia, and time to death (Fig. 1, A and B) and had an 83% fatality rate. Three animals had an expected disease course [group 1 (Normal): NHPs 1, 2, and 12], became quantifiably EBOV-positive as detected by reverse transcriptase (RT)–quantitative polymerase chain reaction (qPCR) on day 6 after exposure, and had a mean time to death of 10.47 days. Four animals had a delayed onset of disease [group 2 (Delayed): NHPs 3, 5, 8, and 10] and had a mean time to death of 13.31 days. In this group, NHPs did not develop quantifiable viremia until days 10 or 12 after exposure, but viral RNA was qualitatively detected (above the limit of detection but below the limit of quantification for the assay) as early as day 3 in NHP 10 and on day 6 in NHPs 3 and 8. NHP 5 did not have any detectable viral RNA until day 10. Three animals had a late onset of disease [group 3 (Late): NHPs 4, 6, and 9] and a mean time to death of 21.42 days. NHP 6 had quantifiable viremia on day 20, but viral RNA was only qualitatively detected on days 6, 10, and 14. NHP 9 first had quantifiable viremia on day 21 but was also qualitatively positive on day 14. NHP 4 had borderline quantitative viremia on day 10, was qualitatively positive on day 14, and then was quantitatively positive again on day 20. Finally, two animals did not become quantitatively RT-qPCR–positive during the time course of the experiment [group 4 (No Response): NHPs 7 and 11] and survived to 41 days, which was the predefined experimental end point. However, viral RNA was qualitatively positive in these animals on days 6 and 10 after exposure, respectively. Survival and time to death did not correlate with the method of intranasal exposure (P = 0.1292).

Fig. 1 Overview of the separation of animals and analysis performed.

(A) Table showing the different animals (rows) and the different days after exposure (columns). White squares signify that the animal was not polymerase chain reaction (PCR)–positive at that time point (n.d., not detectable), blue squares indicate that they were RT-qPCR–positive, and gray squares mean that the RT-qPCR results were below the limit of quantification. The number of plus signs in the gray boxes represents the number of RT-qPCR replicates that were positive. Hashed boxes indicate no data for this time point. (B) Comparison of the four groups and their survival curves. The x axis is the days after exposure (DPE), and the y axis is the percent survival. The green line is the Normal group [nonhuman primates (NHPs) 1, 2, and 12], the blue line is the Delayed group (NHPs 3, 5, 8, and 10), the red line is the Late group (NHPs 4, 6, and 9), and the black line is the No Response group (NHPs 7 and 11). (C) Overview of the RNA quantification analysis performed on the animals. The x axis is the days after exposure, and the y axis is the different animals organized by outcome groups (Normal, green; Delayed, blue; Late, red).

To determine whether physiological differences that correlated to disease outcome were detectable before exposure, we analyzed blood chemistry data, hematology, and soluble proteins collected before exposure (day 0). A one-way analysis of variance (ANOVA) was performed for each of the 39 measured parameters, and a regression coefficient to symptom onset (R2) value was calculated. Only one variable, mean corpuscular volume (MCV), was significantly different between the groups (P < 0.05; fig. S1A). The difference observed in MCV was driven by the Delayed animals, which had a significantly lower value than the Normal group (P = 0.02; fig. S2) that was not significantly different from any other group. A low MCV value, which measures the average size of red blood cells, can signify blood loss or anemia. However, because the values were only significant in one of the group comparisons, we think that it is unlikely to have contributed to the differential disease onset.

To further investigate whether there was a difference observed at the day of infection, we performed principal components analysis (PCA) on the same parameters from day 0 (fig. S1B). We would hypothesize that if the animals were similar across the different response groups before exposure, then one group would not cluster differently from any of the others. There was no apparent clustering within groups, signifying that no group was distinctive from the others. Also, all the groups are centered close to the point (0,0), meaning that none of the groups had large outlier values that could have contributed to the different onset of disease. On the basis of this analysis, we can conclude that it was unlikely that a preexisting physiological difference in the animals contributed to the different disease onsets.

The variable onset of disease in these animals presented a unique opportunity to identify factors that might be associated with different responses to EBOV intranasal challenge. We investigated clinical, biochemical, physiological, and transcriptomic data to identify markers of disease onset within this cohort. An overview of the transcriptomic data available can be seen in Fig. 1C. We used two approaches for transcriptomic analysis—RNA-seq and NanoString. RNA-seq enables a hypothesis-neutral investigation of RNA abundance changes. By contrast, we also used a newly developed NanoString CodeSet that recognizes 769 NHP transcripts, which we have previously demonstrated as effective with studying EBOV/Mak-infected NHPs (15). Although NanoString does not provide the transcriptomic depth of RNA-seq, it is much more rapid compared to RNA-seq and does not have as high of a requirement for RNA quality.

Clinical chemistry changes in symptomatic animals

Typical manifestations of EVD were observed in all nonsurvivors, including maculopapular rash, lymphadenopathy, anorexia, and fever. Indicators of coagulopathy, such as bleeding or hemorrhage, were largely absent in these animals. Some animals displayed respiratory signs including cough, labored breathing, and use of accessory or abdominal muscles for breathing. In addition, all nonsurvivors had alterations in hematological and serum chemistry parameters, regardless of the time to death. These changes mirrored previously reported EBOV infection of NHPs, including increases in blood urea nitrogen, creatinine, calcium, and gamma-glutamyl transferase (Fig. 2). Increases in liver enzymes, such as alanine aminotransferase and aspartate aminotransferase, were also noted, along with an increase in alkaline phosphatase (Fig. 2). Changes in blood chemistry coincided with the development of clinical signs and appearance of viral RNA in plasma; these changes were delayed in those animals with an extended time to death. The two survivors generally did not have any changes in behavior, health, or clinical pathology parameters associated with EVD.

Fig. 2 Blood chemistry data for animals over time.

(A) Blood urea nitrogen (BUN), (B) creatinine (CRE), (C) calcium (CA), (D) gamma-glutamyl transferase (GGT), (E) alanine aminotransferase (ALT), (F) aspartate aminotransferase (AST), and (G) alkaline phosphatase (ALP). Normal, green; Delayed, blue; Late, red; and No Response, black.

Innate immune pathway activation during acute EVD

The varying times to death among the study animals led us to investigate the nature of the host immune response after the onset of symptoms as determined by onset of sustained fever. We performed differential gene expression analysis in the RNA-seq data set and biological pathway analysis with Ingenuity Pathway Analysis (IPA) for each group relative to prechallenge baseline for symptomatic time points (Fig. 3A). All animals showed a similar host response at the end stages of disease as evidenced by up-regulation of the same top 20 signaling pathways based on IPA results (Fig. 3B). The highest level of conservation across the groups was the Normal and Delayed groups at 10 days after exposure and the Late group at 21 days after exposure. Because of the loss of samples from poor RNA quality for day 14 after exposure in the Delayed group, we only had usable samples from one animal remaining at this time point, thus reducing our ability to find statistically significant differentially regulated genes resulting in fewer pathways with enrichment.

Fig. 3 Comparison of leading pathways and upstream regulators.

Overview of the most strongly differentially regulated pathways in the viremic animals. (A) Overview of the samples analyzed (large circles) that were symptomatic (red circles) with the number of biological replicates denoted in the circles. (B) Significantly differentially regulated pathways. The darker purple indicates a strong up-regulation of the pathways (positive z score), and the green means a down-regulation of the pathways (negative z score). The P values represent whether there is a significant number of genes in the pathway that are differentially expressed as determined by Ingenuity Pathway Analysis (IPA) through a Fisher’s exact test. *P < 0.05, **P < 10−3, and ***P < 10−6. (C) Upstream regulators from IPA. Purple denotes up-regulation, and green denotes down-regulation (log fold change). (D) Breakdown of the interferon response to the type of interferon responsive genes. The x axis is the specific types (type I, II, or III only) and overlap between the different responses. The green bars are samples from the Normal group, blue bars are samples from the Delayed group, and red bars are samples from the Late group. The different patterns correspond to different days after exposure.

Enriched pathways included activation of innate immune pathways [role of pattern recognition receptors, interferon signaling, role of retinoic acid–inducible gene I (RIG-I) in antiviral response] and acute phase response signaling and pathways associated with nitric oxide and reactive oxygen species. These pathways were consistent within each group and were also consistent with host gene expression seen during viremic stages of disease after higher-dose exposure to EBOV/Kikwit (22) and other EBOV/Mak challenges (13). We identified two pathways among the top 20 that were down-regulated (peroxisome proliferator–activated receptor signaling and LXR/RXR activation) during infection. These pathways are important in lipogenesis and cell survival and can down-regulate the inflammatory response through inhibition of nuclear factor κB and Toll-like receptor 4 (23, 24). The top 20 regulators of downstream targets were also consistent with an overactivation of the innate immune response and show a strongly conserved host response during the symptomatic stage of disease (Fig. 3, B and C).

To further characterize the interferon response during the symptomatic stage of infection, we used the interferome database (25) to compare type I, type II, and type III interferon responses based on the number of differentially regulated interferon responsive genes. Across all the symptomatic animals there is a strong type I and type II interferon (IFN) response but minimal type III response (Fig. 3D). This is consistent with the upstream regulators where interferon-α (IFNA) and interferon-γ (IFNG) are significantly enriched (P < 10−6) across all the groups.

Cytokine gene expression during acute EVD

Using differentially expressed cytokine genes identified through RNA-seq, we examined whether animals that showed symptoms at different times of infection had differences in their cytokine response during the course of disease. We identified 57 differentially expressed cytokine mRNAs (table S1) and then determined whether the animals showed similar expression based on the groups or on the distinction between symptomatic or asymptomatic. We found that the dominant driver of clustering in this approach was symptom onset (Fig. 4A). By contrast, there was little clustering based on time after EBOV exposure further suggesting that the symptom onset is a more important indicator of disease course.

Fig. 4 Late-state cytokine response conservation.

Analysis of the conservation of the cytokine response. (A) Correlation plot using only the cytokine genes. The rows and columns are the same and represent a single group at a single time point. The dendogram is colored based on a k-means clustering with two centers. The first color bar at the top shows the groups (Normal, green; Delayed, blue; Late, red). The second color bar indicates when the groups become symptomatic with orange for symptomatic and with purple for those animals still in the incubation period. (B) Venn diagram of cluster 2 samples with cytokine genes that are differentially expressed. Spots lacking a value indicate 0 genes.

Of the cytokines in the RNA-seq data set that were differentially regulated, 11 were common to all groups and time points (Fig. 4B). Again, because of the sampling size being reduced to one animal in the Delayed group at day 14 after exposure, fewer differentially regulated genes were identified. Comparing only the Normal group at days 6 and 10 after exposure, Delayed group at day 10 after exposure, and the Late group at day 21 after exposure, the number of conserved cytokines increases to 22. This suggests that the cytokine response is similar in the symptomatic animals. To confirm the host response profiles that we observed with next-generation sequencing RNA-seq, we also used a NanoString CodeSet to profile the expression of 769 NHP genes that has been previously shown to be beneficial as a secondary means for RNA quantification in NHPs infected with EBOV (15). This analysis confirmed the up-regulation of the cytokine genes identified by RNA-seq (fig. S3).

Early and transient transcriptional responses in NanoString data

The NanoString NHP gene expression codeset has the potential for increased sensitivity versus RNA-seq for targeted gene expression profiling (15). Notably, in this data set, we identified a subset of cytokines that were up-regulated at day 3 after exposure in NHPs that showed a delayed course of disease (NHPs 3, 5, 8, and 10).

This up-regulation was consistent across all four animals in this cohort and was not identified as up-regulated in NGS data sets. This finding led us to broaden our analysis of the NanoString data to see whether there were other genes up-regulated at this early time point. Subsequently, we identified 178 genes that were found to be up-regulated on day 3 after exposure in the Delayed animals but had returned to near baseline by day 6 after exposure with the exception of three ISGs (Fig. 5, A and B, fig. S4, and table S2). We performed k-means clustering (k = 3) on the NanoString expression data for all NHPs at all time points (Fig. 5C) using these 178 genes. In the PCA plot, the three clusters are shown. Cluster 1 is made up of most samples. Cluster 2 represents the Delayed animals at day 3 after exposure. One animal from the Normal group (NHP 2), one from the Late group (NHP 9), and one from the No Response group (NHP 11) at day 3 after exposure are also included in this cluster, showing a similar up-regulation at this time point. Three points formed Cluster 3. All three were end-stage samples that were close to death.

Fig. 5 Early and modest immune response at day 3 after exposure in Delayed animals.

Analysis of the NanoString samples with only the genes identified as being up-regulated at an early time point (day 3 after exposure, Delayed). (A) MA plot of the 178 up-regulated genes at day 3 after exposure in Delayed animals. Each point represents a gene. The x axis is the log2 of the base mean counts, and the y axis is the log2 fold change relative to day 0. (B) MA plot for the same genes at day 6 after exposure. (C) Principal components analysis (PCA) of the 178 up-regulated genes across all samples. K-means clustering was performed on all principal components with k = 3 to determine the clustering. The PCA plot is showing the first two principal components. The colors represent the different clusters.

We used PANTHER (26) to determine whether any one set of genes was significantly enriched in these 178 genes. As a reference, we used only those genes included in the NanoString CodeSet to account for potential bias. We found that Gene Ontology (GO) terms associated with cytokines, cytokine response (GO:0005125) and cytokine receptor binding (GO:0005126), were significantly enriched (adjusted P = 8.91 × 10−5 and 6.05 × 10−3, respectively). Cytokine response genes that were up-regulated at day 3 after exposure in Delayed animals included interferons (for example, IFNB1, IFNL2, and IFNL3), proinflammatory cytokines (for example, IL17F, IL6, and IL10), and many T cell recruitment/activation genes (for example, CCL18, IL2, and CCL17). Although there was detection of type I and type III interferons, there was no significant up-regulation of ISGs at this point. This may suggest a modest immune response to the infection in the nasal cavity and that day 3 after exposure was optimal for detection with the Delayed animals. Transcriptional analysis of early events in the nasal cavity would be needed to support this hypothesis.

ISG activation immediately before the onset of EVD symptoms

As the response to infection in the symptomatic stage offered few differences to explain disease onset, we moved to studying events after infection at the time point directly preceding the symptomatic stage (Fig. 6A). Previous studies have suggested that transcriptional analysis can identify host innate immune response genes that are induced after hemorrhagic fever virus infection and that some of these genes can be detected before the onset of symptoms (12, 2729). To determine whether there was any evidence of host-response gene expression changes at early times after exposure, we examined the RNA-seq data for evidence of transcripts that were induced before the display of overt symptoms (fever as determined by telemetry) or quantifiable viremia (determined by RT-qPCR) and that were then sustained through the end of infection. Our analysis identified a number of genes that were induced at or before induction of fever. In challenged animals that did not develop disease, we did not find these genes to be significantly differentially regulated. Many of these “sentinel” genes are classified as ISGs.

Fig. 6 Analysis of the pre-viremic ISG response.

(A) Analysis of the onset of the interferon-stimulated gene (ISG) response. Large circles represent samples that were used in the analysis with red denoting symptomatic samples. The numbers indicate how many samples were used. (B) Comparison of the onset of viremia to induction of the group of ISGs. The days after exposure are on the x axis. The left y axis is the PCR values (black line), and the right is the number of genes annotated as ISGs that are up-regulated (blue line). (C) Similar plot for the Delayed animals, (D) Late animals, and (E) No Response animals.

The transcriptome from animals in the Normal group showed increased accumulation of about 100 ISGs concurrent with the onset of quantifiable viremia and symptoms (Fig. 6B). In animals in the Delayed group, 14 ISG transcripts were significantly up-regulated [log2(fold change) > 1, adjusted P < 0.05] by day 6 after exposure. This represented a gene expression spike that occurred 4 days before the onset of quantifiable viremia and fever (Fig. 6C). The presymptomatic expression of host innate immune genes was even more remarkable in animals that showed a late onset of disease. ISG transcripts were up-regulated at 14 days after exposure, which was 7 days before the onset of quantifiable viremia (Fig. 6D).

These data showed that ISG mRNA up-regulation is a presymptomatic indicator of the development of EVD. This is summarized in Table 1, which shows that the ISG markers are highly expressed by day 6 in the Normal group, coincident with fever and quantifiable viremia, but appear markedly in advance of both fever and quantitative viremia in cases of delayed disease onset. In the Delayed and Late groups, the interferon response preceded a sustained fever by 4 and 5 days, respectively. This presymptomatic response was conserved. Of the four ISGs up-regulated at day 14 after exposure in the Late animals (NHPs 4, 6, and 9), three were identical to ISGs up-regulated at day 6 after exposure in the Delayed animals. All were up-regulated in animals that showed the expected timing of symptom development.

Table 1 Comparison of the ISG response, quantifiable viremia, and fever.

DPE, days post exposure; N.D., not detectable.

View this table:

The early up-regulation of these ISGs was also seen using the NanoString CodeSet, which contained 80 ISGs. Differential expression analysis of genes in the NanoString CodeSet showed that all symptomatic animals had up-regulated ISGs, whereas asymptomatic animals (group 4, No Response) did not (fig. S5). As was seen in our NGS data sets, at day 6 after exposure in the Delayed group and day 14 after exposure in the Late group, there was an up-regulation of IFIT2, IFIT3, IFIT5, OASL, ISG15, and DDX58. Expression of these genes persisted through to the late stage of disease and was conserved between animals in each group (fig. S5).

Conservation of ISG response to human EVD

The ISG response during the acute phase of EVD in NHPs was highly conserved among the groups. Figure 7A shows the fold change of ISGs in Normal animals at day 10 after exposure (x axis) compared to the Delayed animals at day 10 after exposure and the Late animals on day 21 after exposure (y axis). The lines represent linear fits. As can be seen, the ISGs are highly similar in their expression values across the different groups, suggesting that mRNAs that make up the response were highly conserved despite different times to the onset of disease. The conservation of the late-stage animals (Delayed and Late groups) and the early response suggests that the ISG response is up-regulated at similar times and in a similar order. There is a significant correlation of both the Delayed group ISG expression to the Early group ISG expression (Pearson’s r = 0.93, P < 0.0001) and the Late group ISG expression to the Early group ISG expression (Pearson’s r = 0.89, P < 0.0001), demonstrating a high conservation of ISG expression.

Fig. 7 Correlation of NHP and human ISGs.

Analysis of the fold induction correlation with regard to the ISG response in NHPs to human fatalities. (A) Correlation plot of all ISGs that are identified as being differentially regulated at any time point in the NHP data set. The x axis is the log2 fold change for the Normal animals at day 6 after exposure, and the y axis is the log2 fold change for the Delayed animals at day 10 after exposure (red points and line) or the Late animals at day 21 after exposure (blue points and line). The lines are linear fits. (B) Correlation of the human fatal cases (log2 fold change compared to controls) on the x axis and the fold change of the Late group (log2 fold change compared to controls) at day 21 after exposure. The line represents a linear fit of the data. (C) Correlation plot similar to (B) but with only the early ISGs shown. (D) Table of the genes in (C) with their comparison of the log2 fold change in humans and NHPs. Reported P values are from a correlation test.

To determine whether the ISGs up-regulated in NHPs were also up-regulated in clinical patients with EBOV/Mak, we compared the fold induction of the ISGs that were identified in this study to ISGs that showed increased expression in a transcriptional analysis of EBOV-infected patients in Guinea (19) who would go on to succumb to disease. These patients were all symptomatic for disease at the time of sampling and had a mean time from self-reported symptom onset of 5.9 days (19). A correlation plot shows that many of the same ISGs were up-regulated in both NHPs and in human host responses to EBOV infection (Fig. 7B). From this analysis, we identified a small subset of ISGs that could act as early indicators of disease in NHPs and humans: IFI44, IFI44L, IFIT2, IFIT3, IFIT5, ISG15, MX1, and OASL (Fig. 7, C and D). On the basis of their profile of expression in the NHPs and their conservation into the human data set, we suggest that, at a presymptomatic time point, these genes could act as early indicators of viral infection.

Modeling host dynamics of infection

This data set provided a chance to look at modeling the host response relative to the end of the incubation period to determine (i) whether there is a way to align the samples to generally analyze the host response across the three symptomatic groups and (ii) whether we can model the different gene groups to more robustly determine which will be up-regulated first and how we expect these genes to be up-regulated. The gene expression changes across the four groups did not align well when comparing the fold changes observed in the NanoString CodeSet to the days post exposure (Fig. 8A). However, using the telemetry data, we were able to determine a mean time to onset of fever for each of the days after exposure for each group. We assigned 0 as the mean day of onset of sustained fever (>1.5°C above mean temperature). All the sampling days were then either negative (meaning that they happened before the onset of a fever) or positive (meaning that they happened after the onset of a fever). When we looked at the gene expression data using time relative to onset of sustained fever, we found that the samples aligned very well to each other (Fig. 8B). We concluded that, despite differential onset of symptoms and gene expression, the pattern of gene expression relative to fever onset is conserved.

Fig. 8 Host response to infection modeled relative to fever onset.

Alignment and modeling of gene expression across the four NHP groups. (A) Example gene (ISG15) expression (log2 fold change, y axis) in the NanoString data set relative to day after exposure (x axis). (B) Expression of the same example gene (log2 fold change, y axis) relative to the onset of a fever (0) on the y axis. The green line is the Normal animals, the blue line is the Delayed animals, and the red line is the Late animals. (C) Logistic model fit for the expression of ISG15 relative to the onset of fever. The points are the mean log2 fold change for a group at its time relative to fever onset, and the solid line is the model fit for the fold change. (D) Comparison of the early ISG expression (black lines) logistic model fit relative to the onset of fever (x axis) to many cytokine genes (red lines). The dashed line is for a reference of when the gene expression crosses the threshold of a log2 fold change > 2.

Next, we used a logistic model to fit the data to curves. This type of model does limit our analysis to genes that show a pattern of expression consistent with a logistic curve. We fit a single model to a single gene at a time (Fig. 8C). We did not include the time point of day 3 after exposure because we wanted to examine gene expression relative to acute disease, and the day 3 post-exposure spike of cytokines is eliminated by day 6 after exposure. We found that the ISGs identified as early up-regulators preceded expression of cytokine genes in all cases (Fig. 8D). Also, their predicted mean time to detectable expression (log2 fold change > 2) was found to be 4.6 days before the onset of fever (±0.46 days). This means that the early interferon genes are strongly up-regulated 4.5 days before the onset of symptoms. Finally, this response is followed by the cytokine response, which reached higher expression overall than the interferon response (Fig. 8D) but was not significantly expressed earlier than the end of the incubation period. This shows that the cytokine response is concurrent with the onset of fever.

DISCUSSION

A marked aspect of our results is that, despite the different times to onset of disease, the host response was similar across the three symptomatic groups [groups 1 (Normal), 2 (Delayed), and 3 (Late)]. The Delayed and Late groups each had a longer mean time to death than the standard 6 to 10 days. By shifting the sampling days to align the different response groups by time to symptom onset instead of days after exposure, we observed that the host response follows a conserved and predictable pattern of gene expression. We were able to then generate mathematical models of gene expression to show that presymptomatic gene expression can be detected up to 4 days before the onset of fever.

An important implication of this study is the potential use of host markers as presymptomatic indicators of disease. We found that, through modeling the gene data relative to the onset of a fever, the interferon response was detectable up to 4 days earlier. When we corrected for sampling days, ISGs expression was still detectable up to 3 days before the onset of fever. Earlier studies have shown that gene expression changes are detectable in vivo after EBOV exposure (12, 30, 31), but our results here greatly expand the utility of these markers. The appearance of host-based markers up to 4 days before fever suggests that they could be used to track individuals at high risk for infection including hospital personnel caring for EBOV-infected individuals, close contacts of the same, or individuals who have been in an environment where EBOV is present. On the basis of our work here, NanoString-based assays appear to be a feasible approach as a host-based diagnostic, but other multiplexed nucleic acid detection platforms could be equally useful.

The interferon response in symptomatic animals was followed closely by a cytokine response that was more concurrent with the onset of a fever. This cytokine response was similar to what has been seen previously in NHP infections with 100% lethality including cytokines such as IL6, CCL8, and CXCL10 (12). This cytokine response developed after the interferon response and generally reached higher maximum levels of expression than ISGs. Similar to our modeling with the expression patterns of interferon-responsive genes, we found that we could also model the timing of cytokine expression. The cytokines found to be up-regulated are consistent with the cytokine storm that is characteristic of dysregulation of the immune system (32). This is consistent with a strong proinflammatory response in innate immune cells circulating in the blood.

The response observed during the symptomatic phase is consistent with transcriptional observation in human infections. In a patient treated at the National Institutes of Health Clinical Center, the earliest response to viral infection involved many ISGs and cytokine expression (20). Because this patient ultimately survived, these responses eventually subsided. In a Guinea cohort of patients from the 2013–2016 EBOV outbreak in West Africa, one of the gene sets strongly enriched in the fatal cases was ISGs (19). When we compared the expression of the early genes to human fatal cases, a correlation between the genes during the acute phase in both the humans and the NHPs was apparent. This correlation and the conservation of interferon genes at the earliest stages of disease suggest that these genes may be up-regulated early in infection in humans. This introduces the possibility of using the host response as an early diagnostic, and preliminary work to do so has been performed (30). Although the interferon response is not unique to EBOV, identifying other genes with a similar expression pattern may help to isolate potential early markers of infection.

Smaller changes in RNA abundance were detected in our NanoString assays, allowing for the characterization of a modest immune response very early in infection. This response mostly contained proinflammatory cytokines and general immune response genes. This response could be the animals responding to the infection in the nasal cavity, which was the site of initial exposure. The presence of interferon-λ (IFNL) may suggest an immune response in epithelial cells consistent with a modest immune response occurring in the nasal cavity. IFN-λ has been previously shown to be an important factor in protection of respiratory epithelial cells against viral infections (33).

The early immune response detected by NanoString was present in all the Delayed animals, one Late animal, one Normal animal, and one No Response animal. The No Response animal that had a detectable early response at day 3 after exposure did have more instances of elevated body temperature than the other No Response animal as determined by implanted telemetry devices. The combination of increased instances of fever over the course of the study and the presence of the early response suggests that the two No Response animals had different host responses to the initial infection. This early host response was no longer detectable by day 6 after exposure in all the animals except the one Normal animal because that animal had viremia and fever at this time. During this initial immune response, it is likely that myeloid cells such as dendritic cells and macrophages became infected as they are initial targets for EBOV (34) and found to be major contributors to the host response (35). This would contribute to viral dissemination and later transcriptional changes.

One potential limitation of this study is the small sampling size that is a by-product arising from the challenges associated with performing NHP studies in high containment. Although large changes could be detected and described, a more complete analysis of the changes would require additional studies with a similar outcome of disease. Furthermore, we could not make any claims about the route of infection and its relationship to the differential onset of disease. Although there was a separation of animals based on incubation time, we cannot rule out that the Normal or Delayed groups caused a secondary infection in the Late animals that eventually caused them to spike a fever. Nevertheless, it is important to emphasize that the host response to infection after all infection events (including possible secondary events) was highly conserved in the experiment following a predictable pattern. We were unable to find an association between exposure route and survival or exposure route and time to onset of disease. This may become clearer as additional NHP studies with similar outcomes are performed.

Overall, this animal model provided a unique opportunity to study the host response to infection with differential onsets of disease and time to death. Also, the combination of telemetry, blood chemistry data, and gene expression allowed for an in-depth analysis of the exact timing of the host response relative to clinical presentation. This study highlighted that once the disease has started, it follows a predictable pattern over time with interferon genes being up-regulated early, followed by fever and cytokines, and ultimately the end stage of disease, which has excessive immune dysregulation. Despite the different onset of symptoms between animals, the host response, once started, occurred in about the same time frame and can be modeled to create predictions of disease progression in the host.

MATERIALS AND METHODS

Study design

Twelve healthy adult cynomolgus macaques (Macaca fascicularis), equally distributed by sex, were used for this study. All macaques were surgically implanted with T27F-1B radiotelemetry devices (Konigsberg Instruments) at least 12 months before the start of the experiment and were allowed to fully recover. Surgeries were done under anesthesia and, if necessary, with post-operational buprenorphine analgesia to minimize pain.

Animals were randomly assigned to receive a target dose of 100 PFU EBOV/Mak by intranasal exposure either by pipette (n = 6, animals 1 to 6) or by MAD (n = 6, animals 7 to 12). Predefined end-point euthanasia criteria were used to determine end points with the total experiment planned to run for 41 days.

Research was conducted under an Institutional Animal Care and Use Committee–approved animal protocol at the United States Army Medical Research Institute of Infectious Diseases (USAMRIID). This protocol complied with the Animal Welfare Act, the U.S. Public Health Service Policy, and other Federal statutes and regulations relating to animals and experiments involving animals. The facility where this research was conducted is accredited by the Association for Assessment and Accreditation of Laboratory Animal Care, International and adheres to principles stated in the Guide for the Care and Use of Laboratory Animals, National Research Council (NRC), 2011. All experiments were conducted in an animal biosafety level 4 (ABSL-4) laboratory at USAMRIID. Technicians were not blinded to the method of intranasal delivery.

Virus

The virus stock used has been previously described in (15). Briefly, a master stock was made from an isolate from the serum of a 2014 Sierra Leone fatal human case (SL 3864.1; GenBank accession #KR013754.1). The sequence obtained from the clinical specimen SL-3864 was identical to 17 other P0 sequences from the outbreak. Eight of those were from nonsurvivors. Thus, the sequence of SL-3864 is identical to the sequence of lethal viruses. The stock did not have detectable amounts of contaminants, such as mycoplasma, endotoxin, or adventitious agents, based on the assays and techniques used. Particles were examined through electron microscopy and determined to have characteristic shapes of filoviruses with a uniform diameter of about 80 nm.

Challenge

Animals were anesthetized with an IM injection of ketamine before challenge. The pipette group received 250 μl of virus inoculum in each nare dropwise by a pipette. For the MAD exposure, the device (Teleflex Medical) was attached to a 1-ml syringe containing the virus inoculum. The animals’ heads were tilted so that the nostrils were angled upward. Each animal received 0.5 ml of virus inoculum in each nare delivered as a single puff. The head remained tilted upwards for 1 to 2 min after challenge to ensure proper delivery of the challenge material. The challenge dose for all animals was determined by neutral red plaque assay as previously described (36) was 65 PFU.

In-life study activities

Animals were monitored at least once daily for general health and well being. Additional checks were performed when animals began to display signs of illness. Physical examinations and blood collection were performed under anesthesia on days 0 (before challenge), 3, 6, 10, 14, 21, 28, 35, and 41 (end of study) and at terminal. Euthanasia was performed via intracardiac administration of a pentobarbital-based euthanasia solution under deep anesthesia when an animal was deemed moribund (37).

Clinical pathology

Chemistry parameters were evaluated in serum samples using the Piccolo Point-Of-Care Analyzer (Abaxis) and a Piccolo General Chemistry 13 panel. Hematology was assessed on K3 EDTA-treated whole blood using the ADVIA 120 Hematology System (Siemens).

Telemetry monitoring and pathophysiological analysis

The implanted radiotelemetry devices allowed for the continuous monitoring of body temperature, aortic pressure, left ventricular pressure, heart rate, and respiratory rate throughout the baseline and study periods. The telemetry data were captured and analyzed using a NOTOCORD-hem Evolution software platform (version 4.3, Notocord Inc.). Thirty-minute averages of the above listed physiologic parameters were calculated for each subject. Telemetry data obtained during the 5 days before challenge were used to calculate baseline values and provided the average and SD of each 30-min time period of a 24-hour day. Using the average and SD determinations from the baseline physiologic data, values greater or less than 3 SD of the corresponding baseline value were considered significant. Fever was defined as >1.5°C above its corresponding baseline value, and hyperpyrexia was defined as >3.0°C above its corresponding baseline value. The phases of infection were defined as, (i) the incubation period (time from challenge to onset of sustained fever response for more than 2 hours), (ii) the clinical period (time from onset of sustained fever until onset of sustained hypotension, decline in blood pressure for more than 2 hours), and (iii) the terminal period (time from onset of sustained hypotension until death or euthanasia).

Viral RNA isolation and RT-qPCR

Freshly collected plasma samples (100 μl) were combined with 300 μl of TRIzol LS (Thermo Fisher Scientific), and then RNA was extracted manually using either the QIAamp Viral RNA Mini Kit (Qiagen) or with an EZ1 Advanced XL system using the EZ1 Virus Mini Kit (Qiagen). Viral genome equivalents per milliliter (ge/ml) were quantified using an EBOV-specific RT-qPCR assay and a synthetic RNA standard curve, as previously described (38). Alternatively, PFU equivalents per milliliter (PFUe/ml) were extrapolated from the data using an RNA standard curve derived from a viral stock of known infectivity.

RNA isolation from whole blood

Freshly collected whole-blood samples (500 μl) were diluted with an equal volume of molecular biology–grade water (500 μl) and then combined with 3 ml of TRIzol LS. For RNA-seq analyses, RNA was isolated from whole blood in TRIzol LS using the PureLink RNA Mini Kit (Thermo Fisher Scientific), and RNA was evaluated for quality on the Agilent 2200 TapeStation. For NanoString analyses, total RNA was isolated using the miRNeasy Mini Kit (Qiagen) and then quantified using the Qubit RNA HS Assay Kit (Thermo Fisher Scientific).

RNA sequencing

Libraries were generated on the Sciclone G3 Liquid Handling Robot (PerkinElmer) using the TruSeq Stranded Total RNA Library Prep Kit (Illumina). Library quality was evaluated on the Agilent 2200 TapeStation 2200 and quantified by qPCR using the KAPA Complete (Universal) qPCR Kit (Kapa Biosystems) for Illumina libraries. Libraries were diluted to 12 pM, and cluster generation was performed on the Illumina cBot. Libraries were sequenced on the HiSeq 2500 using the paired end 2 × 100–base pair, dual-index format.

NanoString RNA collection and run

Total RNA samples were analyzed using a custom-designed NanoString nCounter gene expression codeset targeting 769 NHP genes. For each gene, a single probeset was designed with predicted cross-hybridization to homologous genes for cynomolgus macaque and rhesus macaque (Macaca mulatta). Probeset-target RNA hybridization reactions were performed according to the manufacturer’s protocol. For each hybridization reaction, 100 ng of total RNA was used or any quantity that was present in a 5-μl aliquot of purified RNA if less than 100 ng. Purified probeset-target RNA complexes from each reaction were processed and immobilized on nCounter Cartridges using the nCounter MAX Prep Station, and transcripts were quantified on the Digital Analyzer (GEN 2).

RNA-seq data processing

Raw sequencing reads were trimmed on the edges for low-quality reads and filtered using the FASTX-Toolkit version 0.0.14. Paired reads were then aligned to the cynomolgus macaque genome version 5.0 (39) using Bowtie2 version 2.2.6 (40) and Tophat version 2.1.0 (41) with default alignment parameters. Finally, count tables were generated using HTSeq Count version 0.5.4 (42) and Python version 2.7.3. Counts were generated using an alignment to a gene and parameter --stranded set to no.

Count tables were then read into R version 3.2.2. Data normalization was performed with the rlog function of DESeq2 (43). Genes were filtered to only include genes whose median count was greater than 1. RNA-seq samples were filtered if the number of aligned reads was less than 1 million mapped reads to genes or if they were determined to be outliers by PCA. In total, six samples were classified as outliers and were not included in the analysis. Finally, differential gene expression analysis was done through DESeq2, and statistical significance was determined if the absolute log fold change was greater than 1 and the Benjamini-Hochberg (BH)–corrected P value was less than 0.05.

Human data sets were collected from BioProject PRJNA352396 and processed according to (19). The data was read into R using DESeq2 to generate a fold change table. Gene conversion between human and macaque Ensembl IDs was done with Ensembl BioMart (44).

NanoString data processing

NanoString raw counts were compiled into a single count table. Normalization of the counts was done using a standard method. First, a global background was calculated using the mean plus 2 SDs of the negative controls. Then, normalization to internal positive controls was performed to adjust for system variation between lanes. This was achieved by first calculating the positive control factor that was the mean of the geometric means of the positive controls of the different lanes. A normalization factor was calculated for each lane by dividing the positive control factor by the geometric mean of each lane. This value was used as a multiplier for each lane’s count values.

To determine genes that would perform as good input normalization factors, we first filtered all the genes to include only those whose count values were above the global background in all samples. We then applied the NormFinder R package (45) on these genes to determine which were most stable in the data set. The top five most stable genes were used for reference gene normalization across the samples to control for input variation of the amount of RNA. This was done similarly to the positive control normalization.

Differential expression analysis of the NanoString data was performed using a one-way ANOVA test on each gene for the different conditions. This was achieved through the lm function in R. P values were adjusted with a BH correction.

Statistical analysis and gene set analysis

PCA was performed using the stats package prcomp function in R. Clustering was determined by k-means clustering with the cluster package (46) pam function in R. Significant blood chemistry data was determined with a one-way ANOVA test. Significance was determined by a P < 0.05. All P value adjustments were done with a BH correction using the p.adjust function in the stats package. Gene set analysis was done using a Fisher’s exact test on a predefined gene set using standard cutoffs for the differential expression analysis (adjusted P < 0.05 and absolute log2 fold change > 1). Pathway analysis was done using the IPA knowledge base. Interferon-type determination was done using the Interferome knowledge base (25). Statistical significance was determined with a hypergeometric test to determine whether the overlap of interferon response type was significant.

Clustering of the samples based on time to death was determined by k-means clustering in R. To determine whether clustering was significant, 10,000 iterations of k-means clustering was performed on the time-to-death data. This was compared to 10,000 iterations of randomly shuffling the animals’ times to death.

To determine a logistic model fit to the data, we first converted all the days after infection into a time point relative to the group mean incubation time. Thus, 0 represents when that group ended the incubation time, which was characterized by a sustained fever of greater than 1.5°C change. For the fitting, the day 3 post-exposure time point was not included because it exhibited early cytokine expression. Also, the nonresponder animals were not included because they did not develop disease. The R stats package was used with the function nls and the following formula for the fit.Embedded Image

Starting parameters were determined as follows: a is the log fold change at the end of the time course, b is the slope of a linear model fit, and c is the intercept of a linear model fit. A deviance cutoff of 20 was used to determine whether the data fit the model well (fig. S6). Plotting of all the data was done in GraphPad Prism and ggplot2.

SUPPLEMENTARY MATERIALS

www.sciencetranslationalmedicine.org/cgi/content/full/10/434/eaaq1016/DC1

Fig. S1. Comparison of animals before exposure.

Fig. S2. Comparison of MCV between groups.

Fig. S3. Correlation of cytokine expression in symptomatic animals.

Fig. S4. Example genes from NanoString in Delayed animals.

Fig. S5. Interferon response in NanoString CodeSet.

Fig. S6. Histogram of deviances from curve fitting.

Table S1. Cytokine expression in RNA-seq.

Table S2. Early expressed genes in the Delayed group.

Table S3. Normalized NanoString counts.

REFERENCES AND NOTES

Acknowledgments: We would like to thank J. J. Bearss for work as a pathologist on this study. Funding: This research was performed with NSF Graduate Research Fellowship under grant no. DGE-1247312 to E.S. J.H.C. was supported through BPS/STP-15-015 and RO1AI1096159. L.A.A., T.M., C.T.-H., C.B., and W.A. were supported through the Defense Threat Reduction Agency Project #CB-10245. S.L.B. was supported through a Chemical and Biological Defense Research Associateship award from the NRC of the National Academies and the Defense Threat Reduction Agency. This research was performed while C.E.A. held an NRC Research Associateship award at the Center for Genome Sciences, USAMRIID. The research described herein was sponsored by the Defense Threat Reduction Agency Program #CB3713 to S.L.B., A.J.G., W.D.P., F.R., J.D.S., S.E.W., and J.M.Z. Author contributions: S.L.B. and A.J.G. designed and executed the NHP study and all associated in-line analyses. W.D.P. and F.R. performed the telemetry data collection and analysis. J.D.S., S.E.W., and J.M.Z. performed the NHP study. C.T.-H., C.B., and W.A. collected NHP samples and/or performed transcriptomic assays. E.N. and G.P. performed RNA-seq. E.S., C.E.A., and G.P. performed pathway analysis. E.S. performed RNA-seq and NanoString analysis. E.S., J.H.C., and L.A.A. wrote the manuscript. S.L.B. and T.M. provided critical review of the manuscript. G.P., J.H.C., and A.J.G. made equal contributions as senior authors. E.S. and S.L.B. made equal contributions as first authors. Competing interests: The authors declare that they have no competing interests. Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the U.S. Army. Data and materials availability: The normalized NanoString count data are available in table S3. The RNA-seq data can be accessed through Gene Expression Omnibus accession number GSE103825.
View Abstract

Navigate This Article