Discussion
We evaluated the performances of recently developed phylodynamic and
coalescent-based analyses using sequence data for RHDV and compared
these results with current understanding of RHDV field epidemiology and
our current understanding of the disease in wild rabbits in Australia.
Phylodynamic analyses
The two phylodynamic analyses (SABDSKY and BDSKY) were very similar,
probably because there was only one sample ancestor in the analysis. In
relative terms, these analyses agreed with our expectations as the
median R e estimates were larger than one for most
of the analysed time windows. However, it was surprising that the
possible range of R e values (i.e. the HPD) in the
first five years since RHDV release included estimates only marginally
larger than one, with a median actually lower than one in some analysis,
in spite of the known rapid spread of the disease in Australia and the
approximate 90% rabbit abundance reduction observed at some sites (Cox
et al., 2013). Furthermore, R e estimates were not
substantially higher than the successive epochs (and possibly not as
high as we would have expected). It is conceivable thatR e rapidly reduced as soon as the first outbreak
wave was over, while in our analysis we are attempting to estimate aR e value that encompasses the first five years
since the release of the virus in Australia, which may force the
analysis to a lower value. We also argue that the limited coverage of
the data that we had access to, possibly contributed to this result. The
data from GenBank was very sparse because we only obtained data from
approximately one individual per outbreak, for a very small fraction of
the total number of outbreaks that have occurred in Australia and, most
likely, our analysis lacks the capacity to estimateR e in a classical epidemiological sense (i.e. the
number of infections created by an infected individual). In other words,
rather than operating at an individual basis, we argue that, given the
resolution of our temporal and spatial sampling, our inference is on
parameters more related to the rate of spread and recovery of local
outbreaks, rather than individuals. In this view,R e could be considered the number of virus field
variants that are generated by a virus lineage; δ represents the rate at
which lineages are not detectable anymore; and the sampling proportion
refers to the proportion of outbreaks detected. Since the infectious
period can be estimated as 1/δ, if our interpretation is correct, this
would be equivalent to the time a variant is detectable in the field. In
our analysis, this would correspond to approximately four months to two
years on average. We recommend further work to evaluate how sampling
regime (i.e. sampling effort) could influence the parameter estimations
and affect these analyses.
An additional challenge that we faced in these analyses was that
geographical heterogeneity is likely to be an important element when
considering RHDV dynamics. In fact, it is believed that rabbit
susceptibility, due to heritable, acquired or abiotic factors, is not
homogeneous throughout Australia, which may causeR e values to vary greatly in different regions
(Henzell et al., 2003). Therefore, not all populations recover equally
from RHDV infections and field data suggest that climatic and other
environmental variables may also influence RHDV transmission (Henzell et
al., 2003). All these factors can make estimating meanR e values that adequately describe the overall
(i.e. nationwide) virus dynamic difficult because neither BDSKY or
SABDSKY cater for such further structuring. It would be particularly
useful if future studies improved our capacity to predict the extent and
the direction of possible biases when model assumptions are not met.
Models that take into account structuring (Kühnert et al., 2016) or more
complex disease dynamics (Volz, 2012) exist, and we attempted the
implementation of one such model. However, these analyses failed. It is
well known that analyses that take into consideration geographical
structuring need adequate sampling within each location (and possibly,
all the relevant demes sampled) and we argue that the lack of a rigorous
geographical coverage of the historic data available to us is the most
likely cause of the difficulties we encountered in the BDMM analysis.
Unfortunately, logistical and resource issues that limit the collection
of samples that would be suitable for these more complex models are a
common occurrence in wildlife research.
Both phylodynamic analyses (BDSKY and SABDSKY) showed a sharp decrease
in R e since 2010 and indeed, inspecting the
phylogenetic tree (Figure S3), we noted that most of post-2010 samples
had very long branches. As mentioned above, the mechanisms of the rabbit
recovery in the last 10 years are not quite clear. Field data indicate
that the frequency of RHDV outbreaks and the prevalence of the disease
in the monitored sites did not change over time (Mutze et al., 2014b,
Wells et al., 2018). Hence, this finding was somewhat unexpected, but
possibly points to an alternative mechanism that caused the recovery of
the rabbit populations. In 2010, weather changes caused by La Niña broke
a prolonged drought that started in 1996 and encompassed the whole
continent, known in Australia as ‘The Millennium Drought’. It is
possible that this drastic change in environmental conditions altered
RHDV epidemiology either directly or by altering vector activities,
reducing the long-distance transmission of the virus. The heavy rainfall
may also have increased the competition between the benign variant
(RCV-A1) and RHDV1 in the areas where this virus was active. It is known
that the former virus’ spread is facilitated in wet conditions (Liu et
al., 2014). Further development of rabbit genetic resistance to the
virus infection could also be a contributing factor, and has been
suggested as a possible explanation for reduced RHDV mortality rates and
recovery of rabbit populations in some parts of Australia {Mutze, 2014
#2683}. We considered RHDV1 sequences up until 2014 because it is
believed that alternative virus variants now known from Australia were
not circulating before this date. For this reason, we judge it highly
unlikely that new variants, especially RHDV2, which caused a substantial
rabbit mortality wave in 2015/16 {Ramsey, 2020 #3137}, would go
undetected and were not competing with RHDV1 between 2010 and 2014.
While all these are possible explanations of the quite drastic reduction
of R e since 2010, and it is possible that field
work focused in sites where RHDV1 remained active rather than being an
unbiased, representative selection across the whole country, these
remain to be proved.
Coalescent based
analyses
The coalescent based analyses (BSP and skyride) were very informative
and generally supported our current knowledge of the RHDV dynamics in
Australia. For example, both analyses correctly detected a sharp initial
increase in the virus population size, which subsequently slowed down
from 1997-1998 onwards. The further increase of the virus population
size between 2007 and 2009, detected by the coalescent-based analyses,
was unexpected. We note that strains isolated in the field in this
period were demonstrated to be of higher virulence compared to the
original strain released in Australia (Elsworth et al., 2014). It has
been hypothesised that virus selection tends to occur in clusters (Gog
and Grenfell, 2002) and it may be possible that this is what we are
witnessing here. Indeed, the VP60 capsid is under selection pressure
from the host immune system and it is possible that obtaining data from
other regions of the genomes could help to provide some evidence to this
end. However, we were not able to confirm this hypothesis with the data
available and it is also possible that other conditions that affect the
virus’ life cycle (such the ones mentioned above in regard to possible
changes in transmission rate) could artificially inflate the virus
population size.