Full-genome sequencing of German rabbit haemorrhagic disease virus uncovers recombination between RHDV (GI.2) and EBHSV (GII.1).

Rabbit haemorrhagic disease virus (RHDV; genotypes GI.1 and GI.2) and European brown hare syndrome virus (EBHSV; genotype GII.1) are caliciviruses belonging to the genus Lagovirus. These viruses pose a serious threat to wild and domestic rabbit and hare populations around the world. In recent years, an expanding genetic diversity has been described within the genus, with recombination events occurring between the different genotypes. Here, we generated and analysed 56 full-genome sequences of RHDV and EBHSV from rabbit and hare livers, collected in Germany between the years 2013 and 2020. We could show that genotype Gl.2 (RHDV-2) almost entirely replaced Gl.1 (classical RHDV) in the German rabbit population. However, GI.1 is still present in Germany and has to be included into disease control and vaccination strategies. Three recombinant strains were identified from rabbit samples that contain the structural genes of genotype Gl.2 and the non-structural genes of genotype Gl.1b. Of special interest is the finding that sequences from two hare samples showed recombination events between structural genes of RHDV Gl.2 and non-structural genes of EBHSV GII.1, a recombination between different genogroups that has not been described before. These findings lead to the assumption that also a recombination of the non-structural genes of RHDV Gl.2 with the structural genes of EBHSV Gll.1 might be possible and therefore increase the potential genetic variability of lagoviruses immensely. Our findings underline the importance of whole genome analysis with next-generation sequencing technology as one of new tools now available for in-depth studies that allow in depth molecular epidemiology with continuous monitoring of the genetic variability of viruses that would otherwise likely stay undetected if only routine diagnostic assays are used.

1. Introduction

Rabbit haemorrhagic disease (RHD) is a highly contagious and fatal disease and therefore a serious threat to native and captive rabbit populations. The RHD virus (RHDV) and the European brown hare syndrome virus (EBHSV) both belong to the family Caliciviridae of the genus Lagovirus (Vinje et al. 2019). The positive-sense, single stranded RHDV genome is 7.4 kb in length and organized in two slightly overlapping open reading frames (ORFs) (reviewed in Abrantes et al. 2012). ORF1 encodes for seven non-structural proteins (p16, p23, helicase, p29, VPg, protease, and RdRp) and the major structural protein VP1/VP60. ORF2 is expressed from subgenomic mRNAs and encodes the structural protein VP2/VP10 (Vinje et al. 2019).

Based on the 2017 attempt to standardize the Lagovirus nomenclature (Le Pendu et al. 2017), RHDV belongs to genogroup GI, which further divides into genotype Gl.1 (former G1−G6), Gl. 2 (RHDV2/b), Gl.3 (non-pathogenic rabbit calicivirus E-1), and Gl.4 (non-pathogenic rabbit calicivirus A-1/E-2). Genotype Gl.1 further splits into the different variants Gl.1a to Gl.1d. Genogroup GII includes genotype GII.1 (EBHSV) and genotype GII.2 (non-pathogenic hare calicivirus HaCV).

The first outbreak of RHDV Gl.1 was described in China in 1984 within a group of rabbits imported from Germany (Liu et al. 1984). Thereupon, the virus rapidly spread throughout the world. It was detected for the first time in Europe (Italy) in 1986 (Cancellotti and Renzi 1991), and the Americas (Mexico) in 1988 (Gregg et al. 1991). In addition, in 1995 RHDV Gl.1 was deliberately released on Wardang Island, South Australia, to be studied as a biological agent for rabbit control (Cooke 2002). Despite strict quarantine measures, the virus reached the mainland in 1995, where it killed more than 95 per cent of the existing wild rabbit population (Mutze, Cooke, and Alexander 1998). RHDV Gl. 2 was first detected in France at the end of summer 2010 (Le Gall-Recule et al. 2011) and subsequently spread worldwide (Hall et al. 2015; Martin-Alonso et al. 2016; Rouco et al. 2019; World Organisation For Animal Health (OIE), 2019). The fast spread and lethal impact on the local rabbit populations has a detrimental effect on the existing ecosystems through bottom-up trophic effects on predators like the Iberian lynx and the Spanish imperial eagle as published for the Mediterranean ecosystems (Delibes-Mateos et al. 2014; Monterroso et al. 2016).

In most cases, RHDV kills the host within 36 h after onset of fever that occurs after an incubation period of up to 72 h (Abrantes et al. 2012). Affected animals show histopathological lesions in the liver, lungs, spleen, heart, and kidney, combined with an acute hepatitis (Fuchs and Weissenbock 1992). Rabbits younger than 8 weeks are protected against the Gl.1 strains, whereas they are not protected and subsequently die in most cases from an infection with Gl.2 (Mikami et al. 1999; Dalton et al. 2012; Marques et al. 2012, 2014). The host range of RHDV Gl.1 is almost exclusively restricted to European rabbits (Oryctolagus cuniculus), whereas RHDV Gl. 2 is also known to infect different hare species like the Italian hare (Lepus corsicanus), mountain hare (Lepus timidus), Cape hare (Lepus capensis mediterraneus), and European brown hare (Lepus europaeus) (Puggioni et al. 2013; Camarda et al. 2014; Hall et al. 2017; Neimanis et al. 2018).

Rapid developments in the field of next-generation sequencing complement the generation of high-quality whole-genomes in the last couple of years. The growing number of complete RHDV genome sequences enables detailed analyses of the genetic diversity of these viruses. For instance, recombinant virus strains that have been discovered since seem to suggest that recombination drive virus evolution (Abrantes, Esteves, and van der Loo 2008; Forrester et al. 2008; Lopes et al. 2015a; Abrantes et al. 2020a). These were shown to comprise non-structural and structural genes from different genotypes, as it was shown for RHDV Gl.1, and RHDV Gl.2 recombinants, GI.4 and RHDV GI.1 recombinants, GI.3 and RHDV GI.1 recombinants and recombinants of non-pathogenic rabbit calicivirus A-1 and RHDV Gl.2 (Almeida et al. 2015; Lopes et al. 2015a, 2019; Hall et al. 2018; Silvério et al. 2018; Mahar et al. 2018b; Abrantes et al. 2020b). Furthermore, triple recombinants have been described that consist of the structural genes of RHDV Gl.2, non-structural genes of GI.1b, and also a non-pathogenic lagovirus (GI.4) (Silvério et al. 2018).

In Germany neither RHD nor EBHS are notifiable diseases, data about the annual cases are not collected centrally and no reference sequences are available from recently sampled virus strains. Therefore, we investigated the proportion of RHDV and EBHSV cases over the years and sequenced a larger number of whole genomes of different strains collected in Germany between 2013 and 2020 to study the genetic variation and occurrence of recombinant virus strains.

2. Results and discussion

Even though a large number of RHDV cases are regularly reported in domestic and wild rabbits in Germany, no contemporary information was available regarding the genetic diversity of the circulating strains. Therefore, we sequenced the complete genomes of 56 lagovirus isolates from different German federal states collected between 2013 and 2020 (Table 1).

Table 1.

Characterization of the liver samples sequenced in this study and details of the assembly data, including reads sequenced, percentage of viral reads, and size of the read subset used for assembly.

	Sample							Reads
Serial number	ID	Species	Year	Location (state)	Genotype by PCR	Variant based on full-length sequence	Processed as described by	Generated	Viral (%)	Subset used for assembly
1	D102-1.L00435/GER-NW/2013	European rabbit (domestic)	2013	NW	Gl.2		Juozapaitis et al. (2014)	3,407,278	46.9	10,000
2	D102-2.L00436/GER-NW/2013	European rabbit (domestic)	2013	NW	Gl.2		Juozapaitis et al. (2014)	3,888,054	36.2	10,000
3	D105-1.L00437/GER-SN/2013	European rabbit (domestic)	2013	SN	Gl.1	a	Juozapaitis et al. (2014)	3,498,060	37.1	10,000
4	D105-2.L00438/GER-SN/2013	European rabbit (domestic)	2013	SN	Gl.1	a	Juozapaitis et al. (2014)	3,761,010	19.7	30,000
5	D107.L00439/GER-SN/2013	European rabbit (domestic)	2013	SN	Gl.1	a	Juozapaitis et al. (2014)	4,387,384	3.8	20,000
6	D24.L00713/GER-RP/2014	European rabbit (domestic)	2014	RP	Gl.2		Juozapaitis et al. (2014)	2,835,154	14.9	10,000
7	D35.L00714/GER-BY/2014	European rabbit (domestic)	2014	BY	Gl.2		Juozapaitis et al. (2014)	2,032,236	34.9	10,000
8	D44.L00716/GER-NW/2014	European rabbit (wild)	2014	NW	Gl.2		Juozapaitis et al. (2014)	2,295,416	14.3	10,000
9	D46-2.L00718/GER-NW/2014	European rabbit (wild)	2014	NW	Gl.2		Juozapaitis et al. (2014)	2,877,674	13.0	10,000
10	D108-2.L00860/GER-NW/2014	European rabbit (domestic)	2014	NW	Gl.2		Juozapaitis et al. (2014)	2,451,222	46.2	10,000
11	D101-2.L00903/GER-TH/2014	European rabbit (domestic)	2014	TH	Gl.2		Juozapaitis et al. (2014)	3,447,616	91.2	30,000
12	D61-2.L00910/GER-NW/2014	European rabbit (domestic)	2014	NW	Gl.2		Juozapaitis et al. (2014)	1,219,162	47.5	10,000
13	D51-1.L00911/GER-NW/2014	European rabbit (domestic)	2014	SH	Gl.2		Juozapaitis et al. (2014)	783,276	99.7	10,000
14	D51-2.L00912/GER-SH/2014	European rabbit (domestic)	2014	SH	Gl.2		Juozapaitis et al. (2014)	929,050	99.5	10,000
15	D123-3.L00914/GER-NW/2014	European rabbit (domestic)	2014	NW	Gl.2		Juozapaitis et al. (2014)	320,212	99.6	10,000
16	D144-2.L01046/GER-NW/2014	European hare	2014	NW	Gl.2		Juozapaitis et al. (2014)	4,765,356	43.7	10,000
17	EI15-1.L03598/GER-TH/2015	European rabbit (domestic)	2015	TH	Gl.2		Wylezich et al. (2018)	635,357	7.8	25,000
18	D66-15.L03609/GER-BB/2015	European rabbit (domestic)	2015	BB	Gl.1	a	Wylezich et al. (2018)	146,992	9.8	25,000
19	D87-2.L03610/GER-RP/2015	European rabbit (wild)	2015	RP	Gl.2		Wylezich et al. (2018)	257,965	0.6	25,000
20	D166.L03611/GER-SN/2015	European rabbit (domestic)	2015	SN	Gl.1	a	Wylezich et al. (2018)	213,399	4.7	25,000
21	D160-1.L03612/GER-NW/2015	European rabbit (wild)	2015	NW	Gl.2		Wylezich et al. (2018)	260,339	4.3	25,000
22	EI15-3.L03599/GER-TH/2016	European rabbit (domestic)	2016	TH	Gl.1	a	Wylezich et al. (2018)	480,753	16.4	25,000
23	EI15-7.L03600/GER-NW/2016	European rabbit (wild)	2016	NW	Gl.2		Wylezich et al. (2018)	445,830	1.7	25,000
24	EI70-1.L03601/GER-NW/2016	European rabbit (domestic)	2016	NW	Gl.2		Wylezich et al. (2018)	601,237	16.0	25,000
25	EI106.L03602/GER-NW/2016	European rabbit (domestic)	2016	NW	Gl.2		Wylezich et al. (2018)	408,754	0.3	75,000
26	EI145-12.L03603/GER-TH/2016	European rabbit (domestic)	2016	TH	Gl.2		Wylezich et al. (2018)	672,442	3.1	25,000
27	EI176-1.L03604/GER-SN/2016	European rabbit (domestic)	2016	SN	Gl.1	a	Wylezich et al. (2018)	606,929	2.2	25,000
28	EI220.L03605/GER-MV/2016	European rabbit (wild)	2016	MV	Gl.2		Wylezich et al. (2018)	549,250	0.7	25,000
29	EI297-1.L03606/GER-SN/2016	European rabbit (domestic)	2016	SN	Gl.1	a	Wylezich et al. (2018)	482,401	13.3	25,000
30	EI327.L03607/GER-BE/2016	European rabbit (domestic)	2016	BE	Gl.2		Wylezich et al. (2018)	574,845	6.9	25,000
31	EI205.L03564/GER-SN/2017	European rabbit (domestic)	2017	SN	Gl.2		Wylezich et al. (2018)	671,729	4.6	25,000
32	EI67.L03566/GER-SN/2017	European rabbit (domestic)	2017	SN	Gl.1	a	Wylezich et al. (2018)	810,644	0.2	125,000
33	EI175-1.L03570/GER-BE/2017	European rabbit (wild)	2017	BE	Gl.2		Wylezich et al. (2018)	941,128	30.1	25,000
34	EI121.L03597/GER-NW/2017	European rabbit (domestic)	2017	NW	Gl.2		Wylezich et al. (2018)	553,836	15.6	25,000
35	EI44.L03567/GER-BB/2018	European rabbit (domestic)	2018	BB	Gl.2		Wylezich et al. (2018)	786,380	15.8	25,000
36	EI125.L03568/GER-NI/2018	European rabbit (domestic)	2018	NI	Gl.2		Wylezich et al. (2018)	1,045,142	7.7	25,000
37	EI126-2.L03569/GER-HE/2018	European rabbit (domestic)	2018	HE	Gl.2		Wylezich et al. (2018)	772,376	1.2	25,000
38	EI117-2.L03571/GER-BE/2018	European rabbit (domestic)	2018	BE	Gl.2		Wylezich et al. (2018)	848,424	5.2	25,000
39	EI73-1.L03576/GER-RP/2018	European rabbit (domestic)	2018	RP	Gl.1	a	Wylezich et al. (2018)	820,213	24.7	25,000
40	EI73-7.L03578/GER-RP/2018	European rabbit (domestic)	2018	RP	Gl.2		Wylezich et al. (2018)	855,556	23.0	25,000
41	EI11-5.L03608/GER-BW/2018	European rabbit (domestic)	2018	BW	Gl.2		Wylezich et al. (2018)	556,690	5.2	25,000
42	EI06-1.L03565/GER-BY/2019	European rabbit (domestic)	2019	BY	Gl.2		Wylezich et al. (2018)	841,740	5.1	25,000
43	EI48.L03572/GER-BE/2019	European rabbit (domestic)	2019	BE	Gl.2		Wylezich et al. (2018)	895,625	1.5	25,000
44	EI53-2.L03573/GER-BY/2019	European rabbit (wild)	2019	BY	Gl.2		Wylezich et al. (2018)	1,326,632	50.3	25,000
45	EI31.L03574/GER-BE/2019	European rabbit (domestic)	2019	BE	Gl.2		Wylezich et al. (2018)	1,059,942	5.0	25,000
46	EI83-2.L03575/GER-BE/2019	European rabbit (wild)	2019	BE	Gl.2		Wylezich et al. (2018)	832,237	13.2	25,000
47	EI17-1.L03577/GER-NW/2019	European hare	2019	NW	Gl.2		Wylezich et al. (2018)	644,099	0.1	600,000
48	EI16-2.L03579/GER-BB/2019	European rabbit (wild)	2019	BB	Gl.2		Wylezich et al. (2018)	1,181,609	2.3	25,000
49	EI129-10.L03595/GER-NI/2019	European hare	2019	NI	Gl.2		Wylezich et al. (2018)	626,138	3.6	25,000
50	EI112-7.L03613/GER-NW/2019	European hare	2019	NW	Gll.1		Wylezich et al. (2018)	250,988	1.4	25,000
51	EI129-20.L03596/GER-NI/2019	European hare	2019	NI	Gll.1		Wylezich et al. (2018)	505,237	0.7	25,000
52	EI104-12.L03475/GER-BY/2019	European hare	2019	BY	Gll.1		Wylezich et al. (2018)	2,565,964	13.2	25,000
53	EI20-1.L03476/GER-NW/2019	European hare	2019	NW	Gll.1		Wylezich et al. (2018)	2,433,281	8.7	25,000
54	EI97.L03477/GER-BY/2019	European hare	2019	BY	Gll.1		Wylezich et al. (2018)	2,863,849	0.1	250,000
55	EI07.L03593/GER-NI/2020	European rabbit (wild)	2020	NI	Gl.2		Wylezich et al. (2018)	755,600	7.6	25,000
56	EI04-2.L03594/GER-TH/2020	European hare	2020	TH	Gll.1		Wylezich et al. (2018)	455,490	1.4	25,000

Sequences were submitted to the European Nucleotide Archive (ENA) under the study number PRJEB38327. Abbreviation for German states: BB, Brandenburg; BE, Berlin; BW, Baden-Württemberg; BY, Bavaria; HE, Hesse; MV, Mecklenburg-Western-Pomerania; NI, Lower Saxony; NW, North Rhine-Westphalia; RP, Rhineland-Palatinate; SN, Saxony; SH, Schleswig-Holstein; TH, Thuringia.

The data on reported Lagovirus cases shown in Fig. 1 were collected by the German Consultant Laboratory for RHDV at the Friedrich-Loeffler-Institut (FLI) between the years 2014 and 2019. Genotypes were determined using differentiating diagnostic RT-qPCR tests. Due to voluntary data acquisition by the German federal states, the set-up is not complete and does not fully represent the situation in Germany. It does, however, reflect very well a trend in the number of cases between the years and predominance of certain genotypes. An overview of Germany including the federal state borders is shown in Supplementary Fig. S1. The total number of reported RHDV (Gl.1 and Gl.2) and EBHSV (Gll.1) cases differs significant by year, with a peak of 425 total cases in 2016 (Fig. 1). In every year, Gl.2 was by far the most common genotype identified (63.3 − 99.4% of all cases), with only a few cases accounting for Gl.1 and even fewer for Gll.1 (Fig. 1). Although the survey data underestimate the overall amount of cases due to in the lack of data submission by certain German federal states, the ratio of different genotypes does suggest the replacement of Gl.1 by Gl.2 in Germany, a situation already described before in France (Le Gall-Recule et al. 2017), Portugal (Lopes et al. 2014), and Australia (Mahar et al. 2018a).

Figure 1.

Number and proportion of reported Gl.1, Gl.2, and Gll.1 cases in Germany between 2014 and 2019. Data were collected by the German Consultant Laboratory for RHDV based on both own test results and voluntary submissions from several federal states, as Lagovirus cases are not reportable in Germany. The figures provide a representative section of the epidemic situation in Germany.

Phylogenetic trees based on the structural- (nucleotide 5296 − 7369; Fig. 2) and non-structural genes (nucleotide 1 − 5295; Fig. 3), respectively, were generated to cluster the sequenced isolates according to their genotypes. Nucleotide numbering is based on the reference sequence NC_001543 (RHDV-FRG, Germany 2000). Based on the tree, ten sequences form a well-supported cluster with genogroup Gl.1a, 35 sequences cluster with genogroup Gl. 2, and 6 sequences cluster with Gll.1, respectively. These finding are in complete accordance with the results of the differentiating diagnostic RT-qPCR tests. Within these genogroups, the sequences show high nucleotide and amino acid identities (Table 2). Therefore, the hypothesis that the steep increase in reported Gl. 2 cases in 2016 was caused by substantial genetic changes of the circulating virus strains could not be confirmed.

Figure 2.

Maximum Likelihood phylogenetic tree based on 205 sequences of the structural genes (nucleotide 5296 − 7369), using the nucleotide substitution model SYM+R4. Sequenced strains fall into genotype Gl.1a, Gl.2, and Gll.1. Horizontal scale bar is proportional to the number of nucleotide substitutions per site. Representative publicly available sequences of the different genotypes Gl.1, Gl.2, Gl.4, Gll.1, and Gll.2 are included and depicted in black. Genotype clusters, which do not include virus variants sequenced in this study, were collapsed and annotated accordingly. Recombinants D51-1.L00911/GER-NW/2014, D51-2.L00912/GER-SH/2014 and EI44.L03567/GER-BB/2018 are depicted in orange, D144-2.L01046/GER-NW/2014 and EI17-1.L03577/GER-NW/2019 in red. All other sequences generated in this study are depicted in blue. Branches with a bootstrap value ≥80 are coloured green. The EBHSV reference strains previously isolated from hares are marked with an asterisk (*).

Figure 3.

Maximum Likelihood phylogenetic tree based on 205 sequences of the non-structural genes (nucleotide 1 − 5295), using the nucleotide substitution model SYM+R5. Sequenced strains fall into genotype Gl.1a, Gl.1b, Gl.2, and Gll. Horizontal scale bar is proportional to the number of nucleotide substitutions per site. Publicly available sequences of the different genotypes Gl.1, Gl.2, Gl.4, Gll.1, and Gll.2 are included. Genotype clusters, which do not include virus variants sequenced in this study, were collapsed and annotated accordingly. Recombinants D51-1.L00911/GER-NW/2014, D51-2.L00912/GER-SH/2014, and EI44.L03567/GER-BB/2018 are depicted in orange, D144-2.L01046/GER-NW/2014and EI17-1.L03577/GER-NW/2019 in red. All other sequences generated in this study are depicted in blue. Branches with a bootstrap value ≥80 are coloured green. The EBHSV reference strains previously isolated from hares are marked with an asterisk (*).

Table 2.

Nucleotide- and amino acid identities of the analysed full-length genomes within the genogroups Gl.1, Gl.2, Gll.1, and the recombinant strains GI.1b−GI.2 and GII.1−GI.2.

Genogroups	Nucleotide identities (%)	Amino acid identities (%)
Gl.1a	95.4 − 99.6	98.0 − 99.9
Gl.2	91.9 − 99.9	96.4 − 100
Gll.1	96.3 − 99.8	98.6 − 99.8
GI.1b-GI.2	96.4 − 100	97.3 − 100
GII.1−GI.2	86.9	95.6

Additional five samples cluster with genogroup Gl.2 based on the structural genes (Fig. 2), but cluster clearly with a different genogroup based on the non-structural genes (Fig. 3). This is a strong evidence for recombination and has been described in lagoviruses before (Lopes, Hubatsch, and Amaris 2015; Hall et al. 2018; Mahar et al. 2018b; Lopes et al. 2019).

Three of the five recombinant strains (D51-1.L00911/GER-NW/2014, D51-2.L00912/GER-SH/2014, and lib03567/GER-BB/2018) cluster with Gl.2 based on the structural genes (Fig. 2), but cluster with Gl.1b and previously identified GI.1b−GI.2 recombinants based on the non-structural genes (Fig. 3). This was also indicated by recombination analyses in recombinant detection program (RDP) (using RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan, and 3Seq), which identifies all three recombinants with strong statistical support (average P-values <0.05) (Table 3). The point of recombination in these chimeric viruses is the junction between non-structural RNA-dependent RNA polymerase and structural VP1 (breakpoint at nucleotide position 5294, Table 3). The three samples were collected in two non-bordering states in 2014 and 2018 and we therefore assume that these recombinants are circulating in Germany until today. Based on the clusters formed, these three recombinants seem to be genetically most closely related to the previously described GI.1b−GI.2 recombinants from Portugal (Lopes et al. 2015a; Silvério et al. 2018). Therefore, we would hypothesize that this GI.1b−GI.2 recombinant emerged from a distinct area, e.g. the Iberian Peninsula, where it has been first described in 2014 (Lopes et al. 2015a) and subsequently spread into Germany. This would be in line with the findings of Lopes et al. (2019), who described GI.1b−GI.2 recombinants found in Africa that resemble recombinants from the Iberian Peninsula. The nucleotide identity to the potential parents is high (96.1%) for the structural genes but relatively low for the non-structural genes (77.9 − 78.1%) (Fig. 4).

Table 3.

Results of the recombination analyses using recombinant detection program (RDP).

Sample	Most likely parental lineage		Breakpoint (99% confidence interval)	Average P value
	Non-structural	Structural		RDP	GENECONV	BootScan	MaxChi	Chimaera	SIScan	3Seq
D51-1.L00911/GER-NW/2014; D51-2.L00912/GER-SH/2014; EI44.L03567/GER-BB/2018	Gl.1b (JX886001)	Gl.2 (KF442964)	5294	1.1×10−38	1.2×10−103	1.8×10−2	1.7×10−18	6.4×10−21	6.2×10−26	8.0×10−26
D144-2.L01046/GER-NW/2014; EI17-1.L03577/GER-NW/2019	Gll.1 (EI104-12.L03475/ GER-BY/2019)	Gl.2 (EI83-2.L03575/ GER-BE/2019)	5297 − 5309	1.4×10−102	2.4×10−168	1.2×10−2	9.7×10−34	1.0×10−37	4.2×10−63	3.1×10−70

To verify potential recombination events, the alignments were screened with a combination of seven different analytical methods. Only recombination events that were detected by three or more methods with highest acceptable P value of 0.05 were considered in this analysis.

Figure 4.

SimPlot results of the recombinant GII.1−GI.2 strains: (a) D144-2.L01046/GER-NW/2014 and (b) EI17-1.L03577/GER-NW/2019; blue line: Gl.2 (EI83-2.L03575/GER-BE/2019), red line: Gll.1 (EI104-12.L03475/GER-BY/2019), grey line: control strain Gl.1a (EF558583). Nucleotide position is depicted on the x-axis (non-structural genes: nucleotide 1 − 5285/95, structural genes: nucleotide 5300/−305 to 7369), similarity in percentage on the y-axis. The following settings were used: window 200 bp, step size 20 bp, GapStrip: On, Kimura distance model. (c) Diagram of the RHDV genome with open reading frames (ORF), non-structural genes (p16, P23, helicase, P29, VPg, protease, and RdRp), and structural genes (VP1 and VP2).

The two other recombinant strains that were isolated from deceased hares (D144-2.L01046/GER-NW/2014 and EI17-1.L03577/GER-NW/2019) stand out in the phylogenetic tree since they cluster with the genogroup Gl.2 based on the structural genes (Fig. 2) but form an isolated cluster in the genogroup GII based on the non-structural genes (Fig. 3). Only nine individual GII.1 whole genomes were available (NCBI database (NCBI 1988; accessed March 2020), March 2020) by the time of our analyses. To allow comparison to recently circulating strains, we sequenced the full genome of a number of Gll.1 strains (n = 6) from Germany, sampled in 2019 and 2020. Nevertheless, these sequences cluster with the already existing Gll.1 strains and do not bridge the gap to D144-2.L01046/GER-NW/2014 and EI17-1.L03577/GER-NW/2019 (Fig. 3). However, this is not surprising since the two recombinants differ by 17 per cent on nucleotide level in the scope of the non-structural genes between each other and therefore we expect a very high diversity within the Gll.1 genotype. Moreover, the high sequence diversity between the non-structural genes of the two recombinants suggests that they arose from two independent recombination events, which is further supported by the divergent position of the two sequences in the structural phylogeny (Fig. 2).

Recombination analysis with RDP supports that both parents of the recombinants are the closest related to German sequences obtained in this study. The closest major parent of both recombinants is a GII.1 strain (EI104-12.L03475/GER-BY/2019) and the closest minor parent is a GI.2 strain (EI83-2.L03575/GER-BE/2019) (Table 3). The nucleotide identity to the potential parents is high (96.1%) for the structural genes but relatively low for the non-structural genes (77.9 − 78.1%). Based on these results, we assume that the two recombinants did emerge in Germany. Nevertheless, a lack of reference sequences (in particular for GII.1) makes a conclusive interpretation of these findings difficult.

The recombination breakpoint of the GII.1−GI.2 recombinants is located at the junction of non-structural and structural genes (Table 3) and therefore in line with other recombination events that were previously found in RHDV recombinants from the Iberian Peninsula, Australia, and Africa (Almeida et al. 2015; Hall et al. 2015; Lopes et al. 2015a, 2019; Mahar et al. 2018a). Recombinants D144-2.L01046/GER-NW/2014 and EI17-1.L03577/GER-NW/2019 were sampled in the same state (100 km air-line distance) five years apart (2014 and 2019), which implies that these virus strains may be stably circulating. Both recombinants D144-2.L01046/GER-NW/2014 and EI17-1.L03577/GER-NW/2019 were isolated from hare samples which is consistent with the finding that the host range of EBHSV and RHDV Gl.2 overlaps in hares (Le Gall-Recule et al. 2017). However, the majority of the hare isolates form a closely related cluster of strains from occasional cases collected early in 2019 and from die-off events towards autumn 2019. In these cases, the molecular analyses of the virus genomes did not substantiate an emergence of EBHSV with enhanced virulence.

The emergence of new virus variants has been widely reported (Abrantes, Esteves, and van der Loo 2008; Almeida et al. 2015; Lopes et al. 2015a, 2019; Hall et al. 2018; Mahar et al. 2018b; Abrantes et al. 2020a) but recombinants between lagovirus genogroups were undetected until now. This shows the importance of whole genome screening to detect the emergence of new virus variants. At this point it might be worth noticing that the finding of the recombination between the structural genes of Gl.2 and the non-structural genes of Gll.1 leads to the assumption that also a recombination of the non-structural genes of Gl.2 with the structural genes of Gll.1 might be possible.

In summary, we describe in this study the replacement of RHDV genotype Gl.1 by Gl.2 in Germany. We provided 56 lagovirus genomes from samples collected between 2013 and 2020, of which five were identified to be recombinants. The non-recombinant isolates from 2013 to 2020 cluster rather homogenously and therefore only the recombinants seem to contribute significantly to the genetic diversity of RHDV in Germany. Among these five recombinants, two isolates stand out in particular, since they appear to represent recombinants of EBHSV GII.1 and RHDV Gl.2. Further research should provide a variety of full-genome sequences of Gll.1 isolates to elucidate the genetic variability and identify possible additional recombination events. Specific real time RT-PCR systems have been developed to facilitate future detection of recombinant Lagoviruses in the course of routine diagnostics. Closer characterization of the recombinants in the respective Lagomorph host will contribute to a better understanding of host spectrum and pathogenicity of the viruses and their possible impact on rabbit and hare populations.

3. Materials and methods

3.1 Samples

All samples analysed in this study were selected from deceased animals submitted to the FLI for confirmatory diagnostics with the suspicion of an RHDV or EBHSV infection. In total, forty-seven rabbit (domestic and wild) and nine hare liver samples were sequenced in this study to generate whole genome sequences. Samples were chosen based on the date of collection (between 2013 and 2020) and the geographic location (covering different regions in Germany). Samples were pre-tested and genotyped (Gl.1 (n = 10); Gl.2 (n = 40); Gll.1 (n = 6)) based on a combination of differentiating RT-qPCRs (Gall et al. 2007; Schirrmeier and König 2019) and the INgezim RHDV ELISA (Eurofins Ingenasa, Madrid, Spain) for the detection of Gl.1, Gl.2, and Gll.1 antigens.

3.2 Nucleic acid extraction, sequencing, and sequence assembly

For RNA extraction, liver samples were homogenized in 1 ml minimum essential medium with the TissueLyser II instrument (QIAGEN, Hilden, Germany) using 5 mm steal beads. Subsequently, RNA was extracted with the QIAamp viral RNA Mini Kit (Qiagen) according to the manufacturer’s instructions.

Samples collected in 2013 and 2014 were processed as described by Juozapaitis et al. (2014) and sequenced on the Illumina platform. Samples collected after 2014 were processed using a modification of the protocol described by Wylezich et al. (2018), and sequenced employing the Ion Torrent S5XL platform. Briefly, extracted RNA was concentrated using the Agencourt RNAClean XP magnetic beads (Beckman Coulter, Fullerton, CA, USA). Subsequently, cDNA was synthesized using a combination of the SuperScript™ IV First-Strand cDNA Synthesis System (Thermo Fisher Scientific, Waltham, MA, USA) and the NEBNext®Ultra™ II Non-Directional RNA Second Strand Synthesis Module (New England Biolabs, Ipswich, MA, USA). The generated cDNA was fragmented to 500 bp using the Covaris M220 Focused-Ultrasonicator (Covaris, Brighton, UK) before Ion Torrent-compatible libraries were generated using the GeneRead L Core Kit (Qiagen) and IonXpress Barcode Adapter (Thermo Fisher Scientific). After size selection, library quality was checked with the Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA) and the library concentration determined with the KAPA Library Quantification kit (Roche, Mannheim, Germany). Using an Ion Torrent S5 XL, libraries were sequenced on an Ion 530 chip in 400-bp mode according to the manufacturer’s protocol (Thermo Fisher Scientific).

For sequence assembly, a random subset of reads was assembled (Newbler v.3.0; 454/Roche). In case fragmented genome sequences were obtained, the contigs were combined into full-genome sequences supported by information from file 454ContigGraph. In case of incomplete 5′ or 3′ termini, additional reads representing the genome termini were identified in the complete datasets by searching a short sequence (50 nucleotides) of the terminal end using BLASTn (Altschul and Lipman 1990). Subsequently, the initial dataset—supplemented with the additional reads—was assembled to yield the whole genome sequence. The complete data sets were finally mapped against the assembled genomes (Newbler v.3.0) to confirm the sequences.

3.3 Phylogenetic analyses

Phylogenetic trees were inferred based on the coding region of the non-structural gene sequences (nucleotide 1 − 5295) and the structural genes (nucleotide 5296 − 7369) of the 56 whole genomes generated in this study, and publicly available genomes of the Lagovirus genotypes Gl.1, Gl.2, GI.3, Gl.4, Gll.1, and Gll.2 (Supplementary Table S1). The final dataset included 205 sequences. The best fitting nucleotide substitution model was determined by ModelFinder (Kalyaanamoorthy et al. 2017), based on the lowest corrected Akaike information criterion value. Using IQtree (Nguyen et al. 2015) and ultrafast bootstrap (Hoang et al. 2018), the maximum likelihood consensus tree was constructed from 100,000 bootstrap trees.

3.4 Recombination analysis

Recombination events were analysed as described by Lopes et al. (2017). To this end, the complete sequences generated in this study were aligned with published sequences from the genotypes Gl.1, Gl.2, GI.3, Gl.4, Gll.1, and Gll.2 using the MAFFT Alignment with the GeneiousPrime software (version 2019.2.3). The final dataset included 205 sequences and the alignment was screened for potential recombination using seven methods (RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan, and 3Seq) implemented in the RDP (v.3.44) (Martin et al. 2010). The following parameters were used: set to linear, Bonferroni correction and 100 permutations. Only recombinations detected by three or more methods with a highest acceptable P value of 0.05 were considered in this analysis. The sequence similarity plotting software SimPlot (v.3.5.1) (Lole et al. 1999) was used for depicting the recombinant strains.

Supplementary data

Supplementary data are available at Virus Evolution online.

Click here for additional data file.