Diversity Analysis of Tick-Borne Viruses from Hedgehogs and Hares in Qingdao, China

ABSTRACT Tick-borne viruses (TBVs) have attracted increasingly global public health attention. In this study, the viral compositions of five tick species, Haemaphysalis flava, Rhipicephalus sanguineus, Dermacentor sinicus, Haemaphysalis longicornis, and Haemaphysalis campanulata, from hedgehogs and hares in Qingdao, China, were profiled via metagenomic sequencing. Thirty-six strains of 10 RNA viruses belonging to 4 viral families, including 3 viruses of Iflaviridae, 4 viruses of Phenuiviridae, 2 viruses of Nairoviridae, and 1 virus of Chuviridae, were identified in five tick species. Three novel viruses of two families, namely, Qingdao tick iflavirus (QDTIFV) of the family of Iflaviridae and Qingdao tick phlebovirus (QDTPV) and Qingdao tick uukuvirus (QDTUV) of the family of Phenuiviridae, were found in this study. This study shows that ticks from hares and hedgehogs in Qingdao harbored diverse viruses, including some that can cause emerging infectious diseases, such as Dabie bandavirus. Phylogenetic analysis revealed that these tick-borne viruses were genetically related to viral strains isolated previously in Japan. These findings shed new light on the cross-sea transmission of tick-borne viruses between China and Japan. IMPORTANCE Thirty-six strains of 10 RNA viruses belonging to 4 viral families, including 3 viruses of Iflaviridae, 4 viruses of Phenuiviridae, 2 viruses of Nairoviridae, and 1 virus of Chuviridae, were identified from five tick species in Qingdao, China. A diversity of tick-borne viruses from hares and hedgehogs in Qingdao was found in this study. Phylogenetic analysis showed that most of these TBVs were genetically related to Japanese strains. These findings indicate the possibility of the cross-sea transmission of TBVs between China and Japan.

Phylogenetic analyses showed that two novel viruses, QDTPV and QDTUV, belong to the family Phenuiviridae (Fig. 4F). In this study, Dabie bandavirus strain QDH05 was clustered with strains isolated previously from Shandong, China (Fig. 4B); it had the closest relationship with a strain of Dabie bandavirus (GenBank accession number KR706567.1) that was isolated from a patient's serum (99.8% nucleotide identity) in China. From the phylogenetic tree, QDTPV was evolutionarily located on a previous branch of MUKV strain HLJ (GenBank accession number YP_009666332.1) and Kuriyama virus (GenBank accession number UXL90891.1) (13, 24) between tick-borne viruses and sandfly/mosquito-borne viruses (Fig. 4C). QDTUV, which had the closest genetic relationship with Toyo virus (13), was assigned to a subgroup of the Kaisodi group ( Fig. 4D) with Silverwater virus (GenBank accession number YP_010086157.1), Huangpi tick virus 2 (GenBank accession number YP_009293590.1) (17), and Kaisodi virus (GenBank accession number AWW17495.1) on the phylogenetic tree, and pairwise similarity analysis based on the full sequences showed up to 90% nucleotide acid identity among four QDTUV strains (H01 and H04 to H06) (Fig. 4F). The M-segment-deficient phleboviruses   (MdPVs), which have only L and S segments, were grouped into three clusters (Group I to III) on the phylogenetic tree (15). OKTV (QDH01 to QDH06, QDR01, QDR02, and QDR04) was assigned to cluster GIII, and pairwise similarity analysis based on the full sequence showed up to 96% nucleotide identity among the currently sequenced OKTVs (Fig. 4F).
(c) CPTV3 of Chuviridae. In this study, five strains of Changping tick virus 3 (QDR01 to QDR04 and QDH03) of Chuviridae were identified in R. sanguineus from hedgehogs and D. sinicus from hares. The virus genome encodes the RdRp, the glycoprotein (G), and the N protein (Fig. 6A). A phylogenetic tree based on the viral RdRp sequences of Chuviridae was constructed. CPTV3 strains QDR01 to QDR04 and QDH03 form a subgroup with CPTV3 (GenBank accession number YP_009177707.1). CPTV3 strains QDR01 to QDR04 and QDH03 were clustered with many other viruses isolated from ticks, including Tacheng tick virus 5 (17) and Wuhan tick virus 2 (17) (Fig. 6B). In addition, the  Table S4 in the supplemental material. Diversity of the Tick Virome Microbiology Spectrum sequence similarity between these five strains and the CPTV3 reference strain was up to 97% (Fig. 6B).

DISCUSSION
This study found 36 strains of 10 viruses from 4 viral families in H. flava, R. sanguineus, D. sinicus, H. longicornis, and H. campanulata from hares and hedgehogs in Qingdao, China, by metatranscriptomics. Among them, three novel viruses, QDTIFV, QDTPV, and QDTUV, were observed. Of the seven known viruses, both HfIFV and OKTV were first detected in China. In addition, we also found some short fragments of suspected new viruses, between 2,000 and 3,000 bp in length. They have the highest identities with Xinjiang tick-associated virus 1 and Ixodes scapularis-associated virus 2,   Table S4 in the supplemental material. respectively (see Table S3 in the supplemental material). Our results showed that there were large numbers of viruses of the Phenuiviridae, Iflaviridae, Nairoviridae, and Chuviridae in the five species of ticks in Qingdao. Our research revealed the diversity and abundance of viruses carried by different tick species (26,27), which enriches the viral diversity of ticks in China.
Phenuiviridae can infect animals, plants, and fungi, which is rare among known virus families (28,29). Phenuiviridae are arboviruses that can replicate in different hosts such as insects, humans, and rice (22,30). Many kinds of Phenuiviridae, such as Dabie bandavirus (3), Rift Valley fever virus (RVFV) (31), Heartland virus (HRTV) (32), and rice stripe virus (33), are highly pathogenic to humans, animals, or plants, and they impose a heavy global burden on human health, animal husbandry and agriculture. SFTSV was first reported in 2011 as a representative virus of the Phenuiviridae (3). In recent years, more than 1,000 cases of SFTSV infection have been diagnosed annually in China (34 Table S4 in the supplemental material.

Diversity of the Tick Virome Microbiology Spectrum
Additionally, South Korea, Japan, and Vietnam have also found human cases of SFTSV infection (35)(36)(37)(38). Currently, many cases of SFTSV infection and studies on ticks carrying SFTSV have been reported in Shandong, China (39), which suggests that SFTSV is widespread in ticks in Shandong, China. Two newly discovered viruses, QDTPV and QDTUV, also belong to the Phenuiviridae. Currently, it is unclear whether the two novel viruses can infect humans or animals to cause new infectious diseases. Many of the known viruses found in this study were first reported in Japan. These findings indicate that the viruses in China and Japan may originate from the same ancestor. Tick-borne viruses are transmitted along with migratory birds across the East China Sea and/or the Sea of Japan. Globally, migratory birds are transporters of tickborne pathogens such as Borrelia, Rickettsia, and Crimean-Congo hemorrhagic fever virus (CCHFV) (40)(41)(42)(43), and tick-borne viruses can be transmitted between humans and livestock through travel and trade. In the future, the detection of these viruses in ticks carried by migratory birds would be powerful evidence. To understand the evolutionary relationships between tick-borne viruses in China and Japan, we performed a phylogenetic analysis on 52 TBVs from China and/or Japan in recent years. The results showed that these viruses were classified into at least nine virus families, including Phenuiviridae, Chuviridae, Rhabdoviridae, Nyamiviridae, Flaviviridae, Nairoviridae, Spinareoviridae, Iflaviridae, and Orthomyxoviridae (Fig. 7). The multiple viruses found in the two countries are distributed on branches of the same virus genus on the phylogenetic tree, indicating a close evolutionary relationship. In our study, HfIFV and OKTV strains were found in China, with up to 99% amino acid identity with strains previously found in Japan. QDTPV and QDTUV were also found in China and showed 80.9% and 90.5% amino acid identities of RdRp with Mukawa virus and Toyo virus previously reported in Japan, respectively. The RdRp protein has the highest level of homology compared to other proteins, indicating that RdRp, as an important enzyme in the process of viral replication that participates in the RNA genome replication process, is very evolutionarily conserved. The diversity of the M segment and S protein enhances the adaptability and immune escape of the virus in the host.
Since only the types of viruses carried by ticks on hares and hedgehogs have been identified, the viruses found in this study are not all species of viruses carried by ticks in Qingdao.
In conclusion, 36 strains of 10 viruses from 4 viral families were found using metatranscriptomics sequencing in this study. Three novel TBVs, QDTIFV, QDTPV, and QDTUV, were found in Qingdao, China. Half of the TBVs in this study are genetically closely related to viruses of the same genus from Japan. These findings provide indirect evidence for the cross-sea transmission of TBVs between China and Japan. Therefore, it is clear that there is a great diversity of viruses in blood-sucking insects of wild animals in China.

MATERIALS AND METHODS
Sample collection. A total of 662 ticks were collected from June to July 2019 in the suburbs of Qingdao, Shandong Province, China, of which 420 were collected from hedgehogs and 242 were collected from hares. Ticks were identified by morphological identification and mitochondrial 16S rRNA sequencing (44). All samples were stored at 280°C throughout the whole process.
RNA library construction and sequencing. The collected ticks were divided into 10 groups according to tick species, and the number of ticks in each group was between 60 and 80. Next, the ticks in each group were washed with phosphate-buffered saline (PBS) three times. One milliliter of Dulbecco's modified Eagle's medium (DMEM) containing 10% fetal bovine serum (FBS) was added to the ticks in each group, and the ticks were ground at a low temperature to homogenize them completely. Pure RNA was obtained by phenol-chloroform extraction and isopropyl alcohol precipitation. The library of RNA was constructed by using the New England BioLabs (NEB) total RNA library building kit according to the manufacturer's instructions. Paired-end (100 bp) sequencing of each RNA library was performed on the Illumina NovaSeq6000 instrument. The sequencing of the qualified library was performed using an Illumina NovaSeq6000 instrument in PE 2Â 150-bp mode.
Metagenomic bioinformatics analyses. Quality control of the next-generation sequencing raw data was performed by using fastp v0.20.0 (45), including adapter trimming and quality filtering. Ribosome fragment elimination was estimated by classifying the reads with SortMeRNA v2.0 (version 4.3.2) (46) against the SILVAdb128 small-subunit (SSU) and large-subunit (LSU) databases (47 The reads were then de novo assembled with metaSPAdes v3.13.1. The assembled contigs were classified into known viral orders and families mapping to the NCBI nt database by BLASTn (parameter of an E value of ,1e25). Predicted viral contigs were mapped against sequences in the NCBI nr database by using diamond software (parameter of an E value of ,1e25) (49). The reads were classified into known viral species using kraken2 (50). At the same time, the species of reads were annotated using the NCBI virus database through the diamond BLASTx program (parameter of an E value of ,1e25).
Validation and annotation of virus genomes. To rule out false-positive results during assembly, specific PCR primers and probes designed according to viral RdRp or polyprotein sequences were used to verify the virus found by NGS in the original sample using reverse transcription-PCR (RT-PCR) (see Table S2 in the supplemental material for details of detection). The potential open reading frames (ORFs) of virus sequences were predicted using ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/). The conserved regions in the sequence were annotated with a CDD search (https://www.ncbi.nlm.nih .gov/Structure/cdd/wrpsb.cgi).
Virus classification. The discovered viruses were classified based on nucleotide and amino acid identities (Table S3). If the nucleotide identity of a new virus species with the whole genome of the reference Diversity of the Tick Virome Microbiology Spectrum virus strain is less than 80% or the amino acid identity with the RNA-dependent RNA polymerase (RdRp) domain of a known virus is less than 90%, the virus is defined as a new virus (27). All new viruses were named "Qingdao tick" (QDT), followed by the names of common viruses classified by them in this study. To distinguish them from other known virus strains, the strains detected in this study are named "QD." Phylogenetic analysis. The nucleotide and amino acid sequences of the viruses from this study were aligned with published reference sequences for the same genus. The nucleotide and amino acid sequences of the alignments were analyzed and compared with representative sequences from GenBank using MAFFT v7.450 (51). Phylogenetic analysis was performed using the full-length nucleotide sequence or the RdRp sequence of the representative members of related viral species or genera. Phylogenetic trees were constructed by the maximum likelihood method based on 1,000 repeated bootstrap replicates using MEGA7 (52). The GenBank accession numbers of the viruses used in this article are detailed in Tables S4 and S5.
Data availability. The nucleotide sequence accession numbers of the viruses detected in this study have been submitted to the international nucleotide sequence database (GenBank) and the National Microbiology Data Center (NMDC) ( Table S3). All sequence reads generated in this project are available in the NCBI Sequence Read Archive under BioProject accession number PRJNA935925. All viral genome sequences have been submitted to the GenBank database under accession numbers OQ513629 to OQ513687 (Table S3).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 7.6 MB.