The exchange of genetic material between individuals is a key driver of evolution and diversification across most branches of life. Segmented viruses can exchange genetic material through reassortment of genomic segments. New viral strains that emerge from reassortments can have greater infection ranges and higher virulence, although concrete examples of the adaptive advantages of reassortants in nature apart from influenza remain rare. We studied here the evolutionary history and consequences of reassortment in Tula orthohantavirus (TULV) in a hybrid zone between evolutionary lineages of its reservoir host, the common vole (Microtus arvalis). Across 58 trapping sites and 127 infected voles, we detected 27 TULV reassortants in a 12.5 km broad zone at the contact of the parental TULV clades, resembling a viral hybrid zone concordant with the hosts’. Phylogenomic analyses revealed three independent reassortment events, but most of the host hybrid zone was dominated by a single strain with a reassorted M-segment, which encodes the surface glycoprotein. We detected clade-specific variation in the glycoprotein’s N-terminal region consisting of five residues, two of which showed evidence of positive selection. In silico 3D modeling of seven glycoproteins confirmed that this N-terminal region has a unique and specific structure for each TULV clade and the dominant reassortants and is the only structurally variable region of the TULV glycoprotein. Our findings suggest that reassortment between the parental TULV clades in the contact region has resulted in a transgressive virus phenotype potentially adapted to hybrid hosts. This demonstrates the potential of zones of hybridization for the emergence of new virus strains with novel evolutionary trajectories.

Homology modeling, structural alignments and molecular dynamics simulations

To assess structural differences in the M-segment glycoprotein among TULV clades and reassortants, we conducted homology modelling and 3D visualization of the TULV M-segment and its N-terminal ectodomain. We selected six M segment sequences which reflect the full diversity spectrum for AA in the N-terminal ectodomain from the Saxony transect, and the original TULV genome Moravia (Kukkonen et al., 1998) as a TULV-EST.S reference. Homology modelling and molecular dynamics (MD) simulations were performed in YASARA Structure version 20 using adapted macros and structures implemented from (Krieger & Vriend, 2014, 2015). Homology models were built using the standard “hm_build.mcr” macro based on the hantavirus glycoprotein structure “6ZJM” template (Serris et al., 2020). Resulting homology models displayed three different subunit configurations (four in Configuration 1, two in Configuration 2, one in Configuration 3). These different configurations likely reflect stochastic variations in fitting or different states of the protein, as the hantavirus glycoprotein can transition through different structural variations (Rissanen et al., 2017). It is unlikely that these variations reflect actual changes in protein folding, as we even observed different configurations between the near identical TULV-EST.S proteins from the Saxony transect and our reference Moravia (4 AA difference).

Pairwise structural alignments, including RMSDs and sequence identities between the homology models were calculated using MUSTANG (Konagurthu et al., 2006). RMSD of atomic locations across the entire protein tetramer (residue 13-1104) were compared to the RMSD of first nine AA of the Gn-subunits ectodomain (N-terminus). Only the results of the alignment for the four proteins in Configuration 1 are displayed in Supp. Table S10, because structural alignments require identical configuration of all subunits for meaningful evaluation. Individual subunits of the three remaining proteins, when aligned, showed similar levels of low RSMDs, but randomness in the predicted configurations artificially inflates the pairwise RMSDs of the complete proteins. This problem can normally be circumvented by using the configuration of a predicted protein for all future predictions. This was not applicable for our set of proteins, as it would also artificially force the flexible N-terminal ectodomain to change shape to match the reference thus eliminating the key variation between proteins. 

MD simulations of all N-terminal regions were performed using the “md.run.mcr” macro and AMBER14 force field, with a simulation duration of 500 ns (Maier et al., 2015). Time simulation steps were set to 1.35 fsec. The simulation box was 'Cube”-shaped and extended at least 10 Å to each side of the model (extension=10), was filled with 0.9% NaCl and the TIP3P water model was used at physiological pH 7.4. Further settings were: temperature at 298K, pressure at 1 bar, density = 0.997, cutoff 8Å- periodic cell boundary and longrange coulomb forces (particle-mesh Ewald). The solute was kept from diffusing and crossing periodic boundaries using the CorrectDrift function (Krieger & Vriend, 2014, 2015). The four flanking amino acids were fixed as anchor point during the MD simulations. To determine structural differences between the regions of interest, root mean square fluctuations (RMSF) were calculated by the “md_analyze” macro. The RMSF is the fluctuation of every heavy atom compared to the mean structure within a simulation cell. The average RMSF of the constituent atoms is used to determine the RMSF per solute residue.