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.