Between 2016 and 2022, 258 wild individuals were collected from 84 locations across the Japanese archipelago (Fig. 1; Suppl. material 1). During the sampling survey, 63 individuals sold in local markets (hereafter marketed individuals) were collected from 26 locations (Fig. 1; Suppl. material 1). The legs of the collected individuals were preserved in 99% ethanol or 99% propylene glycol. DNA was extracted from each sample following the standard protocol of the Wizard Genomic DNA Purification Kit (Promega Corporation, USA).

Figure 1.

Sampling sites per population of T. dichotomus in Japan. Red dots indicate wild populations (W), and blue dots indicate marketed populations (M).

A molecular phylogenetic analysis of T. dichotomus was performed using the mitochondrial DNA COII region, following the PCR conditions described by Yang et al. (2021). The COII region was amplified using the primers F-lue (5´-TCTAATATGGCAGATTAGTGC-3´) and R-lys (5´-GAGACCAGTACTTGCTTTCAGTCATC-3´). The PCR reaction was prepared with 1 µL of DNA sample, 0.1 µL of 20 µM forward primer, 0.1 µL of 20 µM reverse primer, 2.0 µL of 2 mM dNTP, 5.0 µL of KOD FX buffer, and 0.2 µL of KOD FX Neo (TOYOBO, Osaka, Japan), with a final volume of 10 µL. The primers used in the PCR were applied for bidirectional sequencing at Eurofins Genomics (Tokyo, Japan). The collected sequences were registered in GenBank (accession numbers: PX240152–PX240389).

Only wild individuals were used for the molecular phylogenetic analysis. Alignment was performed using CLUSTAL W (Thompson et al. 1994) implemented in MEGA X ver. 10.2 (Kumar et al. 2018). A 685 bp COII region was obtained from all samples. Phylogenetic relationships were estimated using the maximum likelihood method with 1,000 ultrafast bootstrap (UFBoot) replicates. HKY+F was selected as the optimal substitution model based on the Bayesian information criterion (BIC) using ModelFinder in IQ-TREE ver. 2.2.0 (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017). The estimated phylogenetic tree was constructed using FigTree (Rambaut 2018).

Mismatch distribution analysis was performed using mitochondrial DNA (COII region) sequences of T. dichotomus to evaluate the historical demographic dynamics of the populations. The analysis was performed using Arlequin ver. 3.5.2.2 (Excoffier and Lischer 2010). The sudden population expansion model was applied, and the fit between the observed data and the model was assessed using default settings in Arlequin.

The statistical metrics included Tajima’s D and Fu’s Fs. The analysis was performed on two datasets: one including all individuals and the other excluding individuals from Hokkaido and Okinawa Islands. Individuals on Hokkaido Island were excluded because this species is not native there and because the populations were introduced by humans. Individuals on the Okinawa Islands were excluded because they are considered subspecies phylogenetically distinct from the main population of this species.

The multiplexed ISSR genotyping by sequencing (MIG-seq) method was used to detect SNPs in extracted DNA samples (Suyama and Matsuki 2015; Suyama et al. 2022). MIG-seq is a genome-wide SNP analysis approach that amplifies ISSR regions without restriction enzymes. The more than 100 genome-wide SNPs obtained using this method can be applied to studies of genetic differentiation between closely related species and to phylogeographic studies among populations and species (Suyama et al. 2022). Therefore, this method is considered suitable for phylogeographic studies of T. dichotomus. The analytical procedure followed the standard experimental protocol described by Suyama et al. (2022). For the first PCR of genomic DNA, a pair of multiplex ISSR primers (MIG-seq primer set 1) specifically developed as MIG-seq primers was used. The first PCR product was used as the template for the second PCR. Each sample was independently constructed, and an indexed library was created for the Illumina sequencing platform. One µL of each second PCR product was pooled into a single library, purified, and size-selected (300–800 bp) using the Pippin Prep DNA size selection system (Sage Science, Beverly, MA, USA). The size-selected library was sequenced on an Illumina MiSeq system (Illumina, San Diego, CA, USA) using the MiSeq reagent kit ver. 3 (150 cycles, Illumina). Quality filtering of raw sequence data from MIG-seq was performed using Trimmomatic ver. 0.36 (Bolger et al. 2014). The filtered paired-end reads were assembled into contigs based on sequence similarity using the gstacks module in Stacks ver. 2.62. Among 43,162 loci, 99.5% (42,963 loci) were successfully assembled, with a mean contig length of 160.8 bp. The average insert size was 122.8 bp (SD = 35.5), and 98.6% of paired-end reads were successfully aligned (gstacks.log). SNPs shared by more than 10% of all samples (R = 0.6), with a minor allele frequency below 0.05 (min-maf = 0.05) or a heterozygosity above 0.6 (max-obs-het = 0.6), were removed using Stacks ver. 2.65 (Catchen et al. 2011, 2013). Individuals with loci missing over 40% were excluded using TASSEL ver. 5 (Bradbury et al. 2007). After filtering the raw read data, 570 SNPs were detected in 277 samples in the MIG-seq dataset. Raw sequence data of MIG-seq have been deposited in the NCBI Sequence Read Archive under the BioProject accession PRJNA1311727 and will be released upon publication.

To assess spatial genetic structure, we conducted a principal component analysis (PCA) using the 570 SNPs obtained via MIG-seq. PCA was performed using TASSEL ver. 5 (Roweis 1998; Bradbury et al. 2007), and the results were visualized in R ver. 4.2.2 using the ggplot2 package. Two separate analyses were performed: one including all individuals and another excluding individuals from the Okinawa Islands (including Kumejima Island), which are phylogenetically distinct from the Japanese mainland populations.

We further analyzed population structure using STRUCTURE ver. 2.3.4 (Pritchard et al. 2000). Samples were collected from across the Japanese archipelago, including Hokkaido, the Okinawa Islands, and the Goto Islands, such as Yakushima, Tanegashima, and Kumejima. We applied a model-based clustering approach using allele frequencies and an admixture model. STRUCTURE was run with K values ranging from 1 to 10, each replicated 10 times with a burn-in of 10,000 steps and 100,000 MCMC iterations. The optimal K value was determined using STRUCTURE Harvester (Earl and vonHoldt 2012) based on the average log-likelihood and ΔK. The ΔK values were also evaluated using the K algorithm in Structure Selector (Li and Liu 2018), which consistently supported K = 5 as the most suitable number of clusters.

To evaluate the relationship between genetic and geographic distances, we conducted a Mantel test using GenAlex ver. 6.41 (Mantel 1967; Peakall and Smouse 2006). Wild and marketed individuals were analyzed separately, and pairwise codominant genotypic distances were calculated (Smouse and Peakall 1999). The Mantel test was run with 9,999 permutations to assess the correlation between genetic differentiation (Nei’s distance; Nei 1972, 1978) and geographic distance. For this analysis, individuals from the Okinawa Islands (including Kumejima) and Hokkaido were excluded, as they represent genetically distinct or introduced populations.

The genetic diversity indices, including observed heterozygosity (Ho) and allele frequencies, were calculated using the Hierfstat package (Goudet 2005) in the R platform (ver. 3.2.3). We used the glmmTMB package in R (ver. 1.2) to construct a generalized linear mixed model to analyze the relationship between individual heterozygosity and sample characteristics. The response variable was the observed heterozygosity of each individual, calculated as the proportion of heterozygous loci to the total number of loci, and modeled assuming a binomial distribution. The geographic location of each sample, specifically latitude or longitude, was used as an explanatory variable and treated as continuous. To account for the non-independence of samples within the same population, population information was included as a random effect with no specific covariance structure specified. Model fit was assessed using the Akaike information criterion. Individuals from the Okinawa Islands (including Kumejima Island) were excluded from the analysis because they belong to a different subspecies and are phylogenetically distinct from individuals in other Japanese islands (Yang et al. 2021). We also excluded individuals from Hokkaido Island, as they represent introduced populations and are not considered part of the native range.

Phylogenetic analysis revealed that T. dichotomus populations across the Japanese archipelago were divided into two major clades (Fig. 2). The first clade (Japanese mainland Isl. clade) included populations ranging from Hokkaido to Kyushu Islands and their neighboring smaller islands (Sadogashima, Yakushima, Tanegashima, and Fukuejima) (Fig. 2). Only a few obvious subclades were observed in this clade. The second clade (Okinawa and Kumejima Isl. clade) included individuals from Okinawa Island (W83) and Kumejima Island (W84). This clade was clearly different from that of the Japanese mainland populations (UFBoot: 100%).

Figure 2.

Phylogenetic tree of Japanese rhinoceros beetles generated using neighbor-joining analysis based on the maximum likelihood method. Numbers on the branches indicate support values from 1,000 ultrafast bootstrap replicates. Each sample in the tree is labeled with its location name, population number, and individual identification number.

Fu’s Fs values for the populations in the Japanese mainland clade were significantly <0 (Fs = -99.181, P < 0.001). Tajima’s D for these populations was also significantly different from 0 (D = -2.644, P < 0.001). Mismatch analysis revealed no significant deviation from the expectation of the null model, which assumed a population expansion of the Japanese mainland (SSD = 0.00026, P = 0.728; Suppl. material 3: fig S2).

Principal component analysis (PCA) based on 570 SNPs obtained via MIG-seq revealed patterns consistent with the mitochondrial COII phylogenetic analysis. Individuals from the Okinawa and Kumejima Islands formed a distinct cluster in PCA space (Suppl. material 3: fig S3a), corresponding to the clade separation observed in the COII phylogeny (Fig. 2). In contrast, when these individuals were excluded, no distinct clustering was observed among the remaining populations from the Japanese mainland (Suppl. material 3: fig S3b), indicating weak or absent genetic structure. STRUCTURE analysis, conducted using the same SNP dataset for 279 individuals (216 wild and 63 marketed), identified the optimal number of genetic clusters as K = 5 based on the ΔK criterion (Fig. 3; Suppl. material 3: fig S1). For wild populations, weak spatial genetic structure was observed across the Japanese mainland, with Cluster 1 dominant in Hokkaido and northern Honshu, Cluster 2 in western Honshu, Shikoku, and Kyushu, and Cluster 3 primarily in Yakushima and Tanegashima (Fig. 3; Suppl. material 3: fig S4). The Goto Islands showed a distinct cluster (Cluster 4), while Okinawa and Kumejima were represented by Cluster 5 and a mix of Clusters 3 and 5, respectively. Marketed populations exhibited a markedly different structure, with Cluster 1 widely distributed across Japan, including regions where it was less frequent in wild populations, and Clusters 2 and 3 prevalent in the Chugoku and northern Tohoku regions, respectively. To evaluate isolation by distance, a Mantel test was performed for wild and marketed populations separately. In wild populations from the Japanese mainland (excluding Hokkaido, Okinawa, and Kumejima Island), a significant positive correlation was found between geographic and genetic distances (Rxy = 0.237, P < 0.01; Fig. 4). In contrast, no significant correlation was detected in marketed populations (Rxy = 0.016, P = 0.071).

Figure 3.

Pie charts showing cluster distributions of wild and marketed populations by region, based on SNP analysis and STRUCTURE analysis using MIG-seq.

Figure 4.

Relationships between codominant genotypic distances based on data from 570 single-nucleotide polymorphisms and geographic distances at the individual level. a. Gray circles represent relationships based on all wild populations; b. Filled circles represent relationships based on all marketed populations.

In the Japanese mainland (excluding Hokkaido, Okinawa, and Kumejima Islands), the Mantel test revealed a significant positive correlation between geographic and genetic distances in wild populations (Rxy = 0.237, P < 0.01; Fig. 4). Conversely, no significant correlation was observed in marketed individuals (Rxy = 0.016, P = 0.071).

We investigated spatial patterns of genetic diversity using observed heterozygosity values derived from SNPs obtained via MIG-seq. First, we analyzed all wild individuals across Japan (Fig. 5a, b). No significant correlation was found between heterozygosity and latitude (correlation = -0.994, P = 0.1538; Fig. 5a), whereas a significant positive correlation was observed with longitude (correlation = -0.999, P < 0.05; Fig. 5b). To assess whether these patterns were influenced by geographically distinct populations, we reanalyzed the data excluding individuals from the Okinawa and Hokkaido populations (Fig. 5c, d). In this analysis, both latitude (correlation = -0.998, P < 0.05; Fig. 5c) and longitude (correlation = -1.000, P < 0.001; Fig. 5d) showed significant positive correlations with heterozygosity. The model selection results (Suppl. material 2) also supported these findings. The models including longitude as an explanatory variable yielded lower AIC values than the null model, indicating better explanatory power. The dotted lines in each panel represent regression lines estimated using generalized linear models.

Figure 5.

Geographic locations where T. dichotomus individuals were collected in relation to genetic diversity based on 570 single-nucleotide polymorphisms. The dotted lines represent regression lines determined by generalized linear models. a, b. Results for all individuals; c, d. Results excluding Okinawa and Hokkaido populations.

These results suggest that the observed heterozygosity in the studied populations is significantly influenced by latitude and longitude when potential outlier populations, such as those from the Okinawa Islands and Hokkaido Island, are excluded.

The results of the phylogenetic tree based on the COII region indicated a lack of clear genetic differentiation in Japan, except for Okinawa and the Kumejima Islands (Fig. 3). This result supports previous studies (Yang et al. 2021; Weber et al. 2023). Genome-wide SNP analysis using the MIG-seq method showed that weak genetic clines were observed in the Japanese mainland from Hokkaido to the Yakushima Islands. The weak genetic cline observed in the Japanese mainland is likely due to the high dispersal capacity of T. dichotomus. Notably, adult individuals of this species can fly an average distance of approximately 1.3 km per generation, with a maximum recorded distance of 3.6 km (Uchifune 2012). Such high mobility allows T. dichotomus to rapidly colonize new habitats, contributing to its widespread and continuous distribution across Japan.

A generalized linear mixed model analysis revealed higher genetic diversity in the northern and eastern parts of mainland Japan, except for Hokkaido, Okinawa, and Kumejima Islands. The exact factors underlying this higher genetic diversity are unknown. However, possible factors include suitable landscapes for T. dichotomus habitats and genetic disturbances. Historical environmental changes may also influence genetic diversity. For instance, the formation of Satoyama landscapes may have favored species adapted to human-modified environments, allowing them to expand their distribution and maintain genetic diversity. Furthermore, the finding that individuals from Yakushima and Tanegashima Islands belong to the same clade as those from the Japanese mainland suggests that dispersal events to these islands occurred much more recently than to Okinawa and Kumejima Island populations (Weber et al. 2023). Consequently, the genetic differentiation between the Yakushima and Tanegashima Island populations and the Japanese mainland populations remains relatively low. Integrating paleoenvironmental data and historical records is essential to better understand the effects of past environmental changes. Factors related to genetic disturbances are discussed in the following section.

The spatial genetic structure of the marketed and wild populations differed significantly. Isolation by distance was observed in the wild population but not in the marketed population. This is an example of genetic analysis showing that long-distance anthropogenic movement occurred not only in neighboring areas but also in individuals distributed as pets. In fact, there have been reported cases of T. dichotomus being transported long distances, and even on the main island of Okinawa, individuals from mainland Japan are sold (Hamano T., personal communication). This suggests that genetic disturbances can become more serious if sold T. dichotomus are released into the habitat and interbreed with wild beetles.

The STRUCTURE analysis also identified a unique genetic cluster (Cluster 1) in northeastern Japan, distinct from other regions. In general, it is unlikely that insect species with limited dispersal capacity share identical genetic clusters across remote regions unless assisted by anthropogenic movement (Suetsugu et al. 2023). These findings suggest that the release of marketed individuals into the wild may have already taken place, particularly in regions characterized by high genetic diversity.

Several potential ecological and evolutionary consequences may result from such releases. First, historically and geographically unique genetic lineages may be lost. For example, populations from the Goto Islands (e.g., Fukuejima Island) harbor distinctive haplotypes that represent important evolutionary heritage (Hamano et al. 2024). Hybridization with released individuals could lead to their irreversible loss.

Second, the detection of genetically disturbed individuals is extremely challenging, particularly when hybridization occurs within a species with weak reproductive isolation between subspecies or populations (Muranaka and Ishihama 2010). Such cases are difficult to distinguish morphologically, as has been documented in other taxa (Kawamura 2023). Given the high abundance and dispersal capacity of T. dichotomus, genetic disturbance may rapidly expand across wide areas once hybridization begins.

Third, disruption of locally adapted phenology is another major concern. Insects like T. dichotomus often exhibit flexible phenological traits and can adjust to a range of environmental conditions – as seen in the established population in Hokkaido (Kojima et al. 2020). However, if non-native individuals alter the timing of development and activity in local populations, this could erode regionally adapted phenological patterns that have evolved to match local climate conditions. While such flexibility may be beneficial in the short term, it could reduce the population’s ability to respond to future environmental changes, such as abnormal weather or long-term climate shifts, ultimately increasing their vulnerability.

Taken together, these risks highlight the importance of early detection, continuous monitoring, and careful management to prevent irreversible changes to the genetic structure and ecological functioning of native populations.

The increase in Satoyama landscapes during the Holocene may have facilitated the population expansion of T. dichotomus (Sato and Ishikawa 2004). These human-managed secondary habitats likely provided suitable environments for T. dichotomus, supporting population growth and expansion. However, we were unable to estimate the period of change in effective population sizes using demographic analyses such as Bayesian skyline plots or Stairway plots. This is because inaccurate results may be obtained if genetically disturbed individuals are included in the analyses. Despite the high likelihood that the studied wild populations of T. dichotomus contained genetically disturbed individuals, they could not be identified. In addition, human-mediated movements have been observed in other insect species, often resulting in competition with native populations and the potential introduction of new pathogens, which could further complicate conservation efforts (Ishibashi 2007).

This study revealed genetic disturbance–associated risks in pet insects and other commercially valuable species. For example, D. titanus and D. hopei (Lucanidae) exhibited genetic disturbances caused by domestic introductions and hybridization with subspecies (Hosoya and Araya 2010; Tanaka 2017), highlighting anthropogenic impacts on native biodiversity. In the case of T. dichotomus, the differences in genetic composition between wild and marketed populations suggest that large-scale transportation and release have already altered natural genetic patterns. Thus, implementing genetic monitoring and conservation measures is crucial for mitigating these risks. Regular genetic assessments can help detect and manage genetic disturbances before they affect native populations. Public education campaigns are also essential for raising awareness of the ecological consequences of releasing marketed individuals into the wild, and providing clear guidelines on the ethical handling and disposal of pet insects can reduce the likelihood of genetic disturbance.

One limitation of this study is the inability to identify specific genetically disturbed individuals or regions. However, prior studies (Hamano et al. 2024; Yoshitake et al. 2016) suggest that such disturbances are likely to exist. Future studies should use historical specimens to clarify the extent and causes of these genetic disturbances. Historical specimens, such as those used in studies by Waku et al. (2016) and Nakahama and Isagi (2018), provide insights into past gene flow patterns, evolutionary history, and the influence of human activities on species distribution. Furthermore, integrating genetic data with ecological and behavioral studies could provide a more comprehensive understanding of how T. dichotomus interacts with its environment. For example, T. dichotomus larvae are known to grow at different rates depending on latitude (Kojima et al. 2020), suggesting potential ecological consequences of genetic disturbance. These interactions should be explored to inform conservation strategies that address both genetic and ecological dimensions.

In conclusion, this study highlights the differences in the spatial genetic structure between the phylogeographic and marketed populations of T. dichotomus in Japan. In the future, genetic monitoring combined with social education and conservation measures can reduce the risk of genetic disturbance and ensure the long-term survival of this species in its natural habitats.

The authors have declared that no competing interests exist.

No ethical statement was reported.

No use of AI was reported.

This research was supported by the SPRING Research Grant from the University of Hyogo Career Development Program for Young Scientists, a Grant-in-Aid for Scientific Research (19K15856 and 23K13969) from the Japan Society for the Promotion of Science, the 2024 Scholarship from the Idea Environmental and Art Foundation, the 2024 Research Grant from the Hyogo Biological Society, and the 2024 Research Grant from the Kansai Organization for Nature Conservation (KONC).

Tomo Hamano conceptualized and designed the study, conducted field sampling, performed DNA extraction and data analysis, prepared figures and tables, and wrote and revised the manuscript. Yoshihisa Suyama and Ayumi Matsuo performed MIG-seq experiments and contributed to data generation (Investigation). Teruaki Ban and Kohei Watanabe assisted with field sampling. Takeshi Yamasaki and Kazutaka Yamada provided advice on methodology and manuscript preparation. Hiroaki Ishida supervised the research. Naoyuki Nakahama provided training in research methodology, contributed to funding acquisition, manuscript revision, and overall supervision.

Tomo Hamano https://orcid.org/0009-0003-8371-7348

Yoshihisa Suyama https://orcid.org/0000-0002-3136-5489

Teruaki Ban https://orcid.org/0009-0009-9470-3869

Kohei Watanabe https://orcid.org/0000-0002-8761-232X

Takeshi Yamasaki https://orcid.org/0000-0002-2419-188X

Kazutaka Yamada https://orcid.org/0000-0002-4210-6693

Naoyuki Nakahama https://orcid.org/0000-0003-3106-8289

All of the data that support the findings of this study are available in the main text or Supplementary Information.