Research articlesPhylogeographical evidence for historical long-distance dispersal in the flightless stick insect Ramulus mikado
Abstract
Exploring how organisms overcome geographical barriers to dispersal is a fundamental question in biology. Passive long-distance dispersal events, although infrequent and unpredictable, have a considerable impact on species range expansions. Despite limited active dispersal capabilities, many stick insect species have vast geographical ranges, indicating that passive long-distance dispersal is vital for their distribution. A potential mode of passive dispersal in stick insects is via the egg stage within avian digestive tracts, as suggested by experimental evidence. However, detecting such events under natural conditions is challenging due to their rarity. Therefore, to indirectly assess the potential of historical avian-mediated dispersal, we examined the population genetic structure of the flightless stick insect Ramulus mikado across Japan, based on a multifaceted molecular approach [cytochrome oxidase subunit I (COI) haplotypes, nuclear simple sequence repeat markers and genome-wide single nucleotide polymorphisms]. Subsequently, we identified unique phylogeographic patterns, including the discovery of identical COI genotypes spanning considerable distances, which substantiates the notion of passive long-distance genotypic dispersal. Overall, all the molecular data revealed the low and mostly non-significant genetic differentiation among populations, with identical or very similar genotypes across distant populations. We propose that long-distance dispersal facilitated by birds is the plausible explanation for the unique phylogeographic pattern observed in this flightless stick insect.
1. Introduction
The segregation of populations by physical barriers and their dispersal across such obstacles constitute two prominent antagonistic forces shaping the distribution and speciation of organisms [1]. The development of wings in insects is widely regarded as a significant contributor to their prosperity and diversity [2,3], as it has facilitated predator evasion, prey capture and migration [2]. Nevertheless, the loss or reduction of wings in various insect lineages is well documented [4], exerting a profound influence on the biogeographical and speciation patterns of these lineages [5,6]. The loss of flight within a species can restrict dispersal capabilities and foster genetic differentiation among populations, potentially leading to an elevated speciation rate in flightless lineages [5], while the evolution of flight in insects played a crucial role in their early diversification [7].
Consequently, flightless insects offer intriguing models for examining population genetic structure, with numerous species demonstrating substantial genetic differentiation between populations within relatively small geographical distances [5,6]. Phasmatodea (stick insects), encompassing over 3000 extant species of terrestrial herbivores, primarily possess a tropical and subtropical distribution and largely consist of flightless species [8]. Phasmids are recognized for their limited dispersal capacity, with approximately 60% of all phasmid species either displaying significantly reduced wings or entirely lacking wings in their adult form [9]. Furthermore, phasmids with wings may still be incapable of sustained flight [10]. Stick insects have evolved several pivotal adaptations to counterbalance their loss of motility, such as masquerade crypsis and parthenogenesis [11]. The transition from sexual reproduction to parthenogenesis, a form of asexual reproduction, might be correlated with the flightless nature of stick insects, which renders locating mating partners more challenging [11].
Stick insects display remarkable masquerade crypsis as a defensive mechanism against visually hunting avian predators, by morphologically and behaviourally imitating twigs, bark, lichen, moss and leaves [11,12]. Additionally, stick insects use camouflage techniques in their eggs, which often closely resemble plant seeds. The eggs of some stick insect species not only imitate seed appearances but also employ analogous dispersal mechanisms. Numerous stick insect eggs feature a specialized knob-like structure known as a capitulum, which closely resembles the elaiosome of ant-dispersed seeds in both form and function and shares a similar chemical composition [13]. It is hypothesized that both elaiosomes and capitula, being lipid-rich, have evolved to promote ant-mediated dispersal. Furthermore, the considerably hardened egg capsule in Euphasmatodea (including all phasmids except Timema) is regarded as a crucial innovation in phasmid evolution [11]. This allows the eggs to endure potentially damaging falls from the canopy and remain buoyant on seawater for extended periods [11,14].
The sturdiness of phasmid eggs might allow them to remain viable even when contained within gravid female stick insects consumed by avian predators [12,15]. This mechanism is unattainable for many other insect species, as they generally fertilize their eggs immediately before oviposition, using sperm stored in the female's seminal vesicle after copulation. Conversely, numerous stick insects exhibit parthenogenesis, allowing them to produce viable eggs without fertilization [16]. In such instances, predation on gravid female stick insects could facilitate offspring dispersal, similar to the internal seed dispersal seen in plants producing fleshy fruits consumed by frugivorous birds. Nonetheless, it is vital to acknowledge that stick insects have developed a cryptic appearance as a means to avoid predation, rather than actively attracting animals [12]. Furthermore, we note that the study was based on a laboratory feeding experiment, necessitating caution in generalizing the findings to wild populations.
Due to their infrequent occurrences, it is likely challenging to demonstrate that avian predation serves as a factor promoting dispersion in natural settings [1]. One method to indirectly evaluate the importance of historical avian-mediated dispersal by analysing the spatial genetic structures of species with limited mobility, which would otherwise exhibit significant genetic differentiation among populations [1,17,18]. For example, mitochondrial DNA genetic comparisons have implied a minimum of two successful dispersal events between the Pacific and the Atlantic, traversing the digestive systems of birds, as the most plausible explanation for the long-distance dispersal of a certain marine snail [1]. Correspondingly, other investigations have proposed that long-distance dispersal facilitated by birds represents the most parsimonious explanation for the phylogeographic and biogeographic patterns of organisms with limited active dispersal capacities [17–21].
The focus of our study is on Ramulus mikado, a predominantly parthenogenetic stick insect species found commonly in Japan, with only a few documented instances of males [22]. Notably, avian endozoochory has been demonstrated in R. mikado (= R. irregulariterdentatum) eggs in experimental conditions (figure 1) [12], suggesting the possibility of historical avian-mediated dispersal in natural settings. In this investigation, we analyse the population genetic structure of R. mikado to determine the influence of potential geographical barriers on their phylogeographic structure. Fortunately, Japan has undergone extensive research in insect phylogeography, providing a wealth of prior studies for comparison [23,24]. These studies indicated that the division of terrestrial habitats by numerous oceanic straits in Japan has significantly limited insect migration and promoted genetic isolation among populations [23,24]. These studies also showed that the rugged volcanic terrain and diverse altitudes within the Japanese Islands contribute further to isolation [23,24]. Consequently, most insect species exhibit population genetic structures reflecting the geohistory of the Japanese Islands [25–30], while a few species with strong dispersal abilities show weaker genetic differentiation [31–33].
Figure 1. (a) The female adult of the stick insect Ramulus mikado. (b) The brown-eared bulbul (Hypsipetes amaurotis) feeding R. mikado to its chick. (c) H. amaurotis faecal pellets containing intact Ramulus mikado eggs. (d) Intact Ramulus mikado eggs recovered from H. amaurotis faeces. (e) First instar nymph of R. mikado hatched from the excreted egg. Scale bars: 50 mm (a), 100 mm (b) and 2 mm (c–e).
Here, we used a multifaceted molecular approach, incorporating mitochondrial sequences, simple sequence repeat (SSR) markers and genome-wide single nucleotide polymorphisms (SNPs) to determine the phylogenetic relationships among individuals collected throughout the species distribution range. Mitochondrial sequence analysis with maternally inherited markers is anticipated to reveal a genetic structure similar to that of sexually reproducing species, even in predominately parthenogenetic species. SSR markers, which have high mutation rates, are probably appropriate for detecting limited intraspecific genetic variation in parthenogenetic species. Additionally, genome-wide SNP markers can offer more reliable data for estimating kinship relationships among individuals sharing a common parthenogenetic ancestor.
By integrating the results from these three genetic markers, each with distinct advantages, we have examined the intraspecific phylogeographic pattern of parthenogenetic and flightless stick insect species with unparalleled resolution. Consequently, we discovered unique phylogeographic patterns in the flightless stick insect, such as the identification of identical mitochondrial DNA cytochrome oxidase subunit I (COI) haplotypes across significant distances, which do not contradict the hypothesis of passive long-distance genotypic dispersion.
2. Material and methods
(a) Study species, sample collection and DNA extraction
Ramulus mikado is a predominantly parthenogenetic and flightless stick insect widely distributed throughout Japan (figure 1) [22,34]. While the genus Ramulus encompasses over 100 species, R. mikado represents the sole species within this genus in Japan [22,35]. Between 2014 and 2018, we collected 67 R. mikado specimens from two islands in the Japanese Archipelago (figure 2; electronic supplementary material, table S1). Males are rarely documented in this species, and all the specimens we collected were females. Specifically, we obtained 31 individuals from Shikoku Island (SI) and 36 from Honshu Island. Due to the extensive sampling area, we divided Honshu Island into West Honshu (WH) and East Honshu (EH), based on the Itoigawa-Shizuoka Tectonic Line at the westernmost side of the Fossa Magna in central Honshu. The genetic structure of numerous insect groups reflects the separation of eastern and western Japan in the Fossa Magna region, where the current geological structures were formed 0.7–1.0 Ma [24,25].
Figure 2. Map illustrating the sampling locations of Ramulus mikado used in this study. The map was generated using maptools in R software (version 3.6.2; http://cran.r-project.org/) with data sourced from the National Land Numerical Information, Ministry of Land, Infrastructure, Transport and Tourism, Japan (https://nlftp.mlit.go.jp/index.html). Samples collected within a distance of fewer than 10 km were combined.
For data analysis, three populations were defined based on their respective regions: SI, WH and EH. To obtain genomic DNA, R. mikado samples were preserved in absolute ethanol. Subsequently, the Gentra Puregene Tissue Kit (Qiagen) was employed to extract total genomic DNA, following the manufacturer's instructions. This extracted DNA was used for the subsequent molecular analysis.
(b) Mitochondrial DNA analysis
For the analysis of genetic diversity, phylogenetic relationships among haplotypes, and their geographical distribution among three populations, we determined the partial sequences of the mitochondrial DNA cytochrome oxidase subunit I (COI) gene. The COI gene fragment was amplified and sequenced using primers LCO_Phas3 (5′– AAC TCA GCC ATT TTA CTA ATG AAA CG –3′) and HCO_Phas1 (5′– TAT ACT TCT GGA TGA CCA AAA AAT CA –3′), which were designed based on the sequences of related taxon found in the DNA Data Bank of Japan (DDBJ). PCR amplification was conducted in a GeneAmp PCR System 2700 thermal cycler (Applied Biosystems, Foster City, CA, USA). The PCR products were purified by using a High Pure PCR product purification kit (Roche Diagnostics, Mannheim, Germany), and purified products were sequenced directly using an ABI BigDye Terminator Cycle Sequencing Kit v. 3.1 (Applied Biosystems) on the ABI PRISM 3130 Genetic Analyzer (Applied Biosystems). Both strands of the amplified PCR products were sequenced, and electropherograms were assembled using Finch TV (http://www.geospiza.com/finchtv/). The assembled sequences were aligned with all analysed samples using CLUSTAL W [36] with default settings in MEGA-X [37]. We determined the mtDNA haplotype based on this aligned data.
To compare mitochondrial genetic variation among populations, we calculated the number of haplotypes, haplotype diversity and nucleotide diversity using DnaSP 6 [38,39]. We employed analysis of molecular variance (AMOVA) to investigate genetic variance and differences among populations. Overall and pairwise FST values were calculated, and the probability of each pairwise FST value being different from zero was tested based on 999 permutations using Arlequin version 3.5 [40]. To evaluate the genetic relationships among genotypes, we constructed a neighbour-net network using SplitsTree version 4.15.1 [41] and performed principal coordinates analysis (PCoA) using GenAlEx 6.5 [42]. Neighbour-net networks were generated from COI gene sequence alignments based on an uncorrected p-distance matrix. The genetic distance for PCoA was computed employing the maximum composite likelihood model [43] with MEGA-X [37]. The association between geographical and genetic distance was evaluated for (1) all populations, (2) SI and WH populations, (3) WH and EH populations, and (4) each of the three populations using Mantel tests by GenAlEx version 6.5.
(c) Simple sequence repeat analysis
We successfully isolated 13 SSR loci from the nuclear genome of R. mikado and designed corresponding SSR primers (electronic supplementary material, note S1). We employed these newly designed primers to determine genotypes for analysing genetic diversity among the three populations, phylogenetic relationships among multilocus SSR genotypes and their geographical distribution. We performed PCR amplification of the 13 SSR loci using 5 µl reactions with the QIAGEN Multiplex PCR Kit and a fluorescent dye-label protocol [44]. Each reaction contained 10 ng of genomic DNA, 2.5 µl of Multiplex PCR Master Mix, 0.01 µM of forward primer, 0.2 µM of reverse primer and 0.1 µM of fluorescently labelled primer. The amplification protocol consisted of 95°C for 15 min, 33 cycles at 94°C for 30 s, 57°C for 1.5 min and 72°C for 1 min, followed by an extension at 60°C for 30 min. We determined product sizes using an ABI PRISM 3130 Genetic Analyzer and GeneMarker software (SoftGenetics, State College, PA, USA).
We established multilocus genotypes based on 13 SSR loci genotypes. The probability of identical genotypes resulting from sexual reproduction was calculated using GenAlEx 6.5 [42]. For each population, genetic diversity was assessed by evaluating the average number of alleles per locus (A), allelic richness (RS), observed heterozygosity (HO), expected heterozygosity (HE) and the fixation index (FIS). These parameters were calculated using GenAlEx 6.5, except for allelic richness, which was calculated using FSTAT 2.9.3 [45]. To evaluate the genetic relationships among genotypes, we constructed a neighbour-net network using SplitsTree version 4.15.1 [41] and performed PCoA using GenAlEx 6.5. Nei's genetic distance DA [46] among multilocus genotypes calculated by POPULATION 1.2.30 [47] was used for the genetic distance of the network construction and PCoA. Additionally, the association between geographical and genetic distance DA was evaluated for (1) all populations, (2) SI and WH populations, (3) WH and EH populations, and (4) each of the three populations through Mantel tests. Mantel tests were performed using GenAlEx.
(d) MIG-seq analysis
MIG-seq is a recently developed genome-wide genotyping methodology employing a high-throughput sequencing platform [48]. This technique is a microsatellite-associated DNA sequencing approach, a form of reduced representation sequencing that encompasses restriction site-associated DNA sequencing (RAD-seq) [48]. MIG-seq has recently emerged as a potent instrument for population genetics research [49,50]. A MIG-seq library was prepared following the protocol suggested by the development team [51] and sequenced using the MiSeq system (Illumina, San Diego, CA, USA) and MiSeq Reagent Kit v3 (150 cycle). The raw genome-wide SNP data were archived in the DDBJ Sequence Read Archive (DRA, accession number DRA016238).
Upon eliminating primer sequences and low-quality reads [52], 9 558 552 reads (144 827 ± 3256 reads per sample) were acquired from 9 878 722 raw reads (149 678 ± 3360 reads per sample). The Stacks 2.62 pipeline was employed for de novo SNP discovery [53], using the following parameters: minimum depth of coverage required to generate a stack (m) = 3, maximum distance permitted between stacks (M) = 2 and the number of mismatches allowed between sample loci during catalogue construction (n) = 2. SNP sites containing fewer than three minor alleles were filtered, and only SNPs retained by 33 or more samples were extracted. We restrict data analysis to only the first SNP per locus to avoid linkage between SNPs. The SNP filtering for excess heterozygosity was not performed because of the predominantly parthenogenetic reproduction of R. mikado [22]. Ultimately, 980 SNPs were procured for subsequent analyses. During these processes, one sample (EH09) was excluded due to its high missing rate. The nucleotide diversity (ND), observed and expected heterozygosity was assessed by the program populations of Stacks. The number of allelic differences between two individuals was calculated using R package poppr 2.9.4 [54]. Overall and pairwise FST values were calculated, and the probability of each pairwise FST value being different from zero was tested based on 999 permutations using GenAlEx 6.5 [42]. A neighbour-net network was also constructed by employing the uncorrected p-distance matrix and disregarding ambiguous sites, with the use of SplitsTree version 4.15.1 [41]. Moreover, to assess potential genetic structure, a PCoA was performed based on the genome-wide SNPs using R package dartR 2.7.2. [55,56]. The correlation between geographical and genetic distance was assessed for (1) all populations, (2) SI and WH populations, (3) WH and EH populations, and (4) each of the three populations through Mantel tests and GenAlExversion 6.5.
3. Results
(a) Genetic variation and difference among R. mikado populations
Genetic analysis based on mitochondrial COI sequences and nuclear SSR markers revealed a notable accumulation of mutations in R. mikado due to parthenogenetic reproduction. The observed number of alleles at the 13 newly developed loci ranged from 1 to 20 (electronic supplementary material, table S3). The observed heterozygosity for eight loci was either 0 or nearly 0, while the remaining five loci displayed high values ranging from 0.64 to 1.00. From the 67 samples, 55 multi-locus genotypes were identified. The combined random match probability for the 13 loci was 5.53 × 10−6, with a high power for discriminating among individuals. The MIG-seq-based SNP heterozygosity (the highly heterozygous and homozygous loci) corresponds to these SSR heterozygosity patterns (electronic supplementary material, figure S2).
In the mitochondrial COI region sequences, 39 haplotypes were identified, among which 11 haplotypes were present in multiple individuals (table 1). The number of individuals sharing the same haplotype ranged from 2 to 6, and no specific haplotype was predominantly distributed. Intriguingly, some individuals collected from distant sites exhibited the same haplotype, and eight haplotypes were distributed in sites more than 10 km apart (electronic supplementary material, figure S1; table 1). The farthest Hap04 was confirmed from 683 km, and Hap06 was confirmed from 452 km. In the nuclear SSR genotype, the number of individuals exhibiting identical multi-locus genotypes was small, with only 2 or 3 observed, and no widespread specific genotype was observed. Gen34 and Gen30 genotypes were collected from 19 km and 10 km away, respectively (table 1). Two individuals exhibiting the Gen34 genotype also shared the same mitochondrial haplotype. Although the SNP analysis did not obtain the same genotype, individuals with close genetic distances across distant populations have been confirmed. The number of allelic differences between two individuals ranged from 40 (2.04%) to 136 (6.94%). The minimum number of allelic differences between two individuals was observed between EH13 and WH03, which were 259 km apart. The values of haplotype diversity HD and the nucleotide diversity of COI sequences, as well as the allelic richness RS of SSR markers, which are crucial indicators for assessing genetic diversity, were relatively high in the SI and WH populations and low in the EH population (table 2). Nonetheless, no significant differences were observed in the nucleotide diversity of genome-wide SNPs among the populations. AMOVA analysis based on COI haplotypes and SSR genotypes revealed significant genetic differences among the three populations (table 3). The FST values based on mitochondrial sequence data, SSR genotype data, as well as genome-wide SNP data, were all significant among the populations (p < 0.01). Each FST value calculated from mitochondrial sequence data and genome-wide SNP data was significant for every pair of populations, while the pairwise FST value between the EH and WH populations did not exhibit a significant difference when assessed using SSR genotype data (p = 0.218; table 3).
| haplotype or genotype | n | max. dist. (km) | sample ID |
|---|---|---|---|
| mitochondrial DNA | |||
| Hap04 | 4 | 683 | WH12, EH10, EH11, EH12 |
| Hap06 | 3 | 452 | SI26, SI27, WH20 |
| Hap01 | 5 | 179 | EH01, EH02, EH03, EH05, EH06 |
| Hap05 | 6 | 76 | WH19, EH13, EH14, EH15, EH16, EH17, EH18 |
| Hap25 | 5 | 52 | SI07, SI09, SI16, SI17, SI18 |
| Hap34 | 3 | 35 | SI22, SI23, SI24 |
| Hap32 | 2 | 18 | SI19, SI21 |
| Hap20 | 2 | 16 | SI01, SI03 |
| Hap09 | 3 | <1 | WH22, WH23, WH24 |
| Hap03 | 3 | <1 | EH07, EH08, EH09 |
| Hap31 | 2 | <1 | SI14, SI15 |
| nuclear SSR | |||
| Gen34 | 2 | 19 | SI09, SI17 |
| Gen30 | 2 | 10 | SI04, SI05 |
| Gen01 | 3 | <1 | EH01, EH02, EH03 |
| Gen04 | 3 | <1 | EH07, EH08, EH09 |
| Gen07 | 3 | <1 | EH13, EH16, EH17 |
| Gen11 | 2 | <1 | WH22, WH23, WH24 |
| Gen03 | 2 | <1 | EH05, EH06 |
| Gen06 | 2 | <1 | EH11, EH12 |
| population | mitochondrial sequence |
nuclear SSR |
genome-wide SNP |
|||||||
|---|---|---|---|---|---|---|---|---|---|---|
| NH | HD | ND | A | RS | HO | HE | ND | HO | HE | |
| all | 39 | 0.972 | 0.0077 ± 0.0005 | 3.0 | 3.5 | 0.34 | 0.32 | 0.00135 ± 0.00005 | 0.670 | 0.393 |
| SI (n = 31) | 20 | 0.961 | 0.0091 ± 0.0006 | 3.9 | 3.4 | 0.37 | 0.32 | 0.00136 ± 0.00005 | 0.667 | 0.390 |
| WH (n = 18) | 16 | 0.980 | 0.0075 ± 0.0008 | 3.4 | 3.4 | 0.34 | 0.35 | 0.00136 ± 0.00005 | 0.669 | 0.385 |
| EH (n = 18) | 5 | 0.797 | 0.0035 ± 0.0004 | 2.4 | 2.3 | 0.30 | 0.32 | 0.00139 ± 0.00005 | 0.677 | 0.390 |
| mitochondrial DNA | nuclear SSR | genome-wide SNP | |
|---|---|---|---|
| among populations | FST = 0.10*** | FST = 0.04** | FST = 0.012*** |
| SI–WH | FST = 0.11*** | FST = 0.03* | FST = 0.007** |
| SI–EH | FST = 0.12*** | FST = 0.05** | FST = 0.010*** |
| WH–EH | FST = 0.07* | FST = 0.02 | FST = 0.010*** |
(b) Genetic relationship among individuals and their spatial distribution
Our molecular analyses, which employed COI haplotypes, SSR markers and MIG-seq analysis, revealed limited associations between lineage and geographical distribution (figure 3–5), with some exceptions, such as COI haplotypes and genome-wide SNPs consisting exclusively of samples from the SI population (figures 3a,c and 4a,c). However, aside from these few instances, the correspondence between phylogenetic relationships and distribution was not evident. In all analyses of genetic markers, individuals with close genetic relationships were found to be dispersed across different populations. This pattern remained consistent across distinct methods employed in phylogenetic analysis and measures of genetic distance, such as the construction of a most parsimonious network based on COI sequence data, as well as the use of neighbour-joining trees based on Nei's genetic distance [46] and Cavalli-Sforza & Edwards's chord distances [57].
Figure 3. Neighbour-net network of Ramulus mikado samples reconstructed based on (a) the COI haplotypes, (b) nuclear SSR and (c) genome-wide SNP. Branch length denotes the average number of substitutions per site. Figure 4. Principal coordinate analysis plot based on (a) COI haplotypes, (b) nuclear SSR and (c) genome-wide SNP data. Figure 5. Relationships between geographical and genetic distances at the individual level based on (a) COI haplotypes, (b) nuclear SSR and (c) genome-wide SNP data.
The Mantel test results, examining the correlation between genetic distance and geographical distance within specific regions, yielded different outcomes among regions (figure 5). The Mantel test, using all samples, revealed no significant correlation between geographical distance and genetic distance in the mitochondrial sequence and genome-wide SNP data (figure 5a,c). Although a weak yet significant correlation was detected in the SSR data (p < 0.05; figure 5b), the R2 value was low (R2 = 0.041). Notably, genetically very similar individuals are found up to approximately 700 km in the mitochondrial sequence data and up to about 600 km in the SSR data. The combinations of individuals with the highest level of genetic distance between individuals were observed even at a distance of about 100 km in all three genetic markers. Within the samples merging Western Honshu (WH) and SI, ‘isolation by distance' was not observed based on COI haplotypes, SSR and SNP data (electronic supplementary material, figure S3). By contrast, within Eastern Honshu (EH), genetic distance exhibited a significant positive correlation with geographical distance based on all molecular methods. This correlation was particularly pronounced in SSR markers (electronic supplementary material, figure S3, p < 0.05, R2 = 0.678).
4. Discussion
Although the majority of research on avian-mediated dispersal has primarily focused on seed dispersal, birds can disperse a diverse range of invertebrates via endozoochory. Extreme cases include the transportation of insect eggs or larvae, as well as fish eggs, potentially facilitating the colonization of novel habitats [12,58–60]. Nonetheless, it is crucial to recognize that these investigations have primarily showcased passive egg dispersal within laboratory settings. In our comprehensive molecular analysis, spanning different time scales, we reveal evidence suggesting that long-distance dispersal, probably facilitated by avian predation, can influence the distribution and population structure of R. mikado (figure 3–5; electronic supplementary material, S1; table 1).
Our SSR analysis shows high heterozygosity in five loci and zero or near-zero in the remaining eight loci in R. mikado (electronic supplementary material, table S3). This pattern, also observed in genome-wide SNP data (electronic supplementary material, figure S2), suggests the predominance of asexual reproduction, though historical cryptic gene flow cannot be excluded [61]. This assumption is consistent with not only the female predominance in R. mikado [22] but also the non-functionality of the rare males [62]. The pattern of heterozygosity contrasts with terminal fusion automixis or gamete duplication, causing tremendous heterozygosity loss [63,64]. Automixis with central fusion seems most plausible, as SSR analysis revealed mothers and offspring (embryos) mostly shared genotypes, with rare heterozygous to homozygous transitions due to recombination [62]. The mixture of highly heterozygous and homozygous loci probably results from varying recombination probabilities by loci [64].
We identified some differences in COI haplotype frequency and SSR allele distribution among R. mikado individuals (tables 2 and 3). In parthenogenetic species, a single strain or few strains with limited genetic variation tend to have widespread distribution due to efficient parthenogenic reproduction and rapid expansion in a short evolutionary time scale [65,66]. By contrast, the accumulation of genetic variation in R. mikado implies a relatively long parthenogenetic persistence, allowing mutation accumulation. Given insect mitochondrial DNA substitution rates, including the COI region, range from 1.5% [67] to 2.3% [68] per million years, the 0.77% nucleotide diversity of COI sequences probably reflects a history of 0.34–0.51 Myr. Although multiple parthenogenetic origins could account for these differences, SSR marker or genome-wide phylogeny supports a single lineage radiation pattern, inferring a single parthenogenesis origin in R. mikado.
Although R. mikado exhibits a certain degree of genetic variation, only a weak geographical signal was detected among R. mikado individuals with limited active dispersal ability. This is highly unusual, as almost all the Japanese insects (particularly with limited dispersal abilities) exhibit significant genetic differentiation even on a small spatial scale [25–30]. The notion of ‘isolation by distance' serves as a proxy for assessing a species's dispersal ability [17,69]. This concept suggests that when individual dispersal distances are smaller, genetic drift acting on neutral genetic markers eventually leads to a positive correlation between genetic differentiation among locations and the geographical distance that separates them [69]. Notably, no positive correlation between geographical and genetic distances was observed in SI and WH populations, suggesting genotypic dispersal within these sea-separated populations. The detection of COI haplotypes Hap04 and Hap06 at distances of 680 km and 450 km, respectively, implies rapid expansion of these strains over hundreds of kilometres across the Fossa Magna, whose current geological structures had already formed before the origin of R. mikado (0.7–1.0 Ma versus 0.34–0.51 Ma) [24,25], outpacing mutation. The PCoA plots, as well as neighbour-net networks, based on COI and genome-wide SNP data have also revealed some distinct area-specific lineages in the SI population, while other SI lineages mixed with WH and EH populations, suggesting that some genotypes suffer long-distance dispersal. Nonetheless, these patterns also suggest that long-distance dispersal events are infrequent, as the genetic structure between populations remains partially intact, and some regional genetic differentiation is still discernible.
Distinct patterns in the EH population, which is probably a recently formed population following the last glacial period, may also support the infrequency of long-range dispersal. Pollen analysis indicates that coniferous forests predominantly covered this region during the last glacial period, while temperate broadleaf forests, which serve as suitable habitats for R. mikado, did not expand until after the last glacial period [70]. The low genetic variation of mitochondrial DNA and nuclear SSR markers in the EH population (table 2) probably reflects the relatively brief history of this population and the limited number of founders that have arrived from WH and SI populations. Notably, a distinct positive correlation between geographical and genetic distance was observed in the EH population based on SSR markers (electronic supplementary material, figure S3). The isolation-by-distance pattern in the EH population may be attributed to the limitation of long-distance dispersal and the higher mutation rate of SSR markers. Generally, the genetic structure of the EH population might reflect its relatively brief history of distribution expansion and more rapid accumulation of mutations in SSR loci compared to the rate of distribution expansion.
As discussed earlier, the low and mostly non-significant genetic differentiation among R. mikado populations is probably facilitated by long-distance dispersal. Internal dispersal through the digestive tract of stick insect predators is a plausible mechanism for dispersal, particularly considering the demonstrated resilience of their eggs in withstanding passage through avian digestive tracts [12]. Oceanic dispersal may have also contributed to the phylogenetic pattern, as certain phasmid eggs, like those of the Mascarene stick insects, can adhere to branches or leaf surfaces and be transported across oceanic currents [71]. Additionally, some phasmids, like Megacrania, possess sponge-like eggs that can float for long periods without needing to be attached to vegetation [14]. However, these dispersal methods are unlikely in R. mikado since their eggs do not possess either a sponge-like structure or the ability to stick to the branches.
While anthropogenic translocation is another possible explanation [65,72,73], it is highly improbable that ancient anthropogenic translocation occurred tens of thousands or even thousands of years ago, given the limited association between R. mikado and humans. Thus, the somewhat plausible human-mediated dispersal of R. mikado is accidental long-distance transportation with plant seedlings in the past few centuries, coinciding with the advent of steam engines. In such recent human-mediated transportation, genetic variation typically does not accumulate after dispersal, and specific genotypes associated with human transport can be identified [74,75]. However, genetic analysis of R. mikado specimens indicates the presence of certain levels of genetic variations, with several identical or very similar genotypes widely but disjunctly intermixed across long distances. This pattern deviates from the typical patterns observed in recent human-mediated dispersal.
Therefore, internal dispersal via predation provides a plausible explanation for the unexpected phylogeographic pattern. In fact, similar non-significant geographical differentiation can be observed in certain plants dispersed by birds in Japan, where the majority of haplotypes are widespread throughout the country [76,77]. However, we also note that, despite the seed-like appearance of the eggs, the possibility of egg dispersal through granivorous birds seems unlikely, as these birds have evolved to crush seeds in their gizzards, leading to complete digestion of the eggs [78]. Consequently, the predation of gravid females carrying eggs is a probable internal dispersal mechanism [12]. Given that parthenogenetic reproduction facilitates the internal dispersal of eggs within gravid female stick insects [12], similar passive dispersal might be prevalent among parthenogenetic stick insects. Intriguingly, Clitarchus stick insects in New Zealand show distinct genetic structuring in populations with sexual reproduction, whereas populations with parthenogenetic reproduction, which enables internal dispersal, demonstrate a much weaker structure [79]. The contrasting pattern of genetic differentiation may indicate long-distance internal dispersal within a parthenogenetic lineage. Although the ability to colonize without mating partners might also contribute to this pattern [80,81], the sudden increase in dispersal ability in the parthenogenetic lineage also necessitates consideration of other drastic factors.
Various bird species, including the large-billed crow (Corvus macrorhynchos), Eurasian jay (Garrulus glandarius), brown-eared bulbul (Hypsipetes amaurotis), bull-headed shrike (Lanius bucephalus) and eastern buzzard (Buteo japonicus), have the potential to aid in the passive dispersal of R. mikado due to their consumption of stick insects in Japan [12,82]. For instance, some individuals of Hypsipetes amaurotis migrate southward during autumn and winter [83]. Considering its flight speed and gastrointestinal passage time, they could theoretically transport stick insects over several kilometres [12]. While crows do not migrate seasonally, they regularly travel several kilometres between their roosting sites and foraging locations [84]. As large avian species such as crows have longer gastrointestinal retention time, they also contribute to dispersal over several kilometres [85]. Omnivorous mammals might also play a role in stick insect egg dispersal through their digestive tracts. The primate Macaca fuscata and the marten Martes melampus feed on stick insects in Japan, and stick insect eggs are frequently excreted in the faeces of the latter [86,87]. These potential mechanisms of passive dispersal are likely to influence phylogeographic patterns observed in R. mikado.
Overall, we have identified distinct phylogeographic patterns, including instances where identical COI genotypes are found at sites that are geographically distant and disconnected. This supports the idea of passive long-distance dispersal of the genotypes. The question of how organisms with limited active dispersal capabilities achieve extensive distribution has captured curiosity since the time of Darwin [88]. Based on the phylogeographic pattern in conjunction with prior experimental evidence [12], we suggest that R. mikado eggs can potentially survive avian gut passage in the wild, enabling occasional long-distance dispersal. Our finding presumably provides new perspectives by challenging the longstanding notion that insects invariably perish when consumed by predators.
Ethics
This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility
The mitochondrial COI, SSR loci and genome-wide SNP data are available in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive: LC767280–LC767346 (http://getentry.ddbj.nig.ac.jp/getentry/na/LC767280-LC767346), LC623838–LC623855 (http://getentry.ddbj.nig.ac.jp/getentry/na/LC623838-LC623855) and DRA016238 (https://ddbj.nig.ac.jp/resource/sra-submission/DRA016238), respectively.
Supplementary material is available online [89].
Declaration of AI use
We have not used AI-assisted technologies in creating this article.
Authors' contributions
K.S.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, validation, visualization, writing—original draft, writing—review and editing; T.N.: investigation, resources, visualization, writing—review and editing; S.K.H..: formal analysis, investigation, resources, validation, visualization, writing—review and editing; S.F.: investigation, resources, writing—review and editing; K.I.: resources, supervision, writing—review and editing; Y.I.: supervision, writing—review and editing; Y.S.: supervision, writing—review and editing; S.K.: formal analysis, funding acquisition, investigation, project administration, validation, visualization, writing—original draft, writing—review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interest declaration
We declare we have no competing interests.
Funding
The work was funded by the JSPS KAKENHI, grant nos 18K19215 and 21K19108 to K.S.
Acknowledgements
The authors thank Junpei Haga, Katsumi Iwahori, Kota Sakagami, Mineki Yoshida, Naoyuki Nakahama, Osamu Tominaga, Seiki Honda, Shumpei Kitamura, Hajime Ohira and Yuta Mashimo for providing specimens. We are also grateful to Koya Shishido, Kazuki Kurita, Riku Shina, Yuya Nakazawa, Kai Sato, Ryuta Sato, Takako Shizuka, Hidehito Okada, Kazuma Takizawa, Hakuren Kato and Toshihito Takagi for their technical assistance. We also extend our appreciation to Tetsuro Yoshikawa, Shumpei Kitamura, Shota Sakaguchi, Hajime Ikeda, Koji Tojo, Takaya Iwasaki, Kyoko Aoki, Yoshiaki Tsuda, Takeshi Yokoyama, Tatsuya Fukuda, Osamu Miura, Naoyuki Nakahama and Taro Saito for their insightful discussions.
