Comparative mitogenomics elucidates the population genetic structure of Amblyomma testudinarium in Japan and a closely related Amblyomma species in Myanmar

Abstract Ticks are the second most important vector capable of transmitting diseases affecting the health of both humans and animals. Amblyomma testudinarium Koch 1844 (Acari: Ixodidae), is a hard tick species having a wide geographic distribution in Asia. In this study, we analyzed the composition of A. testudinarium whole mitogenomes from various geographical regions in Japan and investigated the population structure, demographic patterns, and phylogeographic relationship with other ixodid species. In addition, we characterized a potentially novel tick species closely related to A. testudinarium from Myanmar. Phylogeographic inference and evolutionary dynamics based on the 15 mitochondrial coding genes supported that A. testudinarium population in Japan is resolved into a star‐like haplogroup and suggested a distinct population structure of A. testudinarium from Amami island in Kyushu region. Correlation analysis using Mantel test statistics showed that no significant correlation was observed between the genetic and geographic distances calculated between the A. testudinarium population from different localities in Japan. Finally, demographic analyses, including mismatch analysis and Tajima’s D test, suggested a possibility of recent population expansion occurred within Japanese haplogroup after a bottleneck event. Although A. testudinarium has been considered widespread and common in East and Southeast Asia, the current study suggested that potentially several cryptic Amblyomma spp. closely related to A. testudinarium are present in Asia.

gions (Takada, 1990). In recent years, A. testudinarium was infesting on grazing cattle and brown bears in the northern parts of Japan such as Aomori and Hokkaido prefectures, which were suspected to be carried by the migratory birds Terada et al., 2019). Amblyomma testudinarium is known as a vector for R. tamurae, a causative agent of spotted fever rickettsiosis in humans (Fournier et al., 2006;Imaoka et al., 2011), and SFTS virus in Japan.
Recently, a novel thogotovirus, named as Oz virus, was isolated from A. testudinarium collected in Ehime prefecture, Japan (Ejiri et al., 2018). Cases of human bites have been reported elsewhere (Isohisa et al., 2011;Nakao et al., 2017). Moreover, the cases of tickassociated rash illness (TARI) showing manifestation of erythema migrans (EM) without association of Lyme disease have been reported in the patients who were bitten by this tick species (Natsuaki et al., 2014).
High-throughput sequencing techniques are very functional in species identification and population genetics of ticks (Burger et al., 2012(Burger et al., , 2013(Burger et al., , 2014. Comparative whole mitochondrial genome (mitogenome) sequence analyses have been used to identify key genetic features among tick species (Kelava et al., 2021;Regilme et al., 2021;Wang et al., 2019). Genetic characterization of complete mitogenome sequences is a core study nowadays due to its highly conserved structure, the high mutation rate, low recombination rate, and maternal inheritance. Recently, two studies reported the complete whole mitogenome of A. testudinarium, which enables the comparative analysis in a higher resolution between tick species (Chang et al., 2020;Nakao et al., 2021) than previously implemented intraspecies phylogenetic analyses targeting partial sequences of 16S ribosomal RNA gene (rDNA), 12S rDNA, and cytochrome c oxidase I gene (COI) (Chao et al., 2017;Yamauchi et al., 2012).
In the present study, we describe the population structure of A. testudinarium based on whole mitogenome sequences obtained from Japan. We investigated the relationships between genetic distances of A. testudinarium samples and the corresponding geographical locations from which they were collected. No geographic clustering was detected among A. testudinarium sequences from Japan, except the samples from Amami island in Kyushu region. Additionally, we characterized a potentially novel Amblyomma species from Myanmar and investigated its genetic relationships with A. testudinarium based on the parsimony informative sites of their mitogenomes.

| Specimen collection and DNA extraction
A total of 39 tick specimens were included in this study. Thirty ticks were collected from 22 collection sites in four geographical blocks (Kinki, Chugoku, Shikoku, and Kyushu) in Japan and nine specimens from Putao town in Myanmar (Table 1, Figure 1). All tick samples from Japan were collected from the vegetations (n = 30), while those from Myanmar were collected from water buffaloes (Bubalus bubalis) (n = 5) and cattle (Bos taurus) (n = 4) ( Table 1). The collected ticks were kept at 4°C until transferred to the laboratory for further examinations. Species were identified morphologically using standard taxonomic keys (Yamaguti et al., 1971) and a related article (Chamuah et al., 2016) under a stereomicroscope. Then, each tick was cut in half with a sterile stainless-steel blade. One half was crushed with stainless-steel beads using a Micro Smash MS-100R (TOMY, Tokyo, Japan) at 2500 rpm, and the DNA was extracted using a blackPREP Tick DNA/RNA Kit (Analytikjena, Germany) as described previously (Thu et al., 2019).

| Mitogenome comparison and phylogeny construction
Our complete mitogenome sequences were imported to Geneious and two rDNA sequences were extracted from each mitogenome sequence and concatenated together.

Sequences of the concatenated 15 mitochondrial genes of
Amblyomma samples collected from Japan and Myanmar along with those obtained from the GenBank were aligned using MAFFT software (Katoh & Standley, 2013). A substitution model was selected using PHYML 3.0 software relying on the Akaike Information Criterion (Guindon et al., 2010). A pairwise identity analysis was conducted by calculating the percentage of residues that are identical between sequences in Geneious software.
Bayesian phylogenetic tree was constructed using BEAST version 1.4. as a cross-platform program for Bayesian analysis of molecular sequences using Markov chain Monte Carlo (MCMC). We modeled the evolution of sequences by using the GTR nucleotide substitution model with discrete gamma-distributed rate variation and assuming a constant evolution rate all over the tree by selecting the strict clock model. We assume that the Bayesian skyline coalescent model was a demographic model in a Bayesian framework. The Bayesian skyline plot analysis was performed using a chain length of 50 million generations sampled every 50,000 MCMC steps with a pre-burn-in of 500,000. Maximum clade credibility (MCC) was selected by using the TreeAnnotator (Drummond & Rambaut, 2007).

| Population genetic structure analyses
The analysis of molecular variance (AMOVA) (Excoffier et al., 2007) in Arlequin software version 3.5.2.2 was used to evaluate the genetic variance among and within the populations of A. testudinarium collected from Amami island in Kyushu region vs other localities in Japan. We adjusted the number of permutations at 1000, and the difference was considered significant when tested at p < 0.05 level based on the calculated fixation indices (F-statistics). F ST estimates the degree of differentiation within the population where the closer F ST is to 0, the greater the extent of allelic fixation or identity within populations (Holsinger & Weir, 2009). F SC estimates the differentiation among populations within the group to which they are assigned.
The closer F SC is to 1, the more heterogeneity among populations exists. In case a strong population genetic structure exists at the population scale being analyzed, F SC should be high relative to F ST .

| Spatial and demographic evolutionary dynamics analyses
The historical expansion dynamics of A. testudinarium populations were estimated by calculating the mismatch distribution based on 15 concatenated mitochondrial genes in Arlequin software (Rogers & Harpending, 1992). We compared the observed and expected distributions of pairwise nucleotide variations between individuals.
Multimodal distribution is thought to be associated with a constant population size, whereas unimodal distribution is representing a sudden expansion. The deviation of observed mismatch from the individual expansion model was assessed by estimating the sum of square deviation (SSD); however, the Harpending's raggedness index (RI) was applied to assess the smoothness of the haplotype frequency distribution (Harpending, 1994) where significant RI indicates a poor fit of the data to the expansion model.
Departure from the neutral model of evolution was evaluated through calculating Tajima's D statistics (Tajima, 1989)  Mantel test was used to test whether there is a geographic separation between A. testudinarium populations from different geographic origins and to evaluate the correlation between genetic variation and geographic distance of the tick population. R function "mantel.rtest" in ape package was used (Wynn et al., 2016). Scatter plot of the correlation was created using the R package ggplot2 (Ginestet, 2011;Villanueva & Chen, 2019).

| Data accessibility
The complete mitogenome sequences of 30 A. testudinarium and 9 Amblyomma sp. were deposited in the DNA Data Bank of Japan (http:// www.ddbj.nig.ac.jp) with accession numbers of LC554761-LC554790 (samples from Japan) and LC633546-LC633554 (samples from Myanmar).

| Features of the sequenced mitochondrial genomes
The length of assembled mitogenomes from 30 Amblyomma specimens from Japan and 9 from Myanmar ranged from 14,830 to 14,839 bp. A total of 37 encoded genes were detected including 13 PCGs, 22 transfer RNA (tRNA) genes, and two rDNAs. Two control regions were identified in the mitogenome. There was no difference in gene rearrangement among the 39 newly sequenced Amblyomma mitogenomes.

| Nucleotide identities based complete mitogenome sequences
Generated nucleotide sequences of 30 collected ticks from Japan

| Haplotyping based on the 16S rDNA of A. testudinarium from Japan
To understand the genetic population structure of A. testudinarium in Japan, we performed a haplotype analysis based on 32 16S rDNA sequences from Japan. The haplotype analysis resolved a total of 16 haplotypes, which were defined by 28 single nucleotide polymorphic sites (SNPs). Median-joining phylogenetic network analysis for the 16 haplotypes revealed that the A. testudinarium population from Japan forms a star-like phylogeny, with the common ancestorial haplotype (H5) in the center surrounded by the remaining haplotypes ( Figure 4). The most frequent three haplotypes in the population from Japan were H2, H5, and H9 that were observed in six, six, and four samples, respectively. A total of five (H1, H2, H5, H7, and H9) were shared between different sampling sites ( Figure 4).

| Population structure of A. testudinarium
The phylogeographic analysis based on the 15 concatenated mitochondrial genes extracted from 32 A. testudinarium from Japan and 9 Amblyomma sp. from Myanmar, where A. javanense was used as an outgroup, showed the separation of Japan population of A.
testudinarium from the Amblyomma sp. sequences from Myanmar.
Interestingly, A. testudinarium sequences from Amami island and two sequences from Miyazaki prefecture were grouped and separated from the remaining Japan populations ( Figure 5). In addition, one A.
testudinarium sequence from Ehime prefecture was separated from the remaining sequences from Japan.  Figure 3b). In addition, the molecular diversity indices showed that A. testudinarium population in Japan is more genetically diverse than the Amblyomma sp. of Myanmar (Table S1). In addition, the nonmetric multidimension scale analysis based on the genetic distances confirmed the geographic separation and genetic distinction of A. testurinarium from Amami island and Miyazaki (two samples collected at site 2), which were clustered together in one group and separated from other samples from Japan ( Figure 6). We also showed the genetic distinction and possible geographical separation of A. testurinarium from Ehime prefecture based on the NMDS analysis.
To statistically test the population genetic structure between and within populations, we used AMOVA to determine whether the genetic structure was influenced by the geographic structure.
Molecular variance results revealed that the genetic variation among-population of A. testudinarium from Amami island (Kyushu) and the other regions in Japan indicated that the population of A. testudinarium is structured with one population present only in Amami (Kyushu), where the among-population genetic variation (57.55%) was higher than that within-population (42.45%).
Significant global F ST values (p < 0.01) indicated significant genetic differentiation between populations from Amami and the remaining localities (p = 0.00098) ( Table 4). The correlation analysis of A.
testudinarium from Japan using Mantel test revealed that there is no significant correlation between the degree of genetic variation based on the pairwise genetic distances and geographic distances with r statistics = 0.11 (p-value = 0.1779 based on 9999 replicates) ( Figure 7).
Results of the mismatch analysis based on pairwise nucleotide differences between 41 nucleotide sequences from Japan (n = 32) and Myanmar (n = 9) revealed that the shape of the mismatch distributions is far from unimodal distribution while analyzing either sequences of

| DISCUSS ION
Phylogeographic structure of tick populations can help enrich the knowledge on the distribution of ticks and tick-borne pathogens (Black et al., 2001). In addition, intraspecies genetic differences  (Burnard & Shao, 2019;Mechai et al., 2013;Sakamoto et al., 2014). The current work is the first to assess the genetic diversity and population structure of A. testudinarium, which is one of the most widely distributed Amblyomma in East and Southeast Asia (Voltzit, 2002) and one of the main ticks that bite humans in Japan (Natsuaki, 2021). Our study is the second to use whole mitogenomes of multiple individuals of the same tick species to assess the genetic diversity based on the regional scale after the report on Ixodes ricinus in Italy (Carpi et al., 2016). Although the earlier study detected intraspecies genetic relationships, population differentiation, and demographic history of I. ricinus, our study provides an advantage of the wide spatial scale over which our samples were collected and phylogenetic tree reconstruction using 30 new A. testudinarium  Mitogenomes have been used as a marker in population genetics, phylogeography, and DNA barcoding (Eimanifar et al., 2018;Fourdrilis et al., 2016). Of note, mitogenomes were proven to provide a better opportunity to understand tick taxonomy than targeting only a single gene (Kelava et al., 2021). Previously, intraspecies phylogenetic relationships was investigated through targeting the partial sequences of 16S rDNA, 12S rDNA, and COI (Chao et al., 2017;Yamauchi et al., 2012). Since different genes have different evolutionary rates, the final conclusions could have been affected by the targeted gene topology of the phylogenetic tree (Liu et al., 2018).
Traditionally, morphological characterization has been used for species identification and differentiation of Amblyomma ticks (Aguilar-Dominguez et al., 2021;Namgyal et al., 2021). Now, a molecular based characterization can provide insights on the genetic variance at the base-pair level and the genetic diversity between and within species of Amblyomma ticks (Burger et al., 2012;Lopes et al., 2016).
An important advantage of targeting tick mitogenomes to study its population genetics and distribution is the significantly smaller size when compared to nuclear genomes (Nadolny et al., 2016), which enables the cost-effective comparative genome analysis. We observed that the gene order of the mitogenome of A. testudinarium from Japan and Amblyomma sp. from Myanmar was identical to that of the other members of genus Amblyomma, which confirms that gene order and composition are highly conserved between Amblyomma species to a common ancestor ( F I G U R E 6 Nonmetric multidimensional scaling (NMDS) ordination plot based on the distance matrix. Red, blue, and cyan ellipsoids represent the grouping of Amblyomma testudinarium from Amami, Miyazaki (two samples collected at site 2), and the remaining Japan localities, respectively TA B L E 4 Analysis of molecular variance (AMOVA) using the concatenated 15 mitochondrial genes sequences extracted from whole mitogenomes of Amblyomma testudinarium populations in Japan These results are supported by a previous report , where the mitogenome of A. testudinarium from Hokkaido, the northernmost main island of Japan, was clustered with the same species from Nara in the Honshu Island in a monophyletic group, but obviously discriminated from other species of Amblyomma ticks as well other ixodid genera (Haemaphysalis, Archaeochroton, Bothriocroton, Dermacentor, Robertsicus, and Rhipicephalus) of ticks.
The results of polymorphic sites comparison showed that relying solely on one single mitochondrial gene for phylogenetic analysis of A. testudinarium might yield incorrect or incomplete conclusions (Table 3). Currently, most studies focusing on population structure and phylogenetic analysis of ticks build the final conclusions based on the results obtained by targeting more than one mitochondrial genes (Gui et al., 2021;Regilme et al., 2021). Furthermore, we showed that A. testudinarium population in Japan is genetically more diverse than the closely related Amblyomma sp. of Myanmar.
In fact, the evolution of mitochondrial genes can be affected by several ecological factors such as endosymbionts and animal host diversity (Kaufman et al., 2018;Sassera et al., 2006). The difference in diversity between A. testudinarium in Japan and Amblyomma sp.
in Myanmar could be attributed to the greater breadth of sampling collection sites in Japan compared to Myanmar where we collected the samples from one site. Unfortunately, we did not investigate the bacterial endosymbionts in our study, so we suggest that more studies are required to correlate the endosymbiont genetics and the mitogenome divergence in A. testudinarium.
Amblyomma testudinarium from China (accession number: MT029329) showed low sequence similarity with those from Japan but relatively high similarity with A. javanense reported from China (Table 2). This was supported by the analysis based on the 16S rDNA obtained from A. testudinarium, which included samples from Taiwan, Japan, and Thailand reported by other studies (Chao et al., 2017;Malaisri et al., 2015). We confirmed that the sequence from China (MT029329) is different from A. testudinarium from Japan, Taiwan, and Thailand ( Figure 2). Collectively, these data suggest a possibility of misidentification of the tick specimen used in the previous report (Chang et al., 2020) or the presence of cryptic species related to A.
javanense in China.
The 16S rDNA sequences of A. testudinarium from Japan are clustered together with those from Taiwan but not with Thailand In general, geographical variations are mostly accompanied with differences in the distribution of tick hosts. The large animals are the primary hosts for A. testudinarium, while larvae and nymphs feed on small mammals, birds, reptiles, and amphibians (Kitaoka, 1974;Takada, 1990;Yamaguti et al., 1971). In Japan, adult A. testudinarium are frequently found on wild boars, which habitats are linked to those of A. testudinarium (Motoi et al., 2013). The rapid expansion of the wild boar habitats (Kodera et al., 2001) may be the reason for no clear geographic substructuring among A. testudinarium populations in Japan ( Figure 6). However, future studies with larger numbers of  (Githeko et al., 2000). We also detected a pattern for the population differentiation of A. testudinarium in accordance with the spatial isolation where specimens collected from Amami island might have its own population in Japan.
This finding could be explained by the complete geographical separation from the Japan mainland. The climate around Amami island is humid subtropical and known as Köppen climate classification (Cfa) with high precipitation throughout the year. Furthermore, the life cycle of A. testudinarium in Amami island is not clear; however, it is possible that this tick species utilizes endogenous animal species of this island as feeding hosts (Chae et al., 2017). Any of these factors may have supported the distinction of A. testudinarium tick population in Amami island.
The algorithms in our study support the expansion history of A. testudinarium in Japan. We detected a relatively high haplotype diversity in Japan, suggesting a strong possibility of sudden expansion (Amzati et al., 2018;Ray et al., 2003;Simonsen et al., 1995).
By contrast, Amblyomma sp. in Myanmar was suggested to have a bottleneck historical event. Although our study did not estimate this sudden expansion and bottleneck times, relatively high diversity of A. testudinarium in Japan within population indicates that A. testudinarium colonies have been existed in Japan for long time. However, our study included samples that were collected from only a single sampling site in Myanmar, and we strongly suggest further investigations using A. testudinarium and the closely related Amblyomma spp.
from all over the country to find a possible explanation of the population equilibrium observed in our study. Our finding that A. testudinarium in Japan has experienced a sudden expansion event agrees with the expected movement of migratory birds that are considered the carrying hosts of ticks in the country (Ishiguro et al., 2000;Kuo et al., 2017). This can be supported by the possible establishment of new colonies in previously unoccupied geographical locations in Japan, such as the recent detection of A. testudinarium in Hokkaido Island for the first time in history .

| CON CLUS ION
Overall, we report whole mitogenome sequences of 30 A. testudinarium from Japan and 9 closely related Amblyomma species from Myanmar. Our study provides the first molecular evidence and convincing sequences to confirm the genetic identity and diversity of Amblyomma sp. in Myanmar. The lack of geographic substructuring among A. testudinarium populations in Japan suggests that geographic distance might not have an impact on the genetic variation.
The data provides information about the sudden expansion event of A. testudinarium in Japan. With the advancement of high-throughput sequencing techniques and bioinformatics, phylogeographic approaches can be used to distinguish A. testudinarium globally based on the genetic data of the whole mitochondrial genomes.
The accumulation of more mitogenome sequences from different geographic regions can help in understanding the phylogeny of Amblyomma ticks present in East and Southeast Asia.

CO N FLI C T O F I NTE R E S T
The authors have no competing interests to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the DNA Data Bank of Japan (http://www.ddbj.nig.ac.jp) with accession numbers of LC554761-LC554790 (samples from Japan) and LC633546-LC633554 (samples from Myanmar).