Tick virome diversity in Hubei Province, China, and the influence of host ecology.

Ticks are important vector hosts of pathogens which cause human and animal diseases worldwide. Diverse viruses have been discovered in ticks; however, little is known about the ecological factors that affect the tick virome composition and evolution. Herein, we employed RNA sequencing to study the virome diversity of the Haemaphysalis longicornis and Rhipicephalus microplus ticks sampled in Hubei Province in China. Twelve RNA viruses with complete genomes were identified, which belonged to six viral families: Flaviviridae, Matonaviridae, Peribunyaviridae, Nairoviridae, Phenuiviridae, and Rhabdoviridae. These viruses showed great diversity in their genome organization and evolution, four of which were proposed to be novel species. The virome diversity and abundance of R. microplus ticks fed on cattle were evidently high. Further ecological analyses suggested that host species and feeding status may be key factors affecting the tick virome structure. This study described a number of novel viral species and variants from ticks and, more importantly, provided insights into the ecological factors shaping the virome structures of ticks, although it clearly warrants further investigation.

Introduction

As one of the major hematophagous ectoparasites of vertebrates, ticks (Acari: Ixodoidea) pose a serious threat to public health worldwide. Direct effects include damage to skin, blood loss, injection of toxins, or even mortality of the hosts, where severe damages are mainly caused by the wide spectrum of pathogens that they harbor, including bacteria, fungi, protozoa, and viruses (Jongejan and Uilenberg 2004; Fang et al. 2015; Getahun et al. 2016). Through feeding different animals, particularly mammals, birds, and terrestrial reptiles, ticks act as efficient vectors of disease transmission, and diverse zoonotic viruses have been identified in ticks. In China, emerging tick-borne viruses have been reported to cause a number of human diseases: severe fever with thrombocytopenia virus (SFTSV) (Yu et al. 2011), Jingmen tick virus (JMTV) (Emmerich et al. 2018; Jia et al. 2019), Alongshan virus (ALSV) (Wang et al. 2019), and the recently identified Songling virus (SGLV) (Ma et al. 2021) and Beiji nairovirus (Wang et al. 2021), highlighting the necessity for routine surveillance of tick-borne viruses.

The advent of high-throughput sequencing (HTS) has greatly speeded up tick viromic studies and has updated our understanding of the virosphere. These viruses showed extremely high variety in genome organization and genetic diversity, with some likely to be the phylogenetic ancestors of vertebrate-infecting members (Li et al. 2015; Shi et al. 2016a,b; Tokarz et al. 2018; Pettersson et al. 2020; Stanojević et al. 2020; Wille et al. 2020b). Exploring the associated factors that affect the tick virome composition and evolution from an ecosystem perspective would recur the ecological events and provide a baseline to assess the probable threats of these viruses to the public (French and Holmes 2020). However, only limited data on virus-tick ecology are available to date. Knap and Avšič-Županc reviewed the studies on tick-borne encephalitis in Slovenia, proposing that forest and agricultural areas were suitable habitats for ticks and were important for tick-borne encephalitis virus (TBEV) establishment and transmission (Knap and Avšič-županc 2015). Benefitting from HTS, Jia and colleagues demonstrated that ecological and geographic factors might have shaped the genetic structure and pathogen composition in several ixodid tick species in China (Jia et al. 2020).

Haemaphysalis longicornis and Rhipicephalus microplus are the two dominant tick species in China. Although they differ from each other in morphological features, life cycle, and host spectrum, both are well known for their vector capacities of diverse pathogens. Haemaphysalis longicornis has a three-host life cycle and parasitizes mammals and birds. It has been reported to harbor lots of viruses, some of which can cause human diseases, including JMTV, TBEV, lymphocytic choriomeningitis virus, Powassan virus (POWV), Nairobi sheep disease viruses (NSDVs), SFTSV, and SGLV (Qin et al. 2014; Gong et al. 2015; Yun et al. 2016; Zhang et al. 2018; Ma et al. 2021), as well as some others involving Yamaguchi virus, Khasan virus, Thogoto virus, Dabieshan tick virus (DBSTV), Huangpi Tick Virus 2, and Imakoba tick virus (Leonova et al. 2009; Li et al. 2015; Luo et al. 2015; Shimoda et al. 2019). Rhipicephalus microplus has a one-host life cycle and is an economically important tick that parasitizes a variety of livestock species, particularly cattle. Likewise, a variety of viruses have also been discovered from R. microplus involving SFTSV, JMTV, Lihan tick virus (LHTV), Wuhan tick virus 1 (WHTV1), Wuhan tick virus 2 (WHTV2), Rhipicephalus associated flavi-like virus, and cattle tick tymovirus-like virus 1 (Barker and Walker 2014; Li et al. 2015; Souza et al. 2018; Sameroff et al. 2019; Gondard et al. 2020).

Herein, we collected H. longicornis and R. microplus ticks in northern Hubei Province, China, where emerging or reemerging viral tick-borne diseases have been reported, including SFTS and Nairobi sheep disease. Analyses of RNA-sequencing data at the viral species level gave more insight into the virus diversity and evolutionary history in ticks in Hubei Province. Further mining the viromic data at different ecological niches demonstrated that host species and feeding status might affect the virome structure in ticks.

Materials and methods

Sample collection

From June to July in 2019, a total of 1,024 ticks were collected from eight towns in northern Hubei Province, China, including Wanhe (n = 90) and Yindian (n = 320) in Suizhou City, and Huahe (n = 190), Yonghe (n = 40), Shunhe (n = 108), Yantianhe (n = 106), Baimiaohe (n = 140), and Shengli (n = 30) in Huanggang City (Fig. 1). The samples have been stored in dry ice before being transported to the laboratory. Tick species was morphologically identified by a trained expert in the field and was subsequently confirmed by amplifying and sequencing a 710-bp amplicon of the mitochondrial cytochrome c oxidase subunit I gene as previously described (Folmer et al. 1994). The 410 ticks captured from Suizhou City were characterized as H. longicornis, with all feeding on goats (Capra hircas) (Fig. 1). Tick species collected in Huanggang City included both H. longicornis and R. microplus, of which the 60 H. longicornis ticks sampled in Shunhe were collected on the grassland (free), 160 R. microplus ticks in Huahe fed on cattle (Bos taurus), and the other 304 H. longicornis and 90 R. microplus ticks all fed on goats (Fig. 1).

Figure 1.

Sampling sites in Hubei Province and the tick species composition. Abbreviations of sampling sites: Wanhe (WH) and Yindian (YD) towns in Suizhou City (SZ), and Huahe (HH), Yonghe (YH), Shunhe (SH), Yantianhe (YTH), Baimiaohe (BMH), and Shengli (SL) towns in Huanggang City (HG). The number of samples at each site was marked. The types of ticks are represented by different colors: blue: H. longicornis fed on goats (H. longicornis—goat), yellow: H. longicornis free (H. longicornis—free), green: R. microplus fed on goats (R. microplus—goat), and orange: R. microplus fed on cattle (R. microplus—cattle).

RNA library construction and sequencing

Ticks were washed using sterile, RNA, and DNA-free phosphate buffered saline (PBS) solution (Gibco) individually, and were pooled (20–30 individuals per pool) based on the collection location and host species (Supplementary Table S1). They were then cut into pieces and homogenized thoroughly in PBS solution, and total RNA was extracted using the RNeasy Plus Mini Kit (QIAGEN). Libraries were constructed with the MGIEasy mRNA Library Prep Kit (BGI) according to the BGI mRNA Library Preparation protocol (https://www.bgi.com). The resultant libraries were quality-controlled using an Agilent 2100 Bioanalyzer (Agilent) and pooled in equal quantity. Paired-end (100 bp) sequencing of each RNA library was performed on the BGISEQ500RS platform (BGI).

Transcriptome analysis and virus discovery

The sequencing reads (raw data) were quality-controlled and preprocessed using Fastp v0.20.0 (Chen et al. 2018). All of the ribosomal RNA (rRNA) reads were removed with Bowtie2 v2.3.3.1 (Langmead and Salzberg 2012). The resultant reads (non-rRNA) were then de novo assembled using Trinity v2.5.1 with default settings (Grabherr et al. 2011). The assembled contigs were compared with the non-redundant nucleotide (nt) and protein (nr) database downloaded from GenBank using Blastn with an e-value cutoff at 1E-1 and Diamond Blast with an e-value cutoff at 1E-5 (Li et al. 2021). All contigs were filtered to remove the host, plant, bacterial, and fungal sequences. Viral contigs were further filtered to remove viruses of bacteria, archaea, plant, and other eukaryotes by searching the Virus-Host DB database (http://www.genome.jp/virushostdb/), as well as retroviruses. For each virus in an individual library, the contigs were further merged using Geneious Prime v2020.0.4 (Biomatters).

The relative abundance of the identified viruses was determined by mapping the reads back to the assembled contigs using Bowtie2 v2.3.3.1 and was calculated as reads per million (RPM) after the removal of rRNA reads.

Viral genome confirmation and annotation

To exclude the probable false positives caused by barcode hopping, the viruses with high positive rates (identified in more than one-third of the libraries) were detected in all of the original samples using nested reverse transcription-polymerase chain reaction (RT-PCR). The integrity of the viral genomes was then verified by comparing the contigs of each virus with their closest reference(s) from GenBank. To further confirm the viral genomes in the original samples, specific RT-PCR primers were designed based on the contigs to amplify the full-length genomes. The positive amplicons were then sequenced using Sanger sequencing (Sangon Biotech), and the full-length viral genomes were assembled. Detailed information on the primers for the detection or whole genome amplification are shown in Supplementary Table S2.

Potential open reading frames (ORFs) in the viral sequences were predicted using ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/) with the minimal length of 100 amino acids (aa) and were further compared with reference sequences. The conserved domains present in the sequences were annotated using Conserved Domain Search Service (CD Search) (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi).

Virus classification

The viruses discovered here were classified and named according to the latest International Committee on Taxonomy of Viruses (ICTV) report of Virus Taxonomy (https://talk.ictvonline.org/ictv-reports/ictv_online_report/) mainly based on the nt and aa sequence identities. If the species demarcation criteria remain unclear within a genus, a novel viral species would also be defined if it holds <80 per cent nt identity across the complete genome or <90 per cent aa identity of the RNA-dependent RNA Polymerase (RdRP) domain with known viruses (Pettersson et al. 2020; Wille et al. 2020a). All the novel viruses were named as ‘Hubei tick’, followed by common viral names according to their taxonomy. For viruses that were classified into established taxonomies, they were marked with ‘Hb-2019’ to distinguish them from described strains.

Phylogenetic analyses

To infer the phylogenetic relationships of the viruses identified here, representative reference sequences were downloaded from GenBank and aligned using MAFFT v7.450 with default settings (Katoh and Standley 2013). The alignment was then trimmed using the TrimAl program (Capella-Gutiérrez, Silla-Martínez, and Gabaldón 2009). The best-fit substitution model for each alignment was selected based on the Akaike Information Criterion of Smart Model Selection (Lefort, Longueville, and Gascuel 2017). Phylogenetic trees were then estimated using the maximum likelihood method in PhyML v3.0 with 1,000 bootstrap replicates (Guindon et al. 2010).

Virome ecology analyses

Ecological and statistical analyses were performed at the viral species level using R v4.0.2 integrated in RStudio v1.3.1093 and were plotted using ggplot2 package (Wickham 2016) v3.3.3. For each library, virome richness, Shannon, and Shannon effective indices (alpha diversity) were measured using the Rhea script sets (Lagkouvardos et al. 2017) and compared between different tick species as well as the feeding status using the Kruskal–Wallis rank sum test. Beta diversity was calculated using the Bray–Curtis dissimilarity matrix using the Vegan package (Oksanen et al. 2007) v2.5-7. Principal-coordinate analysis was performed based on the Bray–Curtis dissimilarity, and additional cluster analysis was performed with Adonis tests (PERMANOVA) using the Phyloseq package (McMurdie and Holmes 2013).

Results

Identified RNA viruses

A total of fifty RNA libraries were constructed and sequenced. After adapter trimming and quality control, 406.2 GB clean data were yielded. After further removing rRNA, ∼2.8 billion non-rRNA reads were obtained. The number of the non-rRNA reads of each library was summarized in Supplementary Table S1. A total of 4,608 viral contigs were obtained by de novo assembly of ∼17.2 million viral reads, which accounted for 0.61 per cent of the total non-rRNA reads. Within each library, the viral reads comprised from 0.004 per cent (library MCSH1) to 12.4 per cent (library HAHH2) of the total non-rRNA reads. After being aligned by Blast and filtered by Virus-Host DB, the viral contigs were finally annotated to twelve viruses belonging to six families: Flaviviridae, Matonaviridae with positive sense single-stranded RNA ((+)ssRNA) genomes, Peribunyaviridae, Nairoviridae, Phenuiviridae, and Rhabdoviridae with negative sense ssRNA ((–)ssRNA) genomes (Table 1).

Table 1.

Viruses identified in the present study.

Classification	Virus (abbreviation)	Genome (bp)	Closet relative (%nt identity)	No.a
Flaviviridae
Jingmenvirus like	Jingmen tick virus (JMTV) Hb-2019	3,110/2,817/2,824/2,790	JMTV SY84 (94.0/95.7/94.4/95.3)	25 (1)
Pestivirus like	Hubei tick flavivirus (HTFV)	17,253	Bole tick virus 4 (78.6)	1 (1)
Matonaviridae
Matonavirus like	Hubei tick hepe-like virus (HTHLV)	5,547	Tick-borne tetravirus-like virus (66.0)	25 (3)
Peribunyaviridae
Peribunyavirus like	Hubei tick peribunyavirus (HTPV)	12,096/4,859/2,180	Wenzhou tick virus TS1-2 (71.1/66.3/68.1)	7 (2)
	Huangpi tick virus 1 (HPTV1) Hb-2019	11,950/5,258/2,492	HPTV1 H124-1 (99.1/97.9/93.3)	2 (2)
Nairoviridae
Orthonairovirus	Nairobi sheep disease virus (NSDV) Hb-2019	12,018/5,053/1,537	NSDV Hubei (96.7/95.2/96.8)	1 (1)
Phenuiviridae
Bandavirus	Severe fever with thrombocytopenia virus (SFTSV) Hb-2019	6,347/3,359/1,709	SFTSV HB2016-13 (99.8/99.8/99.8)	2 (2)
Uukuvirus	Dabieshan tick virus (DBSTV) Hb-2019	6,547/1,766	DBSTV D3 (99.1/98.9)	13 (2)
	Lihan tick virus (LHTV) Hb-2019	6,500/1,807	LHTV LH-1 (98.9/98.5)	8 (2)
Rhabdoviridae
Alphanemrhavirus like	Wuhan tick virus 1 (WHTV1) Hb-2019	10,313	WHTV1 X78-2 (99.9)	13 (3)
Rhabdovirus like	Wuhan tick virus 2 (WHTV2) Hb-2019	11,398	WHTV2 X78-1 (99.1)	13 (3)
	Hubei tick rhabdovirus (HTRV)	11,740	Bole tick virus 2 BL076 (66.5)	7 (3)

Number of libraries positive to the virus, with the number in parentheses indicating the whole genomes confirmed by RT-PCR and Sanger sequencing.

Viral genomic organization and taxonomy

For each of the twelve RNA viruses, the whole genomes (partial 5ʹ- and 3ʹ-termini) of the representative strains were further verified by Sanger sequencing (Table 1 and Supplementary Fig. S1). Of them, four viruses were proposed as novel viral species as they were highly divergent to any of previously identified viruses and were designated as Hubei tick flavivirus (HTFV), Hubei tick hepe-like virus (HTHLV), Hubei tick peribunyavirus (HTPV), and Hubei tick rhabdovirus (HTRV), respectively (Table 1). Another four viruses including JMTV Hb-2019, Huangpi tick virus 1 (HPTV1) Hb-2019, WHTV1 Hb-2019, and WHTV2 Hb-2019 showed the closest relationships with previously described tick-associated viruses that had not been approved by ICTV yet (Table 1). The remaining four viruses were all classified into known species as they showed close relationships and identical genome organizations with ICTV-approved viruses: NSDV Hb-2019, SFTSV Hb-2019, LHTV Hb-2019, and DBSTV Hb-2019, respectively (Table 1).

Virome composition and abundance

Across the fifty libraries, forty-eight libraries possessed by one to five (virus count) species of viruses, except for libraries MCSH1 and MCSH5 that had no virus detected (Fig. 2A). The twelve libraries of R. microplus ticks possessed high virus counts of three to five, with the highest (five) found in the library HAHH8 comprised with JMTV, HTFV, HPTV1, WHTV1, and WHTV2. Except for the library SZWH3 that had four viral species identified, the other thirty-seven libraries of H. longicornis ticks possessed zero to three viral species. Among the twelve viruses (Table 1), both JMTV and HTHLV were identified in twenty-five of the fifty libraries from both tick species sampled from different locations. Interestingly, WHTV1 and WHTV2 were identified in thirteen libraries concurrently with a correlation coefficient of 0.709. However, there were also several viruses only identified in one tick species: HTPV, HPTV1, SFTSV, NSDV, and HTRV in H. longicornis ticks, and HTFV, LHTV, WHTV1, and WHTV2 in R. microplus ticks. It should be noted that both HTFV and NSDV were only detected in one library.

Figure 2.

Viral presence and abundance across the libraries. (A) Virus counts identified in each library. (B) Total and (C) proportions of the RPM (non-rRNA) of the twelve identified viruses shown in different colors.

For each virus, the viral RPM (Fig. 2B) and the proportion to all viral reads in each library were calculated (Fig. 2C). Of them, the libraries of R. microplus ticks fed on cattle collected in Huahe Town (HAHH2–HAHH9) had relatively high viral reads, most of which were of JMTV (Fig. 2B). However, the viral reads of the thirty-eight libraries of H. longicornis ticks were relatively low, and the relatively prominent viruses included DBSTV in Baimiaohe and Yantianhe towns and HTHLV in towns of Shunhe, Wanhe, and Yindian (Fig. 2B and C).

Virus diversity and evolution

(+)ssRNA viruses

Flaviviridae.

There are four approved genera in this family, including Flavivirus, Hepacivirus, Pegivirus, and Pestivirus. The Jingmenvirus group has also been defined recently, comprising a series of related viruses such as JMTV (Qin et al. 2014; Webster et al. 2015; Dinçer et al. 2019; Temmam et al. 2019a), ALSV (Kuivanen et al. 2019; Wang et al. 2019), and Yanggou tick virus (Fig. 3A). JMTV Hb-2019 identified here fell into the Jingmenvirus group with the identity of 94.0 per cent–95.7 per cent (nt) and 97.6 per cent–99.0 per cent (aa) to other JMTVs. In particular, the twenty-five viruses of JMTV Hb-2019 shared 98.7 per cent–100 per cent nt identities to each other and formed an independent clade in all the trees of the four genomic segments of the JMTV (Fig. 3B and Supplementary Fig. S2). Further inspection of the trees revealed that the JMTV Hb-2019 from different sampling locations and hosts was grouped together, indicative of the wide distribution of JMTVs in these regions.

Figure 3.

Phylogenetic analyses of flaviviruses and jingmenviruses. Phylogenetic trees were constructed based on the RdRP protein sequences of representative viruses in the family Flaviviridae (A) and the segment 3 of jingmenviruses (B). In panel A, viruses obtained in ticks here are marked with red-filled circles and highlighted in red, and the closest referenced viruses are also highlighted in bold font. In panel B, the strains described here are shown using the same colors as those in Fig. 1, representing the type of the ticks. The hosts of JMTVs discovered previously are marked as follows: black-filled circles, ticks; grey-filled circles, other arthropods; and grey-filled rectangles, mammals. The accession numbers of the viral sequences used in the trees are shown in Supplementary Table S3.

HTFV was phylogenetically placed into the Pestivirus-like group (defined here) and clustered with Bole tick virus 4 (with nt identities of 78.6 per cent and aa identities of 80.8 per cent) (Fig. 3A), which was first identified in Hyalomma asiaticum ticks in Xinjiang Uygur Autonomous Region, China (Shi et al. 2016b), and later in several other tick species in China, Trinidad and Tobago (Sameroff et al. 2019), and Thailand (Temmam et al. 2019b).

Matonaviridae.

The newly discovered HTHLV fell within the Matonavirus-like group associated with the family Matonaviridae in the order Hepelivirales, sharing the closest yet distant relationship (identities: nt 66.0 per cent and aa 52.1 per cent) to tick-borne tetravirus-like virus strain FI10 that was identified in Dermacentor variabilis in the USA (Tokarz et al. 2014). In addition, it was distantly related to the Rubella virus in the genus Rubivirus (Fig. 4A). The tree generated based on the complete genome sequences of the twenty-five HTHLV viruses displayed that viral sequences detected in the same locations were inclined to cluster together. For instance, the sixteen strains detected in the H. longicornis ticks from goats in Yindian Town (SZYD1–SZYD16) formed a separate clade (Fig. 4B).

Figure 4.

Phylogenetic analyses of hepeliviruses. Phylogenetic trees were constructed based on the RdRP protein sequences of representative viruses in the order Hepelivirales (A) and the HTHLV strains (B). Viruses identified in ticks here are marked with red-filled circles and are formatted in bold red font, and the closest referenced viruses are shown in bold font. The strains described in this study are shown using the same colors as those in Fig. 1, representing the type of the ticks. The accession numbers of the viral sequences used in the trees are shown in Supplementary Table S3.

(–)ssRNA viruses

Peribunyaviridae.

In the RdRP tree, both of the two peribunyaviruses identified here fell within an unclassified clade provisionally designated as Peribunyavirus-like group by us (Fig. 5). HPTV1 Hb-2019 clustered with HPTV1 strain H124-1 with the nt identity of 93.3 per cent–99.1 per cent, and HTPV was distantly related to Wenzhou tick virus strains TS1 and TS2 with low identities of 66.3 per cent–71.1 per cent (nt) and 60.0 per cent–77.4 per cent (aa). These two reference strains were identified in Haemaphysalis ticks that were collected from Huangpi City in Hubei Province and Wenzhou City in Zhejiang Province, respectively (Li et al. 2015). The seven strains of HTPV showed nt identities of 95.3 per cent–99.7 per cent (segment L), 92.1 per cent–99.6 per cent (segment M), and 90.2 per cent–99.8 per cent (segment S) to each other and were grouped together according to sampling sites (Supplementary Fig. S3).

Figure 5.

Phylogenetic tree generated based on the RdRP protein sequences of peribunyaviruses. Viruses identified in ticks here are marked with red-filled circles and are formatted in bold red font, and the closest referenced viruses are shown in bold font. The accession numbers of the viral sequences used in the tree are shown in Supplementary Table S3.

Nairoviridae.

NSDV Hb-2019 was only identified in one library (SZWH2) of the H. longicornis ticks from Suizhou City and clustered together with NSDV strain Hubei (Fig. 6A), which was identified in the same tick species collected from Suizhou City with nt identities of 96.7 per cent (segment L), 96.8 per cent (segment M), and 95.2 per cent (segment S), respectively (Yang et al. 2019). The tree constructed using the complete S gene segment grouped NSDVs into three clades, consistent with the sampling locations: India, Kenya, and China (Fig. 6B).

Figure 6.

Phylogenetic analyses of nairoviruses. Phylogenetic trees were constructed based on the RdRP protein sequences of viruses in Nairoviridae (A) and S segments of known NSDVs (B). In panel A, viruses obtained in ticks here are marked with red-filled circles and highlighted in red, and the closest referenced viruses are shown in bold font. In panel B, the strains described here are shown using the same colors as those in Fig. 1, representing the type of the ticks. The hosts of NSDVs discovered previously are marked as follows: black-filled circles, ticks; and grey-filled rectangles, mammals. The accession numbers of the viral sequences used in the trees are shown in Supplementary Table S3.

Phenuiviridae.

SFTSV, LHTV, and DBSTV identified here fell within the genera Bandavirus and Uukuvirus, respectively (Fig. 7A). The two viruses of SFTSV Hb-2019 found in this study clustered with other SFTSV strains discovered in humans from Hubei province and shared the closest relationship (nt identities 99.8 per cent) to SFTSV strain HBHG30, which was isolated from the serum of a patient (Wu et al. 2020) (Fig. 7B). LHTV Hb-2019 was only detected in the R. microplus ticks, and the LHTV Hb-2019 strains clustered together showing their closest relationship to the LHTV strain LH-1, which was identified in the R. microplus ticks in China (Li et al. 2015). The Chinese clade was paraphyletic to a clade containing strains from Brazil (Souza et al. 2018) and American countries (Sameroff et al. 2019; Gómez et al. 2020; Gondard et al. 2020) (Fig. 7C). DBSTV 2019-Hb was only identified in the H. longicornis ticks from Baimiaohe and the nearby town of Yantianhe, both of which were located in the Dabie Mountain area where the first DBSTV (strain D3) was found (Li et al. 2015). The thirteen DBSTV identified here shared high nt identities of 98.1 per cent–99.9 per cent (Fig. 7D).

Figure 7.

Phylogenetic analyses of phenuiviruses. Phylogenetic trees were constructed based on the RdRP protein sequences (A) of viruses in Phenuiviridae and S segments of SFTSV found in China (B), LHTVs (C), and DBSTVs (D). In panel A, viruses obtained in ticks here were marked with red-filled circles and highlighted in red, and the closest referenced viruses are also highlighted in bold font. In panel B, the hosts of SFTSVs discovered previously are marked as follows: black-filled circles, ticks; grey-filled circles, other arthropods; black-filled rectangles, humans; and grey-filled rectangles, mammals. The strains described here are shown using the same colors as those in Fig. 1, representing the type of the ticks. The accession numbers of the viral sequences used in the trees are shown in Supplementary Table S3.

Rhabdoviridae.

In the phylogenetic tree of the family Rhabdoviridae, WHTV2 Hb-2019 clustered with WHTV2, and they formed a separate clade that is distantly related to other members in the family, while WHTV1 Hb-2019 and HTRV belonged to an unclassified Alphanemrhavirus-like group and showed the closest relationships to WHTV1 and Bole Tick Virus 2, respectively (Fig. 8A) (Li et al. 2015). Interestingly, however, in both phylogenetic trees of WHTV1 and WHTV2, the viral strains detected from the library SZWH3 from H. longicornis ticks showed distinct nt difference from those from R. microplus ticks (Fig. 8B and C).

Figure 8.

Phylogenetic analyses of rhabdoviruses. Phylogenetic trees were constructed based on the RdRP protein sequences of viruses in Rhabdoviridae (A) and full-length genomes of WHTV2 (B) and WHTV1 (C) found in China. In panel A, viruses identified in ticks here are marked with red-filled circles and are formatted in bold red font, and the closest referenced viruses are shown in bold font. In panel B, the strains described in this study are shown using the same colors as those in Fig. 1, representing the type of the ticks. The accession numbers of the viral sequences used in the trees are shown in Supplementary Table S3.

Ecological factors associated with virome structure

The highest virus richness was found in the R. microplus ticks fed on cattle with a median of 4.0 (Fig. 9A), and both the Shannon (Fig. 9B) and Shannon effective indices (Fig. 9C) were the highest in the R. microplus ticks fed on goats, while the lowest values of the three indices were all from H. longicornis ticks. Furthermore, generally, the libraries of R. microplus ticks had higher diversity at the viral species level than those of H. longicornis ticks, indicating that tick species and feeding status may influence the virome structure. For each virus, the abundance (shown in RPM) among different species of ticks (Supplementary Fig. S4A) and sampling locations (Supplementary Fig. S4B) was further compared, and JMTV showed a significant difference in both aspects.

Figure 9.

Comparison of the viral diversity between tick species and feeding statuses. (A) Virome richness, (B) Shannon, and (C) Shannon effective indices. The numbers of libraries of the four types of ticks are as follows: H. longicornis—goat (n = 35), H. longicornis—unfed (n = 3), R. microplus—cattle (n = 8), and R. microplus—goat (n = 4).

Beta diversity analysis described the clustering of the libraries with similar virome composition and abundance. It is evident that the libraries from different types of ticks were not clustered together (df = 3, P = 0.001) (Fig. 10), further demonstrating that tick species and feeding status may affect the virome structure. Notably, for H. longicornis ticks, libraries from the same or nearby sampling locations tended to be clustered together (df = 7, P = 0.001), suggesting the potential association between geolocation and the virome structure of ticks.

Figure 10.

Beta diversity analysis of the viromic composition among libraries.

Discussion

Herein, we performed a viromic investigation of H. longicornis and R. microplus ticks in Hubei Province, China, and identified twelve RNA viruses belonging to six viral families. Of the four novel viruses, three (HTHLV, HTPV and HTRV) were detected in H. longicornis ticks and HTFV was found in R. microplus ticks. Of the other eight viruses, those pathogenic to human or mammals were identified in the same tick species as previously documented, such as NSDV and SFTSV in H. longicornis ticks (Gong et al. 2015; Yun et al. 2016) and JMTV in both species (Qin et al. 2014; Temmam et al. 2019). Notably, HPTV1, DBSTV, WHTV1, and WHTV2 were first discovered in H. longicornis ticks. These results exhibited the virus diversity in the two tick species and expanded the host spectrum of some known viral members, highlighting the necessity of increased surveillance of viral pathogens in ticks in this area.

Phylogenetic analyses provided a comprehensive evolutionary relationship of these tick-associated viruses. Some viruses, including WHTV1 and WHTV2, may co-evolve with their hosts, and HTHLV and LHTV tend to evolve locally as strains from different locations clustered separately. Other viruses, e.g. JMTV, exhibit frequent cross-host transmission even in the relatively isolated mountainous areas. Since the first report of JMTV in R. microplus ticks in China (Qin et al. 2014), JMTV has been subsequently identified in different countries (Ladner et al. 2016; Souza et al. 2018; Dinçer et al. 2019; Guo et al. 2020) and was found to be related to human infection in Kosovo and China, respectively (Emmerich et al. 2018; Jia et al. 2019). Till date, various JMTVs have been discovered in a broad spectrum of hosts involving ticks, cattle, bats, rodents, primates, and humans, with a wide distribution in Asia, Africa, Europe, Central America, and South America. The newly found twenty-five JMTVs in ticks here formed an independent phylogenetic clade, paraphyletic to strains previously identified in ticks and mammals, indicating the increasing genetic diversity of JMTVs.

The transmission and infection of tick-borne viruses mostly occur when feeding and are activated by molecules in tick saliva, such as TBEV and POWV (Labuda et al. 1993; Hermance and Thangamani 2015; Šimo et al. 2017). It is generally assumed that multihost ticks pose more threats to the public due to their capacities to transmit pathogens when changing hosts. However, our study also demonstrated that R. microplus ticks (one-host life cycle) had relatively higher diversity and abundance of viruses, although their vector capacity and relevance to human and animal health deserve further investigation.

As the dominant species and important viral hosts in China, H. longicornis and R. microplus ticks are excellent model species for understanding the ecological factors that influence the virome structure. Herein, the viral composition of the two tick species differed, indicating that tick species may be associated with the virome structure. Meanwhile, the overwhelmingly high abundance of JMTVs in R. microplus ticks fed on cattle suggested a correlation between viral richness and feeding status of ticks (Supplementary Fig. S4A). Considering that viruses could be transmitted among ticks and their mammalian hosts by hemophagia, we conjecture that virome structure might undergo a dynamic change during the life cycles of ticks. However, for most of the sampling sites, only samples of a single type of ticks were collected, and therefore, the correlation between geolocation and virome structure clearly warrants further study. However, our study revealed that multiple factors might be involved in shaping the virome structure and the genetic diversity of viruses, highlighting the complexity of virome composition and the greatly underestimated species and genetic diversity of viruses in nature.

However, there are still many challenges for virus–host ecological studies: (1) viruses are discretely excreted from hosts and heterogeneously distributed, and therefore, a large sampling size is necessary to more accurately characterize the tick-borne virome structures (Plowright et al. 2016); (2) it is evident that there are still a large number of unknown viruses in nature, while the identification and classification of highly divergent viruses remain a big challenge (Geoghegan and Holmes 2017); and (3) many environmental factors can influence the virus–host pattern even at the host population level (Bergner et al. 2020), and thus the same host species sampled from different geographic areas sometimes cannot be compared directly due to the confounding environmental variables.

The novel viral species and variants described here provided new insights into the tick virome structures that exhibited the complex evolutionary history and the potential ecological factors and highlighted the necessity of surveillance of tick-borne viruses. In addition, as the pathogenicity of these tick-borne viruses remains less understood, future studies should include the establishment of the infection and transmission models to develop the prevention and control strategies against tick-associated viral diseases.

Click here for additional data file.