Genetic variability and phylogenetic analysis of Lagovirus europaeus strains GI.1 (RHDV) and GI.2 (RHDV2) based on the RNA-dependent RNA polymerase (RdRp) coding gene

Lagovirus europaeus GI.1 (RHDV-rabbit haemorrhagic disease virus) and GI.2 (RHDV2-rabbit haemorrhagic disease virus 2), family Caliciviridae, genus Lagovirus, are etiological factors of the rabbit haemorrhagic disease (RHD). This small RNA virus is a great model for tracking the variability and evolution of RNA viruses, because it uses an RNA-dependent RNA polymerase (RdRp) to replicate its own genetic material. This polymerase determines the fidelity and the rates of replication and mutation of the virus, conditioning its adaptation to the environment and even to a new host, and thus influencing evolution of the virus. The aim of this study was to determine the genetic variability and phylogenetic relationships of 105 Lagovirus europaeus strains with different genotypes based on the RdRp gene. The strains came from around the world in the years of 1987–2017. The aforementioned group of 105 strains included 14 strains whose RdRp sequences were obtained and analysed in this study, and the rest were retrieved from GenBank: 74 strains classified as genotype GI.1 (RHDV), 14 as GI.2 (RHDV2), 2 strains of Lagovirus europaeus not assigned to any genotype, and a MRCV strain, the sequences of which were collected from GenBank. Among the 14 strains whose RdRp sequences were obtained in this study, the highest variability was presented in the Austrian 237 strain from 2004. The genetic distance between the Austrian 237 strain and the remaining thirteen analysed strains ranged from 0.117 to 0.123 (from 11.7% to 12.3% nucleotide substitutions). The lowest variability, however, was recorded for Hungarian, Czech and Austrian strains. On the phylogenetic tree, the 14 analysed strains were allocated into GI.1c (G2), GI.1d (G3-G5) and GI.1a (RHDVa). Analysis of the genetic variability of the 105 strains of Lagovirus europaeus indicated a growing genetic distance between the strains, both in time and location. Phylogenetic analysis showed a division of the strains into seven groups, dictated by the chronology, geographical location and evolutionary events in the history of the virus, such as mutations and recombinations.


INTRODUCTION
Lagovirus europaeus genotype GI.1 (RHDV) and GI.2 (RHDV2/b) are etiologic factors of Rabbit Haemorrhagic Disease (RHD), belonging to the Caliciviridae family, Lagovirus genus. The disease was recorded for the first time in 1984 in China, in the Jiangsu Province, in rabbits imported from Germany (Liu et al., 1984). In less than a year, the disease had spread across China, causing rabbit deaths on a massive scale (Liu et al., 1984). In Europe, the disease was first reported in Italy in 1986; with further occurrences recorded in Czechoslovakia in 1987, in Germany in 1988, in France and Spain in 1989, and at the turn of 1988 in Poland (Hukowska-Szematowicz et al., 2013). The RHDV strains isolated in the earliest period of RHDV occurrence in Europe (1987)(1988)(1989) were characterised with high homology and low genetic variability among one another, and are referred to as the reference strains, and include the Czech V351 strain from 1987, as well as strains from 1989: German FRG, Spanish AST89, French SD, and Italian BS89. In 1996, the first case of RHDV was recorded (Capucci et al., 1998), which was different from classic RHDV in the aspect of antigenicity, genetic variability and phylogenetic analysis of the gene encoding structural protein VP60, referred to as RHDVa (Capucci et al., 1998;Kerr et al., 2009;Kinnear & Linde, 2010;Fitzner & Niedbalski, 2017). It was also evidenced that RHDVa strains are characterised with higher virulence and mortality rates in rabbits (Capucci et al., 1998). Continuous monitoring points out that in some areas, RHDVa seems to replace RHDV, as is currently being observed in Australia (Mahar et al., 2018b). Phylogenetic surveys have classified the RHDV strains into six genogroups (G1-G6) (Le Gall-Recule et al., 2003), while recently, they have been allocated into four variants of Lagovirus europaeus GI.1 (RHDV): GI.1a (previously G6/RHDVa), GI.1b (previously G1), GI.1c (previously G2), GI.1d (previously G3-G5), according to a newly proposed nomenclature (Le Pendu et al., 2017).
It was also evidenced that RHDV2 causes diseases in other lagomorph species, including the Sardinian Cape hare (Lepus capensis mediterraneus) (Puggioni et al., 2013), Corsican (Italian) hare (Lepus corsicanus) (Camarda et al., 2014), Iberian hare (Lepus granatensis) (Lopes et al., 2014) and European hare (Lepus europaeus) (Hall et al., 2017), and recently even in the mountain hare (Lepus timidus) (Neimanis et al., 2018) and hares in England (Bell et al., 2019), unlike RHDV and RHDVa, which are believed to be characteristic for the European rabbit (Oryctolagus cuniculus). The latest case of infection is the sixth Lepus species to be infected with RHDV2, which suggests that the susceptibility to RHDV2 can be widespread among the Lepus species (Neimanis et al., 2018;Bell et al., 2019). Holms and Grenfell (Holms & Grenfell, 2009) explain that the occurrence of new variants of the virus in the wild is the result of four factors: (1) adaptive evolution of the host's genes that engage in the most thorough interaction with the host's immune response, and react at the "host-pathogen" level, (2) interaction among the locally circulating viruses, (3) possible outbreak dynamics in time and space, and (4) disease spread control methods. Hypothetically, it can be assumed that in the case of RHDV, we are dealing with all four factors. Since the isolation of RHDV more than 30 years ago, the genetic variability and evolution of RHDV have shown a specificity in the RdRp enzyme used by the virus to replicate its own genetic material and to express its genome. The enzyme has no corrective ability to repair the many errors generated during the replication process (Chen et al., 2003). RdRp determines the fidelity, as well as the speed of replication and mutation, and conditions its adaptation to the environment and even to a new host (Chen et al., 2003), thus affecting the evolution of the virus (Smertina et al., 2019). It is also for the above reasons that RdRp may prove to be a more appropriate molecular target to track variability and phylogeny than the capsid protein (Koonim, 1991), as for example, the new RHDV2 genotypes reveal an increasing variability in their genome, which is due to recombination at the edge of structural and non-structural fragments, including RdRp . To date, there has been no research carried out on the variability and phylogeny of RHDV (including RHDV, RHDVa, and RHDV2) based on RdRp. With the emergence of new RHDVa strains in Australia (Mahar et al., 2018b) and RHDV2 in the Macronesian Islands  and Poland (Fitzner & Niedbalski, 2018), and due to the fact that RHDV2 has extended its infectious spectrum to hares, it seems crucial to monitor the variability of this virus for epidemiological reasons.
This study aimed at assessing genetic variability and performing a phylogenetic analysis of a total of 105 strains of Lagovirus europaeus, with RHDV and RHDV2 genotypes originating from around the world in the period of 1987-2017, pursuant to the nucleotide sequence of the RdRp coding gene. The aforementioned group included 14 strains whose RdRp sequences were obtained and analysed in this study (deposited by the author in GenBank), and other strains whose RdRp sequences were retrieved from GenBank: 74 strains classified as genotype GI.1 (RHDV), 14 as GI.2 (RHDV2), 2 strains of Lagovirus europaeus not assigned to any genotype, and a MRCV strain, the sequences of which were collected from GenBank.

MATERIALS AND METHODS
Lagovirus europaeus GI.1 (RHDV) strains. The research included 14 strains of Lagovirus europaeus, which according to the newly proposed classification and nomenclature system for Lagoviruses, are classified in group GI.1 (RHDV) ( Table 1). All the strains originated from research centres in Europe, namely the Central Veterinary Institute, Institute of Debrecen, Hungary; the Collection of Animal Pathogenic Microorganisms, Veterinary Research Institute, Brno, Czech Republic; Friedrich-Loeffler Institute, Institut für Virusdiagnostik, Greifswald, Germany. RHDV strains used for the genetic studies were prepared in the form of a suspension in glycerol.
Isolation of viral RNA. Total RNA from fourteen RHDV strains were isolated from the suspension in glycerol using a High Pure Viral RNA Kit (Roche Diagnostic, Switzerland) according to the manufacturer's recommendations. RNA obtained in this manner was directly used for a reverse transcription reaction, or stored at -80 o C for further analysis.
Reverse transcription (RT) reaction -cDNA synthesis, PCR reaction and primers. A complementary strand of the nucleic acid (cDNA) of the gene encoding RdRp from 14 RHDV strains was obtained in the reverse transcription reaction on a matrix of viral RNA using the Reverse Transcriptase enzyme (RevertAid H Minus M-MuLV Reverse Transcriptase) (Fermentas, Lithuania). The reaction was performed using a reverse transcription kit -First Stand cDNA Synthesis Kit (Fermentas, Lithuania) with the following components: Random Hexamer primer at a concentration of 100 μM, dNTPs at10 mM, RevertAid H Minus M-MuLV Reverse Transcriptase-200U/μl, 5xReaction Buffer, RiboLock RNase inhibitor and water for molecular biology. For the reaction, RNA of the relevant strains (24V,118V, 94VM, V411, 1600VM, 1042, V562, V560, 894V, 948V, 1447V, Wika, 72V, 237) ( Table 1) was used. Before the reaction mixture was prepared, RNA of the 14 analysed RHDV strains was heated for 5 minutes at 65°C, and then stored on ice until the mixture was prepared. The RT reaction was performed in a total volume of 20 μl and included: 1.0 μl of Random Hexamer primer, 6.0 μl of water for molecular biology, 4.0 μl RT buffer, 1.0 μl RiboLock Rnase inhibitor, 2.0 μl dNTPs, 1.0 μl Re-vertAid H Minus M-MuLV Reverse Transcriptase, and 5.0 μl RNA of the relevant strain (Table 1). The RT reaction was conducted in a T-gradient Thermocycler (Biometria, Germany) using the following temperature-time profile: 25ºC -5 minutes, 42ºC -60 minutes, 70ºC -5 minutes and 4ºC -1 minute. cDNA samples were stored at 2-8°C until further analyses. The PCR reaction was prepared in a 25 μl mixture containing: 2.0 μl of the primer pair (1.0 μl each primer: RHDV _RTFor and RHDV _RTRev) (Sigma Aldrich, Germany), 1.0 μl dNTPs (Fermentas, Lithuania), 2.5 μl Taq Plus buffer (Genoplast, Poland), 0.5 μl Taq plus DNA polymerase (Genoplast, Poland), 18.0 μl of water for molecular biology (Eppendorf, Germany), and 1.0 μl cDNA of the relevant RHDV strain (Table 1) added as the last component in the reaction mixture. The PCR reaction was performed in a Mastercycler proS Thermocycler (Eppendorf, Germany), applying the following amplification conditions: preliminary denaturation -94°C -3 min, and then 35 cycles comprised consecutively of: denaturation (94°C -30 s), primer annealing (63°C -1 min), chain elongation (72°C -2 min); final elongation (72°C -8 min), and cooling of the reaction mixture to 4°C. Reaction products were stored at 4°C until further analyses. Primers were designed on the basis of comparison of the nucleotide sequence encoding RdRp of the FRG RHDV strain (M67473) with the sequence of RHDV homologues. For amplification of the RdRp coding sequences (positioned in the genome at 3750-5350), the following primers were used: forward RHDV_RT-For 5'cactggcaagytgttggggttttc3', reverse RHDV_RTRev 5'gttccgggaactgatgctgtgg3'. A gene of length 1600 bp was amplified. Primer synthesis was performed by Sigma-Aldrich (Germany).
Electrophoresis of PCR products, preparative amplification and preparation of the RHDV genes for sequencing. For the purpose of visualisation of the PCR products, electrophoresis was performed in 2% agarose gels (Prona, USA) dyed with ethidium bromide (Fermentas, Lithuania). Molecular weight marker GeneRuler 200 bp (Fermentas, Lithuania) was used for evaluation of the size of the products. Interpretation of the electrophoretic analysis results was completed using a UV visualisation kit (UVP, Germany). After visualisation of the PCR reaction products, a preparative PCR reaction was performed to obtain the appropriate amount of products for sequencing. Conditions and the course of the reaction were identical to the conditions described above. After visualisation of the PCR products on a gel, they were isolated (from the gel) using the Gel-Out kit (A&A Biotechnology, Poland), according to the manufacturer's protocol. The PCR products were then sent to Genomed, Poland, for sequencing.
Comparative analysis of nucleotide sequences of the gene encoding RdRp, creation of the distance matrices and the phylogenetic tree. The obtained 1284 nucleotide-long sequences of the gene encoding RdRp from the 14 RHDV strains were deposited in GenBank (Table 1). Comparative analysis of these nucleotide sequences involved a total of 105 Lagovirus europaeus strains (Tables 1 and 2), as generated by the MEGA 5.0 software, permitting genetic variability to be determined by generation of a distance matrix. Analyses were conducted using the Maximum Composite Likelihood model (Tamura et al., 2004) with the MEGA 5.0 software (Kumar et al., 2018). Variability within the nucleotide and amino acid sequences of RdRp was determined against a reference sequence of the V351 strain (U54983). Supplementary Table S1 and Table S2 (at https://ojs.ptbioch.edu.pl/index. php/abp/) show this variability also with respect to other reference strains from 1989: FRG (M67473), AST89 (Z49271), SD (Z295514), BS89 (X87607). The distance matrix shows the genetic distances for all sequence pairs within the range of 0.000 to 1.000. Genetic distance (manifesting genetic variability) between the sequences of a pair of strains corresponds to the number of nucleotide substitutions (expressed as a percentage) which occurred in such sequences since their divergence from the output sequence. For a clear presentation of the test results, genetic variability is expressed as percentage. Distance matrix tables (Supplementary Material - Table S1 at https://ojs.ptbioch.edu.pl/index.php/abp/) show the values of nucleotide substitutions, which correspond to the percentage of genetic variation. Due to the long strain names in the text, only the first part of the strain name was given. The sequences have been reported to GenBank and accepted, received accession numbers and will be made available in a public database at the time of publication.

P95_PortugalToresNovas1996
Lepus granatensis KJ943791 Due to the long strain names in the text, only the first part of the strain name was given.
The evolutionary history was inferred using the Maximum Likelihood method and General Time Reversible model (Nei & Kumar, 2000), the most statistically appropriate method for the selected data. A discrete Gamma distribution was used to model the evolutionary rate differences among sites. Evolutionary analyses were conducted with MEGA 5.0 (Kumar et al., 2018). Reliability of phylogenetic trees were assessed using the bootstrap method (%). Distance matrices were also developed (using the method stipulated above) for particular groups of the strains generated in the phylogenetic tree in order to determine the intra-group variability. In the last phase, the nucleotide sequences of the 105 RHDV strains were translated to amino acid sequences, and a distance matrix was generated to determine their variability using the MEGA 5.0 software (Kumar et al., 2018).

Genetic variability of Lagovirus europaeus strains
Nucleotide sequences of genes coding RdRp, obtained in this study from 14 Lagovirus europaeus strains, with the GI.1 genotype (RHDV), (Table 1), had a length of 1284 nucleotides and covered positions 3955-5238 in the viral genome, while the full RdRp coding sequence has a length of 1548 nucleotides and occupies positions 3756-5301. Genetic variation of the 14 analysed strains within RdRp indicated occurrence of many polymorphic sites with local mutations in the form of transitions (C↔T, G↔A) and transversions (G↔T, A↔T, C↔G). Transitions were more frequent, comprising over 90% of mutations, involving frequent C↔T changes (over 50%) versus G↔A. Among the 10% of transversions, all mutations occurred with a comparable frequency.

Phylogenetic analysis of the Lagovirus europaeus strains
The image of phylogenetic interdependencies (Fig. 1) among 105 Lagovirus europaeus strains, build on the basis of RdRp, has revealed a clear division of these strains into 7 groups, where the 14 strains analysed in this study were allocated into GI.1c (G2), GI.1d (G3-G5), and GI.1a (RHDVa) groups. Distance matrices generated for groups of Lagovirus europaeus strains (Table S3-S9 at https://ojs.ptbioch.edu.pl/index.php/abp/), isolated pursuant to the topology of the phylogenetic tree (Fig. 1), revealed that variability in RdRp within the groups (intra-group variability) is much lower than variability obtained for all 105 Lagovirus europaeus strains. Lower variability was caused by much lower numbers of nucleotide substitutions in such groups.
Group GI.1b (G1) included 11 strains from the period 1989-2017, namely the last 28 years, and these were the strains isolated in the first period of the disease occurrence in Europe, as well as contemporary strains from continental Portugal and the Portuguese islands. Evolutionary relations in this group clearly point to a chronological distance that has divided the strains into two subgroups. The first has been formed by RHDV from 1989-SD, AST89, Eisenhuettenstadt, with CB137_ Pt from 1995 and CB194_Pt from 2006. While the latter by Portuguese strains from the period of 2014-2017: PSM2, CBMad17-3, CBMad17-2, CBPico17-2, CBPi-co17-1, CBAlgarve14-3. The distance matrix generated to determine intragroup genetic variability (Table S3 at https://ojs.ptbioch.edu.pl/index.php/abp/) revealed variability ranging between 0.2-9.7%. The shortest genetic distance in this group was recorded, mutually, for strains CBMad17-3 from 2017 and PSM2 from 2016 (0% nucleotide substitutions), and then for CBPico 17-1 and CB-Pico17-2 from 2017, and for CBMad 17-2 and CBMad 17-3 from 2017 (in both cases 0.2% nucleotide substitutions), as well as AST89 and Eisenhuettenstadt (0.9%). Portuguese strain CB137Pt showed a genetic distance to strains in its subgroup between 3.1-3.5%, and to other Portuguese strains between 6.8-9.7%.
Group GI.1c (G2) was comprised of 35 strains (including nine strains analysed in this study) from the period of 1987-2014, originating from around the world, including Australia. Evolutionary relations in this group divided the strains into three subgroups. The first subgroup included European strains phylogenetically deriving from the Czech V351 strain from 1987, which included nine strains analysed in this study, as well as the following strains: Italy 90, PD, KGM, MAL, FRG. The second subgroup included strains from New Zealand and Australia. The third subgroup included HB, Korea-90, Mexico89, and SaudiaArabia (Fig. 1). Intragroup variability (Table S4 at https://ojs.ptbioch.edu.pl/ index.php/abp/) ranged from 0% to 8.8%. The shortest genetic distance in this group was recorded for strains NZ54 and NZ61 (0%), 72V and 94VM (0%), V351 and FRG (0.9%). The most genetically distanced from other strains in this group were the four Australian strains from 2013 and 2014, which revealed a genetic distance of 5.7-8.8%.
The image of phylogenetic relations based on RdRp clearly pointed to a new monophyletic group situated between GI.1a-1d and RHDV2, with a high statistical support (bootstrap value 99%). It was created by two strains of the Portuguese origin: P95_PortugalToresNo-vas1996 (host: Lepus granatensis) and P19_PortugalPor-to1994 (host: Oryctolagus cuniculus), which so far have not been assigned to any genotype. The strains (Table S7 at https://ojs.ptbioch.edu.pl/index.php/abp/) revealed a genetic distance from one another of 0.017/1.7% nucleotide substitutions.

DISCUSSION
Genetic variability of viruses depends on the frequency of mutations taking place in the genome, and external factors of the environment where the pathogens are present (Belshaw et al., 2008). The core natural factor that conditions a mutation comes from errors in imprecise replication of nucleic acids by RdRp. This polymerase is incapable of repairing errors occurring during synthesis of descendant RNA strands (Steinhauer et al., 1992). The scale of errors is illustrated by an example where one error occurs per one replication cycle for RNA viruses as compared to one error per 300 cycles for DNA viruses (Smertina et al., 2019). Genetic variability can also result from other mechanisms, such as homological recombination between closely related RNA viruses. The pace of recombination seems to vary and be specific to each virus species, as some recombine very frequently, while others seem to be limited, showing little proof of recombination (Forrester et al., 2008). The emergence of a new genotype or variant in the case of RNA viruses most frequently involves mutations within the RdRp gene or within the capsid genome, particularly the subunit encoding surface epitopes, with the new genotype being more or less virulent or more stable in the environment (Capucci et al., 1998;Belshaw et al., 2008;Bigoraj et al., 2011). Studies on genetic variability of RHDV, among others conducted on the basis of the VP60-coding gene, revealed both -mutations and recombinations (Abrantes et al., 2008;Forrester et al., 2008;Hukowska-Szematowicz et al., 2013;Lopes et al., 2015b;Mahar et al., 2016;Hu et al., 2017;Lopes et al., 2018). It is suggested that the rabbit plague outbreak recorded in China in 1984 was caused by recombination of the genetic material originating from the Angora rabbits imported from Germany. This variability type was also identified for three RHDV strains in the RdRp gene (Forrester et al., 2008). According to some authors, such RHDV variability type as recombination will be increasingly frequent in viral evolution because the recombinants can be characterised with higher pathogenicity and different tropism than the input strains (Abrantes et al., 2008;Forrester et al., 2008;Lopes et al., 2015b), and the phenomenon begins to play an im-portant role in generating biological diversity of RHDV2 (Lopes et al., 2015b).
Analysis of genetic variability in 105 Lagovirus europaeus strains, based on the RdRp gene, indicated that these strains feature an increasing genetic distance, observed both in the time scale and in geographic scale, while RHDV variability occurs both, through mutation and recombination, which is a consequence of employing the existing genetic diversity to develop new genome combinations . The results of the analyses point out that the highest mutation potential in the studied RdRp gene (15.8-16.5%) was recorded for the RHDV2 genotype in Portugal and Spain. The results indicate that these are non-recombinant RHDV2 strains. Smertina and others (Smertina et al., 2019) speculated that RHDV2 in Australia could have acquired a relative fast polymerase, which may explain its higher virulence, and apparent evolutionary success. During 18 months from arrival, RHDV2 had largely replaced the endemic RHDV strains in Australia (Mahar et al., 2018a). Therefore, replication speed may mark viral efficiency because viruses with a higher replication speed can produce more copies of its genome, resulting in a greater number of variants, even if RdRp error rate remains the same (Smertina et al., 2019). Lower mutation potential within the area of the RdRp gene was observed for strains classified as GI.1a (G6) (11.8-13.6%), GI.1c (G2) (1.0-6.8%) and GI.1d (G3-G5) (6.5-8.4%). The mutation potential, however, of the Australian strains in GI.1c (G2) from the period of 2007-2014 totalled 3.4-6.8%, which points to the fact that the mutation coefficient in RdRp had increased in the Australian population of the virus. A very interesting situation was found for some strains (PSM2, CBMad17-3, CBMad17-2, CBPico17-2, CBPico17-1, CBAlgarve14-3/2014) identified in 2014, 2016 and 2017 in GI.1b (G1), which showed variability at 11.1-12.9%, similar to GI.1a (G6/RHDVa) (11.8-13.6%), but phylogenetically were not found in this group. This suggests that these are G1/RHDV2 recombinants, which is in line with the results of previous studies (Mahar et al., 2016;Lopes et al., 2018;Silverio et al., 2018). Genetic variations of 7.0-10.5% of the remaining strains (CB194Pt, CB147 Pt, Eisenhuettenstadt, SD, AST89) in this group, supported by the results of studies by other authors (Abrantes et al., 2012;Lopes et al., 2017), indicate that these are the G1-non-recombinant strains, both within the studied polymerase gene and in the VP60 and VP10 genes . Genetic variability of 14.7% and 14.8% in the RdRp gene was observed in the case of Portuguese strains P95 and P19, which may have been the source of a new variant of polymerase which did not stabilize in the environment. On the other hand, very high genetic variability, which may indicate recombination of genetic material between rabbit viruses, was observed for CBMad-17-1 (37.8%), AUS/NSW/BER-3 (36.3%) and AUS/NSW/BER-2 and AUS/NSW/WAL-1 strains (36.2%).
The matrix of phylogenetic relations (Fig. 1) among the 105 Lagovirus europaeus strains on the basis of RdRp gene sequence has revealed a clear division of the strains into 7 groups. The allocation to the groups depended on chronology, geographic origin, and evolutionary events in the viral history, such as recombination. The location of some strains on the phylogenetic tree, supported by their variability results, testifies to them being recombinants, forming phylogenetically new groups, while recombination seems to be an increasingly frequent event in the viral evolution . Such a situation applies to two strains of the Portuguese origin: P95 and P19. On the phylogenetic tree, the strains are positioned between GI.1a-1d and RHDV2 and have not been assigned to any genotype. Their location on the phylogenetic tree indicates that this is a new group of recombinant strains which "for a moment" appeared in Portugal. Phylogenetic analysis by Lopes and others  showed that in the region of the genome encoding structural proteins VP60 and VP10, the strains were similar to pathogenic strains from the GI.1b group (G1), but in the region of the genome that encodes non-structural proteins they had created a new group which was 13% divergent from known strains. According to Lopes and others , there are no further reports of this strain surviving in the environment. The absence of information could be due to limited sampling . Another important issue is the fact that the P95 strain was isolated from an Iberian hare (Lepus granatensis), which points to breaching the species barrier, although for some unknown reason the infection was not sustained in the environment . Perhaps the occurrence of this strain line can be explained with the survival strategy of the virus through speciation in one of two directions, namely: change of the infected host species or through a change within the same host species. The former strategy which probably occurred in the case of P95 required many changes and was a major adaptation challenge, hence the virus did not survive. One must also account for the fact that perhaps, in the case of the two strains, there was a similar situation as that was observed in the noroviruses (Mahar et al., 2013), where recombination resulted in a new RdRp variant with a higher mutation potential (14.8% and 14.7% nucleotide substitutions, respectively, for P95 and P19). This was supposed to imply a higher genetic variability of RHDV and improved general robustness of the viral population against the selection pressure but, for unknown reasons, this mechanism has failed (Mahar et al., 2013).
The GI.1b (G1) group contains strains that are similar in two ways. On one hand, the groups were created by chronology, as it contains the "classical" strains of RHDV from 1989 and strains CB137 from 1995 and CB194 from 2006, which probably are evolutionarily related to classical strains (Abrantes et al., 2012). On the other hand, this group also includes the G1/RHDV2 recombinant strains currently circulating in the Portuguese islands and replacing the existing G1 strains in the wild rabbit population (Lopes et al., 2015a). It should be noted that the CBAlgarve14-3 strain, in a phylogenetic analysis based on the VP60 gene sequence (Lopes et al., 2015a;Mahar et al., 2016), showed positioning in the RHDV2 genotype, i.e. behaved differently than in our study. This fact underlines the importance of research not only on genes encoding structural proteins such as VP60, but also of non-structural proteins such as RdRp, in order to determine the designation of strains to genotypes. In group GI.1c (G2), strains were characterised with both, a chronological structure: from the period of 1987-1996, and a geographic structure: strains from Australia from the period of 2007-2014, and from New Zealand. Evolutionary relations of Australian and New Zealand strains point to them as being descendants of a European ancestor, the Czech V351 strain, which was artificially released in the mid-1990s in Australia and New Zealand to control the population of the artificially introduced European Rabbit (Cooke & Fenner, 2002;Eden et al., 2015). Furthermore, evolutionary relations in this group point out that all of the remaining strains forming this group derive from the Czech V351 strain. Group GI.1d (G3-G5) was predominantly formed by European strains from the period of 1989-2004, with some of them focusing around the Italian Bs89 strain from 1989. It must be pointed out that this group also included strains which, in prior studies (Forrester et al., 2008), were found to have undergone recombination in the RdRp-coding region, including: Frankfurt5, Frank-furt12, and Wika, but this event did not result in a durable new RdRp variant. Among the strains analysed, 28 originating from the period of 1996-2015 had formed a monophyletic group GI.1a (G6/RHDVa), statistically very well supported, with a variability of 11.8-13.6%, with intragroup variability totalling 0.2-6.1%.
Phylogenetic relations deduced from the RdRp gene had revealed a new phylogenetic group, gathering original non-recombined RHDV2 strains from the period of 2011-2016 circulating in Spain, continental Portugal and the Canary Islands Silverio et al., 2018). The last group in the phylogenetic tree was formed by four strains from Australia, one from Portugal, and by MRCV. Phylogenetic relations in this group suggest that currently the population of Lagovirus is mixed between Europe and Australia, as also reported by other authors (Mahar et al., 2016). Phylogenetic relations in this group and genetic variability suggest that the Portuguese CBMad17-1 strain from 2017 is a recombinant that arose between rabbit caliciviruses and RHDV2, although the occurrence of the strain in Madera still remains a mystery. Positioning of Australian strains AUS/ NSW/BER-3, AUS/NSW/BER-2, and AUS/NSW/ WAL-1 on the phylogenetic tree points to their phylogenetic relations with rabbit calciviruses, and suggests they are recombinants of rabbit caliciviruses and RHDV or rabbit caliciviruses and RHDV or RHDVa, because Australian RHDV strains have been gradually replaced since 2013 with the GI.1a variant (Mahar et al., 2018a).

CONCLUSIONS
The extent and significance of genetic variability through RNA mutation and recombination has not yet been thoroughly investigated, and the study presented here approximates this issue to some extent, which is important in the context of viral evolution, epidemiology and pathogenesis. So far, the main direction of genetic variation and evolution of RHDV which has been present worldwide for over 30 years, has led from RHDV, through RHDVa, to RHDV2. It seems that the recorded variability within the RdRp gene for the analysed lagoviruses, in particular RHDV2, is of major importance for its expansion and survival in the environment.