Sampling was conducted from 6 to 12 May 2021 in the Hiji River, Shikoku Island, Japan. The river originates from Torisaka Pass in Seiyo City, Ehime Prefecture, and flows into the Iyo Sea. The watershed area is 1,210 km². Sixty sites across different environmental conditions were selected for sampling (Fig. 1). We collected three samples (upstream, middle, and downstream) per site, which we processed separately for sorting, DNA extraction, and the first polymerase chain reaction (PCR) but combined into a single sample for the second PCR. Macroinvertebrates were collected using the kick-net method (0.5 mm mesh). The kick-net size was 30 cm in width, 25 cm in length, and 55–60 cm in depth. Sampling effort was distributed across four to five sampling points per site, with approximately 5 min of kick-net sampling per sampling point. The collected macroinvertebrates were preserved in 99.5% ethanol and transported to the laboratory for sorting. Sorted macroinvertebrates that were not processed for morphological identification were stored in 99.5% ethanol at 4 °C until DNA extraction. Environmental parameters, including water temperature, pH, electrical conductivity, and dissolved oxygen, were measured in situ using a portable meter (LAQUA, Horiba Corporation, Japan). Surface water was also collected for the measurement of ammonium nitrogen (NH₄-N), nitrate-nitrogen (NO₃-N), nitrite-nitrogen (NO₂-N), and phosphate-phosphorus (PO₄-P) using the QuAAtro 2-HR system (BeLtech).

We extracted DNA from ethanol-dried macroinvertebrate bulk community samples using the phenol-chloroform-isoamyl alcohol (PCI) method (Suppl. material 1: protocol S1). We used whole individuals and entire bulk community samples for DNA extraction. DNA was quantified using a NanoDrop spectrophotometer (NanoDrop 2000, Thermo Scientific, USA), quality-checked by agarose gel electrophoresis, and diluted to 20 ng/μL. We performed PCR amplification of the mitochondrial cytochrome oxidase subunit 1 (COI) gene using the primer set BF2 (5'-GCHCCHGAYATRGCHTTYCC-3') and BR2 (5'-TCDGGRTGNCCRAARAAYCA-3') with an expected amplicon size of 421 bp (Elbrecht and Leese 2017). A two-step library preparation method was employed following Elbrecht and Steinke (2019), with modifications (Suppl. material 1: protocol S2). We included negative (molecular-grade water) and positive controls (macroinvertebrate bulk-community DNA) in all PCR assays. We quantified the libraries using the KAPA Library Quantification Kit (Kapa Biosystems, USA) in a C1000 Touch Real-Time System thermal cycler (Bio-Rad Laboratories, CA, USA). We then pooled each sample in equimolar amounts (4 nM), cleaned the pooled libraries using SPRIselect (Beckman Coulter Inc., USA) with a left-sided selection ratio of 0.85×, and checked the quality using the Agilent Bioanalyzer 2100 system (Agilent Biotechnologies). High-throughput sequencing was performed on an Illumina MiSeq system with a 300 bp paired-end read kit v3 (Illumina, USA, MS-102-3003) at a final concentration of 8 pM with a 10% Phi-X spike-in.

Raw paired-end reads were quality-checked using FastQC (Andrews 2010) and then demultiplexed via JAMP (Elbrecht et al. 2018) in R v. 4.3.2 (http://github.com/vascoElbrecht/JAMP). After demultiplexing, the reads were processed through the DADA2 pipeline (Callahan et al. 2016), which included primer removal (trimLeft = 20, 20), trimming of forward and reverse reads (truncLen = 260, 220), filtering of erroneous reads (maxEE = 1, 3), denoising, merging of forward and reverse reads, and removal of chimeras. Amplicon sequence variants (ASVs) were further filtered to extract haplotypes. ASVs were filtered for the expected amplicon size (421 bp) using Cutadapt version 4.9 (Martin 2011), and sequences with fewer than 10 reads per sample were removed. We used Multiple Alignment of Coding Sequences (MACSE) version 2.07 (Ranwez et al. 2011; Ranwez et al. 2020) to assess the presence of internal frameshifts, stop codons, and deletions. ASVs were then translated into protein sequences using Geneious Prime version 2023.0.1 (https://www.geneious.com). We performed pseudogene filtering by comparing the open reading frame (ORF) to the arthropod amino acid COI gene hidden Markov model (HMM) profile using an HMM scan analysis (http://hmmer.org/) (Porter and Hajibabaei 2021). We discarded sequences with a bit score below the first quartile (Q1) – 1.5× interquartile range (IQR) as putative pseudogenes (Porter and Hajibabaei 2021). ASVs that passed these filtering steps were considered haplotypes.

In this study, operational taxonomic units (OTUs) were used to define species boundaries. Haplotype sequences were clustered into OTUs with 97% identity using the cluster_smallmem command in VSEARCH v. 2.21.1 (Rognes et al. 2016). Because OTUs can group haplotype sequences from multiple species (Martin et al. 2021), we screened all OTUs for signals of mixed-species ancestry. Initial screening was conducted based on the maximum pairwise p-distance. OTUs with a maximum pairwise p-distance of < 0.03 were considered single species, whereas those with a maximum pairwise p-distance of > 0.03 were considered potentially mixed species. OTUs with a maximum pairwise p-distance > 0.03 were further confirmed for mixed-species signals based on two criteria: genetic distance and phylogeny. The genetic distance method involved testing the modality of pairwise p-distance distribution of all haplotypes within the OTU using a histogram and excess mass-based multimodality testing according to Ameijeiras-Alonso et al. (2019), implemented in the multimode package v. 1.5 in R (Ameijeiras-Alonso et al. 2021). The phylogeny method was based on the Poisson Tree Process (PTP) (https://species.h-its.org/ptp/) (Zhang et al. 2013) (Suppl. material 1: protocol S3). We considered OTUs with a maximum pairwise distance of > 0.03 to be single species if the pairwise distance of all haplotypes within the OTU showed a unimodal distribution based on a histogram and excess mass-based multimodality testing and demonstrated a monophyletic clade with strong maximum likelihood (ML) and Bayesian support values (≥ 75) based on PTP. Otherwise, OTUs that did not meet both criteria were classified as mixed species.

For taxonomic assignment, haplotype sequences were first queried against the Barcode of Life Database (BOLD v. 5, http://www.boldsystems.org/; accessed October 2025) using a ≥ 97% identity cutoff. Sequences that did not return a taxonomic assignment using the BOLD database were then queried against an in-house database using the --usearch_global command in VSEARCH v. 2.21.1 (Rognes et al. 2016), with a ≥ 97% identity threshold. The in-house database consisted of Sanger sequences of morphologically identified adult crane flies and was deposited in the DNA Data Bank of Japan (accession numbers LC898117–LC898137). Sequences that did not have taxonomic identification using BOLD or the in-house database were further searched against the National Center for Biotechnology Information (NCBI) non-redundant database (accessed October 2025) using the megablast option in Geneious Prime version 2023.0.1 (https://www.geneious.com) at a ≥ 97% similarity threshold and a minimum query coverage of ≥ 97%. Final identification in each database was based on the highest percent identity. If multiple hits were obtained with equal percent identity but with different species identification, we downgraded to genus-level identification (Arulandhu et al. 2017). Haplotype sequences without significant hits and those assigned to non-insect species or terrestrial insect species were excluded from subsequent analysis. When two or more OTUs share the same species identification, which may indicate cryptic species, we distinguish them by appending a number to the scientific name (e.g., Cladotanytarsus vanderwulpi 1, Cladotanytarsus vanderwulpi 2).

OTUs detected in at least five sites were used for population genetic analysis to ensure a minimum of five data points for the correlation analysis. More frequently detected and less frequently detected OTUs in the catchment were determined by the number of sites. We classified OTUs identified at ≥ 17 sites (75^th percentile) as more frequently detected and those at ≤ 6 sites (25^th percentile) as less frequently detected. Haplotype read counts were then transformed into presence/absence data. Haplotype richness was based on the total number of unique haplotypes of an OTU in a sample. Global ɸ-statistics (GɸST) among populations within an OTU were calculated using analysis of molecular variance (AMOVA) (Weir and Cockerham 1984; Excoffier et al. 1992; Weir 1996) according to differences in the allelic content of haplotypes (pairwise difference) within and among populations (Excoffier et al. 1992). Pairwise ɸST was also computed using the dissimilarity of allelic content of the haplotypes (pairwise difference) between populations. We also determined the average pairwise ɸST per OTU. Both GɸST and pairwise ɸST were computed using Arlequin version 3.5.2.2 (Excoffier and Lischer 2010).

To understand the relationship and location of haplotypes within an OTU in the catchment, we analyzed their phylogeography and haplotype network. We used PopART (Bandelt et al. 1999) to create the minimum spanning network. We classified haplotype network structures as “star-like,” “complex star,” “reciprocally monophyletic,” or “complex mutational” (Jenkins et al. 2018). Populations with few haplotypes and networks that did not fit these types were described as showing “no pattern.” We defined ancestral haplotypes as those at the network center and found more often in the catchment. For phylogeographic analysis, we mapped all haplotypes within an OTU using QGIS version 3.32.0-Lima (https://qgis.org/en/site).

We used the ggscatter function (ggpubr v. 0.6.0, Kassambara 2023) to perform a Kendall correlation test on the number of sites, ɣ-genetic diversity, average haplotype richness per site, GɸST, and average pairwise ɸST. We excluded OTUs with GɸST < –1 from correlation analyses employing GɸST because these values may indicate a statistical artifact resulting from the use of haplotype presence/absence data in AMOVA. Because we transformed haplotypes to equal frequencies within a population, the degrees of freedom for computing within-population genetic variation may be very low, potentially affecting the variance estimate, particularly for OTUs in which most sampling sites harbor a single haplotype. We created correlation plots with ggplot2 (Wickham 2016) and ggpubr. In the average pairwise ɸST analysis, we excluded OTUs that were genetically homogeneous (one haplotype per species). For GɸST, we excluded OTUs that were either genetically homogeneous or monomorphic (one haplotype per site).

α- and β-SGDCs were analyzed using presence/absence data for OTUs and haplotypes. Specifically, for the α-SGDC, haplotype richness was correlated with community OTU richness. OTU richness was measured as the number of unique OTUs in the community. Kendall correlation tests were performed using the ggscatter function of the ggpubr package (version 0.6.0; Kassambara 2023). For the β-SGDC, population pairwise ɸST between sites was correlated with the dissimilarity of the community composition (Sørensen dissimilarity) between the same sites where the OTU was detected using Mantel tests with the Kendall method and 9,999 permutations in vegan version 2.6-4 (Oksanen et al. 2022). For OTUs with a small number of occurrence sites, vegan version 2.6-4 (Oksanen et al. 2022) conducted an exact permutation test using all possible data arrangements (Berry et al. 2011). In all SGDC analyses, we included only sites in which the species of interest was detected. The p-values were adjusted using the Benjamini-Hochberg false discovery rate (FDR) correction.

We used female wing length as an estimate of insect dispersal ability. Female wing length was categorized from one to eight (1 = shortest or < 5 mm; 8 = longest or ≥ 50 mm) according to the DISPERSE database (Sarremejane et al. 2020). We included OTUs identified at the genus level, along with the female wing-length data, in DISPERSE in the analysis. However, for OTUs within Diptera, the data were based on subfamilies and tribes available in DISPERSE. The longest female wing-length category for each OTU was correlated with the number of sites, average α-genetic diversity, average pairwise ɸST, and GɸST using the Kendall method implemented in the ggscatter function of the ggpubr package version 0.6.0 (Kassambara 2023).

We used Mantel tests to analyze the correlation between geographic distance (river and straight-line distances) and pairwise ɸ-statistics for each OTU. River distance is defined as the waterline distance along river corridors. In contrast, straight-line distance is the air distance across the landscape. Mantel statistics were implemented in the vegan package version 2.6-4 (Oksanen et al. 2022) using the Kendall method, with 999 permutations. For OTUs with a small number of occurrence sites, vegan conducted an exact permutation test using all possible data arrangements (Berry et al. 2011). We used the same analysis to test the correlation between geographic distance and dissimilarity in the community composition among sites where the OTU was detected. We excluded genetically homogeneous populations from the analysis. The p-values were adjusted using the Benjamini-Hochberg FDR correction.

We correlated the mean pairwise differences of environmental parameters with the number of sites, ɣ-genetic diversity, average pairwise ɸST, and GɸST for all 59 OTUs selected for population genetic analysis using the Kendall method in the ggscatter function of ggpubr (v. 0.6.0; Kassambara 2023). We adjusted p-values with the Benjamini-Hochberg FDR correction. For each OTU, the mean pairwise differences for each environmental parameter were calculated from the environmental distance matrix using all sites where the OTU was detected. We excluded genetically homogeneous OTUs from the average pairwise ɸST analysis and both genetically homogeneous and monomorphic OTUs from the GɸST analysis.

We used the capscale function in the vegan package version 2.6-4 (Oksanen et al. 2022), with 999 permutations for univariate distance-based redundancy analysis (dbRDA), to examine the correlation of eight environmental parameters, namely, pH, conductivity, dissolved oxygen, NH₄-N, NO₃-N, NO₂-N, PO₄-P, and elevation, with population pairwise ɸST in each OTU. The significance of the capscale model was tested using the anova.cca function. The p-values were adjusted using the Benjamini-Hochberg FDR correction. Univariate dbRDA was used because the less frequently detected OTUs in the catchment had fewer sites than environmental parameters. However, for 16 OTUs with a substantially higher number of sites (≥ 16 sites) than environmental parameters (8 parameters), an additional analysis using ordiR2step with forward selection and 999 permutations in the vegan package version 2.6-4 (Oksanen et al. 2022) was conducted to verify the results. The same analysis was performed to correlate the environmental parameters with the dissimilarity in the community composition among sites where the OTU was detected.

The sequencing run yielded 6,932,659 demultiplexed reads. We obtained 586,062 reads, comprising 1,879 ASVs, after primer removal, quality filtering, denoising, merging, and chimera removal. After length and abundance filtering, the number of reads was reduced to 526,468, yielding 1,423 ASVs. None of the ASVs contained internal frameshifts, stop codons, or deletions according to MACSE analysis. After pseudogene filtering based on HMM scan analysis, 1,297 ASVs, hereafter referred to as haplotypes, were retained. Of these, 918 haplotypes had taxonomic identification at least at the family level. After excluding non-insect and terrestrial insect sequences, we obtained 889 haplotypes and clustered them into 198 OTUs. Six of 198 OTUs (3%) were verified to be a mixture of haplotypes from two or more species. We excluded them from downstream population genetic and community analyses because they may affect the biological interpretation of the results. Of 60 sites, 59 were deemed successful. Site 9 was excluded from the downstream analysis because of the limited number of reads that matched public and in-house databases.

We considered 59 OTUs detected at more than five sites in the catchment. These OTUs belong to seven orders, namely Coleoptera (four OTUs), Diptera (17 OTUs), Ephemeroptera (20 OTUs), Megaloptera (one OTU), Odonata (five OTUs), Plecoptera (four OTUs), and Trichoptera (eight OTUs). We found 15 more frequently detected and 19 less frequently detected OTUs in the catchment (Table 1). Haplotype richness ranged from one to 31. GɸST ranged from –0.954 to 0.537 after excluding OTUs with GɸST < –1, while average pairwise ɸST ranged from –0.716 to 1 (Table 1).

We found a positive correlation between the number of sites and haplotype richness using 59 OTUs detected at more than five sites (Fig. 2a) and all 194 OTUs (Suppl. material 2: fig. S1). A positive correlation was also observed between the number of sites and average haplotype richness per site using 59 OTUs with more than five sites (Fig. 2b) and 117 OTUs present in at least two sites (Suppl. material 2: fig. S1). However, we did not find a correlation between the number of sites and the average pairwise ɸST or GɸST, even after removing OTUs with GɸST < –1 (Fig. 2c, Suppl. material 2: fig. S1).

Haplotype network analysis revealed that 17 OTUs displayed star-like networks, 13 had complex stars, three had complex mutational patterns, one had a reciprocally monophyletic pattern, four had genetically homogeneous patterns, and 21 had no pattern (Fig. 3, Suppl. material 3; Suppl. material 4: table SS1). In addition to network structures, phylogeographic analysis revealed that 16 OTUs exhibited one high-frequency haplotype in the catchment, six OTUs showed two high-frequency haplotypes in the catchment, and 14 had an ancestral haplotype that was frequently detected in the catchment (Fig. 3, Suppl. material 3; Suppl. material 4: table SS1).

We found no significant α-SGDC in any of the 48 OTUs examined. However, one more frequently detected OTU in the catchment (Rhyacophila sp.) exhibited a statistically significant β-SGDC (R = 0.2639, p = 0.0378) (Suppl. material 4: table S3).

We found a positive correlation between the number of sites and the longest female wing-length category (Fig. 4). However, we did not find a significant correlation between the longest female wing-length category and haplotype richness (R = 0.23, p = 0.071), average pairwise ɸST (R = 0.044, p = 0.74), or GɸST (R = 0.031, p = 0.82), even after removing OTUs with GɸST < –1 (R = 0.01, p = 0.93) (Suppl. material 2: fig. S3).

We found no significant correlation between geographic distance, population genetic differentiation, and dissimilarity in community composition using Mantel tests in all OTUs examined (Suppl. material 4: table S2). Furthermore, our results also showed no significant correlation between the number of sites and the mean pairwise difference in environmental parameters. However, as the number of sites increased, we found that the mean pairwise difference in environmental parameters converged toward the mean of the whole dataset (Fig. 5). Except for NH₄-N, NO₂-N, and elevation, the mean pairwise differences in environmental parameters also showed a positive correlation with haplotype richness (Fig. 5). No correlation was found between the mean pairwise difference in all environmental parameters and the average pairwise ɸST (Suppl. material 2: fig. S4). The mean pairwise difference in pH showed a marginally significant correlation with GɸST (R = 0.28, p = 0.05) (Suppl. material 2: fig. S5). However, the relationship was no longer significant after excluding OTUs with GɸST < –1 from the analysis (Suppl. material 2: fig. S5). Because pH was the only one of 10 environmental parameters to show a marginal correlation with GɸST in the initial analysis, and this result was not robust after removing OTUs with extreme negative GɸST, the two analyses may consistently suggest a lack of correlation between mean pairwise differences in environmental variables and GɸST.

Univariate dbRDA indicated that conductivity explained variation in the genetic structure of OTU3 (Onychogomphus viridicostus) and OTU60 (Amphinemura bulla). PO₄-P accounted for the genetic variation of OTU44 (Hexatoma sp.), while pH captured the genetic variation of OTU69 (A. kawamurai) (Suppl. material 4: table S4). Similarly, ordiR2step showed that conductivity, PO₄-P, and pH explained the genetic variation in OTU7 (Afronurus sp.), OTU26 (Antocha bifida), and OTU13 (Mataeopsephus japonicus), respectively (Suppl. material 4: table S5). At the community level, univariate dbRDA showed that elevation explained the most variation in community composition (eight OTUs). NO₃-N accounted for the variation in the community composition of three OTUs. This trend was also confirmed by ordiR2step: elevation explained variation in most community compositions (12 of 16 OTUs, 75%), followed by NO₃-N (three of 16 OTUs, 18.75%) (Suppl. material 4: table S5). We did not observe a correlation between genetic structure or community composition and environmental parameters for OTUs less frequently detected in the catchment.

We found a positive correlation between the number of sites and haplotype richness (ɣ-genetic diversity). However, no correlation was observed between the number of sites and the average pairwise ɸST or GɸST. This partially supports our hypothesis (1) that species with higher occurrence in the catchment show greater β- and ɣ-genetic diversity. The positive correlation between species occurrence and ɣ-genetic diversity may reflect asynchronous extinction-colonization dynamics across habitats, consistent with metapopulation theory. These dynamics are often shaped by habitat quality and spatial heterogeneity (Sutcliffe et al. 1997; Kindvall 1999; Franzén and Nilsson 2010; Cosentino et al. 2012). Species with higher occurrence in the catchment may experience more spatially heterogeneous patterns of local extinction and recolonization, which may promote genetic structuring among habitats. Founder effects during colonization, genetic drift, and restricted gene flow may in turn increase β-genetic diversity among populations (García-Verdugo et al. 2010; Cosentino et al. 2012; Angst et al. 2022) and consequently lead to higher ɣ-genetic diversity at the catchment scale.

In addition to spatial asynchrony in extinction and colonization influenced by habitat heterogeneity, other mechanisms may also explain the observed patterns. Because previous studies demonstrated that species abundance is positively correlated with the proportion of sites it occupies within a region (Gaston et al. 2000; Ten Caten et al. 2022), we expect that species more frequently detected might support metapopulations with larger local population sizes in the catchment, diminishing the effects of genetic drift and increasing ɣ-genetic diversity.

However, we did not find a correlation between the number of sites and β-genetic diversity. These results did not support our claim in hypothesis (1) that species more frequently detected in the catchment exhibit high β-genetic diversity due to spatial asynchrony in extinction and colonization. One possible explanation is that, assuming species with stronger dispersal ability can expand their range within the catchment and occupy more sites, the effects of spatial asynchrony in extinction and colonization may have been masked by gene flow in those species occupying more sites. This is consistent with García-Verdugo et al. (2010), who found that gene flow can counteract the differentiation produced by extinction and colonization.

We found no significant α-SGDC in any of the 48 taxa examined, and only one taxon (Rhyacophila sp.) exhibited a significant positive β-SGDC, thus failing to support our hypothesis (4), which predicts that species with higher occurrence in the catchment will show a higher probability of α- and β-SGDC resulting from environmental heterogeneity. A plausible explanation is that different factors act at both levels of biodiversity. This is in agreement with previous studies, which reported null SGDC due to different mechanisms shaping species and genetic diversity (Taberlet et al. 2012; Seymour et al. 2016).

Although we did not find a correlation between the mean pairwise difference in environmental parameters and the number of sites, we observed that as the number of sites increased, the mean pairwise difference in environmental parameters converged toward the mean of the entire dataset. However, based on elevation data, the observed convergence pattern is generally consistent with null expectations, providing limited support for our hypothesis (2), which postulates that environmental conditions are more variable in species frequently detected in the catchment compared to those with rare occurrences. Nevertheless, these findings may suggest that species found across more sites tend to sample broader environmental heterogeneity in the catchment than those with restricted occurrence. In addition, we also found a positive correlation between mean pairwise differences in environmental variables and haplotype richness (ɣ-genetic diversity). Assuming that β-genetic diversity contributes more strongly to ɣ-genetic diversity than α-genetic diversity, we anticipate that an increase in ɣ-genetic diversity potentially reflects an increase in β-genetic diversity. Hence, we predict that β-genetic diversity will increase with environmental heterogeneity, likely driven by spatial asynchrony in extinction and colonization in metapopulations (Sutcliffe et al. 1997; Franzén and Nilsson 2010; Angst et al. 2022). Contrary to our expectation, the number of sites did not correlate with β-genetic diversity, suggesting that gene flow seems to offset the genetic variation arising from asynchronous extinction and colonization dynamics in the metapopulation. The observed absence of correlation between geographic distance and genetic differentiation may further reinforce our findings on gene flow, as previous studies have shown that aquatic insect populations often lack isolation-by-distance, in which gene flow overwhelms drift (Miller et al. 2002; Phillipsen et al. 2015).

Although we found a positive correlation between female wing length and the number of sites, our results show no correlation between female wing length and β-genetic diversity, contradicting our hypothesis (3), which proposed that β-genetic diversity decreases with species dispersal ability, estimated from female wing length, due to gene flow. Taken together, these results may imply that although a species with strong dispersal ability can colonize and recolonize more sites, gene flow may not be determined solely by species’ morphological dispersal ability (e.g., female wing length). It may be influenced by other species-specific factors, such as behavior and demographic and evolutionary histories. Previous studies have shown that some species of aquatic insects, despite strong dispersal abilities, may preferentially select nearby sites that are favorable for oviposition (Monaghan et al. 2001). Other species display a lack of genetic differentiation among drainages due to their past evolutionary history (e.g., historical gene flow; Monaghan et al. 2002). However, our inference on the relationship among female wing length, dispersal ability, gene flow, and β-genetic diversity is necessarily indirect, as none of these parameters were directly measured in this study. For a more robust interpretation of these relationships, future studies should gather and integrate morphological, dispersal, and genetic data.

At the community level, spatial structuring appears to be driven by species sorting, as the majority of species showed that the elevation gradient accounts for variation in community composition. Previous research indicated that elevation reduces the functional diversity of macroinvertebrates caused by enhanced environmental filtering in high-altitude areas (Zhang et al. 2025). The decoupling between the two levels of biodiversity, that is, gene flow potentially promoting genetic homogenization in metapopulations and species sorting possibly driving divergence in the metacommunity, may explain the absence of SGDC in this study.

Our results demonstrate the effectiveness of haplotype-level metabarcoding for analyzing multiple parallel populations. This approach enabled us to compare genetic diversity, population genetic structure, and SGDC across species inhabiting the same watershed. This interspecific comparison provides evidence that species-specific evolutionary traits are strongly associated with the extent of each species’ distribution. Because we can capture genetic signatures across many species with varying traits, behavior, and demographic histories at a large spatial scale, our approach is likely to benefit metapopulation and metacommunity research. This overcomes limitations of single-species analysis (Akçakaya et al. 2007) and may open opportunities for multiscale assessment of metapopulation and metacommunity (Malagon et al. 2025). Therefore, the application of haplotype-level metabarcoding is critical for biodiversity conservation.

However, there are limitations to haplotype-level metabarcoding that warrant further improvements in future research. First, there are concerns about potential biases arising from reliance on the number of sequence reads for estimating allele frequencies. This bias can be significantly influenced by PCR amplification, next-generation sequencing depth, and specimen size. In this study, we transformed haplotype sequence reads into presence/absence data to remove bias, and the only other source of variation is mutational differences among haplotypes (Azarian et al. 2020). However, to increase the likelihood of obtaining statistically significant population differentiation, future studies should consider using replicate samples (Azarian et al. 2020). Additionally, we recommend that future studies explore whether approaches such as quantitative MiSeq (qMiSeq) sequencing (Ushio et al. 2018), which can estimate the relative abundance of each species in species-level metabarcoding, can also be adapted to estimate allele frequencies at the haplotype level. In qMiSeq, by incorporating a known quantity of short DNA fragments (internal standard DNA) from haplotypes of species absent from the study area into all samples during PCR, the number of reads from non-standard DNA may be adjusted subsequently. However, estimating precise animal biomass and abundance based on the measured or calculated copy numbers of mitochondrial DNA (mtDNA) markers is still challenging due to the high variation in mtDNA copy numbers per input biomass (Wolany et al. 2026).

In future studies, environmental DNA (eDNA) may also be used for eDNA haplotyping. This approach allows less invasive, time-saving, and cost-effective sampling. It enables the collection of a greater number of samples, reduces the effort required for sample replication, and allows a broader comparison of intraspecific variation across taxa and regions. However, reliable data require careful attention to the challenges of environmental samples, including eDNA degradation, transport, and ecology, as well as PCR inhibition (Jane et al. 2015; Yoon et al. 2025).