Since its introduction to control overabundant alien rabbits (Oryctolagus cuniculus), the highly virulent Rabbit Haemorrhagic Disease Virus (RHDV) has caused regular annual disease outbreaks in Australian rabbit populations. Although initially reducing rabbit abundance by 60%, continent-wide, experimental evidence has since indicated increased genetic resistance in wild rabbits that have experienced RHDV-driven selection. To identify genetic adaptations, which explain the increased resistance to this biocontrol virus, we investigated genome-wide SNP (single nucleotide polymorphism) allele frequency changes in a South Australian rabbit population that was sampled in 1996 (pre-RHD genomes) and after 16 years of RHDV outbreaks. We identified several SNPs with changed allele frequencies within or in proximity of genes that have roles potentially important for increased RHD resistance. Many of the identified genes are known to be involved in virus infections or immunity, or had previously been identified as being  differentially expressed in healthy vs. acutely RHDV-infected rabbits. Furthermore, we show in a simulation study that the allele/genotype frequency changes cannot be explained by drift alone, and that several candidate genes had also been identified as being associated with surviving RHD in a different Australian rabbit population. Our unique dataset allowed us to identify candidate genes for RHDV resistance that have evolved under natural conditions, and over a time span that would not have been feasible to study in an experimental setting. Moreover, it provides a rare example of host genetic adaptations to virus-driven selection in response to a suddenly emerging infectious disease.

Rabbit samples had been collected in 1996 from Gum Creek sheep station and abutting Flinders Ranges National Park in northern South Australia by CSIRO, Canberra, (CSIRO ethics approval 95/96-21) as part of a serological study to follow the epidemiology of RHDV. Samples taken included liver tissue for virus-capture ELISA (N = 59) to detect early stages of infection as well as eye lenses and body weight for age estimation (Cooke et al., 2000). The tissues used here were from rabbits which were too young to be fatally affected during the initial outbreak (Robinson et al., 2002) and those which had not yet been infected. They are termed the ‘pre-RHDV selection’ group.

In 2012, sixteen years after the first collection, more samples were collected by PIRSA, Adelaide (PIRSA ethics approval 11/09). Rabbits (N = 53) were collected in the same way by shooting at night with a 0.22 calibre rifle from a vehicle equipped with a 100 W spotlight. The location of each rabbit shot was recorded by GPS (Garmin 72, Schaffhausen, Switzerland) to nearest 10 m. A liver sample was collected from each rabbit using a sterile scalpel and stored in a labelled 5 ml Sarstedt vial and frozen. All samples were stored at -18o C. The 53 rabbits sampled in 2012 (post-RHDV selection) had been subject to approximately 15-16 generations under RHDV selection as every rabbit that is recruited into the breeding population would have had been previously exposed to RHD, and RHDV outbreaks had occurred regularly every year (Mutze et al., 2015, and additional unpublished serological data).

DNA was extracted using the Qiagen’s Gentra Puregene mouse tail kit (Hilden, Germany) according to the manufacturer’s protocol. Sixteen barcodes with a length of 4-9 nucleotides were designed using the GBS BARCODE GENERATOR (v.2.0, Deena Bioinformatics). GBS libraries were constructed using the restriction enzyme PstI (New England Biolabs) following a published protocol (Elshire et al., 2011) with the adjustment that we initially constructed an individual library for each rabbit (N=108). All individual libraries were quality controlled for fragment size distribution and absence of adapter dimers on an Agilent 2100 Bioanalyzer using the High Sensitivity Kit (Agilent, Waldbronn, Germany). We prepared seven sequencing libraries by always pooling 16 individual libraries at equal molarities. Sequencing (paired-end 100-bp reads) was conducted on an Illumina HiSeq2000 machine by the Beijing Genomic Institute (BGI), China.