To date, there is no standardised classification of the subspecies of Puumala virus based on quantitative criteria. The current classification consists of several lineages such as North-Scandinavian (N-SCA), Finnish (FIN), Russian (RUS), Latvian (LAT), Alpe-Adrian (ALAD) and Central European (CE) (Sironen et al., 2001; Dekonenko et al., 2003; Castel et al., 2019; Blinova et al., 2022; 2023). Some sequences from the same lineage may be genetically more distant from each other than sequences from different lineages, as there are currently no strict thresholds for the determination of the lineages. Furthermore, some new sequences cannot be clearly assigned to one of the current lineages (Blinova et al., 2023). Here, we propose a threshold for genetic distance, to separate the lineages unambiguously. On the other hand, we show that such a quantitative criterion can be used to divide PUUV into distinct genotypes.

Phylogenetic reconstruction is only reliable for sequences without horizontal gene transfer, and for viruses, this means that only parts of the genome that do not undergo recombination can be used. Recombination occurs when at least two viruses co-infect the same host cell and exchange genetic segments. Different types of viral recombination are distinguished by the structure of the crossover site. Homologous recombination occurs at the same site in both parental strands, whereas non-homologous or illegitimate recombination occurs at different sites of the genetic fragments involved. Reassortment occurs in viruses with segmented genomes that can exchange complete segments (Pérez-Losada et al., 2015). Since PUUV has three segments, it is important to understand whether its genome undergoes reassortment and homologous recombination before conducting any studies on its phylogeny.

Homologous recombination is an important aspect of viral evolution in the context of phylogenetic analysis, as different parts of the genome can evolve independently. Since Puumala viruses have three RNA molecules in their virion they can also exchange whole genome segments.

Homologous recombination has been reported for hantaviruses. E.g., it occurs in Tula virus and Dobrava-Belgrade virus which are the relatives of PUUV (Sibold et al., 1999; Klempa et al., 2003). In addition, numerous reassortment events have been found in other hantaviruses (Kim et al., 2016; Lee et al., 2017). Due to their host specificity, there can only be one species of hantavirus in a host species, so that hantaviruses, unlike other bunyaviruses, can generally only undergo reassortment events within the same viral species (Klempa, 2018). Reassortment is well described for Puumala orthohantavirus. In the studies in Finland, a high frequency of reassortment, 19.1%–32%, was observed. The reassortment events within Finnish lineage, between Finnish and North Scandinavian groups have been found (Razzauti et al., 2008; 2009). The authors estimated the contributions of homologous recombination and reassortment and concluded that recombination is rarer than reassortment. In PUUV from Central Europe, sequences from Slovakia were clustered with sequences from Bohemia (Czech Republic) and Bavaria (Germany) only for M, but not for S segments. This means that the M segments of these viruses are closely related, but the S segments are not. The S segment of the Slovakian viruses was related to the S segments of PUUV from Belgium and France. This fact may indicate reassortment events (Razzauti et al., 2009). On the other hand, this study was performed for partial genomes, so these observations can also be explained by the presence of homologous recombination.

The presence of reassortment events means that the phylogenies for the different segments do not match. At the same time, different regions within the same segment may also have different evolutionary histories. Therefore, the segment without signs of recombination should be used for phylogenetic analysis. On the other hand, the segments are not evenly sequenced. Using the most frequently sequenced segments makes it possible to comprehensively uncover the phylogeny of PUUV.

Only viruses with three fully sequenced segments were selected for recombination and reassortment analysis (total number of sequences = 50). The GenBank database was searched for viruses for which all three segments were available. Then these genomic sequences were joined together with strings of 700 N’s between segments for each PUUV. In addition, only the protein-coding regions were used for the analysis. These sequences were aligned using the Muscle algorithm in MEGA software (Kumar et al., 2018) and then recombination and reassortment events were searched for using the online web application Shiny ¹ (Vakulenko et al., 2021). For each pair of genomic regions, the distances between all possible pairs of sequences were calculated and then plotted on a pairwise nucleotide distance comparison plot (PDCP) (Vakulenko et al., 2021). If no recombination or reassortment has taken place, the distance between two sequences in one region should correlate with the distance in another region. If recombination or reassortment has occurred, these distances may not correlate. Without recombination and reassortment, the distances between pairs of sequences for each pair of genomic regions should follow a linear relationship, so a high value of root mean square error (RMSE) indicates recombination or reassortment. A pairwise distance deviation matrix (PDDM) was then calculated for all possible pairs of genomic regions using a sliding window with a length of 700 nt and a step of 50 nt. The PDDM thus summarises the RMSE of all PDCPs. Red areas on the PDDM indicate a high RMSE value and thus the fact that the corresponding genome regions could originate from different ancestors. Blue areas indicate a lower probability of recombination or reassortment.

To detect recombination and reassortment event GiRaF (Nagarajan and Kingsford, 2011) and RDP4 (Martin et al., 2015) software was also implemented. GiRaf used as an input phylogenetic trees that were obtained using Bayesian algorithms by the MrBayes software (Huelsenbeck and Ronquist, 2001). In RDP4 seven methods were used to find possible recombination and reassortment events (RDP, GENECONV, Chimaera, MaxChi, BootScan, SiScan, and 3Seq). When looking for potential recombination inside each one of the segments with RDP4, we used alignments of sequences for that particular segment, and when analyzing reassortment events between different segments, we used a full-genome alignment with segments separated by N’s, as mentioned before.

The visualisation of the similarity was created using the SimPlot software (Lole et al., 1999). It shows the comparison between a set of sequences and a reference sequence.

The trees for the recombination and reassortment analysis were created using the maximum likelihood algorithm with 100 bootstraps with the MEGA software and visualised with the Figtree software (Rambaut, 2020).

All 514 currently available (as of June 30th, 2023) Puumala S segment sequences were obtained from the GenBank. Only sequences that contain a complete nucleocapsid coding fragment were used in the study (n = 505). A phylogenetic tree was built using the maximum likelihood method and branches were named according to currently accepted systematics.

Protein and nucleotide pairwise distances were calculated using MEGA software. Scatter plots in nucleotide distances/protein distances coordinates were built using R Studio [ggplot2 library (Wickham, 2016)]. CD-HIT (Fu et al., 2012) was used for sequence clustering to visualise the results obtained by rationally reducing the sample size. Then another maximum likelihood phylogenetic tree was built from the clustered sequences using IQ-TREE with ultrafast bootstrap analysis (n = 1,000) (Nguyen et al., 2015). The protein and nucleotide sequences are available at ² .

The S segment was the most represented in the GenBank segment (as of June 30th, 2023 number of sequences = 514) since it has the shortest length and it’s sequencing is the most cost-effective. L segment was the least represented segment (number of sequences = 65). RDP4 software was used to determine whether the S segment undergoes recombination events. No event was found that would have been supported by more than two methods out of the seven used.

Since sequences of other segments are not available for most viruses and the S segment doesn’t seem to undergo any recombination events, the phylogenetic analysis was based on the nucleocapsid protein gene located in the S segment (n = 505). The distributions of pairwise nucleotide and amino acid distances were plotted to visualize the sequence space occupied by PUUV (Figure 1). Each pair of sequences in a dataset was represented by a dot with coordinates reflecting the nucleotide and amino acid sequence distance between them. It should be noted that without strict quantitative criteria, the taxonomic classification may be ambiguous and may change as new distinct sequences are discovered. Each dot was coloured according to whether the virus pair consisted of inter- or intragroup viruses. The current classification of PUUV lineages was used in Figure 1A to define the inter- or intralineage assignment of the virus pairs (Figure 1A). At the same time, PUUVs were clustered based on the threshold of 10% of nucleotide sequences. These groups of sequences were used for the inter- or intracluster assignment of the virus pair in Figure 1B.

Two different groups of virus pairs on the plot could be clearly distinguished (Figure 1). One group of virus pairs has less than 10% nucleotide difference in the S segment sequence, while the other group differed in more than 10% of the nucleotides in the S segment sequence. Viruses from different groups always differed in more than 10% of the nucleotides, while some sequence pairs from the same group also differed by more than 10% of nucleotides (Figure 1A, cyan dots to the right of the threshold line). Indeed, sequence pairs from different lineages can be phylogenetically closer than a pair from one group. For example, #KX815394 (LAT lineage) and #JQ319167 (FIN lineage) are more similar (12.97% of nucleotides were different) than #MG812440 (CE lineage) and #AF294652 (CE lineage) (16.5% of nucleotides were different). Thus, in some cases, intralineage diversity is higher than interlineage diversity.

To test the validity of the “10% nucleotide difference” criterion usage as a quantitative measure for classifying Puumala virus, we divided the viruses into 21 groups (Figure 2). Remarkably, several lineages (CE, RUS, N-SCA, and FIN) were subdivided into multiple groups specific to the isolated regions according to this criterion (Table 1). At the same time, other lineages (ALAD, S-SCA, and DAN) correspond to the novel taxon. This means that the nucleotide distance between the coding regions of the S segments of two sequences can be the criterion for an unambiguous classification of PUUV into the subspecies taxon. At the same time, such grouping of PUUVs has its limitations (cyan dots to the right of the “10% nucleotide difference” criterion in Figure 1B). #JN657228 (LAT lineage) has a nucleotide difference of less than 10% with the sequences MT580935 and KX815394, but more than 10% with other representatives of the LAT lineage.

The accession numbers of all S segment PUUV sequences assigned as “lineages” and “genotypes” are provided in Supplementary Table S1.

The concatenated sequences of protein-coding regions of all three segments were used (number of sequences = 50). The results of reassortment and recombination analysis are shown on the pairwise distance deviation matrix (PDDM) in Figure 3. The coordinates of three genes in the final alignment with separations were: 1–1,302 for N gene, 2003–5,449 for G gene and 6,150–12620 for RdRp gene ³ .

This plot uses red colour to mark areas of high phylogenetic inconsistency among the genomic regions, while blue indicates low phylogenetic inconsistency. A high RMSE value suggests potential reassortment among the most divergent viruses in the respective genomic regions.

For the next part of the analysis, an alignment without the separation between the segments was used ⁴ . We visualised PDCPs for the following pairs of regions: S segment (1–1,302 nt in alignment) against the M segment (1,303–4,752 nt in alignment) (Figure 4B) and M segment (1,303–4,752 nt in alignment) against L segment (4,753–11,223 nt in alignment) (Figure 4C). Each sequence pair in a dataset was represented by a dot with coordinates reflecting the nucleotide sequence distance between the different genomic fragments. The density of possible distance values was indicated by colour. The PDCP for the S segment and the M segment showed that some viral pairs had phylogenetic incongruence (Figure 4B, circled) that could be explained by a reassortment event. Indeed, the S segments of these viruses differed in about 6% of the nucleotides in the S segment and in about 14% of the nucleotides in the M segment. To obtain a negative control, we concatenated the S segment and the M segment, selected odd or even sites from this artificial sequence and calculated the pairwise distances for this data set (Figure 4A).

Two concatenated viral genomic sequences that were off trend (Figure 4B, circled): “strain K27” (#MZ673552-MZ673554) and “isolate 989” (#OL343606, OL343584, and OL343562) were selected for further analysis. The similarity plot for this pair of sequences and three control sequences: “isolate Klyukva” (#MZ580951-MZ580953), “strain Osnabruck V29” (#MN639737-MN639739) and Muju virus “isolate 11-1” (#JX028271-JX028273) was shown in Figure 5A. The sequence of “isolate 989” was used as a reference. Phylogenetic trees of these five sequences were constructed using the S segment and the M segment (Figure 5B). Both trees showed high bootstraps but different topologies. The S segment of the sequence “isolate 989” was close to the S segment of the sequence “strain K27,” while the M segment of these viruses was more diverse.

The S segment of strain K27 has 94.5% identical nucleotides with the segment of isolate 989, while its M and L segments have 86.9% and 87.4% identical nucleotides respectively. The most plausible explanation for this fact is a reassortment event in which the S segment was transmitted from a precursor virus whose descending M and L segments are currently unknown. However, this is speculation, and without direct evidence, the drastic change in the proportion of identical nucleotides in the different segments can also be explained by the regional differences in mutation rates. According to the similarity plot, the different regions of the M segment have varying percentages of identical nucleotides (Figure 5A), and the phylogenetic trees of these regions also have different topologies ⁵ . In addition, conflicting tree topologies were observed for the phylogenetic trees for these M-segment fragments generated for all PUUVs from the dataset ( ⁶ , red arrows). These observations can be explained by possible recombination events in the M-segment.

In addition, GiRaF and RDP4 were used. GiRaF found no potential reassortment events, and RDP4 confirmed the reassortment event described in this study with three of the seven methods used ⁷ . To further validate the results, phylogenetic trees for each segment were built. Preliminarily all of the segments were analyzed by RDP4 for possible recombination events. Recombinant sequence “isolate 836” (#OL343604, OL343582, OL343560) was omitted from the analysis. The phylogenetic trees ⁸ for the S segment and the L segment have the same topology for the clades containing sequences “isolate 989,” “K27” and “Klyukva” that were previously mentioned, but for the M segment the topology is different (red arrows on the phylogenetic tree). Another potential reassortment event was discovered during the same analysis between clades that contain the sequences “Suo/renal,” “Sotkamo,” and “Mg7/2008,” respectively (blue arrows). The topologies for the S and M segment are identical for those clades, and the one for the L segment is different from them. Interestingly, an event between the same clades was discovered by RDP4 as well, although it suggests that the M segment reassorts, not L.