Correlative Gene Expression to Protective Seroconversion in Rift Valley Fever Vaccinates

Abstract

Rift Valley fever Virus (RVFV), a negative-stranded RNA virus, is the etiological agent of the vector-borne zoonotic disease, Rift Valley fever (RVF). In both humans and livestock, protective immunity can be achieved through vaccination. Earlier and more recent vaccine trials in cattle and sheep demonstrated a strong neutralizing antibody and total IgG response induced by the RVF vaccine, authentic recombinant MP-12 (arMP-12). From previous work, protective immunity in sheep and cattle vaccinates normally occurs from 7 to 21 days after inoculation with arMP-12. While the serology and protective response induced by arMP-12 has been studied, little attention has been paid to the underlying molecular and genetic events occurring prior to the serologic immune response. To address this, we isolated RNA from whole blood of vaccinated calves over a time course of 21 days before and after vaccination with arMP-12. The time course RNAs were sequenced by RNASeq and bioinformatically analyzed. Our results revealed time-dependent activation or repression of numerous gene ontologies and pathways related to the vaccine induced immune response and its regulation. Additional bioinformatic analyses identified a correlative relationship between specific host immune response genes and protective immunity prior to the detection of protective serum neutralizing antibody responses. These results contribute an important proof of concept for identifying molecular and genetic components underlying the immune response to RVF vaccination and protection prior to serologic detection.

Fig 1. Serologic and perturbed gene analysis of cattle vaccinated with arMP-12.

(A) Log10 plaque reducing neutralizing titers (Log10(PRNT₈₀)) of test animals prior to and after vaccination with arMP-12. A dilution of 1:80 (log10 = 1.903) was used as the cutoff to determine protection from WT RVFV challenge based on prior published work. (B) Schematic of the immunologic status of test animals by day of experiment: unvaccinated (gray), vaccinated, below PRNT₈₀ cutoff (light red), vaccinated, first day at or above PRNT₈₀ cutoff (dark green), vaccinated, above PRNT₈₀ (light green). (C) Summary count of significantly perturbed genes normalized against unvaccinated samples occurring in at least 50% of animals by day post inoculation, with a |Z-score ≥2.24| and a false discovery rate (FDR) ≤0.05.

Fig 2. Dynamic Bayesian Gene Group Activation (DBGGA) Pathway analysis.

(A) Summary table of cellular pathways significantly activated or repressed (Bayesian score ≥|2.24|) over the time course of the experiment as determined by the DBGGA tool. (B-C) Read out and heat map of significantly perturbed pathways during days post inoculation (p.i.). Red color indicates activation, green color indicated repression. Intensity of color represents amplitude of perturbation. Analyses are oriented to the Early Phase (days 2–6 p.i.) (B) or the Late Phase (days 7, 10, 14, 21 p.i.) (C) to identify significant perturbations at different periods during the vaccine response.

Fig 3. Dynamic Bayesian Gene Group Activation (DBGGA) Gene Ontology (GO) Term analysis.

(A) Summary table of significantly perturbed GO terms and genes described by day post inoculation. Only GO terms and their genes with Bayesian score ≥|2.24| are included in analysis. (B-C) Heat maps of perturbed GO terms described by day post inoculation (p.i.) and identified from the Early Phase (B) and the Later Phase (C). Red color indicates activation, green color indicates repression. Intensity of color represents amplitude of perturbation. The list of GO terms shown represents a subset of all perturbed terms that were selected as being most relevant to innate and adaptive immune response.

Fig 4. Serologic and perturbed gene expression in time-shifted data.

(A) Log10 plaque reducing neutralizing titers (Log10(PRNT₈₀)) of test calves prior to and after vaccination with arMP-12 with data synchronized to the time point that each animal reached the minimum threshold (T = 0; 1:80, log10 = 1.903). (B) Schematic illustrating the immunologic status of test animals in time shifted analysis: unvaccinated (gray), vaccinated, below PRNT₈₀ cutoff (light red), vaccinated, first day at or above PRNT₈₀ cutoff (dark green), vaccinated, above PRNT₈₀ (light green). (C) Count of significantly perturbed genes normalized against unvaccinated samples occurring in at least 50% of animals by time pre-serum neutralization, with a |Z-score ≥2.24| and a false discovery rate (FDR) ≤0.05.

Fig 5. DBGGA pathway analysis on time-shifted data.

(A) Summary table of significantly activated or repressed cellular pathways (Bayesian score ≥|2.24|) on time shifted data. Perturbation of cellular pathways and constituent genes are reported. (B-C) Heat map representation of significantly perturbed pathways at times pre-serum neutralization (pSN). Red color indicates activation, green color indicates repression. Intensity of color represents amplitude of perturbation. Analyses are oriented to the Early Phase (time –6, –5, –4, –3 pSN) (B) or the Late Phase (time –2, –1, 0, 1 pSN) (C) to identify significant pathway perturbations at different periods during the vaccine response.

Fig 6. DBGGA GO term analysis on time-shifted data.

(A) Summary table of GO Terms and component gene perturbation by time pre-seroconversion. Only GO terms and their genes with Bayesian score ≥|2.24| are included in analysis. (B-C) Heat maps of perturbed GO terms described by time pre-seroconversion and identified from the Early Phase (time –6, –5, –4, –3 pSN) (B) and the Later Phase (time –2, –1, 0, 1 pSN) (C). Red color indicates activation, green color indicates repression. Intensity of color represents amplitude of perturbation. The list of GO terms shown represent a subset of all perturbed terms that were selected as being most relevant to innate and adaptive immune response.

Fig 7. Sliding window correlation (SWC) to identify pathways associated with serum neutralization titers.

(A) Visualization of sliding window correlation approach. (i) Hypothetical Log10(PRNT₈₀) data taken from time points capturing neutralizing antibody levels during time period in which animals reached threshold for protection (log10 = 1.903) (dark green bars). (i-ii)The trajectory of the Log10(PRNT₈₀) data is applied to other time points prior to serum neutralization (light green bars) to identify pathway with complementary or antithetical trajectories (A-ii for hypothetical Pathway Bayesian Z-score data, dark blue bars). (iii) Graphical representation of the R correlation coefficient value between Log10(PRNT₈₀) and Pathway Bayesian Z-score. (B) Pathways correlated to PRNT₈₀ values at the time of seroconversion listed by pathway at incremented window times (each consisting of three time points), with significant correlation p-values (i), or plotted as log(PRNT₈₀) vs normalized Bayesian Z-score (ii) for six selected pathways having highest correlations.

Fig 8. Sliding window correlation (SWC) to identify GO terms associated with serum neutralization titers.

(A) Immunologically important GO terms correlated with log10(PRNT₈₀) listed by time window with associated p-value indicated significance of correlation. (B) Correlated GO terms plotted as log10(PRNT₈₀) vs Normalized GO term DBGGA Bayesian z-scores. The GO term Viral life cycle had a negative correlation with seroconversion at the windowed time points –6, –5, and –4 pSN (upper panel), while Positive regulation of T cell cytokine production shows a positive correlation with seroconversion at the windowed time points –4, –3, and –2 pSN (lower panel).

Fig 9. Summary of DBGGA Pathways or GO terms and constituent genes with high correlation to protective antibody response.

(A) A selected set of DBGGA scored pathways or GO terms most correlated with protective antibody response employed for learning the dynamic Bayesian network model for inferencing future protective immune response. (B) A subset of immunologically relevant genes associated with pathways and GO terms listed in (A) that are highly correlative to serologic protection as determined by PRNT₈₀. Various time windows are represented, with time –6, –5, –4 pSN window being the furthest from the designated time of serologic protection.

Fig 10. Prediction of protective immunity by the Dynamic Bayesian Network (DBN).

(A) The Dynamic Bayesian network model representing the learned gene regulatory network structure. Each node represents a continuous variable with a Gaussian distribution. The DBN was trained with the time course data associated with each gene. Relationships between gene nodes, or between the gene node and the PRNT₈₀ node, are represented by arrows as a directed acyclic graph. (B) Example inferencing results on a selected set of gene expression data taken at time window -5, -4, -3 pSN to predict the PRNT₈₀ outcome. Graph indicates strong serologic response and log10(PRNT₈₀) levels sufficient for protection given hypothetical gene expression data. (C) Performance table for model experiment in (B) show high specificity and sensitivity of the model predictions in a time window up taken at –5, -4, -3pSN. Data sets for true negative testing were randomly selected from other genes at the same time points as the true positive gene data sets. Data from five individual animal subjects were employed for the model evaluation.