Comparative population genetic structure of two Ixodidae ticks (Ixodes ovatus and Haemaphysalis flava) in Niigata Prefecture, Japan

doi:10.21203/rs.2.18669/v1

PDF

Research

Comparative population genetic structure of two Ixodidae ticks (Ixodes ovatus and Haemaphysalis flava) in Niigata Prefecture, Japan

https://doi.org/10.21203/rs.2.18669/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 11 Dec, 2019

You are reading this latest preprint version

Background

Ixodid tick species such as Ixodes ovatus and Haemaphysalis flava function as important vectors of tick-borne diseases in Japan. The study of the genetic patterns of tick populations can reveal information regarding the spread of tick-borne disease. We hypothesized that I. ovatus and H. flava have different population genetic structure because of their host mobility in different tick life stages despite sharing of hosts.

Methods

Samples (n = 1 to 77) were collected in 29 (I. ovatus) and 17 (H. flava) sampling locations across Niigata. In this study, we used genetic structure at two mitochondrial loci (cox1, 16S rRNA gene) to infer gene flow patterns of I. ovatus and H. flava from Niigata Prefecture, Japan.

Results

For I. ovatus, pairwise F_STand analysis of molecular variance (AMOVA) analyses of cox1 sequences indicated significant among-population differentiation. This was in contrast to H. flava, for which there were only two cases of significant pairwise differentiation and no overall structure. A Mantel test revealed isolation by distance and there was positive spatial autocorrelation of haplotypes in I. ovatus cox1 and 16S sequences, but non-significant results were observed in H. flava in both markers. We found three genetic groups (China 1, China 2 and Japan) in the cox1 I. ovatus tree. Newly sampled I. ovatus grouped together with a published I. ovatus sequence from northern Japan and were distinct from two other I. ovatus groups that were reported from southern China.

Conclusions

The three genetic groups in our data set suggest the potential for cryptic species within the lineage. While many factors can potentially account for the observed differences in genetic structure, including population persistence and large-scale patterns of range expansion, we propose that differences in the mobility of hosts of tick immature stages (small mammals in I. ovatus; birds in H. flava) may be driving the observed patterns.

Parasitology

gene flow

host mobility

Mantel test

Spatial autocorrelation

cryptic species

Tick-borne diseases are a public health concern and their control is often challenging due to the fact that it involves a complex transmission chain of vertebrate hosts and ticks interacting in a changing environment [1]. In order to understand the vector-host relationship, studies on the spatial distribution and population genetic structure of ticks are needed. Population genetic studies have helped in understanding the dispersal patterns and potential pathogen transmission in ticks [2]. The direction, distance and gene flow patterns between tick populations can be inferred from population genetic studies. For example, if high gene flow is observed there might be a greater chance of colonizing new areas or of re-colonizing areas following successful vector control programs [3]. Therefore, estimates of genetic structure and gene flow in tick populations may give insights regarding the spread of tick-borne disease.

Host mobility can affect the genetic patterns of tick populations, although its effects are not consistent. Several studies have reported low levels of gene flow in ticks with less mobile hosts (e.g. smaller mammals) and high levels of gene flow in ticks with highly mobile hosts [2]. For example, Amblyomma americanum and A. triste (Ixodidae, Acari) exhibited high gene flow across various spatial scales (137, 000 km²to 2.78 million km²) and this was attributed to their hosts’ dispersal capabilities (large mammals and birds) [4-6], while low gene flow was observed in Amblyomma dissimile populations and this was attributed to the low mobility of its hosts (small mammals, reptiles and salamanders) [7]. In contrast, some studies have reported low gene flow despite the high mobility of the host, for example in Ixodes scapularis [8] and Ornithodoros coriaceus [9] which can be attributed to the host mobility.

In Japan, tick-borne diseases are an increasing public health concern, affecting humans and animals [10]. A total of 8 genera of ticks have been recorded in Japan, composed of 47 species: 43 Ixodidae species and 4 Argasidae species [11]. Out of 47 species, 21 species are known to parasitize humans [12]. Among these 21 species, Ixodes ovatus, the main vector of Lyme disease [13] and Haemaphysalis flava, a vector of severe fever with thrombocytopenia syndrome (SFTS) and Japanese spotted fever (JSF), were reported in Japan [10,14]. A previous study found that the hosts of adult I. ovatus were mainly hares and can also be large mammals (e.g. cows and horses), and the hosts of immature forms were small rodents [15]. H. flava adult host preference are cows, dogs, horses, wild boar, bear, and deer, while immature forms only parasitized birds [15]. Despite the similarity of I. ovatus and H. flava in which blood meal from different vertebrate host is taken in each life stage (larva, nymph and adult) [16], comparative population genetic studies on two Ixodid tick species such as I. ovatus and H. flava from the same sampling location remain nonexistent. Thus, our study is important because there are no previous studies about dispersal movement of two important tick species: I. ovatus and H. flava, we conclude that low gene flow might be observed in I. ovatus due to the lower mobility of hares in the adult stage and small rodents in its immature stage while H. flava will exhibit high gene flow because of its avian mediated dispersal in its immature stage.

In addition to tick population structure studies, genetic studies may also reveal the presence of cryptic species, where morphologically identified individuals might represent more than one species [17]. Previous studies on Rhipecephalus appendiculatus [18], I. holocylus [19] and I. ovatus [20] have revealed the presence of cryptic species based on the clustering of haplotypes in a phylogenetic tree, and concluded that morphological criteria for species differentiation alone are equivocal and that genetic analysis is important. A putative I. ovatus species complex was identified based on the presence of 3 distinct haplotype clusters in both cox1 and 16S rRNA genes; 2 groups in China and one in Japan that included haplotypes from North America [20].

Here, we studied the population genetic structure of I. ovatus and H. flava and also tested for the presence of cryptic species using mitochondrial DNA sequences of cox1 and the 16S rRNA gene. We hypothesized that I. ovatus and H. flava may have differences in their population genetic structure despite some overlap in their adult host preference because of the differences in mobility of the immature tick hosts. We predicted that H. flava exhibits high gene flow patterns due to avian-mediated dispersal of immature ticks, in contrast to I. ovatus which uses small mammals as hosts during its immature stages.

Study site, collection, sampling and identification

From April 2016 until November 2017, ticks were collected by the standard flagging methods [21] across Niigata Prefecture in Japan and a total of 29 sampling locations were surveyed for H. flava and I. ovatus ticks (Additional File 1. Table S1). Altitude per location ranged from 8 -1402 m/a.s.l. and the geographic distance between sites ranged from 4.28 - 247.65 km. Collected ticks were stored in microcentrifuge tubes with ethanol at 4⁰C. A total of 2,103 individual tick samples were collected. Identification of sex, developmental stage and species was performed using a compound microscope and identification keys of [15].

DNA extraction, PCR amplification and sequencing

Genomic DNA (I. ovatus n=320; H. flava n=220) from each identified tick was extracted using the Isogenome DNA extraction kits (NIPPON GENE Co. Ltd., Tokyo, Japan) following the manufacturer’s recommended protocol. Prior to DNA extraction, each tick was washed with alcohol and PBS solution. DNA quality was checked using a NanoDrop™ 2000 Spectrophotometer (Thermo Scientific™). Two mitochondrial genes were amplified by polymerase chain reaction (PCR): Cytochrome Oxidase 1 (cox1) (658 base pairs) using the primer pairs LCO-1490 (5'- GGTCAACAAATCATAAAGATATTGG -3') and HCO1-2198 (5'- AAACTTCAGGGTGACCAAAAAATCA-3') [22] and the 16S rRNA gene (16S) (407 base pairs) with the following primer pairs 16S+1 (5'CTGCTCAATGAATATTTAAATTGC-3') and 16S – 1 (5' -CGGTCTAAACTCAGATCATGTAGG-3’) [23]. All PCR amplifications of both cox1 and 16S were performed with a final volume of 10 µl with 1µl of genomic DNA. The PCR reaction for both markers were composed of the following: 10x Ex Taq buffer, 25mM MgCl₂, 2.5mM dNTP, 10μm of forward and reverse primers and 5 units/μl of TaKaRa Ex Taq™ (Takara Bio Inc.). The cox1 PCR amplification were as follows: an initial denaturation of 94⁰C for 2 minutes, denaturation of 94⁰C for 30 seconds, annealing of 38⁰C for 30 seconds, extension of 72⁰C for 1 minute for 30 cycles and final extension of 72⁰C for 10 minutes. While the 16S PCR amplification followed the protocol of [23] with some modifications, initial denaturation of 94⁰C for 3 minutes, denaturation of 94⁰C for 30 seconds, annealing of 50⁰C for 40 seconds, extension of 72⁰C for 40 seconds for 30 cycles and final extension of 72⁰C for 5 minutes. PCR products were purified using the QIAquick 96 PCR Purification Kit (Qiagen) in accordance with the manufacturer’s instructions, and sequenced by Eurofin Genomics, Inc., Tokyo, Japan.

Sequence data analysis

We assembled forward and reverse reads for each individual using CodonCode Aligner ver 1.2.4 software (https://www.codoncode.com/aligner/). No ambiguous bases were observed and the low quality bases were removed in the start and end of the reads. Multiple alignment of sequences was performed using the MAFFT alignment online program with default settings (https://mafft.cbrc.jp/alignment/server/). To ensure sequence quality and correct species identification, we checked the similarity of the sequences against reference sequences from GenBank using BLASTN (Basic Local Alignment Search Tool – Nucleotide, https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome). The quality of the final aligned sequences was checked in Mesquite ver 3.5 [24] wherein the protein coding genes were translated into amino acids to confirm the absence of stop codons.

Population genetic analysis

We analyzed the final sequences of both H. flava and I. ovatus per genetic marker separately using DNASp ver 6.12.03 [25] and Arlequin ver 3.5.2.2 [26] for the following parameters: number of haplotypes, average number of polymorphic sites, average number of nucleotide differences, gene/haplotype diversity and nucleotide diversity. The population genetic structure among and within sampling locations in Niigata Prefecture was calculated by analysis of molecular variance (AMOVA) performed in Arlequin with 1023 permutations. The degree of molecular genetic differentiation between sampling locations defined as populations was assessed by calculation of the pairwise F_ST values in Arlequin.

To determine if the genetic differentiation was influenced by geographical distance or altitudinal differences among populations, we performed Mantel Test in GenAlEx ver 6.51b2 [27]. Two tests per species and marker were conducted. First we compared pairwise genetic (pairwise F_ST values) and geographical distances (km) and second we compared genetic distance (F_ST values) with altitudinal differences (m/a.s.l.). The geographic distances were obtained from the GPS coordinates (latitude and longitude) recorded during the sampling. All Mantel tests were assessed using 9999 permutations for the significance of the correlation.

The spatial component of genetic variation was further assessed by spatial autocorrelation [28] using GenAlEx [27]. The autocorrelation coefficient (r) was computed from the geographic distance between sampling locations and the genetic distance (pairwise F_ST values). This measures the genetic similarity between individual pairs within an identified distance class. The size of the distance class influences the estimation of the spatial autocorrelation. We identified the appropriate distance class based on the observed distribution of pairwise geographic distance between the sampling locations (data not shown). The distance class sizes used were the following: 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300, and 320 kilometers. Bootstrap estimates of r and random permutations were set at 9999 each for the test for significance. The upper and lower error bars in the correlogram bound the 95% confidence interval of r value as identified by bootstrap resampling. The upper and lower confidence limits bound the 95% confidence interval of the null hypothesis of no spatial structure in the data set.

The genetic relationships among the sampling locations (i.e., populations) were calculated by unweighted pair group method with the arithmetic mean (UPGMA) cluster method using the APE package [29] and R program [30] to create dendrogram using the genetic distance matrix (pairwise F_ST values) generated from GenAlEx.

Phylogenetic analysis

We constructed maximum likelihood (ML) gene trees for cox1 and 16S haplotype sequences of I. ovatus and H. flava using PhyML ver 3.1 [31] under the default settings. We applied HKY and GTR nucleotide substitution models for cox1 and 16S, respectively, as suggested by jModelTest ver 2 [32]. Additional sequences from China (MH208506, MH208512, MH208514, MH208522, MH208515-19, MH208524, MH208531, MH208574, MH208577, MH208579, MH208681-87, MH208689-93, MH208706, KU664519 [20]), Japan (Hokkaido, Yamanashi and Aomori) (AB231670, AB819241, AB819243, AB819244 [33-34]) and the USA (U95900; [35]) were included. We used Ixodes canisuga as an outgroup (KY962023, KY962074; [36]).

The number of successfully sequenced and analyzed H. flava ranged from 1 to 36 individuals per study site, with more successfully sequenced for cox1 (220/223; 98.7%) than for 16S (172/223; 77.1%). The number of I. ovatus ranged from 1 to 36 individuals per study site and, as with H. flava, more were successfully sequenced for cox1 (307/320; 95.9%) than for 16S (284/320; 88.8%) (Additional File 1. Table S1). There were 63 and 60 cox1 haplotypes in H. flava and I. ovatus, respectively, and 40 and 24 16S haplotypes (Table 1). Haplotype diversity ranged from 0.442 to 0.852 and nucleotide diversity ranged from 0.001 to 0.006 for both markers and species (Table 1).

AMOVA results revealed a high among-population divergence (41.54%) in I. ovatus cox1 sequences, whereas both cox1 and 16S markers of H. flava indicated no significant genetic variation among populations (Table 2). Global F_ST values were highest in I. ovatus with values 0.415 in cox1 and 0.074 in 16S sequences. Pairwise F_ST values (0.500 to 1.00) of I. ovatus based on cox1 showed significant genetic differences between most pairs of populations such as between locations 18 and 25 and between locations 28 and 25 (Additional File 2 Table S2). No significant differences were observed in 16S (Additional File 3 Table S3). Pairwise F_ST values from the cox1 marker of H. flava showed no significant genetic differences except between locations 18 and 15 and 18 and 25 (Additional File 4 Table S4). No significant genetic differentiation was also observed in 16S marker of H. flava except in between locations 3 and 20 and 4 and 20 (Additional File 5 Table S5).

The Mantel test of I. ovatus showed evidence of isolation by geographic distance in both markers: (cox1 r = 0.197, p = 0.008; 16S r = 0.293, p = 0.011) (Figures 1a, b). In contrast, there was no evidence for isolation by altitudinal difference (cox1 r = -0.066, p=0.225; 16S r = -0.023, p = 0.577). Both Mantel tests were not significant for H. flava (distance: cox1 r= -0.145, p= 0.23; 16S r = -0.049, p= 0.33; altitude: r = 0.092, p = 0.30; 16S r = 0.217, p= 0.06). The spatial autocorrelation test in I. ovatus revealed positive and significant correlation between genetic and geographic distance from 20 – 40 km, and 60 – 80 km, and 20 – 50 km and 120 – 160 km from cox1 and 16S, respectively (Figure 1c, d).

The UPGMA cluster dendrogram from the pairwise F_STvalues from the cox1, I. ovatus (Figure 2) revealed 2 genetic clusters across the 29 sampling locations, where the cluster 2 included populations from the more western sites including the island (location 24 and 25) while cluster 1 mostly included populations in the northern and southern part of Niigata Prefecture (Figure 3). No evident genetic clustering on the dendrogram was observed from the I. ovatus 16S sequences and in both cox1 and 16S H. flava sequences (Additional File 6 Figure S1; Additional File 7 Figure S2; Additional File 8 Figure S3).

The cox1 gene tree of 60 I. ovatus haplotypes indicated three groups: group 1 of the published sequences from western China, one large group 2 containing the 58 haplotypes and 2 divergent haplotypes (Hap 60 and Hap 59), and group 3 from western China (Figure 4). The 2 divergent haplotypes were singletons. One was found in location 29 and the other in location 7. The published haplotype from Hokkaido [33] occurred within our large group and the published sequences from western China did not. The I. ovatus 16S gene tree (24 newly sequenced haplotypes) (Additional File 9 Figure S4) also recovered published haplotypes from Japan (Yamanashi and Aomori; [34]) as well as from North America [35] within the same cluster as our new sequences. The cox1 haplotype sequences of H. flava were similar to reference sequences from China (Additional File 10 Figure S5) [unpublished manuscript] while the 16S haplotype sequences (Additional File 11 Figure S6) showed similarities with Japan (Kagoshima, Fukui Aomori, Fukui, Yamanashi, Kagawa and Ehime [34]) and China [unpublished manuscript] reference sequences from Genbank.

Contrasting population genetic structures between H. flava and I. ovatus

Our results supported our hypothesis that H. flava and I. ovatus display contrasting population genetic structure. The noticeable high gene flow in H. flava may be a result of the higher level of mobility of its hosts such as birds in tick’s immature stage and the more pronounced genetic structure and lower gene flow in I. ovatus may be due to the restricted ability of its small mammalian hosts at the immature stage to disperse further. 28 species of birds were previously reported as hosts of H. flava, mostly from the order Passeriformes known as migratory birds [37]. The host for adults of I. ovatus included large mammals such as hares while immature forms are found on small mammals such as small rodents and rats [15]. The movement of small mammals is limited because of environmental and ecological factors such as difficulty in crossing roads and higher possibility of road kill due to vehicle collision and or traffic [38], which may explain the genetic structuring in I. ovatus. The host preference of the adults H. flava are large mammals (e.g. cows, horses, wild boar, deer, bear) and hares, dogs, rodents while the immature forms parasitize birds [15]. Mammals and birds may have wide habitats ranges that may allow maintaining the gene flow of H. flava between the locations in Niigata, as previously observed in Amblyomma americanum populations [5,6,39], and Ixodes ricinus [40]. Birds present the greatest potential for farther dispersal of H. flava that could reduce its genetic structuring [41] that supports the results of our study.

Several alternative factors such as tick behavior, biology and ecology can also affect the tick dispersal movement. For example, long movement patterns, resistance against unfavorable weather conditions causing high survival in their hosts and low host specificity [4] are possible reasons for the low genetic structuring in H. flava populations. It could also be influenced by ecological factors such as changes in land use and cover, wildlife and forestry management as observed in the previous study of [42] in the Ixodes ricinus tick populations from Europe. Previous study on two tick species: Hyalomma rufipes and Amblyomma hebraum also displayed contrasting genetic pattern despite the overlap sharing of highly mobiles hosts such as birds [43] due to the capacity of the immature stages of ticks to survive off the hosts in various habitat conditions [44-46]. Population structure can also be influenced by assortative mating (e.g. Ixodes ricinus) wherein genetically similar individuals tend to mate, rather than random mating resulting in increased genetic divergence [47]. Our study does not have supporting data to test these alternative factors thus it is suggested in future studies to further analyze these factors that could play a role in the tick’s dispersal movement.

Based on previous studies [48-49], isolation by distance suggests that the populations are in equilibrium of migration and genetic drift at the spatial scale of Niigata Prefecture. Our results were congruent with a previous study of Amblyomma ovale that showed positive relationship of geographic and genetic distance at the spatial scale of the whole of Brazil whose main hosts for the immature stage are small rodents while tapir and jaguar in its adult stage [50]. In contrast, H. flava populations showed an absence of isolation by distance. In circumstances when geographic distance can’t explain the genetic structure of tick populations, other factor such as ecological variables (e.g. land use) maybe affecting the genetic structuring within tick populations [7, 40, 51].

A significant and positive spatial autocorrelation of I. ovatus in both cox1 and 16S sequences was observed within the geographic range of 20 – 40 km and 60- 80 km in cox1 sequences and 20 – 50 km and 120 – 160 km range in the 16S sequences. We can infer that the positive r values at the maximum geographic rages of 60-80 km in cox1 and 120 - 160 km in 16S I. ovatus sequences indicate that gene flow is occurring in the wide geographic range. The lack of significance at 40-60 km in cox1 may suggest patch structure indicating genetic similarity or variation [28] as observed on limpet species, Siphonaria raoulensis from New Zealand that also displayed recovery of positive r values at greater distances [52]. However, no autocorrelation was observed in H. flava at any distance class, suggesting that gene flow occurring has no specific pattern favoring certain distances [53-55]. Since the sampling geographic range of our study was within Niigata Prefecture and distances >220 km were not included in our sampling, we suggest that future studies collect samples at a wider geographic scale to more precisely determine the geographic distance at which autocorrelation exists.

16S I. ovatus sequences showed genetic differentiation (pairwise F_ST= 0.5-0.6) in few pairs of sampling locations as compared to the cox1 sequences which may be due to lower variability of the marker thus providing inadequate genetic variation for population genetic analysis, as previously observed in Amblyomma ovale [50] and Rhipicephalus microplus [56-57]. The low variability of the 16S in I. ovatus may not be sufficient to infer intraspecific relationship and is further supported by its low (nd = 0.001) nucleotide diversity as compared in cox1 I. ovatus (nd = 0.004). Thus, it may also explain the grouping observed in the cox1 I. ovatus tree which is absent in the 16S I. ovatus sequences. Despite of these, we used 16S mitochondrial gene because this gene mutates at a rate that is informative for species-level phylogenetic and broad biogeographic inferences [2]

The negative pairwise F_STvalues obtained in H. flava sequences do not have any biological meaning. These results from the unbiased estimator allowing a negative F_STwhen the sample size is small. [58-60]. It can also be noted that we found significant pairwise F_STvalues in 2 sampling location in H. flava cox1 sequences probably due to the fact that both locations have only one tick sample. Although, some sampling sites in this study have small sample size (n < 10), the cox1 and 16S sequences provided valuable insight into the genetic structure and gene flow patterns of the present-day tick populations. Thus, for future research, we recommend to increase the sample size per sampling site.

Species complex formation in I. ovatus cox1 sequences

Our cox1 I. ovatus gene tree showed Japanese individuals to form a group that was distinct from haplotypes from southern China [20]. Despite the low gene flow we found in cox1 I. ovatus, haplotypes found from Niigata were closely related to the published sequence from Hokkaido, Japan which may indicate that these ticks originated from a diverse set of geographical locations in Japan which might be transported by its hosts or is undergoing recent population expansion from northern Japan (Hokkaido) to south (Niigata) or vice versa. We found three groups (China 1, Japan and China 2) and 2 slightly divergent cox1 haplotypes (Hap 60 and 59) in the Japan group of I. ovatus in Niigata. Considering two or more cryptic species can be concealed in one morphologically described species [59], the occurrence of the 3 groups and the divergent haplotypes suggests that I. ovatus may be a species complex. Previous studies have also observed species complexes in Ixodes and Rhipicephalus [19-20,57;62] suggesting that morphological criteria for species differentiation alone are equivocal and that genetic analysis is important.

Two Ixodid ticks, H. flava and I. ovatus exhibited contrasting patterns of genetic structure in Niigata Prefecture, Japan. We suggest that these differences are due to differences in host mobility of the immature life stages of the ticks. The more pronounced population genetic structure of I. ovatus might be influenced by the restricted movement of its small mammalian host in its immature stage, whereas the lack of structure in H. flava may be due to its avian-mediated movement in its immature stage. Even though I. ovatus populations were genetically structured within Niigata, a published haplotype from Hokkaido was also found, indicating that widespread dispersal is possible. The occurrence of three groups and the divergent cox1 haplotypes in I. ovatus emphasizes the need of additional methods in determining the existence of species complex in I. ovatus populations in Japan.

flava – Haemaphysalis flava

ovatus – Ixodes ovatus

AMOVA – Analysis of molecular variance

UPGMA - unweighted pair group method with arithmetic mean

ML – maximum likelihood

F_ST– fixation index

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

All data generated or analyzed during the study are included in this published article and its supplementary additional files. All the newly generated sequences are available in the Genbank database under the accession numbers _______.

Competing interests

The authors declare that they have no competing interests.

Funding

This work is funded by the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (16K00569, 19K21996, 19H02276), the Sumitomo Electric Industries Group Corporate Social Responsibility Foundation and the Research Unit Program of Ehime University. The research was partially supported by the German Academic Exchange Service (DAAD, Programm Projektbezogener Personenaustausch Japan 2018, project 57402018).

Authors’ contributions

MAFR, MG, MM, KW conceptualized and the designed the experiment. MS, TT, RA, SI and MOS designed the sampling collection, collected and identified the tick samples of this study. MD and KT conducted the molecular analyses. MAFR, MG and KW performed the data analysis and wrote the manuscript. All authors read and approve the manuscript.

Acknowledgements

The authors would like to thank Niigata Prefectural Office for their valuable help during the tick sampling collection. Thank you to Dr. Thaddeus Carvajal and Joeselle Serrana for their valuable suggestions in the population genetic analyses in this study and Micanaldo Francisco for the technical assistance in constructing the sampling map of this study.

Dantas-Torres F, Chomel BB, Otranto D. Ticks and tick-borne diseases: a One Health perspective. Trends in parasitology. 2012 Oct 1;28(10):437-46.
Araya-Anchetta A, Busch JD, Scoles GA, Wagner DM. Thirty years of tick population genetics: a comprehensive review. Infection, Genetics and Evolution. 2015 Jan 1;29:164-79.
Mccoy K. The population genetic structure of vectors and our understanding of disease epidemiology. Parasite. 2008Aug;15(3):444–8.
Guglielmone AA, Nava S, Mastropaolo M, Mangold AJ. Distribution and genetic variation of Amblyomma triste (Acari: Ixodidae) in Argentina. Ticks and tick-borne diseases. 2013 Sep 1;4(5):386-90.
Mixson TR, Lydy SL, Dasch GA, Real LA. Inferring the population structure and demographic history of the tick, Amblyomma americanum Linnaeus. Journal of Vector Ecology. 2006 Jun;31(1):181-93.
Trout RT, Steelman CD, Szalanski AL. Population genetics of Amblyomma americanum (Acari: Ixodidae) collected from Arkansas. Journal of medical entomology. 2014 Dec 1;47(2):152-61.
Lampo M, Rangel Y, Mata A. Population genetic structure of a three-host tick, Amblyomma dissimile, in eastern Venezuela. Journal of Parasitology. 1998 Dec 1;84:1137-42.
Qiu WG, Dykhuizen DE, Acosta MS, Luft BJ. Geographic uniformity of the Lyme disease spirochete (Borrelia burgdorferi) and its shared history with tick vector (Ixodes scapularis) in the northeastern United States. Genetics. 2002 Mar 1;160(3):833-49.
Teglas MB, May B, Crosbie PR, Stephens MR, Boyce WM. Genetic structure of the tick Ornithodoros coriaceus (Acari: Argasidae) in California, Nevada, and Oregon. Journal of medical entomology. 2005 May 1;42(3):247-53.
Yamaji K, Aonuma H, Kanuka H. Distribution of tick-borne diseases in Japan: past patterns and implications for the future. Journal of infection and chemotherapy. 2018 Jul 1;24(7):499-504.
Fujita H. Identification of immature ticks in Japan. Acari and Emerging/Reemerging Disease. 2007:53-68.
Okino T, Ushirogawa H, Matoba K, Hatsushika R. Bibliographical studies on human cases of hard tick (Acarina: Ixodidae) bites in Japan (7) cases of unidentified tick infestation. Kawasaki Medical Journal. 2010;36(2):127-41.
Miyamoto K, Nakao M, Fujita H, Sato F. The ixodid ticks on migratory birds in Japan and the isolation of Lyme disease spirochetes from bird-feeding ticks. Medical entomology and zoology. 1993 Dec 15;44(4):315-26.
Yu XJ, Liang MF, Zhang SY, Liu Y, Li JD, Sun YL, Zhang L, Zhang QF, Popov VL, Li C, Qu J. Fever with thrombocytopenia associated with a novel bunyavirus in China. New England Journal of Medicine. 2011 Apr 21;364(16):1523-32.
Yamaguti N, Tipton VJ, Keegan HL, Toshioka S. Ticks of Japan, Korea, and the Ryukyu islands. Brigham Young University Science Bulletin, Biological Series. 1971;15(1):1.
Herrmann C, Gern L. Search for blood or water is influenced by Borrelia burgdorferi in Ixodes ricinus. Parasites & vectors. 2015 Dec;8(1):6.
Fegan M, Prior P. How complex is the Ralstonia solanacearum species complex. Bacterial wilt disease and the Ralstonia solanacearum species complex. 2005 Jan 1;1:449-61.
Kanduma EG, Mwacharo JM, Githaka NW, Kinyanjui PW, Njuguna JN, Kamau LM, Kariuki E, Mwaura S, Skilton RA, Bishop RP. Analyses of mitochondrial genes reveal two sympatric but genetically divergent lineages of Rhipicephalus appendiculatus in Kenya. Parasites & vectors. 2016 Dec;9(1):353.
Song S, Shao R, Atwell R, Barker S, Vankan D. Phylogenetic and phylogeographic relationships in Ixodes holocyclus and Ixodes cornuatus (Acari: Ixodidae) inferred from COX1 and ITS2 sequences. International journal for parasitology. 2011 Jul 1;41(8):871-80.
Li LH, Zhang Y, Wang JZ, Li XS, Yin SQ, Zhu D, Xue JB, Li SG. High genetic diversity in hard ticks from a China-Myanmar border county. Parasites & vectors. 2018 Dec;11(1):469.
Ginsberg HS, Ewing CP. Comparison of flagging, walking, trapping, and collecting from hosts as sampling methods for northern deer ticks, Ixodes dammini, and lone-star ticks, Amblyomma americanum (Acari: Ixodidae). Experimental & applied acarology. 1989 Sep 1;7(4):313-22.
Folmer O, Black, M, Hoeh W, Lutz R and Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology. 1994;(3):294-9.
Tian Z, Liu G, Xie J, Yin H, Luo J, Zhang L, Zhang P, Luo J. Discrimination between Haemaphysalis longicornis and H. qinghaiensis based on the partial 16S rDNA and the second internal transcribed spacer (ITS-2). Experimental and applied acarology. 2011 Jun 1;54(2):165-72.
Maddison WP, Maddison DR. Mesquite: a modular system for evolutionary analysis. 2011. See http://mesquiteproject. org. 2011.
Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular biology and evolution. 2017 Sep 18;34(12):3299-302.
Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary bioinformatics. 2005 Jan;1:117693430500100003.
Peakall RO, Smouse PE. GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular ecology notes. 2006 Mar;6(1):288-95.
Smouse PE, Peakall RO. Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity. 1999 May;82(5):561-73.
Paradis E, Schliep K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2018 Jul 17;35(3):526-8.
R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2016.
Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Systematic biology. 2003 Oct 1;52(5):696-704.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nature methods. 2012 Aug;9(8):772.
Mitani H, Takahashi M, Masuyama M, Fukunaga M. Ixodes philipi (Acari: Ixodidae): phylogenetic status inferred from mitochondrial cytochrome oxidase subunit I gene sequence comparison. The Journal of parasitology. 2007 Jun 1:719-22.
Takano A, Fujita H, Kadosaka T, Takahashi M, Yamauchi T, Ishiguro F, Takada N, Yano Y, Oikawa Y, Honda T, Gokuden M. Construction of a DNA database for ticks collected in Japan: application of molecular identification based on the mitochondrial 16S rDNA gene. Medical Entomology and Zoology. 2014 Mar 25;65(1):13-21.
Norris DE, Klompen JS, Black WC. Comparison of the mitochondrial 12S and 16S ribosomal DNA genes in resolving phylogenetic relationships among hard ticks (Acari: Ixodidae). Annals of the Entomological Society of America. 1999 Jan 1;92(1):117-29.
Hornok S, Sándor AD, Beck R, Farkas R, Beati L, Kontschán J, Takács N, Földvári G, Silaghi C, Meyer-Kayser E, Hodžić A. Contributions to the phylogeny of Ixodes (Pholeoixodes) canisuga, I.(Ph.) kaiseri, I.(Ph.) hexagonus and a simple pictorial key for the identification of their females. Parasites & vectors. 2017 Dec;10(1):545.
Yamauchi T, Takeno K. Record of Haemaphysalis flava Neumann (Acari, Ixodidae) from Zoothera dauma (Latham)(Aves) in Kyushu, Japan. Esakia. 2000 Mar 31;40:99-101.
Haemig PD, Waldenström J, Olsen B. Roadside ecology and epidemiology of tick-borne diseases. Scandinavian journal of infectious diseases. 2008 Jan 1;40(11-12):853-8.
Reichard MV, Kocan AA, Van Den Bussche RA, Barker RW, Wyckoff J, Ewing SA. Sequence variation of the ribosomal DNA second internal transcribed spacer region in two spatially distinct populations of Amblyomma americanum (L.)(Acari: Ixodidae). Journal of Parasitology. 2005 Apr;91(2):260-4.
Casati S, Bernasconi MV, Gern L, Piffaretti JC. Assessment of intraspecific mtDNA variability of European Ixodes ricinus sensu stricto (Acari: Ixodidae). Infection, Genetics and Evolution. 2008 Mar 1;8(2):152-8.
Mclain DK, Wesson DM, Oliver Jr JH, Collins FH. Variation in ribosomal DNA internal transcribed spacers 1 among eastern populations of Ixodes scapularis (Acari: Ixodidae). Journal of Medical Entomology. 1995 May 1;32(3):353-60.
Medlock JM, Hansford KM, Bormane A, Derdakova M, Estrada-Peña A, George JC, Golovljova I, Jaenson TG, Jensen JK, Jensen PM, Kazimirova M. Driving forces for changes in geographical distribution of Ixodes ricinus ticks in Europe. Parasites & vectors. 2013 Dec 1;6(1):1.
Cangi N, Horak IG, Apanaskevich DA, Matthee S, Das Neves LC, Estrada-Peña A, Matthee CA. The influence of interspecific competition and host preference on the phylogeography of two African Ixodid tick species. PLoS One. 2013 Oct 9;8(10):e76930.
Needham GR, Teel PD. Water balance by ticks between bloodmeals. Morphology, physiology, and behavioral biology of ticks/editors, John R. Sauer and J. Alexander Hair. 1986.
Cumming GS. Host distributions do not limit the species ranges of most African ticks (Acari: Ixodida). Bulletin of Entomological Research. 1999 Apr;89(4):303-27.
Estrada-Peña A. Ticks as vectors: taxonomy, biology and ecology. Revue scientifique et technique. 2015 Apr 1;34(1):53-65.
Kempf F, De Meeûs T, Arnathau C, Degeilh B, Mccoy KD. Assortative pairing in Ixodes ricinus (Acari: Ixodidae), the European vector of Lyme borreliosis. Journal of medical entomology. 2014 Oct 29;46(3):471-4.
Yaegashi S, Watanabe K, Monaghan MT, Omura T. Fine-scale dispersal in a stream caddisfly inferred from spatial autocorrelation of microsatellite markers. Freshwater Science. 2014 Feb 4;33(1):172-80.
Watanabe K, Monaghan MT, Takemon Y, Omura T. 1775201. Dispersal ability determines the genetic effects of habitat fragmentation in three species of aquatic insect. Aquatic conservation. 2010;20(5):574-9.
Bitencourth K, Amorim M, de Oliveira SV, Voloch CM, Gazêta GS. Genetic diversity, population structure and rickettsias in Amblyomma ovale in areas of epidemiological interest for spotted fever in Brazil. Medical and veterinary entomology. 2019 Jun;33(2):256-68.
Krakowetz CN, Dergousoff SJ, Chilton NB. Genetic variation in the mitochondrial 16S rRNA gene of the American dog tick, Dermacentor variabilis (Acari: Ixodidae). Journal of Vector Ecology. 2010 Jun;35(1):163-73.
Wood AR, Gardner JP. Small spatial scale population genetic structure in two limpets (Scutellastra kermadecensis and Siphonaria raoulensis) endemic to the Kermadec Islands, New Zealand. Marine Ecology Progress Series. 2007;349:159-70.
Pardanani N, Mather TN. Lack of spatial autocorrelation in fine-scale distributions of Ixodes scapularis (Acari: Ixodidae). Journal of medical entomology. 2004 Sep 1;41(5):861-4.
Kitron U, Kazmierczak JJ. Spatial analysis of the distribution of Lyme disease in Wisconsin. American journal of Epidemiology. 1997 Mar 15;145(6):558-66.
Nicholson MC, Mather TN. Methods for evaluating Lyme disease risks using geographic information systems and geospatial analysis. Journal of Medical Entomology. 2014 Sep 11;33(5):711-20.
Low VL, Tay ST, Kho KL, Koh FX, Tan TK, Lim YA, Ong BL, Panchadcharam C, Norma-Rashid Y, Sofian-Azirun M. Molecular characterisation of the tick Rhipicephalus microplus in Malaysia: new insights into the cryptic diversity and distinct genetic assemblages throughout the world. Parasites & vectors. 2015 Dec;8(1):341.
Burger TD, Shao R, Barker SC. Phylogenetic analysis of mitochondrial genome sequences indicates that the cattle tick, Rhipicephalus (Boophilus) microplus, contains a cryptic species. Molecular Phylogenetics and Evolution. 2014 Jul 1;76:241-53.
Weir BS, Cockerham CC. Estimating F‐statistics for the analysis of population structure. evolution. 1984 Nov;38(6):1358-70.
Bitencourth K, Voloch CM, Serra‐Freire NM, Machado‐Ferreira E, Amorim M, Gazêta GS. Analysis of Amblyomma sculptum haplotypes in an area endemic for Brazilian spotted fever. Medical and veterinary entomology. 2016 Sep;30(3):342-50.
Vial L, Durand P, Arnathau C, Halos L, Diatta G, Trape JF, Renaud F. Molecular divergences of the Ornithodoros sonrai soft tick species, a vector of human relapsing fever in West Africa. Microbes and infection. 2006 Sep 1;8(11):2605-11.
Kwon-Chung KJ, Bennett JE, Wickes BL, Meyer W, Cuomo CA, Wollenburg KR, Bicanic TA. The Case for Adopting the" Species Complex" Nomenclature for the Etiologic Agents of Cryptococcosis. Sphere 2017; 2 (1): e00357-16. DOI: https://doi. org/10.1128/mSphere.:00357-16.
Xu G, Fang QQ, Keirans JE, Durden LA. Molecular phylogenetic analyses indicate that the Ixodes ricinus complex is a paraphyletic group. The Journal of parasitology. 2003 Jun 1:452-7.

Table 1. Summary of cox1 and 16S haplotype distribution, variability, genetic diversity of Ixodes ovatus and Haemaphysalis flava populations in Niigata Prefecture

Marker	Species	n	h	s	k	hd±SD	nd±SD
cox1	I. ovatus	307	60	65	2.728	0.852 ± 0.012	0.004 ± 0.000
	H. flava	220	63	60	1.472	0.789 ± 0.026	0.002 ± 0.000
16S	I. ovatus	284	24	22	0.699	0.442 ± 0.037	0.001 ± 0.000
	H. flava	172	40	49	2.447	0.835 ± 0.024	0.006 ± 0.000

Abbreviations: n sample size; h number of haplotypes; s number of polymorphic sites; k mean number of nucleotide differences; hd±SD haplotype diversity ± Standard Deviation; nd±SD nucleotide diversity ± Standard Deviation

Table 2. Analysis of molecular variance (AMOVA) using cox1 and 16S of Ixodes ovatus and Haemaphysalis flava populations

Marker	Species		df	ss	vc	pv	F_ST
cox1	I. ovatus	Among	28	190.38	0.5803 Va	41.54	0.4154**
	I. ovatus	Within	278	227.04	0.8167 Vb	58.46	0.4154**
	H. flava	Among	16	12.42	0.0039 Va	0.53	0.0053
	H. flava	Within	203	148.73	0.7327 Vb	99.47	0.0053
16S	I. ovatus	Among	22	14.04	0.0260 Va	7.39	0.0739**
	I. ovatus	Within	261	84.87	0.3252 Vb	92.61	0.0739**
	H. flava	Among	12	23.6	0.0642 Va	5.21	0.0521**
	H. flava	Within	159	185.58	1.1672 Vb	94.79	0.0521**

Abbreviations: df degrees of freedom; ss sum of squares; Va, Vb and Vc are associate covariance components; pv percentage variation; *P<0.05 ** P<0.01

Additional File 1.xls Table S1. Summary of Haemaphysalis flava and Ixodes ovatus collected from the different locations of Niigata Prefecture and its corresponding sample number

Additional File 2.xls Table S2. Pairwise comparison of genetic differentiation (F_ST) calculated for all I. ovatus using the cox1 gene

Additional File 3.xls Table S3. Pairwise comparison of genetic differentiation (F_ST) calculated for all I. ovatus using the 16S gene

Additional File 4.xls Table S4. Pairwise comparison of genetic differentiation (F_ST) calculated for all H. flava using the cox1 gene

Additional File 5.xls Table S5. Pairwise comparison of genetic differentiation (F_ST) calculated for all H. flava using the 16S gene

Additional File 6.docx Figure S1. An unweighted pair group method with the arithmetic mean (UPGMA) dendrogram of I. ovatus based on the pairwise genetic distance (F_ST) of 16S among 23 sampling locations across Niigata Prefecture, Japan.

Additional File 7.docx Figure S2 An unweighted pair group method with the arithmetic mean (UPGMA) dendrogram of H. flava based on the pairwise genetic distance (F_ST) of cox1 among 17 sampling locations across Niigata Prefecture, Japan.

Additional File 8.docx Figure S3 An unweighted pair group method with the arithmetic mean (UPGMA) dendrogram of H. flava based on the pairwise genetic distance (F_ST) of 16S among 17 sampling locations across Niigata Prefecture, Japan.

Additional File 9.docx Figure S4 Phylogenetic tree of 16s haplotype sequences of I. ovatus inferred by maximum likelihood analysis with evolutionary model GTR with sequences from China (MH208506, MH208512, MH208514 – MH208519, MH208522, MH208524, MH208531, MH208574, MH208577, MH208579, KU664519), Japan (AB819241, AB819242, AB819243, AB819244) and US (U95900). An outgroup (Ixodes canisuga) was also included.

Additional File 10.docx Figure S5 Phylogenetic tree of cox1 haplotype sequences of H. flava inferred by maximum likelihood analysis with evolutionary model HKY with sequences from China (KY021800- KY021807; KY021810-KY021819, KY003181, JQ62588 - JQ625889; JQ737097; JF758632). An outgroup (Ixodes canisuga) was also included.

Additional File 11.docx Figure S6 Phylogenetic tree of 16S haplotype sequences of H. flava inferred by maximum likelihood analysis with evolutionary model GTR with sequences from China (KC844858 -KC844867; KC844880-KC844882, M696720, KP324926, KX450279) and Japan (AB19177- AB819192). An outgroup (Ixodes canisuga) was also included.

SupplementaryFiles.zip

PDF

Version 1

posted 11 Dec, 2019

You are reading this latest preprint version

Comparative population genetic structure of two Ixodidae ticks (Ixodes ovatus and Haemaphysalis flava) in Niigata Prefecture, Japan

Status:

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusion

List Of Abbreviations

Declarations

References

Tables

Additional Files Description

Supplementary Files

Status:

Version 1