All the work involving infectious MARV was performed following standard operating procedures (SOPs) approved by the Rocky Mountain Laboratories (RML) Institutional Biosafety Committee (IBC) in the maximum containment laboratory at the RML, Division of Intramural Research, National Institute of Allergy and Infectious Diseases, National Institutes of Health. 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 United States Department of Agriculture and was approved by the RML Animal Care and Use Committee (ACUC). Procedures were conducted by trained personnel under the supervision of veterinary staff on animals anesthetized with ketamine. 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/). Animals were housed in adjoining individual primate cages that enabled social interactions, under controlled conditions of humidity, temperature, and light (12-h light:12-h dark cycles). Food and water were available ad libitum. Animals were monitored and fed commercial monkey chow, treats, and fruit at least twice a day by trained personnel. Environmental enrichment consisted of commercial toys, music, and video. Endpoint criteria based on clinical score parameters as specified and approved by the RML ACUC were used to determine when animals were humanely euthanized.

Vero E6 cells (mycoplasma negative) were grown at 37°C and 5% CO₂ in Dulbecco’s modified Eagle’s medium (DMEM) (Sigma-Aldrich, St. Louis, MO) containing 10% fetal bovine serum (FBS) (Wisent Inc., St. Bruno, Canada), 2 mM L-glutamine, 50 U/mL penicillin, and 50 μg/mL streptomycin (all supplements from Thermo Fisher Scientific, Waltham, MA). VSV vaccines and challenge virus were propagated in Vero E6 cells using DMEM supplemented with 2% FBS, L-glutamine and penicillin/streptomycin. VSV-MARV expressing the MARV-Angola glycoprotein (GP) (8) and VSV-EBOV expressing the EBOV-Mayinga GP (15) were used for IM vaccination. MARV-Angola was obtained from the Public Health Agency of Canada (16) (GenBank accession number KY047763 (17)), propagated on Vero E6 cells, titered and stored in liquid nitrogen. All viruses were confirmed by sequencing and mycoplasma testing revealed no contaminants.

Sixteen male and female cynomolgus macaques (Macaca fascicularis) 5-7 years of age and 4.0-6.7 kg in weight were used for this study. Three groups of cynomolgus macaques (n=4 per group) were vaccinated with a single intramuscular (IM) injection of 1x10⁷ plaque forming units (PFU) VSV-MARV at 14, 7 or 3 days before challenge. Control animals were IM-vaccinated with 1x10⁷ PFU VSV-EBOV (n=1 each day -14 and -7, n=2 day -3; n=4 total). Challenge was performed by IM injection of 1,000 PFU of MARV (confirmed by back-titration) as previously described (8). Clinical exams including a blood draw on anesthetized animals were conducted on days -14, -7, -3, 0, 3, 6, 9, 14, 21, 28, 35 and 42 post-challenge. The animals were observed at least twice daily for clinical signs of disease according to a RML ACUC-approved scoring sheet and humanely euthanized when they reached endpoint criteria. The study ended on day 42 post-challenge when all surviving animals were humanely euthanized.

The total white blood cell, neutrophil, lymphocyte, and platelet counts were determined from EDTA blood with the IDEXX ProCyte DX analyzer (IDEXX Laboratories, Westbrook, ME). Serum biochemistry including aspartate aminotransferase (AST), alkaline phosphatase (ALP), glucose, creatinine, and total bilirubin was analyzed using the Piccolo Xpress Chemistry Analyzer and Piccolo General Chemistry 13 Panel discs (Abaxis, Union City, CA).

Viremia was determined from EDTA WB samples using Vero E6 cells (mycoplasma negative). To this end, cells were seeded in 48-well plates the day before titration. At the day of titration, blood samples were thawed, and 10-fold serial dilutions were prepared in DMEM without supplements. Media was removed from cells and triplicate wells were inoculated with each dilution. After one hour, DMEM supplemented with 2% FBS, penicillin/streptomycin and L-glutamine was added and cells were incubated at 37°C. Cells were monitored for cytopathic effect (CPE) and 50% tissue culture infectious dose (TCID₅₀) was calculated for each sample employing the Reed and Muench method (18).

Post-challenge NHP sera were inactivated by gamma irradiation (5 MRad) (19) and removed from the maximum containment laboratory according to RML standard operating procedures (SOP) approved by the RML IBC. Antibody titers were determined using ELISA kits based on recombinant soluble MARV-Angola GPδTM (Alpha Diagnostics, San Antonio, TX). Serum samples were diluted 1:200 for assessment in the ELISA. The ELISA was performed, and titers were calculated as per manufacturer’s instructions.

Serum samples were diluted 1:2 in serum matrix for analysis using the Milliplex Non-Human Primate Magnetic Bead Panel as per manufacturer’s instructions (Millipore, Burlington, MA). Concentrations for IL-1Ra, IL-8, IL-10, MCP-1, IFN-γ, IL-6, MIP-1α, IL-2 and IL-15 were determined for all samples. Values below the limit of detection of the assay were assigned the value of 1.

RNA quality and concentration were assessed using an Agilent 2100 bioanalyzer. A total of 100 ng was used for library preparation. Total RNA samples were treated with 1 U of DNase I (New England BioLabs, Ipswich, MA) at 37°C for 10 min and cleaned with RNAClean XP beads (Beckman Coulter, Brea, CA). Ribosomal RNAs were depleted using rRNA removal beads (Illumina, San Diego, CA). The remaining RNA was fragmented at 94°C for 8 min to yield a median fragment size of 155 nucleotides (nt). Libraries were prepared using the TruSeq stranded total RNA library prep kit according to the instructions provided by the manufacturer (Illumina, San Diego, CA). Final library quality was assessed using an Agilent high-sensitivity DNA kit. Libraries with unique barcoded adaptors were pooled and sequenced in a single end read (1x100) on the Illumina NextSeq 500 platform.

Data analysis was performed using the RNA-Seq workflow module of the systemPipeR package (20). RNA-Seq reads were demultiplexed, quality filtered, and trimmed using Trimmomatic with an average phred score cutoff of 30 and minimum length of 50 bp. Quality reports were generated with the FastQC function. Trimmed reads were mapped to the Macaca fascicularis reference genome (Macaca_fascicularis.Macaca_fascicularis_5.0.dna.toplevel.fa) using HISAT2 and the corresponding gene annotation (Macaca_fascicularis.Macaca_fascicularis_5.0.94.gtf) from Ensembl. Uniquely mapped reads were counted using summarizeOverlaps in strand-specific mode. Normalization and statistical validation of differentially expressed genes (DEGs) was performed using EdgeR package’s pair-wise function (21). 0 DPV data were used as the reference. DEGs were defined as protein coding genes with human homologues with at least a 2-fold change in expression, a multiple hypothesis Benjamini-Hochberg false discovery rate (FDR) corrected P value less than 0.05 and an average of at least 5 read per kilobase of transcript per million mapped reads (RPKM).

Temporal gene expression patterns and signatures that distinguish vaccine groups, survivors and non-surviving, and negative controls were analyzed using maSigPro, which is a two-way regression-based approach that finds a set of statistically significant DEGs for the entire time course (22). Only protein-coding genes with human homologs and an average of at least 5 read per kilobase of transcript per million mapped reads (RPKM) were included in this analysis.

DEGs were first mapped to human homologs using BioMart (Ensemble Gene 94). Only protein-coding genes with human homologs were included for further analysis. The functional enrichment of DEGs was assessed using Metascape (23).

Heatmaps, Venn diagrams, bar graph and volcano were generated using R packages VennDiagram, dplyr, and ggplot2. Line graphs were generated using GraphPad Prism V8 (San Diego, CA).

Clinical data were examined for statistical significance using Prism version 9 (GraphPad, San Diego, CA). Survival curves were analyzed with Mantel-Cox test and values representing groups were analyzed by two-way ANOVA with Tukey’s multiple comparisons. Statistical analysis of the maSigPro data was carried out using Prism version 8 (GraphPad, San Diego, CA). Significance was determined using a one-way ANOVA with a Dunnett’s multiple-comparison test. Statistically significant differences are indicated as follows: p<0.0001 (****), p<0.001 (***), p<0.01 (**), and p<0.05 (*).

The minimum time between vaccination and challenge for protection with a single IM dose of 1x 10⁷ PFU VSV-MARV was determined by vaccinating groups of 4 NHPs at 14, 7 and 3 days prior to challenge. The control group consisted of 4 VSV-EBOV-vaccinated NHPs that were vaccinated with a single IM dose of 1x 10⁷ PFU on day -14 (n=1), day -7 (n=1) and day -3 (n=2). On day 0, all NHPs were challenged with 1,000 PFU MARV by IM injection. The control animals developed signs of MVD and were humanely euthanized 6- and 8-days post-challenge (DPC) when they reached IACUC-approved endpoint criteria ( Figure 1A ). One NHP in the day -3 vaccine group developed clinical signs of MVD as demonstrated by the increase in clinical score ( Figure 1B ) and was euthanized 7 DPC ( Figure 1A ). Of the remaining three animals in the day -3 group, one developed moderate signs of MVD, two developed very mild signs of MVD, and all three recovered ( Figure 1B ). None of the NHPs in the day -14 or -7 vaccine groups developed signs of MVD ( Figure 1B ). Only the 5 NHPs that succumbed to disease developed hallmarks of MVD including thrombocytopenia ( Figure 1C ), high titer viremia ( Figure 1D ), and increased levels of AST ( Figure 1E ) and ALP ( Figure 1F ). Other parameters examined in WB and serum after challenge demonstrated changes to abnormal levels for cell populations and metabolites for the 5 NHPs that succumbed to MVD ( Figure S1 ). Similarly, only these 5 NHPs developed cytokine levels indicative of the disease-associated cytokine storm ( Figure S2 ).

VSV-based vaccines have been shown to primarily confer protection by humoral responses (15, 24). Therefore, antigen-specific immune responses were determined in serum samples collected from all NHPs throughout the study. Levels of MARV GP-specific IgM peaked in the day -14 group 11 days post-vaccination (DPV; -3 DPC) before dropping by DPC 0 ( Figure 2A ). MARV GP-specific IgM was detected in the day -7 group 3 DPC increasing slightly 6 DPC dropping 9 and 14 DPC ( Figure 2A ). In the day -3 group MARV-GP specific IgM was detected 6 DPC increasing slightly 9 and 14 DPC ( Figure 2A ). Levels of MARV GP-specific IgG were detected 11 DPV (day -14 group), 3 DPC (day -7 group) and 9 DPC (day -3 group) ( Figure 2B ) coinciding with peak levels of MARV GP-specific IgM ( Figure 2A ). For the day -14 vaccination group the MARV challenge served as a boost as demonstrated by the increase in MARV GP-specific IgG 6 and 9 DPC ( Figure 2B ).

Longitudinal transcriptional changes in WB in response to VSV-MARV and VSV-EBOV vaccination were assessed using RNA sequencing (RNA-Seq). To that end, transcriptomes of all available WB samples collected on 3, 7, 11, and 14 DPV were compared to those observed 0 DPV. Gene expression changes were evident as early as 3 DPV with VSV-MARV with 88, 73, 483 and 828 DEGs detected 3, 7, 11, and 14 DPV, respectively ( Figure 3A ). The largest number of DEGs was detected 14 DPV ( Figure 3A ). DEGs enriching to gene ontology (GP) terms associated with innate immunity to viral infection such as “defense response to virus” and “response to interferon alpha” were detected 3 DPV and remained differentially expressed throughout 14 DPV ( Figure 3B ). Notable DEGs in this group include interferon stimulated genes (ISG; ISG15, BST2, IRF7, IFIH1) as well as genes that play a role in antigen presentation (HLA-B), cytokine production (MIF) and antiviral response (RSAD2, CARD9, DDX58, DDX60, and DHX58) ( Figure 3C ). DEGs important for lymphocyte activation and adaptive immunity (e.g., GO terms “lymphocyte activation” and “lymphocyte mediated immunity”) were detected 11 DPV and included genes important for T cell activation and effector functions (CD3D, CD8B, IL2RA, IL2RB1, GZMB and CD274). DEGs important for B cell activation and antibody secretion (MZB1, BATF, IGHM, IGHA1) were upregulated 14 DPV ( Figures 3B, C ).

Samples collected 3, 7, and 11/14 DPV following VSV-EBOV vaccination were also compared to 0 DPV. Samples collected 11 and 14 DPV were combined because of the small number of animals. We detected 179, 10, and 576 DEGs at 3, 7, and 11/14 DPV, respectively ( Figure 3D ). As described for VSV-MARV above, DEGs important for innate immune response to virus were detected as early as 3 DPV and enriched to GO terms “type I interferon signaling pathway” and “defense response to virus” ( Figure 3E ). These DEGs included genes important for anti-viral defense (OAS1, OAS2, OASL, MX1); interferon stimulated genes (IFIT2, IFIT3, IFI44, ISG15); and inflammation (S100A8, S100A9, CXCL10) ( Figure 1F ). DEGs associated with lymphocyte activation and B cell proliferation were detected 11/14 DPV (GO terms “regulation of B cell proliferation” and “lymphocyte proliferation”) ( Figure 3E ) including ARG2, CD40, IL12B, TNFSF13B and TNFRSF13C ( Figure 3F ).

Next, we examined gene expression changes post MARV challenge in the day -14 and -7 groups. All comparisons were made relative to 0 DPV in a pairwise manner in each group. Limited transcriptional changes were detected with the exception of 9 DPC in the day -7 group ( Table 1 ). DEGs detected in the day -14 group at 3 and 6 DPC mapped to GO terms associated with viral life cycle, leukocyte migration, and lymphocyte activation ( Figure 4A ). DEGs important for lymphocyte activation include genes involved in T cell activation (CD3E, decreased 9 DPC), and B cell activation and antibody secretion (IGLL1, IGHG2, IGHD, and IGHA1; increased 9 DPC) ( Figure 4B ). DEGs detected 14-42 DPC played a role in “oxidative phosphorylation” and “ATP metabolic process” including (ATP5F1D, ATP5F1E, COX2, NDUFA2) ( Figure 4A ).

As described above, transcriptional changes detected in the day -7 group following MARV challenge enriched to GO processes associated with viral life cycle, myeloid cell differentiation, and apoptotic signaling pathway ( Figure 4C ). DEG that enriched to GO term “viral life cycle” included genes that play a role in viral budding (VPS28, VPS37B, GRK2), and Interferon-stimulated genes and chemoattractant (ISG15, ISG20, CXCL8, CCL5). DEGs important for adaptive immunity are detected 9 DPC with DEGs enriching to “T cell activation”, and “B cell activation” including DEGs that played a role in antigen presentation as well as T and B cell activation (CD74, HLA-DRB5, HLA-C, CD247, CD19, CD79A, IGHA1, and CD27) ( Figure 4D ). Interestingly, DEGs enriching to GO terms associated with lymphocyte activation were down regulated 9 DPC ( Figure 4D ). Finally, DEGs detected 14-42 DPC played a role in “myeloid cell differentiation” including ITGA2B, SLC4A1, STAT5B, H3C11, and EPB42 ( Figure 4C ).

To better understand the temporal patterns of gene expression changes and signatures that distinguish the day -7 group, day -14 group, and negative controls, we used maSigPro (22). This approach provides a set of statistically significant DEGs for the entire time course rather than at each time point. This analysis identified 4 clusters with temporal expression changes that were significantly different between the groups ( Figure 4E ). Expression of genes in cluster 1 (1088 genes) either decreased 3 and 6 DPC (day -14 group) or remained unchanged (day-7 group) in vaccinated animals but increased in the controls 6 and 9 DPC ( Figure 4E ). Functional enrichment indicated that these genes were associated with inflammatory processes such as “response to wounding” (e.g., F2, SERPINC1, APOH, ARRB1) and “cytokine metabolic process” (e.g., GDNF, IL34, VEGFC, CXCL6) ( Figures 4F , S3A ). Genes in cluster 3 (320) had similar expression patterns as genes in cluster 1 (unchanged in day -7 group, decreased at 3 DPC in the day -14 group, significantly increased in controls). These genes enriched to “lymphocyte apoptotic process” (FAS, BTK, CRKL, LGALS9, LGALS14, STAT3) and “regulation of cell adhesion” (ITGA5, PAK1, ADAM19) ( Figures 4F and S3A ). Expression of genes in cluster 2 (905) increased 9 DPC in both day -7 group and controls but decreased in in day -14 group 3 DPC. These genes mapped to GO terms “response to virus” (e.g., DHX9, IRF7, MX2, OAS2, STAT2, TLR8), “autophagy” (e.g., MAPK8, CASP1, ANXA7, VPS13A, VPS13C, VPS35), and “leukocyte activation” (e.g., CD44, CD53 CD59, HLA-DRA, CD46, FCER1G) ( Figures 4F and S3A ). Finally, expression of genes in cluster 4 (325) significantly decreased 3 DPC in the day -14 group, 6 and 9 DPC in the controls, and 9 DPC in the day -7 group. Those genes enriched to GO terms associated to “lymphocyte activation” (CD3D, CD27, IGHD, IGLC1, CSK, TNFRSF1B) and “antigen processing and presentation” (HLA-B, HLA-C, HLA-DMB, HLA-E, CTSB) ( Figures 4F and S3A ).

VSV-MARV vaccination 3 days before challenge led to a mixed outcome with 3 animals surviving and 1 animal succumbing to MARV infection. We first carried out pairwise analysis of the transcriptional profiles of the samples in animals that survived challenge relative to 0 DPV. DEGs were detected at each DPC relative to 0 DPV peaking at 6 DPC ( Table 1 ). DEGs that enriched to GO terms “viral transcription” (e.g., POLR2E, POLR2F, RPL8, RPS27A, and RPS28); “myeloid leukocyte activation” (SERPINB1, MYD88, CD14, TLR4, and TLR8); and several other immune related GO terms (MX1, OAS1, OAS2, STAT2, S100A8, S100A12, STAT1, HERC5, and TGFB1) were upregulated at 3 and 6 DPC ( Figures 5A and S3B, C ). Genes important for antigen processing and presentation, and T and B cell activation were detected 6, 9, and 14 DPC, respectively ( Figure 5A ). DEGs involved in “antigen processing and presentation” included MHC/HLA molecules (HLA-B, HLA-C, HLA-DQB1) and components of the proteosome (PSMA3) ( Figure S3C ). Genes encoding HLA molecules were mostly downregulated while genes encoding proteosome subunits were upregulated. DEGs mapping to “T cell activation” and “B cell activation” include CD4, CD8B, CCL5, TNFSF13B, CD274, IL2RA, CD19, CD27, IGHA1, and IGHD the expression of which peaked 9 DPC ( Figure 5B ).

To uncover gene expression signatures that distinguish survival versus fatal outcome, we analyzed transcriptional changes 0 DPV, 0, 3, 6, and 9 DPC using maSigPro with all NHPs vaccinated on day -3. Animals were divided into survivors (n=3), non-survivor (n=1), and controls (n=2) for the analysis. This approach identified 4 clusters of genes with temporal expression changes that were significantly different between the groups ( Figure 5C ). Genes in clusters 1 (277), 2 (240), and 4 (109) were upregulated in both the negative controls as well as the animal that succumbed to disease. Genes in cluster 1 play a role in natural killer cell activation (e.g. VAMP7, SLAMF7, and SNX27) and integrin activation (RAP1B, ITGB1BP1, and P2RY12) ( Figures 5C, D and S3D ). Genes in cluster 2 genes enriched to cell cycle G2/M phase transition (PSMB5, PSMC1, KHDRBS1, and CCNB1); antigen processing and presentation (KLC1 and DCTN2); and signaling (MAP3K20, KIDINS220, PARK7) ( Figures 5C, D and S3D ). Cluster 4 genes mapped to GO terms including neutrophil activation (PNP, CD63, and FCAR); positive regulation of phagocytosis (FCER1G, FPR2, and PTX3); and regulation of innate immune response (IRF7, TRIM21, IL18RAP, and SOCS3) ( Figures 5C, D and S3D ). Genes in cluster 3 were significantly upregulated only in the negative control animals at 9 DPC. Those genes enriched to cell cycle checkpoints (CDK1 and CDK5RAP2) and translation (EIF2D and EIF4B) ( Figures 5C, D and S3D ).

To further capture any transcriptional differences between the 3 surviving animals and the animal that succumbed to MARV challenge, we compared the WB transcriptional landscape between these animals at 0 DPV. This comparison identified 55 upregulated genes in the surviving relative to non-surviving animal on the day of vaccination some of which play a role in myeloid leukocyte activation (CD14, GRN, HEXIM1, and SIGLEC5); regulation of MAPK cascade (INPP5K, CEACAM1, and MINK1); and regulation of cell adhesion (ST3GAL4, ADAM8, and LAMA5) ( Figures S3E, F ).