Release of marketed individuals increases the risk of genetic disturbance in the pet insect Trypoxylus dichotomus

303 Release of marketed individuals increases the risk of genetic disturbance in the pet insect Trypoxylus dichotomus Tomo Hamano1, Yoshihisa Suyama2, Ayumi Matsuo3, Teruaki Ban4,5 , Kohei Watanabe6, Takeshi Yamasaki7,8 , Kazutaka Yamada7,8 , Hiroaki Ishida7,8, Naoyuki Nakahama7,8 1 The Graduate School of Human Science and Environment, University of Hyogo, 1-1-12 Shinzaike-honcho, Himeji, Hyogo 670-0092, Japan 2 Kawatabi Field Science Center, Graduate School of Agricultural Science, Tohoku University, 232-3 Yomogida, Naruko-onsen, Osaki, Miyagi 989-6711, Japan 3 GENODAS Inc., Urbannet Sendai-Chuo Building, 4-4-19 Chuo, Aoba-ku, Sendai, Miyagi 980-0021, Japan 4 Laboratory of Entomology, Obihiro University of Agriculture and Veterinary Medicine, Japan, c/o Tanaka-shinden, Matsudo, 270-2255, Japan; 5 Natural History Museum and Institute, 955-2 Aoba-cho, Chuo-ku, Chiba, Chiba 260-0852, Japan 6 Ishikawa Insect Museum, 3 Inu, Yawata-machi, Hakusan, Ishikawa 920-2113, Japan 7 Institute of Natural and Environmental Sciences, University of Hyogo, 6 Yayoigaoka, Sanda, Hyogo 669-1546, Japan 8 Museum of Nature and Human Activities, Hyogo, 6 Yayoigaoka, Sanda, Hyogo 669-1546, Japan Corresponding author: Tomo Hamano ([email protected]) Copyright: © Tomo Hamano et al. This is an open access article distributed under terms of the Creative Commons Attribution License (Attribution 4.0 International – CC BY 4.0). Research Article Abstract Genetic disturbance can be caused by the release or escape of individuals with different genetic characteristics into wild habitats, risking impacts on native biodiversity. The risk of genetic disturbance in pet insects due to release and escape is particularly common because a wide variety of affordable pets are available on the market. Trypoxylus dichotomus (Coleoptera, Scarabaeidae), the Japanese rhinoceros beetle, is a renowned pet insect in Japan and thus is a suitable target species for studying genetic disturbances in pet insects. However, the detailed spatial genetic structure and genetic disturbances of this species in Japan remain unclear. Here, we estimated the genetic diversity and spatial genetic structure of wild and marketed individuals using mitochondrial DNA sequences and genome-wide single-nucleotide polymorphisms (SNPs) obtained via MIG-seq. Using MIG-seq, 570 SNPs were obtained, revealing a weak yet significant spatial genetic structure in the Japanese archipelago. Although significant isolation by distance (IBD) was observed in wild individuals, no significant IBD was observed in marketed individuals. Comparisons between wild and marketed individuals revealed clear differences in spatial genetic structure. These findings highlight the risks of releasing marketed individuals into the wild owing to their artificial long-distance migration. Our results provide valuable insights into the genetic disturbance of human-mediated distribution and underscore the need for informed management practices to protect native biodiversity. Key words: Coleoptera, genetic distance, MIG-seq, pet insects, release, Scarabaeidae Introduction Genetic disturbance indicates the disruption of the original genetic structure of native wild populations caused by the artificial introduction and breeding of native and artificially introduced species (Rhymer and Simberloff 1996; Frankham 2010; Nakahama et al. 2022). In recent years, this problem has been recognized as a driver of biodiversity loss (Dufresnes et al. 2016; Nakahama et al. 2021). Academic editor: Nathan Havill Received: 23 May 2025 Accepted: 26 August 2025 Published: 3 October 2025 Citation: Hamano T, Suyama Y, Matsuo A, Ban T, Watanabe K, Yamasaki T, Yamada K, Ishida H, Nakahama N (2025) Release of marketed individuals increases the risk of genetic disturbance in the pet insect Trypoxylus dichotomus. NeoBiota 101: 303–320. https://doi.org/10.3897/ neobiota.101.159665 NeoBiota 101: 303–320 (2025) DOI: 10.3897/neobiota.101.159665 Advancing research on alien species and biological invasions A peer-reviewed open-access journal NeoBiota 304 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Genetic disturbances alter the original genetic composition of native wild organisms, sometimes leading to species or population extinction (Nakahama et al. 2021; Kato et al. 2024). Genetic disturbance is mainly caused by the release of individuals raised for purposes such as stocking, aquaculture, pets, or horticulture (Iguchi 2009; Nakao 2017; Nakahama et al. 2021). Human activities have facilitated the artificial introduction of organisms, which extends beyond their natural dispersal capacity and increases the risk of invasion into new habitats (Chiba et al. 2022). Therefore, phylogeographic or population genetic studies are needed to understand actual genetic disturbances and conserve the original genetic structure of wild organisms. Pet insects are at high risk of genetic disturbance because they are easy to purchase, and a wide variety is available on the market (Goka and Kojima 2004; Hosaka et al. 2017). Because of the popularity of rhinoceros and stag beetles in Japan, many have been widely marketed and bred (Goka et al. 2004; Hosaka et al. 2017). These pet insects include many endangered taxa and populations that require conservation. For example, the Okinawan rhinoceros beetle, Trypoxylus dichotomus takarai, and the stag beetle Dorcus hopei binodulosus have been designated as Data Deficient (DD) and Vulnerable (VU), respectively, by the Red List (Ministry of the Environment Japan 2020), owing to excessive collection, hybridization among subspecies, and habitat loss. However, studies exploring pet insect-induced genetic disturbances in Japan remain lacking. T. dichotomus (Coleoptera, Scarabaeidae) is widely distributed in East Asia, including China, Japan, Korea, Vietnam, Myanmar, Laos, India, and Thailand (Nagai 2006, 2007; Satoru 2014; Adachi 2017). In Japan, this species mainly inhabits forests at altitudes below 1,500 m (Unno 2006), particularly in Satoyama landscapes – traditional rural areas shaped by long-standing interactions between people and nature that support rich biodiversity (Takeuchi et al. 2016). Within Japan, five subspecies of the rhinoceros beetle are distributed as follows: the Japanese mainland subspecies T. d. septentrionalis (Kono 1931), the Okinawa subspecies T. d. takarai (Kusui 1976), the Kumejima subspecies T. d. inchachina (Kusui 1976), the Kuchinoerabujima subspecies T. d. tsuchiyai (Nagai 2006), and the Yakushima and Tanegashima subspecies T. d. shizuae (Adachi 2017). T. dichotomus is a popular pet owing to its large size and distinctive horns; however, genetic disturbance caused by introducing them artificially or through escape into non-native habitats has become a concern (Kusui 1976; Muranaka and Ishihama 2010; Yang et al. 2021). In fact, on Hokkaido Island, Japan, where T. dichotomus was not originally distributed, introduced individuals were recorded in 1936 because of human activities (Kida 2003). A few individuals with different genetic characteristics from their original habitats have also been reported on the Korean Peninsula and the Japanese Goto Islands (Yang et al. 2021; Hamano et al. 2024). The positive correlation between growth rate during the larval stage and latitude suggests that genetic disturbances from this species negatively impact wild populations (Kojima et al. 2020). No population genetic analysis of T. dichotomus has been conducted in Japan, and the current status of genetic disturbances remains unknown. Although phylogeographic studies of T. dichotomus in East Asia have been conducted, the sample size was insufficient to estimate the population genetic structure within Japan (Yang et al. 2021). Therefore, this study examined the genetic disturbance of T. dichotomus in Japan by conducting population genetic analysis of wild and marketed individuals 305 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus across Japan. We conducted this study using genome-wide SNP analysis (multiplexed ISSR [inter-simple sequence repeat] genotyping by sequencing [MIG-seq]) and mitochondrial DNA sequences. This study further assessed the risk of genetic disturbance from the release and escape of marketed individuals by comparing their spatial genetic structures. Materials and methods Collection of study samples 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). Phylogenetic analysis based on mitochondrial DNA 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 306 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus 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) (b) 307 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus 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. MIG-seq analysis 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 pairedend 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. Population genetic structure 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 308 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus 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. Genetic diversity 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. Results Phylogenetic analysis based on mitochondrial DNA 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%). 309 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus 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. Fukuoka_W67 Kochi_W71 (N=3) Nigata_W31 Akita_W19_4 Yamagata_W25 Kyoto_W51 Aomori_W18_3 Fukuoka_W66_4 Kyoto_W51_2 Kochi_W71 Nigata_W31 Yamagata_W29 Kochi_W71_3 Fukui_W35_8 Miayazaki_W77_2 Tokushima_W63 (N=2) Kanagawa_W45 Osaka_W61_2 Hokkaido_W02, W11 Akita_W19 Tottori_W44 Fukui_W35_30 Osaka_W57 Nigata_W31_2 Hokkaido_W03, W05, W06, W09, W10 Hokkaido_W08 (N=3) Fukui_W35_31 Aomori_W18 Fukui_W35 (N=4) Tottori_W47 Chiba_W48 (N=2) Nagasaki_W72 Kanagawa_W45_2 Iwate_W20 Yamagata_W22 Chiba_W41 Kagoshima_W81 Nigata_W28 Nigata_W28_5 Nigata_W28 (N=2) Yamagata_W29_4 Kochi_W71_6 Chiba_W42 Osaka_W74 Hyogo_W56 Nagsaki_W47_6, W74_5, W75 Nagasaki_W7 (N=2) Kagoshima_W80 Hokkaido_W02 Gunma_W34 Nagasaki_W73_3 Nagasaki_W74 (N=4) Nagsaki_W74_4 Miyagi_W27 (N=4) Miyagi_W26 (N=2) Chiba_W37_4 Miyagi_W26 (N=3) Chiba_W50 Tottori_W44_5 Hokkaido_W16 Fukui_W35_29 Hyogo_W56_5 Fukui_W35 (N=2) Osaka_W58 (N=2) Osaka_W53 Saitama_W36_2 Miyagi_W27_2 Hyogo_W56_3 Miyazaki_W77 Saitama_W36 Kyoto_W51_3 Kochi_W71_5 Kochi_W65 Osaka_W61 Tokyo_W40 Nagano_W33 Hokkaido_W14 Hokkaido_W07 Kochi_W70 Fukuoka_W66 (N=2) Hokkaido_W11_2 Hyogo_W52 Nagasaki_W72_2 Nagasaki_W66 Chiba_W39 Hokkaido_W12 Yamagata_W24 Yamanashi_W43 Kochi_W68 Hokkaido_W01 (N=2) Nagano_W38 Kagoshima_W79 Miyazaki_W78 Miyagi_W27_4 Fukuoka_W67_2 Fukushima_W32_3 Tottori_W47_2 Nigata_W28 (N=4) Akita_W19 (N=2) Yamaguchi_W60 Miyazaki_W76 Hokkaido_W02 (N=2) Nigata_W31_3 Nagasaki_W75_2 Miyagi_W26 (N=3) Aomori_W1 (N=2) Hokkaido_W15 (N=2) Nigata_W30 Hokkaido_W05 (N=2) Tottori_W44 (N=4) Yamagata_W29 (N=4) Hyogo_W56 (N=3) Osaka_W59 (N=2) Yamagata_W23 Wakayama_W62 Hyogo_W55 Hokkaido_W08_2 Hokkaido_W16 (N=2) Kagoshima_W82 Chiba_W37 (N=3) Chiba_W49 Hokkaido_W13 Osaka_W58 (N=3) Hokkaido_W04 (N=2) Fukushima_W32 (N=2) Kagoshima_W80 (N=3) Kochi_W69 (N=4) Tokushima_W63 (N=6) Ehime_W64 (N=6) Hokkaido_W14 (N=6) Fukui_W35 (N=22) Akita_W21 Nagasaki_W73 (N=5) Okinawa_W84 Okinawa_W83 (N=4) 0.005 100 100 Japanese mainland Isl. clade Okinawa and Kumejima Isl. clade 310 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Demographic analysis based on mitochondrial DNA 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). Spatial genetic structure based on MIG-seq 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). 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). Genetic diversity based on MIG-seq 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 311 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Figure 3. Pie charts showing cluster distributions of wild and marketed populations by region, based on SNP analysis and STRUCTURE analysis using MIG-seq. Wild Marketed Cluster. 5 Cluster. 2 Cluster. 1 Cluster. 4 Cluster. 3 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. 0 200 400 600 800 1000 1200 0.0 500.0 1000.0 1500.0 2000.0 Genetic distance Co - dominant genotypic distance Geographic distance (km) (a) Wild individuals Mantel test, P = 0.002 0 200 400 600 800 1000 1200 0.0 500.0 1000.0 1500.0 2000.0 2500.0 Genetic distance Co - dominant genotypic distance Geographic distance (km) (b) Marketed individuals Mantel test, P = 0.071 318 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Nagai S (2006) A new species and new subspecies of the genus Trypoxylus from Asia and a new subspecies of the genus Beckium from New Guinea (Coleoptera, Scarabaeidae, Dynastinae). Gekkan-Mushi 428: 13–17. Nagai S (2007) Nihon no Kabutomushi Daizukan. BE-KUWA 22: 8–29. Nakahama N, Isagi Y (2018) Recent transitions in genetic diversity and structure in the endangered semi-natural grassland butterfly Melitaea protomedia in Japan. Insect Conservation and Diversity 11: 330–340. https://doi.org/10.1111/icad.12280 Nakahama N, Asai T, Matsumoto S, Suetsugu K, Kurashima O, Matsuo A, Suyama Y (2021) Detection and dispersal risk of genetically disturbed individuals in the endangered wetland plant Pecteilis radiata (Orchidaceae) in Japan. Biodiversity and Conservation 30: 1913–1927. https:// doi.org/10.1007/s10531-021-02174-y Nakahama N, Hanaoka T, Itoh T, Kishimoto T, Ohwaki A, Matsuo A, Kitahara M, Usami S, Suyama Y, Suka T (2022) Identification of source populations for reintroduction in extinct populations based on genome-wide SNPs and mtDNA sequence: A case study of the endangered subalpine grassland butterfly Aporia hippia (Lepidoptera, Pieridae) in Japan. Journal of Insect Conservation 26: 121–130. https://doi.org/10.1007/s10841-022-00369-4 Nakao R (2017) Current status of genetic disturbance in wild medaka (Oryzias latipes species complex) in Japan. Nippon Suisan Gakkaishi 83: 235–235. https://doi.org/10.2331/suisan. WA2353-4 Nei M (1972) Genetic distance between populations. The American Naturalist 106: 283–292. https://doi.org/10.1086/282771 Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89: 583–590. https://doi.org/10.1093/genetics/89.3.583 Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ (2015) IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution 32: 268–274. https://doi.org/10.1093/molbev/msu300 Peakall R, Smouse PE (2006) GenAlEx 6: Genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes 6: 288–295. https://doi.org/10.1111/j.14718286.2005.01155.x Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959. https://doi.org/10.1093/genetics/155.2.945 Rambaut A (2018) FigTree v1.4.4. https://tree.bio.ed.ac.uk/software/figtree/ Rhymer JM, Simberloff D (1996) Extinction by hybridization and introgression. Annual Review of Ecology and Systematics 27: 83–109. https://doi.org/10.1146/annurev.ecolsys.27.1.83 Roweis S (1998) EM algorithms for PCA and SPCA. Advances in Neural Information Processing Systems 10: 626–632. Sato Y, Ishikawa R (2004) The world of plants from the Sannai-Maruyama site: a DNA archaeological perspective. SHOKABO, Tokyo. [in Japanese] Satoru T (2014) A new subspecies of Trypoxylus dichotomus (Coleoptera, Scarabaeidae, Dynastinae) from China. Smouse PE, Peakall ROD (1999) Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity 82: 561–573. https://doi.org/10.1038/sj.hdy.6885180 Suetsugu K, Nozaki T, Hirota S, Funaki S, Ito K, Isagi Y, Suyama Y, Kaneko S (2023) Phylogeographical evidence for historical long-distance dispersal in the flightless stick insect Ramulus mikado. Proceedings of the Royal Society B, Biological Sciences 290: 20231708. https://doi.org/10.1098/ rspb.2023.1708 Suyama Y, Matsuki Y (2015) MIG-seq: An effective PCR-based method for genome-wide single-nucleotide polymorphism genotyping using the next-generation sequencing platform. Scientific Reports 5: 16963. https://doi.org/10.1038/srep16963 319 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Suyama Y, Hirota S, Matsuo A, Tsunamoto Y, Mitsuyuki C, Shimura A, Okano K (2022) Complementary combination of multiplex high-throughput DNA sequencing for molecular phylogeny. Ecological Research 37: 171–181. https://doi.org/10.1111/1440-1703.12270 Takeuchi K, Ichikawa K, Elmqvist T (2016) Satoyama landscape as social–ecological system: Historical changes and future perspective. Current Opinion in Environmental Sustainability 19: 30–39. https://doi.org/10.1016/j.cosust.2015.11.001 Tanaka Y (2017) Notes on non-native stag beetle species (Coleoptera, Lucanidae) observed in Itami City, Hyogo Prefecture, Japan. Itakon 5: 31–33. Thompson JD, Higgins DG, Gibson TJ (1994) CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Research 22: 4673–4680. https://doi. org/10.1093/nar/22.22.4673 Uchifune T (2012) Subsequent report on movement of a Japanese rhinoceros beetle (Coleoptera, Scarabaeidae) in the Miura Peninsula in 2011. Scientific Reports of Yokosuka City Museum 59: 31–32. Waku D, Segawa T, Yonezawa T, Akiyoshi A, Ishige T, Ueda M, Ogawa H, Sasaki H, Ando M, Kohno N, Sasaki T (2016) Evaluating the phylogenetic status of the extinct Japanese otter on the basis of mitochondrial genome analysis. PLOS ONE 11: e0149341. ttps://doi.org/10.1371/ journal.pone.0149341 Weber JN, Kojima W, Boisseau RP, Niimi T, Morita S, Shigenobu S, Gotoh H, Araya K, Lin C-P, Thomas-Bulle C, Allen CE, Tong W, Lavine LC, Swanson BO,. Emlen DJ (2023) Evolution of horn length and lifting strength in the Japanese rhinoceros beetle Trypoxylus dichotomus. Current Biology 33: 4285–4297. https://doi.org/10.1016/j.cub.2023.08.066 Yang H, You CJ, Tsui CKM, Tembrock LR, Wu ZQ, Yang DP (2021) Phylogeny and biogeography of the Japanese rhinoceros beetle, Trypoxylus dichotomus (Coleoptera: Scarabaeidae) based on SNP markers. Ecology and Evolution 11: 153–173. https://doi.org/10.1002/ece3.6982 Yoshitake H, Hosoya T, Yamada R (2016) Rhinoceros beetles collected aboard the ferry ‘Toshima’. Sayabane, new series 23: 47. Supplementary material 1 Sampling information and genetic diversity indices of wild and marketed populations of Trypoxylus dichotomus Authors: Tomo Hamano, Yoshihisa Suyama, Ayumi Matsuo, Teruaki Ban, Kohei Watanabe, Takeshi Yamasaki, Kazutaka Yamada, Hiroaki Ishida, Naoyuki Nakahama Data type: xlsx Explanation note: The file provides detailed sampling information and genetic diversity indices (e.g., nucleotide diversity, heterozygosity) for wild and marketed populations of Trypoxylus dichotomus analyzed in this study. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.101.159665.suppl1 320 NeoBiota 101: 303–320 (2025), DOI: 10.3897/neobiota.101.159665 Tomo Hamano et al.: Genetic disturbance of T. dichotomus Supplementary material 2 Pairwise FST values among wild and marketed populations of Trypoxylus dichotomus Authors: Tomo Hamano, Yoshihisa Suyama, Ayumi Matsuo, Teruaki Ban, Kohei Watanabe, Takeshi Yamasaki, Kazutaka Yamada, Hiroaki Ishida, Naoyuki Nakahama Data type: xlsx Explanation note: This file provides pairwise FST values among wild and marketed populations of Trypoxylus dichotomus, showing the degree of genetic differentiation between populations. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.101.159665.suppl2 Supplementary material 3 Optimal number of clusters (K) in STRUCTURE analysis Authors: Tomo Hamano, Yoshihisa Suyama, Ayumi Matsuo, Teruaki Ban, Kohei Watanabe, Takeshi Yamasaki, Kazutaka Yamada, Hiroaki Ishida, Naoyuki Nakahama Data type: docx Explanation note: Optimal number of clusters (K) in STRUCTURE analysis; Mismatch distribution for wild samples based on 685 bp of mitochondrial COII gene; Principal component analysis (PCA) plots based on 570 SNPs obtained via MIG-seq; Results of STRUCTURE analysis for K values ranging from 2 to 5. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/neobiota.101.159665.suppl3