The discovery of three new hare lagoviruses reveals unexplored viral diversity in this genus.

Our knowledge of mammalian viruses has been strongly skewed toward those that cause disease in humans and animals. However, recent metagenomic studies indicate that most apparently healthy organisms carry viruses, and that these seemingly benign viruses may comprise the bulk of virus diversity. The bias toward studying viruses associated with overt disease is apparent in the lagoviruses (family Caliciviridae) that infect rabbits and hares: although most attention has been directed toward the highly pathogenic members of this genus-rabbit haemorrhagic disease virus and European brown hare syndrome virus-a number of benign lagoviruses have also been identified. To determine whether wild European brown hares in Australia might also carry undetected benign viruses, we used a meta-transcriptomics approach to explore the gut and liver RNA viromes of these invasive animals. This led to the discovery of three new lagoviruses. While one was only detected in a single hare, the other two viruses were detected in 20 per cent of all animals tested. All three viruses were most closely related to other hare lagoviruses, but were phylogenetically distinct from both known viruses and from each other, indicating that lagoviruses have circulated for longer than previously assumed. Their evolution was also characterised by complex recombination events. Mapping mutations onto the lagovirus phylogeny revealed no amino acid changes that were consistently associated with virulence phenotype. Overall, our study points to extensive unsampled diversity in this genus, such that additional metagenomic studies are needed to fill gaps in the lagovirus phylogeny and better understand the evolutionary history of this important group of mammalian viruses.

1. Introduction

Although viruses likely infect all cellular organisms (Koonin, Dolja, and Krupovic 2015; Zhang, Shi, and Holmes 2018), their true diversity is underappreciated (Shi, Zhang, and Holmes 2018). Recent metagenomic studies have demonstrated that a wealth of viruses exist in apparently healthy species, and that their characterisation is central to understanding virus evolution and ecology (Li et al. 2015; Shi et al. 2016, 2018;). The study of such ‘background’ viruses in vertebrates is of particular interest as these may have the capacity to jump species boundaries and sometimes evolve new levels of virulence (Calisher et al. 2006; Domingo 2010; Watson et al. 2014; Geoghegan and Holmes 2018).

The genus Lagovirus (Caliciviridae) comprises both benign and highly pathogenic viruses (Capucci et al. 1996; Strive, Wright, and Robinson 2009; Le Gall-Recule et al. 2011; Cavadini et al. 2016; Le Pendu et al. 2017; Droillard et al. 2018). Lagoviruses are positive-sense, single stranded RNA viruses that infect members of the Leporidae family of mammals (i.e. rabbits and hares) (Wirblich et al. 1994; Strive, Wright, and Robinson 2009; Abrantes et al. 2012). Lagovirus genomes are approximately 7.5 kb in length and are made up of two open-reading frames (ORF), encoding a polyprotein that is proteolytically cleaved to produce the non-structural proteins and the major capsid protein, VP60, as well as a minor structural protein (Wirblich et al. 1994; Le Gall et al. 1996; Wirblich, Thiel, and Meyers 1996). Like other caliciviruses, lagoviruses possess a subgenomic RNA that is collinear with the 3’ end of the genomic RNA and encodes the structural genes (Meyers, Wirblich, and Thiel 1991; Wirblich et al. 1994; Clarke and Lambden 1997; Bull et al. 2005; Coyne et al. 2006). Viruses in this genus have been classified into two proposed genogroups: GI encompasses all viruses related to rabbit haemorrhagic disease virus (RHDV, GI.1), while GII includes viruses related to European brown hare syndrome virus (EBHSV, GII.1) (Le Pendu et al. 2017). Pathogenic viruses of the Lagovirus genus, such as RHDV, RHDV2 (GI.2), and EBHSV, primarily affect the liver, causing massive hepatic necrosis usually associated with high mortality, while benign viruses exhibit an intestinal tropism (Fuchs and Weissenbock 1992; Capucci et al. 1996; Strive, Wright, and Robinson 2009; Le Gall-Recule et al. 2011).

In Australia, extensive work has been done to characterise rabbit lagoviruses due to their importance in rabbit biocontrol. Indeed, rabbits are a major pest species in Australia, and since the mid-1990s RHDV has been used as a biocontrol for overabundant rabbit populations, and is still released periodically (Cooke and Fenner 2002; Abrantes et al. 2012). European rabbits and European brown hares (Lepus europaeus [L. europaeus]) were both successfully introduced into Australia in the 1800s and eventually reached plague densities, causing agricultural and ecological damage (Rolls 1984; Fenner 2010; Stott 2015). While rabbits successfully colonised most of the continent except the wet tropics and extremely arid zones (Myers and Parker 1965), hares did not spread as far, and were primarily located in the south-east of Australia (Stott 2015). For uncertain reasons, the hare population declined rapidly at the start of the 1900s (Rolls 1984), although they are still considered as pests in parts of Australia today (Agriculture Victoria, Victoria State Government 2017; Department of Primary Industries, Government of South Australia 2017).

Following its initial release, a reduced effectiveness of RHDV was noted in south-eastern temperate regions of Australia, and serological data suggested the existence of a related benign virus (Nagesha et al. 2000). Subsequently, a benign lagovirus, RCV-A1 (GI.4) was isolated (Strive, Wright, and Robinson 2009), and shown to confer a degree of cross-protection against RHDV, potentially interfering with biocontrol (Strive et al. 2010, 2013). Additionally, a number of seemingly benign lagoviruses were detected in Europe (RCV-E1, GI.3) and New Zealand (RCV-A1, GI.4) (Capucci et al. 1996; Le Gall-Recule et al. 2011; Le Pendu et al. 2017; Nicholson et al. 2017; Lemaitre et al. 2018). In the last decade, additional lagoviruses have appeared in Australia, including RHDVa-Aus (GI.4eP–GI.1a) (Mahar et al. 2018a), RHDV2 (GI.2) (Hall et al. 2015), and a recombinant of these two variants (GI.4eP–GI.2, [RdRp-capsid genotype]) (Le Pendu et al. 2017; Hall et al. 2018). RHDV2 has also been detected in hares (L. europaeus) in Australia (Hall et al. 2017), and although these infections were likely the result of spill-over from sympatric rabbit populations (Mahar et al. 2018b), recent reports suggest hare-to-hare transmission of RHDV2 may occur (Neimanis et al. 2018). Apart from RHDV2, no other lagoviruses have been detected in hares in Australia. In Europe, L. europaeus is affected by the pathogenic lagovirus EBHSV (Duff et al. 1994; Frolich et al. 2001; Billinis et al. 2005; Syrjala, Nylund, and Heinikainen 2005; Le Gall-Recule et al. 2006; Frolich et al. 2007), as well as RHDV2, which is highly virulent in multiple hare species (Puggioni et al. 2013; Camarda et al. 2014; Le Gall-Recule et al. 2017; Velarde et al. 2017; Mahar et al. 2018b; Neimanis et al. 2018). Two additional, presumably benign, hare caliciviruses from Europe (denoted GII.2) have recently been reported (Cavadini et al. 2016; Le Pendu et al. 2017; Droillard et al. 2018).

Less emphasis has been placed on the investigation of viruses in hares than to those affecting rabbits, likely because hares are not commercially farmed and have a comparatively lower impact as a pest species in Australia. However, exploration of the hare virome is of importance, not just for understanding the biological diversity of lagoviruses, but for understanding aspects of the evolution of virulence (Geoghegan and Holmes 2018). We therefore aimed to characterise the RNA virome of healthy hares in Australia, and in so doing increase our understanding of RNA virus evolution.

2. Materials and methods

2.1 Sample collection

Liver and duodenum were collected post-mortem from apparently healthy (shot) hares and frozen at −20 °C. Samples were taken from thirty hares in Hamilton, Victoria (VIC) (−37.717168, 142.026622) over two nights—30 June 2016 and 23 May 2017—and twelve hares from Mulligan’s Flat, Australian Capital Territory (ACT) (−35.166512, 149.164925) in December 2012, February 2016, or May–July 2016 (Table 1). Samples were collected as part of a routine vertebrate pest control programme and lagovirus surveillance studies. All work was carried out according to the Australian Code for the Care and Use of Animals for Scientific Purposes with approval from the institutional animal ethics committee (ESAEC 12-15 and CWLAAEC 16-02).

Table 1.

Prevalence of new hare lagoviruses and sequencing library details.

	Sample collection		Sequencing librarya		Virus detection by PCR in duob
Name	Locationc	Date	Liver	Duodenum	HaCV-A1	HaCV-A2	HaCV-A3
JM-1	Ham	30 June 2016	N/A	N/A	−	−	−
JM-2	Ham	30 June 2016	N/A	JM-2-duo	+	−	−
JM-3	Ham	30 June 2016	N/A	N/A	−	−	−
JM-4	Ham	30 June 2016	N/A	N/A	−	+	−
JM-5	Ham	30 June 2016	N/A	N/A	−	−	−
JM-6	Ham	30 June 2016	N/A	N/A	−	+	−
JM-7	Ham	30 June 2016	N/A	N/A	−	−	−
JM-8	Ham	30 June 2016	N/A	N/A	−	−	−
JM-9	Ham	30 June 2016	N/A	N/A	−	+	−
JM-10	Ham	30 June 2016	N/A	N/A	−	−	−
JM-11	Ham	30 June 2016	N/A	N/A	−	+(L)	−
JM-12	Ham	30 June 2016	N/A	N/A	−	+	−
JM-13	Ham	30 June 2016	N/A	N/A	−	−	−
JM-14	Ham	30 June 2016	N/A	N/A	−	−	−
JM-15	Ham	30 June 2016	N/A	N/A	−	−	−
JM-16	Ham	30 June 2016	N/A	N/A	+	−	−
JM-17	Ham	30 June 2016	N/A	N/A	−	−	−
JM-18	Ham	30 June 2016	N/A	N/A	−	−	−
JM-19	Ham	30 June 2016	N/A	N/A	−	−	−
JM-20	Ham	30 June 2016	N/A	N/A	−	+	−
JM-22	Ham	23 May 2017	Ham1-L	Ham1-D	+	−	−
JM-24	Ham	23 May 2017	Ham1-L	Ham1-D	+	+	−
JM-26	Ham	23 May 2017	Ham1-L	Ham1-D	−	+/−	−
JM-27	Ham	23 May 2017	Ham2-L	Ham2-D	+/−	−	−
JM-29	Ham	23 May 2017	Ham2-L	Ham2-D	+	−	−
JM-30	Ham	23 May 2017	Ham3-L	Ham3-D	+	−	−
JM-31	Ham	23 May 2017	Ham3-L	Ham3-D	−	−	−
JM-34	Ham	23 May 2017	Ham3-L	Ham3-D	−	+(L)	−
JM-35	Ham	23 May 2017	Ham4-L	Ham4-D	+	−	−
JM-40	Ham	23 May 2017	Ham4-L	Ham4-D	+(L)	−	−
MF-01	MF	20 December 2012	N/A	N/A	−	−	−
MF-02	MF	20 December 2012	N/A	N/A	−	+	−
MF-07	MF	20 December 2012	MF1-L	MF1-D	−	−	−
MF-22	MF	20 December 2012	MF1-L	MF1-D	−	−	−
MF-137	MF	3 February 2016	MF1-L	MF1-D	−	−	−
MF-148	MF	5 May 2016	MF2-L	MF2-D	−	−	−
MF-149	MF	5 May 2016	MF2-L	MF2-D	−	−	−
MF-150	MF	9 June 2016	MF3-L	MF3-D, MF-150	−	−	+
MF-151	MF	9 June 2016	MF3-L	MF3-D	−	−	−
MF-152	MF	9 June 2016	MF3-L	MF3-D	−	−	−
MF-155	MF	7 July 2016	MF4-L	MF4-D	−	−	−
MF-156	MF	27 July 2016	MF4-L	MF4-D	−	−	−

N/A, not applicable—RNA sequencing was not performed on these samples.

duo, duodenum; +, positive; −, negative; +/−, weak positive; (L) weak positive in liver RNA.

Ham, Hamilton, Victoria; MF, Mulligan’s Flat, Australian Capital Territory.

2.2 RNA isolation

RNA was extracted from 20–30 mg of tissue using the Maxwell 16 LEV simplyRNA tissue kit and extraction robot (Promega, WI, USA) as per the manufacturer’s instructions.

2.3 cDNA synthesis and PCR for detection of diverse lagoviruses

First-strand cDNA was prepared using Invitrogen SuperScriptTM IV Reverse Transcriptase (Thermofisher Scientific, MA, USA) according to the manufacturer’s instructions using 5 μl of RNA and 500 ng of Oligo(dT)(18mer) or 10 μM CaVuniR specific primer (Supplementary Table S1). For cDNA prepared for 3’ end amplification, 10 μM of primer GV270 (Eden et al. 2013) was used. Duodenum samples (n = 38) were screened for the presence of lagoviruses using a universal lagovirus PCR as described previously (Strive, Wright, and Robinson 2009), and positive amplicons were confirmed by Sanger sequencing at the Australian Cancer Research Foundation (ACRF) Biomolecular Resource Facility (BRF) in Canberra, ACT.

2.4 Initial amplification of HaCV-A3 for Sanger sequencing

Regions of the genome of HaCV-A3 (sample MF-150) were initially amplified using EBHSV primers EBHSV_VP60_01728R and EBHSV_VP60_0813F (Lopes et al. 2014), or specifically designed broadly-reactive primers (Supplementary Table S1). PCRs were conducted using Invitrogen Platinum Taq Polymerase High Fidelity kit (Thermofisher Scientific) according to the manufacturer’s protocol, using 1.5 μl of diluted cDNA (1:2) as template in a 25 μl reaction. Positive amplicons were sequenced at ACRF-BRF.

2.5 RNA sequencing

2.5.1 RNA library construction and sequencing

RNA from two hare duodenum samples that tested positive in the universal lagovirus PCR (MF-150 and JM-2), as well as RNA from the liver and duodenum of ten additional hares from Hamilton VIC, and ten hares (including MF-150) from Mulligans Flat ACT, were selected for RNA sequencing (Table 1). Freshly extracted RNA was treated using Invitrogen TURBO DNase (Thermofisher Scientific) and further purified and concentrated using the RNeasy MinElute cleanup kit (Qiagen, Hilden, Germany). RNA was quantified using the Qubit RNA Invitrogen Broad-range Assay kit with the Qubit Fluorometer v3.0 (Thermofisher Scientific), and further quantified and assessed for quality using the Agilent RNA 6000 nano kit and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). JM-2 and MF-150 duodenum RNA were each submitted for sequencing as a single library, while the remaining RNA samples were pooled in equal proportions by location and tissue type into pools of two–three individuals, totalling eighteen libraries (Table 1). MF-150 duodenum RNA was sequenced individually, and was also included in the Mulligan’s Flat duodenum pools, as the virus loads initially detected in MF-150 duodenum RNA appeared very low. Library preparation and sequencing was carried out at the Australian Genome Research Facility (Melbourne) using the TruSeq total RNA library preparation kit (Illumina, CA, USA) with host rRNA depletion using the Illumina Ribo-Zero-Gold rRNA removal kit (Epidemiology). Paired-end sequencing (100 bp) was performed on the HiSeq 2500 sequencing platform.

2.5.2 Contig assembly and annotation

Reads were trimmed using Trimmomatic (Bolger, Lohse, and Usadel 2014) and assembled into contigs de novo using Trinity (Grabherr et al. 2011). Abundance (as expected counts) was estimated for each contig using the RSEM tool (Li et al. 2010), an alignment-based quantification method implemented in Trinity. BLASTn and DIAMOND BLASTx were used to compare Trinity contigs to the National Center for Biotechnology Information (NCBI) nucleotide (nt) database (e-value cut-off 1 × 10−10) and non-redundant protein (nr) database (e-value cut-off 1 × 10−5), respectively. Results were filtered and contigs that had a viral hit for either BLAST search were retained. Virus-host associations were allocated using the Virus-Host database (https://www.genome.jp/virushostdb). All reads were mapped to host rRNA (rabbit rRNA sequences were used as hare rRNA sequences were not available) using Bowtie2 (Langmead and Salzberg 2012) to quantify remaining host rRNA reads. The rabbit host rRNA target index was generated from a complete Oryctolagus cuniculus (O. cuniculus) 18s rRNA reference sequence obtained from GenBank (accession NR_033238) and a near complete O. cuniculus 28s rRNA sequence obtained from the Silva high quality ribosomal database (Quast et al. 2013) (accession GBCA01000314).

2.5.3 Lagovirus genome assemblies

The lagovirus contigs from all libraries, in addition to Sanger sequences obtained for JM-2 and MF-150 duodenum samples, were aligned, and contigs with overlapping regions were merged using the Geneious assembler (Kearse et al. 2012). Four merged contigs were generated with lengths 7,364 nt, 5,588 nt, 1,628 nt, and 1,375 nt. The three longest merged contigs were generated from contigs from the Hamilton, VIC libraries, while the shortest contig was compiled from Sanger sequences from MF-150 duodenum (Mulligan’s Flat, ACT). To generate library-specific lagovirus contigs, the consensus sequences of these four merged contigs were used as reference sequences and reads from each individual library were mapped to these reference sequences using Bowtie2 (Langmead and Salzberg 2012). Consensus sequences were extracted from the library-specific contigs. To obtain more of the genome sequence of the lagovirus detected in MF-150, reads and contigs assembled from the MF-150 duodenum library and the Mulligan’s Flat library containing MF-150 duodenum (MF3-D) were aligned to an EBHSV reference sequence (KC832839.1/EBHSV/SWE/O4021-9/1982). Two contigs and six reads aligned to EBHSV in regions upstream of the MF-150 Sanger sequence already obtained, allowing the design of primers to amplify and sequence across some missing regions.

2.6 Genome confirmation and extension PCRs

Further primers were designed to amplify missing parts of the newly identified hare calicivirus genomes, to merge the 5,588 nt and 1,628 nt contigs, and confirm the genome sequence for HaCV-A1 (JM-29 duodenum) by amplicon sequencing (Supplementary Table S1). Primer GV271 (Eden et al. 2013) was used for 3’ end amplification from within the polyA tail. PCRs were conducted using Invitrogen Platinum Taq Polymerase High Fidelity kit according to the manufacturer’s protocol, using 2.4 μl of cDNA as template in a 40 μl reaction. DNA libraries were prepared and sequenced using either Illumina Miseq technology as described previously (Eden et al. 2015; Mahar et al. 2016) or Sanger sequencing conducted at ACRF-BRF.

2.7 Hare calicivirus screening PCRs

To determine the prevalence of each of the new hare caliciviruses, hare duodenum RNA was individually screened for each new lagovirus. Specific primer sets were designed based on the sequence of the three new lagoviruses to enable detection of each virus; HaCV-A1, HareCaV1_F6.2, and HareCaV1_R6.5 (330 bp amplicon); HaCV-A2, HareCaV2_F5.2, and HareCaV2_R5.4 (213 bp amplicon); HaCV-A3, HareCaV4_F5.5, and HareCaV4_R5.9 (408 bp amplicon) (Supplementary Table S1). Liver RNA was also screened from hares for which a product was amplified from the duodenum RNA. RT-PCR was conducted using the One-Step Ahead RT-PCR kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions using 1 μl of RNA diluted 1:10 in nuclease free water in a 10 μl reaction volume.

2.8 Identification of conserved domains and potential ORFs

The NCBI Conserved domains tool (Marchler-Bauer et al. 2017) was used to check for the presence of conserved functional domains, and the ExPASY translate tool (https://web.expasy.org/translate/) and the Geneious ORF prediction tool were used to identify ORF. The location of ORF 1 in all three viruses was inferred due to its size (largest ORF) and through sequence similarity with other lagoviruses (accession NC_002615, NC_001543). The location of ORF 2 (in HaCV-A1 and HaCV-A2) was chosen based on the location of putative termination upstream ribosomal binding site (TURBS) motifs (identified manually), as translation re-initiation in caliciviruses tends to occur within 12–24 nt of the TURBS structure (Luttermann and Meyers 2014), and only one potential ORF fit this criteria. The likely cleavage fragments of the ORF 1 polyprotein were inferred from sequence homology with EBHSV and RHDV (accession NC_002615, NC_001543).

2.9 Recombination analyses

The RDP, GENECONV, MAXCHI, and BOOTSCAN methods within the RDP4 package (Martin et al. 2015) were used to screen for recombination on a data set containing the three new sequences plus eighteen non-recombinant lagovirus sequences. A P-value of 0.05 represented a significant result for all tests, and putative recombination events were considered to be those detected by at least two methods and confirmed manually. A pairwise identity plot was generated with a sliding window of 30 nt. To confirm recombination events, we inferred phylogenetic trees on sections of the alignment either side of the putative recombination breakpoint using a maximum likelihood (ML) approach as described below. Significant evidence for recombination was reported as cases of clear phylogenetic incongruence with strong (i.e. >70%) bootstrap support.

2.10 Phylogenetic analysis

The lagovirus genome sequences identified here were aligned with twenty-eight (twenty-four complete genomes, four capsid sequence only) sequences available on GenBank, representing the known diversity of lagoviruses, using MAFFT (Standley and Katoh 2013). ML phylogenetic trees were inferred for both the RdRp (1,548 nt, twenty-seven sequences) and capsid (1,704 nt, thirty-one sequences) using PhyML (Guindon et al. 2010) employing the GTR+Γ+I model of nucleotide substitution (as selected using jModelTest v2.1.6) (Darriba et al. 2012) with four rate categories, an estimated proportion of invariant sites and gamma distribution parameter. Topology searching used a combination of nearest-neighbor interchange and subtree pruning and regrafting branch-swapping. Branch support was estimated using 1,000 bootstrap replicates using the same ML procedure as described above.

2.11 Detection of virulence determinants and selection analysis

RdRp and capsid protein gene data sets containing all available lagovirus sequences were obtained from NCBI and aligned with the HaCV-A1 sequences described here, using MAFFT (RdRp n = 628, capsid n = 1,443). Sequences were visually inspected to identify amino acid sites that differed between benign and virulent viruses. Where information on the tissue tropism and virulence status was unavailable, viruses named ‘EBHSV’, ‘RHDV’, or ‘RHDV2’ were assumed to be virulent with a liver tropism, while those named ‘RCV’ or ‘HaCV’ were assumed benign with intestinal tropism (Le Pendu et al. 2017).

To screen for positive selection at individual sites, branches or entire genes, the data sets prepared for the phylogenies were used, along with a curated version of the nucleotide alignments of all available RdRp and capsid lagovirus sequences. For the curated lagovirus data sets, incomplete sequences were removed and the capsid data set was trimmed to the length of HaCV-A3 (1,704 nt). To make computational analyses viable, duplicate sequences were removed from the capsid data set and sequences were randomly sampled to produce final data sets of RdRp n = 347, capsid n = 500. ML phylogenies were inferred as described above. Selection analyses were conducted using the Datamonkey online server (Li et al. 2018), employing user-defined trees and the aBSREL (Smith et al. 2015) (branch, small data sets only), RELAX (Murrell et al. 2014) (branch), FEL (Kosakovsky Pond, and Frost 2005) (site-specific), and BUSTED (Murrell et al. 2015) (gene-wide) methods. For the branch and gene-wide methods, selection analyses were run on the entire phylogenies as well as on the designated virulent or benign lineages separately. Selection results that were not consistent between the small and large data sets were not considered.

3. Results

3.1 Initial PCR for lagovirus detection

As an initial screen to identify diverse lagoviruses a broad-range universal lagovirus RT-PCR (Strive, Wright, and Robinson 2009) was used to analyse hare duodenum samples (n = 38) collected in two locations Hamilton, VIC and Mulligan’s Flat, ACT. Lagoviruses were detected in two hare duodenum samples: MF-150 and JM-2. This was confirmed by Sanger sequencing, and initial phylogenetic analysis of the ∼300 nt sequences indicated that these viruses were distinct from known lagoviruses and each other (data not shown).

3.2 RNA sequencing

RNA sequencing (i.e. ‘meta-transcriptomics’ (Shi et al. 2016)) was used to characterise the complete genomes of the JM-2 and MF-150 duodenum RNA and a selection of other hare duodenum and liver RNA pooled into sixteen libraries (Table 1). An aggregate of 792,023,978 reads were obtained for all libraries, 605,624,790 (76%) of which did not map to host rRNA, averaging 38,086,664 non-rRNA reads per library. No viral contigs were assembled from any of the liver libraries or most of the Mulligan’s Flat duodenum libraries. The four Hamilton duodenum pooled libraries, one Mulligan’s Flat duodenum pooled library (MF3-D), and the JM-2 duodenum library together had a total of fifty-six viral contigs. These contigs matched six different viruses in a BLAST analysis; EBHSV, Hare calicivirus, Hubei partiti-like virus 54, Hubei partiti-like virus 49, Mammalian orthoreovirus 1, and Mammalian orthoreovirus 3. Contigs with highest identity to lagoviruses (either EBHSV [GII.1] or hare calicivirus [GII.2]) were present in all of these libraries and, on average, made up 86 per cent (71–100%) of all viral contigs in each library. However, lagovirus reads were in very low abundance overall, comprising less than 0.003 per cent of the non-rRNA transcriptome in each library. The fifty lagovirus contigs had an average nucleotide identity of only 84.7–89.9 per cent to the top BLAST result, indicating potential new viruses. The remaining viral contigs (those other than lagoviruses) were short (range 214–483 nt), and all had very low abundance of less than twenty reads.

3.3 Lagovirus genome assembly and annotation

The near complete genomes of two new lagoviruses, termed Hare calicivirus Australia-1 (HaCV-A1) and Hare calicivirus Australia-2 (HaCV-A2), were assembled from RNA sequencing libraries, Ham-2D (mean coverage 22.8X) and Ham-1D (mean coverage 12.8X), respectively, and extended through amplicon or Sanger sequencing (Fig. 1A). The partial genome sequence of a third new virus, Hare calicivirus Australia-3 (HaCV-A3), was obtained by Sanger and amplicon sequencing of MF-150 duodenum using a combination of EBHSV primers (Lopes et al. 2014), specifically designed broadly-reactive primers, and primers based on RNA sequencing reads from library MF-3D. Excluding the polyA tail, and primer inferred sequence, consensus sequences of 7,386 nt, 7,412 nt, and 4,570 nt were obtained for HaCV-A1, HaCV-A2, and HaCV-A3, respectively (Fig. 1). The 3’ untranslated region (UTR) was 64 nt for HaCV-A1 and 77 nt for HaCV-A2, and was not obtained for HaCV-A3. The complete 5’ end was not obtained, but only 14 nt were likely missing from the start of HaCV-A1 and 24 nt likely missing from HaCV-A2 based on inference from other lagoviruses.

Figure 1.

Genome structure and pairwise identity of new viruses. (A) The region of the genome sequenced for each new virus is represented schematically. ORFs are represented by coloured arrow bars (blue, HaCV-A1; red, HaCV-A2; green, HaCV-A3). Conserved protein domains detected using the NCBI conserved domains search tool are indicated by dark grey boxes (RdRp, RNA-dependent RNA polymerase; DUF840, lagovirus protein of unknown function). The likely cleavage fragments/peptides of the ORF 1 polyprotein, inferred from sequence homology with EBHSV and RHDV, are indicated by the light grey arrow bars (2C-like, 2C-like RNA helicase; VPg, genome-linked viral protein; 3C-like, 3C-like proteinase; RdRp, RNA-dependent RNA polymerase; VP60, major capsid protein). Amino acids flanking the likely cleavage sites in the polyprotein are indicated at the junction between the peptides using amino acid one-letter identifiers. A broken appearance at either end of the arrow bars indicates incomplete sequence for that ORF or peptide. The 3’ UTR and polyA tails (A(n)) are indicated where sequence was obtained. Genome numbering at regular intervals is indicated above each schematic. The sequence coverage for the initial assemblies of HaCV-A1 and HaCV-A2, from Ham-2D and Ham-1D libraries respectively, is shown beneath the genome schematics. Coverage is shown for each position along the length of the genome. The regions obtained/confirmed by Sanger sequencing for HaCV-A2 are indicated by the red bars underneath the coverage plot. (B) Pairwise nucleotide identity (y-axis) according to genome position (x-axis) is plotted using a sliding window of 30 nt. The plot was generated in the RDP4 program from an alignment trimmed to the length of HaCV-A3. A clear cross-over between the blue line (identity between HaCV-A3 and HaCV-A1) and the red line (identity between HaCV-A3 and HaCV-A2) suggests that HaCV-A3 is a recombinant between parental viruses related to HaCV-A1 and HaCV-A2. The cross-over event occurs near the junction of the RdRp and capsid and is indicated by a black arrow.

All three viruses have the same genome organisation as other lagoviruses. HaCV-A1 and HaCV-A2 have two ORFs—one encoding a polyprotein containing the non-structural genes and capsid gene, and one encoding the minor structural protein (Fig. 1A). The polyprotein encoded by ORF 1 is likely 2,332 amino acids in length for HaCV-A1 (based on start codon in the primer inferred sequence) and is likely the same for HaCV-A2, although the start of the coding sequence was not obtained for the latter virus. The likely cleavage products resulting from post-translational processing of the ORF 1 polyprotein (and cleavage sites) were inferred from sequence similarity with EBHSV and RHDV. These included peptides for which conserved domains were identified (RNA helicase, peptidase C37/3C-like proteinase, RNA-dependent RNA polymerase (RdRp), calicivirus capsid protein, and DUF840, a lagovirus protein of unknown function), as well as the genome-linked viral protein, VPg, which binds to the 5’ end of calicivirus RNA molecules (Meyers, Wirblich, and Thiel 1991); and three proteins with unknown function (Fig. 1A). Although much of the HaCV-A3 genome sequence is missing, the cleavage sites and predicted products for the region obtained match those of HaCV-A1 and HaCV-A2 (Fig. 1A). Potential TURBS motifs were identified at positions 6,904–6,908 (motif 2*), 6,910–6,916 (motif 1), 6,959–6,963 (motif 2) of the HaCV-A1 partial genome sequence and at the equivalent location in the HaCV-A2 partial genome sequence; 6,907–6,911 (motif 2*), 6,913–6,919 (motif 1), and 6,962–6,966 (motif 2). For both viruses, based on the putative location of the TURBS motifs (Royall and Locker 2016), the second ORF is likely to overlap with the first ORF by 8 nt, as seen for EBHSV (Wirblich et al. 1994), and encode a protein of 113 amino acids.

To validate our genome assembly, we confirmed the HaCV-A1 genome sequence by amplicon sequencing of an individual sample (JM-24). There were only ten nucleotide differences between the amplicon-based sequencing approach and the de novo RNA sequencing assembly approach. This level of variation is expected given that amplicon sequencing was conducted on a single sample while RNA sequencing was conducted on a pool of three samples (Ham-2D), all of which were positive for HaCV-A1.

3.4 Evolutionary relationships of hare lagoviruses

The three new hare lagoviruses were genetically diverse, with HaCV-A1 sharing only 74 and 77 per cent nt identity across sequenced regions with HaCV-A2 and HaCV-A3, respectively, while HaCV-A2 and HaCV-A3 share 78 per cent nt identity. Notably, HaCV-A3 was most similar to HaCV-A2 in the RdRp gene and most similar to HaCV-A1 in the capsid region, suggestive of recombination.

Phylogenetic analysis revealed that the three new hare lagoviruses form a strongly supported monophyletic group with other hare-specific lagoviruses (GII)—EBHSV and hare caliciviruses—in both the capsid and RdRp gene trees (Fig. 2). However, all three viruses are clearly distinct, separated both from each other and other known lagoviruses by relatively long branches with strong bootstrap support (Fig. 2). Indeed, there is approximately the same phylogenetic distance between HaCV-A1 and HaCV-A2 as among all known rabbit lagoviruses. According to the proposed classification guidelines (Le Pendu et al. 2017) it is likely that each of these viruses represent a new genotype within the GII genogroup of the Lagovirus genus, provisionally termed GII.3 (HaCV-A1), GII.4 (HaCV-A2), and GII.5 (HaCV-A3).

Figure 2.

Phylogenetic analysis of new lagoviruses. Maximum likelihood phylogenies of the (A) RdRp gene (n = 27; 1,548 nt) and (B) capsid gene (n = 31; 1,704 nt) were inferred for the three new lagoviruses along with representative members of the genus Lagovirus. The accession number of sequences obtained from GenBank is shown in the taxa labels. Trees were mid-point rooted for clarity only, and branch support was estimated using 1,000 bootstrap replicates, which are shown at the major nodes. The taxa names of the three benign hare viruses reported in this study are bolded. The two major clades in each phylogeny are labelled according to proposed genogroup and prototypical host. Proposed genotypes are indicated in taxa labels (those without a genotype label are unclassified). The branches within and leading to virulent clades are coloured red, while branches within and leading to benign clades are coloured blue. The branch to MRCV, a virulent virus that only caused local outbreaks and apparently died out, is coloured orange.

In addition, there is marked phylogenetic incongruence between the RdRp and capsid protein gene trees among the hare caliciviruses (Australian and European). In the capsid gene tree (Fig. 2B), HaCV-A2 clusters most closely with the two European benign hare caliciviruses (GII.2), while the other two new viruses cluster together and share a more recent common ancestor with the pathogenic hare lagovirus, EBHSV. In the RdRp phylogeny (Fig. 2A), the two European hare caliciviruses cluster most closely with EBHSV and the three viruses discovered here are more distant.

3.5 Lagovirus recombination

To determine the cause of the incongruence between the RdRp and capsid gene phylogenies we performed a more detailed analysis of this putative recombination event. This revealed that HaCV-A3 is likely to be a recombinant of viruses related to HaCV-A1 and HaCV-A2 (RDP pairwise distance plot; Fig. 1B). The estimated location of the putative breakpoint was at position 2,909 in the HaCV-A3 sequence (99% CI: 2, 827–3,183), nineteen nucleotides downstream of the RdRp/capsid putative cleavage site (Fig. 1B), which is the equivalent of position 5,307 on the reference EBHSV genome sequence (accession NC_002615). Recombination at this location essentially divides the genome into a region encoding the non-structural proteins and a second region encoding the structural proteins (Fig. 1B). Phylogenetic analysis on regions either side of the breakpoint strongly supported recombination, with HaCV-A3 clustering with HaCV-A2 in the non-structural genes tree and clustering with HaCV-A1 in the capsid tree, with robust bootstrap support (>70%). This incongruence is captured in the RdRp and capsid phylogenies presented in Fig. 2. However, given the substantial diversity of potential parent sequences and lack of sampling in this clade, it was difficult to determine the evolutionary history of recombination events with certainty, and an end breakpoint could not be identified. Hence, HaCV-A1 or HaCV-A2 may be the recombinant, or more than one recombination event may have occurred among these and related GII viruses. Importantly, the regions flanking the putative recombination breakpoint were amplified in a single amplicon for each of the Australian hare caliciviruses, excluding the possibility of incorrect assembly.

3.6 Identification of virulence determinants and selection pressures

We manually explored amino acid sequence alignments of all available lagovirus capsid and RdRp sequences to identify putative virulence determinants within these proteins that may be associated with the pathogenic lagovirus lineages (Fig. 2). Two sites were identified in the capsid protein—289 and 455 in EBHSV capsid sequence NC_002615.1—that differed between all benign sequences and virulent sequences, with the exception of one RCV sequence in each case that contained the putative virulence mutation. A selection analyses conducted on the RdRp and capsid data sets revealed neither sites nor branches under positive selection associated with differences in virulence among the lagoviruses sampled here.

3.7 Prevalence of new hare lagoviruses

Screening PCRs amplifying a short region of the capsid gene were designed for each new virus, and forty-two duodenum samples from healthy hares from both locations were tested for the presence of these viruses. HaCV-A1 and HaCV-A2 were both found at a prevalence of 30 per cent in Hamilton, VIC, and one hare was infected with both viruses (Table 1). HaCV-A1 was not detected in Mulligan’s Flat, ACT, while HaCV-A2 was detected in 1/12 hares in this location, making it the only virus of the three newly identified lagoviruses to be detected at both locations (Table 1). HaCV-A3 was detected in only one hare duodenum, MF-150, and was likely to be present at a very low concentration, as viral contigs could not be assembled from a total of 39,993,840 reads from RNA sequencing of this sample. For PCR-positive duodenum samples, liver samples from the same individuals were also screened for the presence of the new lagoviruses. HaCV-A1 was detected in one liver sample and HaCV-A2 was detected in two liver samples (Table 1), although the amplicons were faint, indicating a lower viral abundance compared to the duodenum or possible contamination during sample collection.

4. Discussion

Using a bulk RNA sequencing approach we identified three new viruses in European hares present in Australia: HaCV-A1, HaCV-A2, and HaCV-A3. Prior to this, the only lagovirus detected in hares in Australia was RHDV2 (Hall et al. 2017), which is primarily a rabbit virus, and phylogenetic evidence suggests that RHDV2 infection in hares in Australia occurred as a result of transient spill-over events (Mahar et al. 2018b). While HaCV-A3 was only found in one animal, the other two new viruses were both detected in almost one third of hares tested at the Hamilton, VIC site, during both sampling periods (one year apart). This suggests that similar to the non-pathogenic rabbit calicivirus RCV-A1 (GI.4) (Liu et al. 2014), these two viruses may be prevalent in certain populations, although more extensive screening is needed.

The three new viruses were all members of the genus Lagovirus (family Caliciviridae). This genus comprises both virulent and benign viruses that infect hares (Lepus) and rabbits (O. cuniculus) (Strive, Wright, and Robinson 2009). The virulent viruses, RHDV (GI.1), RHDV2 (GI.2), and EBHSV (GII.1), have a liver tropism and are associated with necrotic hepatitis often resulting in fatality (Chasey and Duff 1990; Fuchs and Weissenbock 1992; Le Gall-Recule et al. 2013), while the benign viruses, RCV-A1 (GI.4), RCV (unclassified), RCV-E1 (GI.3), and the French hare calicivirus (GII.2), have an intestinal tropism (Capucci et al. 1996; Strive, Wright, and Robinson 2009; Le Gall-Recule et al. 2011; Droillard et al. 2018). All three new viruses discovered here were found in low abundance in the duodenum of apparently healthy hares, consistent with the intestinal tropism observed for benign lagoviruses (Capucci et al. 1996; Strive, Wright, and Robinson 2009; Droillard et al. 2018). HaCV-A1 or HaCV-A2 was also detected in the liver of three of these apparently healthy hares, although at objectively lower levels, suggesting that the site of replication is likely to be the intestine.

There was evidence of recombination between the three new viruses, confirming its role in generating lagovirus genetic diversity (Lopes et al. 2015; Hall et al. 2018). The location of the putative breakpoint is near the junction of the RdRp and capsid protein genes. Recombination in this region results in chimeric viruses with non-structural genes derived from one virus and structural genes derived from another (Bull et al. 2005; Coyne et al. 2006; Lopes et al. 2015; Hall et al. 2018). This is a recombination hotspot in caliciviruses, with similar events observed between members of the GI lineage of lagoviruses, where several recombinants between RHDV (GI.1), RHDV2 (GI.2), and RCV-A1 (GI.4) have been reported (Lopes et al. 2015; Hall et al. 2018). The sequence in the RdRp/capsid junction is highly conserved in caliciviruses, and is predicted to form a stem loop structure that may facilitate a pause in replication and a subsequent template switch (Bull et al. 2005; Coyne et al. 2006). In addition, several caliciviruses possess a subgenomic RNA encoding the structural genes which may serve as a secondary template for re-initiation of RNA synthesis (Meyers, Wirblich, and Thiel 1991; Clarke and Lambden 1997; Bull et al. 2005; Coyne et al. 2006). Due to the diversity between the three new Australian hare caliciviruses and related European hare caliciviruses, and general under-sampling of the Lagovirus GII clade, it is difficult to establish which of the new viruses is the recombinant and which is the parent, although our analysis suggested that HaCV-A3 was the most likely recombinant. Indeed, the pattern of phylogenetic incongruence between the RdRp and capsid gene trees may mean that several recombination events have occurred among these and related GII strains (i.e. the two European ‘GII.2’ hare caliciviruses, which also show incongruence), although analysis is complicated by the substantial sequence divergence, under-sampling of potential parental sequences, and/or older events being over-written by more recent recombination events (Martin et al. 2015).

The viruses identified here were most closely related to the only two known hare specific lagoviruses (GII), EBHSV (GII.1), and the recently reported European hare calicivirus (denoted GII.2) (Cavadini et al. 2016; Le Pendu et al. 2017; Droillard et al. 2018). The three new viruses are strikingly distant from the previously characterised hare lagoviruses and from each other, and each would constitute a new genotype of the Lagovirus genus (Le Pendu et al. 2017), provisionally GII.3 (HaCV-A1), GII.4 (HaCV-A2), and GII.5 (HaCV-A3). The addition of these three viruses has therefore greatly increased the phylogenetic diversity of this genus. The relatively large genetic distance both among these viruses and between the hare and rabbit viruses likely reflects a lack of sampling, with the possibility that lagoviruses may in fact infect a more diverse range of taxa, such as Sylvilagus sp. It should be noted that the two European hare caliciviruses, both tentatively denoted GII.2 (Cavadini et al. 2016; Le Pendu et al. 2017), cluster together, but likely exhibit sufficient diversity (85% nucleotide identity in the capsid protein gene) to be classified as different genotypes (Fig. 2).

Historically, research efforts have been focussed towards the highly virulent viruses due to their apparent impacts, although non-virulent lagoviruses are increasingly being reported (Capucci et al. 1996; Strive, Wright, and Robinson 2009; Le Gall-Recule et al. 2011; Le Pendu et al. 2017; Nicholson et al. 2017). Indeed, with additional comprehensive sequencing studies, the Lagovirus genus may prove to be largely comprised of asymptomatic viruses. These viruses may provide an extensive gene pool for recombination or even the possible emergence of additional virulent lagoviruses like RHDV2. Although mid-point rooted, the lagovirus phylogenies generated here are compatible with the idea that the common ancestor of these viruses was benign, and that high virulence evolved independently in EBHSV, RHDV, and RHDV2 (Fig. 2) (Mahar et al. 2016). That antibodies to lagoviruses were detected in the serum of healthy rabbit populations before the emergence of virulent variants supports this view (Liu et al. 2014). However, sequence comparisons did not identify a universal virulence determinant among the three virulent lagovirus lineages, suggesting that each virulence acquisition may have involved a unique mutation(s).

Overall, our data suggests that the Lagovirus genus is likely substantially under-sampled and that the true diversity of the hare (and rabbit) virome is underestimated. Additional sequencing of RNA viromes of hares and other lagomorph species is clearly needed to fill in the gaps in the lagovirus phylogeny and those of other RNA viruses, in turn providing broad-scale insights into RNA virus evolution and ecology.

Click here for additional data file.