Hong Kong SAR, China is densely populated with ~7 million inhabitants, peaking at >57,000 people per km², and averaging ~7,000 individuals per km² across the territory. While urban density is exceptionally high, approximately 70% of the territory is undeveloped, with 40% currently designated as country parks. The 30% of the territory that is developed includes a mosaic of land use that has shifted considerably over the years and reflects a tiered intensity of human impact. More specifically, 22% of the Hong Kong SAR is classified as high-density urban cores, 3% as rural, low population density settlements, and 5% as fragmented agricultural land. The unique land use structure offers a sharp natural to urbanization gradient within a compact geographic area.

We sampled from 30 Hong Kong streams to assess the biodiversity of mayflies across the natural-to-urbanized gradient (Fig. 1). Three urbanization categories were defined. 1) Natural (low disturbance): located within protected country parks; characterized by intact riparian buffers (woodland/shrubland/grassland) and minimal impervious cover. 2) Peri-urban (moderate disturbance): sites with low-density residential or agricultural matrices while retaining a high degree of ecological environments and connectivity; subject to moderate nutrient loading and hydrological alteration. 3) Urbanized (high disturbance): sites exhibiting urban stream syndrome (Walsh et al. 2005), characterized by extensive impervious surfaces, channelization, and flashy hydrographs. Sites were selected to encompass the full natural-to-urbanized land-use gradient across Hong Kong while maximizing geographic representation and accessibility.

Mayfly specimens were collected from 22 sites in September 2022 (10 natural, 8 peri-urban, 5 urban) and 30 sites in November 2023, including the same 22 sites sampled in September 2022, along with 8 additional sites (3 natural, 4 peri-urban and 1 urban) (Fig. 1). We prioritized head water sites that were accessible while ensuring we encompassed the full natural-to-urbanized land-use gradient across Hong Kong. At each sampling event, mayfly nymphs were collected using a standardized 3-minute kick sampling protocol, with a D-frame kick net (500 μm mesh) across all microhabitats present (Furse et al. 1981; Seymour et al. 2021). Bulk specimens were preserved in 95% ethanol and counted and identified to morphospecies at the University of Hong Kong’s eDNA & eEcology Lab.

DNA was extracted from 2–5 representative specimens per morphospecies using the Universal Genomic DNA kit (CW2298, CoWin Biotech, China), following the manufacturer’s protocol. For each specimen, 4–6 legs or the thorax (for smaller specimens) were used. We aimed to generate 3–5 barcodes for all morphospecies identified, but were restricted to 1 or 2 barcodes for 3 groups due to limited abundances or poor quality DNA extraction (see Suppl. material 1: table SS2 for a complete breakdown of specimens and barcodes generated). The 658 bp cytochrome c oxidase subunit I (COI) barcode region was amplified using LCO1490/HCO2198 primers (Folmer et al. 1994). If samples failed to amplify with the standard primers, we used our custom primers (MayF1/R1; MayF2/R2; Table 1), which were designed for this study by modifying the standard primers to bind to other target groups (e.g., Han et al. 2026a). The MayF1/MayR1 and MayF2/MayR2 primers increased overall success by 21.4%, enabling amplification in previously recalcitrant taxa (Afronurus alces: 100%, n = 4; Bungona longisetosus: 40%, n = 10; Cloeon navasi: 43%, n = 7) and improving sequence quality for samples that generated poor standard sequences (19.8% cases with incomplete traces using standard primers while zero using modified primers). All PCR reactions (total volume = 25 µl) contained 12.5 µl 2 × Hieff PCR Master Mix (YEASEN, China), 6.5 µl water, 0.5 µl of each primer (10 µM), and 5 µl of template. Thermal cycling conditions: 94 °C for 2 min; 35–41 cycles of 94 °C for 30 s, 49–50 °C for 30 s, 72 °C for 50 s; 72 °C for 10 min (ProFlex PCR System, Applied Biosystems, USA). PCR cycles were only increased after failed amplification at 35 cycles as a means to account for low concentration DNA samples without adversely affecting amplicon quality (Mitchell 2015; Vierna et al. 2017). Products were visualized on 1.2% agarose gels and sequenced bidirectionally using Sanger sequencing at GeneWiz Inc. (Suzhou, China).

We used an integrative morphology and molecular approach to ensure robust species delimitation (Zhou et al. 2009; Cardoni et al. 2015; Han et al. 2026a). Mayfly nymphs were initially sorted and identified to the lowest possible taxonomic level using a stereo microscope and regional taxonomic keys (Tong 2001). To resolve challenging groups, particularly within the families Baetidae, Heptageniidae, and the genera Choroterpes and Habrophlebiodes, subsets of specimens were dissected to determine taxonomic differences. Dissected mouthparts, antennae, legs, gills, and paraprocts were mounted on slides using Hoyer’s medium (Phygene Biotech, China) and examined under high magnification (Flexacam i5, Leica, Germany). For specimens that could not be confidently assigned to a known species due to a lack of descriptive literature or diagnostic features, interim morphospecies codes were assigned to allow for downstream molecular assessment of putative species. Representative specimens of each morphospecies were photographed dorsally, ventrally, and laterally using a M205C stereo microscope (Leica, Germany).

COI barcode sequences were edited in BioEdit v7.2.5 (Hall 1999) to generate consensus sequences, which were subsequently aligned using MUSCLE (Edgar 2004). Barcode sequences and associated specimen photos and trace files were then uploaded to the Barcode of Life Data Systems (BOLD) (dataset DS-HKMAYFLY) (Ratnasingham and Hebert 2007) for primary molecular grouping and Barcode Index Number (BIN) assignment.

To complement the BOLD framework and cross-check our morphological assignments against a broader public repository, sequences were also queried against the GenBank (NCBI) database via BLAST (Camacho et al. 2009) to determine taxonomic proximity. In this workflow, a conservative 95% similarity BLAST threshold was utilized as a diagnostic baseline to flag broad taxonomic alignment (e.g., confirming order- or genus-level placement) and to screen for potential database gaps or contamination (Padial et al. 2010). Six potential outcomes were addressed:

1) Concordance: where molecular matches (>95% similarity to a reference sequence) were consistent with morphological identification, the species name was confirmed. 2) Morphological re-evaluation (intra-ordinal conflict): if a high-similarity molecular match (>95%) conflicted with the initial morphological ID but remained within Ephemeroptera, the voucher specimens were re-examined. The molecular ID was only accepted if morphological features could be re-interpreted to support it. 3) Contamination detection (extra-ordinal conflict): if a high-similarity match identified the specimen as a non-target organism (out of the order Ephemeroptera), the sequence was flagged as contamination and removed from the dataset. 4) Database gaps (low similarity, low divergence): if similarity to reference databases was <95% but the specimens within the morphospecies group showed zero or negligible genetic divergence, the morphological identification was retained, and the record was treated as a new barcode for a known species or a regional endemic not yet represented in global databases. 5) Cryptic diversity (low similarity, high divergence): if the top-hit match was <95% and high genetic divergence was observed within a single morphospecies group, we hypothesized the presence of cryptic lineages. These were assigned unique interim names and further validated via neighbor-joining (NJ) tree clustering. 6) Reciprocal refinement (systematic ambiguity): in instances where morphological diagnosis was constrained by specimen condition (e.g., underdeveloped diagnostic features) and molecular database queries yielded low similarity (<95%), taxonomic assignments were resolved through internal sequence comparison. By analyzing the genetic cohesion within our own sample clades, we were able to either: (a) confirm the identity of ambiguous specimens, or (b) correct initial “morphospecies” misassignments through the detection of zero or negligible internal divergence.

In cases where sequences failed to amplify after multiple attempts, the identification defaulted to the morphological assignment. This integrated approach ensures that our taxonomic library provides a curated system based on concordant evidence.

To validate the efficacy of the COI library for species delimitation, we employed four complementary analytical approaches. First, sequences were assigned barcode index numbers (BINs) within the BOLD System using the RESL (Refined Single Linkage) algorithm (Ratnasingham and Hebert 2013). The RESL algorithm uses an initial 2.2% divergence threshold as a heuristic baseline to generate primary sequence clusters. The RESL then refines these clusters into Operational Taxonomic Units (OTUs) through subsequent graph theoretic analyses, including Markov clustering and silhouette score evaluations, rather than relying on a rigid, universal species boundary (Ratnasingham and Hebert 2013). Second, genetic distances were calculated using the Kimura-2-Parameter (K2P) model (Kimura 1980) in MEGA 12 (Kumar et al. 2024). We quantified mean and maximum pairwise intraspecific genetic divergences, as well as mean and minimum pairwise interspecific genetic divergences (distance to the nearest neighbor). These metrics were used to assess the presence of a barcode gap for species delimitation, defined as the separation between the maximum internal variation of a species and the minimum distance to its closest relative. Third, to visualize the concordance between molecular clusters and morphological identifications, a Neighbor-Joining (NJ) tree was constructed using the K2P model with 1,000 bootstrap replicates (Saitou and Nei 1987). The NJ tree included the 146 assigned barcodes. Among the 146 BIN-assigned records, three sequences were shorter than 500 bp (Suppl. material 1: table SS2) and therefore do not meet the BOLD barcode compliance criteria. However, we included these short barcodes in NJ tree analyses because they exhibited high sequence quality, with no ambiguous bases or evidence of contamination. The short sequences also perfectly aligned with, and showed genetic congruence to, their complimentary longer compliant barcode sequences assigned to the same BINs. The NJ tree showed that the short sequences clustered tightly with their corresponding BINs, supporting their taxonomic placement and analytical value, therefore verifying that the shorter sequences did not introduce noise or inconsistencies into the analysis. The NJ tree was visualized and annotated using iTOL v7 (Letunic and Bork 2024), providing a means to evaluate species monophyly and assess potential instances of paraphyly or polyphyly, which serves as a diagnostic check for specimens that may need re-examining or for potential contaminant screening. Lastly, a conservative 2.2% standard divergence threshold was used as a benchmark for potential cryptic diversity or taxonomic revision in our dataset. Arthropoda rarely exceeds a 2.2% divergence threshold, hence the chosen cutoff, unless there are potential cryptic lineages or a flawed taxonomic classification system (Hebert et al. 2003).

To objectively characterize the natural-to-urbanized gradient, site classification was determined by land use composition analysis within a 100 m radius of each sampling reach. A 100 m radius was chosen because it is wide enough to capture surrounding variability linked to local surface runoff, while also allowing enough space to prevent overlapping data between nearby sites (Morley and Karr 2002; Seymour et al. 2016). The land use classification analysis utilized the 2023 Land Utilization Raster Grids dataset provided by the Hong Kong Planning Department (Planning Department 2024). Land use classification was based on the proportional composition of land-cover types within a 100 m buffer surrounding each sampling site. Percent cover values for woodland, shrubland, grassland, agriculture, rural settlement, and urbanized surfaces were first verified to sum to 100% per site. Sites were then assigned to one of three categories (natural, peri-urban, urbanized) using a rule-based decision framework. Specifically, sites with ≥25% urbanized cover were classified as urbanized (Hill et al. 2000). Sites with ≥80% combined natural cover (woodland, shrubland, and grassland) were classified as natural. Remaining sites were classified as peri-urban if they contained ≥10% combined agricultural and rural settlement land use or >5% urban cover. Sites not meeting any threshold criteria were conservatively assigned to the peri-urban category, reflecting fragmented or mixed land use conditions (Ahogle et al. 2025).

At each sampling site, substrate composition (percentage of silt, sand, boulder, and gravel) was quantified (Seymour et al. 2021) during the November 2023 sampling period. Physicochemical variables, including acidity (pH), salinity (conductivity), dissolved oxygen (DO), and total dissolved solids (TDS) were measured in situ using a multiparameter meter probe (U-50 HORIBA, Japan). Elevation was recorded at each site using GPS. Land cover classification included natural, peri-urban, and urban, as defined above.

All statistical analyses were conducted in R v4.5.1 (R Core Team 2025). Gamma diversity was defined as the total number of taxa recorded across all 30 sites. Site-level diversity was assessed using species richness, calculated as the number of taxa recorded per sampling event. To evaluate inventory completeness and estimate regional species richness, species accumulation curves (SAC) and the Abundance-based Coverage Estimator (ACE) were calculated using the ‘vegan’ package (Oksanen et al. 2025). Total regional richness was estimated by pooling abundance data across all sampling events using the ACE metric. To visualize the rate of species discovery as a function of sampling effort, a sample-based species accumulation curve was generated using the random method with 1000 permutations (Suppl. material 1: fig. S1).

Given the overdispersion in the data, environmental effects on richness were tested using a negative binomial generalized linear model (Zuur et al. 2009). Models were fitted using the R package 'glmmTMB' (McGillycuddy et al. 2025). Species richness was initially modeled as a function of scaled elevation, urbanization intensity, peri-urban intensity, natural intensity, pH, dissolved oxygen, conductivity, total dissolved solids, gravel coverage, sand coverage, silt coverage, and boulder coverage. Predictor variables were centered and standardized prior to analysis. Collinearity was assessed using the 'performance' package (Lüdecke et al. 2021), with a standard threshold of Variance Inflation Factor (VIF) > 5, resulting in total dissolved solids, peri-urban intensity, boulder coverage, and silt coverage being removed. We first assessed whether to account for seasonal effect of sampling by including site as a random effect using a negative binomial mixed effects model, however, the site variance was near zero (σ² < 0.001), indicating that the random effect was not improving the model fit and caused a singular fit, so we reverted to the more parsimonious negative binomial glm. We used a backward model selection process to select the most parsimonious final model (Borcard et al. 2011). Model significance was evaluated using Wald z-tests, and model fit was assessed using the Akaike Information Criterion (AIC).

Indicator taxa associated with each urbanization category were identified using the Indicator Value (IndVal) approach via the R package 'indicspecies' (Cáceres and Legendre 2009). IndVal quantifies the strength of association between individual taxa and predefined site groups by combining two complementary components: specificity (the probability that a surveyed site belongs to the target group given the presence of the taxon) and fidelity (the probability of finding the taxon within sites of that group). The product of these two components yields the indicator value for each taxon, with high values reflecting taxa that occur frequently and predominantly within a single habitat category. Indicator analysis was performed for the three urbanization classes, natural, peri-urban, and urbanized, to detect taxa that were significantly associated with each class. Statistical significance was evaluated by permutation testing at α = 0.05. Significant taxa were considered habitat indicators and were used to characterize assemblage differences across the urbanization gradient.

A total of 12,155 mayfly nymphs were collected from 30 stream sites during September 2022 and November 2023. Initial taxonomic analyses resulted in 31 morpho-taxonomic groups, which formed the baseline barcoding focus. Subsequently, integrated taxonomic analysis identified 45 morphospecies across 28 genera and 6 families (Table 2). Of the 207 specimens selected for DNA extraction, 154 COI sequences were obtained, yielding an average of 3.85 sequences per morphospecies (range: 1–12). DNA amplification was unsuccessful for three taxa: Nigrobaetis facetus, Thraulus bishopi, and Baetis cf. iIlliesi. While we attempted to re-extract and sequence these groups, we were unsuccessful due to limited available material.

Of the total barcode sequences, 146 (94.8%) were successfully assigned to 53 BINs. The 10 remaining sequences failed assignment due to high proportions of ambiguous bases (>1% N-count; n = 8) or represented highly divergent singletons awaiting further taxonomic validation (n = 2). Overall, this first regional DNA barcode library for mayflies in China includes 146 barcodes (mean length: 626 bp; range: 409–658 bp; Suppl. material 1: table SS2) from our 45 integrative taxonomy morphospecies groups. Of these, 39 were successfully barcoded (32 BOLD confirmed species), and clustered into 53 unique BIN assignments (Table 2). In total, 5 of the 45 morphospecies groups are represented by a single barcode (singleton). Notably, 34 of the 53 BINs were unique records in the Barcode of Life Data System (BOLD). The database is publicly available, with digital micrographs, geographic metadata, and sequence trace files uploaded to BOLD to facilitate future taxonomic revisions (dataset DS-HKMAYFLY).

Comparisons with public records in GenBank revealed that 38 sequences (26.4%) yielded a high-confidence molecular match with globally available sequences (>95% similarity). Approximately 41% of the sequences (59 sequences across 21 species) represent taxa previously absent from global databases.

A species accumulation curve indicated that sampling effort was sufficient and approached an asymptote (Suppl. material 1: fig. S1). The Abundance-based Coverage Estimator (ACE) predicted a regional species pool of approximately 52 species, compared to the observed gamma diversity of 45 morphospecies recorded across the 30 sampling sites. This suggests that our survey captured approximately 86.5% of the estimated total regional richness.

Rigorous taxonomic re-evaluation resulted in four new records and seven putative new lineages for Hong Kong: Arcobaetis cf. gracilentus, Baetidae sp. HK10, Baetis cf. maculosus, Caenis sp. HK06, Electrogena sp. HK08, Procloeon sp. HK04, Procloeon sp. HK05, Procloeon cf. chinensis, Procloeon cf. strumiferum and Tenuibaetis sp. HK09. These findings expand the existing regional checklist by 14%, with a 65% species overlap (Suppl. material 1: table S3).

The mean genetic divergence across the dataset was 25.4%. A distinct barcoding gap was observed: the maximum intraspecific genetic divergence averaged 1.4% (range: 0–9.4%), while the minimum interspecific genetic divergence averaged 22.4% (range: 12.2–34.8%). No barcode sharing was detected between species, and no overlap occurred between maximum intraspecific genetic distances and distances to the nearest neighbor. Neighbor-joining (NJ) analysis (partial deletion, 90% coverage, 572 bp), which is not meant to be a definitive phylogenetic assessment, recovered monophyletic clusters for all morphospecies groups with high (>99%) bootstrap support. While the families Baetidae, Ephemerellidae, and Heptageniidae indicated monophyletic groupings, Ephemeridae, Leptophlebiidae, and Caenidae suggested large paraphyletic clustering. Within the species-rich Baetidae (23 species), only the genera Labiobaetis, Bungona, and Baetis did not form monophyletic clusters (Fig. 2).

Overall, 89% of species groups were successfully delimited, with three taxa exceeding the 2.2% baseline threshold: Bungona fusina (2.8%), Paegniodes cupulatus (4.7%), and Isca purpurea (8.2%). Excluding BINs with elevated threshold values (i.e. BINs: BOLD:AGW9771, BOLD:AGX7348, BOLD:AGV1965), divergence values ranged between 1.0–2.3%. Subsequent morphological re-evaluation of these voucher specimens revealed no consistent diagnostic characteristics between these deeply divergent lineages.

Species richness declined with increasing urban coverage, peaking at 26 species and declining to 0 (mean = 9) (Figs 3, 4). The final negative binomial GLM results showed that richness decreased significantly with increasing urbanization (β = -0.308, p = 0.007) and conductivity (β = -0.427, p = 0.002), while elevation had a marginally non-significant negative effect (β = -0.210, p = 0.051) (Fig. 4).

Indicator species analysis identified three taxa significantly associated with habitat categories. One taxon (Choroterpes elliptica HK01), a species not officially described in the literature, was associated with natural sites, two taxa (Tenuibaetis pseudofrequentus and Caenis lubrica) were associated with peri-urban sites, and no species were indicative of urbanization. Overall, the indicator analysis highlighted a diverse assemblage of mayfly taxa in peri-urban sites suggesting that these habitats serve as an important buffer between natural and urbanized sites (where biodiversity declined significantly) (Figs 4, 5).