Tracking the Origin of a Rabbit Haemorrhagic Virus 2 Outbreak in a Wild Rabbit Breeding Centre in Portugal; Epidemiological and Genetic Investigation

It is known that RNA viruses’ adaptability is the result of their extreme mutation frequency, as many genomic variants are created at a high generation rate [2]. Among them, non-synonymous mutations at one or more sites may alter the RNA virus phenotypes and virulence [3,4]. Indeed in RHDV infections, like for other viruses, antigenic drift (minor changes in surface generated through point mutations) render hosts susceptible if at least a portion of the population has not been previously exposed to the new antigen(s) [5].

It is known that RNA viruses' adaptability is the result of their extreme mutation frequency, as many genomic variants are created at a high generation rate [2]. Among them, non-synonymous mutations at one or more sites may alter the RNA virus phenotypes and virulence [3,4]. Indeed in RHDV infections, like for other viruses, antigenic drift (minor changes in surface generated through point mutations) render hosts susceptible if at least a portion of the population has not been previously exposed to the new antigen(s) [5].
The mean evolution rate for RHDV based on the analysis of complete capsid gene sequences was estimated in 5.48 to 7.7 × 10 -4 substitutions/site/year [6,7]. It has been suggested that the emergence of virulence occurred in the first quarter of the 20 th century in a different host until it jumped to Oryctolagus spp. in China, where the disease was first reported in 1984 [7][8][9]. A similar process of lagovirus species jump was proposed for the emergence of RHDV2, although its putative host was not yet identified [8].
Although genetically related with RHDV, RHDV2 emergence in 2012 in Portugal had a tremendous impact in wild rabbit populations and was a major drawback for the re-introduction of the captive-bred Iberian lynx (Lynx pardinus), the most endangered wild felid in the world.
Wild rabbits are the main component of Iberian lynx diet and a keystone species on the equilibrium of the Iberian ecosystem [10]. Endemic to the Iberian Peninsula [11], the Iberian lynx population size dramatically decreased during the 20 th century [12]. In 2002 was declared 'critically endangered' (CR) by the International Union for the Conservation of Nature (IUCN) [13,14]. Major efforts (Action Plan for the Conservation of the Iberian Lynx in Portugal; Iberian Lynx exsitu Conservation Programme; projects under the European LIFE Programme [Lince Moura/Barrancos -LIFE06 NAT/P/000191; Habitat Lince Abutre -LIFE08 NAT/P/000227; Iberlince -LIFE10 NAT/ES/000570], among others) have been made in the last decade to reverse the decline of this species. As a result, in 2015, the Iberian lynx species status was shifted to "Endangered" (EN), due to the increase of the populations in Spain. However, the RHDV2 emergence in Portugal in 2012 almost decimated the entire wild rabbit population in particular niches of mainland [15,16] and Azores [17]. A comparable scenario was also observed in Spain [18]. Although both Oryctolagus cuniculus subspecies are equally affected by RHDV2 [15,18], the susceptibility of the subspecies O. cuniculus algirus acquires special importance considering its restricted distribution to southwest [15] and the dramatic implications of RHDV2 outbreaks on its frail conservation status [16,18]. O. cuniculus algirus is a key prey species for several carnivores, including the Iberian Lynx and the Iberian imperial eagle (Aquila adalberti), both emblematic and endangered species in Portugal, and a population downsize may lead to a series of major ecological and economic problems [19].
To counteract the effects of RHDV2 on the decrease of wild rabbit populations in Portugal, several measures were implemented among which restrictive hunting regulations and captive breeding of wild rabbits for subsequent return to wildlife. Rabbits restocking efforts have been recently reinforced to repopulate depleted areas in the sequence of RHDV2 outbreaks as a measure to fasten self-sustainable population establishment. Wild rabbits are breed in captivity units and set free in areas assigned either for cinegetic activity or for the Iberian lynx reintroduction. For the latter purpose, most units are located in Moura and Barrancos municipalities.
Since its first detection in Portugal in 2012 [15], RHD has been observed in some of those wild rabbit breeding units. However, to date, none of these past outbreaks were characterized and publically reported.
Here, we describe an RHDV2 outbreak in a wild rabbit captive breeding unit, built in 2004 for the wild rabbit repopulation of the Noudar Nature Park, one of the possible sites for future conservation translocations of the Iberian lynx. Phylogenetic analyses were carried out to assess the variability and putative relation of two RHDV2 Barrancos-2016 strains with other strains obtained previously in the unit, the first time in late December 2012, as well as with strains from the South, envisaging to unravel possible links between different outbreaks.

Samples and inquiry
A total of 13 victimized rabbits collected between March and April 2016, were investigated at Instituto Nacional de Investigacao Agraria e Veterinaria (INIAV). An inquiry was carried out to gather epidemiologic information on the wild rabbit captive breeding unit regarding its exact location and surroundings, extension of the parks, road accesses and animal imports (dates and animal number).

Estimation of density and mortality in 2016
Direct counts are difficult since these animals spend most of the day hidden in their warrens, making indirect counting methods preferable. Rabbits' abundance was therefore estimated based on the intake of commercial dry rabbit food, supplied weekly to each park, assuming it was the main food source and that the average consume per rabbit per week is about 1 kg.
Mortality rates were estimated comparing the animal density inferences, calculated as described above, before and after disease outbreaks. In the unit, enumeration methods [20] (collection of cadavers) are used quite rarely for mortality estimations since they can contribute to bias results and underestimate mortality numbers as a consequence of pre-emergence death (before litters' emerge from the nest) and predation [21].

Anatomopathological examination
Necropsies were carried out at the University of Evora. Liver and lung samples were collected and sent to INIAV for histopathological examination. Samples were fixed in 10% buffered formalin and embedded in paraffin by standard procedures. Five micrometer-thick sections were stained with haematoxylin and eosin (H & E) prior light microscopy examination.

Virological examination
Liver and lungs samples were homogenized with phosphate buffered saline (PBS) and clarified at 3,000 g for 5 min. DNA and RNA were extracted from 200 µl of the clarified supernatant in a BioSprint 96 nucleic acid extractor (Qiagen, Hilden, Germany) according to the manufacturer's instructions. RHDV2-RNA was assessed by RT-qPCR [22] using the One Step RT-PCR kit (Qiagen, Hilden, Germany). Screening for classical RHDV strains was performed by conventional RT-PCR followed by sequencing analysis of the amplicons obtained with primers RC-9 and RC-10 [23], also using the One Step RT-PCR kit (Qiagen, Hilden, Germany). Myxoma virus was investigated by qPCR as described by Duarte et al. [24] with the FastStart TaqMan Probe Master Kit (Roche, Roche Diagnostics GmbH, Manheim, Germany). Cq values above 40 were considered negative.
Sequencing was carried out using a BigDye TM Terminator cycle sequencing kit (Applied Biosystems, Foster City, CA, USA) and the nucleotide sequences determined on an automated 3130 Genetic Analyzer system (Applied Biosystems, Foster City, CA, USA). The complete vp60 sequences obtained from two animals (strains Barrancos1/PT16 and Barrancos2/PT16) were submitted to GenBank and attributed the accession numbers KX132812 and KX132813. Nucleotide alignments were performed with Clustal Omega [25].
The percent nucleotide variability among strains based on the vp60 gene (Table 2) and the genetic divergence between populations, inferred by the nucleotide differences and substitutions between viral populations, was calculated using the DnaSP software (Version 5.10.01) [26], in order to further explore the genetic relationships among strains. Strains from the same geographic origin and year of collection were included in the same population.

Phylogenetic analysis
The appropriated substitution model was determined resourcing to R software (R Development Core Team, 2009). The GTR model [27] with gamma-distributed rate variation across sites showed the lowest BIC and AICc values and was subsequently used to infer phylogenetic relationships using ML analysis. Robustness of the tree nodes was assessed by bootstrapping 1000 times.
For the Bayesian analysis the CLUSTAL Omega results were converted to the NEXUS format using Mesquite software [28]. The phylogenetic tree was obtained with a Bayesian inference of phylogeny throughout the MrBayes version 3.1.2 software that uses the Markov chain Monte Carlo simulation technique to approximate the posterior probabilities (PP) of trees [29,30]. MrBayes analysis was performed using the GTR model (nst=6) with gamma-shaped rate variation with a proportion of invariable sites (rates=invgamma). The analysis was run for 10 6 generations (ngen=10 6 ) with four chains of temperature (nchains=4), and each chain was sampled every 10 th generations (samplefreq=10).
The graphical representation and edition of the phylogenetic trees were performed with FigTree v1.3.1.
The Network software (version 5.0) was used to reconstruct a network using haplotypes (vp60 gene) by median joining [31].

Results
Data collected from the inquiry revealed that the unit affected by RHDV2, located in Barrancos municipality near the Spanish border, encompasses one quarantine and four reproduction parks. The five parks are sited in holm oak woodlands, only accessible by sand-clay roads, and lay on a 999 hectare extension area delimited by two riversides: Southeast and within the national territory by the Múrtega River, and Northwest by the Ardila River, in the Spanish border ( Figure 1). Parks 1 (P1) and 3 (P3) are located closer to a busy road and therefore more exposed to human activity. In contrast, parks 2 (P2) and 4 (P4) lay on a more protected area, given their proximity to the riverside and, in the P4 case, also to a closed ore mine more distant from the roads.
The first animals (100 adults, 80 females and 20 males from subspecies O. cuniculus algirus) were introduced in 2004 (by the time of P1 construction) and originated from captive-bred wild populations from Montemor-o-Novo, located 136 Km northwest. Introduction of additional O. cuniculus algirus rabbits from other breeding parks was performed regularly (Table 1). From the second half of 2013 afterwards, following the first RHDV2 outbreak a higher number of foreign bred rabbits, mostly originated from the Spanish Granada Province, were introduced into the parks ( Table 1). The last introduction took place in 2015. The actual animal density is estimated in approximately 120 animals per hectare (≈1500 rabbits over a 12.5 hectare area). Rabbits are vaccinated against myxoma virus and classical RHDV and identified with a tattoo or ear tag before release. No other prophylactic measures besides a coccidiostatic administration through the feed are carried out.
Mortality by RHDV2 had been observed in the parks since 2012 (Table 1).
In 2016, the first RHDV2 fatalities occurred in March, affecting rabbits of all ages from P1, P2 and P3. No casualties were registered in P4. Until April 2016, the mortality rate was estimated at approximately 8.67% (130/1500). No more casualties were registered from April onwards.
A total of 13 animals, corresponding to approximately 10% of the estimated casualties, were necropsied and epistaxis was observed in all animals. Haemorrhagic tracheitis, hepatic congestion and friable liver were also observed. In two animals, histopathology revealed congestion of the lungs and necrotic lesions in the liver.
All liver samples (n=13) were positive to RHDV2 by RT-qPCR [22] and negative to myxoma virus by qPCR [24]. Sequencing analysis of vp60 386 bp long fragments ruled out the concomitant presence of classical RHDV strains.
The complete vp60 nucleotide sequences (1740 bp long) obtained from two animals (KX132812 and KX132813) showed to be identical (100% similarity  The 2016 strains were compared with four strains from Barrancos obtained in previous years from dead wild rabbits collected in the same unit, two in January 2013 (KF442963 and KF442964) and other two in February 2014 (KM115675 and KM115676), as well as with two strains from Mértola (KM115712 and KM115713), a municipality located 81.22 km southwest from the unit, obtained in January 2014. The percent of nucleotide similarity among strains based on the vp60 gene (Table 2) and the nucleotide differences and substitutions between viral populations showed concordant values (Table 3). Besides the six SNPs already mentioned above, the Barrancos-2016 strains differed from the Barrancos-2013 and from the Barrancos-2014 strains in 26 and 37 additional residues, accounting for a total of 32 and 43 mutations, respectively (Table 3). Interestingly, the Barrancos-2016 strains showed a slightly higher nucleotide similarity and lower nucleotide differences with the Barrancos-2013 strains than with Barrancos-2014 (Tables 2  and 3). The nucleotide variability found between the Barrancos-2013 and Barrancos-2014 strains (Table 3), involved 27 nucleotide variations. Regarding strains obtained in the same year in the Barrancos municipality, a similarity of 99.94% and 100% was observed. Among the oldest (2013 and 2014) and the 2016 strains, respectively ( Table  2). No elations could be taken regarding the cumulative nucleotide variation of RHDV2 since the collection of the biological materials in each year occurred always within a short period of one month.
When comparing strains from the two municipalities, the Barrancos-2016 sequences (KX132812 and KX132813) showed 97.93% of similarity with the Mértola-2014 strains (Table 2), from which they differed in 31 nucleotides (Table 3)

Discussion and Conclusions
Immediately before release, the wild rabbits raised in the unit    are vaccinated against myxoma virus and classical RHDV. However, RHDV2 vaccination, a growing practice in the domestic rabbit industry, is not practiced in this unit. Vaccination of captive wild rabbits against RHDV2 may provide temporary protection against infection and acquired immunity may increase survival of vaccinated rabbits if contact with the field strains occurs seven to 10 days after vaccination up to one year [32]. Vaccination is crucial for the establishment of herd immunity, although it can also be used as an effective post-exposure tool, in particular situations, for disease control according to the OIE [32]. The difficult feasibility of vaccine administration prior return to the wild, and the limited immunity induced by a single boost raises doubts about its contribution to a faster reestablishment of natural wild rabbit populations. Since RHDV2 vaccines are not freely available on the market, their use requires special licenses from the Veterinarian Authorities. Vaccination of wild captive rabbits against RHDV2 was exceptionally authorised by Portuguese National Authority for Animal Health (DGAV) in a few other wild rabbit breeding units.
For the revival of wild populations, the development of natural immunity against RHDV2 is crucial. Many factors will affect infection outcome at the population level, such as strain virulence [9,33], herd immunity due to previous infections [4,34] population density, season of the year [7,35], anthropogenic changes in the environment [4], landscape [9] and viral transmission between populations [4,9,36].
The geographic features of the unit, circumscribed by rivers, may offer a natural barrier to animal movement which is expected to slow down virus dissemination within the parks, similarly to what was described for RHDV [7,37], giving the virus more time to evolve [9]. This might explain the genetic diversity found between the strains from this unit obtained across this 4-years period. Also, the rapid transmission of a pathogen will affect the local spread of the disease reducing the number of susceptible individuals [4]. Hence, RHDV2 moderate virulence compared to RHDV, constitutes a selective advantage of this virus [9], as it was demonstrated that moderately virulent RHDV-related strains can invade wild rabbit populations with greater efficiency than highly virulent or non-pathogenic strains [9,33].
Interestingly, despite all the other parks in the unit were affected by RHDV2, there were no casualties recorded in P4. Since the virus first detection in the unit in 2012, animals from P4 seem to be passing unharmed. Although the factors behind the lower incidence of infection in P4 are not fully clarified, it is likely that its geographic isolation may have protected rabbits more efficiently from infection. However, the higher levels of humidity near the riverside probably also favour mosquitoes, which are often implicated in disease transmission as mechanic vectors for RHDV [38], therefore puzzling this uneven occurrence of disease.
Also, an artificial boundary by electric fencing delimits each park adding up to the protection from predators given by several natural shelters located within the parks. Yet, the involving forestland provides suitable habitats for many free-ranging species including wild rabbits and small predators and scavengers such as badgers (Meles meles), mongooses (Herpestes ichcneumon), weasels (Mustela nivalis), foxes (Vulpes vulpes), stone martens (Martes foina), as well as birds of prey, often seen in the area, whose proximity with the parks may favour the direct or indirect contact with captive wild rabbits. Some predators and scavengers species can remove small animal carcasses from the site of death to other locations [39] and, concurring with that, a few specimens collected during this outbreak were damaged, probably due to scavenging. In fact, consumption of RHDV-victimized rabbits may be at the origin of the RHDV antibodies detected in red foxes and other predators and scavengers living in sympatry with affected rabbit populations [40][41][42]. These serological evidences pointed to the possibility of other species, apart from rabbits, being involved in the epidemiology and persistence of the disease [43].
Given the complex eco-epidemiology of RHDV2 depicted above, while efficient herd immunity is not achieved, serious difficulties may be posed to the establishment of self-sustainable wild rabbits populations in Baixo-Alentejo.
The phylogenetic analysis did not allow clarifying the origin of the infection since the Barrancos-2016 strains did not group with any other strain presently known, showing however to be more closely related with strains from Azores (supported by a PP of 0.85) (Figure 2). When looking to viruses originated in the same district, the Bayesian tree showed that Barrancos-2016 strains are more closely related with the Mértola-2014 strains than with other strains from Barrancos (pp of 0.99). The Barrancos-2013 strains clustered with strains from the South and Centre of Portugal mainland, while the two Barrancos-2014 strains grouped independently in a more ancestral branch, contradicting the temporal structure with regards to the position in the tree. However, these results were further clarified by the phylogenetic network analysis performed ( Figure 3) where missing intermediates from which the Barrancos 2013, 2014 and 2016 strains may have evolved, present a chronological relative correct position. Based on the phylogenetic analyses, the Barrancos-2016 could represent an independent introduction of RHDV2 from unknown origin. However, it is worthy of mention that the unit is located in a preserved area were hunting activity is not allowed, therefore movement of hunters, which could be implicated in the virus dissemination from other locations, did not account for the occurrence of the outbreak.
Despite disease reports in the last five years suggest that RHDV2 may have circulated endemically in the region with epizootic occurrences, mortality was not observed in the immediate surrounding areas of the parks. Nonetheless, the possible involvement of the local wild rabbit population in the outbreak couldn't be proven or excluded, since the outside area was not systematically investigated prior the outbreak, a limitation of this study. The new virulent Barrancos-2016 strains could have evolved from previous circulating strains through cumulative mutations. The phylogenetic similarity with other strains isolated from Portugal, even though with no apparent geographical or epidemiological link, is more in favour of a common infection source rather than a less probable convergent evolution.
Finally, regarding rabbit translocations, the last introduction of animals in the park took place in January 2015, 14 months before the 2016 fatalities. It is unlikely that the virus was introduced from the Spain park supplier, where no records of RHDV2 infection were reported. Furthermore, during the at least 15 day quarantine carried out before translocation, no clinical signs of disease were observed.
The mortality rate estimations can contain some bias, a limitation inherent to the used methodology and related to the fact that in P2 and P3, in contrast with P1, rabbits have access to some seasonal natural food source, like acorn and grass, being also more vulnerable to predation, particularly by birds of prey. Nevertheless, the mortality rate estimated in 2016 was comparably lower than that observed in former years, namely in the outbreak between December 2012 and January 2013 (≈200 victimized animals) when the virus was first detected in the unit. At the time, the centre encompassed two parks with 200 and 300 animals, were a mortality rate of 42% and 39% (average 40.5%) was registered. This lower mortality rate compared with other reported outbreaks may be the reflex of a gradual RHDV2 immunity development in wild rabbit populations and its passive transmission to the offspring, a promising sign towards the establishment of host-virus equilibrium.