Distinct transcriptional responses to fatal Ebola virus infection in cynomolgus and rhesus macaques suggest species-specific immune responses.

Ebola virus (EBOV) is a negative single-stranded RNA virus within the Filoviridae family and the causative agent of Ebola virus disease (EVD). Nonhuman primates (NHPs), including cynomolgus and rhesus macaques, are considered the gold standard animal model to interrogate mechanisms of EBOV pathogenesis. However, despite significant genetic similarity (>90%), NHP species display different clinical presentation following EBOV infection, notably a ∼1-2 days delay in disease progression. Consequently, evaluation of therapeutics is generally conducted in rhesus macaques, whereas cynomolgus macaques are utilized to determine efficacy of preventative treatments, notably vaccines. This observation is in line with reported differences in disease severity and host responses between these two NHP following infection with simian varicella virus, influenza A and SARS-CoV-2. However, the molecular underpinnings of these differential outcomes following viral infections remain poorly defined. In this study, we compared published transcriptional profiles obtained from cynomolgus and rhesus macaques infected with the EBOV-Makona Guinea C07 using bivariate and regression analyses to elucidate differences in host responses. We report the presence of a shared core of differentially expressed genes (DEGs) reflecting EVD pathology, including aberrant inflammation, lymphopenia, and coagulopathy. However, the magnitudes of change differed between the two macaque species. These findings suggest that the differential clinical presentation of EVD in these two species is mediated by altered transcriptional responses.

Introduction

Ebola virus (EBOV) is a single-stranded RNA virus and a member of the Filoviridae family. Infection with EBOV causes Ebola virus disease (EVD), which is characterized by excessive inflammation, coagulopathy, lymphopenia, and apoptosis that ultimately results in organ failure and death [1]. In the absence of antivirals and vaccines, case fatality rates range from 40% to 90% [1-3]. The recent 2013–2016 West Africa epidemic that incurred over 28,000 reported infections and the recent outbreaks of EBOV variants in the Democratic Republic of the Congo demonstrate a continued need for ongoing research into the mechanisms of EBOV pathogenesis [3-5].

EBOV preferentially infects antigen presenting cells notably dendritic cells (DCs) and monocytes/macrophages, which are believed to be major drivers of pathology [6,7]. Infection of monocytes/macrophages contributes to the massive induction of pro-inflammatory cytokines and chemokines that recruit additional myeloid cells to the site of infection, promote neutrophil-mediated immunity, and contribute to lymphocyte apoptosis [8-10]. In contrast, infection of DCs results in reduced expression of co-stimulatory molecules and impaired antigen presentation which in turn disrupts the mobilization of adaptive immunity [10-13]. This includes T cell-mediated and humoral immunity, both which are critical for controlling infection and conferring vaccine-mediated protection [2,14,15]. Fatal EVD cases are associated with a robust and sustained adaptive systemic cytokine storm and dramatic loss of lymphocytes, while survivors show early control of innate and immune responses [16-21].

Nonhuman primates (NHPs) are the gold standard animal model for EBOV pathogenesis studies as they accurately recapitulate pathobiology of EVD observed in humans [22-24]. Historically, cynomolgus macaques have been used to study response to EBOV infection and assess vaccine efficacy, while rhesus macaques are primarily utilized in therapeutic studies [25]. Both species are equally susceptible to EBOV infection; however, differences in the duration/severity of EVD and histopathological presentation exist. Importantly, time to euthanasia following infection with a lethal dose of EBOV-Mayinga or EBOV-Makona is shorter in cynomolgus macaques (5-6 days) compared to rhesus macaques (∼8 days) [26-29]. This is accompanied by earlier changes in inflammatory mediators, earlier onset of viremia and greater severity of cellular stress and coagulopathy (e.g. elevated BUN, CRE, and AST; reduced PLT) in rhesus macaques [26-29]. Interestingly, viral titres in lymphoid organs at the time of euthanasia are comparable between the two species despite the delayed viremia [26,27]. Other studies reported small variations in the kinetics of blood and serum markers of EVD between the two macaque species although extensive comparisons have not been conducted [27,28,30]. These differences in disease progression are not unexpected given that rhesus and cynomolgus macaques exhibit stark differences following infection with influenza A virus (IAV), severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2), simian immunodeficiency virus (SIV), and simian varicella virus (SVV)[31-41]. However, no study to date has compared the host transcriptional responses to EBOV infection in these two species to elucidate molecular mechanisms explaining differences in disease pathogenesis.

To address this gap in knowledge, we compared transcriptional signatures in whole blood samples obtained from cynomolgus and rhesus macaques following infection with EBOV-Makona Guinea C07 infection. The EBOV-Makona C07 isolate was obtained early in the 2013–2016 West Africa EBOV epidemic [26,27]. Our analysis determined a shared core of antiviral, innate immune and adaptive immune genes induced in both species during infection that reflected canonical EVD pathology. In addition, we identified a more robust expression of genes that play roles in inflammation, apoptosis, coagulopathy and regulation of T lymphocytes in rhesus macaques. In contrast, expression of genes related to B cell-mediated immunity were more frequently dysregulated in cynomolgus macaques.

Materials and methods

Study cohorts

Historical whole blood RNA samples from rhesus macaques infected with Makona C07 (n = 3), as described in [26], were obtained from Dr. Marzi. Complete description of this study can be found in the original publication [26]. Briefly, whole blood samples were collected at days 0, 4 and 6 from rhesus macaques (n = 3) following intramuscular inoculation with 1,000 plaque-forming units (PFU) of EBOV-Makona C07 divided equally between both caudal thighs 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. RNA was extracted using the QIAmp viral Mini RNA kit (Qiagen) and shipped to the Messaoudi Laboratory [26]. Clinical signs of disease as well as blood chemistry values were reported in [26].

RNA sequencing data were obtained from our previous study ([29], PRJNA398558). Similar to the study with rhesus macaques, whole blood samples were collected at days 0, 4 and 6 from cynomolgus macaques (n = 5) following intramuscular inoculation with 1,000 PFU of EBOV-Makona Guinea C07 divided equally between both quadriceps. This study was conducted in the BSL4 at the University of Texas Medical Branch at Galveston (UTMB). Clinical signs of disease as well as blood chemistry values were reported in [29].

Both RML and UTMB studies obtained EBOV-Makona Guinea C07 from the same source (PMID24738640) and passaged only 1–2 times on Vero cells. Both labs confirmed the sequence of the virus before use and determined identicality with the published sequence (KJ660347).

Both RML (OLAW# A4149-01) and UTMB (IBC# 201-1581, OLAW# A3314-01) are AAALAC accredited institutions. All procedures followed standard operating procedures (SOPs) approved by the respective 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 studies were approved by the Institutional Animal Care and Use Committees (IACUC). 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 IACUC. All efforts were made to ameliorate animal welfare and minimize animal suffering in accordance with the Weatherall report on the use of nonhuman primates in research (https://royalsociety.org/policy/publications/2006/weatherall-report/).

Library generation, sequencing, and Bioinformatic analysis

Quality and quantity of whole blood RNA obtained from rhesus macaques were determined using an Agilent 2100 Bioanalyzer. cDNA libraries were constructed using the TruSeq stranded total RNA LT-LS kit. RNA was treated with RNase H and DNase I following depletion of ribosomal RNA (rRNA). Adapters were ligated to cDNA products and the subsequent ∼300 base pair (bp) amplicons were PCR-amplified and selected by size exclusion. Quality and concentration of the cDNA libraries were assessed using Agilent 2100 Bioanalyzer and subjected to 150 bp single-end sequencing using the Illumina NovaSeq platform.

Preliminary data analysis was performed with RNA-Seq workflow module of systemPipeR, developed by Backman et. al. RNA-Seq reads were demultiplexed, quality-filtered and trimmed using Trim Galore (average Phred score cut-off of 30, minimum length of 50 bp). FastQC was used to generate quality reports. Reads were aligned using hisat2 to the reference genome Macaca mulatta (Macaca_mulatta.Mmul_8.0.1.dna.toplevel.fa) and the Macaca_mulatta.Mmul_8.0.1.97.gtf was used for annotation. Raw expression values (gene-level read counts) were generated using the summarizeOverlaps function and normalized (read per kilobase of transcript per million mapped reads, rpkm) using the edgeR package. Statistical analysis with edgeR was used to determine DEGs meeting the following criteria: genes with median rpkm of ≥5, a FDR corrected p-value ≤ 0.05 and a log2fold change ≥ 1 compared to day 0.

Construction, sequencing and bioinformatics analysis of cDNA libraires from blood samples obtained from cynomolgus macaques were previously described in detail in [29]. Briefly, cDNA libraries were constructed using the TruSeq stranded total RNA LT-LS kit and subjected to 100 bp (rhesus) single-end sequencing using the Illumina NextSeq platform. This data set is available through NCBI SRA under project PRJNA398558. Analysis was carried out exactly as described above with the exception of using reference genome Macaca fascicularis (Macaca_fascicularis.Macaca_fascicularis_5.0.dna.-toplevel.fa) and annotation file Macaca_fascicularis.Macaca_fascicularis_5.0.94.gtf.

To identify commonly regulated genes during infection with adjustments for species-specific differences (Figure 3), edgeR was performed with GLM capabilities. To identify DEGs between infected rhesus macaques at d6 and infected cynomolgus macaques at d4, the GLM approach was used with makeContrasts function to account for average species-related effects. Sparse partial least squares discrimination analysis (SPLS-DA) was performed using the mixOmics R package for classification and validation. Models were initially built with 10 components and validated with three-fold cross-validation and 10 repetitions.

Figure 3.

Cynomolgus and rhesus macaques share a core of EVD-related genes. (A) Volcano plot of global gene expression changes shared between infected cynomolgus and rhesus macaques. DEGs (average rpkm ≥5) are denoted in red. Exemplar DEGs are labelled. (B) GO network depicting functional enrichment of DEGs expressed by both infected rhesus and cynomolgus using Metascape. Clustered nodes of identical colour correspond to one GO term. Node size represents the number of DEGs associated with the GO term. Grey lines represent shared interactions between GO terms, with density and number indicating the strengths of connections between closely related GO terms. Heatmaps representing DEGs enriching to GO terms (C) “response to virus,” “regulation of production of type I interferon,” “interferon-gamma-mediated signaling pathway,” (D) “lymphocyte activation,” (E) “myeloid cell activation,” “myeloid cell differentiation,” (F) “apoptotic signaling pathway” depicted in part B. Where GO terms consisted of more than 60 DEGs, only 60 are represented. Each column represents the median of the normalized rpkm of samples. Range of colours is based on scaled and centred rpkm values of the represented DEGs. Red represents upregulated; blue represents downregulated. DEGs expressed to a greater extent in cynomolgus (orange) or rhesus (green) macaques are coloured (two-tailed T test).

Functional enrichment of DEGs was performed using Metascape to identify relevant GO biological process terms. Heatmaps, bubbleplots, Venn diagrams and violin plots were generated using R packages ggplot2 and VennDiagrams. GO network plots were generated in Cytoscape (Version 3.5.1). Graphs were generated using GraphPad Prism software (version 8).

Statistical analysis

Differences in transcript abundance between the two macaque species was determined with a two-tail independent T test using the genefilter R package. PERMANOVA was performed using the vegan R package to assess contributions of variance from each principal component in Figure S1. *p-value≤0.05, **p-value≤0.01, ***p-value≤0.001.

Click here for additional data file.