Spatial Variation of Acanthophlebia cruentata (Ephemeroptera), a Mayfly Endemic to Te Ika-a-Māui—North Island of Aotearoa, New Zealand

Simple Summary Patterns of variation within a widespread species can provide evidence of population history. Adjacent stable populations with gene flow show clinal genetic divergence forming a pattern of isolation by distance. Populations that have grown due to an increase in potential habitat such as forest cover since the last glacial maximum will have low genetic variation showing patterns of range expansion. The mayfly Acanthophlebia cruentata of Aotearoa, New Zealand, is widespread in streams of North Island but absent from the cooler South Island. Mayfly nymphs are restricted to streams but adults fly, facilitating gene flow among catchments. We detected higher genetic diversity at lower latitudes of this mayfly’s range compared to most of its distribution, concordant with predictions of limited forest cover in New Zealand during Pleistocene glacial periods. A signature of recent range expansion was observed in the higher latitudes. Despite initial observation suggesting mayfly size correlated with latitude, we found sex, elevation and sampling date were significant predictors of size, and some size variation is also explained by three regional groups based on haplotype distribution. Abstract The mayfly Acanthophlebia cruentata of Aotearoa, New Zealand, is widespread in Te Ika-a-Māui North Island streams, but has never been collected from South Island despite land connection during the last glacial maximum. Population structure of this mayfly might reflect re-colonisation after volcanic eruptions in North Island c1800 years ago, climate cycling or conceal older, cryptic diversity. We collected population samples from 33 locations to estimate levels of population genetic diversity and to document phenotypic variation. Relatively low intraspecific haplotype divergence was recorded among mitochondrial cytb sequences from 492 individuals, but these resolved three geographic-haplotype regions (north, west, east). We detected a signature of isolation by distance at low latitudes (north) but evidence of recent population growth in the west and east. We did not detect an effect of volcanic eruptions but infer range expansion into higher latitudes from a common ancestor during the last glacial period. As judged from wing length, both sexes of adult mayflies were larger at higher elevation and we found that haplotype region was also a significant predictor of Acanthophlebia cruentata size. This suggests that our mitochondrial marker is concordant with nuclear genetic differences that might be explained by founder effect during range expansion.


Introduction
Aotearoa, New Zealand, hosts a rich freshwater invertebrate fauna with high levels of endemism [1]. The distribution of freshwater species provides important regional biodiversity reflective of the country's varied environmental conditions and dynamic geological history [2,3]. Stream invertebrate composition is considered a key indicator of Insects 2022, 13, 567 3 of 19 AD232 Hatepe eruption in the Taupo Volcanic Zone Taupo [28] distributed some 65 km 3 of ash and rock over about 3000 km 2 [29], and had a major impact on nearby biology [30]. It is expected that abrupt changes in sediment, temperature and chemistry associated with volcanic fallout would have been deleterious for freshwater life [31]. Population recovery after such events is typically characterised by lowered genetic diversity and evidence of recolonisation from elsewhere [21,32,33]. Alternatively, it is plausible that A. cruentata comprises cryptic diversity that would be indicated by distinct regional lineages. expect a phylogeographic signature of southward population expansion. Another plausible source of population disruption is volcanism [25]. The high energy AD232 Hatepe eruption in the Taupo Volcanic Zone Taupo [28] distributed some 65 km 3 of ash and rock over about 3000 km 2 [29], and had a major impact on nearby biology [30]. It is expected that abrupt changes in sediment, temperature and chemistry associated with volcanic fallout would have been deleterious for freshwater life [31]. Population recovery after such events is typically characterised by lowered genetic diversity and evidence of recolonisation from elsewhere [21,32,33]. Alternatively, it is plausible that A. cruentata comprises cryptic diversity that would be indicated by distinct regional lineages.

Materials and Methods
Specimens of Acanthophlebia cruentata were collected at small gravel streams in catchments with native vegetation, either as nymphs by hand searching stones or kick-sampling with a net, or, on emergence by UV trapping [11,36]. Specimens were preserved in ethanol for subsequent analysis. For analysis of genetic variation, we targeted locations  [19], (B) [34], (C) [35]. Here, we examine the phylogeographic structure and adult phenotypic variation within Acanthophlebia cruentata across its full spatial range, to examine these possibilities.

Materials and Methods
Specimens of Acanthophlebia cruentata were collected at small gravel streams in catchments with native vegetation, either as nymphs by hand searching stones or kick-sampling with a net, or, on emergence by UV trapping [11,36]. Specimens were preserved in ethanol for subsequent analysis. For analysis of genetic variation, we targeted locations for which published sequence data for this species were not available, particularly in the southern and northern extents of the species' distribution. For examination of phenotypic variation of adults, we selected locations to ensure spatial coverage, but also representative altitudes across the species' distribution ( Figure 2, Table 1).

Genetic Variation
Ethanol preserved specimens (nymphs, subimagos or imagos) were used for whole genomic DNA extraction using a salting-out method [37]. Tissue was macerated and incubated with 10 µL of 10 mg/mL proteinase-K in 600 µL of TNES buffer (20 mm ethylenediaminetetraacetic acid, 50 mM Tris, 400 mM NaCl, 0.5% sodium dodecyl sulphate) at 55 • C for 1-4 h. 10% 5M NaCl was added and the extractions shaken vigorously for one minute before centrifuging at 16,000 g for 5 min. The supernatant was removed and precipitated with an equal volume of cold 100% ethanol. DNA was collected by spinning and washed with 70% ethanol, then air-dried and dissolved in water.
Mitochondrial DNA sequences were obtained using primers that target a portion of the cytochrome-b gene (cytb) homologous with published data for A. cruentata [25]. Polymerase chain reactions (PCR) used the primers CB-J-10612 and CB-N-10920 [38] in 10 or 20 µL reactions (200 µm dNTPs, 2 mm MgCl 2 , 0.5 U Taq), treated to 36 cycles of 94 • C for 30 s, 54 • C for 30 s, 72 • C for 90 s with an initial denaturation of 94 • C for 2 min and a final extension at 72 • C for 5 min. We also generated sequences for a portion of the mtDNA cytochrome oxidase subunit 1 (cox1) for a subsample of Acanthophlebia cruentata to provide comparative 'DNA-barcode' information, using PCR primers LCO1490 and HCO2198 [39].
DNA sequences were checked for ambiguities, translated and aligned in Geneious [40]. We identified distinct haplotypes including those (#1-#34) recorded by Smith et al. (2006) [25] and continued the same numbering system (#35-#52, #54-#84), for a total of 83 haplotypes. We examined the spatial distribution of haplotypes across the range of the species. Population sample locations were mapped in R [41] using the maps package [42] and evolutionary relationships among haplotypes were inferred with a minimum spanning network [43] implemented in PopArt [44]. During relatively slow range expansion, the populations at the leading edge might be characterized by increasingly derived sequences due to mutations accumulating on external branches [45]. However, evidence of recent recolonisation would be seen as the sharing of the same haplotypes among source and putative sink populations.
We used spatial Principal Component Analysis [46], implemented with the R package adegenet [47] to visualise the non-random geographic distribution of genetic variation using simple distances between DNA sequences. Spatial PCA (sPCA) utilises spatial coordinates via a connection network to explore non-random distribution of genetic variation among population samples, and uses a Monte-Carlo simulation to assess significance of structure. A useful feature of this approach is that it does not assume population groupings a priori. Based on this analysis, population samples were then grouped for analysis to identify signal consistent with hypothesised prehistoric partitioning. We used ARLEQUIN v3.5.2.2 [48] for AMOVA with population sample data grouped into three regional sets (north, west, east). southern and northern extents of the species' distribution. For examination of phenotypic variation of adults, we selected locations to ensure spatial coverage, but also representative altitudes across the species' distribution ( Figure 2, Table 1).

Figure 2.
Locations in North Island, New Zealand of Acanthophlebia cruentata population samples. Black circles-previously reported locations [25], black squares-this study. Abbreviated codes refer to locations listed in Table 1. Degrees of latitude are given on the left; north is towards the equator. Inset: Male subimago Acanthophlebia cruentata from Pukenui Forest © Olly Ball and Steve Pohe Collection. Black circles-previously reported locations [25], black squares-this study. Abbreviated codes refer to locations listed in Table 1. Degrees of latitude are given on the left; north is towards the equator. Inset: Male subimago Acanthophlebia cruentata from Pukenui Forest © Olly Ball and Steve Pohe Collection. Table 1. Sampling of Acanthophlebia cruentata in North Island, New Zealand ordered by latitude. Shaded entries relate to population samples reported here, and others are as given by Smith et al., (2006) [25]. N is the number of individuals sequenced from that location. Latitudes and longitudes are in decimal degrees from the New Zealand Geodetic Datum 2000 system. Region is the inferred grouping of population samples. Range expansion is expected to leave a signature of increasingly divergent allelic frequencies and reduced genetic variation at the leading edge [44,49]. Southward population expansion from LGM forest refugia would result in a southward cline of reduced genetic diversity, while post-Tāupo eruption re-colonisation would result in lower diversity in central North Island with higher diversity either side. We calculated haplotype and nucleotide diversity using DnaSP v6.0 [50] and ARLEQUIN v3.5.2.2 [48] and examined population sample nucleotide diversity (π) with respect to sample location (latitude).

Site
Range expansion is usually associated with population growth and therefore we predicted a signature of southward demographic expansion if this mayfly species was restricted to higher latitudes in Northland during the LGM. We calculated neutrality statistics and population size changes for the entire intraspecific sample and for regional subsets (north, west, east). We looked for evidence of demographic expansion and population size constancy considering the frequency of pairwise genetic distances among haplotypes [51,52].
Genetic structure resulting from gene flow rather than range expansion would leave a signature of isolation by distance (IBD). To estimate population differentiation we calculated population sample pairwise Φ ST [53] with Arlequin, and tested to see whether estimates were significantly greater than 0 using 1000 sample permutations. We compared pairwise Φ ST with linear distances between sampling sites, obtained by triangulation of New Zealand Map Grid coordinates, and used a mantel test implemented in R with the package vegan [54] to determine whether there was a signature of IBD.

Size Variation
We compared 1064 adult Acanthophlebia cruentata from 33 sites across the species range to document intraspecific phenotypic variation, using forewing length as a proxy for individual size [55]. A single forewing from each specimen was mounted under a coverslip in a syracuse dish with 80% ethanol and measured to the nearest 0.1 mm from the apex of the wing to the point of attachment at the basal plate with an eyepiece graticule through a Leica stereomicroscope (Model: MZ95) at 8-10× magnification [11]. Sampling of winged mayflies took place in one flight season between November 2013 and February 2014. All measurements were made by the same person ( Figure 3). nsects 2022, 13, x FOR PEER REVIEW 7 of 2 central North Island with higher diversity either side. We calculated haplotype and nu cleotide diversity using DnaSP v6.0 [50] and ARLEQUIN v3.5.2.2 [48] and examined pop ulation sample nucleotide diversity (π) with respect to sample location (latitude). Range expansion is usually associated with population growth and therefore we pre dicted a signature of southward demographic expansion if this mayfly species was re stricted to higher latitudes in Northland during the LGM. We calculated neutrality statis tics and population size changes for the entire intraspecific sample and for regional sub sets (north, west, east). We looked for evidence of demographic expansion and population size constancy considering the frequency of pairwise genetic distances among haplotype [51,52].
Genetic structure resulting from gene flow rather than range expansion would leav a signature of isolation by distance (IBD). To estimate population differentiation we cal culated population sample pairwise ΦST [53] with Arlequin, and tested to see whether es timates were significantly greater than 0 using 1000 sample permutations. We compared pairwise ΦST with linear distances between sampling sites, obtained by triangulation o New Zealand Map Grid coordinates, and used a mantel test implemented in R with th package vegan [54] to determine whether there was a signature of IBD.

Size Variation
We compared 1064 adult Acanthophlebia cruentata from 33 sites across the specie range to document intraspecific phenotypic variation, using forewing length as a proxy for individual size [55]. A single forewing from each specimen was mounted under a co verslip in a syracuse dish with 80% ethanol and measured to the nearest 0.1 mm from th apex of the wing to the point of attachment at the basal plate with an eyepiece graticul through a Leica stereomicroscope (Model: MZ95) at 8-10× magnification [11]. Sampling of winged mayflies took place in one flight season between November 2013 and February 2014. All measurements were made by the same person ( Figure 3). Analysis used R with the package ggplot2 [56] to plot density distributions of mal and female wing lengths by sample location. Wing length was analysed with linear mixed effects models using R packages lme4 [57] and lmerTest [58]. Wing length was log-trans formed to remove positive skew and heteroscedasticity in residuals. Adult body size (in cluding wing length) in mayflies is expected to be temperature dependent [59] so the co variates, elevation (100 m units), latitude (100 km units) and time of season (week afte start of sampling on 22 November) were included. Preliminary plots suggested non-linea relationships of wing length with these covariates so quadratic effects were added. Sex Analysis used R with the package ggplot2 [56] to plot density distributions of male and female wing lengths by sample location. Wing length was analysed with linear mixed effects models using R packages lme4 [57] and lmerTest [58]. Wing length was log-transformed to remove positive skew and heteroscedasticity in residuals. Adult body size (including wing length) in mayflies is expected to be temperature dependent [59] so the covariates, elevation (100 m units), latitude (100 km units) and time of season (week after start of sampling on 22 November) were included. Preliminary plots suggested non-linear relationships of wing length with these covariates so quadratic effects were added. Sex, stage (imago, subimago) and region (based on haplotype groups, see Figure 4) were fixed factors and sampling location a random effect (R code in Supplementary Material).   Supplementary Table S1). This genetic variation is also depicted by interpolations of the first two global axes of spatial PCA from these data revealing overall trends in spatial structure (right). Inset is late instar A. cruentata nymph.

Diversity
We obtained mtDNA cytb sequences from 308 Acanthophlebia cruentata individuals, from 33 population samples spanning the known range from North Cape to Wellington, North Island Aotearoa-New Zealand (Figure 1). These complemented published data representing 185 individuals from 19 other locations from part of this area [25]. The maximum population sample size for the genetic analysis was 13 sequenced individuals with a mean sample size of 9.7 and mode of 10 (Supplementary Table S1). From two sites (T4  Supplementary Table S1). This genetic variation is also depicted by interpolations of the first two global axes of spatial PCA from these data revealing overall trends in spatial structure (right). Inset is late instar A. cruentata nymph.

Diversity
We obtained mtDNA cytb sequences from 308 Acanthophlebia cruentata individuals, from 33 population samples spanning the known range from North Cape to Wellington, North Island Aotearoa-New Zealand (Figure 1). These complemented published data representing 185 individuals from 19 other locations from part of this area [25]. The maximum population sample size for the genetic analysis was 13 sequenced individuals with a mean sample size of 9.7 and mode of 10 (Supplementary Table S1). From two sites (T4 Waingongoro, KH Khandallah) we had only a single specimen (each with a haplotype common in an adjacent location sample), these data were excluded from most analyses. This gave a data set of 492 individuals in the cytb alignment of 280 bp with 57 variable positions, comprising 83 haplotypes (49 of them new to the present study) with an overall nucleotide diversity (π) of 0.01034 (Table 2). These data were characterised by high haplotype diversity (Hd: 0.861), but low genetic divergence (mean K2P: 0.93%). New GenBank accessions ON227131-ON227179 (Supplementary Table S1). More than 50% of all specimens from which we obtained cytb data (254/492) had one of two common haplotypes (#16 and #27), while more than three quarters of all haplotypes (85.7%) were recorded in just one or two individuals (Supplementary Table S1). One or other of the common haplotypes (#16 or #27) were observed in 80% of population samples; sixteen eastern samples with #27, and twenty-four western samples with #16 ( Figure 4). Only one population sample (BP; Bay of Plenty) had both of the common haplotypes. Population samples lacking both of the common haplotypes came from Northland and Little Barrier Island (Figure 4). In the north only one sample (NC; North Cape, n = 8) had individuals with a haplotype very similar to the common haplotype (#80, #81 both one-step from #16; Figure 4). North of Whangaparāoa ( Figure 4) the cytb haplotypes differed from one another by up to eight substitutions (2.86%), but few of these haplotypes are recorded further south. One rare exception was haplotype #45 recorded in individuals in the north (NC, n = 1), Coromandel Peninsular (FB, n = 1) and Waikato (TA, n = 1).
Interpolation of the spatial distribution of the cytb genetic variation using two global sPCA axes revealed three regions with distinct genetic diversity: north, west and east ( Figure 4). An imaginary line running north northeast-south southwest through central North Island (Lake Taupo) would fall between our population samples with one or other of the two common haplotypes and their derivatives (haplotype #16 in the west, haplotype #27 in the east). AMOVA indicated that more than 50% of genetic variation could be explained by the regional groups north, west and east, the remainder was explained by variation within groups (28%) and within population samples (22%), reflecting the high haplotype diversity.

Range Expansion
Genetic diversity as measured by the number of mtDNA haplotypes per individual within population samples (h) was generally higher in the north and northwest. For example, samples of 12 individuals from Whangarei (WG) and Waikato (TA) each had seven or eight haplotypes, while some similar sized samples from the east had a single haplotype (e.g., MO, KI; Table 2).
Population sample genetic diversity (h) was lowest at the southern limit of this mayfly's range (Table 3). Population sample nucleotide diversity (π) showed a similar pattern across North Island, New Zealand (Table 2; Figure 5A). The average π for the eleven northern population samples was 0.0105, for West 0.0073 (23 population samples) and the lowest estimate of π was in the East (0.00123; 15 population samples; Table 3) reflecting this uneven distribution of genetic variation.
Demographic analyses applied to the total cytb data set returned significant negative values for all statistics (Tajima D −1.94599 sig p < 0.05, Fu & Li's D* −7.18167 sig p < 0.02, Fu & Lis F* −5.624 sig p < 0.02, Fu's Fs −98.926), which, under the assumption of functional neutrality of the marker, is a signature of rapid recent population increase. However, when subdivided into regional population sample groups (North, West, East as in Figure 4), statistical support for expansion in the North was less compelling; Tajima's D (−1.54881) was not significantly negative (p > 0.10) whereas the West (Tajima's D −1.76828, p < 0.05) and East (Tajima's D −2.35176, p < 0.01) were (Table 3). Viewing the data in a different way, a ragged mismatch distribution for the total data was dominated by the northern sample and indicative of a stable population whereas unimodal distributions as in the east are consistent with range expansion.  sequence diversity for all and regional subsets (see Figure 4). Table 3. Summary statistics of regional variation among partial mtDNA cytochrome-b sequences from Acanthophlebia cruentata mayflies in New Zealand (π, nucleotide diversity; k, average number of nucleotide differences; S, number of segregating sites; h, number of haplotypes; Hd, haplotype diversity), and 'neutrality' statistics that can indicate population size change. Significantly negative values of Tajima's D (bold) imply population expansion under the assumption of locus neutrality. R2 represents the relationship between singleton mutations and average nucleotide difference with a smaller value (bold) consistent with a scenario of population expansion. A larger value (bold) of the Raggedness statistic r indicates a poor fit with the assumption of constant population size. Regional groupings as shown in Figure 4 and Table 1. Analysis seeking evidence of isolation by distance (IBD) using pairwise population sample differentiation (Φ ST ) also indicated regional-specific processes. A significant correlation between geographic and genetic distances was found when all population samples were analysed (p = < 10 −6 ) but with a mediocre coefficient (r = 0.4146). However, we found strong statistical support (p = 0.00004) and a higher correlation coefficient (r = 0.7164) within the northern population samples ( Figure 6). In contrast, analyses of geographic and genetic distances within the west and east population samples (see Figure 4) gave marginal (West, p = 0.01214, r = 0.2467) or no signal of IBD (East, p = 0.3196, r = 0.05115). In combination with the demographic analyses we infer longer standing population structure in the north, but recent, rapid population expansion elsewhere that is consistent with mismatch distributions.  Figure  4).   DNA sequence from a fragment at the 5 -end of cox1. An alignment of 585 bp contained 27 variable sites and we identified 34 distinct haplotypes (GenBank ON442277-ON442310). Among those haplotypes the average Jukes-Cantor distance was 0.00728 (0.73%) and the maximum was 0.014 (1.4%).

Timeframe
Mean DNA sequence divergence within Acanthophlebia cruentata mitochondrial haplotypes was low (K2P: cytb 0.93%, cox1 0.73%). Applying substitution rates inferred for insect intraspecific mtDNA molecular evolution of 0.0285-0.0792 substitutions per site per million years [60,61] to the median of these mtDNA distances (0.074) we can infer that the common ancestor of the diversity observed within A. cruentata existed approximately 38,000-50,000 years ago. Though by no means a robust time calibration, this coincides with the last glacial maximum [62] and is consistent with a southward range expansion of this species when forest habitat would have been extending and there was land connection between North and South Islands [19].

Wing Length-Adult Size
Forewing lengths of 1064 individuals (imagos and subimagos) from 33 sites were measured and analysed taking into account latitude, elevation, instar and time of sampling, and lineages defined from cytb sequence data ( Figure 4). Density distributions for 21 population samples with >10 individuals were graphed showing smaller size of males (Range: 7.4-12.5 mm, mean 9.96 mm, n = 533) compared to females (Range: 8.1-13.8 mm, mean 10.87 mm, n = 417) but different mean sizes across the species range (Figure 7). Data from all 33 sites were used for statistical model fitting. Attempts to fit a linear mixed effects model with interaction terms among elevation, latitude and season of sampling were abandoned as these variables are strongly confounded due to geographical constraints and the temporal pattern of sampling. In addition, the combinations of these interactions could equally be covered by the location random effect. Females were generally larger than males (0.9 mm or approximately 9%) but the response to other predictor variables was also stronger in females resulting in many significant interactions. To simplify interpretation, models that were fitted separately for males and females are reported here. abandoned as these variables are strongly confounded due to geographical constraints and the temporal pattern of sampling. In addition, the combinations of these interactions could equally be covered by the location random effect. Females were generally larger than males (0.9 mm or approximately 9%) but the response to other predictor variables was also stronger in females resulting in many significant interactions. To simplify interpretation, models that were fitted separately for males and females are reported here. After stepwise removal of non-significant effects the final model included the random location effect (55% of variance), quadratic elevation, and haplotype region (Table  4). Male wing length is greatest at 500 m elevation where it is about 19% longer than at sea level. Locations in the northern haplotype region (Figure 4) are generally lower elevation than elsewhere in North Island but in addition to the elevation effect, northern males were 7% smaller than those in the eastern region. After stepwise removal of non-significant effects the final model included the random location effect (55% of variance), quadratic elevation, and haplotype region (Table 4). Male wing length is greatest at 500 m elevation where it is about 19% longer than at sea level. Locations in the northern haplotype region (Figure 4) are generally lower elevation than elsewhere in North Island but in addition to the elevation effect, northern males were 7% smaller than those in the eastern region.
The results of model fitting and simplification (Table 4) were similar in females but with the addition of a weak quadratic seasonal effect. Wing length was greatest early in the season and declined with increasing slope to 5% smaller by the end of sampling 13 weeks later. The elevation and regional effects on wing length were stronger in females compared to males. Female wing length plateaued at 460 m asl, 31% longer than at sea level ( Figure 8). Females in the northern region were 15% smaller than those in the eastern region. The random location effect accounted for 58% of the variance, similar to that for males. No effect of latitude was detected in either sex. The results of model fitting and simplification (Table 4) were similar in females but with the addition of a weak quadratic seasonal effect. Wing length was greatest early in the season and declined with increasing slope to 5% smaller by the end of sampling 13 weeks later. The elevation and regional effects on wing length were stronger in females compared to males. Female wing length plateaued at 460 m asl, 31% longer than at sea level ( Figure 8). Females in the northern region were 15% smaller than those in the eastern region. The random location effect accounted for 58% of the variance, similar to that for males. No effect of latitude was detected in either sex.

Discussion
The monotypic genus Acanthophlebia might have contained cryptic taxa reflective of a long history of diversification in New Zealand, but we did not find that. Instead we found low genetic divergence within Acanthophlebia cruentata sampled throughout its range in North Island, New Zealand; from North Cape in the Far North to Wellington in

Discussion
The monotypic genus Acanthophlebia might have contained cryptic taxa reflective of a long history of diversification in New Zealand, but we did not find that. Instead we found low genetic divergence within Acanthophlebia cruentata sampled throughout its range in North Island, New Zealand; from North Cape in the Far North to Wellington in the south, Taranaki in the west and East Cape. Within the large sample of A. cruentata mean sequence divergence was less than 1% (cytb 0.93%, cox1 0.73%), which is at the lower end of the intraspecific range (1.2-1.4%) reported for New Zealand mayflies, stoneflies and caddisflies [63]. Other New Zealand mayflies have not been sampled at this scale (see [64]), but unpublished DNA sequences for a selection of New Zealand species sampled across their ranges revealed intraspecific variation among cox1 sequences ranging from 0.6% (Arachnocolous phillipsi) to 3.9% (Zephlebia versicolor). This contrasts with much higher sequence divergence within a number of endemic terrestrial insects in the same landscape; average mtDNA distances of~8% in the wētā Hemideina crassidens [65] and Hemiandrus pallitarsus [66] suggest large sustained populations of these forest species.
Mayflies are habitat sensitive aquatic stream invertebrates whose populations are likely to be responsive to geophysical processes that abruptly alter local conditions [67]) and so it is not surprising that we found population structure consistent with very recent environmental changes in North Island (Figure 1) [23]. We detected isolation by distance in the northern part of the Acanthophlebia cruentata but a signature of population expansion in the south. Despite low divergence among the cytb sequences the haplotype diversity revealed a pattern of regionally-restricted variants delimitating three sections of the North Island. Population samples north of Whangaparāoa contained region-exclusive haplotypes and relatively high nucleotide diversity per location. In contrast, to the south, genetic diversity was lower and dominated by two common haplotypes. Here, population samples fell into east and west groups, each comprising a common haplotype and its derivatives. This pattern of haplotype distribution is typical of recent range expansion scenarios.
Overall we found a pattern of population expansion [44,68] in central and southern North Island, contrasting with a more stable population history in Northland. The prominent east-west structure implies two parallel range expansions that likely reflect occupation of habitat in recent geological time. The directionality of range expansion within the east and west regions cannot be inferred, although the northern population samples within the western group have genetic diversity levels equivalent to the northern population samples. Within both east and west there is statistical evidence of population expansion but not isolation by distance ( Figure 6). It is probable that the phylogeographic structure of Acanthophlebia cruentata reflects changing habitat opportunities during Pleistocene climate cycling, and we found that the level of genetic diversity indicated a common ancestor about 38,000-50,000 years ago. This suggests climate during the last (Otiran) ice age (~115,000-11,700 years ago) reduced population size of this mayfly, and the impact may date to the coldest phases associated with ice advances since~65,000 years ago [69]. All putative genetic diversity accumulated in previous interglacials or prior to the Pleistocene has been erased, leaving this species on a relatively isolated phylogenetic lineage with shallow clade depth. As a result we find no signal for cryptic diversity or population structuring associated with land formation in southern North Island [19]. Similarly, we do not find strong evidence that volcanic eruptions at Taupo influenced local mayfly population structure (see sPCA Figure 4) despite deep ignimbrite deposits [70]. This implies a rapid reoccupation of streams around the crater from surviving populations to the east and west, although we note few unique alleles in population samples from sites near Taupo (KK-Kakaho #47; OM-Omarowa #51: Supplementary Table S1).
In New Zealand, a pattern of high intraspecific mtDNA diversity in Northland and less to the south is found in the tree wētā Hemideina thoracica with a concordant decline in nucleotide diversity (π) southwards [21,71]. A number of other New Zealand forest species have a genetic signature of northern refugia, suggesting a similar response to climate shifts, including the stick insect Clitarchus hookeri [27], and giraffe weevil Lasicrhynchus barbicornis [72]. Genetic signatures of population expansion were not resolved in the cicada Kikihia cutora but three mtDNA clades were identified with distributions that closely match the north, west, east locations of the genetic groups within A. cruentata [73]. Interestingly, levels of intraspecific haplotype diversity are higher in all these examples, suggesting that population reductions during periods of inhospitable conditions were more intense in Acanthophlebia cruentata than these other endemic insects. Given the possibility that the time of population expansion for these different species was not concordant we cannot conclude that the processes involved were the same (Cf [73]).
Intraspecific size variation in stream insects generally correlates with sex and water temperature associated with latitude and elevation, with smaller winged, small bodied adults (of each sex) at warmer (northern) latitudes and lower elevation [59,74]. This is largely the case for Acanthophlebia cruentata, with smaller wings observed from warmer streams, but after accounting for sex and elevation, our models found no support for an effect of latitude. Adult size, measured here as forewing length, is greater in females and they also show stronger responses to other predictors of wing length than males. Wing length was greater at higher elevation as expected from the thermal equilibria hypothesis [59], but unexpectedly declined at the highest elevation sites. Maximum wing length was observed at a similar elevation for males and females (500 m and 460 m) but there are only five sample locations above this elevation (max 775 m) so it could be simply a 'plateau' effect. In females there was a significant decline in wing length through the flight season which is also expected from the thermal equilibrium hypothesis as later emerging mayflies would be exposed to higher temperatures during nymphal growth. The three geographical groups (regions defined by cytb haplotypes) are strongly confounded with these proxies for temperature. The populations in the east were sampled later in the flight season and also included sites at higher elevation, especially compared to the northern region; nevertheless haplotype region was a significant predictor of Acanthophlebia wing length even in competition with strong proxies for temperature and after accounting for sampling location effects. This significant effect of haplotype region in our data, hinting at a genetic difference between Acanthophlebia cruentata in the north, west and east, might be explained by founder effect on functional genes [75].

Conclusions
The mayfly Acanthophlebia cruentata has, in the context of endemic insects of Aotearoa-New Zealand, shallow lineage diversity. Nevertheless, we are able to identify from the distribution of genetic variation a spatial pattern of late Pleistocene age. We find signal for two, parallel, range expansions within the mtDNA data and an indication that adult mayfly size is linked to distinct, though shallow, regional haplotype clusters. Changing habitat availability associated with Pleistocene climate cycling provides the parsimonious explanation for this shallow population history but it is plausible that the Acanthophlebia lineage itself has a shallow ancestry in Aotearoa-New Zealand.