Extensive diversity of RNA viruses in ticks revealed by metagenomics in northeastern China

Background Ticks act as important vectors of infectious agents, and several emerging tick-borne viruses have recently been identified to be associated with human diseases in northeastern China. However, little is known about the tick virome in northeastern China. Methods Ticks collected from April 2020 to July 2021 were pooled for metagenomic analysis to investigate the virome diversity in northeastern China. Results In total, 22 RNA viruses were identified, including four each in the Nairoviridae and Phenuiviridae families, three each in the Flaviviridae, Rhabdoviridae, and Solemoviridae families, two in the Chuviridae family, and one each in the Partitiviridae, Tombusviridae families and an unclassified virus. Of these, eight viruses were of novel species, belonging to the Nairoviridae (Ji’an nairovirus and Yichun nairovirus), Phenuiviridae (Mudanjiang phlebovirus), Rhabdoviridae (Tahe rhabdovirus 1–3), Chuviridae (Yichun mivirus), and Tombusviridae (Yichun tombus-like virus) families, and five members were established human pathogens, including Alongshan virus, tick-borne encephalitis virus, Songling virus, Beiji nairovirus, and Nuomin virus. I. persulcatus ticks had significant higher number of viral species than H. japonica, H. concinna, and D. silvarum ticks. Significant differences in tick viromes were observed among Daxing’an, Xiaoxing’an and Changbai mountains. Conclusions These findings showed an extensive diversity of RNA viruses in ticks in northeastern China, revealing potential public health threats from the emerging tick-borne viruses. Further studies are needed to explain the natural circulation and pathogenicity of these viruses.


Introduction
Ticks are obligate haematophagous ectoparasites of wild and domestic animals as well as humans, and there are approximately 900 tick species worldwide, of which many can transmit pathogenic agents, including viruses, bacteria, and protozoa [1]. Several tick-borne viruses (TBVs) are associated with serious diseases in humans and animals, such as tick-borne encephalitis virus (TBEV), Crimean-Congo hemorrhagic fever virus (CCHFV) and Nairobi sheep disease virus (NSDV) [2]. In recent decades, the incidence and geographical distribution of tick-borne viruses have an increasing tendency, highlighting the public health importance of these arboviruses [3]. Due to the application of next generation sequencing (NGS) in recent years, many novel viruses have been identified in different tick species of different regions worldwide [4][5][6][7]. To our knowledge, the known identified TBVs include hundreds of viral members of at least 12 genera in 9 families of two orders as well as other unassigned members [8].
The northeastern region has the richest forest resources in China, mainly concentrated in the Daxing'an mountain (DXAM), Xiaoxing'an mountain (XXAM) and Changbai mountain (CBM), resulting in abundant tick populations. In addition to TBEV, several emerging tickborne viruses that may infect humans have been discovered in northeastern China, such as ALSV [16], SGLV [17], BJNV [18], and JMTV [15]. Here, using metagenomic analysis, we found extensive diversity of RNA viruses in ticks in northeastern China, revealing potential public threats from the emerging TBVs. Further studies are needed to explain the natural circulation and pathogenicity of these viruses.

Sample collection
From April 2020 to July 2021, the questing ticks were collected by the flagging method in Northeastern China, and the blood-sucking ticks were only collected from cattle in Shulan (Fig 1). These ticks were identified to species as described elsewhere [29,30]. The questing ticks were pooled to 10 ticks per tube based on the tick species and sampling sites, and Dermacentor silvarum ticks collected from cattle were grouped and pooled to 1, 5, or 10 ticks per tube based on the engorgement levels (fully engorged, partially engorged, or unengorged). All the ticks were stored at -80˚C until use.

RNA library construction and sequencing
After washing with 75% ethanol and RNA/DNA-free water, pooled ticks in tubes were added with 500 μL Dulbecco's modified Eagle's minimum essential medium (DMEM) and two stainless steel beads (3 mm diameter), and crushed using the Tissuelyser (Jingxin, Shanghai, China) at 70 Hz for 2 min. The lysates were centrifuged at 12000 rpm for 10 min at 4˚C, and the supernatant was further pooled for library construction according to the collection sites and species (S1 Table). After digested with micrococcal nuclease (NEB, USA) in 37˚C for 2 h, the pooled samples were used for viral RNA extraction with the TIANamp Virus RNA kit (TIANGEN, Beijing, China).
The extracted RNA was subjected to metagenomic sequencing at Tanpu Biological Technology (Shanghai, China). Briefly, the RNA from each pool was used for library preparation with the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer's instructions. After adapter ligation, 10 cycles of PCR amplification were performed for target enrichment. The libraries were pooled at equal molar ratio, denatured and diluted to optimal concentration, and sequenced with an Illumina NovaSeq 6000 System. [32,33]. After being compared with the nonredundant nucleotide (nt) and protein (nr) database downloaded from GenBank using BLAST+ v2.10.0, the assembled contigs were filtered to remove the host and bacterial sequences. The relative abundance of the identified viruses was determined by mapping the reads back to the assembled contigs using Bowtie2 v2.3.3.1.

Viral genome confirmation and annotation
The assembled contigs were compared with NCBI nucleotide and viral refseq database using BLAST (V2.10.0+), and used as a reference for designing specific primers to confirm and analyze the sequences of terminal ends, using the nested reverse transcription-polymerase chain reaction (RT-PCR) and the rapid amplification of cDNA ends (RACE) as described elsewhere [17,18]. Potential open reading frames (ORFs) in the viral sequences were predicted using ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/).

Virus classification
All the viruses identified in this study were classified according to the latest International Committee on Taxonomy of Viruses (ICTV) report of virus taxonomy (https://talk.ictvonline. org/ictvreports/ictv_online_report/). A novel viral species should be satisfied with one of the  following conditions as described before [24], namely, (i) <80% nucleotide (nt) identity across the complete genome; or (ii) <90% amino acid (aa) identity of the RNA-dependent RNA polymerase (RdRp) domain with the known viruses. All novel viruses were named as the collection sites that the virus was first identified, followed by common viral names according to their taxonomy. All the viral strains would be marked with 'Northeastern (NE)' to distinguish them from the virus strains identified in other studies.

Phylogenetic analyses
To confirm the phylogenetic relationships of the viruses discovered in this study, representative reference viral sequences were retrieved from the GenBank database (S2 Table), and aligned using ClustalW available within MEGA 7.0. Phylogenetic analyses were conducted with the aligned sequences using the maximum-likelihood method in MEGA version 7.0 with the best-fit substitution model for each alignment [34]. A bootstrapping analysis of 1000 replicates was conducted in the analysis, and the bootstrap values more than 70 were shown in the trees.  Table). The collection sites were distributed in Ji'an (n = 100), Dunhua (n = 382), and Shulan (n = 300) in Jilin Province, Fangzheng (n = 300), Mudanjiang (n = 156), Yichun (n = 253), and Tahe (n = 316) in Heilongjiang Province, and Songling (n = 224) in Inner Mongolia Autonomous Region (Fig 1, S1 Table). Of these, five were located in Changbai Mountain (CBM), two in Daxing'an Mountain (DXAM), and one in Xiaoxing'an Mountain (XXAM) (Fig 1, S1 Table).

Identified RNA viruses
A total of 24 RNA libraries were constructed and sequenced, resulting in 245.8 GB clean data and *0.5 billion non-rRNA reads (S1 Table). Totally, 2,059 viral contigs were obtained by de novo assembly from *0.93 million viral reads that accounted for 0.2% of the total non-rRNA reads. Within each library, the viral reads ranged from 0.02% (library SL3) to 0.53% (library JA1) of the total non-rRNA reads. After being aligned by Blast and filtered by virus-host database, the viral contigs were finally annotated to 22 viruses, belonging to the 8 viral families, including Flaviviridae, Rhabdoviridae, Nairoviridae, Phenuiviridae, Chuviridae, Partitiviridae, Tombusviridae, and Solemoviridae, and one unclassified viral species (Table 1).

Virome composition and abundance
Across the 24 libraries, 22 possessed 1-11 viral species, with the exception of the libraries FZ1 and YC1, which had no virus detectable (Fig 2A).  2B).

Virus diversity and evolution
Flaviviridae. According to the latest ICTV report of virus taxonomy, Flavivirus, Hepacivirus, Pegivirus, and Pestivirus are the approved genera in family Flaviviridae. While in recent years, a flavivirus-like group, the Jingmenvirus group, including a series of genetically related viruses, such as Jingmen tick virus [14] and ALSV [16], have been discovered (Fig 3A). In this study, ALSV was identified in the I. persulcatus library (TH4) from Tahe in DXAM, and clustered together with ALSV strain H3 isolated from tick-bitten patients in NE China (Fig 3B) [16], but different from the isolates detected in I. persulcatus in Russia and I. ricinus in Finland (Fig 3B), with nt identities of 97.0-98.5% (S4 and S5 Tables) [35,36].
Bole tick virus 4 (BLTV4)-NE strains were phylogenetically grouped into the pestivirus-like group, with nt identities of 77.5-83% and RdRp aa identities of 94.6-96.9% to other BLTV4 . All the viruses obtained in ticks here were highlighted in red. In panel A, closest referenced viruses were highlighted in bold font, and the virus genera or groups were marked with different colors of background. In panel B, ALSV isolated from humans were marked with blue. ALSV were highlighted with light green background. The host tick species of the viruses and the countries that the viruses discovered were also labeled. In panel C, TBEV isolated from CBM and XXAM were marked with blue-filled circles, while strains found in DXAM were marked with red-filled circles. In panel D, the strain names of BLTV4 were labeled. The host tick species and the countries that BLTV4 discovered were also marked. ALSV, Alongshan virus; TBEV, Tick-borne encephalitis virus; BLTV4, Bole tick virus 4. The accession numbers of the viral sequences are shown in S2 and S3 Tables.
Nairoviridae. In the phylogenetic tree of the family Nairoviridae, SGLV-NE, together with SGLV, formed a separate clade from other viral members in the genus Orthonariovirus, including Ji'an nariovirus, Tacheng tick virus 1, and Tamdy virus (Fig 4A). SGLV-NE was identified in two libraries TH1 and TH2 of the H. concinna ticks from Tahe in DXAM, and clustered together with SGLV strains isolated from tick-bitten patients (Fig 4B), with nt identities of 92.4-99.2% (S8 and S9 Tables) [17].
JANV, genetically related to SGLV with nt identities of 70.7-73.5%, was a novel identified nairovirus (S8 and S9 Tables). Four libraries were detected JANV-positive, including DH1 and JA libraries of H. japonica ticks in CBM, and MDJ1 and YC2 of H. concinna ticks in XXAM (Fig 4B).
BJNV and YCNV, belonging to an unclassified Norwavirus-like group, showed close relationships to Norway nairovirus 1 and Grotenhout virus (Fig 4A), with nt identities of 75.4-79.6% (Segment L and S, S10 and S11 Tables). BJNV was identified in seven I. persulcatus tick libraries in all the three regions, suggesting the wide distribution of the virus in NE China. All the BJNV NE strains were clustered together with BJNV strains identified in tick-bitten patients and I. persulcatus tick in NE China and GSTV detected in Russia (Fig 4C), with nt identities of 96.2-100% (Segment L and S, S10 and S11 Tables).
YCNV was only detected in three libraries (FZ3, YC3, and YC4) of I. persulcatus ticks from Fangzheng and Yichun in CBM and XXAM. Although YCNV showed nt identities of 82.3- 83.9 (>80) (Segment L and S, S10 and S11 Tables) with BJNV, the virus formed a separate clade, with aa (RdRp) identities of 87.7-88.5% (Fig 4C, S10 and S11 Tables), indicating that YCNV was a novel viral species that may be different from BJNV.
Phenuiviridae. MKWV and MJPV, together with STPV and OTPV identified here, fell within the genera Phlebovirus and Ixovirus of the family Phenuiviridae, respectively (Fig 5A). MKWV NE strains were clustered together with the strain MKW73 identified in I. persulcatus ticks in Japan, with nt identities of 92.4-100% for L segment, 90.7-100% for M segment, and 93.5-99.9% for S segment (Fig 5B, S12 and S13 Tables). Of the six MKWV NE strains, five were detected in the I. persulcatus ticks from Shuangzi, Dunhua, Yichun, and Tahe, while one was identified in the D. silvarum ticks collected from cattle in Shulan.

PLOS NEGLECTED TROPICAL DISEASES
Extensive diversity of RNA viruses in ticks in northeastern China MJPV was only identified in two libraries in the I. persulcatus ticks from Fangzheng and Mudanjiang in CBM, which was distantly related to Kuriyama virus and Mukawa virus strains, with low nt identities of 75.8-77.4% for L segment, 67.6-68.7% for M segment, and 69.8-70.8% for S segment (Fig 5B, S12 and S13 Tables) [42].
STPV and OTPV NE strains found in this study were clustered with STPV and OTPV strains discovered in I. persulcatus ticks from Karelia in Russia, with nt identities of 98.0-99.2% and 98.6-99.1%, respectively (Fig 5C, S14 and S15 Tables). The two viruses formed separate clades from each other, with nt identities of 56.4-57.3%, and were identified in paired in I. persulcatus ticks in six libraries from three collection sites, including Yichun in XXAM, and Tahe and Songling in DXAM (Fig 5C, S14 and S15 Tables).
Chuviridae. Nuomin virus (NUMV) and Yichun mivirus (YCMV) identified here fell within the Mivirus and Nigecruvirus genera in the Chuviridae family, respectively (Fig 7A). NUMV NE strains were detected in all the ten I. persulcatus libraries in the three regions in NE China, while no libraries from other tick species were identified NUMV positive, indicating the wide distribution and specificity of host tick species of the virus. The viral strains of NUMV found in this study were clustered with other NUMV strains discovered in humans from NE China and Lesnoe mivirus isolated from the I. persulcatus ticks in Russia, showing nucleotide identities of 96.6-99.5% with each other (S17 Table).
YCMV was a novel virus identified in this study, and clustered together with Blacklegged tick chuvirus 2 found in I. scapularis ticks from USA, with nt identities of 68.5%. YCMV was only detected in the I. persulcatus ticks from Dunhua and Yichun in CBM and XXAM ( Fig  7C), and shared high nt identity of 99.2%.
Plant-related viruses. Several viral species closely related to plant viruses were found in this study. The Jilin partiti-like virus 1 (JPLV1) identified here fell within an unclassified clade, provisionally designated as Deltapartitivirus-like group related to the Partitiviridae family ( Fig  8A), and was closely related to Norway partiti-like virus 1 strains in Norway with high identities of nt 91.3-91.9% (S18 Table). The newly discovered Fangzheng tombus-like virus (FTLV) fell within the Tombusviridae-like group associated with family Tombusviridae (Fig 8B), sharing a close relationship (identity: nt 73.4-73.6% and aa 80.6-81.3%) to Upmeje virus strain OTU9.IU20 that has been identified in Ixodes uriae in Sweden (S19 Table) [43].
Three Solemoviridae-like viral species, including Ixodes scapularis associated virus 1 (ISAV1), Xinjiang tick associated virus 1 (XTAV1), and Jilin luteo-like virus 2 (JLLV2), were found in this study, and fell within the Sobemo-like virus group associated with the family Solemoviridae (Fig 8C). In the aa RdRp phylogenetic tree, ISAV1-NE was clustered with Norway luteo-like virus 3 identified in Ixodes Ricinus ticks in Norway and Ixodes scapularis associated virus 1 identified in I. scapularis ticks in USA, with nt identities of 81.5-85.7% (S20 Table). The XTAV1-NE strains were clustered together with XTAV1 strains in Xinjiang with nt identities of 95.9-96.2% (S20 Table). The JLLV2 strains were detected in the I. persulcatus ticks and showed a close relationship with Norway luteo-like virus 2 identified in Ixodes ricinus ticks in Norway, with nt identities of 87.7-88.2% (S21 Table).
Ixodes scapularis associated virus 3 (ISAV3) is still unclassified to date, but the strains identified in this study showed close relationship with ISAV3 and ISAV4 discovered in Ixodes scapularis ticks from USA with nt identities of 86.7-88.2% (Fig 8D, S22 Table) [44,45].
There were significantly differences in tick viromes among CBM, XXAM, and DXAM. In this study, ALSV, TBEV, and THRV2-3 were only detected in DXAL, while MJPV, YCNV, YCMV, and FTLV were identified in CBM and XXAL. Nairovirus SGLV was detected in DXAM, while its closely related virus JANV was detected in CBM and XXAM; Rhabdovirus THRV1 tended to evolve into two genotypes: one was distributed in DXAM, the other distributed in CBM and XXAM. Additionally, the genetic differences of TBEV also supports the occurrence of geographical barriers of tick-borne viruses in northeastern China.
Segmented flaviviruses have recently been reported as emerging tick-borne viruses, with a wide distribution in Asia, Africa, Europe, Central America, and South America [46,47]. Of them, Jingmen tick virus (JMTV) and ALSV are associated with the febrile illness in tick-bitten patients in northeastern China. JMTV has been found in various vertebrates, including cattle, sheep, rodents, and non-human primate, showing that the virus cocirculates between ticks and mammals [46]; JMTV and ALSV have also been detected in mosquitoes [14,16]; however, the transmission modes of these viruses remain to be investigated. In this study, only I. persulcatus tick was tested ALSV-positive in DXAM, where ALSV patients have been recently found [16], suggesting potential public health risk of ALSV infection and limited distribution of the emerging tick-borne virus in northeastern China. The pathogenicity of segmented flaviviruses in humans and animals needs to be further verified through animal infection models.
Interestingly, co-feeding transmission might affect the viromes in ticks collected from animals. In this study, two viral species, THRV1 and MKWV, identified in D. silvarum collected from cattle (Shulan), were mainly found in H. japonica and H. concinna, and I. persulcatus ticks, respectively. However, the two viruses were not detected in the questing D. silvarum ticks in Dunhua adjacent to Shulan. As cattle can act as hosts of ticks, such as D. silvarum, H. japonica, H. concinna, and I. persulcatus ticks, and co-feeding transmission may be an efficiency way of viral transmission in ticks [52][53][54]. D. silvarum collected from cattle here may contract the viruses from Haemaphysalis sp. or I. persulcatus ticks by co-feeding transmission, but not the actual vector of these viruses. Thus, it is suggested to investigate tick virome PLOS NEGLECTED TROPICAL DISEASES diversity using questing ticks instead of ticks collected from animals. Moreover, further studies should be focus on the vector competence of ticks for the transmission of the identified viruses, which may further confirm the roles of different tick species for virus transmission.
The Bunyavirales order includes important human pathogens, such as Crimean-Congo hemorrhagic fever virus in the Nairoviridae family, and SFTSV and Valley fever virus in the Phenuiviridae family, whose genome are negative single-stranded RNA of small (S), medium (M), and large (L) segments, encoding structural nucleoprotein (NP), glycoprotein precursor (GPC), and RNA-dependent RNA polymerase (L) proteins, respectively. Recently, several bisegmented viruses without the M gene have been found in nairoviruses (such as Gakugsa tick virus and Norway nairovirus 1) and phleboviruses (such as Tacheng tick virus 2) [48][49][50][51]. We also identified four viruses, including YCNV and BJNV within the Nairoviridae family, and STPV together with OTPV in the Phenuiviridae family. However, we did not find the M gene encoding the glycoproteins. The possible reason may be the substantially genetical diversity of these viruses from the reference viruses [18,19].
We also found several viral species closely related to plant viruses, including JPLV1 in the Partitiviridae family, FTLV in the Tombusviridae family, and ISAV1, XTAV1, and JLLV2 in the Solemoviridae family. Compared with viruses in the Flaviviridae, Nairoviridae and Phenuiviridae families, plant-related viruses identified here may have low pathogenicity to humans or mammals.
There are some limitations to the present study. Although some tick-borne viruses associated with diseases in humans or mammals, including TBEV, ALSV, SGLV, BJNV, and NUMV, were detected in this study, some other pathogenic viruses, such as severe fever with thrombocytopenia syndrome virus (SFTSV) [55,56], lymphocytic choriomeningitis virus (LCMV) [57], Nairobi sheep disease virus (NSDV) [21], and Jingmen tick virus (JMTV) [15] identified in previous studies, were not detected here, indicating that larger tick samples and wider sampling sites are necessary to characterize the tick-borne viruses in NE China. Some novel viruses had close relationship with TBVs of public health significance. For example, JANV and YCNV were genetically related to SGLV and BJNV, respectively, which have been shown to be associated with febrile diseases, indicating the potentially pathogenic to humans and animals of these novel viruses. As in the Phenuiviridae family, however, although MKWV does not detected in tick-bitten patients, it can grow in a human-derived cells and mice, and its nonstructural proteins can suppress the anti-innate immune responses [42], suggesting the potential public health significance of MKWV and its closely related virus (MJPV), and the necessity of epidemiology studies on tick-bitten patients, livestock, and even wild animals.

Conclusions
These findings showed an extensive diversity of RNA viruses in ticks in northeastern China, revealing potential public health threats from the emerging tick-borne viruses. Further studies are needed to explain the natural circulation and pathogenicity of these viruses.
Supporting information S1