A novel genotype of Hantaan orthohantavirus harbored by Apodemus agrarius chejuensis as a potential etiologic agent of hemorrhagic fever with renal syndrome in Republic of Korea

Background Orthohantaviruses, causing hemorrhagic fever with renal syndrome (HFRS) and hantavirus cardiopulmonary syndrome, pose a significant public health threat worldwide. Despite the significant mortality and morbidity, effective antiviral therapeutics for orthohantavirus infections are currently unavailable. This study aimed to investigate the prevalence of HFRS-associated orthohantaviruses and identify the etiological agent of orthohantavirus outbreaks in southern Republic of Korea (ROK). Methodology/Principal findings We collected small mammals on Jeju Island during 2018–2020. We detected the Hantaan virus (HTNV)-specific antibodies and RNA using an indirect immunofluorescence assay test and reverse transcription-polymerase chain reaction on Apodemus agrarius chejuensis (A. chejuensis). The prevalence of anti-HTNV antibodies among rodents was 14.1%. A total of six seropositive mouse harbored HTNV RNA. The amplicon-based next-generation sequencing provided nearly full-length tripartite genomic sequences of six HTNV harbored by A. chejuensis. Phylogenetic and tanglegram analyses were conducted for inferring evolutionary relationships between orthohantaviruses with their reservoir hosts. Phylogenetic analysis showed a novel distinct HTNV genotype. The detected HTNV genomic sequences were phylogenetically related to a viral sequence derived from HFRS patient in southern ROK. Tanglegram analysis demonstrated the segregation of HTNV genotypes corresponding to Apodemus spp. divergence. Conclusions/Significance Our results suggest that A. chejuensis-borne HTNV may be a potential etiological agent of HFRS in southern ROK. Ancestral HTNV may infect A. chejuensis prior to geological isolation between the Korean peninsula and Jeju Island, supporting the co-evolution of orthohantaviruses and rodents. This study arises awareness among physicians for HFRS outbreaks in southern ROK.

Emerging orthohantaviruses pose a significant public health threat worldwide. According to the Korea Disease Control and Prevention Agency, approximately 400 HFRS cases are recorded annually in the Republic of Korea (ROK), with a mean mortality rate of 1-4% [7]. In 1976, HTNV harbored by the striped field mice (Apodemus agrarius), was first identified in the ROK [8]. In addition, five types of hantaviruses have been characterized: SEOV carried by Rattus norvegicus and R. rattus; Soochong virus (SOOV) carried by A. peninsulae; Muju virus carried by Myodes regulus; Imjin virus carried by Crocidura lasiura, and Jeju virus carried by C. shantungensis [9][10][11][12][13]. Previous studies have reported genomic and phylogeographic relationships between the patients with HFRS and A. agrarius collected in endemic regions [14]. However, most epidemiological surveillance studies have been performed in the northern ROK, which limits the understanding of the correlation between HFRS clinical cases and their reservoirs in southern areas [15][16][17][18]. Although a few studies have demonstrated the prevalence of orthohantavirus in the southern ROK, epidemiological investigations of orthohantavirus outbreaks in humans and rodents in the southern ROK were limited due to the lack of available viral sequences [19][20][21][22][23]. On Jeju Island, an administrative island in southern ROK, approximately 18 HFRS cases have been reported from 2011 to 2019 [7]. Based on the incidence of HFRS, one or more distinct, rodent-borne orthohantaviruses have been suspected to exist on the island. However, the etiological agent of HFRS on Jeju Island still remains unknown.
In this study, we report the nearly full-length genome sequences and characterize a genetically distinct orthohantavirus, harbored by the subspecies A. agrarius chejuensis (A. chejuensis) captured on Jeju Island. The phylogenetic association of A. chejuensis-borne HTNV with viral sequences from HFRS patient suggests that A. chejuensis-borne HTNV may be a potential etiological agent of HFRS on Jeju Island, ROK. This study provides significant insights into the genetic diversity and molecular evolution of HTNV in the ROK. We highlight the necessity of continued surveillance and large-scale investigations for understanding the evolutionary diversification, transmission, and pathogenicity of orthohantaviruses in the southern areas, ROK.

Trapping for small mammals on Jeju Island, ROK
Rodents and shrews were captured on Jeju Island, ROK from 2018 to 2020 (Fig 1).

Whole-genome sequencing of HTNV
Using multiplex PCR-based next-generation sequencing (NGS), the nearly whole-genome sequence of HTNV was obtained from the lung tissues of A. chejuensis captured on Jeju Island. where small mammals were captured on Jeju Island from 2018 to 2020. The red circles indicate the HTNV RNA positive areas: a, Ora-dong; b, Bongseong-ri in Jeju-si, respectively. The black circle represents the regions where no HTNV RNA was identified: c, Ara-dong; d, Aewol-eup; e, Haengwon-ri; f, Youngsu-ri. The Quantum Geographical Information System (QGIS) 3.10 for Mac was used to create the map.

Phylogenetic and co-evolutionary analyses of HTNV
The phylogenetic analysis of the HTNV L, M, and S segments showed well-supported geographical clusters (Fig 3). The L segment of HTNV from Jeju Island formed an independent geographic and distinct cluster with all other viral genomes collected from mainland ROK and China. The M segment of HTNV from Jeju Island shared a common ancestor with the HTNV from mainland ROK. The S segment of HTNV from Jeju Island formed an independent genetic lineage with other HTNV, SOOV, and Amur virus (AMRV). In addition, the phylogenetic analysis of the partial HTNV L, M, and S segments revealed well-supported geographical clusters (Fig 4). In all segments, A. chejuensis-borne HTNV from Jeju Island formed homologous genetic clades with CUH15-126 obtained from the patients with HFRS in Jangheung-gun (Jeollanam Province) and distinguished from other A. agrarius-borne HTNV (Gyeonggi and Gangwon Provinces). Tanglegram analysis showed the segregation of representative rodentborne hantaviruses according to the species of their reservoir hosts (Fig 5). The phylogenetic positions of HTNV from A. chejuensis and other hantaviruses mirrored the co-evolutionary relationships of their Rodentia hosts, except for Tula and Lúxī virus.

Genetic reassortment analysis of HTNV
Recombination Detection Program 4 (RDP4) and Graph incompatibility based Reassortment Finder (GiRaF) analyses were performed to evaluate the possibility of genetic reassortment events. The genetic exchange event was undetectable among HTNV genomic sequences in this study.

Discussion
Hantaviruses have specifically evolved along with one or a few closely related rodent species [4,5,24]. The genus Apodemus (family Muridae, order Rodentia) consists of at least 20 small rodent species and is mainly widespread in the Palearctic area [25]. Molecular data based on the rodent mitochondrial DNA cytochrome b gene and the interphotoreceptor retinoid binding protein gene in nuclear DNA revealed that Apodemus spp. diverged approximately 8-10 million years ago [26]. In Asia, A. agrarius was first identified as the primary rodent reservoir of HTNV [8]. A. peninsulae was also a rodent reservoir for SOOV (closely related to AMRV) distributed in ROK, China, and Russia [10,27,28]. The striped field mouse (A. agrarius coreae), the most abundant rodent species, are widely distributed in the Korean peninsula, whereas A. agrarius chejuensis is mainly found on Jeju Island, the southernmost island of ROK [29].     The present study proposes that A. chejuensis-borne HTNV is a genetically distinct genotype due to the clear natural host species classification and co-divergent evolution of HTNV and Apodemus spp. after the geographic segregation between Jeju Island and the Korean peninsula. The species demarcation criterion of the International Committee for Taxonomy of Viruses (ICTV) is currently a 7% difference in amino acid sequences of G N /G C and N proteins. The differences in amino acid sequences of the RdRp, glycoproteins, and N proteins of HTNV strains harbored by A. chejuensis and A. agrarius did not exceed 5%. Taken together, A. chejuensis-harbored HTNV is considered Hantaan orthohantavirus species, including A. agrarius-borne HTNV and A. peninsulae-borne SOOV (closely related to AMRV). Continued surveillance and large-scale investigations should be conducted to understand the evolutionary diversification and geographic distribution of orthohantavirus and their natural reservoir hosts in ROK. Rodent-borne orthohantaviruses pose a critical public health threat worldwide. In the ROK, approximately 400 HFRS cases are reported annually, with 33 mortality cases over a decade [7]. Approximately 70% HFRS cases are caused by A. agrarius-borne HTNV and the remaining by SEOV carried by R. norvegicus or other unidentified agents. Jeju Island includes urban and rural regions with a population of approximately 700,000. It is a popular tourism destination for international tourists to be listed as a Biosphere Reserve (2002), World Natural Heritage Site (2007), and Global Geopark (2010). Overall, visitation to the site has been sustainably increased by an influx of international tourists due to global reputations [38]. The island has been an endemic HFRS region since the first clinical case was officially identified in 2005. After five HFRS cases were recorded from 2001 to 2009, 18 clinical cases have been reported in the recent decade on Jeju Island [7]. Despite the distribution of the rodent reservoir population, HTNV and SEOV have not yet been detected on Jeju Island. In the present study, the genetically distinct HTNV was detected in A. chejuensis captured on Jeju Island. Phylogenetic analyses of orthohantaviruses delineated a clear molecular epidemiological link between the patients with HFRS and their reservoir hosts collected at the infection sites [39]. In particular, active surveillance with targeted rodent trapping in the sites suspected of orthohantavirus outbreaks associated with HFRS allowed for phylogenetic distinction within a relatively short distance (about 5 km) to track the precise infectious agents and potential exposed sites [14]. Our previous studies demonstrated a phylogenetic relationship of partial HTNV M segment genome between four US soldiers diagnosed with HFRS and HTNV-positive A. agrarius collected at the very likely infection sites [40]. The phylogenetic analysis showed that HTNV derived from A. chejuensis genetically linked CUH15-126, which is the only hantaviral sequence collected from the patient with HFRS in southern ROK (Jeollanam Province) [41]. These results indicate that A. chejugensis-borne HTNV may be a potential etiological agent in clinical HFRS cases in southern ROK.
The subdivision of DOBV genotypes (Dobrava, Kurkino, Saaremaa, and Sochi) was differentiated by the characteristics in their phylogeny, rodent reservoirs, and geographical distribution [42]. Despite their high genetic homology, different DOBV genotypes induce HFRS with varied severities. The case-fatality rate (CFR) of clinical cases due to the Dobrava genotype harbored by A. flavicollis was 10-12% [43][44][45]. The CFR of HFRS cases caused by the Sochi genotype harbored by A. ponticus was approximately 6% [46]. Furthermore, the CFR range for orthohantavirus outbreaks associated with the Kurkino genotype carried by A. agrarius was 0.3-0.9% in European Russia. A. agrarius-borne SAAV genotype infection appears to be subclinical [47,48]. Thus, the transmission, pathogenicity, and disease risk assessment of HFRS by distinct A. chejugensis-borne HTNV in humans should be determined.
Genomic reassortment plays an important role in the evolutionary mechanism by which segmented RNA viruses may shuffle genome segments to generate new viral strains [49]. Genetic exchange of hantaviruses occurs more frequently between genetically closed strains than between distant strains. The incongruence of phylogenetic patterns in the tripartite segments showed the possibility of genetic reassortment among hantaviruses [46,50]. In this study, the phylogenies of A. chejuensis-borne HTNV L, M, and S segments revealed differing levels of incongruence in phylogenetic positions, which indicates the differential evolution of each segment. The L segment of the HTNV shared common ancestors with other HTNV lineages collected from mainland ROK and China. The M segment formed a clade with HTNV strains in the ROK. The S segment showed the independent geographic and distinct cluster with Hantaan orthohantavirus, including HTNV, SOOV, and AMRV. The phylogenetic pattern of A. chejuensis-borne HTNV may be a configuration compatible with genetic reassortment. However, recombination and reassortment were considered insignificant by RDP4 and GiRaF, respectively. To better understand the complexity of orthohantavirus evolution, future studies should focus on continuous sample collection from reservoir hosts and analysis of hantaviral sequences.
Sexual maturity and aggression are associated with viral transmission among natural host reservoirs [51,52]. The prevalence of orthohantaviruses in rodent populations plays a critical role in understanding hantaviral disease in HFRS-endemic areas. The seroprevalence of orthohantaviruses in rodents demonstrated a higher prevalence of infection in males than females [16,[53][54][55][56][57]. Hinson et al. reported that wounded male rats were more frequently infected with SEOV than non-wounded male rats [58]. In this study, the prevalence of HTNV antibodies and RNA was higher in males (7/35, 20%; 5/7, 71.4%, respectively) than those in females (2/29, 6.9%; 1/2, 50%, respectively). Our results showed a gender-specific prevalence of HTNVinfected A. chejuensis, supporting that males may contribute to the spread of viral infection via competition with other rodents. However, the number of samples (n = 64) was limited in this study. Thus, continuous large-scale epidemiological surveillance should be conducted to clarify the pattern of male-preferential prevalence of HTNV-infected A. chejuensis.
In conclusion, a genetically distinct HTNV was found in A. chejuensis captured on Jeju Island, ROK. The phylogeny of orthohantaviruses and their rodent hosts demonstrated that A. chejuensis may be infected with ancestral HTNV prior to geological isolation. Our results suggest that A. chejuensis-borne HTNV may be a potential etiological agent of HFRS in southern ROK. The discovery of a novel HTNV genotype provides important insights into the evolutionary history, genetic diversity, gender-specific prevalence of orthohantaviruses and a clue to account for HFRS cases that cannot be attributed to A. agrarius-borne HTNV infection in the ROK.

Ethics statement
Small mammals were handled in accordance with the Ethical Guidelines of the Korea University Institutional Animal Care and Use Committee (KUIACUC #2019-4 and 2019-171). Captured rodents and shrews were euthanized via cardiac puncture under alfaxalone-xylazine anesthesia. The experiment was performed in an animal biosafety level 3 facility at Korea University.

Sample collection
Small mammals were captured using Sherman live traps (8 × 9 × 23 cm; H. B. Sherman; Tallahassee, FL, USA) on Jeju Island, ROK from 2018 to 2020. A total of 100 traps were placed in unmanaged grasses, herbaceous vegetation, farmlands, mountains, fields, and forests for 1-3 days (S5 Table). The captured small mammals were identified morphologically and transported to Korea University in accordance with the safety guidelines. Multiple tissues of the rodents and shrews were dissected aseptically and the sera were separated by centrifugation for 5 min at 4˚C. The tissue samples were stored at −80˚C until use.

Mitochondrial DNA (mtDNA) analysis
Total DNA was extracted from liver tissues using a High Pure PCR Template Preparation Kit (Roche; Basel, Switzerland). The cytochrome b gene of mtDNA was amplified using universal primers [59]. The sequences were deposited in GenBank (accession numbers: MW219774-MW2197749).

Indirect immunofluorescence antibody (IFA) test
The sera and heart fluids of small mammals were used to detect anti-HTNV immunoglobulin G (IgG). The sera or heart fluids were diluted 1:32 and 1:2, respectively, in phosphate buffered saline (PBS) and added to acetone-fixed Vero E6 cells infected with HTNV. The slide was incubated at 37˚C for 30 min. The cells were washed with PBS and distilled water (D.W.). The slides were treated with fluorescein isothiocyanate-conjugated anti-mouse (for A. chejuensis and T. triton) or mouse/rat (for C. shantungensis) IgG (ICN Pharmaceuticals; Laval, Quebec, Canada). The plates were incubated at 37˚C for 30 min and washed with PBS and D.W. The virus-specific fluorescence was evaluated using a fluorescent microscope (Axio Scope; Zeiss; Berlin, Germany).

Reverse transcription-polymerase chain reaction (RT-PCR)
Total RNA was extracted from lung tissues of rodents and shrews with TRI Reagent Solution (Ambion; Austin, Texas, USA). cDNA was synthesized from 1 μg of total RNA using the High Capacity RNA-to-cDNA kit (Applied Biosystems; Foster City, CA, USA) with OSM55 (5 0 -TAG TAG TAG ACT CC-3 0 ). The PCR conditions were previously described [18]. All oligonucleotide primer sequences used in this study are shown in S6 Table. Quantitative PCR (qPCR) qPCR was performed in a 10 μL reaction mixture containing 1 μg of total RNA. The cycling conditions consisted of denaturation at 95˚C for 10 min, followed by 40 cycles at 95˚C for 15 s, and 60˚C for 1 min using SYBR Green PCR Master Mix (Applied Biosystems) on a QuantStudio 5 Real-Time PCR System (Applied Biosystems).

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 instructions. The composition consisted of 12.5 μL 2× Uh pre-mix, 2.0 μL each primer mixture, 10 μL cDNA template, and 10.5 μL D.W. in 25 μL reaction mixture. The PCR cycling conditions were previously described [18]. The PCR amplicons were prepared using the TruSeq Nano DNA LT sample preparation kit (Illumina; San Diego, USA) according to the manufacturer's instructions and previously described study [60]. NGS was performed on a MiSeq benchtop sequencer (Illumina) with 2 × 150 bp using a MiSeq reagent kit v2 (Illumina). The output files were imported and analyzed by CLC Genomics Workbench version 7.5.2 (CLC Bio; Cambridge, MA, USA). The sequences of HTNV strains were deposited in GenBank (accession numbers: MW219756-MW219773).

Rapid amplification of cDNA ends (RACE) PCR
To obtain the termini 3 0 and 5 0 genome sequences of HTNV, RACE PCR was conducted using a 3 0 -and 5 0 -RACE System for Rapid Amplification of cDNA Ends, Version 2.0 (Invitrogen; Carlsbad, CA, USA), according to the manufacturer's specifications.

Phylogenetic analysis
The tripartite genomic sequences of HTNV were aligned using the Clustal W algorithm in Lasergene version 5 (DNASTAR Inc.; Madison, WI, USA). For nearly full-length HTNV, phylogenetic analysis was performed using the best-fit GTR+G+I (for L and M segments) and GTR+G (for S segment) substitution models of evolution by the maximum-likelihood (ML) method in MEGA 7. Partial phylogenetic trees were generated by the ML method under the optimal evolutionary models T92+I (for L segment) and T92+G (for M and S segments). The topologies were assessed by bootstrap analysis for 1,000 iterations.

Co-divergence analysis
A tanglegram algorithm was generated to compare different phylogenetic links matching between HTNV strains and their hosts. The auxiliary lines in the center connect between the phylogenetic trees. The complexity between dendrograms of phylogenies was reduced to the maximum before the full reconciliation analysis. The method was implemented in the R package (dendextend) [61].

Genetic reassortment analysis
Alignment analysis of the concatenated HTNV L, M, and S segment open reading frame regions were performed using RDP, GENECONV, MAXCHI, CHIMAERA, 3SEQ, BOOTSCAN, and SISCAN in the RDP4 software [62]. A p-value of recombination and reassortment events under 0.05 was considered statistically significant. Genetic events were likely to indicate the occurrence of a genetic exchange when p-value was under 0.05, and the RDP recombination consensus score (RDPRCS) was between 0.4 and 0.6. All parameters of RDP4 were left at the default.
The GiRaF was performed to confirm reassortment events [63]. Nucleotide alignments of the tripartite segments were used as an input source for MrBayes [64]. The best nucleotide substitution models were determined by using MEGA 7. 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 default settings. The process was repeated ten times, with ten independent MrBayes based tree data per segment.

N-glycosylation site analysis
For N-glycosylation site prediction, the full-length amino acid sequences were submitted to the NetNGlyc 1.0 server (https://services.healthtech.dtu.dk) and the parameter was set to default [65].  119-1,080 nt). The phylogenetic tree was generated by the maximum likelihood method using MEGA 7.0. A bootstrap support value >70% is described at the nodes.

Supporting information
Apodemus mice are color-coded according to geographic regions (red, mainland ROK; blue, Jeju Island of ROK; green, China). (TIF)