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 40C. 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 MgCl2, 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 940C for 2 minutes, denaturation of 940C for 30 seconds, annealing of 380C for 30 seconds, extension of 720C for 1 minute for 30 cycles and final extension of 720C for 10 minutes. While the 16S PCR amplification followed the protocol of [23] with some modifications, initial denaturation of 940C for 3 minutes, denaturation of 940C for 30 seconds, annealing of 500C for 40 seconds, extension of 720C for 40 seconds for 30 cycles and final extension of 720C 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 FST 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 FST values) and geographical distances (km) and second we compared genetic distance (FST 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 FST 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 FST 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]).