Phylogeographic diversity and hybrid zone of Hantaan orthohantavirus collected in Gangwon Province, Republic of Korea

Background Hantaan orthohantavirus (Hantaan virus, HTNV), harbored by Apodemus agrarius (the striped field mouse), causes hemorrhagic fever with renal syndrome (HFRS) in humans. Viral genome-based surveillance at new expansion sites to identify HFRS risks plays a critical role in tracking the infection source of orthohantavirus outbreak. In the Republic of Korea (ROK), most studies demonstrated the serological prevalence and genetic diversity of orthohantaviruses collected from HFRS patients or rodents in Gyeonggi Province. Gangwon Province is a HFRS-endemic area with a high incidence of patients and prevalence of infected rodents, ROK. However, the continued epidemiology and surveillance of orthohantavirus remain to be investigated. Methodology/Principal findings Whole-genome sequencing of HTNV was accomplished in small mammals collected in Gangwon Province during 2015–2018 by multiplex PCR-based next-generation sequencing. To elucidate the geographic distribution and molecular diversity of viruses, we conducted phylogenetic analyses of HTNV tripartite genomes. We inferred the hybrid zone using cline analysis to estimate the geographic contact between two different HTNV lineages in the ROK. The graph incompatibility based reassortment finder performed reassortment analysis. A total of 12 HTNV genome sequences were completely obtained from A. agrarius newly collected in Gangwon Province. The phylogenetic and cline analyses demonstrated the genetic diversity and hybrid zone of HTNV in the ROK. Genetic exchange analysis suggested the possibility of reassortments in Cheorwon-gun, a highly HFRS-endemic area. Conclusions/Significance The prevalence and distribution of HTNV in HFRS-endemic areas of Gangwon Province enhanced the phylogeographic map for orthohantavirus outbreak monitoring in ROK. This study revealed the hybrid zone reflecting the genetic diversity and evolutionary dynamics of HTNV circulating in Gangwon Province. The results arise awareness of rodent-borne orthohantavirus diseases for physicians in the endemic area of ROK.

Phylogeographic analysis has become an essential tool for the public health surveillance and molecular epidemiology of infectious diseases when used for tracing the sources of epidemic infections [7]. Recently, the phylogenetic association between patients with HFRS and natural reservoirs demonstrated the putative infection location of HTNV [8]. Active surveillance in HFRS-endemic areas identified the infectious source of HTNV by real-time next-generation sequencing (NGS), epidemiological interview, and targeted rodent trapping [9]. Emerging orthohantavirus infections may occur at any time through contaminated urine, feces, or saliva in rodent-infested areas. To ascertain geographic prevalence and disease risk assessment of orthohantavirus in HFRS-endemic areas, ROK, genetic and molecular epidemiological studies on small mammals have consistently been conducted for decades [8,[10][11][12][13][14][15][16][17]. Most studies have demonstrated the serological prevalence and genetic diversity of orthohantaviruses collected from HFRS patients or rodents in Gyeonggi Province [11][12][13][14][15][16]. In Gangwon Province, an administrative province in northeast ROK, approximately 371 HFRS cases have been reported from 2001-2019 [18]: Cheorwon-gun and Hwacheon-gun are highly HFRS-endemic areas. However, the continued surveillance and genetic study of HTNV remain to be investigated owing to the lack of virus genome sequences obtained from small mammals in Gangwon Province.
Hybrid zones are geographic regions in which two well-divergent taxa meet and mix [19]. The areas take various shapes from large scales of overlap to narrow contact regions or mosaic zones [20]. The particular zones are considered stable over evolutionary time [21]. The regions are characterized by variation in the genotype frequency across hybridization areas. A cline is formed by the frequencies of different genetic or phenotypic traits across hybrid regions. The contact zones continued by maintaining a balance between the homogenizing effect of dispersal and the diversifying effect of natural selection [22,23]. Hybrid zones were discovered in a variety of wild organisms [24][25][26][27]. Several studies have investigated the occurrence of hybrid zone between pathogens and their hosts in nature. The contact area plays a role in genetic diversity and speciation process in beak and feather disease viruses and their reservoirs [28]. Murine cytomegalovirus and their hosts formed a spatial contact in a hybrid zone maintained by natural selection [29]. Recently, a hybrid zone analysis revealed two distinct Tula virus (TULV) lineages were co-circulated in a geographical area in which two different evolutionary clades in the common vole (Microtus arvalis) interact and interbreed [30]. However, to the best of our knowledge, the hybrid zone of orthohantavirus in the ROK remains unknown.
In this study, epidemiological surveys of small mammals demonstrated the geographic prevalence of HTNV in Gangwon Province during 2015-2018. Targeted-enriched NGS elicited whole-genome sequences of HTNV acquired at the new expansion sites. The phylogeny of HTNV showed highly divergence and possible genetic exchanges of tripartite genomes in nature. In addition, the cline analysis suggested a hybrid zone of HTNV in Cheorwon-gun and Hwacheon-gun. This study provides significant insights into the phylogenetic diversity and evolutionary dynamics of orthohantaviruses in HFRS-endemic areas, ROK.

Ethics statement
All trappings of small mammals were carried out in accordance with the Ethical Guidelines of the Korea University Institutional Animal Care and Use Committee (KUIACUC #2016-49). Captured animals were euthanized via cardiac puncture under isoflurane anesthesia. To minimize hazards from potentially infected animals, workers wore thick rubber gloves and protective clothes. Collected rodents were placed in double plastic bags and transported to Korea University. Gloved hands were washed with a suitable disinfectant [31]. The experiment was conducted in the biosafety level 3 facility at Korea University.

Mitochondrial DNA analysis
Total DNA was extracted from liver tissues using a High Pure PCR Template Preparation Kit (Roche, Basel, Switzerland). To verify a precise species of small mammals, the cytochrome b gene of mitochondrial DNA was amplified by universal primers; forward: 5'-CGA AGC TTG ATA TGA AAA ACC ATC GTT G-3' and reverse: 5'-CTG GTT TAC AAG ACC AGA GTA AT'-3' [33]. The HTNV-positive A. agrarius was phylogenetically confirmed by mitochondrial DNA cytochrome b gene.

Indirect immunofluorescence antibody (IFA) test
Sera from A. agrarius were diluted 1:32 in phosphate buffered saline (PBS) [34]. The sera were added to wells of acetone-fixed Vero E6 cells infected with HTNV. The slides were incubated at 37˚C for 30 min. After washing with PBS and distilled water (D.W.), fluorescein isothiocyanateconjugated goat anti-mouse immunoglobulin G (IgG) antibody (ICN Pharmaceuticals, Laval, Quebec, Canada) was added. The plates were incubated at 37˚C for 30 min and then washed with PBS and D.W. The pattern of virus-specific fluorescence was evaluated to be an indication of HTNV infection using a fluorescent microscope (Axioscope, Zeiss, Berlin, Germany).

RNA extraction and reverse transcription-polymerase chain reaction (RT-PCR)
Total RNA was extracted from lung, liver, kidney, and spleen tissues using a FastPrep-24 5G Instrument (MP Biomedicals, USA) with TRI Reagent Solution (Ambion, Austin, Texas).  Reverse transcription was conducted with 1 μg of total RNA using a High Capacity RNA-to-cDNA kit (Applied Biosystems, Foster City, CA, USA) with random hexamers and OSM55 (5 0 -TAG TAG TAG [37,39]. First and second RT-PCR was performed at 94˚C for 5 min, followed by 6 cycles of denaturation at 94˚C for 30 s, annealing at 37˚C for 40 s, and elongation at 72˚C for 1 min; then, 32 cycles of denaturation at 94˚C for 30 s, annealing at 42˚C for 40 s, and elongation at 72˚C for 1 min and a final cycle at 72˚C for 5 min (ProFlex PCR System, Life Technology, CA, USA). PCR products were purified using a MinElute PCR purification kit (Qiagen, Hilden, Germany) or a QIAquick Gel Extraction Kit (Qiagen). Sequencing was performed in forward and reversed directions of each PCR product using a BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) on an automated sequencer (ABI 3730XL DNA Analyzer, Applied Biosystems). The whole-genome sequences of HTNV strains deposited in GenBank (Accession numbers: MT012546-MT012581).

Real-time quantitative PCR (RT-qPCR)
RT-qPCR was performed at 95˚C for 10 min, followed by 40 cycles at 95˚C for 15 s and 60˚C for 1 min using a Power SYBR Green PCR Master Mix (Applied Biosystems) on a QuantStudio 5 Real-Time PCR System (Applied Biosystems). The primer sequences of HTNV S segment included a forward primer (5'-TTA TTG TGC TCT TCA TGG TTG C-3') and a reverse primer (5'-CAT CCC CTA AGT GGA AGT TGT C-3') [40].

Multiplex PCR-based next-generation sequencing (NGS)
cDNA was amplified using HTNV-specific primer mixtures and Solg 2X Uh-Taq PCR Smart mix (Solgent, Seoul, Republic of Korea) according to the manufacturer's instruction. The enrichment was performed in 25 μL of reaction mixtures containing 12.5 μL of 2X Uh premix, 2.0 μL of each primer mixture, 10.5 μL of D.W., and 1.0 μL of DNA template. Multiplex PCR was performed by a cycle at 95˚C for 15 min, then 40 cycles and/or 25 cycles at 95˚C for 20 s, 50˚C for 40 s, 72˚C for 1 min, and a cycle at 72˚C for 3 min.
DNA libraries were prepared using a TruSeq Nano DNA LT sample preparation kit (Illumina, San Diego, USA) according to the manufacturer's instructions. To obtain size-selected amplicons, cDNA templates were mechanically sheared using the M220 focused ultrasonicator (Covaris, Woburn, MA, USA). The cDNA amplicons were prepared by size-selection, A-tailing, and ligation with indexes and adaptors. The enrichment reaction contained 5 μL of PCR primer cocktail and 20 μL of enhanced PCR mixture. The quality of libraries was evaluated by a bioanalyzer using an Agilent DNA 1000 chip kit (Agilent Technologies, Santa Clara, USA). NGS was performed using MiSeq reagent V2 (Illumina) with 2×150 bp of MiSeq benchtop sequencer (Illumina).

Rapid amplification cDNA ends (RACE) PCR
The 3' and 5' end sequences of viral genomes were acquired by rapid amplification cDNA end (RACE) PCR using a SMARTer RACE 5'/3' Kit (Clontech Laboratories Inc., Mountain View, CA, USA) according to the manufacturer's specifications. Both 3' and 5' ends of HTNV strains were empirically filled up by incomplete complementary sequences [41].

Phylogenetic analysis
Whole-genome sequences of HTNV were aligned by the Clustal W algorithm (Lasergene program version 5, DNASTAR Inc. Madison, WI). The phylogenetic trees of HTNV were generated using the best-fit General Time Reversible (GTR) +gamma (G) +invariable (I) (for L and M segments) and T92+G (for S segment) models of evolution. Support for the topologies was assessed by bootstrapping for 1,000 iterations. Model optimizations were calculated for each data set, followed by the calculation of pairwise genetic distances between HTNV strains using MEGA 7.0 [42].

Cline analysis
Geographic clines were inferred along one-dimensional transect axes crossing the contact areas for two diverged HTNV clades using the HZAR package [43]. The software includes functions for fitting molecular genetic or morphological data from hybrid zones to classic equilibrium cline models using the Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm. HTNV strains were analyzed by performing 10 6 generations of an MCMC sampling after 10 5 burn-in iterations. Input files for cline analysis include information on locality distance, genotype frequency, and the number of samples. The hybridization frequency of trapping sites was displayed starting from the western-most locality in Gyeonggi Province. Genotype frequency data objects were created using the function hzar.doMolecular-Data1DPops in the software. The hzar.plot.obsData function plots data with mean frequencies for molecular clines and mean values for morphological clines.

Genetic reassortment analysis
The graph incompatibility based reassortment finder (GiRaF) was performed to confirm reassortment events [44]. Nucleotide alignments of HTNV tripartite segments were used as an input source for MrBayes [45]. The best nucleotide substitution models were determined using MEGA 7.0. As input for the software, 1,000 unrooted candidate trees were estimated using the GTR+G+I substitution model, the burn-in 50,000 iterations (25%), and sampling every 200 iterations. These trees were used to simulate the phylogenetic uncertainty for segments; the parameters of the GiRaF were at the default settings. The default value of the confidence threshold was 0.7 for the data set; all events described using GiRaF were at over 0.9 confidence levels [46]. The process was repeated 10 times, with 10 independent MrBayesbased tree data per segment.

Multiplex PCR-based NGS of HTNV from A. agrarius captured in Gangwon Province
The coverage of genome sequences of 15 HTNV strains were 94.7-99.8% for L segments, 97.4-99.6% for M segments, and 99.2% for S segments (Table 3). In total, whole-genome sequences of 12 HTNV strains were nearly obtained in Cheorwon-gun, Chuncheon-si, Hwacheon-gun, and Yanggu-gun. The whole-genome sequencing of Aa17-338 showed low coverages of the L and M segments corresponding to the lowest viral load. We acquired the 3' and 5' termini sequences of Aa17-421 using RACE PCR to complete the whole-genome sequencing of HTNV.

Hybrid zone analysis
A cline analysis inferred the change in the population frequency of diverged clades, defined by Gyeonggi and Gangwon Provinces, along geographic transects ( Table 4). The frequency variation was estimated according to each segment, and this analysis showed that the pattern of transition between HTNV clades in the L and S segments was spatially similar (Fig 3). The transition pattern of the M segment was heterogeneous and significantly steeper compared to other segments. Geographic contacts were found between two phylogenetic lineages of HTNV in the ROK. In the L and S segments, hybrid zones were observed in Cheorwon-gun and Hwacheon-gun (64.3-89.4 km from Paju-si) whereas contact regions were detected only in Hwacheon-gun in the M segment. The range of the hybrid zone is approximately 20-25 km and  Fig 4). Analysis of the genotype frequency for HTNV L, M, and S segments was divided into three clades; HTNV strains from Paju-si, Pocheon-si, Yeoncheon-gun, and Chuncheon-si showed the frequency of clade was Ada a of 1 (Genotype of Gyeonggi) as referred to the Clade I. The clade frequency of HTNV strains from Yanggu-gun belonged to  Ada b of 1 (Genotype of Gangwon), as referred to the Clade II. Finally, Cheorwon-gun and Hwacheon-gun were selected as a hybrid zone between the clades I and II.

Genetic reassortment of HTNV
The occurrence of genetic exchange between segments (L and M, L and S, and M and S) was confirmed using the GiRaF program. The GiRaF analysis demonstrated that the reassortment event was detected between the M and S segments of HTNV in Aa15-56 and Aa15-58. The reassortment events were found in 100% of independent GiRaF runs by 10 times, and all events were supported over 0.9 confidence levels. However, additional genetic exchanges were not identified between L and the other two segments.

Discussion
Epidemiological prevalence and genetic diversity of orthohantaviruses in nature reservoirs play a critical role in understanding of hantaviral diseases in HFRS-endemic areas. In 2005, HTNV genome sequences acquired from four HFRS patients at training sites near the Demilitarized Zone, ROK, showed an epidemiological link with viral sequences obtained by targeted rodent trapping at six training sites where patients had exercised [13]. These results help early diagnosis and prevention of HFRS patients by establishing a database to infer epidemiological and emergent dynamics of hantaviral genomes sequenced in endemic regions [47]. To intensify the resolution of phylogeographic map of orthohantaviruses, whole-genome sequences of HTNV were recovered from small mammals in Gangwon Province during 2015-2018. Additional genome sequences from three regions, Cheorwon-gun, Chuncheon-si, and Hwacheongun in Gangwon Province, were obtained to enhance the geographical data available for surveilling HFRS incidence by genomic epidemiology (Fig 2). Here, HTNV strains were first detected in Chuncheon-si. Thus, there should be vigilance for potential human orthohantavirus infections in the region. Hybrid zones are geographic regions in which genetically divergent populations contact and mix [23]. The areas allow to observe interactions with two different genotypes and the process of speciation [21,48]. The two diverged clades of PUUV were circulated in a single host lineage, and inter-lineage reassortments of PUUV were detected in northern Finland [26]. The hybrid zone of TULV demonstrates correlation between the distribution of two phylogenetic lineages in TULV and their host clades in the European common vole (M. arvalis) [30]. The spatial transition between TULV lineages was narrower than clades of their reservoirs. In this study, spatial contact regions of HTNV were investigated through the cline analysis. The results revealed a remarkable hybrid zone at two sites (Cheorwon-gun and Hwacheon-gun) in Gangwon Province and showed the spatial separation and sequence divergence across genome segments of HTNV (Fig 3). Although the impact of fitness differences between individuals in the local population cannot be ignored, one might expect the hybrid zone to move in favor of the fitter genotype [49][50][51][52]. The genotype frequency (Ada a of 1) of the M segment in Cheorwon-gun indicated the M segment is completely compatible with the clade I. In Hwacheongun, the HTNV M segment showed a selected genetic phenotype of the clade I with the Ada a of 0.667 and Ada b of 0.333. These results demonstrated that the M segment of HTNV may be preferentially compatible with the clade I in the hybrid zone. The genetic lineage of HTNV L segment in Cheorwon-gun showed higher frequency of Ada b (Genotype of Gangwon) compared with the M and S segments. The results showed that HTNV L segment may be compatible with the L segment depending on the geographic origin.
Segmented RNA viruses have the capacity to exchange genome segments by genetic reassortment [53]. Several orthohantaviruses, including the PUUV, Sin Nombre virus (SNV), and HTNV, have been described with the intra-and inter-lineage reassortment between closely related variants [26,[54][55][56]. In the previous study, reassortment analysis reported that a genetic exchange in HTNV might occur in nature [18]. The genetic exchange of SNV was reported on the M segment in nature and in vitro [57,58]. Andes virus also has high level molecular diversity in the M segment, resulting in five diverged clades related to geographic origins in the South Americas [59]. These reports consistently showed that the possible exchanges of the M segment had lower genetic compatibility requirements and higher genetic tolerance than the L and S segments [49,60,61]. However, a homologous association between the L and M segments was observed in PUUV [62]. Genetic reassortments of HTNV strains (2/11; 18.2%) in the hybrid zone was detected using GiRaF analysis. The GiRaF analysis demonstrated that Aa15-56 and Aa15-58 from Cheorwon-gun were reassortants that have phylogenetic heterogeneity between the M and S segments. However, additional reassortants were not found between the L and the other two segments. The genetic compatibility and homologous relationships among hantaviral genome segments remain to be studied.
This study led us the hypothesis that genetic exchange rates of HTNV may be associated with the hybrid zone. Co-circulation of different lineages within hybrid zones presented the opportunity to confer genetic exchanges and dominance of the virus genome between these two genetic groups. The genetic exchanges confer viral characteristics including fitness, transmission, and pathogenesis [63]. The genetic reassortments of different segments promoted the evasion of host immunity and the occurrence of epidemics in influenza A virus and rotavirus A [64][65][66]. Continued reassortment of human immunodeficiency virus gave rise to genetic variants with transmissibility and altered virulence in humans [67]. A genetic exchange of the M segment of DOBV occurred between low pathogenic DOBV-Aa and highly pathogenic DOBV-Af [68]. In this study, HTNV in Cheorwon-gun showed genome organization compatible with reassortments between Gyeonggi and Gangwon Provinces. The cline analysis demonstrated that Cheorwon-gun and Hwacheon-gun were a hybrid zone among HTNV population in northern areas, ROK. However, some questions remain to be further investigated: 1) whether the genetic exchange of hantaviruses associates with the co-existence of different genetic lineages; 2) whether the biological consequence of the reassortment affects the pathogenicity in humans.
In summary, the prevalence and distribution of HTNV in HFRS-endemic areas of Gangwon Province enhance the phylogeographic maps for the monitoring and response to orthohantavirus outbreaks in the ROK. The hybrid zone analysis reveals hybridization of HTNV strains and the reassortment analysis suggests a natural occurrence of HTNV genetic exchange in the ROK. The finding provides significant insights into the genetic diversity and evolutionary dynamics of HTNV circulating in Gangwon Province. This report increases awareness of the rodent borne-orthohantavirus disease for physicians in the endemic area of ROK.
Supporting information S1 Table. Trapping sites of small mammals colleted in Gangwon Province. (DOCX)

S1 Fig. Immunofluorescence on Hantaan virus (HTNV)-infected Vero E6 cells. (A)
The representative HTNV-infected sample (Aa15-69) was shown in the S1 Fig. (B) The negative control was confirmed by PBS. (C) The positive control was an image of Vero E6 cells infected with HTNV 76-118 (the prototype of orthohantavirus