Epidemiology, Genetic Characterization, and Evolution of Hunnivirus Carried by Rattus norvegicus and Rattus tanezumi: The First Epidemiological Evidence from Southern China

Hunnivirus is a novel member of the family Picornaviridae. A single species, Hunnivirus A, is currently described. However, there is limited information on the identification of Hunnivirus to date, and thereby the circulation of Hunnivirus is not fully understood. Thus, the objective of this study was to investigate the prevalence, genomic characteristics, and evolution of rat hunnivirus in southern China. A total of 404 fecal samples were subjected to detection of Hunnivirus from urban rats (Rattus norvegicus and Rattus tanezumi) using PCR assay based on specific primers targeted to partial 3D regions, with the prevalence of 17.8% in Rattus norvegicus and 15.6% in Rattus tanezumi. An almost full-length rat hunnivirus sequence (RatHuV/YY12/CHN) and the genome structure were acquired in the present study. Phylogenetic analysis of the P1 coding regions suggested the RatHuV/YY12/CHN sequence was found to be within the genotype of Hunnivirus A4. The negative selection was further identified based on analysis of non-synonymous to synonymous substitution rates. The present findings suggest that hunniviruses are common in urban rats. Further research is needed for increased surveillance and awareness of potential risks to human health.

In 1965, the identification of the first hunnivirus was detected from sheep cell cultures in Northern Ireland. Subsequently, bovine and ovine hunniviruses were determined in fecal samples of cattle and sheep in central Hungary from 2008 to 2009 [1]. Hunnivirus was then proposed as a member of a novel genus in the family Picornaviridae. At the request of the ICTV Executive Committee, the name was modified from 'hungarovirus' in July 2013 to the name hunnivirus, derived from these two countries, Hungary and Northern Ireland (https://www.picornaviridae.com/hunnivirus/hunnivirus.htm, accessed on 20 May 2021). Between 2014 and 2016, the hunniviruses were identified in commensal Norway rats (Rattus norvegicus) in New York City [2] and detected in rodents from several representative regions of China [3]. A novel hunnivirus, provisionally designated feline hunnivirus (FeHuV), has been recently discovered in a cat with diarrhea in Guangdong, China, indicating that pets have the potential to spread members of this genus [4].
In the last few years, the hunniviruses have gradually been identified worldwide in different animals, including cattle, sheep, rats, and cats [1,2,4], which raises interesting questions regarding the global distribution, pathogenesis, and potential risks to human health. However, little is known about these questions, which are thereby challenging to answer. Herein, this study investigated the prevalence, genetic characteristics, and evolution of Hunnivirus carried by Rattus norvegicus and Rattus tanezumi in southern China.

Phylogenetic Analysis of Partial 3CD Region Sequences
Twelve representative sequences of the partial 3CD regions (1550 nt) were randomly selected for further PCR amplification and sequencing within the seventy-one positive samples (GenBank accession numbers MW417230-MW417241). Our sequences that were isolated from two different rodent species captured in different geographic locations shared 94.6-99.3% nucleotide identities and 97.2-99.6% amino acid identities with one another. Besides, they displayed the closest relatedness (94.0-95.8% nucleotide and 97.0-99.2% amino acid identities) with a Norway rat hunnivirus sequence (rat/NYC-E21/USA/2012) previously reported in the New York City (accession no. KJ950971.1) [2]. Maximum-likelihood phylogenetic analysis based on the 1550-nt 3CD nucleotide sequences was carried out using other hunnivirus reference sequences in GenBank. The result demonstrates that our sequences clustered along with the rat hunniviruses previously found in Asia and America and formed a tight cluster in a monophyletic branch with a hunnivirus sequence Wencheng-Rt38-2 isolated from Rattus tanezumi. These sequences also shared a common root with the hunnivirus sequences identified in felines and rats, but were related more distantly to bovine and ovine hunniviruses ( Figure 1). Maximum-likelihood phylogenetic analysis based on the 1550-nt 3CD nucleotide sequences was carried out using other hunnivirus reference sequences in GenBank. The result demonstrates that our sequences clustered along with the rat hunniviruses previously found in Asia and America and formed a tight cluster in a monophyletic branch with a hunnivirus sequence Wencheng-Rt38-2 isolated from Rattus tanezumi. These sequences also shared a common root with the hunnivirus sequences identified in felines and rats, but were related more distantly to bovine and ovine hunniviruses ( Figure 1).

Figure 1.
Phylogenetic tree of Hunniviruses based on the partial 3CD nucleotide sequences. The tree was generated by the maximum-likelihood method based on the General Time Reversible (GTR) model (gamma-distributed with invariant sites (G + I) and partial deletion) with 1000 bootstrap replicates, and the statistics values > 70% are displayed above the tree branches. Red font represents the sequences of rat hunniviruses identified in the present study.

Genomic Analyses of the Near-Complete Sequence
An almost full-length rat hunnivirus sequence (RatHuV/YY12/CHN, 7282 nt) was successfully obtained from Rattus norvegicus (accession no. MW417242). This sequence included a partial 5′ untranslated region (UTR) of 481 nt, one complete ORF of 6705 nt that encoded a potential polyprotein of 2234 amino acids, and a partial 3′ UTR of 96 nt. Further, we predicted the genome organization and potential cleavage sites for our near-complete sequence. It was found that the L protein was 210 nt (70 aa) in length, whereas the complete P1, P2, and P3 regions were 2337 (779 aa), 1827 (609 aa), and 2331 nt (776 aa) long, respectively. The detailed genome structure of RatHuV/YY12/CHN is shown in Figure  2A. In addition, we presented the predicted RNA secondary structure of 5′UTR by RNAfold, evidencing that the 481-nt 5′UTR sequence can construct a stable secondary structure with the minimum free energy of −164.80kcal/mol ( Figure 3). Phylogenetic tree of Hunniviruses based on the partial 3CD nucleotide sequences. The tree was generated by the maximum-likelihood method based on the General Time Reversible (GTR) model (gamma-distributed with invariant sites (G + I) and partial deletion) with 1000 bootstrap replicates, and the statistics values > 70% are displayed above the tree branches. Red font represents the sequences of rat hunniviruses identified in the present study.

Genomic Analyses of the Near-Complete Sequence
An almost full-length rat hunnivirus sequence (RatHuV/YY12/CHN, 7282 nt) was successfully obtained from Rattus norvegicus (accession no. MW417242). This sequence included a partial 5 untranslated region (UTR) of 481 nt, one complete ORF of 6705 nt that encoded a potential polyprotein of 2234 amino acids, and a partial 3 UTR of 96 nt. Further, we predicted the genome organization and potential cleavage sites for our nearcomplete sequence. It was found that the L protein was 210 nt (70 aa) in length, whereas the complete P1, P2, and P3 regions were 2337 (779 aa), 1827 (609 aa), and 2331 nt (776 aa) long, respectively. The detailed genome structure of RatHuV/YY12/CHN is shown in Figure 2A. In addition, we presented the predicted RNA secondary structure of 5 UTR by RNAfold, evidencing that the 481-nt 5 UTR sequence can construct a stable secondary structure with the minimum free energy of −164.80 kcal/mol ( Figure 3).     The findings of BLASTn analysis indicated that our complete polyprotein sequence shared the highest nucleotide and amino acid identities with the rat/NYC-E21/USA/2012 sequence, with values of 91.2 and 92.7%, respectively. In contrast to other rat-derived hunnivirus sequences, RatHuV/YY12/CHN shared 70.3-79.5% nucleotide identities and 71.5-79.1% putative amino acid identities with two Vietnamese sequences and one Chinese sequence. Interestingly, RatHuV/YY12/CHN had a higher sequence similarity with the FeHuV identified in Guangdong, China, at the nucleotide (85.1%) and putative amino acid (86.5%) levels. However, it shared a relatively low nucleotide/amino acid identity with bovine (67.6/70.1%) and ovine (67.7/69.9%) hunnivirus sequences. The similarities of each functional region were further compared between RatHuV/YY12/CHN and other reported hunniviruses in GenBank ( Table 2). We found the highest nucleotide (97.2%) and amino acid (97.4%) identities presented in the rat/NYC-E21/USA/2012 L and P2 regions, respectively. In comparison, the lowest nucleotide (63.0%) and amino acid identities (45.9%) occurred in the P2 and L regions of bovine hunnivirus sequence, respectively. Overall, RatHuV/YY12/CHN shared 61.2-92.9%, 63.5-93.1%, 70.4-97.4%, and 81.2-96.9% amino acid identities with the rat-derived hunniviruses at the L, P1, P2, and P3 regions, respectively. Taken together, the genetic diversity is existent in the hunniviruses from different host species. Additionally, a similarity plot analysis was conducted to further analyze the genetic characteristics of our near-complete nucleotide sequence, with the bovine hunnivirus sequence BHUV1/HUN/2008/JQ941880.1 used as the out-group sequence and the Norway rat hunnivirus sequence rat/NYC-E21/USA/2012/KJ950971.1 used as the query sequence. The findings indicate that RatHuV/YY12/CHN exhibited relatively high similarities to the query sequence in the 5 UTR, L, VP4, VP3, 2B-2C, and from the 3A to 3D regions, while different similarities were presented in the VP2, VP1, and 2A regions, and the sequence connecting 2B and 2C. Moreover, the P1 coding region was referred to as having a higher range of genetic variability ( Figure 2B).

Phylogenetic Analyses and Negative Selection during the Evolution of Hunnivirus
The maximum-likelihood phylogenetic tree was constructed to investigate the phylogenetic relationship between RatHuV/YY12/CHN and other Hunnivirus A strains based on the complete P1 coding regions. The analytical results demonstrated that RatHuV/YY12/ CHN clustered with the Chinese FeHuV sequence (accession no. MF953886.1), within the genotype of Hunnivirus A4 (Figure 4). These P1 sequences in the phylogenetic tree were then aligned for selective pressure analyses. Despite no positively selected sites, the strong negative selection on the non-synonymous sites was detected in SLAC, FEL, and REL methods [5], with the negatively selected sites of 319, 519, and 340, respectively. phylogenetic tree were then aligned for selective pressure analyses. Despite no positively selected sites, the strong negative selection on the non-synonymous sites was detected in SLAC, FEL, and REL methods [5], with the negatively selected sites of 319, 519, and 340, respectively. To further characterize the evolution of Hunnivirus and detect possible differences in selective pressure between phylogenetic branches, we obtained the annotated complete coding regions (CDSs) in the hunnivirus reference sequences in GenBank. Then, we reconstructed the phylogenetic tree using the synonymous sites ( Figure 5). The analytical findings revealed that our sequence RatHuV/YY12/CHN clustered into a single group with the Chinese feline and rat hunnivirus sequences. The mean dN/dS values were calculated using the GA Branch method for each branch, indicating that all the phylogenetic branches were under strong negative selection (0 < dN/dS < 1). To further characterize the evolution of Hunnivirus and detect possible differences in selective pressure between phylogenetic branches, we obtained the annotated complete coding regions (CDSs) in the hunnivirus reference sequences in GenBank. Then, we reconstructed the phylogenetic tree using the synonymous sites ( Figure 5). The analytical findings revealed that our sequence RatHuV/YY12/CHN clustered into a single group with the Chinese feline and rat hunnivirus sequences. The mean dN/dS values were calculated using the GA Branch method for each branch, indicating that all the phylogenetic branches were under strong negative selection (0 < dN/dS < 1).

Discussion
In 1965, the hunnivirus was first discovered in sheep cell cultures. However, it was not reported in cattle, sheep, rats, and cats until the last few years [1,2,4]. The detailed data related to hunnivirus in different host species and different geographic locations are quite limited. To the best of our knowledge, this study represents the first study to specifically

Discussion
In 1965, the hunnivirus was first discovered in sheep cell cultures. However, it was not reported in cattle, sheep, rats, and cats until the last few years [1,2,4]. The detailed data related to hunnivirus in different host species and different geographic locations are quite limited. To the best of our knowledge, this study represents the first study to specifically investigate the prevalence of hunniviruses in fecal samples from Rattus norvegicus and Rattus tanezumi in southern China. We investigated 404 fecal samples of which 71 (17.6%) were positive for hunnivirus using self-designed primers. The mean detection rates of Hunnivirus in Rattus norvegicus and Rattus tanezumi were 17.8 and 15.6%, respectively, similar to the prevalence of rat hunnivirus (16%, 21/133) in the USA and bovine hunnivirus 1 (15%, 4/26) in central Hungary [1,2]. Nevertheless, information on the global distribution of Hunniviruses remains to be determined due to the lack of epidemiological evidence. Further studies are warranted to investigate the geographic distribution and host species of Hunnivirus worldwide, thereby confirming variations in prevalence among different host species.
In previous studies, the hunniviruses were serendipitously identified when detecting other picornaviruses using the primers UNIV-kobu-F/UNIV-kobu-R. These primers were designed to screen three prototypes of kobuviruses (Aichi virus, bovine kobuvirus, and porcine kobuvirus) based on the 3D conserved viral regions [6]. To date, numerous kobuviruses have been determined worldwide among various host species using this primer pair, including human [7], cattle [8,9], pig [10], cat [11,12], dog [13], deer [14], and wolf [15]. Similarly, an epidemiological study of rat kobuvirus isolated from murine rodents in Guangdong, China, was previously reported by our research team [16]. However, this primer pair was considered more specific for picornaviruses than for kobuviruses [1,6], confirming the universal applicability in identifying novel picornaviruses [4]. For instance, the bovine hunnivirus was initially identified in a cattle sample but further found in three additional samples after using specific hungarovirus screening primers (Hungaro-3D-R/F) [1]. Likewise, our present study found that one 215-nt-long nucleotide sequence had no similarity to kobuviruses available in the GenBank database but shared 83.2-87.9% nucleotide identity with five hunnivirus sequences by BLAST.
To address this issue, we designed generic hunnivirus screening primers (RHuV-F/RHuV-R) based on the 3D conserved region to determine rat hunniviruses from additional fecal samples and consequently, expected PCR products of amplicon size (254 nt) for hunnivirus were identified. We confirmed the low sensitivity of the primers UNIV-kobu-F/UNIV-kobu-R that may reflect the detection rate of hunnivirus in rats [1]. Additionally, we also screened kobuvirus for hunnivirus-positive specimens, of which 15 specimens were co-infected with kobuvirus. This finding was consistent with a prior study demonstrating that kobuviruses frequently mixed infections with other pathogens [17]. Hence, further investigation should be carried out to determine whether the co-detection or multiple virus detection is equally common for hunniviruses and other pathogens.
Maximum-likelihood phylogenetic analysis based on the partial 3CD regions suggested that the nucleotide sequences of rat kobuvirus isolated from Rattus norvegicus and Rattus tanezumi clustered tightly together in a clade according to the species, but regardless of being derived from different geographical locations (Figure 1). Moreover, our sequences formed a group with other rat-derived hunnivirus sequences identified in China, Vietnam, and the USA. Interestingly, they also clustered with a Chinese FeHuV (accession no. MF953886.2), and the majority shared relatively high nucleotide (95.1-96.6%) and amino acid identities (97.0-99.2%) with the Chinese FeHuV sequence [4]. What is more, the phylogenetic relationships based on P1 regions and complete coding regions between our sequence and other hunnivirus sequences (Hunnivirus A1-A9) illustrated that our rat hunnivirus sequence shared a common root with the Chinese FeHuV sequence (Figures 4 and 5). The combined findings suggest that the cross-species transmission of hunnivirus possibly occurs in rats and other animal species, which would result in a sustained threat to public health due to the increased contact between humans and these animals. Since less attention has been focused on the potential risk of cross-species transmission for hunniviruses, it is of significance to conduct serological investigations and assess the pathogenicity posed to human health.
The almost complete gene sequence of RatHuV/YY12/CHN (7282nt) was successfully acquired in the present study. Except for the 5 and 3 UTR regions, the genetic analyses suggested that the polyprotein structure of RatHuV/YY12/CHN could be cleaved into 12 viral proteins, including L, P1 (VP4, VP2, VP3, and VP1), P2 (2A-2C), and P3 (3A-3D). The results of a comparative sequence analysis for each functional region show that RatHuV/YY12/CHN had a high amino acid identity with FeHuV in the complete polyprotein gene and different functional regions, especially in the P3 region (97%). Interestingly, the P1 structural protein showed fewer nucleotide and amino acid identities (mostly < 70%) between RatHuV/YY12/CHN and other hunnivirus reference sequences, demonstrating that genetic diversity is exhibited and the P1 region might be the most variable motif in the Hunnivirus genome (Table 2). Such highly variable motifs encoding the significant viral capsid protein might determine the pathogenicity and antigenicity for picornavirus [18]. Thereby, we used SLAC/FEL/REL methods in the HyPhy package when calculating the dN/dS values at each P1 sequence site, identifying strong negative selection on the non-synonymous sites. Taking SLAC results as an example, a total of 319 negatively selected sites and no positively selected site were reported with a statistical significance of 0.05. Additionally, the mean dN/dS ratio (0.0170) represented the average 98.3% of nonsynonymous mutations which were eliminated by negative selection during the evolution of virus.
A branch-specific analysis using the HyPhy package (GA Branch analysis) was subsequently generated because the dN/dS ratio between branches in the phylogenetic tree may vary when hunnivirus sequences in the alignment come from different host species ( Figure 5). The findings demonstrated evidence of complex and variable selective pressures on the hunnivirus sequences, with the strong support of negative selection along all branches representing the evolution in each host species. Therefore, the negative selection might be responsible for the evolution of P1 and even the entire coding regions of hunniviruses, which deserve putatively functional studies in the future.
This study predicted the positions of protease-cleavage sites in the rat hunnivirus viral proteins by sequence alignment with known cleavage sites in hunnivirus sequences, which contained E/G, A/D, Q/G, P/T, and G/P (Figure 2A). The cleavage sites between VP4 and VP2 (A/D), 2A and 2B (G/P), 3B and 3C (Q/G), and 3C and 3D (Q/G) are conserved among the feline, bovine, ovine, and several rat-derived hunniviruses, in line with a recent study reporting in detail on the protease-cleavage sites of hunniviruses [4]. We considered that understanding and feature analysis concerning the protein cleavage sites would be helpful for the in-depth investigations into the mechanism of protein cleavage [19]. Furthermore, recent studies have illustrated that functional non-coding RNAs have critical impacts on defects functionally related to various diseases [20,21]. Several studies have revealed the predicted RNA secondary structure of picornaviruses, while little is known currently on this structure in hunniviruses. In this study, the analysis of predicted RNA secondary structure showed that the 5 UTR region of our sequence can construct a stable secondary structure with the minimum free energy of −164.80kcal/mol (Figure 3). Since the functions of non-coding RNAs are highly associated with their structures, understanding the structures could clarify the functions of non-coding RNAs. Thus, more 5 UTR sequence data should be obtained to determine the genetic features of Hunnivirus.
The number of recently determined gastrointestinal viruses is expanding rapidly. To date, many picornaviruses are associated with diarrhea among humans and other animals. A previous review summarized the viruses that have been described in gastroenteritis cases by next-generation sequencing platforms, such as human Aichi virus [22,23], Cosavirus [24], and Salivirus [25][26][27]. The evidence combined with the epidemiological studies and statistical analysis showed that kobuviruses might be a potential causative agent for diarrhea disease in calves and pigs [28,29]. Although a feline hunnivirus described previously was found in a cat with diarrhea, the ability of hunniviruses to cause diarrhea disease remains to be determined. Therefore, identifying hunniviruses from different host species with and without diarrhea is needed to explore the potential association between the prevalence of Hunnivirus and diarrhea disease in the future.

Sample Collection
A total of 404 fresh fecal samples were collected between October 2015 and September 2017 from rats with unknown health status captured close to human residences using live traps in five different areas in southern China, including Xiamen in Fujian Province, Malipo in Yunnan Province, Yiyang in Hunan Province, and Guangzhou and Maoming in Guangdong Province ( Figure 6). Individual fresh stools were immediately placed in RNase-free tubes with 700 µL phosphate-buffered saline (PBS; 0.3% homogenate) and kept frozen at −80 • C until further use.

Nucleic Acid Extraction and cDNA Synthesis
The thawed stool specimens were fully resuspended in PBS and centrifuged at 8000× g for 10 min at 4 °C to collect the supernatants. Following the manufacturer's instructions, viral nucleic acid was extracted from 200 μL of each supernatant using the MiniBEST Viral RNA/DNA Extraction Kit (TaKaRa, Kusatsu, Japan). The obtained RNA was reverse transcribed to synthesize cDNA using a Transcriptor First Strand cDNA Synthesis Kit (Roche, Basel, Switzerland). The cDNA was directly used as the template for polymerase chain reaction (PCR) or stored at −20 °C.

Nucleic Acid Extraction and cDNA Synthesis
The thawed stool specimens were fully resuspended in PBS and centrifuged at 8000× g for 10 min at 4 • C to collect the supernatants. Following the manufacturer's instructions, viral nucleic acid was extracted from 200 µL of each supernatant using the MiniBEST Viral RNA/DNA Extraction Kit (TaKaRa, Kusatsu, Japan). The obtained RNA was reverse transcribed to synthesize cDNA using a Transcriptor First Strand cDNA Synthesis Kit (Roche, Basel, Switzerland). The cDNA was directly used as the template for polymerase chain reaction (PCR) or stored at −20 • C.

PCR Detection for Hunnivirus
During the present study, a hunnivirus sequence was serendipitously found in the fecal sample. Thus, we designed a pair of specific primers for rat hunnivirus to determine additional hunniviruses from the collected fecal samples (RHuV-F: 5 -TGGTGACCGGACTGAT GGACCC-3 , corresponding to nucleotides [nt] 6291-6312 of the sequence rat/NYC-E21/USA/2012 and RHuV-R: 5 -TCAGTTCAGCATGCAGCACCGG-3 , corresponding 2020). RNA secondary structure prediction was generated using RNAfold 2.4.13 (http: //rna.tbi.univie.ac.at/, accessed on 23 September 2020). Phylogenetic tree construction was generated by the maximum-likelihood method with 1000 bootstrap replications in MEGA v6.0 and visualized in FigTree v1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/, accessed on 8 September 2020). The model selections were based on the results of "Find best DNA/Protein models" in MEGA v6.0. A similarity plot analysis of the almost full-length genome was conducted using SimPlot 3.5.1 software.

Selective Pressure Analyses
The signatures of selection from homologous sequences were inferred by estimating the relative rates of non-synonymous (dN) and synonymous (dS) substitutions. To detect site-specific selective pressure on the codon alignment of P1 regions, we implemented the HyPhy package on the Datamonkey Webserver (www.datamonkey.org, accessed on 25 May 2021) [5] using three codon-based methods. The methods included single-likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), and random effects likelihood (REL). Further, the genetic algorithm (GA) Branch analysis was used to vary selection pressure between different phylogenetic branches using the non-synonymous and synonymous substitutions ratio (dN/dS). It noted that the sequence alignments came from translated protein sequences and then mapped aligned residues to codons to avoid introducing frameshifts and preserve codons.

Data Summary
The partial 3CD nucleotide sequences and near-complete genome of rat hunnivirus sequence have been lodged within the GenBank database under accession numbers MW417230 to MW417242.

Conclusions
In conclusion, the present study is expected to be the first to provide an epidemiological manifestation for rat hunnivirus from Rattus norvegicus and Rattus tanezumi in southern China, indicating that hunniviruses are common in rats close to humans' residences. Our findings from genetic analyses suggest that genetic diversity is displayed among the hunniviruses from different host species. Furthermore, the almost complete genome of RatHuV/YY12/CHN was successfully sequenced in this study. The phylogenetic trees based on the P1 coding regions demonstrate that the RatHuV/YY12/CHN belongs to the genotype of Hunnivirus A4. Strong negative selection was detected based on non-synonymous to synonymous substitution rates, providing novel insights into the viral evolution. Even though the list of bioinformatical analyses performed in the current study is supposed to be incomplete, the present results will be helpful to understand the limitation of the epidemic, molecular characteristics, and evolution of hunnivirus in rats in China. These provide an excellent beginning to dive even deeper into viral analysis. Therefore, there is a high need for research focused on the epidemiology, geographic distribution, heterogeneity, pathogenesis, transmission route, and the potential risk of cross-species transmission of rat hunnivirus.

Institutional Review Board Statement:
The research protocol has been approved by the Animal Ethics and Welfare Committee of the School of Public Health, Southern Medical University and met the guidelines for the Rules for the Implementation of Laboratory Animal Medicine (1998) from the Ministry of Health, China.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated or analyzed during this study are included in this published article. Access to raw data can be acquired by contacting the corresponding author via email.