Ticks feeding on patients who sought treatment at Mudanjiang Forestry Central Hospital in Heilongjiang Province, China, were collected between May and August 2016. The study was approved by the ethics committee of the hospital (Mudanjiang Forest Central Hospital 2011-03) in accordance with the medical research regulations of China. The species and developmental stage were identified by an entomologist (Yi Sun). Adult I. persulcatus ticks were included in the study, and the species was later confirmed by analyzing sequences of the mitochondrial cytochrome c oxidase subunit I gene. We further divided the feeding ticks into three groups: flat (ticks that had imbibed small amounts of blood and had not started to enlarge), partially fed (ticks that had begun to enlarge but were not fully engorged), and fully fed (ticks that had fully engorged and were almost ready to detach). All samples were captured alive and thoroughly surface-sterilized (two successive washes in 70% ethanol for 30 s each) before storage at −80°C for DNA/RNA extraction.

Extraction of total DNA and RNA from ticks was performed using the AllPrep DNA/RNA Mini Kit (Qiagen, Valencia, CA, United States) with modifications. Briefly, ticks were homogenized in RLT solution under liquid nitrogen. The homogenate was then incubated at 55°C for 10 min with proteinase K (Qiagen) and centrifuged for 30 s at 15,000 × g. The homogenized lysate was transferred to an AllPrep DNA spin column and centrifuged for 30 s at 8,000 × g. The AllPrep DNA spin column was used for later DNA purification, and the flow-through was used for RNA purification as per the manufacturer’s instructions. Ten ticks were pooled together to extract RNA/DNA, and each group comprised three pools.

The extracted RNA was used for transcriptome sequencing (RNA-seq) after RNA quantification and quality checking. The rRNA was removed using RiBo-Zero Gold rRNA removal reagents (human/mouse/rat) (Illumina). Then, the sequencing library was prepared following the Illumina standard protocol. Paired-end (2 × 150 bp) sequencing of the RNA library was performed on an Illumina HiSeq 4000 platform at Novogene Tech (Beijing, China).

RNA sequencing reads were de novo assembled using the Trinity program (Grabherr et al., 2011). Assembled contigs were subject to BLASTx against all non-redundant (nr) database downloaded from GenBank, and the threshold E value was set to 1e-5. Putative viral contigs were further merged by high-identity overlaps using the SeqMan program of Lasergene package v7.1 (DNAstar, Madison, WI, United States). Original reads were aligned to the contigs again using Bowtie2 (Langmead and Salzberg, 2012) to complete the remaining gaps, and the assembly was verified in the Integrated Genomics Viewer (Thorvaldsdóttir et al., 2013).

The highly conserved RNA-dependent RNA polymerase (RdRp) gene was used to construct family-level phylogenies comprising the viruses that were highly divergent. Predicted viral proteins of the RdRp genes in this study were aligned with reference proteins of the same viral families using the E-INS-i algorithm in MAFFT version 7 (Katoh and Standley, 2013). Ambiguously aligned regions were removed using TrimAl (Capella-Gutiérrez et al., 2009). The WAG + Γ model was identified as the best-fit amino acid substitution model using Prot-Test 3.4 (Darriba et al., 2011), which was used for maximum likelihood (ML) phylogeny reconstruction with bootstrap tests (1,000 replicates) in PhyML version 3.0 (Guindon and Gascuel, 2003). The ML trees were visualized with FigTree v1.4.2.

A single de novo assembly across all samples was generated using Trinity (Grabherr et al., 2011), then reads from each sample were separately aligned back to the single Trinity assembly for downstream analyses of the differential expression of transcripts. Transcript abundance containing RNA-Seq fragment counts for each transcript (or gene) across each sample was estimated using the alignment-based abundance estimation method RSEM (Li et al., 2009). The trimmed mean of M-values normalization method (TMM) were subsequently used to normalize the transcript abundance. Differentially expressed (DE) transcripts or genes were identified by Fisher’s oexact test running “edge” R packages (Robinson et al., 2010; McCarthy et al., 2012), and then those transcripts or genes that had p-values of at most 5e-2 and were at least 2^1.5-fold differentially expressed were extracted to generate a sample correlation matrix heatmap and a DE gene vs. samples heatmap using Trinity downstream analysis script. The DE genes shown in the heatmap were grouped into gene clusters with similar expression patterns by cutting the hierarchically clustered gene tree at 60% from the tip. The functional annotations of the DE transcripts or genes were performed using Trinotate (Bryant et al., 2017) based on Pfam database. We classified these annotated transcripts or genes based on the COG database using eggnog mapper (Huerta-Cepas et al., 2017).

A total of 90 feeding adult I. persulcatus ticks were collected from patients. We sequenced nine pools of ticks (10 ticks/pool) by transcriptomic survey. Pools TG1, TG2, and TG13 comprised flat ticks; pools TG3, TG4, and TG14 comprised partially fed ticks; and pools TG5, TG6, and TG15 comprised fully fed ticks. We analyzed cox1 genes from RNA-seq data to identify the source of previous bloodmeals of I. persulcatus and found five mammal host sources, namely, Tonatia saurophila, Mus musculus, Ovis aries, Sus scrofa, and Homo sapiens (Figure 1). The host cox1 gene sequences assembled in the study had over 98.5% identity to the respective reference sequences (Supplementary Table 1). The most abundant cox1 sequence was from Homo sapiens, accounting for 1,570 reads across all samples. Intriguingly, one bat species, Tonatia saurophila, was detected with a read count of 48 in one of the partially fed tick pools (Figure 1 and Supplementary Table 1).

We performed nine RNA-seq runs on each of the nine pools of I. persulcatus ticks, generating 72 GB of total sequence data. Among the nine resultant RNA-seq libraries, library TG13 (flat ticks) had the most abundant viral reads, followed by library TG4 (partially fed ticks), while libraries TG6 and TG15 (fully fed ticks) had the fewest viral reads (Figure 2). Through a BLASTx search with the assembled genome sequences, we discovered 14 known or novel viruses belonging to 11 viral families (Table 1). All of which had fully or almost fully intact genomes with undisrupted reading frames. These included, with proposed names of new viruses in brackets, one Mononegavirales-like virus (I. persulcatus mononegavirales-like virus HLJ), one novel chuvirus (I. persulcatus chuvirus HLJ), one novel rhabdovirus (I. persulcatus dimarhabdovirus HLJ), three distant relatives of Totiviridae (I. persulcatus totivirus 1 HLJ, I. persulcatus totivirus 2 HLJ and I. persulcatus totivirus 3 HLJ), one novel Partitiviridae (I. persulcatus picobirna virus HLJ), one novel orbivirus (I. persulcatus orbivirus-like HLJ), one Luteoviridae-like virus (I. persulcatus luteo-like virus HLJ), four distant relatives of Bunyavirales (I. persulcatus phlebovirus 1 HLJ, I. persulcatus phlebovirus 2 HLJ, I. persulcatus nairovirus HLJ, I. persulcatus bronnoya virus HLJ), and Jingmen tick virus (JMTV). Reads from I. persulcatus nairovirus HLJ were most abundant, followed by reads from I. persulcatus dimarhabdovirus HLJ (Figure 2). These viruses shared 36–99% amino acid identities compared with previously reported viruses, and the query coverage varied from around 30 to 90% (Table 1).

We identified three viruses in the Mononegavirales. Among them, I. persulcatus mononegavirales-like virus HLJ was closely related to deer tick mononegavirales-like virus, which had been previously reported in Heilongjiang Province, China (Meng et al., 2019; Figure 3A). Ixodes persulcatus dimarhabdovirus HLJ clustered with Long Island tick rhabdovirus and formed a monophyletic clade with other rhabdoviruses (Figure 3A). Ixodes persulcatus chuvirus HLJ had 83% amino acid similarity with blacklegged tick chuvirus-2 but only 57% coverage of the L segment sequence. We obtained the nearly complete genomes of these three viruses. The genome of I. persulcatus mononegavirales-like virus HLJ, I. persulcatus dimarhabdovirus HLJ, and I. persulcatus chuvirus HLJ were 10,865, 10,365, and 11,440 nucleotides (nt) long, and all encoded four predicted ORFs.

Three distant relatives of Totiviridae (I. persulcatus totivirus 1 HLJ, I. persulcatus totivirus 2 HLJ, I. persulcatus totivirus 3 HLJ) had extremely low amino acid identity (37–40%) with previously described viruses. Their RdRp sequences were highly divergent from the existing clades and fell into the unclassified groups on the phylogenetic tree (Figure 3C). This indicates the distinct phylogenetic position of these newly identified viruses from the currently known viruses.

Three members of the Bunyavirales, I. persulcatus phlebovirus 1 HLJ, I. persulcatus phlebovirus 2 HLJ, and I. persulcatus nairovirus HLJ showed high amino acid identities of 75, 87, and 99% with previously reported blacklegged tick phlebovirus 3, Norway phlebovirus 1, and Beiji nairovirus, respectively, based on the RdRp gene. The fourth member of the Bunyavirales, I. persulcatus bronnoya virus HLJ exhibited 38% amino acid identity with Bronnoya virus (88% sequence coverage of RdRp gene). The phylogenetic tree of Bunyavirales indicated that the two phleboviruses and the nairovirus clustered with previously reported Chinese strains and formed well-supported monophyletic groups closely related to the “classic” phleboviruses and nairoviruses (Figure 3D). It is suggested that they may share a single common ancestor for this gene. The newly discovered Bronnoya virus RdRp sequence clustered with the previously reported Bronnoya virus.

Ixodes persulcatus luteo-like virus HLJ, I. persulcatus picobirna virus HLJ, and JMTV found in this study exhibited 87, 94, and 99% amino acid identities to previously reported viruses with, respectively, 98, 91, and 93% sequence coverage (Figures 3B,E,F).

Ixodes persulcatus orbivirus-like HLJ revealed 40% amino acid similarity across 94% of the VP1 gene sequence of the orbivirus Skunk River virus (MK100569.1), and it showed a phylogenetic position distinct from other orbiviruses in the phylogenetic tree (Figure 3G).

We next detected the differential expression responses of ticks at two different stages of feeding: flat (pools TG1, TG2, and TG13) and fully fed (pools TG5, TG6, and TG15) to investigate the effect of the bloodmeal. We first generated a single Trinity assembly across all pools and then calculated the abundance separately for each pool. The subsequent differential expression analysis showed that a total of 1,302 genes were differently expressed (p < 0.05, logFC > 1.5). Most of the DE genes were negatively correlated (Supplementary Figure 1). As shown in Figure 4A, most of the DE genes (71.4%) were downregulated in fully fed ticks; of these, ATP synthase c-subunit, a putative cuticular protein, ribosomal protein L6, and the putative apoptosis-promoting RNA-binding protein TIA-1/TIAR were the most downregulated (Supplementary Table 2). The gene encoding a putative secreted salivary gland peptide was most upregulated at engorgement, followed by a putative cystatin, putative coiled-coil domain-containing protein, and putative secreted glycine-rich protein (Supplementary Table 2). One antimicrobial peptide, lysozyme, was downregulated in fully fed ticks (Supplementary Table 3).

All DE genes shown in the heatmap (Figure 4A) were grouped into clusters of genes with similar expression patterns. A total of 251 transcripts were clustered in the upregulation category (Figure 4B), while 930 transcripts were clustered in the downregulation category (Figure 4C). According to the taxonomy classification based on the COG database, 18.7% of DE genes in the upregulated category were related to post-translational modification, protein turnover, and chaperones, and 12.7% of upregulated DE genes displayed functions of translation, ribosomal structure, and biogenesis (Supplementary Figure 2A). Most of the DE genes in the downregulated category were related to functions of nucleotide transport and metabolism (18.2%) and cytoskeleton (18.3%) (Supplementary Figure 2B).

Reactive Oxygen Species (ROS) produced by the mosquito midgut play a role in maintenance of redox homeostasis (Citelli et al., 2007; Hernandez et al., 2019), and previous literatures also suggested that the balance between positive and negative aspects of ROS in ticks should be carefully maintained (Hernandez et al., 2019). To achieve redox homeostasis, the ticks have a complex anti-oxidant system composed of enzymes such as glutathione peroxidase, glutathione-S transferase, and so on (Hernandez et al., 2018, 2020; Kusakisako et al., 2018, 2020). We found that antioxidant genes including glutathione peroxidase, glutathione S-transferase, peroxidase, peroxiredoxin, and protein kinase C were all downregulated in the fully fed ticks compared with the flat ticks (Supplementary Table 3).