During the outbreak season of E. onukii (June to October 2017), individuals were collected from 59 sites encompassing all 22 tea-growing provinces in China (Figure 1A; Supplementary Table S1). This comprehensive sampling approach facilitated the inclusion of a wide range of climate conditions, elevation gradients, and soil pH variations. Furthermore, the sampled locations were characterized by intricate topography, including towering mountains (e.g., Qinling Mountain), expansive plateaus (e.g., Tibetan Plateau), large rivers (e.g., Yangtze River), as well as nearby seas. This meticulous sampling scheme enabled the researchers to disentangle the influences of IBD, IBE, and IBR on the population genetic connectivity of E. onukii.

Field sampling was systematically conducted in old tea plantations situated in the prominent and historically significant tea-cultivating counties of each province. In each plantation, a minimum of 50 adult individuals were meticulously sampled from five to ten randomly selected plots. These selected plots were spaced at distances ranging from 20 to 100 m, and kept away from the edges of the tea plantation to prevent unintended sampling of other leafhopper species. To preserve the collected leafhoppers, individual specimens were placed in vials containing absolute ethanol. Subsequently, these samples were transported back to the laboratory and then stored at −80°C for further analyses.

Following a thorough morphological examination under a stereomicroscope (Qin et al., 2015) to validate their taxonomic identity, individual DNA extractions were performed with the TIANamp Micro DNA Kit (TIANGEN, Beijing, China). Subsequently, the extracted DNA was utilized for amplifying three mtDNA genes (COI, COII, and 16S rRNA), which successfully differentiated species of Empoascini leafhoppers or exhibited variation among E. onukii populations in previous studies (Demichelis et al., 2010; Fu et al., 2014; Zhou, 2014). Following the primers and protocols described by Demichelis et al. (2010), Zhou (2014), and Fu et al. (2014), the amplifications were conducted with a BIO-RAD C1000 TOUCH Thermal Cycler (Bio-Rad, Hercules, CA, USA). The resulting PCR products were then subjected to the analysis with an ABI 3730XL capillary sequencer (Applied Biosystems, Foster City, CA, USA). Each of the successfully amplified sequences was blasted with the complete mitochondrial genome (Liu et al., 2017) of E. onukii available in the NCBI database, and multiple sequence alignments were performed separately for each of the three genes through MAFFT (Katoh and Standley, 2013). The sequence data obtained from 1518 individuals sampled from 57 sites were utilized for subsequent analyses. However, individuals from the ES (Enshi) and RC (Rongcheng) sites were excluded from the analysis due to poor sequencing quality. To maintain consistent sequence lengths and ensure reliable alignment, the poorly aligned ends of the sequences were filtered, resulting in sequence lengths of 763 bp, 679 bp, and 435 bp for COI, COII, and 16S rRNA, respectively.

To obtain the most influential variables in predicting the mitochondrial genetic divergence of leafhoppers, we conducted a comprehensive analysis that simultaneously examined the impacts of IBD, IBE, and IBR.

For IBD (geographic distance, Supplementary Table S2), we employed the “Vincenty (ellipsoid)” method (Vincenty, 1975) implemented in the R package ‘geosphere’ (Hijmans et al., 2022) to calculate the shortest geographic distances between sites.

For IBE, we utilized climatic variables, elevation, and soil pH. The soil pH was particularly selected due to its significant association with tea plant growth (Han et al., 2018), which, in turn, can influence the productivity, distribution, and evolutionary dynamics of the leafhopper populations. For climatic and elevation factors, we obtained data of the current 19 bioclimatic variables (Supplementary Table S3; Fick and Hijmans, 2017) and elevation at a resolution of 30 s (approximately 1 km) from https://worldclim.org. Soil pH data at a resolution of 250 m (at a depth of 30 cm) were acquired from http://data.isric.org specifically for the pH analysis. The soil pH data were subsequently transformed to a resolution of 30 s in the R package ‘raster’ (Hijmans et al., 2023) for analyses. To mitigate collinearity or autocorrelation issues, we initially examined the correlations among these 21 variables in the R package ‘enmtools’ (Warren et al., 2010). Variables with correlations (|r|) below 0.7 were remained, while for the two variables with correlations (|r|) above 0.7, the one known to be important determinants of insect growth was selected for the following analyses. The selection ultimately resulted in the retention of soil pH and five bioclimatic variables (Supplementary Table S2): Isothermality, Max Temperature of Warmest Month, Precipitation Seasonality, Precipitation of Warmest Quarter, and Precipitation of Coldest Quarter. Through applying the R package ‘raster’, we extracted the environmental variation associated with each sampling site from these six layers at a resolution of 30 s.

For IBR, we generated seven resistance surfaces (Supplementary Table S2) to represent the hypothesized resistance of specific landscape features to leafhopper gene flow. These surfaces were created for mountains, rivers, elevation, habitat suitability of E. onukii, and availability of tea plants. Rivers and mountains were selected due to their known influence as barriers to gene flows in previous studies, particularly the Yangtze River, Huaihe River, Qinling Mountain, and Tibetan Plateau (e.g., Wei et al., 2015). For the rivers, we obtained a shape file of rivers from https://www.naturalearthdata.com and assigned a value of 1 to rivers and 0 to land, with the rivers defined as barriers, so as to generate a resistance surface (“riv_res”). For the mountain resistance surfaces, we utilized the K1 and K3 mountain resources from https://rmgsc.cr.usgs.gov/gme/. Three types of resistance surfaces were generated as follows: (a) through utilizing the K3 resource, we defined mountains above 300 m as barriers and employed this criterion as a threshold to assign low versus high resistance (assigning values of 0 and 1 to mountains below and above 300 m, respectively) (“mou300_res”); (b) with the K3 resource, we defined mountains above 900 m as barriers and used this criterion as a threshold to assign low versus high resistance (assigning values of 0 and 1 to mountains below and above 900 m, respectively) (“mou900_res”); (c) based on the K1 resource, we defined all mountains as barriers, with resistance increasing proportionally to mountain height, as well as assigning normalized values ranging from 0 to 1 (“mouclass_res”). In addition to investigating elevation dissimilarity in IBE, we generated a resistance surface to model the potential impact of topographic complexity on dispersal routes and gene flow patterns. Through the application of the same digital elevation data, we defined resistance values based on elevation, assigning normalized values ranging from 0 to 1 to generate a surface (“ele_res”). Finally, habitat suitability and tea plant availability were selected as resistance surfaces because these factors tended to restrict the potential dispersal routes and gene flows among populations. To obtain resistance surfaces related to habitat suitability and tea plant availability, we employed MaxEnt 3.4.1 (Phillips et al., 2006) to construct ecological niche models for E. onukii (“ensuit_res”) and C. sinensis (“teasuit_res”) as detailed in Li et al. (2023). From the seven generated resistance surfaces at a resolution of 30 s (ele_res, riv_res, mou300_res, mou900_res, mouclass_res, ensuit_res, and teasuit_res; Supplementary Table S2), we calculated resistance-distance values between each pair of the 57 locations through applying Circuitscape 5.0 (McRae et al., 2016) implemented in Julia, based on circuit theory. In the analysis, rivers, mountains, and elevation were considered as dispersal barriers, and thus regarded as “resistance” in the Circuitscape analysis, where higher surface values indicated increased limitation rates of leafhopper gene flows. Regarding habitat suitability and tea plant availability, we assumed that E. onukii could be more likely to occur (or to exist) in its suitable habitats and in habitats suitable for its host plant (tea). In addition, we assumed a stepping-stone dispersal model for E. onukii, based on its limited flight capabilities (Bian et al., 2014), and considering the observed patterns of isolation by geographic distance among the sampled populations. Consequently, high values of habitat suitability and tea plant availability were regarded as facilitators of dispersal and were employed as “conductance” parameters in the Circuitscape analysis. Greater surface values of these two factors indicated higher rates of gene flows in E. onukii.

Genetic diversity was firstly evaluated for each population through calculating the number of haplotypes (H) and haplotype diversity (h) in DnaSP 6.12.1 (Rozas et al., 2017), and estimating nucleotide diversity (π) with ARLEQUIN 3.5.2 (Excoffier and Lischer, 2010).

Previously, population differentiation was examined through the calculation of pairwise F_ST values (F_ST) in ARLEQUIN 3.5.2, with coalescent tree construction and phylogenetic relationship analysis conducted (Li et al., 2022). In this study, to further capture discrete genetic variation and explore potential correlations with geographical and landscape factors, we additionally calculated pairwise Prevosti’s genetic distance (P_GD; Prevosti et al., 1975; Kamvar et al., 2014) and Edwards’ genetic distance (E_GD; Cavalli-Sforza and Edwards, 1967) using the R package ‘adegenet’ 2.0.0 (Jombart et al., 2010), based on the concatenated three mitochondrial genes. P_GD has been recognized as an absolute genetic distance (Prevosti et al., 1975), while E_GD has been estimated as a Euclidean distance and has demonstrated robust performance in the presence of missing data (Bishop et al., 2021). Notably, F_ST differed from P_GD and E_GD in that it exclusively considered the data pertaining to only the two populations concerned, rather than encompassing the entire pool of population data simultaneously (Weir, 1990). Hence, all three distinct measures of genetic distance were estimated and subsequently utilized as response variables associated with spatial data.

Principal coordinates analysis (PCoA) was then performed using GenAlEx 6.503 (Peakall and Smouse, 2012) to better understand how the populations genetically grouped. Analyses were conducted based on the generated matrices of the P_GD, E_GD, and F_ST, respectively.

To identify non-IBD patterns in genetic differentiation across China, we employed the Estimated Effective Migration Surface (EEMS; Petkova et al., 2016) based on matrices of the P_GD, E_GD, and F_ST, respectively. The EEMS can model genetic variation by considering migration rates in relation to geographic factors, which can generate an estimated effective migration surface to reveal the regions where population divergence and gene flows deviated from the expected patterns of the IBD. By utilizing EEMS as an exploratory tool, we aimed to identify regions of the landscape that could potentially serve as dispersal corridors, exhibiting higher-than-expected gene flows, as well as those acting as biogeographic barriers which might impede gene flow. To capture the explicit and continuous genetic variations, we partitioned the landscape into demes and modeled the dispersal of individuals between neighboring demes. The EEMS analysis was run with a deme size of 1,200. Each analysis consisted of three independent chains, running for 5 × 10⁶ Markov chain Monte Carlo (MCMC) iterations, after a burn-in of 1 × 10⁶ iterations. The thinning was performed at intervals of 5,000 iterations, and distinct starting seeds were adopted for each run. To ensure convergence, posterior plots from independent runs were compared with the R package ‘rEEMSplots’ (Petkova et al., 2016). Subsequently, outputs from the three runs were combined and visualized.

In landscape genetics analyses, it has been recommended to employ various approaches and models to ensure robust conclusions (Balkenhol et al., 2015). Thus, to elucidate the relationship between dissimilarity in leafhopper mitochondrial variation among sites and the environmental distance (IBE), as well as the degree of isolation between these sites (IBD and IBR), we adopted a comprehensive approach which combined Mantel tests, multiple regression models, and machine learning methods.

To gain a better understanding of the interaction between the top predictor variable (mountain resistance) and genetic connectivity, we constructed generalized linear models (GLM) between populations genetic diversity (H, h, and π) and elevation using the R package ‘stats’ (R Core Team, 2022).

The genetic diversity (Supplementary Table S1) varied among populations, with H ranging from 7 (MH, BF) to 34 (TS), h from 0.639 (YX) to 1 (DY), and π from 0.0014 (MH) to 0.0102 (BS). The mean values of H, h, and π across populations were 16.842 (SD = 5.441), 0.898 (SD = 0.097), and 0.005 (0.001), respectively.

The estimations of P_GD, E_GD, and F_ST between sampling sites ranged from 0.02 to 0.30, 0.07 to 0.42, and −0.03 to 0.55, respectively, with the mean values of 0.07 (SD = 0.05), 0.16 (SD = 0.06), and 0.10 (SD = 0.12), respectively (Supplementary Figure S1).

The YX population was genetically unique in all of the three Principal Coordinate Analyses (PCoA; Figure 2). Additionally, results showed that most of the populations from the mountainous regions of Wester China (e.g., MH, PE, TC, PJ, MT, QR, MY, BS, DY) and some populations located at the South China region (e.g., BF, WZS, HZ, NC) tended to be genetically segregated from each other and also from other populations.

The EEMS analysis conducted on P_GD, E_GD, and F_ST (Supplementary Figure S3) consistently revealed the notable reductions in gene flows, across the majority of the populations inhabiting the mountainous regions of Western China, as well as the islands of Taiwan and Hainan. Additionally, similar patterns were observed among the populations situated along the middle and lower reaches of the Yangtze River.

The results from the simple Mantel tests indicated a significant correlation between the geographic distance and all three genetic distance measures. When controlling for the impact of geographic distance through applying the partial Mantel tests, the variables significantly associated with all three genetic distance measures consisted of isothermality, max temperature of warmest month, mou300_res, mouclass_res, and ele_res (Table 1).

The MRM analyses (Table 2) utilizing P_GD, E_GD, and F_ST yielded similar results, with the final model including only geographic distance and mouclass_res after the backward selection procedure. These two predictors significantly explained 27.38% of the variation in P_GD, 32.08% in E_GD, and 35.39% in F_ST.

The GDM models (Table 3; Supplementary Tables S4–S6) combining all possible predictor variables significantly explained 37.71%, 38.95%, and 42.65% of the variation in P_GD, E_GD and F_ST, respectively. The variables contributing significantly to the full models remained consistent irrespective of the genetic indices, including geographic distance and mouclass_res. The GDM models incorporating IBR variables demonstrated similar predictive power to the full model, while those exclusively utilizing IBD, IBE or both tended to be insignificant and accounted for a smaller proportion of the variation. Among models incorporating IBR variables, mouclass_res consistently emerged as the main predictor, effectively explaining a substantial portion of the observed genetic dissimilarity. In models formed solely by IBD, IBE, or both, the variables of importance were geographic distance, isothermality, and max temperature of warmest month, though their significance was not consistently observed.

The RDA models (Supplementary Figure S3; Supplementary T able S7) incorporating all variables and the optimal models with selected minimal variables exhibited significant associations with the three genetic distance measures. Moreover, the optimal models, in particular, accounted for a considerable proportion of the genetic variance, explaining 29.34% in P_GD, 28.49% in E_GD, and 13.32%in F_ST. When controlling for geographic location (longitude) with partial RDAs, the environmental variables included in the optimal models consistently predicted nearly the same amount of genetic variance as the full optimal models (22.19% for P_GD, 20.88% for E_GD, and 8.36% for F_ST). In contrast, when controlling for environmental factors, the percentage of variation explained by geographical location was substantially lower (4.64% for P_GD, 5.35% for E_GD, and 2.50% for F_ST), while still statistically significant. The variables involved in the optimal models varied depending on the specific genetic indices, with PCNM4, PCNM2, and longitude being consistently present for all three indices.

The RF models demonstrated a substantial R-squared value of 0.63 for P_GD, 0.65 for E_GD, and 0.69 for F_ST, which indicated that the predictor variables accounted for 63% of the genetic variance in P_GD, 65% in E_GD, and 69%in F_ST. While the increase in node purity analyses consistently highlighted the prominent role of mouclass_res in the RF models (Figure 1B), irrespective of the genetic indices, the projected maps vividly illustrated diminished genetic connectivity across the mountainous region of Western China (Figures 1C–H). For assessing the performance of the models through leave-one-out cross-validation, the mean RMSE for model evaluation (Supplementary Figure S4A) across all 57 runs was 0.03 (SD = 0.00), 0.04 (SD = 0.00), and 0.06 (SD = 0.00) for P_GD, E_GD, and F_ST, respectively. Furthermore, the spatial evaluation (Supplementary Figure S4B) yielded RMSE values of 0.03 (SD = 0.02), 0.04 (SD = 0.02), and 0.07 (SD = 0.04) for P_GD, E_GD, and F_ST, respectively.

The GLM analyses (Figure 3; Supplementary Table S8) demonstrated that both the number of haplotypes (H; r = –0.350, R^{2 =} 0.106, P = 0.008), haplotype diversity (h; r = –0.365, R^{2 =} 0.117, P = 0.007) significantly decreased with increasing elevation.

Previous studies have extensively documented the significance of long-standing landscape features, such as mountains and rivers, as prominent barriers to gene connectivity and drivers of genetic variation in various organisms (e.g., Noguerales et al., 2016; Trense et al., 2021). In the case of small-sized insect species characterized by limited dispersal capacity, strong associations between biogeographic barriers and genetic differentiation have been expected. Within our study region, the Yangtze River, as well as the complex topography of the Qinling Mountains and Tibetan Plateau, which are characterized by deep valleys and high mountains, have been previously recognized as effective impediments to gene flow in diverse species (e.g., Chilo suppressalis in Meng et al., 2008; Grapholita molesta in Wei et al., 2015; Pinus tabuliformis in Xia et al., 2018). Through our EEMS analysis, we observed distinct patterns of reduced gene flow on islands, along the course of the Yangtze River, and in the mountainous regions of Western China. Building upon these findings, in conjunction with our spatial clustering-based population structure analysis of E. onukii (Li et al., 2022), we postulated that landscape features might play a crucial role in driving genetic differentiation among E. onukii populations in China.

Our comprehensive analyses employing multiple models have consistently indicated the potential prominence of IBR in driving population divergence within E. onukii. Although the resolution of the raster layer representing river resistance (at 30-s resolution, approximately 1,000 m) may limit our ability to detect a direct association between functional connectivity and river permeability in this study, it is noteworthy that resistances associated with mountainous emerged as significant and predominant drivers of genetic variation in E. onukii. Additionally, our GLM analysis revealed that elevation was significantly and negatively associated with the number of haplotypes and haplotype diversity. And although the analysis did not detect a significant correlation between nucleotide diversity and elevation, populations at higher elevation tended to exhibit slightly higher nucleotide diversity. Nucleotide diversity is known as the average number of nucleotide differences per site in pairwise comparisons among DNA sequences, while haplotype diversity (also defined as gene diversity) refers to the probability that two randomly sampled alleles differ (Nei, 1987; de Jong et al., 2011). These lines of information suggest that the E. onukii populations at higher elevations show fewer haplotypes but greater nucleotide sequence difference from other populations, pointing towards genetically distinct populations with smaller effective size at high mountains compare those within low elevations. While this pattern observed here may be attributable to that the populations in high mountains more frequently experience sharp population contraction and strong environment selection due to the more variable climate there, it could also be a phenomenon formed because gene flow is restricted by mountain ridges as had previously been observed in other organisms (e.g., Rana luteiventris; Funk et al., 2005). Hence, mountains may impede genetic connectivity among E. onukii populations, aligning with the situations observed in other insect species that inhabit regions characterized by pronounced landscape heterogeneity (e.g., Noguerales et al., 2016).

Interestingly, the significant role of IBR in driving genetic divergence is not a universal outcome observed in numerous other organisms. For example, in Tetragonisca angustula bees across southern Brazil, the population divergence is primarily associated with ecological conditions and geographic distance rather than physical barriers (Francisco et al., 2017). Similarly, for 13 snake species across western North America, environmental heterogeneity rather than biogeographic barriers has been identified as the main driver of the population divergence (Myers et al., 2019). These findings suggest that species-specific traits play a crucial role in shaping population structure (Zamudio et al., 2016), and landscape genetic patterns are closely linked to life-history traits, dispersal capacities, and habitat preferences of each species (Myers et al., 2019). In the case of the leafhopper species we studied, the observed isolating mechanism of mountains in its populations could be attributed to its limited long-distance migration and weak dispersal ability. E. onukii serves as a small-sized insect, with adults measuring only about 3 mm in length (Qin et al., 2015). Although long-distance migration with the assistance of air streams has been documented in certain small insects such as Nilaparvata lugens (Hu et al., 2019), no such behavior has been observed in E. onukii. Due to their small size, E. onukii adults tend to be able to actively fly for only a few kilometers. Indoor tethered flight mill tests have confirmed the sedentary nature of E. onukii, as they tend to feed or perch behind tea leaves for most of their time and fly distances of less than 6 m (Bian et al., 2014). Furthermore, the species’ short generation time (9–15 generations per year; Qin et al., 2015) may also contribute to the observed patterns. It is widely recognized that species with short generation times are more likely to retain genetic signals left by historical dispersal routes interacting with habitat barriers and corridors (Noguerales et al., 2016). Additionally, the species’ limited feeding preferences likely augmented the pronounced isolating effects of mountains on its populations. E. onukii exhibits a marked preference for tea trees as its primary host plants (Shi et al., 2015). Feeding on other plants is sparsely reported, while elevation gradients have been evidenced as an important factor influencing the distribution of tea plants. Consequently, these host–plant associations may further limit population connectivity among mountain barriers (Noguerales et al., 2016). In addition to these species-specific traits, the successful detection of mountain resistances in this study may also be relevant with the marker type we employed, which will be discussed further in the following section.

Collectively, it is crucial to emphasize the significance of understanding species-specific traits in elucidating the landscape features that can promote gene flow and functional connectivity among populations, though determining the exact traits that contribute to these processes can be challenging. Furthermore, in theory, our findings can provide insights for E. onukii management strategies. Our multiple landscape genetics models demonstrated that mountains can serve as effective barriers to dispersal for E. onukii, which indicated a high likelihood of low population connectivity among mountainous regions. The potential for high isolation among mountains suggested that internal management efforts can be effective and highlight the demand of focusing on and controlling human-mediated dispersal from outside regions into these areas. Additionally, the presence of low-connected populations among mountains presented an opportunity to explore innovative pest control methods and to locally eradicate the leafhoppers. Moreover, the knowledge of mountains acting as dispersal barriers can aid in identifying regions with a low risk of reinvasion when planning spatially explicit management strategies. This information can also guide the placement of treatment technologies to block human-mediated dispersal from outside regions, such as the inter-regional transportation of tea plants. In fact, the recognition of isolated tsetse fly populations through population genetics has previously been applied to successfully plan pest control efforts in Senegal, which has led to the eradication of tsetse flies in specific areas and the opening of new agricultural zones (Bishop et al., 2021).

In contrast to the climate-related pattern observed in our previous study on microsatellites (Li et al., 2023), where ecology was proved to play a more prominent role than biogeography in population diversification, our mitochondrial analyses indicated that geography and gene flow were of more significance than ecology and selection in driving population differentiation. The combination of microsatellites and mtDNA sequences has been widely employed in population genetics studies across various species (e.g., Bonnin et al., 2023). However, a recent debate has emerged regarding the incongruence of the results obtained from these two markers and the preferred marker for population structure analysis or long-term monitoring of species in management programs (Muñoz-Valencia et al., 2023). Certain researchers have raised concerns about the utility of single-locus mtDNA and argued that limited mtDNA sequence data may not provide sufficient information to determine spatial genetic structure and a species’ true dispersal potential (e.g., Teske et al., 2018). On the other hand, others have reported that mtDNA still serves as an effective tool for broad-scale geographic studies, and simple mtDNA studies can serve as powerful initial surveys to gain insights into phylogeographic patterns and demographic history across environments. These insights can then guide more detailed investigations of selection and local adaptation in populations using nuclear genomic approaches (e.g., Bowen et al., 2014; Fonseca et al., 2023).

Despite the challenge of mitochondrial-nuclear microsatellite discordance and the criticism of single-locus studies, it has been widely acknowledged that these two markers evolve at different rates and can provide insights into different evolutionary events over time. Specifically, the evolutionary rates of mtDNA and microsatellites have been observed to be 5–10 and 100–1,000 times faster, respectively, than those of single-copy nuclear DNA (Wan et al., 2004). In other words, microsatellites evolve 10–200 times faster than mtDNA. The rapid mutation rate of microsatellites can lead to issues of saturation and homoplasy, where genetic information about past events may be effectively lost. As a result, microsatellites are more informative for studying recent and/or contemporary genetic mutations and are suitable for investigating patterns of genetic diversity and gene flow over shorter timescales. Comparatively, mtDNA retains more information about historical changes in genetic variation, making it the optimal option for resolving spatial patterns of gene flow and colonization over long timescales. These differences in evolutionary rates and the information captured by each marker may be the reason underlying the incongruent findings commonly reported among studies utilizing different markers, and can also explain the comparative results found between this study and our previous investigation based on microsatellites. The present distribution and connectivity of E. onukii populations across the landscape are the combined result of historical and contemporary forces. Long-standing natural barriers, such as mountains, represent historical mechanisms that have naturally and geographically isolated E. onukii populations, which has led to evolutionary processes such as changes in genetic diversity, population differentiation, and potentially even speciation. In addition to these long-term landscape features that have hindered population dispersal and exhibited their DNA signatures, contemporary factors, such as recent climate changes, can also impact the movement and gene connectivity of individuals among populations, leaving their marks on genetic structure. Thus, our mtDNA analyses successfully identified the prominent role of mountains as barriers to E. onukii population divergence, while the previous microsatellite study (Li et al., 2023) based on the same sampling system revealed that the spatial population structure of E. onukii has been primarily shaped by climatic differences. In light of these findings, it is evident that unrelated phenomena can leave their marks on genetic structure, and these different genetic signatures can be captured by diverse types of molecular markers. Hence, there is a demand of the combination of multiple markers to obtain a comprehensive understanding of intraspecific diversification and to facilitate the management of a species’ evolutionary history and contemporary interconnectedness.

Undoubtedly, both microsatellites and mtDNA genes are commonly utilized to investigate neutral processes and may not fully capture the genomic-wide variations associated with environmental adaptation (Álvarez et al., 2021). However, the identification of candidate adaptive loci is widely recognized as a valuable approach for gaining insights into potential adaptive divergence among populations and the capacity for local adaptation (Zimmerman et al., 2020). Although both marker types can generally inform us about population-level genetic differentiation and diversity, addressing questions related to local adaptation and evolutionary independence necessitates the inclusion of non-neutral loci. Notably, genomic-wide evidence can benefit the exploration of the candidate genes responsible for adapting to ecological gradients in natural populations, which consisting of not only responses to heat, but also precipitation, cold, and other environmental factors under strong selection pressures (Waldvogel et al., 2018). Considering E. onukii, conducting further studies through the application of the genomic-wide approaches would be valuable for elucidating the underlying mechanisms driving adaptive evolution in this species.