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.