First Detection and Circulation of RHDV2 in New Zealand

Rabbit haemorrhage disease virus 2 (RHDV2) is a highly pathogenic lagovirus that causes lethal disease in rabbits and hares (lagomorphs). Since its first detection in Europe in 2010, RHDV2 has spread worldwide and has been detected in over 35 countries so far. Here, we provide the first detailed report of the detection and subsequent circulation of RHDV2 in New Zealand. RHDV2 was first detected in New Zealand in 2018, with positive samples retrospectively identified in December 2017. Subsequent time-resolved phylogenetic analysis suggested a single introduction into the North Island between March and November 2016. Genetic analysis identified a GI.3P-GI.2 variant supporting a non-Australian origin for the incursion; however, more accurate identification of the source of the incursion remains challenging due to the wide global distribution of the GI.3P-GI.2 variant. Furthermore, our analysis suggests the spread of the virus between the North and South Islands of New Zealand at least twice, dated to mid-2017 and around 2018. Further phylogenetic analysis also revealed a strong phylogeographic pattern. So far, no recombination events with endemic benign New Zealand rabbit caliciviruses have been identified. This study highlights the need for further research and surveillance to monitor the distribution and diversity of lagoviruses in New Zealand and to detect incursions of novel variants.

The RHDV2 taxonomic designation comprises several distinct variants that share the genotype GI.2 structural (S) gene coding sequences but differ in the genotype of their NS gene sequences [28,29].These variants are described by the nomenclature [x]P-[y], where [x]P refers to the genotype of the polymerase (representative of the NS sequence, based on the nomenclature for human noroviruses [1]), and -[y] refers to the S sequence genotype.Reported RHDV2 variants include the GI.1bP-GI.2[30][31][32], GI.4eP-GI.2[33], GI.4cP-GI.2[28], GI.3P-GI.2[6,29,[34][35][36][37][38][39][40][41][42][43], and GII.1P-GI.2[39] variants.While all reported RHDV2 variants can lethally infect hares and young rabbits (unlike RHDV1), they appear to differ in epidemiological fitness [28].However, the precise biological mechanisms underpinning this are yet to be elucidated.European rabbits are an introduced species in both New Zealand and Australia that cause enormous damage to native ecosystems and agricultural industries [45].In New Zealand, rabbits compete with livestock for pasture and provide a stable food source for other invasive mammals and carriers of bovine tuberculosis, such as ferrets, stoats, and feral cats [46,47].Rabbits graze on many vulnerable native plant species and damage erosion-prone soils through burrowing and digging activity.In Australia, rabbit caliciviruses have been used to help manage feral rabbit populations since 1995.Subsequently, in 1997, RHDV1 was illegally imported and released in New Zealand [48].Following the official approval of RHDV1 in New Zealand, this variant was observed to reduce rabbit numbers at release sites, making it an effective measure for rabbit population control there [49].
Prior to the detection of RHDV2 in 2018, three lagovirus variants were known to be present in New Zealand.Two variants serve as biological controls to manage the population of pest rabbits, RHDV1 (V351 Czech, genotype GI.1cP-GI.1c),introduced in 1997 [48], and RHDV1a (K5, genotype GI.1aP-GI.1a),present since 2018 [50].Additionally, a benign European rabbits are an introduced species in both New Zealand and Australia that cause enormous damage to native ecosystems and agricultural industries [45].In New Zealand, rabbits compete with livestock for pasture and provide a stable food source for other invasive mammals and carriers of bovine tuberculosis, such as ferrets, stoats, and feral cats [46,47].Rabbits graze on many vulnerable native plant species and damage erosionprone soils through burrowing and digging activity.In Australia, rabbit caliciviruses have been used to help manage feral rabbit populations since 1995.Subsequently, in 1997, RHDV1 was illegally imported and released in New Zealand [48].Following the official approval of RHDV1 in New Zealand, this variant was observed to reduce rabbit numbers at release sites, making it an effective measure for rabbit population control there [49].
Prior to the detection of RHDV2 in 2018, three lagovirus variants were known to be present in New Zealand.Two variants serve as biological controls to manage the population of pest rabbits, RHDV1 (V351 Czech, genotype GI.1cP-GI.1c),introduced in 1997 [48], and RHDV1a (K5, genotype GI.1aP-GI.1a),present since 2018 [50].Additionally, a benign enterotropic rabbit calicivirus (genotype GI.4P-GI.4)has circulated subclinically in wild rabbit populations since the early 1990s [51].This non-pathogenic variant can induce partial immunological cross-protection to lethal RHD, reducing the effectiveness of rabbit biocontrol programs in regions with high seroprevalence.Then, in April 2018, RHDV2 was detected in New Zealand for the first time [50].In this work, we describe the field and molecular epidemiological characteristics of this incursion.

Tissue Sample Collection
Wild and domestic rabbit liver samples (n = 195) were collected from various locations across the North and South Islands of New Zealand between 2017 and 2021.Samples were collected opportunistically from the carcasses of rabbits found dead by council staff and members of the public.No animal ethics approvals are required for sampling rabbits that are found dead in New Zealand [50].Approximately 100-200 mg of liver tissue was subsampled into 1 mL of RNAlater buffer (ThermoFisher Scientific, Waltham, MA, USA), transported at ambient temperature to Manaaki Whenua-Landcare Research molecular laboratory in Lincoln, and then stored at −80 • C until analysis.

Virus Detection
Total RNA was extracted from 20-50 mg of liver tissue using the iNtRON easy-spin™ [DNA free] Total RNA Extraction Kit (iNtRON Biotechnology, Seongnam, Kyonggi-do, Republic of Korea).cDNA was synthesised from the extracted RNA using SuperScript IV Reverse Transcriptase and Oligo d(T) 20 (Life Technologies Corporation, Carlsbad, CA, USA).Samples were screened for the presence of lagovirus RNA with a multiplex PCR using KAPA2G Robust HotStart PCR Kit (Kapa Biosystems, Inc., Wilmington, MA, USA), as previously described [33] (Table 1).PCR products were visualized on a 2% agarose gel stained with SYBRsafe (ThermoFisher Scientific).

Virus Sequencing
Total RNA (20 µL) from lagovirus-positive liver samples (n = 52) was precipitated with a 0.1 volume of 3M sodium acetate and 3 volumes of 100% ethanol for transportation, and RNA pellets were later resuspended in nuclease-free water prior to PCR amplification.RNAs were reverse-transcribed using Superscript IV Reverse Transcriptase (Thermofisher Scientific) and an Oligo d(T) 18 primer (Sigma-Aldrich, St. Louis, MI, USA), according to the manufacturer's directions.Overlapping primer pairs (Table 2) specific for GI.3P-GI.2were used to amplify the viral genome in five fragments, and libraries were prepared and sequenced using Illumina MiSeq technology (MiSeq reagent kit v2 300-cycles) (Illumina, San Diego, CA, USA), as described previously [52,53].

Phylogenetic and Phylogeographic Analyses
The quality of raw sequencing data was assessed using FastQC (v.0.11.8).Adapter and primer sequences, as well as low-quality regions, were trimmed using Trimmomatic (v.0.38) [55], and overlapping reads were merged using FLASh (v1.2.11) [56].Cleaned sequence reads from each isolate were mapped to the rabbit hemorrhagic disease virus-FRG complete genome (GenBank accession M67473.1)using the Geneious mapper as implemented in Geneious Prime 2021.1.1,and consensus sequences were extracted.
Representative global lagovirus sequences were selected by retrieving all near-complete lagovirus sequences from the NCBI Nucleotide (nt) database (as of 8 September 2021) and identifying clusters with a 95% sequence identity cut-off using CD-HIT-EST [57].Sequences without sufficient temporal data were removed from the dataset.Newly sequenced New Zealand RHDV2 genomes were aligned with these representative global sequences using the FFT-NS-2 algorithm within MAFFT v7.450 [58].Maximum likelihood (ML) phylogenies (n = 179) were estimated separately for the NS and S gene regions (nucleotide positions 28-5304 and 5305-7382, respectively, based on GenBank accession M67473.1)using IQTREE v1.7.0b7 [59], with the best-fit model as determined by ModelFinder [60].Branch support was estimated using 1000 ultrafast bootstrap replicates [61] and 1000 replicates of the SH-aLRT test [62].Phylogenies were rooted at the midpoint between the genogroup I and II clades.We were not able to amplify the 5 ′ -end of three of the New Zealand RHDV2 sequences (OM372659, OM372625, and OM372627); NS sequence phylogenies were estimated both with and without these sequences, and we observed no significant changes in topology; therefore, analyses were conducted with these sequences included.
To reveal any geographic clustering and structure, a maximum likelihood phylogeny was also inferred separately for the 52 New Zealand RHDV2 NS sequences, and the tips were matched to their location on a map of New Zealand.

Bayesian Evolutionary Analysis
All New Zealand RHDV2 sequences and global lagovirus sequences were aligned, and phylogenetic trees were inferred as described above.Additionally, a phylogenetic tree was constructed as described above for a near-full genome alignment.A linear regression of root-to-tip distances against sampling times, as implemented in TempEst [63], was used to determine sufficient temporal signal.The alignment of non-structural genes was found to have the strongest temporal signal (correlation coefficient = 0.53) and was thereafter used for all subsequent analyses.A time-scaled phylogeny was inferred by using a Bayesian Markov chain Monte Carlo (MCMC) approach.To assess the most appropriate clock prior (strict versus uncorrelated log-normally distributed (UCLD)) and tree prior (Gaussian Markov random field Bayesian Skyride model versus constant size coalescent versus exponential coalescent), we used the marginal likelihood estimations (MLEs) (path sampling/stepping-stone sampling) as implemented in BEAUti (v1.10.4) [64].Substitution model GTR + F + I + G4 was used for all MLEs.The best MLE values were reached with a strict clock and exponential coalescent with a chain length of 20 million and were subsequently used for a Bayesian Evolutionary Analysis Sampling Tree (BEAST) analysis (v1.10.4).The analysis was run twice to convergence (ESS > 200) to confirm consistency.TreeAnnotator (v1.10.4), available in the BEAST package, was used to create a maximum clade credibility (mcc) tree from the Bayesian phylogenetic inference results with a 10% burn-in.

RHDV2 Was First Detected in New Zealand in 2017
Routine sampling of dead rabbits was previously conducted in New Zealand to monitor the impacts of the existing biocontrols, RHDV1 (since 1997), and the more recently introduced RHDV1a-K5 (since March 2018).Between 2017 and 2020, 195 rabbit liver samples were tested for the presence of lagoviruses, the majority of those samples in 2018 (n = 146), using a multiplex RT-PCR.Samples were collected from across New Zealand, with most of them originating from the South Island's Otago region (n = 105) (Figure 2).
prior (strict versus uncorrelated log-normally distributed (UCLD)) and tree prior (Gaussian Markov random field Bayesian Skyride model versus constant size coalescent versus exponential coalescent), we used the marginal likelihood estimations (MLEs) (path sampling/stepping-stone sampling) as implemented in BEAUti (v1.10.4) [64].Substitution model GTR+F+I+G4 was used for all MLEs.The best MLE values were reached with a strict clock and exponential coalescent with a chain length of 20 million and were subsequently used for a Bayesian Evolutionary Analysis Sampling Tree (BEAST) analysis (v1.10.4).The analysis was run twice to convergence (ESS > 200) to confirm consistency.TreeAnnotator (v1.10.4), available in the BEAST package, was used to create a maximum clade credibility (mcc) tree from the Bayesian phylogenetic inference results with a 10% burn-in.

RHDV2 Was First Detected in New Zealand in 2017
Routine sampling of dead rabbits was previously conducted in New Zealand to monitor the impacts of the existing biocontrols, RHDV1 (since 1997), and the more recently introduced RHDV1a-K5 (since March 2018).Between 2017 and 2020, 195 rabbit liver samples were tested for the presence of lagoviruses, the majority of those samples in 2018 (n = 146), using a multiplex RT-PCR.Samples were collected from across New Zealand, with most of them originating from the South Island's Otago region (n = 105) (Figure 2).RHDV2 was first detected in New Zealand in a liver sample from a wild adult European rabbit found dead on 6 April 2018 in the Awatere Valley, Marlborough, in the northern point of the South Island (Figure 2) [78].However, following this detection, December 2017 (Figures 2 and 3), preceding the original April 2018 sample collection date.RHDV2 was subsequently detected in 71 samples in 2018 (March, April, June, and August to November), in 4 samples in 2019 (April, June, and December), and in 9 samples in 2020 (January, April, October, and December) across New Zealand's North and South Islands (Figures 2 and 3).Mixed testing results suggesting co-infection with multiple lagoviruses were observed in three samples in 2018 and one sample in 2019 (Figure 2B), which were positive for both the RHDV2 and RHDV1-V351 Czech variants.
retrospective testing of previously collected samples identified RHDV2 in two wild rabbit samples collected from a single farm from the Bay of Plenty region in the North Island in December 2017 (Figures 2 and 3), preceding the original April 2018 sample collection date.RHDV2 was subsequently detected in 71 samples in 2018 (March, April, June, and August to November), in 4 samples in 2019 (April, June, and December), and in 9 samples in 2020 (January, April, October, and December) across New Zealand's North and South Islands (Figures 2 and 3).Mixed testing results suggesting co-infection with multiple lagoviruses were observed in three samples in 2018 and one sample in 2019 (Figure 2B), which were positive for both the RHDV2 and RHDV1-V351 Czech variants.Amongst the New Zealand RHDV2 sequences, there was a clear pattern of increasing genetic diversity over time (Figure 5), further supporting a single introduction with ongoing transmission.Several 2020 sequences have comparatively long branch lengths, most likely suggesting substantial undetected spread in New Zealand rabbit populations.A clear phylogeographic signal was also observed despite the heavy sampling bias towards the South Island.Most North Island sequences clustered close to the root of the New Zealand phylogeny, correlating with the earliest detections in the North Island.Sequences from in infected (leading to shedding large amounts of virus into the environment), indicate that fomite transmission, mechanical transmission via insect vectors [80], and/or the movement of infected carcasses via scavengers [31,40,[81][82][83] are likely key contributors to long-distance spread.
While RHDV1 has been endemic in New Zealand since 1997 [48], New Zealand remained free of RHDV2 until its first detection in 2018, with subsequent retrospective detection to 2017 as described in this study.Here, we show that a single incursion event, probably into the North Island in 2016 or 2017, led to the establishment of RHDV2 in New Zealand's wild and domestic rabbit populations.While the earliest detection was from December 2017, there was limited sampling of wild rabbits in New Zealand for molecular testing prior to 2017.Phylodynamic analysis supports a recent incursion, with a most likely incursion date between March and November 2016.It is likely that RHDV2 underwent cryptic spread for approximately one year before the initial detection.This is not overly surprising since the sampling of wild rabbit carcasses is opportunistic and biased.Carcasses can be difficult to recover when infected rabbits die in inaccessible locations, carcasses may be taken by scavengers, and someone must be willing to collect and submit a sample for analysis.
Since the incursion of RHDV2 into New Zealand, phylogenetic analysis suggests ongoing spread throughout New Zealand rabbit populations and at least two spread events from the North to the South Island.The North and South Islands are separated by the Cook Strait, a distance of 23 km at its narrowest point.Transmission across the Strait may have been human-mediated, via fomites, through the self-directed or passive/windassisted movement of mechanical vectors such as flies, or perhaps via raptors or scavenging birds.Other authors have suggested similar routes for the incursion of lagoviruses onto and/or between islands [31].The limited scope of the available epidemiological data precludes a more comprehensive analysis of the spread mechanisms of RHDV2 in New Zealand.Importantly, in response to the incursion of RHDV2 in New Zealand, Agricultural Compound and Veterinary Medicines approved the importation of the bivalent Filavac VHD K C+K vaccine for the immunisation of domestic rabbits against RHDV1 and RHDV2.This enables pet owners, breeders, and meat producers to protect domestic rabbits from both RHDV1 and RHDV2 virus variants.Previously, only the monovalent Cylap RHDV1 vaccine was available in New Zealand.
Elsewhere in the world, RHDV2 has rapidly replaced previously dominant RHDV1 variants [39,84,86].Likewise, in New Zealand, RHDV2 was the only variant detected in 2020, with the pathogenic RHDV1 Czech and K5 variants only detected until the end of 2019 (Figure 2B).However, it must be noted that the number of samples analysed in 2020 was relatively small compared with previous years.A large number of samples were submitted in 2018 in association with the coordinated national release of the biocontrol variant RHDV1-K5 and a public awareness campaign.However, sample submission declined following the release.
In contrast with Australian observations, no recombination events between RHDV2 and circulating benign rabbit caliciviruses have been identified in New Zealand so far.While multiple recombination events were detected in Australia in the 2-year timeframe

Figure 2 .
Figure 2. Stacked bar charts of tested samples from 2017 to 2020.Samples are coloured by (A) the location of collected samples and (B) testing results by RT-PCR per year.

Figure 2 .
Figure 2. Stacked bar charts of tested samples from 2017 to 2020.Samples are coloured by (A) the location of collected samples and (B) testing results by RT-PCR per year.RHDV2 was first detected in New Zealand in a liver sample from a wild adult European rabbit found dead on 6 April 2018 in the Awatere Valley, Marlborough, in the northern point of the South Island (Figure2)[78].However, following this detection, retrospective testing of previously collected samples identified RHDV2 in two wild rabbit samples collected from a single farm from the Bay of Plenty region in the North Island in

Figure 4 .
Figure 4. Maximum likelihood phylogenies of New Zealand and representative global RHDV2 sequences.All near-complete lagovirus genomes were downloaded from GenBank, and representative sequences spanning the known diversity were derived using CD-HIT-EST.These were aligned with the New Zealand RHDV2 sequences.Phylogenies were inferred for the (A) non-structural and (B) structural genes independently.Taxa are coloured by country of origin, and the genogroup, genotype, and variant type are indicated by clade labels.Strong (100) and moderate (80-99) clade support is indicated by the filled dark grey and light grey circles, respectively, at internal nodes.The scale bar represents substitutions per site.The clade for New Zealand samples (this study) has been collapsed and is shown as a big, red-filled dot.The collapsed New Zealand clade for the non-structural (A) genes is shown in Figure 5.

Figure 4 .
Figure 4. Maximum likelihood phylogenies of New Zealand and representative global RHDV2 sequences.All near-complete lagovirus genomes were downloaded from GenBank, and representative sequences spanning the known diversity were derived using CD-HIT-EST.These were aligned with the New Zealand RHDV2 sequences.Phylogenies were inferred for the (A) non-structural and (B) structural genes independently.Taxa are coloured by country of origin, and the genogroup, genotype, and variant type are indicated by clade labels.Strong (100) and moderate (80-99) clade support is indicated by the filled dark grey and light grey circles, respectively, at internal nodes.The scale bar represents substitutions per site.The clade for New Zealand samples (this study) has been collapsed and is shown as a big, red-filled dot.The collapsed New Zealand clade for the non-structural (A) genes is shown in Figure 5.

Figure 5 .
Figure 5. Maximum likelihood phylogeny and phylogeographic analysis of New Zealand RHDV2 non-structural genome sequences.A phylogeny was inferred for the 52 non-structural New Zealand RHDV2 sequences.Taxa names are coloured by region of sample collection, and branches are highlighted by collection year.Tip points represent whether samples were collected in the North or the South Island.Strong and moderate clade support is indicated by the filled dark grey and light grey circles, respectively, at internal nodes.The scale bar represents substitutions per site.Lines connect the sample location to the location on the New Zealand map.

Figure 5 .
Figure 5. Maximum likelihood phylogeny and phylogeographic analysis of New Zealand RHDV2 non-structural genome sequences.A phylogeny was inferred for the 52 non-structural New Zealand RHDV2 sequences.Taxa names are by region of sample collection, and branches are highlighted by collection year.Tip points represent whether samples were collected in the North or the South Island.Strong and moderate clade support is indicated by the filled dark grey and light grey circles, respectively, at internal nodes.The scale bar represents substitutions per site.Lines connect the sample location to the location on the New Zealand map.

Table 2 .
Primers used in this study.