Integrated Jingmenvirus Polymerase Gene in Ixodes ricinus Genome

Members of the jingmenviruses group have been found in arthropods and mammals on all continents except Australia and Antarctica. Two viruses of this group were isolated from patients with fever after a tick bite. Using a nested RT-PCR assay targeting a jingmenvirus polymerase gene fragment, we screened ticks collected in seven regions of Russia and found that the abundant jingmenvirus-positive were of Ixodes ricinus species, with the prevalence ranging from 19.8% to 34.3%. In all cases, DNase/RNase treatment suggested that the detected molecule was DNA and subsequent next generation sequencing (NGS) proved that the viral polymerase gene was integrated in the I. ricinus genome. The copy number of the integrated polymerase gene was quantified by qPCR relative to the ITS2 gene and estimated as 1.32 copies per cell. At least three different genetic variants of the integrated polymerase gene were found in the territory of Russia. Phylogenetic analysis of the integrated jingmenvirus polymerase gene showed the highest similarity with the sequence of the correspondent gene obtained in Serbia from I. ricinus.

The viruses of this group have a genome segmented into four segments [1]. Segment 1 encodes a nonstructural protein (NSP1), which has homology to the flavivirus NS5 protein (RNA-dependent RNA polymerase). Segment 3 encodes a nonstructural protein that resembled the NS2b-NS3 complex of flaviviruses. The other two segments encode nonflavivirus-related proteins and appear to be unique for this virus group.
In the first publication describing the Jingmen tick virus, the authors only found evidence of the RNA form of the virus [1]. Similarly, when Maruyama et al. [17] discovered the Mogiana tick virus in Brazil they checked for the presence of viral DNA, but found only viral RNA. Almost all subsequent studies of JMV were performed as a viral RNA search, with no attention to viral DNA as no evidence was indicating its significance. There are two exceptions: the study describing the Takachi virus [4], where authors had confirmed the absence of DNA forms of the virus, and the work on RNA viruses of ixodid ticks in the French Antilles [18]. The authors of this work showed the presence of JMV-related DNA fragments in Rhipicephalus microplus from Guadeloupe and Martinique. All four segments of the Jingmen tick virus were found as DNA in tick samples [18]. However, the authors concluded that the detection of the viral DNA may be the result of residual RNA contamination [18].
Mostly, the presence of DNA sequences of non-retroviral RNA viruses in genomes of eukaryotic cells arise after the viral infection of a host's germline cells and the integration of viral genes or genomes into a host's chromosomes [19,20]. The term endogenous viral elements (EVEs) is usually used for integrated viral sequences in the animal genome. A number of non-retroviral EVE sequences belonging to several virus families were found in various animal genomes [19,20]. The search for EVEs is carried out mainly by bioinformatics/computational methods. Therefore, in regards to arthropods, the largest number of EVEs were found in Aedes, Culex, and Anopheles mosquitoes, since their genomes are better annotated [21][22][23][24][25][26]. Mosquito genomes most often contained EVEs of negative-sense RNA viruses (Rhabdoviridae, unassigned mononega-like viruses) [27]. Complete genomes of the hard ticks (Ixodidae) have only been annotated for a few species, and as a result, there are only a few works that have investigated EVEs in ticks. A number of EVEs were found in the genomes of Ixodes scapularis and I. ricinus; and phylogenetic analysis revealed that most of them are negative-sense RNA virus-derived sequences [27,28].
In this paper, we report that the JMV polymerase gene has been integrated into the I. ricinus genome.

Study Design
This study began from a tick survey in the Moscow region. Using system targeting of the JMV polymerase gene, we detected a high prevalence of RT-PCR-positive ticks. This led us to the hypothesis that the DNA form of the virus was integrated into the tick genome or the genome of a tick's endosymbiont. To test this hypothesis, we performed DNase and RNase treatment and re-tested all samples for viral DNA (without a reverse transcription step) ( Figure 1). Furthermore, we conducted high-throughput sequencing of several JMV-positive tick samples on HiSeq Illumina to illustrate the presence of sequences mapping to both the viral JMV polymerase gene and tick DNA. Subsequently, we also expanded the study by analyzing ixodid ticks collected in different regions of Russia for the presence of JMV DNA and/or RNA. To avoid confusion, we have created a flowchart that demonstrated the experimental strategy used in the study (Figure 1).

The Study Areas
The list of species included I. ricinus, I. persulcatus, Ixodes trianguliceps, D. reticulatus, D. marginatus, D. nuttalli, H. concinna, and Haemaphysalis japonica, which were collected in seven regions of Russia ( Figure 2). More detailed information including the number of collected ticks and georeferenced data are presented in Table S1. Additionally, we used a collection of ticks published earlier [3].
Russia for the presence of JMV DNA and/or RNA. To avoid confusion, we have created a flowchart that demonstrated the experimental strategy used in the study (Figure 1).

The Study Areas
The list of species included I. ricinus, I. persulcatus, Ixodes trianguliceps, D. reticulatus, D. marginatus, D. nuttalli, H. concinna, and Haemaphysalis japonica, which were collected in seven regions of Russia ( Figure 2). More detailed information including the number of collected ticks and georeferenced data are presented in Table S1. Additionally, we used a collection of ticks published earlier [3].

The Study Areas
The list of species included I. ricinus, I. persulcatus, Ixodes trianguliceps, D. reticulatus, D. marginatus, D. nuttalli, H. concinna, and Haemaphysalis japonica, which were collected in seven regions of Russia ( Figure 2). More detailed information including the number of collected ticks and georeferenced data are presented in Table S1. Additionally, we used a collection of ticks published earlier [3].

Tick Sampling and Processing
Ticks were collected by flagging on transects and tick developmental stage, sex, and species were identified using morphological characteristics [29,30]. Additionally, species identification was confirmed by sequencing the COI gene fragment using primers described previously [31].
The ticks were handled in two ways: approximately half of the samples from the Moscow region (325 from 676 ticks) were processed individually (single tick per tube). For the remaining half (351 from 676), we pooled the ticks as follows: three adult ticks per pool for D. reticulatus and five adult ticks per pool for I. ricinus. Nymphs were studied separately Viruses 2022, 14,1908 4 of 17 (without pooling). I. ricinus ticks from other regions of Russia were tested in pools, while the other ixodid species were tested separately (without pooling). The ticks were washed in 70% alcohol and twice in 0.15 M NaCl, and then homogenized with TissueLyser LT (Qiagen, Hilden, Germany) in the following volumes of 0.15 M NaCl solution: 500 µL for Dermacentor ticks, 300 µL for adult ticks of Ixodes and Haemaphysalis genera, and 200 µL for I. ricinus nymphs. Pools were homogenized in 500 µL of 0.15 M NaCl.
Total nucleic acids were extracted from 100 µL of tick suspension using RIBO-prep kit (AmpliSens, Moscow, Russia) following the manufacturer's instructions. Extracted nucleic acids were split into two aliquots that were treated with DNase I (Thermo Fisher Scientific, Vilnius, Lithuania) and Rnase A (Thermo Fisher Scientific, Vilnius, Lithuania), respectively. We used EcoRI endonuclease (Thermo Fisher Scientific, Vilnius, Lithuania) to distinguish between single-or double-stranded DNA. The nuclease treatment was carried out according to the manufacturer's instructions. Purified RNA was DNase I treated and used to check for the presence of viral RNA.

PCR Assays
For PCR screening, JMV RNA was reverse transcribed using Reverta-L kit (Central Research Institute of Epidemiology, Moscow, Russia) following the manufacturer's instructions. We developed nested PCR assays for the JMV polymerase gene to enable amplification of a 425-bp fragment in the first PCR round and a 233-bp fragment in the second PCR round. The reaction mix for each round contained a set of primers designed to amplify all related JMV, available in GenBank on April 2020 ( Table 1). The detailed description of PCR conditions and reaction mix compositions is provided in Table S2. An aliquot of PCR product (5 µL) was run on a 3% agarose gel and visualized with UV illumination after staining with 0.5 µg/mL ethidium bromide. Table 1. Oligonucleotide sequences of the primers used in PCR amplification and sequencing.

Name
Nucleotide Sequence ( Table 1). The PCR was conducted with the conditions detailed in Table S2.

Quantitative PCR
To quantify the copy number of the integrated JMV polymerase gene, three qPCR assays were developed (Table 2) targeting the JMV polymerase gene, the genome integration site (a fragment of I. ricinus genome where the JMV polymerase gene was inserted), and the internal transcribed spacer 2 (ITS2) fragment of I. ricinus. Estimated limits of detection for the qPCR assays were approximately 1000 copies per ml. All reactions were performed in duplicates. Quantitative PCR amplification was performed with RotorGene Q instrument (Qiagen, Hilden, Germany). Reaction mix compositions and PCR conditions are detailed in Table S2. Tick-borne encephalitis virus (TBEV) RNA was detected using a commercial multiplex PCR kit (AmpliSens TBEV, B. burgdorferi s.l., A. phagocytophilum, E. chaffeensis/E. muris-FL; Central Research Institute of Epidemiology, Moscow, Russia) according to the manufacturer's instructions. Tick-borne encephalitis virus strain Absettarov were obtained from Chumakov Institute of Poliomyelitis and Viral Encephalitides collection. There is no information on the number of ITS2 copies per genome of I. ricinus. However, such data exist for the closely related I. scapularis. According to NCBI, the genome of I. scapularis contains 442 copies of the ITS2 region (Genome assembly ASM1692078v2), hence one diploid tick cell contains 884 copies of the ITS2 region. We used this number as a divisor to calculate tick cells number in each sample: Further, we used the copy number of the integrated JMV polymerase gene and the estimated number of tick cells to calculate the copy number of the integrated JMV polymerase gene per cell:

Copy number insertion per cell =
Copy number insertion per sample N tick cells 2.6. Sequencing and Cloning 2.6.1. Next Generation Sequencing A total of six JMV-positive ticks collected in the Moscow region were used for the library preparation. Naïve RNA/DNA templates without any nuclease treatment were reverse-transcribed using Reverta-L kit (Central Research Institute of Epidemiology, Moscow, Russia). Then, 20 µL of first strand cDNA product was immediately taken into second strand cDNA synthesis using NEBNext ® Ultra™ II Non-Directional RNA Second Strand Synthesis Module, which is a part of NEBNext ® Ultra II RNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, USA). DNA fragmentation was performed in the Focused-ultrasonicator ME220 (Covaris, Woburn, MA, USA) according to the man- ufacturer's instructions with parameters for obtaining 500-bp fragments. The remaining steps of library preparation were completed using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, USA). Size selection of the final libraries to~400-500 bp library fragment size was performed using SPRIselect (Beckman Coulter, Danvers, MA, USA) solid-phase reversible immobilization beads. The quality and fragment length distribution of the obtained libraries were evaluated with Agilent Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany). Sequencing was performed with HiSeq PE Rapid Cluster Kit v2 and HiSeq Rapid SBS Kit v2 (500 cycles) on HiSeq 2500 (Illumina, San Diego, CA, USA).

Sanger Sequencing
Purified PCR products were sequenced bidirectionally using the BigDye Terminator v1.1 Cycle Sequencing kit (Thermo Fisher Scientific, Austin, TX, USA) on an Applied Biosystems 3500xL Genetic Analyzer (Applied Biosystems, Foster City, CA, USA).

Cloning
We obtained a long PCR product of 2361 bp for two JMV-positive samples using the MS-F-875 (5 -TCGGATAGGCTGGAGACTCA-3 ) and MS-R-3100 (5 -TGACGAGATGCCGGTT CCCGA-3 ) primers. These primers flank the overlapping region of JMV polymerase gene and an Ixodes' genomic sequence. More detailed data on reaction mix composition and PCR conditions are provided in Table S2. PCR product (5 µL) was run on a 1.5% agarose gel and visualized using UV illumination after staining with 0.5 µg/mL ethidium bromide. The fragment was directly ligated and cloned in pGEM-T vector (Promega, Madison, WI, USA). Positive clones were screened by PCR with primers flanking the vector and inserts and verified by Sanger sequencing using three pairs of primers (Table S3). For each sample, at least 10 clones were selected for sequencing.

Data Analysis
The 95% confidence intervals of the prevalence for the unpooled ticks were estimated using exact Clopper-Pearson methods in R in the package 'PropCIs' [33]. Pooled tick prevalence for the fixed pool size and confidence limits were estimated using a Bayesian approach and a Gibbs sampler [34,35]. Alpha and beta values for the prior distributions of true prevalence and test sensitivity and specificity were calculated based on unpooled data. The simulation was run on the EPITOOLS web platform for 25,000 iterations, with 5000 iterations discarded to allow for convergence [35]. For a part of the analyzed ticks, pooled into pools of variable size (1 to 23 ticks per pool), we calculated the prevalence using maximum likelihood estimator with the assumption of 100% test sensitivity and specificity [36]. The calculations were performed on the EPITOOLS web platform (https: //epitools.ausvet.com.au, accessed on 26 August 2022).
The nucleotide sequences obtained by sequencing were aligned, compared, and analyzed using ClustalW and BLAST. Gblocks (version 0.91b, Castresana J., Heidelberg, Germany) was used to remove divergent or ambiguously aligned regions [37]. Mega X package [38] was used for phylogenetic analyses. Phylogenetic trees were constructed using the maximum-likelihood method with the general time-reversible model [39] and the Tamura 3-parameter model [40], with bootstrap analysis based on 1000 replicates.
Raw reads of high-throughput sequencing were filtered with Trimmomatic [41] using the SLIDINGWINDOW:4:20 MINLEN:40 parameters. Genome assemblies were completed by using SPAdes 3.14 [42] with default parameters. Flaviviridae sequences were selected by calling BLASTn [43] on assembled contigs using all of the available Flaviviridae nucleotide sequences available in the nucleotide database.

High Prevalence of Jingmenviruses Positive Ticks in Moscow Region
A total of 676 ticks were collected in the Moscow region, 325 of them were analyzed individually, and 351 ticks were pooled (Table 3). Total nucleic acids extracted from ticks and tick pools were reverse transcribed and screened for the JMV fragment with the nested PCR assay. We found 66 PCR-positive specimens ( Table 3), 49 of which were confirmed by sequencing with the nested PCR primers (see Table 1, 'Multiplex Mix In'), and/or Mos-Seq-F-2380, Mos-Seq-R-2925 (Table 1). The remaining 17 positive specimens were confirmed with qPCR assay for the integrated JMV polymerase gene (Table 2).  (Table 3). All samples of D. reticulatus yielded negative results in the PCR assay. JMV cDNA was found in adult I. ricinus ticks (both females and males) as well as in questing nymphs and larvae. The sequencing of positive samples with the JMTV-In-F/JMTV-In-Rev primers (Table 1) confirmed their JMV group origin and demonstrated high similarity to the nucleotide sequence of Serbia-24-28 (MT822179) obtained from I. ricinus in Serbia [44].

RNA Virus or DNA Form?
To test the hypothesis that we detected DNA fragments derived from JMV and not the viral RNA, we used 10 JMV-positive ticks and the TBEV strain Absettarov as a positive control of the RNA virus. Total nucleic acids were treated with either DNase I or Rnase A and then tested by qPCR or qRT-PCR using the primers IRJPG-NS5-F, IRJPG-NS5-Rev, IRJPG-NS5-Probe (Table 2), and commercial kit for TBEV detection. Expectedly, the TBEVpositive control yielded positive signals in the qRT-PCR for both naïve RNA/DNA and DNase I-treated, but not RNase-treated samples (Table 4). In comparison, for JMV-positive ticks the results were quite different. All 10 tested specimens were negative for viral RNA (after DNAase I treatment), thus revealing that all 10 tested JMV-positive ticks contained the DNA form of the JMV polymerase gene.
Encouraged by the results, we re-analyzed all ticks collected in the Moscow region looking for viral DNA by the nested PCR assay without reverse transcription. DNA fragments of the JMV polymerase gene were detected in all previously JMV-positive specimens. Additionally, the treatment of samples with EcoRI resulted in signal loss in the PCR assay indicating that the detected JMV polymerase gene fragment is double-stranded DNA.

. NGS and Further Confirmation by Sanger Sequencing
We have performed high-throughput sequencing of cDNA from six samples of JMVpositive I. ricinus ticks. De novo genomic assemblies of all these samples contained contigs of the JMV polymerase gene, and one sample yielded a contig in which a large fragment of the coding region of the JMV polymerase gene (2733 bp) was connected downstream to an unknown sequence of 1273 bp (Figure 3, Specimen 264). The BLASTx algorithm comparing translated nucleotide sequences showed similarity of the unknown sequence to an uncharacterized I. scapularis protein LOC115330668 (GenBank accession number XP_042148831, 62% identity) and a hypothetical protein HPB47_003995 from I. persulcatus (GenBank accession number KAG0419594, 58% identity) ( Figure 3). The other 639 bp of the unknown sequence did not map to any sequences in GenBank by either BLASTn or BLASTx algorithms.  Encouraged by the results, we re-analyzed all ticks collected in the Moscow region looking for viral DNA by the nested PCR assay without reverse transcription. DNA fragments of the JMV polymerase gene were detected in all previously JMV-positive specimens. Additionally, the treatment of samples with EcoRI resulted in signal loss in the PCR assay indicating that the detected JMV polymerase gene fragment is double-stranded DNA.

NGS and Further Confirmation by Sanger Sequencing
We have performed high-throughput sequencing of cDNA from six samples of JMVpositive I. ricinus ticks. De novo genomic assemblies of all these samples contained contigs of the JMV polymerase gene, and one sample yielded a contig in which a large fragment of the coding region of the JMV polymerase gene (2733 bp) was connected downstream to an unknown sequence of 1273 bp (Figure 3, Specimen 264). The BLASTx algorithm comparing translated nucleotide sequences showed similarity of the unknown sequence to an uncharacterized I. scapularis protein LOC115330668 (GenBank accession number XP_042148831, 62% identity) and a hypothetical protein HPB47_003995 from I. persulcatus (GenBank accession number KAG0419594, 58% identity) (Figure 3). The other 639 bp of the unknown sequence did not map to any sequences in GenBank by either BLASTn or BLASTx algorithms.  Subsequently, the overlapping of the JMV polymerase gene and the unknown sequence was confirmed by PCR and Sanger sequencing using primers flanking the 3 -end of the JMV polymerase gene and the Ixodes' gene (Table S3, MF2835 and MR3100). Ten positive ticks were sequenced, confirming the overlapping (Table S4). Furthermore, we obtained 2361 bp amplicons that cover a fragment of the JMV polymerase gene (1950 bp) and the unknown sequence on the 3 -end (411 bp) and cloned them in the pGEM-T vector (see primers MS-F-875 and MS-R-3100 in the Methods section). Sequencing of the clones also confirmed the overlapping of the JMV polymerase gene and an Ixodes' gene (Table S4).
Both I. scapularis and I. persulcatus ticks belong to the subgenus Ixodes Latreille, 1795, which is also true for I. ricinus [29,45]. Therefore, presence of a hypothetical protein sequence from these two Ixodes species allows us to assume that the genome of I. ricinus also contains the same or a very similar gene. To verify this, we designed a qPCR assay for the detection of the integration site in the I. ricinus genome ( Table 2, primers Ir-3100-F, Ir-3100-R, and Ir-3100-probe) and used it to screen the I. ricinus collection. As a result, we found the integration site in all tested specimens of I. ricinus, both JMV-positive and JMVnegative. Furthermore, we found the integrated JMV polymerase gene in the uninfected I. ricinus cell line IRE/CTVM19 [46]. Thus, we confirmed our hypothesis that the integration site is a part of the tick genome, and the DNA form of the JMV polymerase gene is inserted into the I. ricinus genome. Interestingly, no other integrated segments of JMV were not found after the analysis of high-throughput sequencing data.

Quantifying the Integrated Jingmenvirus Polymerase Gene in Ixodes Ricinus Genome
We developed a qPCR assay to estimate the copy number of the integrated JMV polymerase gene relative to the I. ricinus gene with a known copy number (Figure 4). Currently, the genome of I. ricinus has not been annotated. At the same time, there is a datum on the number of nuclear ribosomal genes for a closely related I. scapularis-(442 copies per genome). Therefore, we decided to use the ITS2 fragment as a reference gene in this assay. A total of 18 I. ricinus ticks were individually tested by qPCR (as naïve DNA without reverse transcription) including eight nymphs, five females, and five males. Two samples with the highest values of the copy number were excluded from consideration as outliers exceeding the interquartile range by 1.5-fold (Table S5). As a result, the median copy number of the integrated JMV polymerase gene was calculated to be 1.32 copies per cell (Table S5). to the fragment of the Ixodes scapularis genome (GenBank accession number XM_042293990). The light blue box represents the fragment with an unknown sequence. The horizontal axis shows the length of the nucleotide sequences.
Subsequently, the overlapping of the JMV polymerase gene and the unknown sequence was confirmed by PCR and Sanger sequencing using primers flanking the 3′-end of the JMV polymerase gene and the Ixodes' gene (Table S3, MF2835 and MR3100). Ten positive ticks were sequenced, confirming the overlapping (Table S4). Furthermore, we obtained 2361 bp amplicons that cover a fragment of the JMV polymerase gene (1950 bp) and the unknown sequence on the 3′-end (411 bp) and cloned them in the pGEM-T vector (see primers MS-F-875 and MS-R-3100 in the Methods section). Sequencing of the clones also confirmed the overlapping of the JMV polymerase gene and an Ixodes' gene (Table S4).
Both I. scapularis and I. persulcatus ticks belong to the subgenus Ixodes Latreille, 1795, which is also true for I. ricinus [29,45]. Therefore, presence of a hypothetical protein sequence from these two Ixodes species allows us to assume that the genome of I. ricinus also contains the same or a very similar gene. To verify this, we designed a qPCR assay for the detection of the integration site in the I. ricinus genome ( Table 2, primers Ir-3100-F, Ir-3100-R, and Ir-3100-probe) and used it to screen the I. ricinus collection. As a result, we found the integration site in all tested specimens of I. ricinus, both JMV-positive and JMV-negative. Furthermore, we found the integrated JMV polymerase gene in the uninfected I. ricinus cell line IRE/CTVM19 [46]. Thus, we confirmed our hypothesis that the integration site is a part of the tick genome, and the DNA form of the JMV polymerase gene is inserted into the I. ricinus genome. Interestingly, no other integrated segments of JMV were not found after the analysis of high-throughput sequencing data.

Quantifying the Integrated Jingmenvirus Polymerase Gene in Ixodes Ricinus Genome
We developed a qPCR assay to estimate the copy number of the integrated JMV polymerase gene relative to the I. ricinus gene with a known copy number (Figure 4). Currently, the genome of I. ricinus has not been annotated. At the same time, there is a datum on the number of nuclear ribosomal genes for a closely related I. scapularis-(442 copies per genome). Therefore, we decided to use the ITS2 fragment as a reference gene in this assay. A total of 18 I. ricinus ticks were individually tested by qPCR (as naïve DNA without reverse transcription) including eight nymphs, five females, and five males. Two samples with the highest values of the copy number were excluded from consideration as outliers exceeding the interquartile range by 1.5-fold (Table S5). As a result, the median copy number of the integrated JMV polymerase gene was calculated to be 1.32 copies per cell (Table S5).

Screening for Jingmenviruses and the Integrated Jingmenvirus Polymerase Gene in Ixodes Ricinus from Different Regions of Russia
A total of 140 pools of I. ricinus (451 ticks) from six regions were tested using two qPCR assays: (1) qPCR specific for the integrated JMV polymerase gene (Table 2); and (2) nested PCR using cDNA as a template, which was obtained after DNase I treatment of naïve RNA/DNA ( Table 5). The collection sites of these ticks are widely distributed across the Russian part of the I. ricinus geographic range, covering the northern, eastern, and southern parts of the global range of the species (Figure 2). DNA fragments of the JMV polymerase gene were detected in samples from each region. Due to the pool size variability of the samples (from 1 to 10 for adult ticks and from 5 to 23 for nymphs), we estimated the prevalence of ticks with the integrated JMV polymerase gene using the maximum likelihood method, assuming 100% test sensitivity and specificity. The observed prevalence varied from 10.1% to 36.7% (Table 5). Two screened samples were PCR-positive for the Alongshan virus and sequencing revealed the presence in these samples of both Alongshan virus RNA and DNA of the integrated jingmenvirus polymerase gene. However, it is important to note, that both Alongshan-positive samples were pools of five ticks per sample. Therefore, we cannot conclude that the Alongshan virus was found in the same tick that also had the integrated JMV polymerase gene.

Screening for Jingmenviruses and the Integrated Jingmenvirus Polymerase Gene in Other Ixodid Tick Species
A total of 412 ticks representing eight species collected in the different regions of Russia were analyzed by qPCR specific for the integrated JMV polymerase gene ( Table 2) and nested PCR on cDNA obtained after the DNase I treatment of naïve RNA/DNA (Table S6). This screening revealed five specimens positive for the Alongshan virus and two specimens positive for the Yanggou tick virus. These findings were described in our previous paper [3]. Here, we analyzed the samples for the presence of the DNA fragments of the integrated JMV polymerase gene, which were not detected in any of the ticks (Table S6).

Genetic Variability of the Jingmenvirus Polymerase Gene Fragment in the Studied Ticks
Phylogenetic analysis of the JMV polymerase gene based on the 2730 bp fragments showed that the discovered integrated DNA (I. ricinus Russia Moscow-264 isolate) is the most similar (p-distance = 0.012) to the sequence of the JMV polymerase gene obtained in Serbia from I. ricinus (GenBank accession number MT822179, Figure 5). Additionally, GenBank contains a similar sequence (GenBank accession number JXMZ02142294, Figure 5), which was annotated in the genomic DNA of the I. ricinus strain (Charles River) obtained by whole genome shotgun sequencing (p-distance = 0056). As a result, the integrated JMV polymerase sequence together with the sequences obtained from Serbian ticks form a separate cluster in the JMV group and shows a close relationship to the JMV polymerase gene of the Alongshan virus ( Figure 5). which was annotated in the genomic DNA of the I. ricinus strain (Charles River) obtained by whole genome shotgun sequencing (p-distance = 0056). As a result, the integrated JMV polymerase sequence together with the sequences obtained from Serbian ticks form a separate cluster in the JMV group and shows a close relationship to the JMV polymerase gene of the Alongshan virus ( Figure 5). The genetic variability of the integrated JMV polymerase gene in ticks is also noteworthy. Sequencing of the 3′ end region of the gene (509 bp fragment) from 35 ticks from the Moscow region and 24 ticks from other regions of Russia revealed the presence of at least three different genetic variants ( Figure 6). The first genovariant, represented in our sample by one specimen collected in the Moscow region (id 100), appears to have been the earliest deviation into a separate branch and has a 23-29 nt difference to the other genovariants (p-distances-0.045-0.057). The second genovariant was present in seven samples and is removed from the others by a distance of 17-23 nt (p-distances-0.033-0.045). This genovariant was isolated from ticks collected in the Moscow, Kaliningrad, and Voronezh regions (Figure 2). Two of the samples simultaneously contained the first and the second genovariants. One of these samples was a pool of five ticks, and the second sample was an adult female of I. ricinus. The genetic variability of the integrated JMV polymerase gene in ticks is also noteworthy. Sequencing of the 3 end region of the gene (509 bp fragment) from 35 ticks from the Moscow region and 24 ticks from other regions of Russia revealed the presence of at least three different genetic variants ( Figure 6). The first genovariant, represented in our sample by one specimen collected in the Moscow region (id 100), appears to have been the earliest deviation into a separate branch and has a 23-29 nt difference to the other genovariants (p-distances-0.045-0.057). The second genovariant was present in seven samples and is removed from the others by a distance of 17-23 nt (p-distances-0.033-0.045). This genovariant was isolated from ticks collected in the Moscow, Kaliningrad, and Voronezh regions (Figure 2). Two of the samples simultaneously contained the first and the second genovariants. One of these samples was a pool of five ticks, and the second sample was an adult female of I. ricinus.
The third genovariant in the clade was the most abundant and included 51 sequences obtained from the Moscow region, the Republic of Karelia, the Republic of Tatarstan, Voronezh, Kaliningrad, Belgorod, and Ulyanovsk regions of Russia. Moreover, a sequence from GenBank obtained from I. ricinus in Serbia (accession number MT822179) also belongs to this cluster ( Figure 6). Notably, the ticks from the Moscow region in all three clusters were collected on the same path and the distance between the collection sites did not exceed 1 km. Viruses 2022, 14, x FOR PEER REVIEW 12 of 17 The third genovariant in the clade was the most abundant and included 51 sequences obtained from the Moscow region, the Republic of Karelia, the Republic of Tatarstan, Voronezh, Kaliningrad, Belgorod, and Ulyanovsk regions of Russia. Moreover, a sequence from GenBank obtained from I. ricinus in Serbia (accession number MT822179) also belongs to this cluster ( Figure 6). Notably, the ticks from the Moscow region in all three clusters were collected on the same path and the distance between the collection sites did not exceed 1 km.

Discussion
In our study, we clearly demonstrate the integration of the JMV polymerase gene fragment into the genome of I. ricinus from seven regions of Russia. A number of other ixodid species inhabiting Russia were free from the integrated JMV polymerase gene, but additional studies of different tick species in an expanded geographic range are needed to validate a more general assumption. Based on the phylogenetic study of the hard ticks [47], the most promising strategy could be to look at the closely related species of the genus Ixodes: I. acuminatus, I. redikorzevi, I. persulcatus, and I. pavlovskyi.
The possible integration of viral sequences into the host genome was previously investigated by other authors working with viruses of the JMV group. For example, when the Jingmen tick virus was first discovered, the presence of DNA forms of this virus in the ticks' samples was excluded based on the absence of an amplicon in PCR (without reverse transcription) of total DNA [1]. In the work reporting the first detection of the Mogiana tick virus-related nucleic acid sequences in R. microplus ticks in Brazil, the authors checked whether these sequences are of RNA or DNA origin [17]. A comparison of the amplification efficiency starting from either DNA or cDNA isolated from ticks showed that the virus-related sequences were in the form of RNA, and the DNA form was absent [17].

Discussion
In our study, we clearly demonstrate the integration of the JMV polymerase gene fragment into the genome of I. ricinus from seven regions of Russia. A number of other ixodid species inhabiting Russia were free from the integrated JMV polymerase gene, but additional studies of different tick species in an expanded geographic range are needed to validate a more general assumption. Based on the phylogenetic study of the hard ticks [47], the most promising strategy could be to look at the closely related species of the genus Ixodes: I. acuminatus, I. redikorzevi, I. persulcatus, and I. pavlovskyi.
The possible integration of viral sequences into the host genome was previously investigated by other authors working with viruses of the JMV group. For example, when the Jingmen tick virus was first discovered, the presence of DNA forms of this virus in the ticks' samples was excluded based on the absence of an amplicon in PCR (without reverse transcription) of total DNA [1]. In the work reporting the first detection of the Mogiana tick virus-related nucleic acid sequences in R. microplus ticks in Brazil, the authors checked whether these sequences are of RNA or DNA origin [17]. A comparison of the amplification efficiency starting from either DNA or cDNA isolated from ticks showed that the virus-related sequences were in the form of RNA, and the DNA form was absent [17]. Discovering the Guaico Culex virus, Ladner et al. also interrogated the nature of the genomic material using nuclease treatment (DNase I and RNase A) and confirmed that all five genome segments were single-stranded, positive-sense RNA, and DNA forms were absent [48]. In the work describing the Takachi virus (JMV group), authors checked the presence of DNA forms by PCR without reverse transcription [4]. In subsequent works, only NGS sequencing of the virus was carried out and the presence/absence of DNA forms was not accounted for [2,7,9,10]. In these works, library preparation was done after DNase treatment, therefore, even if DNA forms were present in the samples, they would have been lysed with nucleases. The presence of DNA forms of JMVs was detected in R. microplus from the Antilles, and both viral RNA and DNA of all four segments of the viral genome were found in 1-8% of JMV-positive samples [18]. Because of this key difference, we cannot argue that the authors also found insertions into the tick genome. The authors themselves questioned the validity of their findings of the JMV DNA form and suggested contamination as a possible explanation.
On the other hand, in the study of ter Horst et al. [28], genomes of arthropods, including I. ricinus and I. scapularis, were screened for the presence of EVEs. In the genome of I. scapularis and I. ricinus, 387 and 168 viral sequences were found, respectively [28]. One of the sequences found in the I. ricinus genome (GenBank accession number: JXMZ02142294) corresponds to the polymerase gene of the Jingmen tick virus. Phylogenetic analysis showed that this sequence is closely related to the Serbian sequence MT822179 [44] and to the JMV polymerase gene inserted in the I. ricinus genome detected in our study ( Figure 5). Importantly, the preparation of Serbian I. ricinus samples did not involve DNase treatment. Based on the genetic similarity of the Serbian sequences to the one we describe, the fact that they were both found in I. ricinus, and the presence of only the JMV polymerase gene and that no other viral genome segments were found, we can conclude that the sequences (GenBank accession numbers MT822179, MT822180) deposited by Stanojević et al. [44] come not to the Jingmen tick virus directly, as the authors pointed out, but are probably also an integrated DNA of the JMV polymerase gene.
Our phylogenetic analysis showed that the integrated fragment of the JMV polymerase gene is variable and forms at least three genovariants sufficiently removed from each other (range of p-distances 0.033-0.057). There are two possible explanations for the genetic heterogeneity of the insertion: (1) the integration event occurred once and the observed variability have been accumulated by the DNA form; (2) observed variability had existed in the ancestral virus population before the integration, and integration events occurred several times with the different genovariants. It is known that the mutation rate for virus-derived sequences integrated into a eukaryotic genome is lower than for the RNA viruses [49], so the observed genetic divergence had most likely been accumulated before integration into the I. ricinus genome. This suggests that insertion into the I. ricinus genome has occurred several times, and these insertion events happened in germline cells. One of the insertion's genovariants was widespread in the I. ricinus area. Therefore, we can conclude that the insertion into the I. ricinus genome of this variant occurred long enough to spread within the species range. Our study showed that in the Moscow region the prevalence of ticks with this insertion was approximately 27-35%. The fact that not all studied I. ricinus had this insertion suggests that the insertion event has occurred after I. ricinus split into a separate species. Thus, the insertion of the JMV polymerase gene could potentially become a convenient marker for studying the genetic structure and phylogeography of I. ricinus populations.
While screening our I. ricinus collection, we found two specimens that were PCRpositive for both the JMV-integrated polymerase gene and the Alongshan virus. Both Alongshan-positive samples were pools of five ticks, preventing us from making a conclusion that the virus could successfully replicate in ticks with the integrated JMV polymerase gene. However, we also found the integrated gene in uninfected cell culture IRE/CTVM19. This cell culture has previously been used to successfully isolate the Alongshan and Yanggou tick viruses [3,8]. This suggests that the integrated JMV polymerase gene does not significantly affect the ability of the Alongshan virus and the Yanggou virus to replicate in tick cell culture. It is possible that virus-derived DNA of ancestral JMV could be involved in the antiviral immune response of host cells and lead to increased tolerance (low damage for host cell and high virus load). It was shown that during RNA-virus infection, Aedes mosquitos generate virus-derived DNA that is essential for viral tolerance [50]. A similar mechanism of tolerance to arboviruses could potentially exist in ticks even before the integration event.
Viruses 2022, 14, 1908 14 of 17 We were unable to pinpoint the location of the insertion in the I. ricinus genome because the tick genome has not yet been annotated. The insertion was found in ticks of both sexes, which allows us to conclude that the JMV polymerase gene is not integrated into the sex Y chromosome. The assessment of the copy number of the insertion showed that one tick cell contains, on average, 1.5 copies of the insertion. Therefore, we hypothesize that only one copy of the JMV polymerase gene was inserted in the haploid genome of I. ricinus. Respectively, in a diploid genome, this number would be increased to two copies per cell and lower values do not contradict the single insertion hypothesis. Indeed, the insertion could be present in both the biallelic variant and in the single allelic variant (hemizygote) of the tick's genome. A monoallelic variant is possible because not all individuals in the population have this insertion (particularly in the Moscow region population where the prevalence of ticks with insertion was around 27-35%). The offspring of two ticks, one with the insertion and one without, will inherit a monoallelic variant of the insertion. With that, the copy number estimation in our work has a certain limitation in that it is based on the data for the I. scapularis, not I. ricinus, and therefore this issue requires deeper investigation. Regarding the copy number of non-retroviral integrated RNA virus sequences, in Aedes mosquitoes, up to 71% of the virus-derived sequences had at least one additional copy in the mosquito's genome [27]. In contrast, in the genome of I. scapularis, only 18% of all virus-derived sequences had an additional copy [27]. Our results about the copy number of the integrated fragment are thus consistent with what has been observed for I. scapularis.
Russo et al. [27] showed that among the non-retroviral integrated RNA virus sequences in I. scapularis, 81% is accounted for the genes encoding RNA-dependent RNA polymerase (RdRp gene), and the remaining 19% are from structural genes (nucleoprotein, glycoprotein, and matrix regions). In the case of I. scapularis, when the high frequency of RdRp integration was observed, the suggested explanation was an unusual template preference for the endogenous reverse transcriptases and integrases, which could facilitate the integration of virus-derived sequences in I. scapularis [27]. We cannot directly extrapolate the results of the study by Russo et al. [27] to our findings. However, we can speculate that for the endogenous reverse transcriptases there could also exist a template preference for the JMV polymerase gene of segmented positive-sense RNA viruses.
The presence of the reverse transcriptase enzyme is required to make a DNA copy from the viral RNA template. Presumably, non-retroviral integrated RNA virus sequences use the reverse transcriptase encoded by the retroelements abundant in eukaryotic genomes [19]. Indeed, some studies of EVE in insect genomes showed that retrotransposons are involved in reverse transcription and the integration of viral genomic fragments [50][51][52][53][54]. Particularly, for the Aedes mosquito genome, there is clear evidence for the association of RNA-virusderived sequences and long-terminal repeat (LTR) retrotransposons [52,55,56]. For ixodid ticks, Russo et al. [27] suggested that the reverse transcriptase activity of non-long-terminal repeat (non-LTR) retrotransposons, particularly I, L1, and R1 elements, might be responsible for the appearance of RNA virus-derived sequences.

Conclusions
In conclusion, we would like to highlight the following main results of our work. Populations of I. ricinus ticks contain a JMV polymerase gene, which was integrated into the genome after separation of I. ricinus as a species. Insertion has occurred several times and long ago enough for ticks with the insertion to spread widely within the species range.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/v14091908/s1, Table S1: Georeferenced data of collected ticks; Table S2: PCR mixture and PCR conditions for self-developed PCR assays; Table S3: Oligonucleotide primers used for sequencing of cloned amplicon, which overlapped the flavi-NS5-like protein gene and the fragment of host genome; Table S4: Data on Ixodes ricinus ticks that were PCR-positive on the Integrated Jingmenvirus polymerase gene; Table S5: Quantitative PCR assessment of the copy number of the inserted flavi-NS5-like protein gene in I. ricinus ticks; Table S6: The results of nested-PCR screening on JLV of other ixodid ticks (not I. ricinus).