1. Introduction

Vaccination programs are one of the most effective means of controlling infectious diseases and with the development of oral vaccines and bait delivery systems, the elimination of diseases circulating in wildlife populations has become a realistic possibility. The large-scale oral rabies vaccination (ORV) campaigns that have eliminated fox-mediated rabies from Western Europe and North America and substantially reduced disease incidence in central Europe are pre-eminent examples for the success of such control programs.

ORV programs in foxes are aimed at increasing herd immunity in the target population using oral rabies vaccines distributed into the environment. Over the past four decades several oral rabies vaccines, mainly live replication-competent attenuated rabies virus vaccines, have been successfully used in ORV campaigns. In Canada for example, the ERA-BHK21vaccine virus strain, a derivative of the cell culture adapted vaccine virus strain Street Alabama Dufferin (SAD), was the only live attenuated vaccine deployed in fox ORV campaigns. In Europe, with the exception of a recombinant vaccinia virus expressing the rabies virus glycoprotein, all constructs have been based on live attenuated rabies virus strains, derived from the SAD Bern original (SAD Bernorig) vaccine virus strain, a successor of the ERA strain. While all these vaccines have been highly efficient in fox rabies control, the first generation of SAD-derived vaccines demonstrated residual pathogenicity in non-target species particularly in rodents. Although several cases of vaccine virus-induced rabies were observed even in species other than rodents over the course of vaccination campaigns in a number of countries, including Germany, Austria, Slovenia, Romania, Poland, and Canada, such cases were without epidemiological relevance.

While previous analyses using high-throughput sequencing approaches revealed substantial genetic heterogeneity within commercial SAD-derived oral rabies virus vaccines, the sub-consensus genetic heterogeneity of viruses isolated from these vaccine virus-induced rabies cases on the contrary revealed nearly clonal genotypes, indicating the presence of a strong bottleneck during infection.

In this study, we attempted to further analyze and elucidate the mechanisms of genetic selection using a combined deep sequencing, variant, and haplotype analysis of a large set of vaccine-induced rabies cases that included additional SAD- and ERA-induced rabies cases from Poland and Canada alongside vaccine virus batches. With this comprehensive approach, we could demonstrate the utility of this approach for identity and genetic stability (revision to virulence) testing of vaccines with a heterogeneous genetic background. Additionally, we were interested to know whether viruses from vaccine-induced rabies cases differ from virulent field rabies viruses (RABV) at a population level.

2.1. Samples

In this study, original brain material of additional ERA-BHK21 and SAD Bernorig-derived vaccine-induced rabies cases from Canada and Poland, respectively, were analyzed (Table 1). In addition, new batches of the attenuated rabies virus vaccine SAD Bern (Lysvulpen) were obtained and included in the analysis (Table 2). In order to check for validity of the PCR based and Illumina generated sequence reads for the use of the methods to be applied, we included datasets of two case samples (C/CAN/1994/MB and C/CAN/1996/B) as well as a dataset of the ERA/lot16 vaccine batch already sequenced on an IonTorrent PGM platform from a previous study.

2.2. Nucleic Acid Extraction, Sample Processing, and Sequencing

Samples were processed according to two different methodical approaches. Methodology 1 corresponds to the published protocol of Wylezich and colleagues. Briefly, for sample disintegration prior to nucleic acid extraction, 20 mg tissue material or 200 µL of the vaccine virus preparation from the blister were cryofractured using the cryoPREP CP02 (Covaris, Brighton, United Kingdom). The pulverized material was lysed with pre-heated (56 °C) buffer AL (Qiagen, Hilden, Germany). A threefold volume of Trizol was added and RNA extracted using the RNeasy Mini Kit (Qiagen) followed by on-column DNase I digestion (Qiagen). Subsequently, the obtained RNA was converted into double stranded cDNA using a combination of SuperScript™ IV First-Strand cDNA Synthesis System (Invitrogen/Thermo Fisher Scientific, Waltham, MA, USA) and NEBNext® Ultra™ II Non-Directional RNA Second Strand Synthesis Module (New England Biolabs, Ipswich, MA, USA). The generated cDNA was fragmented using the Covaris M220 (Covaris) and subsequently converted to Ion Torrent compatible libraries using the GeneRead L Core Kit (Qiagen) and IonXpress Barcode Adaptors (Thermo Fisher Scientific, Waltham, MA, USA) followed by size selection as described. After quality control (Agilent 2100 Bioanalyzer; High Sensitivity DNA Kit, Agilent Technologies, Santa Clara, CA, USA) and quantification (KAPA Library Quantification Kit—Ion Torrent PGM Uni; KAPA Biosystems/Roche, Basel, Switzerland), the libraries including the Ion S5 Calibration Standard were sequenced on an Ion 530 chip with an Ion S5 XL System (Thermo Fisher Scientific) in 400 bp-mode according to the manufacturer’s instructions.

Methodology 2 followed sample preparation as described before. In short, RNA was prepared using a protocol combining Trizol and a semi-automated MagMax™ (Applied Biosystems/Thermo Fisher Scientific, Waltham, MA, USA) extraction system. Subsequently, dsDNA for input into library preparation was generated by reverse-transcription polymerase chain reaction (RT-PCR), yielding overlapping amplicons of the RABV genome (Tables S1 and S2, Supplementary File). The generated amplicons were individually quantified for each sample utilizing a Qubit fluorometer (ThermoFisher Scientific) and pooled equimolarly followed by library generation using the Nextera XT kit (Illumina, San Diego, CA, USA). Sequencing was performed with an Illumina MiSeq instrument generating 250 base paired end reads.

2.3. Population-Based Analysis, Sequence Assembly, and In-Depth Variant Analysis

Data analysis was performed essentially as described before with the respective published data included (Tables S3–S6, Supplementary File). After trimming for quality of reads (all datasets), the dataset sizes were adjusted to match the size of the published datasets. All sequence reads were mapped along a strict majority consensus sequence and its reverse complement (454 Sequencing System Software v3.0; Roche) that was previously generated and introduced via a MAFFT (multiple sequence alignment using fast Fourier transform) alignment of published SAD-related sequences. Using R (version 3.6.0) and RStudio (version 1.2.1335) in combination with additional packages (Table 3), base frequencies of the forward and reverse mapping were calculated for each position of the reference sequence. Subsequently, the data were corrected for strand-bias and pairwise Manhattan distances were calculated for each population. After non-linear multidimensional scaling, the distances were displayed in 2-dimensional plots. For reasons of simplification, previously sequenced in vitro selected attenuated rabies vaccine strains, i.e., P5/88, VA1, SAG2, SAD B19CS, and SAD B19P1 (Table S3, Supplementary File), were excluded to achieve a more convenient presentation of the SAD Bernorig-derived vaccine-induced rabies cases analyzed in this study.

Full genome sequences of all samples listed in Table 1 and Table 2 were obtained by de-novo assembly of full or partial matching reads from the forerun mapping (454 Sequencing System Software v3.0; Roche) and submitted to the European Nucleotide Archive (ENA) under the study accession PRJEB35810.

For in-depth variant analysis, strand-bias corrected data of base frequencies were used as the basis for comparison of single nucleotide variants between vaccine virus strains and their induced rabies cases. Here, all nucleotide variants with a variant frequency of at least five percent (Table S7, Supplementary File) were considered. For sample datasets generated with both IonTorrent and Illumina (C/CAN/1994/MB—C/CAN/1994N6762M, C/CAN/1996/B—C/CAN/1996N5367, and ERA/lot16), variants were taken into account if they were present in both datasets and their mean variant frequency was equal or above five percent. In parallel, fully or partially mapped sequence reads were analyzed by Geneious Prime (2019.2.3; build 2019-09-24) for an additional variation analysis using standard settings and a minimum variant frequency threshold of 0.05 to confirm positions and frequencies of single nucleotide variants found in the strand-bias corrected dataset. Noteworthy, the populations determined for samples C/CAN/1994/MB—C/CAN/1994N6762M and C/CAN/1996/B—C/CAN/1996N5367 from Ion Torrent and Illumina data were nearly identical, putting the related viruses in close proximity in the distance analysis.

3.1. Viruses from Vaccines Form Product Specific Clusters in Population Analysis

Overall, the distance plot generated from the population data segregates the different oral rabies vaccine viruses into two major clusters, i.e., ERA-BHK21- and SAD Bernorig-derived vaccines (Figure 1), confirming previous observations. When additional SAD Bern (Lysvulpen) batches were included, the SAD Bernorig-derived cluster appeared to form two sub-clusters. While the majority of SAD B19 and SAD Bern vaccine batches clearly separated into these sub-clusters, one batch of each group was dislocated. Namely, batch B/SAD/B19/958 clustered with SAD Bern (Lysvulpen) vaccine batches, whereas vaccine batch B/SAD/Bern/0213 was shifted to the vaccine-induced rabies cases (Figure 1) as previously described.

Population analysis of the viruses isolated from vaccine-induced rabies cases revealed two distinct patterns. Those viruses that were associated with SAD Bernorig-derived oral rabies vaccines (including the newly analyzed Polish case C/POL/2017/B) did not cluster with the related vaccines as shown before. In contrast, the viruses from the ERA vaccine-induced rabies cases all clustered closely together with the ERA-BHK21 vaccine batch (B/ERA/lot16, Figure 1). To further elucidate the observed differences, we conducted an in-depth analysis of the relations between viruses from SAD Bernorig- and ERA vaccine-induced rabies cases to their progenitor vaccines.

3.2. ERA-Derived Subclusters Comprise Unique Patterns of Base Exchanges

The ERA-BHK21 vaccine virus population displayed a high genetic heterogeneity represented by diverse single nucleotide variants distributed across the N, P, G, and L genes (Figure 2), supporting previous findings. Interestingly, the in-depth variant analysis revealed distinct dependencies of the ERA vaccine-induced cases from vaccine variants, with three specific combinations of single nucleotide variants (SNVs) being selected from the progenitor population at the consensus level (Figure 2). Most of these differences were derived from single nucleotide variants (SNVs) already present in the progenitor ERA-BHK21 vaccine virus strain, while a few additional mutations at the consensus level were of spontaneous origin. While all but one of the latter mutations appeared to be unique, the selected SNVs could be found across several viruses from ERA vaccine-induced rabies cases. Thus, the analyses of viruses from the ERA vaccine-induced cases confirmed distinct combinations of SNVs. According to these combinations, three groups of viruses from ERA vaccine-induced rabies cases could be defined (Table 4, Figure 2).

Each of these groups was characterized by a unique set of differences at the consensus level ranging between three and nine genome positions that were not present in any of the other groups (Table 4, Figure 2). This suggests the presence of at least three replication-competent virus haplotypes in the progenitor ERA-BHK21 vaccine virus strain that were selected in vivo by a genetic bottleneck in affected animals. However, although the limited number of ERA vaccine-induced rabies cases hampers a statistical analysis, there seems to be no spatio-temporal nor species-specific correlation among the three groups of viruses from ERA vaccine-induced rabies cases.

With 8–14 differences at the consensus level, sequence diversity of the ERA vaccine-induced rabies cases and their progenitor vaccine virus strain was relatively high compared to the viruses of the SAD Bernorig rabies cases. This observation seemingly contradicts results from the population-based analysis where the ERA vaccine-induced cases clustered closely together with the ERA-BHK21 vaccine batch. However, differences at the consensus level of the viruses of ERA vaccine-induced rabies cases were already included in the progenitor vaccine population (ERA-BHK21) and likewise, in some variants, the former consensus sequence was present as a minor variant. For instance, in each of two ERA vaccine-induced rabies cases (C/CAN/1992N8944 and C/CAN/1994N6762M) one additional single nucleotide variant (A10840[G,A] and C3085[T,C]) appeared in the viral populations; in both cases, the progenitor consensus was maintained as a minor variant (Table 5).

3.3. SAD Bernorig-Derived Cases are Defined by High Sequence Identities Relative to Their Related Vaccines

Strikingly, in contrast to the ERA vaccine-induced rabies cases, for the SAD Bernorig-induced rabies cases no SNVs already present in the progenitor vaccines were selected, even though a higher number of variants were detected in the SAD Bernorig-derived vaccine virus strain (Nmean ≈ 47, see Table S8, Supplementary File) than in the ERA-BHK21 strain (N = 31, see Figure 2 and Table S7). Also, the consensus sequences of viruses from the SAD Bernorig vaccine-induced rabies cases showed a high sequence identity to their progenitor vaccine virus strain at the consensus level, ranging from zero to four nucleotide exchanges, confirming previous studies.

3.4. Known Antigenic Sites Are Rarely Affected by Amino Acid Exchanges in Vaccine-Induced Cases

To elucidate whether the induction of disease in these vaccine-induced rabies cases was a result of mutations at known antigenic sites, all sequences were screened for nucleotide exchanges at the respective genome positions. The analysis revealed no nucleotide exchanges between the SAD Bernorig-derived vaccine virus strains and their induced rabies cases in any of the known antigenic sites, including SNVs on these particular positions. Furthermore, none of the antigenic sites were affected by mutations in the ERA vaccine-induced rabies cases, except for one amino acid (AA) exchange Thr178Ile at an interferon-antagonist motif in the P-gene found in all members of group 2 (Table 4). Whether this single AA exchange has any influence on pathogenicity remains speculative, as the entire motif contains eleven AAs and only the deletion of segments or the complete motif was studied. On the other hand, if pathogenicity is increased by this presumed haplotype of group 2, more vaccine-induced rabies cases of this variant would be expected.

3.5. Viruses from Vaccine-Induced Cases Can Clearly Be Distinguished from Field RABV Viruses

When selected field RABV viruses were included in the Manhattan distance matrix analysis, the attenuated vaccine virus strains and viruses from their vaccine-induced rabies cases were clearly separated from them by high distances (Figure 3). Furthermore, dimensions of the plot increased substantially compared to Figure 1, resulting in the condensation of minor differences and the formation of a single tight cluster.

4. Discussion

Although various studies have focused on oral rabies vaccine-induced cases, only a few studies have comprehensively investigated those cases and their progenitor vaccine viruses at a population level. The latter showed a high genetic heterogeneity in different vaccine virus strains as opposed to viruses from vaccine-induced rabies cases, indicating a substantial loss in viral population diversity. In this study, we could identify two different vaccine virus-dependent patterns in the genetic markup of viruses from vaccine-induced rabies cases. This was only possible by including further samples of ERA and SAD Bern vaccine-induced rabies cases in the in-depth variant analysis (Table 1). Hence, in total 20 datasets from viruses of 15 cases of vaccine-induced rabies from Europe and Canada as well as 15 related batches of the ERA-BHK21, SAD Bern, and SAD B19 vaccine strains that were collected between 1991 and 2017 were considered in this analysis.

While viruses from ERA vaccine-derived cases seem to descend from different replication-competent virus haplotypes that were already present within the progenitor ERA-BHK21 vaccine virus strain (Figure 2), the absence of SNV combinations in viruses from SAD Bernorig vaccine-induced rabies cases strongly suggests that only one replication-competent haplotype was present in those cases, despite a higher number of variants. One reason could be the passage history. Both ERA-BHK21 and SAD Bernorig vaccine virus strains descend from the same derivative of the cell culture-adapted field virus strain, namely SAD, but were independently further developed by multiple in vitro passaging on various cell types (see for a more comprehensive overview on this topic). Genetic as well as phylogenetic analyses, however, suggest that the ERA-BHK21 vaccine virus remained closer to its progenitor SAD strain. On the other hand, the SAD Bernorig vaccine virus strain was further developed into the commercial vaccine virus strains SAD Bern and SAD B19. Eventually, increased cell culture adaptation of SAD Bern and SAD B19 vaccine virus strains may have selected for only a single replication-competent haplotype. In parallel, this adaptation towards efficient replication in cell culture is favoring the development of defective interfering particles that provide the background for the high genetic diversity as found by the variant analyses (Figure 1). Unfortunately, our analysis was restricted to viruses for which sufficient sequencing data was available, and thus, could not cover the SAD P5/88 vaccine-induced rabies cases.

As attenuated live vaccines consist of replication-competent virus particles, a major concern in their application is the potential reversion to virulence as demonstrated by vaccines for the severe acute respiratory syndrome coronavirus (SARS-CoV), peste des petits ruminants virus (PPRV), the avian metapneumovirus (aMPV) and Rift Valley fever virus. For the attenuated oral rabies vaccines investigated here, our analyses on alterations in antigenic sites showed no evidence of a reversion to virulence of viruses from vaccine-induced rabies cases either as result of the production process nor of the in vivo replication in affected animals. This is also supported by the population-based analysis conducted, which clearly separated field strains from vaccines and their induced cases (Figure 3). Rather, other factors including the individual immune status or genetic predisposition of the infected animals most likely have contributed to the development of disease. This is supported by the observation, that none of these sporadic cases resulted in onward transmission or was of epidemiological relevance.

The results of this study clearly demonstrate the benefits of utilizing deep sequencing followed by a comprehensive stepwise data analysis focusing on consensus, population, and variant level. While each of the analyses may lead to biased conclusions, only the combination allows for a holistic assessment, and thus should be standard in similar scenarios.