Identification of Novel Viruses in Amblyomma americanum, Dermacentor variabilis, and Ixodes scapularis Ticks

The incidence of tick-borne disease is increasing, driven by rapid geographical expansion of ticks and the discovery of new tick-associated pathogens. The examination of the tick microbiome is essential in order to understand the relationship between microbes and their tick hosts and to facilitate the identification of new tick-borne pathogens. Genomic analyses using unbiased high-throughput sequencing platforms have proven valuable for investigations of tick bacterial diversity, but the examination of tick viromes has historically not been well explored. By performing a comprehensive virome analysis of the three primary tick species associated with human disease in the United States, we gained substantial insight into tick virome diversity and can begin to assess a potential role of these viruses in the tick life cycle.

IMPORTANCE The incidence of tick-borne disease is increasing, driven by rapid geographical expansion of ticks and the discovery of new tick-associated pathogens. The examination of the tick microbiome is essential in order to understand the relationship between microbes and their tick hosts and to facilitate the identification of new tick-borne pathogens. Genomic analyses using unbiased high-throughput sequencing platforms have proven valuable for investigations of tick bacterial diversity, but the examination of tick viromes has historically not been well explored. By performing a comprehensive virome analysis of the three primary tick species associated with human disease in the United States, we gained substantial insight into tick virome diversity and can begin to assess a potential role of these viruses in the tick life cycle. KEYWORDS I. scapularis, ticks, virome T icks (order Ixodida) are hematophagous arthropods that are second only to mosquitos as vectors of human pathogens worldwide. The ability of ticks to transmit a wide range of microbial pathogens, combined with their promiscuous feeding and geographical range expansion, makes them a substantial threat to animal and human health (1). Several recent studies have described novel tick-transmitted pathogens, suggesting that the full spectrum of tick-borne pathogens has not been identified (2)(3)(4)(5).
Genomic analyses using unbiased high-throughput sequencing (UHTS) platforms have proven valuable for investigations of microbially diverse populations, and recently these analyses have been extended to include ticks. The characterization of the tick microbiome is fundamental to understanding the relationship between microbes and their tick hosts and may facilitate the discovery of novel pathogens. UHTS has enabled the examination of tick viromes, a component of the tick microbiome that historically has not been well explored (6)(7)(8)(9)(10)(11)(12). In 2014, we reported the discovery of nine new viruses that were detected in Ixodes scapularis (the blacklegged tick), Dermacentor variabilis (the American dog tick), and Amblyomma americanum (the lone star tick) collected at a single site in New York State (6,13). Several of these viruses were also identified in a recent study of I. scapularis ticks from Pennsylvania (14). These investigations have provided an insight into the breadth of viral genetic diversity found in ticks in the northeastern United States. In this study, we expanded on our previous work by examining a greater population of these three tick species collected from multiple geographically diverse sites in New York, Connecticut, and Virginia. blacklegged tick phleboviruses 1 and 2 (BTPV-1 and -2) (6). All three of these viruses were present in every I. scapularis pool examined in this study. We also identified one additional phlebovirus-like virus with sequence similarity to BTPV-1 and -2, provisionally designated blacklegged tick phlebovirus 3 (BTPV-3). We obtained the complete coding sequence for both the L and S segments of this virus. However, as was the case with BTPV-1 and -2, a third M-like segment could not be identified for BTPV-3. The S segment of BTPV-3 contains a 1,404-nucleotide (nt) putative nucleocapsid open reading frame (ORF), whereas the L segment encodes a 6704-nt putative RNA-dependent RNA polymerase (RdRp). The amino acid sequences of BTPV-3 N and L proteins were approximately 30% and 50% similar to the S and L proteins of BTPV-1 and -2. Phylogenetically, BLTPV-3 clusters with BTPV-1, BTPV-2, and Norway phlebovirus 1, a virus recently identified in Ixodes ricinus ticks in Europe (12) (Fig. 2). BLTPV-3 sequences were present in 65% of I. scapularis pools. We also examined 23 individual ticks from one of the positive pools, with 8 ticks testing positive, indicating this virus has a lower prevalence than BTPV-1 and -2.
Three pools (8%) of I. scapularis adults collected at sites 4, 9, and 10 on Long Island (Table 1) were positive for a highly divergent Bunyavirales-like virus that we provisionally named Laurel Lake virus (LLV). The genome of LLV contained three segments: a 1.2-kb segment (RNA1) encoding a putative nucleocapsid, a 2.5-kb segment (RNA2)   I. scapularis ticks from New London, CT, collected in 2010. Only two ticks were positive-both from New York. These results, combined with LLV being present in Ͻ10% of UHTS pools, suggest this virus is rare in I. scapularis.
We also identified sequences of two distinct contigs from I. scapularis with homology to the RdRp of newly discovered unclassified viruses that cluster within Bunyavirales. We provisionally named these Ixodes scapularis-associated viruses 5 and 6 (ISAV-5 and -6). The sequence of ISAV-6 was most similar to Rhizoctonia solani negativestranded virus 4 (GenBank accession no. KP900923), while ISAV-5 clustered with Fusarium poae negative-stranded virus 2 (accession no. LC150619) (Fig. 2)   of fungi, although other segments were not identified. ISAV-5 and -6 were both present in the same two I. scapularis pools. We subsequently tested the 24 individual ticks that comprised one of the positive pools and detected both viruses in the same two ticks. Our data suggest that both viruses are likely acquired simultaneously by I. scapularis.
Rhabdoviridae. We identified four new viruses with sequence homology to viruses in the family Rhabdoviridae: two in I. scapularis and two in D. variabilis ( Table 3). The two I. scapularis-borne viruses were provisionally named blacklegged tick rhabdovirus 1 (BLTRV-1) and New Kent County virus (NKCV). BLTRV-1 was present in two pools of adult ticks (one each from Connecticut and New York) and phylogenetically clustered with rhabdoviruses recently identified by UHTS in insects (Fig. 3). We tested 22 individual ticks that comprised one of the BLTRV-1-positive pools by reverse transcription-PCR (RT-PCR), with only one tick testing positive. NKCV was identified in two of the three pools of I. scapularis ticks collected in Virginia, both from New Kent County. None of the Connecticut and New York pools were positive for this virus. Analysis of the NKCV genome revealed that it is similar in structure and sequence to rhabdoviruses classified within the genus Ephemerovirus. NKCV encodes a nonstructural glycoprotein and 5 accessory proteins upstream of the polymerase, a feature of other ephemeroviruses (Fig. 4). Phylogenetic analysis of the L protein revealed that NKCV clustered with Kotokan, Koolpinyah, and Yata viruses, which make up one of the two clades within Ephemerovirus (Fig. 3).
The two rhabdoviruses identified in D. variabilis were provisionally named American dog tick rhabodoviruses 1 and 2 (ADTRV-1 and -2). ADTRV-1 was present in two ADTRV-2 sequences were present in three out of eight D. variabilis pools sequenced. However, due to insufficient sequence coverage, we only assembled the nucleocapsid and the polymerase portions of the genome. Analysis of the amino acid sequence of both proteins revealed that ADTRV-2 sequences were most similar to those of Taishun tick virus and Huangpi tick virus 3, which were identified in ticks from China (Fig. 3).
Little is known about this group of viruses, although they likely will also be classified into a distinct genus within Rhabdoviridae.
Chuviruses. We previously reported the discovery of a virus in I. scapularis with homology to Mononegavirales-like RdRp. Here, we report the complete sequence of this virus, provisionally named Suffolk virus (SUFV) (accession no. NC028243). Analysis of the amino acid sequence and genome structure indicated that SUFV is similar to a group of recently identified arthropod-associated viruses that were assigned to a provisional family designated Chuviridae (7). SUFV was present in I. scapularis ticks in all three geographic areas surveyed, including 11 out of 28 pools from Long Island, 17 out of 18 pools from Connecticut, and 2 out of 3 pools in Virginia.
In addition, we identified two other putative chuviruses, designated blacklegged tick chuvirus 2 (BLTCV-2) and lone star tick chuvirus 1 (LSTCV-1). BLTCV-2 was detected in a single I. scapularis pool from Virginia. LSTCV-1 was present in Ͼ90% of pools from New York and was also present in the one sequenced pool of A. americanum from Virginia.
Analysis of complete genomes revealed all three viruses have a circular genome of approximately 10.8 to 11.5 kb in length that includes four putative ORFs (Fig. 5A). ORFs 1, 2, and 3 encode the RdRp, a putative glycoprotein, and a putative nucleoprotein, respectively. Due to a lack of homology to proteins outside of Chuviridae, the potential function of ORF4 is unknown. Phylogenetically, chuviruses cluster into two distinct groups (Fig. 5B). Clade I includes viruses with circular genomes, and clade II contains viruses with linear or circular multisegmented genomes. All three viruses identified in our study cluster within clade I, although they are more closely related to viruses identified in China than to each other.
Nodaviridae. A new nodavirus, named lone star tick nodavirus 1 (LSTN-1), was present in three A. americanum pools. Nodavirus genomes typically consist of a bipartite, linear, single-stranded positive-sense RNA genome (15). RNA1 encodes the polymerase, and RNA2 typically encodes a polyprotein posttranslationally cleaved into two capsid components. Two genera of nodaviruses are recognized: Alphanodavirus, associated with arthropods, and Betanodavirus, isolated from fish and other marine animals. Phylogenetically, LSTN-1 clusters with arthropod viruses and is the first noda- Characterization of the Tick Virome virus discovered in ticks (Fig. 6). LSTN-1 contains an atypical structure of RNA2, with two ORFs instead of a single polyprotein. The 528-nt 5=-proximal ORF has limited similarities to viral capsid proteins of nodaviruses, while ORF2 has no similarity to any known protein.
Sobemo-and Luteo-like viruses. Recent UHTS studies of invertebrates revealed the existence of a large diversity of Sobemo-like or Luteo-like viruses present across a wide range of arthropod families (7,10). These viruses typically consist of a 3-to 3.5-kb genome that contains two ORFs. We identified four viruses that have sequence homology to viruses within this group. At least one virus was present in pools from each tick species, with one novel virus in A. americanum (lone star tick-associated virus 1), one novel virus in D. variabilis (American dog tick-associated virus 1), and two viruses in I. scapularis we identified previously (Ixodes scapularis-associated viruses 1 and 2) (  Long Island, while one POWV-positive pool of ticks originated from site 1 in Connecticut. All viruses were genetically homogenous (99% nucleotide identity) to the LI-1 strain we identified in 2013. POWV was not present in any of the 10 pools of nymphs or in the three I. scapularis pools of adults from Virginia. By PCR, we screened an additional 70 adult ticks from Virginia that were also negative.
Other viruses. We detected two viruses (90% amino acid identity) in I. scapularis and D. variabilis (Ixodes scapularis-associated virus 4 and American dog tick-associated virus 2) that were similar to viruses identified in Drosophila (Table 3) (21). In lone star ticks, we obtained the complete coding sequence of a new densovirus and totivirus. We also identified several viruses in I. scapularis that were present in only a single pool. These included an avian-like circovirus, an Ilarivirus, a partitivirus, and two picorna-like viruses. All five viruses displayed Ͻ65% amino acid similarity to the next closest virus within their group.

Characterization of the Tick Virome
Since most of the surveyed sites were clustered within New York and Connecticut, we sampled ticks from Virginia to assess variations in the tick virome due to geography. All of the highly prevalent I. scapularis-borne viruses were also found in Virginia. However, despite testing only three pools, we identified two viruses (BLTCV-2 and NKCV) that were not present in any of the 46 I. scapularis pools from New York and Connecticut. Thus, we expect that exploration of additional sites throughout the range of I. scapularis would identify additional viruses.
Transmission. We next examined whether vertical transmission from the female ticks to their eggs may account for the high prevalence of SBV, BLTPV-1, and SUFV in I. scapularis. cDNA generated from 20 individual replete I. scapularis females was tested by RT-PCR for the presence of each virus. We detected a high prevalence (Ͼ50%) of each virus in individual ticks, including nine infected with all three viruses (Table 4). We then examined a portion of the egg mass or larvae from four of these ticks for evidence of these viruses. We detected at least one of the viruses in the eggs/larvae from each of the corresponding parental ticks, confirming that these viruses were transovarially transmitted.

DISCUSSION
The increasing incidence of tick-borne diseases presents a major challenge to public health worldwide. Improving our understanding of the diversity of the tick microbiome will help guide the identification of emerging tick-borne agents. The principal aim of our study was to catalogue viruses present in I. scapularis, D. variabilis, and A. americanum, the primary tick species associated with transmission of human pathogens in the eastern part of the United States. Through this work, we sought to gain insight into the virome diversity of these important tick species, in order to assess any potential role of these viruses in the tick life cycle.
Based on our previous findings and UHTS analyses of other arthropods, the high diversity of viruses uncovered in this study was not unexpected (7,10,21,22). Recent virome studies have revealed that arthropods and other invertebrates are hosts to a wide range of highly divergent viruses. Members of some of these viral groups have been found across a large number of different arthropods but have not been detected  in vertebrates. Although the roles of these viruses in the arthropod life cycles are largely unknown, they are not likely to be horizontally transmitted to other hosts. Although some of the viruses characterized in this study are phylogenetically related to human or animal pathogens, we currently do not have evidence suggesting their potential for horizontal transmission. We also cannot exclude the possibility that some of the identified viral sequences may represent viruses that do not replicate in the tick itself. Instead, these sequences may originate from viruses parasitizing other hosts present within the tick, such as a fungus or nematode, or may represent remains of a virus ingested during a blood meal. The majority of the samples analyzed in our study were comprised of I. scapularis pools. We focused our work on this tick species because of its role as the principal vector of Borrelia burgdorferi, the agent of Lyme disease (23). With an estimated 300,000 annual cases of Lyme disease in the United States, I. scapularis has a substantial impact on public health, and coinfections have the potential to modulate the clinical course of acute Lyme disease or posttreatment persistent illness (24,25). Additionally, sampling success was a factor in the final pooling. At the majority of our collection sites, I. scapularis and/or A. americanum was abundant, whereas D. variabilis was not collected in large numbers or was absent.
A striking feature of the I. scapularis virome was the presence of a large diversity of Bunyavirales-like sequences. Seven of the eight Bunyavirales-like viruses detected in our study were identified in I. scapularis. Four of these (BTPV-1, BTPV-2, BTPV-3, and SBV) share common features. Each virus was present in ticks from all three states and was highly prevalent in I. scapularis, with Ͼ25% of individual adult ticks infected. These viruses appeared to lack the glycoprotein-encoding M segment that prototypical bunyaviruses utilize for cellular attachment. We first reported bunya-like viruses that lacked the M segment in our initial tick virome study, and similar viral sequences have subsequently been identified in a wide array of ticks from Europe and Asia (7,10,12). We also found that transovarial transmission may play a role in maintaining these viruses in their tick hosts. Based on these common features, we propose that BTPV-1, BTPV-2, BTPV-3, and SBV are viral symbionts of I. scapularis. Ticks are known to host a wide array of symbionts, and the presence of transovarially transmitted bacterial symbionts, including many closely related to human pathogens, has been described in all three of the tick species surveyed here (26). Viral symbionts have recently been suggested to exist in several types of arthropods, including ticks (11). Future studies examining horizontal transmission of these viruses will be necessary to confirm our hypothesis.
Due to the high prevalence of these putative viral symbionts, we identified a much greater number of viruses in I. scapularis than A. americanum. Despite examining a similar number of adult ticks (839 I. scapularis versus 720 A. americanum ticks), I. scapularis ticks were hosts to three times more viruses than A. americanum. We found that most A. americanum pools were rarely positive for more than one virus, whereas I. scapularis pools had, on average, Ͼ4 viruses. The primary contributing factor to this difference was the high prevalence of five I. scapularis viruses (SBV, BLTPV-1 and -2, SUFV, and ISAV-1), present in Ͼ75% of pools from Connecticut and New York. All five viruses were also present in I. scapularis from Virginia.
Our study identified four putative novel members of the family Rhabdoviridae, a large group of single-stranded negative-sense RNA viruses (15). Rhabdoviruses infect a broad range of vertebrates, invertebrates, and plants (27). Arthropods are hosts to a wide variety of rhabdoviruses and are known vectors for several pathogenic vertebrate viruses. Based on sequence identity and genome structure, only one of the four viruses identified in our study, the I. scapularis-borne NKCV, can be grouped into an International Committee on Taxonomy of Viruses (ICTV)-recognized genus, Ephemerovirus. Eight other ICTV-recognized ephemerovirus species are recognized. These viruses have been isolated from mosquitos, biting midges, or ruminants from subtropical and tropical regions of Africa, Australia, and Asia (27). Thus, NKCV represents the first tick-associated ephemerovirus and the first to be discovered in North America. Al-though ephemeroviruses have not been associated with human disease, Kontonkan and bovine fever emphemeroviruses have been implicated in bovine disease in Africa (28,29).
The phylogeny of the remaining three viruses is less clear. DTRV-2 and BLTR-2 are both similar to a wide array of viruses identified through UHTS of diverse arthropod species, and little can be discerned about their potential for transmission. DTRV-1 forms a unique cluster with Farmington virus, a virus capable of infecting mammalian cells that was isolated from an unidentified wild bird in Connecticut in 1969 (20). Further studies of DRTV-1, including determination of host tropism, may help clarify if this virus presents a risk for human infection.
POWV is currently the only tick-borne virus associated with human disease circulating in the northeastern United States. POWV was present in Ͻ10% of I. scapularis pools, in agreement with surveillance studies that typically show low prevalence of POWV (Ͻ5%) in individual ticks (30)(31)(32)(33). We detected POWV in pools from both New York and Connecticut, but not in ticks from Virginia. To our knowledge, the virus has not been shown to be present in the southern United States.
In summary, we performed a comprehensive UHTS analysis of tick viromes that substantially expands the current knowledge of viruses associated with I. scapularis, D. variabilis, and A. americanum. The identification and characterization of these viruses will enable future work focused at determining the role of these viruses in the tick life cycle.

MATERIALS AND METHODS
Tick collections and nucleic acid extractions. Ticks (I. scapularis, D. variabilis, and A. americanum) were collected in 2015 and 2016 at 14 sites in New York, Connecticut, and Virginia ( Fig. 1; Table 1). Specimens were kept frozen at Ϫ80°C until processing. All tubes, wash reagents, beads, and media were sterilized by UV radiation prior to processing.
Ticks were sorted by species, pooled (mean n ϭ 22 ticks per pool) (Table 1), and washed (once with 70% ethanol and three times with H 2 O). Each pool was homogenized using the TissueLyser II (Qiagen) in 1 ml of viral transport medium (VTM) with four 5-mm stainless steel beads in an Eppendorf Safe-Lock tube. Homogenates were centrifuged at 8,000 rpm for 1 min to precipitate cellular debris; the supernatant was removed and then purified through a 0.45-m-pore filter (Millipore). Nonencapsulated nucleic acids were digested with RNase A (15 min at room temperature) followed by treatment with Turbo DNase and Benzonase (30 min at room temperature). Total nucleic acids (TNAs) from 300 l of nuclease-treated homogenate were extracted using the easyMAG extraction platform (BioMérieux) and eluted in 40 l, with 11 l used in a reverse transcriptase reaction (Superscript III; Invitrogen) prior to library preparation. For six I. scapularis adult pools, we used a modified method for homogenization, in which individual ticks (n ϭ 20) were washed as described above and individually homogenized in 50 l of VTM. Twenty microliters of the homogenate from each tick was pooled, and the entire pooled homogenate was filtered and nuclease treated as described above. TNAs from 300 l of the pooled homogenate were extracted and eluted in 40 l. The remaining 30 l of homogenates from individual ticks were extracted individually and used for single-agent reverse transcription-PCR (RT-PCR).
Library preparation. Double-stranded DNA was mechanically sheared to an average length of 200 nt and purified. Sequencing was performed on the Illumina HiSeq 4000 system (Illumina, San Diego, CA) using the Hyper Prep kit (KAPA Biosystems, Boston, MA). The demultiplexed FastQ files were adapter trimmed using the cutadapt program (v1.8.3) (34). Adapter trimming was followed by generation of quality reports using FastQC software (v0.11.5), which were used to determine filtering criteria based on the average quality scores of the reads, presence of indeterminate nucleotides, and homopolymeric reads (35). The reads were quality filtered and end trimmed with PRINSEQ software (v0.20.3) (36). Host background levels were determined by mapping filtered reads against a tick reference database using Bowtie2 mapper (v2.2.9) (37). The host-subtracted reads were de novo assembled using the MIRA (4.0) assembler (38). Contigs and unique singletons were subjected to homology search using Megablast against the GenBank nucleotide database. Sequences that showed low or no homology at the nucleotide level were subjected to a BLASTX homology search against the viral GenBank protein database. Sequences from viral BLASTX analysis were subjected to a second round of BLASTX homology search against the complete GenBank protein database to correct for biased E values and taxonomic misassignments. For some viruses present at a low abundance, we only obtained interspersed reads and no contigs. We then used RT-PCR on cDNA from the virus-positive pool to fill in gaps in the sequence and to generate long genomic fragments. We also employed PCR to confirm that the novel viral RNA sequences originated from authentic viral nucleic acid and not endogenous viral elements. For all authentic viral RNA sequences, amplification products were obtained with RT-PCR, but not with PCR.
Analysis of transovarial transmission. We received 20 I. scapularis adult female ticks that previously fed to repletion on a rabbit. Engorged females were placed in a single individual tube, maintained at 21 Ϯ 2°C and Ͼ90% relative humidity, and allowed to oviposit. When completed, the individual ticks were removed and placed in 1.5-ml Eppendorf tubes, homogenized, and TNAs were extracted from each tick. After cDNA synthesis, each sample was tested for the presence of different viruses by RT-PCR. A portion of egg mass from each tick was placed into a separate 1.5-ml Eppendorf tube. TNAs from the combination of the egg mass and newly hatched larvae were extracted, and cDNA was synthesized and tested by RT-PCR for the presence of the same viruses as the female ticks to provide evidence for vertical transmission.
Phylogenetic analysis. All alignments were performed using MUSCLE alignment software in Geneious 10.1.3 (39). Phylogenetic trees were constructed with MEGA 6.06, and the robustness of each node was determined using 1,000 bootstrap replicates (40). For rhabdoviruses, amino acid sequences of the complete RdRp from 157 viruses were aligned, and the phylogenetic relationship was determined using the maximum likelihood (ML) method employing an LG ϩ G model with nearest-neighbor interchange (NNI). For chuviruses, phylogenetic analysis was performed on a 331-aa conserved RdRp fragment from 38 chuvirus-like amino acid sequences, using the ML method employing an LG ϩ G ϩ I model with NNI. For bunyaviruses, phylogenetic analysis was performed on a 386-aa RdRp fragment consisting of the premotif through motif E from 46 viruses within Bunyavirales. The phylogenetic relationship was determined using the ML method employing an LG ϩ G ϩ F model with NNI. For nodaviruses, a 606-aa fragment of the conserved region of the RdRp from 66 putative nodaviruses was used to generate a maximum likelihood tree with the Whelan and Goldman model with discrete gamma distribution (ϩG) and using the nondefault amino acid frequencies (ϩF) of the model. Accession number(s). Accession numbers newly reported here are shown in bold in Table 3.