Introduction

Applications of coalescent-based analytical techniques to molecular data allow the reconstruction of demographic changes over time of the studied organisms (Drummond et al., 2005). These techniques have been successfully applied to reconstruct the demographic history of pathogens as well as higher order organisms (Raghwani et al., 2019, Pacioni et al., 2015). In parallel, recent efforts have aimed to extend phylogenetic birth-death models to estimate parameters relevant for epidemiological studies. In these models, the rate of transmission and recovery can be estimated allowing calculation of the basic reproductive number (the number of secondary infections per each infectious case in a naïve population, R 0) (Stadler et al., 2012). Further developments of these models aimed to relax limiting assumptions (e.g. that the sampled infected individuals either die or are removed) while co-inferring evolutionary (i.e. phylogeny) and epidemiological processes (e.g. Gavryushkina et al., 2014, Stadler et al., 2013). Investigations on these research topics constitute the currently very active field of research termed “phylodynamics”.
Viruses are ideal targets for these models as they have very high mutation rates, allowing detection of changes in population dynamics over a very short period (Drummond et al., 2003). The utility of these models for recovering epidemiological processes has been investigated with simulation studies (Stadler et al., 2013) as well as observational comparisons, mostly with the current understanding of diseases affecting humans (e.g. Stadler et al., 2014). Their applications to wildlife diseases have been limited and their accuracy has been broadly assumeda priori . However, several aspects of wildlife disease investigations present additional challenges compared to human infectious disease research. Knowledge of the epidemiology of the pathogen under study is often limited (Wobeser, 2007), making the selection of biologically appropriate models difficult. Compared to human diseases, the sampling regimes in wildlife are generally limited both temporally and spatially. Also, the lack of clear definition of host’s epidemiological compartments (groups of individuals with similar characteristics in respect to the pathogen: e.g. susceptible, infected, etc.) due to knowledge gaps, or inadequate sampling presents further challenges. Models that do not consider this additional complexity may fail to accurately represent dynamics of infectious diseases. Furthermore, it is difficult to predict the possible biases of parameter estimates when analysing data with these simplified models. Hence, there is a need for further investigation of phylodynamic models to determine their utility for analysing the epidemiology of wildlife diseases. The analysis of Rabbit Haemorrhagic Disease Virus (RHDV) data from Australia provides a unique opportunity to validate the application of these methods to monitor virus dynamics within a wildlife host. A strong molecular clock-likeness (i.e. the fact that the accumulation of mutations is strictly related with time) of RHDV sequence data has been demonstrated by a number of previous studies (Eden et al., 2015b, Kovaliski et al., 2014, Eden et al., 2015a). This property makes RHDV particularly suitable for the analysis of the heterochronous sequence data that were available for this study. Subsequent to the escape of the virus from quarantine field trials in the mid-1990s, monitoring sites were set up throughout the country to detect the arrival of the virus in different regions. This monitoring allowed a detailed reconstruction of the temporal-spatial spread of the epidemic (e.g. Bruce et al., 2004, Kovaliski, 1998). We, consequently, have a strong prior knowledge of what dynamics we should be able to retrieve from our phylodynamic models.
Since their introduction in 1859, European rabbits (Oryctolagus cuniculus ) have had a devastating impact on agricultural production and biodiversity in Australia, with competition and land degradation by rabbits listed as key threatening processes under theCommonwealth’s Environmental Protection and Biodiversity Act (1999) . Prior to the occurrence of Rabbit Haemorrhagic Disease in Australia, rabbit damage was reported to cost Australian agriculture up to $600 million annually. More recent estimates place the cost of rabbits on agricultural production losses at approximately $200 million annually, with an additional $25 million spent each year on management and research costs (Cox et al., 2013). In addition to physical control techniques (e.g. shooting, poisoning, warren ripping and harbour destruction), the myxoma virus (Myxomatosis cuniculi ) and RHDV constitute the most important landscape scale control strategies for rabbits in Australia, having resulted in excess of $70 billion in economic savings since the 1950s (Cooke et al., 2013). The RHDV (Czech strain), commonly identified as RHDV1, escaped from Wardang Island, South Australia, in October 1995, where its efficacy as a biological control agent was being tested under natural conditions. Soon after the virus escaped to the mainland, monitoring programs that looked at both rabbit abundance and individuals’ serological status, were established in multiple sites across all States and Territories. Through this monitoring, it was demonstrated that the virus became established nationwide within 12 months (Kovaliski, 1998) and it was hypothesised that the rapid spread of the virus was partly due to flies feeding on infected carcasses (Henning et al., 2005, McColl et al., 2002). In experimental trials, death of susceptible rabbits normally occurs within two to three days post infection (Elsworth et al., 2014) and at some field monitoring sites, population reductions of more than 90% were recorded (Mutze et al., 1998, Bruce and Twigg, 2005). However, based on subsequent monitoring, it also became apparent that the efficacy of RHDV was highly variable with average abundance reductions of 28% in high rainfall areas compared to 67% in lower rainfall areas throughout Australia (Bruce et al., 2004, Cooke et al., 2002, Cox et al., 2013).
At some locations, serological analyses detected antibodies cross reacting to RHDV in rabbits prior to the arrival of RHDV, suggesting that RHDV effectiveness was likely impeded by the presence of a related, non-pathogenic calicivirus (Cooke et al., 2002, Robinson et al., 2002a, Nagesha et al., 2000), which was later identified (Strive et al., 2009, Strive et al., 2010, Strive et al., 2013).
Since approximately 2005, rabbit abundance has recovered by 5-10 fold compared to initial post-RHDV1 levels, likely due to the decreasing efficacy of RHDV as control agent (Mutze et al., 2014a). However, the exact mechanisms involved are not well understood. While Elsworth et al. (2014) demonstrated that field RHDV variants collected between 2007 and 2009 were associated with increased case fatality rates in Australian wild rabbits from the same location, it is generally accepted that other mechanisms may have facilitated the recovery of rabbits. Young rabbits have a physiological resistance to lethal infection until at least 4 weeks of age (Fordham et al., 2014, Abrantes et al., 2012, Robinson et al., 2002b). Mechanisms that enable rabbit populations to recover may therefore include a shift in the age classes that are infected (Mutze et al., 2014b, Elsworth et al., 2014), which results in a lower fatality rate, variability of the doses of infection (Wells et al., 2018), or possibly a combination of both. Further elements that may influence RHDV1 epidemiology are the host population structure, existing (age-related) epidemiological compartments within each rabbit sub-population and the presence/absence of localised myxoma disease epidemics (Barnett et al., 2018).
We aim to analyse available RHDV molecular data to reconstruct the virus dynamics from the first RHDV epidemic in Australia to 2014. As we have a relatively good understanding of the disease dynamics within the first 10 years (1995-2005) since the initial virus release in Australia, we can validate the applied phylodynamic models with field data. A secondary aim is to retrieve additional epidemiological parameters such as effective reproductive number (the average number of secondary infections per each infectious case, R e) with a particular focus on the detection of recent changes (post 2005) in naturally occurring strains.
Based on general epidemiological theory, we expect the virus to evolve to maximise R e, with a trade-off in virulence. However, field data suggest that RHDV has not evolved to become more attenuated but appears to be co-evolving with rabbits to maintain very high levels of virulence (Elsworth et al., 2014) - probably because flies that have acquired RHDV from rabbit carcasses are one of the main means for disease transmission over long distances , as opposed to live, diseased animals, which keeps the selection pressure for high levels of virulence that produces dead bodies for effective transmission (Schwensow et al., 2014, Elsworth et al., 2014). Concurrently, as the virus established in Australia, various levels of acquired immunity in animals surviving RHDV infection would have led to a reduction in the pool of susceptible individuals. Our initial hypothesis, based on the available evidence, is that we would expect R e to be much higher than 1 during the initial outbreak in Australia because of the naïve conditions of Australian rabbit populations (with the exception of a small proportion of rabbits with cross-immunity due to the non-pathogenic calicivirus). As a result of a substantial reduction of the host density, with the number of susceptible hosts further reduced by acquired immunity in animals surviving virus infection, we would expect R e to decline, although staying above 1 because RHDV1 successfully established itself as endemic in Australia. It is harder to anticipate what the results of the analysis may be in regard to RHDV1 dynamics between 2005 and 2014, when field data indicated a drastic recovery of most rabbit populations, with one possibility being that R e continued to slowly decrease below one or stayed stable.