Discovery and Genetic Characterization of Novel Paramyxoviruses Related to the Genus Henipavirus in Crocidura Species in the Republic of Korea

Paramyxoviruses, negative-sense single-stranded RNA viruses, pose a critical threat to human public health. Currently, 78 species, 17 genera, and 4 subfamilies of paramyxoviruses are harbored by multiple natural reservoirs, including rodents, bats, birds, reptiles, and fish. Henipaviruses are critical zoonotic pathogens that cause severe acute respiratory distress and neurological diseases in humans. Using reverse transcription-polymerase chain reaction, 115 Crocidura species individuals were examined for the prevalence of paramyxovirus infections. Paramyxovirus RNA was observed in 26 (22.6%) shrews collected at five trapping sites, Republic of Korea. Herein, we report two genetically distinct novel paramyxoviruses (genus: Henipavirus): Gamak virus (GAKV) and Daeryong virus (DARV) isolated from C. lasiura and C. shantungensis, respectively. Two GAKVs and one DARV were nearly completely sequenced using next-generation sequencing. GAKV and DARV contain six genes (3′-N-P-M-F-G-L-5′) with genome sizes of 18,460 nucleotides and 19,471 nucleotides, respectively. The phylogenetic inference demonstrated that GAKV and DARV form independent genetic lineages of Henipavirus in Crocidura species. GAKV-infected human lung epithelial cells elicited the induction of type I/III interferons, interferon-stimulated genes, and proinflammatory cytokines. In conclusion, this study contributes further understandings of the molecular prevalence, genetic characteristics and diversity, and zoonotic potential of novel paramyxoviruses in shrews.

and 'Daeryong virus' from the first discovered location of origin (Mount Gamak and Daeryong, respectively). GAKV, DARV, and rodent-borne henipavirus (MojV) appear to be distinct from previously identified bat-borne henipaviruses (Hendra virus, NiV, Cedar virus, and Kumasi virus), suggesting that they may be novel species within the genus Henipavirus. Infection with GAKV elicited the expression of type I/III interferons (IFNs), IFN-stimulated genes (ISGs), and proinflammatory cytokines in A549 cells. Therefore, this report demonstrates the molecular prevalence, genetic characteristics and diversity, and virus-host interactions of paramyxoviruses in shrews, ROK.

Ethics Statement
The US Forces Korea (USFK) approved the animal trapping procedure by USFK Regulation 40-1 "Prevention, Surveillance, and Treatment of Hemorrhagic Fever with Renal Syndrome". The trapping, experiments, and handling of animals were conducted under a protocol approved by the Korea University Institutional Animal Care and Use Committee (KUIACUC, #2016-0049).

Animal Capture and RNA Extraction
Small mammals were captured during 2017-2018 using Sherman traps (8 cm × 9 cm × 23 cm; H. B. Sherman, Tallahassee, FL, USA). A total of 115 individuals belonging to Crocidura species were collected from the previous study [30]. Shrews were identified using morphological criteria and polymerase chain reaction (PCR), when required. Serum, brain, lung, spleen, kidney, and liver tissues were collected aseptically and frozen at −80 • C until use. Total RNA was extracted from the kidney tissues of shrews using TRI Reagent Solution (AMBION Inc., Austin, TX, USA) according to the manufacturer's instructions.

Metagenomic Next-Generation Sequencing (NGS) Using Sequence-Independent, Single-Primer Amplification (SISPA)
cDNA was generated from the RNA extracted from shrew-borne paramyxovirusinfected cells and kidney tissues using SISPA method was carried out according to a protocol described previously [30].

NGS for Illumina MiSeq
The cDNA library was prepared using the TruSeq Nano DNA LT Sample Preparation Kit (Illumina, San Diego, CA, USA), and the amplicon was size-selected, A-tailed, ligated with indexes and adaptors, and enriched. The library was verified for quality and size with Agilent 2100 Bioanalyzer (Agilent High Sensitivity DNA kit, Agilent Technologies) and sequenced with MiSeq reagent V2 (Illumina) using MiSeq benchtop sequencer (Illumina) with 2 × 150 bp.

NGS for Illumina HiSeq
Libraries were prepared by using the NEBNext Ultra II Directional RNA-Seq Kit (New England Biolabs, Ipswich, MA, USA). Additionally, mRNA was isolated using the Poly(A) RNA Selection Kit (Lexogen, Vienna, Austria). The mRNAs were synthesized for cDNA and sheared according to the manufacturer's instructions. The fragments were size-selected, A-tailed, ligated with indexes and adaptors, and enriched by PCR. The libraries were evaluated for the mean fragment size by the Agilent 2100 bioanalyzer (DNA High Sensitivity Kit). Quantification was performed by using a library quantification kit and StepOne Real-Time PCR System (Life Technologies). Libraries were conducted as paired-end 100 sequencing by the HiSeq X10 system (Illumina).

Extraction of Paramyxoviral Genome Sequences from NGS Data
The reads were trimmed with Trimmomatic (v.0.36) to remove adapter sequences for the metagenomic approach [44]. Aligned reads against the host genomic sequences were excluded using Bowtie2 (v.2.2.6) [45]. The complete mitochondrial sequence of the species on the National Center for Biotechnology Information (NCBI) RefSeq was utilized as a host reference [46]. The remaining reads were filtered for quality by FaQCs (v0.11.5), and de novo assembly was conducted by using SPAdes (v3.11.1) [47,48]. The assembled contigs were subsequently analyzed in the NCBI RefSeq database consisting of complete viral genomes (updated in June 2019) by nucleotide-Basic Local Alignment Search Tool (BLASTn) (v2.6.0). In addition, for the reference mapping strategy, adaptor and index sequences of reads were trimmed, and low-quality sequences were filtered using the CLC Genomics Workbench v.7.5.2 (CLC Bio, Cambridge, MA, USA). The genome sequences of Hendra virus (HeV), NiV, and MojV were used in the reference mapping method. The consensus genomic sequences of GAKV and DARV were determined by combining viral contigs extracted from CLC analysis and de novo assembly.

Reverse Transcription-Polymerase Chain Reaction (RT-PCR) Screening of Samples from Crocidura Species Individuals for Paramyxovirus
RNA was used for cDNA synthesis using a high-capacity RNA-to-cDNA kit (Applied Biosystems, Foster City, CA, USA). First and nested PCRs were conducted in a 25 µL reaction mixture containing 2.5 U of Ex Taq DNA polymerase (TaKaRa BIO Inc., Shiga, Japan), 1.5 µg of cDNA, and 10 pM of pan-paramyxoviral primers [49]. PCR condition and purification were carried out according to a protocol described previously [30]. Sequencing was performed in both 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).

Rapid Amplification of cDNA Ends (RACE) PCR
The 3 and 5 terminal genome sequences of paramyxovirus were determined by using a SMARTer ® RACE 5 /3 Kit (Takara Bio), according to the manufacturer's specifications. PCR products were purified by the LaboPass PCR Purification Kit (Cosmo Genetech). Sequencing was conducted in both directions of each PCR product using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) on an automated sequencer (Applied Biosystems).

Mitochondrial DNA (mtDNA) Analysis
DNA was extracted from liver or lung tissues using a High Pure PCR template preparation kit (Roche, Basel, Switzerland). The cytochrome b gene of mtDNA was identified by universal primers [50].

Phylogenetic Analysis
The viral genomic sequences were aligned with MUSCLE algorithm in MEGA v.7.0 [51] and trimmed using the Lasergene program v.5 (DNASTAR, Madison, WI, USA). Phylogenetic analyses were conducted by the maximum likelihood method according to the best-fit substitution model (GTR+G+I) using MEGA 7.0. Support for the topologies was evaluated by bootstrapping for 1000 iterations. Additionally, the Markov chain Monte Carlo (MCMC) method and the BEAST package (v.1.10.4) as the Bayesian inference method were implemented to infer Bayesian phylogenetic trees [52]. The MCMC chain length was established to 100 million states by sampling every 50,000 states. The results of parameters produced sufficient sample sizes (ESS > 200). Maximum clade credibility trees were generated using TreeAnnotator (v.1.10.4) and FigTree (v.1.4.0).

Virus-Host Co-Divergence Analysis
To compare different phylogenetic links matching between paramyxoviruses and their hosts, a tanglegram algorithm was performed. The auxiliary lines in the center connect the phylogenetic trees. The complexity between dendrograms of phylogenies was diminished to a maximum before complete reconciliation analysis. The method was implemented using the R package "dendextend" [53].

Virus Isolation
Tissue homogenates of kidneys from GAKV-positive samples were inoculated on Vero E6 cells. The inoculum was removed after one and a half hours of adsorption. The culture was maintained with 5.5 mL of DMEM containing 5% FBS at 37 • C in a 5% CO 2 incubator. Cytopathic effects were inspected daily using an inverted microscope.

Plaque Assay
Vero E6 cells were prepared at a density of 1.5 × 10 6 cells per well onto 6-well plates. The cells were washed twice with phosphate-buffered saline (PBS) and inoculated serially with 10-fold diluted viruses. A 1:1 medium-melting-point agarose and overlay medium mix was added on the cells after 90 m of incubation at 37 • C with constant shaking. The agarose overlay was discarded at day 6. The plaques were stained with 0.1% crystal violet in 10% formaldehyde for visualization.

Electron Microscopy
GAKV-inoculated Vero E6 cells were collected at 5 d postinfection and fixed with 2% paraformaldehyde and 2.5% glutaraldehyde in 0.1 M PBS (pH 7.4). Ultrathin sections were placed on 400-mesh square copper electron microscopy grids (Electron Microscopy Sciences) and observed by using a transmission electron microscope (Model H-7650; Hitachi, Japan).

In Vitro Infection
A total of 1 × 10 6 cells were infected with GAKV at a multiplicity of infection of 0.02. The samples were then collected at 2, 6, 12, 18, and 24 h postinfection (hpi). The mRNA genes of type I/III IFNs, ISGs, and proinflammatory cytokines were examined by using the primers (Table S1).

Statistical Analysis
Statistical analyses were conducted as shown in figures using the GraphPad Prism v.5.00 for Windows (GraphPad Software, San Diego, CA, USA; www.graphpad.com (accessed on 1 September 2021)).

Isolation of Novel Henipaviruses in Shrews
GAKV was isolated from the kidney tissues of Cl17-32 using a cell culture-based method. The first isolate of GAKV was confirmed by passaging twice for 14 days postinoculation. GAKV particles were observed using a transmission electron microscope ( Figure 2A). In addition, the number of infectious GAKV particles was 1 × 10 6 PFU/mL, as quantified using the plaque assay ( Figure 2B).

Nearly Whole-Genome Sequencing of Shrew-Borne Henipaviruses Using NGS
We performed host rRNA depletion on SISPA-based MiSeq of the kidney tissue of Cl17-46, resulting in two contigs (602 and 17,892 nt in length) with significant similarities in the genomic sequence of paramyxoviruses. The NGS of Cl17-46 generated 245,846 reads, and the number of viral genome sequences was 3608 (Table S5). To obtain whole-genome sequences of shrew-borne paramyxovirus, RNA-Seq (HiSeq) was performed on samples from the kidney tissues of Cl17-46 and Cs17-65 and on cells infected with the paramyxovirus isolated from Cl17-32. The genomic sequences of henipaviruses were identified across 23 contigs (408-18,382 nt in length) with significant similarities. The number of viral genome reads was 14,779,024, 6497, and 50,987 from Cl17-32, Cl17-46, and Cs17-65, respectively (Table S6). Nearly complete genome sequences of GAKV and DARV were obtained except for both 3 and 5 terminal ends. The genomic sequences of GAKV and DARV have been deposited in GenBank (accession number: MZ574407-MZ574409).

Nearly Whole-Genome Sequencing of Shrew-Borne Henipaviruses Using NGS
We performed host rRNA depletion on SISPA-based MiSeq of the kidney tissue of Cl17-46, resulting in two contigs (602 and 17,892 nt in length) with significant similarities in the genomic sequence of paramyxoviruses. The NGS of Cl17-46 generated 245,846 reads, and the number of viral genome sequences was 3608 (Table S5). To obtain wholegenome sequences of shrew-borne paramyxovirus, RNA-Seq (HiSeq) was performed on samples from the kidney tissues of Cl17-46 and Cs17-65 and on cells infected with the paramyxovirus isolated from Cl17-32. The genomic sequences of henipaviruses were identified across 23 contigs (408-18,382 nt in length) with significant similarities. The number of viral genome reads was 14,779,024, 6497, and 50,987 from Cl17-32, Cl17-46, and Cs17-65, respectively (Table S6). Nearly complete genome sequences of GAKV and DARV were obtained except for both 3′ and 5′ terminal ends. The genomic sequences of GAKV and DARV have been deposited in GenBank (accession number: MZ574407-MZ574409).  The N, M, F, G, and L genes encode a single protein, whereas the P gene encodes multiple accessory proteins: the viral phosphoprotein, C protein, and V/W protein. GAKV and DARV possessed a putative RNA editing site (TTAAAAAAGGCA) of the P gene, a conserved motif sequence (YTAAAARRGGCA) in the members of the genus Henipavirus. The 3′ leader and 5′ trailer sequences were a length of 52 nt and 25 nt, respectively. Tables S7 and S8 show the start, stop, and intergenic region (IGR) sequences of GAKV and DARV in shrews.

Phylogenetic Analysis
The phylogenetic analysis of all 26 currently recognized Orthoparamyxovirus species elucidated that GAKV and DARV belong to the genus Henipavirus and are closely related to the henipaviruses (Figure 4). In addition, the amino acid sequences (N, P, M, F, G, and L) of GAKV and DARV showed the phylogenetic shape comparable to the viral RNA genome sequences among major paramyxoviruses ( Figure S1). The genetic cluster of

Phylogenetic Analysis
The phylogenetic analysis of all 26 currently recognized Orthoparamyxovirus species elucidated that GAKV and DARV belong to the genus Henipavirus and are closely related to the henipaviruses (Figure 4). In addition, the amino acid sequences (N, P, M, F, G, and L) of GAKV and DARV showed the phylogenetic shape comparable to the viral RNA genome sequences among major paramyxoviruses ( Figure S1). The genetic cluster of GAKV and DARV shared a common ancestor with MojV with an amino acid, similarity to the L gene of 62.9% and 74.1%, respectively (Table S9).
GAKV and DARV shared a common ancestor with MojV with an amino acid, similarity to the L gene of 62.9% and 74.1%, respectively (Table S9).

Co-Phylogeny of Paramyxovirus and Host
Cophylogeny mapping demonstrated the segregation of paramyxoviruses according to the subfamily of their reservoir hosts, using consensus trees based on the L protein amino acid sequences (

Co-Phylogeny of Paramyxovirus and Host
Cophylogeny mapping demonstrated the segregation of paramyxoviruses according to the subfamily of their reservoir hosts, using consensus trees based on the L protein amino acid sequences ( Figure 5). The phylogenetic positions of GAKV, DARV, and shrew paramyxovirus (Shrew PV) mirrored the phylogenetic relationships of the shrew species. In contrast, the phylogenetic positions of viruses from rodents such as Bank Vole virus (BaVV) and Pohorje Myodes paramyxovirus 1 (PMPV-1) from Myodes glareolus, Mount Mabu Lophuromys virus 1 (MMLV-1) and Mount Mabu Lophuromys virus 2 (MMLV-2) from Lophuromys machangui, and PAPV-1 and PAPV-2 from Apodemus agrarius, were not directly matched, as per the tanglegram between the virus and host species.

Discussion
We characterized two novel paramyxoviruses, GAKV and DARV, identified in Crocidura species, ROK. GAKV and DARV showed a high positive rate: GAKV was found in five regions, while DARV was detected in one region. The genomic structure (3′-N-P-M-F-G-L-5′) of these viruses is representative of the genus of henipaviruses, including HeV, NiV, MojV, CedV, and KV. The phylogenies of GAKV and DARV shaped independent

Discussion
We characterized two novel paramyxoviruses, GAKV and DARV, identified in Crocidura species, ROK. GAKV and DARV showed a high positive rate: GAKV was found in five regions, while DARV was detected in one region. The genomic structure (3 -N-P-M-F-G-L-5 ) of these viruses is representative of the genus of henipaviruses, including HeV, NiV, MojV, CedV, and KV. The phylogenies of GAKV and DARV shaped independent genetic lineages corresponding to the paramyxovirus species distinctive criterion (the L gene amino acid distance of <0.82 for the Orthoparamyxovirinae) [26]. These observations suggest that GAKV and DARV are genetically distinct novel shrew-borne paramyxoviruses within the genus Henipavirus in the family Paramyxoviridae.
Cophylogenetic relationships of zoonotic viruses and natural reservoir hosts play a critical role in understanding the evolutionary process for coadaptation, spillover, and host sharing [60][61][62]. Multiple natural reservoirs harbor Orthohantavirus species, and in most cases, it is a specific relationship of "one virus-one host" [63]. The congruent relationships between whole groups of hantaviruses and (sub)families of hosts proposed a view of the long-term coevolution of these viruses with their hosts [64,65]. However, recent studies, based on cophylogenetic reconciliation and estimation of evolutionary rates and divergence times, demonstrate that local host-specific adaptation and preferential host switching account for the phylogenetic similarities between viruses and their mammalian hosts [66,67]. Our cophylogenetic analysis indicates that shrew-borne henipaviruses and their host species have coadapted together differently from the phylogenies of rodent-borne paramyxoviruses. Thus, continuous and large-scale investigations will provide insights into the evolutionary history and dynamics of paramyxoviruses, harbored by shrews and rodents, in the family Orthoparamyxoviridae.
Triggering of innate antiviral genes modulates viral infectious diseases in humans and animals [68,69]. The expression of antiviral genes, Mx1, Rsad2/Viperin, Isg15, and OAS1, was highly induced during early NiV and HeV infections in ferrets [70]. These robust innate antiviral responses were insufficient to prevent viral dissemination and tissue damage in vivo [71]. The isolate of GAKV from C. lasiura was evaluated for the innate antiviral response in human lung epithelial cells. We found that GAKV enabled replication and elicited the rapid induction of type I/III IFNs, ISGs, and pro-inflammatory cytokines in A549 cells. These results demonstrated that the novel paramyxovirus may infect and induce pro-inflammatory responses via the human respiratory tract. The characterization and effect of innate antiviral responses to shrew-borne henipaviruses in humans await for further studies.
In conclusion, this study reports two generically distinct novel paramyxoviruses, GAKV in C. lasiura and DARV in C. shantungensis. Phylogenetic inference and genomic characterization demonstrated that these viruses belong to henipaviruses within the family Paramyxoviridae. Infection with GAKV showed the upregulation of multiple human innate antiviral genes in vitro. Therefore, these results facilitate understanding of the molecular epidemiology, genetic characteristics and diversity, and zoonotic potential of shrew-borne paramyxoviruses, ROK. These observations raise awareness and caution among virologists and physicians about the henipa-like paramyxoviruses in shrews.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/v13102020/s1, Figure S1: Phylogenetic trees using the coding sequences (CDSs) of N, P, M, F, G, and L genes of Gamak virus (GAKV) and Daeryong virus (DARV); Figure S2: Phylogenetic analyses of the paramyxoviruses using partial genome sequences of L gene; Table S1: Sequences of primers for quantitative reverse transcription-polymerase chain reaction; Table S2: Characteristics of Crocidura species infected with Gamak virus (GAKV) and Daeryong virus (DARV) and the nucleotide position of GAKV and DARV RNA obtained in this study; Table S3: Molecular prevalence of Gamak virus (GAKV) by region, sex, weight, and season in C. lasiura, the Republic of Korea from 2017 to 2018; Table S4: Molecular prevalence of Daeryong virus (DARV) by region, sex, weight, and season in C. shantungensis, the Republic of Korea from 2017 to 2018; Table S5: Genome coverage of consensus sequences for Gamak virus using next-generation sequencing (SISPA-based MiSeq) method; Table S6: Genome coverage of consensus sequences for Gamak and Daeryong viruses using nextgeneration sequencing (RNA-Seq based HiSeq) method; Table S7: Sequences of intergenic regions (IGRs) and transcriptional stop and start signals of Gamak virus; Table S8: Sequences of intergenic regions (IGRs) and transcriptional stop and start signals of Daeryong virus; Table S9: Amino acid similarities of GAKV and DARV with other viruses in Paramyxoviridae.

Informed Consent Statement: Not applicable.
Data Availability Statement: All the data generated for this publication have been included in the current manuscript.