iSeq 100 for metagenomic pathogen screening in ticks

Background Ticks are blood-sucking ectoparasites that play a pivotal role in the transmission of various pathogens to humans and animals. In Korea, Haemaphysalis longicornis is the predominant tick species and is recognized as the vector of pathogens causing various diseases such as babesiosis, borreliosis, rickettsiosis, and severe fever with thrombocytopenia syndrome. Methods In this study, the targeted high-throughput sequencing of the 16S rRNA V4 region was performed using the state-of-the-art sequencing instrument, iSeq 100, to screen bacterial pathogens in H. longicornis, and the findings were compared with those using conventional PCR with specific primers. Microbiome analyses were performed with EzBioCloud, a commercially available ChunLab bioinformatics cloud platform. ANOVA-Like Differential Expression tool (ALDEx2) was used for differential abundance analysis. Results Rickettsia spp. were detected in 16 out of 37 samples using iSeq 100, and this was confirmed using a PCR assay. In the phylogenetic analysis using gltA and ompA sequences of the detected Rickettsia, the highest sequence similarity was found with ‘Candidatus Rickettsia jingxinensis’ isolate Xian-Hl-79, ‘Ca. R. jingxinensis’ isolate F18, and ‘Ca. R. longicornii‘ isolate ROK-HL727. In the microbiome study, Coxiella AB001519, a known tick symbiont, was detected in all 37 tick samples. Actinomycetospora chiangmaiensis was more abundant in Rickettsia-positive samples than in Rickettsia-negative samples. Conclusions In this study, iSeq 100 was used to investigate the microbiome of H. longicornis, and the potentially pathogenic Rickettsia strain was detected in 16 out of 37 ticks. We believe that this approach will aid in large-scale pathogen screening of arthropods to be used in vector-borne disease control programs. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-021-04852-w.


Background
Ticks are blood sucking ectoparasites that play a pivotal role in the transmission of a variety of pathogens to humans and animals [1,2]. Ticks harbor numerous bacterial, protozoal, and viral pathogens and inject these into the host's body when sucking blood, causing serious infectious diseases [3,4]. Tickborne diseases such as anaplasmosis, ehrlichiosis, borreliosis, babesiosis, and rickettsiosis represent emerging threats to public and animal health worldwide [5,6]. Each tick species has preferred environmental conditions, which determine the geographic distribution of the tick and, consequently, the risk areas for tick-borne diseases [7].
For pathogen detection, PCR and culture-based tests are commonly used. PCR assays are comparatively expensive and time-consuming and have poor operational characteristics, limiting their utility in identifying the cause of infectivity [20]. Culture-based approaches are limited as many infectious agents require special culture conditions [20]. Moreover, in the case of pathogens transmitted from arthropods, many are unculturable.
The main advantage of next-generation sequencing (NGS) is that a wide array of known and unknown pathogens can be identi ed simultaneously, without the need for designing individual speci c bacterial primers [21]. The NGS approach is faster and more cost-effective than conventional methods (such as pathogen-speci c PCR or cultivation) in cases where the sample number or the number of potential pathogens to be detected is high.
NGS can also provide quantitative data on the concentration of organisms in a sample by counting the sequence reads [21]. iSeq100, a recently released instrument, is capable of processing a targeted NGS library with accuracy comparable to that of MiSeq, the most widely used instrument in microbiome studies [22,23]. iSeq100 also features shortened sequencing work ow preparation time and a total sequencing run time of under 24 hours, as it allows the sequencing of a lower read count than the MiSeq [22]. The reduced capital and maintenance costs are also advantages of iSeq100 over MiSeq [22].
In this study, we used iSeq100 to detect bacterial pathogens in H. longicornis and compared the pathogen detection rates to those of pathogen-speci c PCR test. In addition, we identi ed the strain of the detected potential pathogens using PCR with speci c primers.

Sample collection
Ticks were collected from the vegetation by agging in Wonju (Gangwon-do Province; 37.389545, 127.801770) in July 2015. Species identi cation of collected ticks was performed by examination under a dissecting microscope according to the method of Yamaguti et al. [24]. DNA was extracted from each tick using a NucleoSpin DNA Insect kit (Macherey-Nagel, Düren, Germany) according to the manufacturer's instructions, and stored at -20 ˚C until use.

Illumina sequencing and bioinformatics
Total DNA was extracted using the NucleoSpin DNA Insect Kit (Macherey-Nagel, Düren, Germany) following the instructions of the manufacturer. The 16S rDNA V4 region was ampli ed by PCR using primers of primers of 515F (5 -TCGTCGGCAGCGTCAGATG TGTATAAGAGACAGGTGCCAGCMGCCGCGGTAA-3 ) and 806R (5 -GTCTCGTGGG CTCGGAGATGTGTATAAGAGACAGGGACTACHVGGGTWTCTAAT-3 ). A limited-cycle (eight cycles) ampli cation step was performed to add multiplexing indices and Illumina sequencing adapters. Mixed amplicons were pooled, and the sequencing was performed with the Illumina iSeq100 sequencing system according to the manufacturer's instructions, utilizing an Illumina iSeqTM 100 i1 Reagent v2 kit (San Diego, California, USA).

Bioinformatics and statistics
Processing raw reads started with a quality check and the ltering of low quality (< Q25) reads by Trimmomatic ver. 0.321 [25]. After a quality control pass, sequence data were merged using the fastq_mergepairs command of VSEARCH version 2.13.42 with default parameters. Primers were then trimmed with the alignment algorithm of Myers & Miller [26] at a similarity cut-off of 0.8. Non-speci c amplicons that do not encode 16S rRNA were detected by nhmmer [27] in the HMMER software package ver. 3.2.1 with hmm pro les. Unique reads were extracted, and redundant reads were clustered with the unique reads by derep_fulllength command of VSEARCH2. The EzBioCloud 16S rRNA database [28] was used for taxonomic assignment using the usearch_global command of VSEARCH2, followed by more precise pairwise alignment [26]. Chimeric reads were ltered on reads with < 97% similarity by referencebased chimeric detection using the UCHIME algorithm [29] and the non-chimeric 16S rRNA database from EzBioCloud. After chimeric ltering, reads that were not identi ed to the species level (with < 97% similarity) in the EzBioCloud database were compiled, and the cluster_fast command (2) was used to perform de novo clustering to generate additional operational taxonomic units (OTUs). Finally, OTUs with single reads (singletons) were omitted from further analysis. All subsequent analyses were performed with EzBioCloud, a commercially available ChunLab bioinformatics cloud platform for microbiome research (https://www.ezbiocloud.net/). The reads were normalized to 2,500 to perform the analyses. We computed the Shannon index [30] and performed principal coordinate analysis (PCoA) [31]. The Wilcoxon rank-sum test was used to test the difference in the number of OTUs and the Shannon index between two groups. Linear discriminant analysis effect size (LEfSe) analysis was used to identify signi cantly different taxa among the groups [32].
A BLAST search was used to compare the obtained sequence of the gltA gene and ompA gene to those available in GenBank (USA). The obtained sequences were compared for similarity to sequences deposited in GenBank using BLAST. Gene sequences, excluding the primer regions, were aligned using the multisequence alignment program in Clustal Omega (Cambridgeshire, UK) [35], and phylogenetic analysis was performed using MEGA-X software [36].

Results
The microbiome composition of ticks analyzed with iSeq100 In the high-throughput sequencing of the 16S rRNA gene of 39 H. longicornis tick samples using iSeq100, the average reads assigned to bacteria were 7,285 assigned to 72 species (OTUs). At the species level, all samples were dominated by Coxiella AB001519_s (5.48-89.51% of the total community, average: 42.10%) (Fig.1). Williamsia maris group (0.04-41.83% of the total community, average: 13.83%), which was the second most abundant group, was also detected in all samples. The Rickettsia rickettsii group was detected in 18 out of 39 tick samples.

Detection of rickettsiae using PCR
For the same set of samples, only 14 samples were detected as positive for rickettsiae using nested PCR of the 17 kDa surface protein gene, but four samples were negative. Positive samples in PCR analysis were also positive in NGS using iSeq (Fig. 1, Table 1).

Alpha and beta diversity of Rickettsia-positive and Rickettsia-negative groups
To investigate whether Rickettsia infection can affect the composition of the tick microbiome or if the tick microbiome can affect the Rickettsia infection, we compared the structure and relative abundances of bacteria in the Rickettsia-positive group with those in the Rickettsia-negative group. For the Rickettsia-positive and Rickettsia-negative samples, the average reads assigned to bacteria were 7,012 reads to 73 species (OTUs) and 7,604 reads assigned to 71 species, respectively. The number of OTUs (bacterial species richness) were not signi cantly different between the Rickettsia-positive (n = 18) and Rickettsia-negative (Fig. 3a) groups. However, the Shannon index (bacterial diversity) was signi cantly different between the two groups (Fig. 3b). In the PCoA, the samples were organized well according to Rickettsia positivity (Fig. 3c). To assess the beta diversity, we performed permutational multivariate analysis of variance (PERMANOVA). Rickettsia positivity was a signi cant factor in determining the tick microbiome composition (P = 0.000).
Composition in the Rickettsia-positive and Rickettsia-negative groups Next, we analyzed the composition of bacteria between each sample in the two groups to determine whether the relative abundances of bacteria were different in the Rickettsia-positive and Rickettsia-negative groups. Coxiella AB001519_s accounted for 44.82% and 35.96% of the total reads in the Rickettsia-positive and Rickettsia-negative groups, respectively (Fig. 3). There was no difference in the relative abundances of other species that accounted for more than 1%, on average, of the total reads between the Rickettsia-positive and Rickettsia-negative groups.

LEfSe analysis
To identify signi cant differences in bacterial abundance between the groups, we performed LEfSe analysis (Fig. 4). The bacterial species with the highest linear discriminant analysis (LDA) scores in the Rickettsia-positive group were the Actinomycetospora chiangmaiensis group (LDA = 3.68) and the Rickettsia rickettsii group (LDA = 4.57). In the Rickettsia-negative group, Coxiella AB001519 (LDA=4.81) was the species with the highest LDA score.

Discussion
Ticks can transmit bacterial, parasitic, and viral pathogens (often zoonotically) and often harbor more than one agent simultaneously [37]. Thus, obtaining broader information about pathogens in ticks is important from the perspective of proper diagnosis and treatment of these diseases [38]. In this study, we demonstrated that the recently released NGS instrument, iSeq100, is useful for screening of bacteria in ticks. NGS approaches have the ability to identify a wide range of known or unknown pathogens or discover new organisms from a single test [21] without the need to design speci c primers for each pathogen. This method makes it possible to identify pathogens immediately, not only in ticks, but also in arthropods that serve as vectors and reservoirs for pathogens, such as mosquitoes, tsetse ies, and sand ies.
Considering this advantage, in this study, we screened the tick-borne bacteria using iSeq100 and found that 18 of the 31 ticks harbored pathogens of the R. rickettsii species. For an accurate taxonomic classi cation of the species of the detected Rickettsia spp., we used conventional PCR to amplify gltA and ompA sequences of the Rickettsia spp. and compared them with sequences deposited in GenBank using BLAST. In the phylogenetic analysis performed using MEGA-X software, the sequence similarity to sequences of Ca. R. jingxinensis isolate Xian-Hl-79, Ca. R. jingxinensis isolate F18, and Ca. R. longicornii isolate ROK-HL727 was 100 % for gltA, indicating a close relationship between rickettsial isolates from H. longicornis from Korea with those from other East Asian countries. The close clustering of the Chinese and Korean strains of Rickettsia spp. may indicate a close epidemiological link between these strains.
In this study, we used iSeq100 and PCR with speci c primers for screening tick-borne bacteria, and more Rickettsia-positive samples were detected with the iSeq100 method. In iSeq100 analysis, 18 samples showed positive results, while in the pathogen-speci c PCR method, 14 samples were positive. However, sequencing using conventional PCR remains essential for identifying the speci c strain of the pathogen species.
Ixodid ticks (e.g. H. longicornis, H. ava, Ixodes persulcatus, and I. nipponensis) in Asia have the potential to be primary vectors/repositories of rickettsiae of medical and veterinary importance [16]. In 2006, the rst case in Korea of R. japonica was isolated from a spotted fever patient [14]. Recent studies show a high prevalence of the emerging pathogen, Rickettsia raoultii, in canine ticks [43].
In our study, the sequences of gltA and ompA are identical or highly homologous to those of Ca. R. jingxinensis isolate Xian-Hl-79, Ca. R. jingxinensis isolate F18, and Ca. R. longicornii isolate ROK-HL727.
Ca. R. jingxinensis, a novel Rickettsia species in Rhipicephalus microplus and H. longicomis ticks, was rst discovered in China (Shenyang and Wuhan) [44][45][46] and subsequently reported in Korea (Chungnam, Jeonbuk, and Gwangju) [47,48]. Many associated Ca. R. jingxinensis sequences have been deposited in GenBank. Of these, a gltA sequence (KU853023) was recovered from a patient, suggesting its potential pathogenicity to humans [49]. The pathogenicity of Ca. R. longicornii is yet to be determined. However, 99.6% identity was detected with the Ca. R. longicornii ompA sequence of an isolate from rodent spleen tissue obtained in Korea, and 99.9% identity was detected with the Ca. R. longicornii gltA sequence of an isolate from a human blood sample obtained in China [19,49]. These results suggest that Ca. R. longicornii has the potential to infect mammalian hosts, including humans [50].
Coxiella AB001519, known as a Coxiella-like symbiont of H. longicornis, was detected in all 39 samples [57]. A bacterium belonging to group Coxiella was reported as the primary endosymbiont of Amblyomma americanum ticks and was found to improve the reproductive health of the ticks [58-60]. Coxiella AB001519 was rst identi ed in Japan in a phylogenetic association with the Coxiella-like endosymbiont of H. longicornis, and the Coxiella-like endosymbiont presenting more than 99% homology with Coxiella AB001519 in Thailand, China, and Korea was found in H. longicornis [61][62][63].
In the present study, the abundance of Actinomycetospora chiangmaiensis found to be higher in Rickettsia spp.

Conclusions
In conclusion, in this study, we demonstrated that iSeq100 can rapidly and economically screen pathogens in ticks. This approach enables large-scale pathogen screening of arthropods for vector-borne disease control programs.  Table   Table 1 The microbiome composition of each tick at the species level (n = 39). Species accounting for more than 1% of total reads are shown. Tick samples 1 to 18 (red) were found to harbor the Rickettsia rickettsii group (relative abundance > 0.1%). Samples positive for Rickettsia spp. in the nested PCR are marked with a clover (♣) Figure 2 Phylogenetic trees based on partial sequences of gltA (a) and ompA (b) genes. Sequences of the Rickettsia spp. detected in the present study were aligned with those retrieved from the GenBank database. Phylogenetic analysis was performed using MEGA-X software. The tree is drawn to scale, with branch lengths (below the branches) in the same units as those of the taxonomic distances used to infer the phylogenetic tree. The average microbiome composition of Rickettsia-positive ticks (n = 18) and Rickettsia-negative ticks (n = 21) at the species level. Species accounting for than 1% of total reads are shown.