The southern Great Escarpment forms an 800-km-long stretch of mountain complexes extending from the Nuweveldberge in the west to the Eastern Cape Drakensberg in the east. Ancient erosional features divide the mountains into five main blocks that range in altitude from 1 600–3 000 m a.s.l., making the area interesting for study of dispersal-limited groups.

Nymphs of A. sudafricanum, D. capensis, and D. natalensis were collected from 21 rivers in the Eastern Cape Great Escarpment, relating to 12 study areas within the Escarpment and non-Escarpment sites (Table 1). An additional six rivers were sampled for A. sudafricanum in lower-altitude (non-Escarpment) areas in the Eastern Cape and KwaZulu-Natal (Table 1). All specimens were preserved in 80% ethanol.

A related species of Baetidae, Baetis rhodani Pictet, was used as the outgroup for phylogenetic analyses, and relevant sequence data (Rutschmann et al. 2014) were obtained through Genbank (Benson et al. 2012) for both cytochrome c oxidase subunit I (COI) (KP438135 and KP438160) and 16S rRNA (16S) (KP438109 and KP438119) gene regions.

DNA was extracted using the Invisorb Spin Tissue Mini Kit following manufacturer’s protocol (Invitek, Berlin, Germany). Extraction was non-destructive, using internal body digestion, which ensured the preservation of the exoskeleton for future morphological analysis (housed in the Albany Museum, Makhanda, South Africa, along with additional material that is stored in the collection, listed under the GEN catalogue.)

Two mitochondrial gene regions were amplified: cytochrome c oxidase subunit I (COI) and small subunit ribosomal RNA (16S). A 528-bp section of the COI regions of D. natalensis and D. capensis was successfully amplified with the standard ‘universal’ primer pair, LCO1490 and HCO2198 (Folmer et al. 1994), which worked with only limited initial success with A. sudafricanum. A new forward primer (5’–GGT GGA TGG GCA GGA ATG GTA GGA–3’) was designed and used with HCO2198 to successfully sequence the remaining samples of A. sudafricanum. The 16S region was amplified with the primer pair 16Sar (5’–CGC CTG TTT ATC AAA AAC AT–3’) and 16Sbr (5’–CCG GTC TGA ACT CAG ATC ACG T–3’) (Palumbi 1996). However, these primers proved problematic for the Demoreptus samples, and this region is thus excluded from subsequent analyses for this taxon.

The polymerase chain reaction (PCR) was performed in a 50 μl volume using the following thermal regime: 95 °C for 5 min, 35 cycles of 95 °C for 45 s, 50 °C for 45 s, and 72 °C for 90 s, followed by a final extension period of 72 °C for 5 min. PCR amplifications were checked for the presence of amplified PCR products by gel electrophoresis (0.5% agarose gel stained with SYBR green) and viewed with a UV-transilluminator. Successful PCR products were cleaned up using the Invisorb PCRapace® Quick purification kit (Invitek, Berlin, Germany) and cycle-sequenced in both directions using the primers used for amplification, the ABI Big Dye Sequencing kit v.3.1 (following manufacturer’s instructions (Applied Biosystems)), and a ABI Genetic Analyzer 3500 (Applied Biosystems).

Sequence trace files were assembled and edited using Sequencher 3.0 (DNA sequence analysis software, Gene Code Corporations, Ann Arbor, MI USA, http://www.genecodes.com). The sequences were then aligned in MEGA v.6 (Tamura et al. 2013) using the ClustalW algorithm and subsequently each non-synonymous mutation was manually cross-checked in the original trace files.

Each gene was tested for substitution saturation using plots of transitions and transversions against F84 distance in DAMBE v7.0.58 (Xia et al. 2003, Xia and Lemey 2009, Xia 2017).

Congruence between the COI and 16S datasets was assessed using the partition homogeneity test (PHT) in PAUP* (Swofford 2002) with 1000 replicates to verify that the gene sections could be combined for analysis.

Bayesian Inference (BI) analyses were conducted with MrBayes v.3.1.2 (Huelsenbeck and Ronquist 2001) using the GTR+I+G model since it is the most complex model, allowing the nesting of simpler models that could be estimated through the Bayesian sampling. Each analysis comprised two independent runs with random starting trees and four chains (three heated and one cold) each, sampled every 200 generations for 20 million generations per run. The cumulative sample sizes were plotted against the likelihood scores and tree length using Tracer v1.7.0 (Rambaut et al. 2018), to ascertain when the analysis reached stationarity after the first 10% of the trees were discarded as burn-in. The analysis was run on the CIPRES Science Gateway (Miller et al. 2010) using default parameters for variables not mentioned above.

Maximum likelihood (ML) analyses were conducted with 2 000 bootstrap replicates using the GARLI (Genetic Algorithm for Rapid Likelihood Inference) on XSEDE via the CIPRES (Cyberinfrustructure for Phylogenetic Research) Science Gateway v3.3 (Miller et al. 2010), which is supported by the San Diego Supercomputer Center (SDSC) and the University of California (UC San Diego). Models of molecular evolution for each dataset were selected using the Akaike information criteria (AIC) as implemented by jModeltest 2.1.6 (Darriba et al. 2012) (Table 2). The COI and COI+16S ML phylograms were compared and presented using Phylo.io software (Robinson et al. 2016).

Parsimony analyses were performed in PAUP* version 4.0b10 (Swofford 2002) using the heuristic search option with 100 random addition replicates. A search with TBR (Tree Bisection and Reconnection) branch-swapping was used to find the approximate length of the shortest trees, with one tree kept with each random addition. To investigate nodal support, all of the trees from this search were used as starting trees for a second heuristic search with MAXTREES set to 5 000. The Demoreptus analysis used FASTBOOTSTRAP with 10 000 replicates.

Molecular diversity was investigated using the COI datasets. The number of variable sites (S), number of haplotypes (Hap) and haplotype diversity (Hd), Nucleotide diversity (p) and neutrality tests (Tajima’s D and Fu’s F_S) were calculated in DNAsp (Rozas et al. 2017). Population structures within each species were estimated using one-level Analyses of Molecular Variance (AMOVA) in ARLEQUIN ver. 3.5.2 (Excoffier and Lischer 2010). F_ST (fixation index) values were calculated between localities to determine whether putatively conspecific populations differed significantly in their genetic composition. For all AMOVA analyses (listed in Table 4), a priori groups were defined by each site where the insects were collected. Haplotype networks were illustrated with a median-joining network (MJN) algorithm (ε = 0) (Bandelt et al. 1999) using the software PopART v. 1.7 (Leigh and Bryant 2015) to analyse haplotype genealogy.

COI sequences (649 bp) were obtained from 86 individuals and 16S sequences (542 bp) obtained from 59 individuals of A. sudafricanum. Shorter (528 bp) COI sequences were obtained for 24 Demoreptus individuals (for D. natalensis, N = 12; D. capensis, N = 11; unidentified, N = 1). DNA characteristics for each gene dataset are summarised in Table 2. The Partition Homogeneity Test for incongruency (Swofford 2002) showed that the combined COI and 16S datasets were not significantly incongruent (P = 0.3890) and could therefore be combined for analysis. The COI+16S molecular dataset consisted of 59 specimens and 1191 nucleotides including the outgroup. No evidence of saturated substitution was found for either gene (data not shown).

The parsimony analyses’ results are summarised in Table 2. Phylogenetic analyses of the habitat generalist A. sudafricanum consistently retrieved six distinct clades and an unresolved grade of specimens (referred to hereafter as the “widespread grade”) from a wide range of sites for all analyses (Fig. 1). The tree comparison shows that the relationships between these clades in the analysis of the CO1 data and the CO1+16S data sets were consistent, with improved support for the combined dataset. The well-supported clades did not conform to the mountain blocks described in Table 1 and included specimens from across these ranges. The clades were roughly separated into overlapping geographic groups: Southern Montane, Stormberg/Barkly East, and Eastern Cape, while more restricted geographic areas included KwaZulu-Natal and Eastern Cape Drakensberg, the latter clade showing a longer stem branch compared to other clades (Fig. 1).

Figure 1. 

Bayesian inference phylograms of A. sudafricanum for gene markers COI (left) and COI + 16S (right). Support for major nodes shown in the order Bayesian Inference / Maximum Parsimony / Maximum Likelihood (BI/ML/MP). Bars next to clades refer to distinct clades that are colour-coded according to the study areas found within that clade (see colour legend), except for the widespread grade which is designated by a solid black line. Branches bearing outgroups have been omitted to save space and their position is depicted by a dashed line.

The phylogenetic analysis of the habitat-restricted D. capensis and D. natalensis clearly indicated strong genetic structure corresponding to geographic location (Fig. 2). Both species had genetically distinct populations with strong support from parsimony, Bayesian and maximum likelihood analyses. The clades found for D. capensis and D. natalensis were more closely aligned with the mountain ranges described in Table 1 and appear to be site-restricted, apart from one instance where individuals of D. natalensis from Rhodes and Barkly Pass fell into the same well-supported clade. For D. capensis the Rhodes clade had a long branch and is clearly distinct from the other clades, which is noteworthy considering the close geographic proximity to Barkly East, which formed a separate well-supported clade nested with other clades (Fig. 2). This pattern was not apparent in D. natalensis, where the Barkly East specimen was clearly separate from the D. natalensis clade; morphological re-examination suggests that it does not belong to any described Demoreptus species.

Figure 2. 

Bayesian inference phylogram of Demoreptus spp for the COI gene marker. Support for major nodes is shown in the order Bayesian Inference / Maximum Parsimony / Maximum Likelihood (BI/ML/MP). Bars next to clades refer to distinct clades that are colour-coded according to the study areas found within that clade (see colour legend). Baetis rhodani Pictet was used as the outgroup.

MJN analysis collapsed the 86 A. sudafricanum COI sequences into 60 haplotypes (Fig. 3, Table 3), 45 of which were singletons or private haplotypes. Haplotype 17 was the most abundant (N = 8) and occurred in three of the 12 study areas (Fig. 4), which included non-Escarpment Grahamstown and two main mountain Escarpment blocks (Sneeuberg and Winterberg–Amathole). Haplotype 20 was next-most-abundant (N = 4) and exclusive to non-Escarpment Makhanda (= Grahamstown). Haplotype 10 (N = 3) was found in one non-Escarpment site (KwaZulu-Natal) and one main mountain Escarpment block (Eastern Cape Drakensberg, in two study areas: Barkly Pass and Maclear; Figs 3, 4). Haplotypes were clustered according to a broad geographical structure, which correspond to clades from the phylogenetic analyses (Fig. 1). The numerous missing mutational steps in the haplotype network (Fig. 3) suggest that more sampling is needed for some clusters, particularly between sites that are separated by long sampling gaps (for example, the non-Escarpment sites). Other clusters that are separated by numerous missing intermediates could represent cryptic species or relict lineages that have re-joined the metapopulation (Hinojosa et al. 2019) (encircled with dashes in Fig. 3). The divergent Haplotype 27 from the Kamdebooberg did not cluster with the other haplotypes from the same area and may represent such a lineage. The widespread grade showed little geographic structure, and all haplotypes from Stormberg (Hap 51, 52 and 53) grouped together exclusively, otherwise all other sites are mixed.

Figure 3. 

Median-joining network of A. sudafricanum based on COI haplotypes generated in this study. The network was estimated using the median-joining algorithm in PoPArt v.1.7 with epsilon = 0. Each circle represents a different haplotype and the size of a circle correlates with the number of individuals assigned to that haplotype. Only haplotypes found in more than one sample are numbered. Colours indicate the geographic origin of sequences; black dots indicate unsampled or extinct haplotypes.

The MJN analyses for D. capensis retrieved eight haplotypes (Hd = 0.9273, S = 101), six of which were singletons and D. natalensis retrieved six haplotypes (Hd = 0.8636, S = 21) including three singletons (Fig. 5). Haplotypes were largely site-restricted for both species with the exception of Haplotype 1 (N = 3) for D. capensis and Haplotype 4 (N = 3) for D. natalensis (study areas Barkly Pass and Rhodes), which are both found in the Eastern Cape Drakensberg main mountain block in the Escarpment (Fig. 4). Both networks show many missing mutational steps between haplotypes grouped by locality, which could result from undersampling or haplotype filtering.

Figure 4. 

Distribution of A. sudafricanum, D. natalensis and D. capensis COI haplotypes across the study area. The map shows the study areas defined in Table 1, and the pie charts indicate the haplotype composition of the population from each area. Each colour represents a shared haplotype found across the study area; private haplotypes (singletons found in the samples from one particular population and are absent in the samples from other populations) are represented as clear sections within the pie charts.

Nucleotide diversities (P_i) are reported in Table 3 and are not interpreted further because the small sample sizes for Demoreptus spp. make the estimates imprecise. Neutrality tests (Tajima’s D and Fu’s F_S) were not significant for A. sudafricanum, D. capensis or D. natalensis indicating that the nucleotide patterns of variation are consistent with the neutral theory of evolution. Fu’s F_S statistic for the widespread grade of A. sudafricanum was negative (F_S = –11.544) and significant (P < 0.02), indicating a recent population expansion (Table 3). The Fu’s F_S statistics for D. capensis and D. natalensis were positive, indicating a deficiency of alleles as expected from a population bottleneck, but they were not significant and need a larger sample size to confirm these results.

The AMOVA results for A. sudafricanum revealed that over all localities, 52.33% of the total variance was explained by variation among populations (df = 10, Va = 12.073) while 47.67% (df = 75, Vb = 10.998) was explained by variation within populations (Table 4). A similar result was found with the widespread grade of A. sudafricanum, with 39.43% of the total variance explained by variation among populations (df = 5, Va = 2.238) and 60.57% (df = 28, Va = 3.438) explained by variation within populations. In contrast, the AMOVAs for the habitat-restricted species, D. capensis and D. natalensis, indicated a higher proportion of variation among populations: 94.83% (df = 4, Va = 24.950) and 95.39% (df = 4, Va = 5.423), respectively (Table 4). The total variance explained by variation within populations was only 5.17% (df = 6, Vb = 1.361) for D. capensis and 4.61% (df = 7, Vb = 0.262) for D. natalensis.

Figure 5. 

Median-joining networks of D. natalensis and D. capensis based on COI haplotypes generated in this study. The network was estimated using the median-joining algorithm in PoPArt v.1.7 with epsilon = 0. Each circle represents a different haplotype and the size of a circle correlates with number of individuals belonging to that given haplotype. Only haplotypes found in more than one sample are numbered. Colours indicate the geographic origin of sequences; black dots indicate unsampled or extinct haplotypes.

The measure of population differentiation due to genetic structure (F_ST) was much lower for A. sudafricanum compared to the Demoreptus species (Table 4). The widespread grade for A. sudafricanum had a very low F_ST value of 0.39 while D. natalensis and D. capensis had very high F_ST values of over 0.94 (Table 4).

These results help to illuminate some of the evolutionary processes occurring in mayfly species and highlight the effect of habitat-specificity on haplotype diversity and partitioning within a species. While all three species have qualitatively similar levels of dispersal potential in terms of flight, they show differences in gene flow, suggesting that other processes, such as species-specific habitat requirements, may contribute to genetic population structure. These results have implications for the conservation of riverine organisms, the reintroduction of locally extinct taxa and the rehabilitation of disturbed environments (Jump et al. 2009, Razgour et al. 2019).

In South Africa, it is legislated that catchments are used as management units (Republic of South Africa 1998). Previous studies on the genetic population structure of winged aquatic insects in South Africa have further supported the use of catchments as units for conservation (Wishart and Hughes 2003, Wishart et al. 2003, Price et al. 2010). The results found here for D. capensis and D. natalensis further highlight the genetic distinctiveness of populations between catchments, further corroborating the value of using catchments in conservation, management and legislative frameworks. These genetically distinct populations form an important component in the evolutionary legacy of a species. Therefore, the development of inter-basin water transfer schemes poses a threat to both D. capensis and D. natalensis and many other species by potentially connecting historically isolated and genetically distinct populations (Snaddon and Davies 1998, Davies et al. 2000)

In addition, dispersal among adjacent catchments has implications for the recovery of lotic systems following disturbance (Wishart and Davies 2003, Bellingan et al. 2019, Razgour et al. 2019). These factors should be considered in the development of strategies for the conservation of aquatic biodiversity (Wishart 2000, Thieme et al. 2007, Castello et al. 2013, Bellingan et al. 2019), and most particularly for high altitude catchments.

This study highlights the importance for future studies on community structure, biodiversity, and biomonitoring, where the taxonomic accuracy of species identification is crucial (Hajibabaei et al. 2016). The identification of possible cryptic species in A. sudafricanum and new species of Demoreptus affect the field of aquatic research in South Africa. Mayflies form a very important component of applied aquatic biology, particularly biomonitoring, the presence of cryptic taxa is being discovered at an increasing rate and poses challenges for some aquatic ecosystem monitoring methods. With bioassessment methods gaining increasing popularity, a detailed understanding of commonly collected species will aid in further development of assessment methods and clarify species identification (Delić et al. 2017, Suh et al. 2019). In addition, a deeper understanding of evolutionary processes and gene flow with regard to commonly occurring mayfly taxa contributes to broader research on ecosystem functioning and environmental processes. The utility of DNA barcoding for elucidating such phenomena is already proven (Jackson and Resh 1998, Plaisance et al. 2009, Raupach and Radulovici 2015) and widely used, with new technologies allowing for the rapid assessment of biodiversity using DNA metabarcoding (Pavan-Kumar et al. 2015, Elbrecht et al. 2017, Daravath et al. 2018, Alvarez-Yela et al. 2019). This approach to rapid biodiversity assessment has the potential to revolutionise and streamline management and conservation practices by providing detailed data for informed decision- and policy-making.