High diversity and local endemism in Aotearoa New Zealand's groundwater crustacean fauna

Abstract We used DNA barcoding to assess the diversity and distribution of New Zealand's groundwater amphipods and isopods (Crustacea) and to determine whether biodiversity and endemism within tectonically active New Zealand are similar to those of more tectonically stable continents. Sixty‐five wells were sampled in seven aquifers across four regions within the North and South islands of New Zealand, and resident invertebrates were morphologically identified and then assessed using sequencing of the mitochondrial DNA cytochrome c oxidase subunit one (COI) gene. Invertebrates were found in 54 wells. Of the 228 individual amphipods and isopods found in 36 of the wells, 154 individuals were successfully sequenced for COI (68% success rate) from 25 wells, with at least one well in each aquifer containing sequenced individuals. Of the 45 putative species identified using Barcode Index Numbers (BINs), 30 BINs (78% of all taxa and 83% of amphipods) were previously unrecorded. Substantial morphologically cryptic, species‐level diversity was revealed, particularly within the amphipod Family Paraleptamphopidae. Similarly, one isopod taxon morphologically identified as Cruregens fontanus was assigned to five well‐separated BINs based on COI sequences. Endemism appeared high, with all taxa regionally endemic; 87% of species were restricted to one aquifer and more than 50% restricted to one well. Non‐saturated species accumulation curves indicated that, while additional sampling may increase the range of some currently identified taxa, additional range‐restricted taxa are also likely to be discovered. Patterns of diversity and short‐range endemism were similar to those found elsewhere, including locations which are more tectonically stable. The predominance of local endemism within New Zealand's groundwater fauna suggests that land‐use activities and groundwater extraction require careful evaluation to minimize threats to groundwater biodiversity.


| INTRODUC TI ON
Groundwater systems are often viewed as lifeless conduits of subsurface water flow (sensu Hancock & Boulton, 2008). However, research over the last few decades has identified a rich diversity of groundwater fauna (the stygofauna), which provide important ecosystem services . Stygofaunal communities are typically dominated by invertebrates and are characterized by high levels of biodiversity, particularly Crustacea (Danielopol et al., 2000;, and by endemism over small spatial scales (Boulton, 2020;Hancock & Boulton, 2008).
Logistical difficulties in sampling groundwater ecosystems (Larned, 2012) and the often cryptic morphology of stygofauna (Bradford et al., 2010;Danielopol & Pospisil, 2001;Finston et al., 2007) have meant that biodiversity inventories of subterranean ecosystems are severely lacking in many locations (Ficetola et al., 2019;. Investigating the spatial scales of endemism within groundwater ecosystems is a critical step in understanding the implications of increasing threats, such as water abstraction and contaminant infiltration, as well as the efficacy of different management policies and practices (Boulton, 2020;Mammola et al., 2019).
Endemism over relatively small spatial scales appears to be high in most groundwater systems (Danielopol et al., 2003). For example, DNA sequencing of 14 nominal, widespread species indicated more than 50 morphologically cryptic amphipod lineages (Trontelj et al., 2009), most with highly restricted spatial distributions. Fortyone percent of the stygobitic (obligate groundwater dwellers) species found across six European regions were reported from areas <500 km 2 ) and ranges of <200 km were common (Trontelj et al., 2009). Some taxa were even restricted to a single cave or sampling location (Gibert & Deharveng, 2002). (Amphipoda), which lacks any pigmentation or eyespots, reflecting its subterranean existence.
For stygobionts, low dispersal capabilities, coupled with geographical isolation over evolutionary time scales, will result in genetic divergences among populations, leading to small-scale or short-range endemism (Harvey, 2002;Harvey et al., 2011). Shortrange endemism has been identified in areas that are tectonically stable including Australia (e.g., Hancock & Boulton, 2008) and North America  and where physical barriers, such as glacial deposits or catchment shrinkage (due to aridity), subdivide aquifers and isolate populations. New Zealand's active tectonic environment, with relatively rapid uplift, subsidence, erosion, and deposition (Brown, 2001), provides a contrasting setting for stygofaunal evolution compared to that of other continents and could suggest an alternative to the emerging paradigm of short-range endemism within stygofaunas.
Much of New Zealand's groundwater resides within alluvial aquifers underlying extensive plains comprising relatively young, unconsolidated, and often highly porous matrices, resulting in high hydraulic conductivities and high interstitial water velocities (Close et al., 2002;Pang et al., 1998). This might be expected to facilitate stygofaunal movement within or between aquifers. However, the country's mountainous terrain is also likely to provide physical barriers that would facilitate short-range endemism. For example, New Zealand's spring and spring-stream hydrobiid snails include several examples of allopatric, short-range endemic species (Haase, 2008).
Previous work on New Zealand's stygofauna has largely relied on morphological identifications of taxa, with three families of amphipods and two families of isopods conspicuously present (e.g., Fenwick, 2001;. This suggests either a smaller number of widespread taxa or a larger number of more restricted, morphologically cryptic taxa. Molecular markers, such as the mitochondrial DNA cytochrome c oxidase subunit one gene (COI), are particularly helpful for identifying morphologically conservative taxa, including Crustacea (Costa et al., 2007;Hogg et al., 2006;Watson et al., 2015). For example, molecular studies have revealed that several European stygofaunal species, previously considered widespread within karstic environments of southern and western Europe, actually comprise several morphologically cryptic taxa, each confined to single locations or catchments, with geographic ranges comprising single or multiple localities spanning no more than c. 180 km (e.g., Ferreira et al., 2007;Lefébure et al., 2006Lefébure et al., , 2007. Similarly, morphologically cryptic, subterranean stygofauna (e.g., amphipods, isopods, and water beetles) inhabiting groundwater calcretes in the arid Yilgarn region of Australia are actually endemic to single calcrete aquifers, with some ranges smaller than a few square kilometers (Cooper et al., 2007).
Here, we assess stygofaunal diversity across New Zealand using COI gene sequences. We focus on amphipod and isopod crustaceans, as they generally dominate the stygofaunal assemblages of shallow alluvial aquifers (Gibert & Deharveng, 2002) and compare stygofaunal diversity and endemism within New Zealand to more tectonically stable continents.
F I G U R E 1 A typical groundwater (stygobitic) crustacean, Paracrangonyx sp. (Amphipoda), which shows the lack of pigmentation or eyespots. The head is to the lower right of the photo. Photo credit: N. Boustead 2 | MATERIAL S AND ME THODS

| Study sites
Sampling locations were stratified hierarchically, across: (1) North and South islands of New Zealand; (2) regions within the South Island; (3) aquifers within regions; and (4) elevation within an aquifer ( Figure 1). Sampling was focused on larger alluvial aquifers in four regions along the drier, east coasts of both islands. The largely north-south orientation of New Zealand and its associated mountain ranges, in conjunction with predominantly westerly weather patterns, results in orographic precipitation on the west coast and a drier east coast. Further, we focused on alluvial aquifers to reduce variability in invertebrate communities potentially caused by geographical or hydrogeological differences. These aquifers generally extend from foothills to the coast and have high potentials to store and transport groundwater Tschritter et al., 2017). We collected from the uppermost aquifer at each location.
Candidate aquifers were identified using a two-dimensional aquifer map (Ministry for the Environment, 2015), which was generated using data from White (2001) and updated by Moreau and Bekele (2015). The location of major aquifers corresponds with a more recent map using finer-scale GIS data . Given the unknown quantity of water exchange between aquifers at a small scale, we conservatively assigned wells to the major aquifers identified by both maps and named in Ministry for the Environment (2015). We caution that the "aquifers" identified in this report may comprise two or more smaller aquifers that are variously hydrologically connected.
We assigned a group of wells within the Moutere Valley aquifer to the adjacent Motueka River Terraces aquifer ( Figure 1) because our sampling sites were close (440-630 m) and likely hydrologically connected to the Motueka River, which traversed both aquifers.
Also, because the Waimea Plains are underlain by multiple major and minor aquifers with some hydrological inter-connectivity (White, 2001), we assigned our sampling wells in this area to a composite "Waimea Plains" aquifer ( Figure 1). There was no established name for the aquifer beneath the Southland wells, so we named this after the nearby Mataura River. Sampling locations within aquifers were restricted to existing wells where sampling equipment could be deployed. Candidate wells were identified with help from local groundwater monitoring agencies (regional and district councils) and had been installed for multiple purposes, including water quality monitoring and research, and for water abstraction.
Sixty-two wells were sampled once across the four regions and two islands (Figure 2). Due to logistical constraints of finding and accessing suitable wells, the number of wells sampled varied between aquifers. However, at least two wells (and a maximum of 14 wells) were sampled within each aquifer. Where possible, these wells were located across a range of elevations within an aquifer (Appendix 1). We also included invertebrates collected during other sampling excursions from three Canterbury wells (Central Plains aquifer) because they contained specimens that complemented those from the current sampling program. This resulted in a total of 65 wells sampled. We use the term "well" to include both traditional wells (installed by excavation) as well as drilled bore holes. Sampling was undertaken between May 12, 2017, and December 1, 2017, with one well sampled on March 1, 2018.

| Invertebrate sampling
We used two main sampling methods to maximize capture rates.
Firstly, we pumped 60 or 100 L (depending on well flow rate) of water from the screened (i.e., open to the aquifer via slots or perforations) section of the well through a 200µm mesh collecting bag. Neoprene flanges or inflatable packers were used to restrict pumping to the screened section of well. For shallower wells (water table < c. 8 m, 48 wells 77% of wells), we used a Bou-Rouche pump (Malard et al., 2002). A pneumatic Bennett pump (Bennett Sample Pumps Inc., Amarillo, Texas; pumping rate c. 30 L/min. on average) was used for 11 deeper wells. Secondly, a plankton net (64 or 100 µm mesh, depending on the amount of suspended sediment) with a flexible rim was folded into a weighted bailer, lowered to bottom of the well, bounced to suspend sediment and any associated stygofauna, and retrieved slowly to filter the entire water column.
The bailer was then retrieved, and its contents added to the contents of the net haul. This procedure was repeated three times. Contents of repeat net and bailer collections from each well were pooled.
Sampling methods were modified for sampling very large and very small wells (casing internal diameters <50 and >500 mm diameter, 5% of wells) where the flanges or packers could not be used. Three samples were also collected coincidentally with other activities (e.g., well conditioning or purging) and employed other field methods and filtered larger volumes of groundwater. In all cases, stygofauna were collected and concentrated using a 200µm mesh bag. All samples were preserved in the field with 100% ethanol, chilled, transported to the laboratory, and stored in the dark at −20°C until needed for further processing. All equipment was washed thoroughly after sampling each well and air-dried between regions to avoid transferring any specimens between wells.
In the laboratory, the contents of each sample were concentrated on a 250µm sieve, sorted under a stereomicroscope into separate vials for each recognizable taxon, and stored in the dark at −20°C in 100% ethanol. Amphipods and isopods were identified as far as practicable, based on whole specimen morphology (dissection was avoided to retain material for DNA analyses) using existing literature and guides to the New Zealand stygofauna (Appendix 2).

| Physical and chemical parameters
Two 250 ml water samples were collected in acid-washed bottles from each well between pumping and plankton net sampling to determine dissolved organic carbon content and nutrient concentrations, respectively (nutrients were dissolved reactive phosphorus (DRP), nitrite-nitrogen (NO 2 -N), nitrate-nitrogen (NO 3 -N), ammoniacal nitrogen (NH 4 -N), and total dissolved nitrogen (TDN) and total dissolved phosphorus (TDP)). Dissolved oxygen concentration was measured in a five-liter container of gently pumped well water using a TPS WP-82 meter (TPS Pty Ltd, Brisbane), and conductivity (µS/ cm), temperature, and pH were measured in situ using a TPS WP-81 meter (TPS Pty Ltd, Brisbane). Well depth and water column depth were measured in situ and information on well diameter and casing material extracted from local council databases.

| DNA analyses
Individuals were photographed and loaded into single wells on 96well microplates for processing at the Canadian Centre for DNA Barcoding (CCDB). Total DNA was extracted from specimens using a glass fiber plate method (Ivanova et al., 2006. Following DNA extraction, residual cuticular material for each specimen was depos- Polymerase chain reaction (PCR) amplification of the mitochondrial cytochrome c oxidase subunit I (COI) gene region used the primer pairs LepF1 and LepR1 (Hebert et al., 2004) and LCO490 and HCO2198 (Folmer et al., 1994) according to CCDB standard protocols (Ivanova & Grainger, 2007). Successfully amplified products progressed to cycle sequencing using BigDye ™ v3.1 terminator chemistry (Applied Biosystems ™ ). Products were then cleaned using a semi-automated AutoDTR ™ method (EdgeBio ® ) before being sequenced in forward and reverse directions on an ABI 3730xl DNA Analyzer (Applied Biosystems ™ ) using the same primers used for PCR amplification.
Specimen images, collection data, raw trace files, and edited sequences were all uploaded to and are available on the Barcode of Life Datasystems (BOLD) database  (http://dx.doi.org/10.5883/DS-GDWMS) and cross-referenced to GenBank (accession numbers OK072722-OK072875). Barcode Index Numbers (BINs; Ratnasingham & Hebert, 2013) assigned by BOLD were used to delineate putative species based on the sequence data (Milton et al., 2013).

| Data processing and statistical analyses
The amphipod sequences were aligned in Geneious Prime 2020.0.4 (Kearse et al., 2012, https://www.genei ous.com) using MUSCLE (2015)). Note that wells in the Moutere Valley aquifer were assigned to the Motueka River Terraces aquifer due to a surface water connection nearby (the Motueka River). Waimea is a composite label for multiple aquifers that have some degree of hydrological connectivity and Mataura is an unofficial aquifer name generated for this project. See Appendices 1 and 2 for well details (Edgar, 2004) and trimmed to 462 bp. A Maximum-likelihood (ML) phylogenetic tree was generated in MEGA7, (Kumar et al., 2016) with We calculated species accumulation curves using the abundance of BINs within each well to investigate how sampling effort (number of wells) affected the diversity of BINs both at the national and at regional scales. If most species within an area are collected, an asymptote in cumulative species richness is expected as sample number increases. The mean, median, and variance (quartiles and range) of species richness estimates for each additional well (1-n wells with genetic sequences) were calculated from 100 permutations with wells added in random order using the package Vegan in the statistical program R (R Development Core Team, 2017).

F I G U R E 2 Well locations across four regions of New Zealand with labeled aquifers (from Ministry for the Environment
Spearman rank correlation was used to assess whether the diversity of BINs detected within a well was positively correlated with the number of sequenced individuals. To investigate potential range restriction and, thus, endemism at different scales, the spatial occurrence of individual BINs between aquifers, regions, and islands was also assessed. Separate one-way analyses of variance (ANOVA) were used to assess whether wells that contained specimens with successful COI sequences differed in physical and chemical parameters from wells that either contained no stygofauna or from which successful sequences were not generated. These parameters included well depth, spot measurements of water temperature, conductivity and pH, and nutrient concentrations (DRP, TDN, TDP, nitrate-N, nitrite-N, ammoniacal-N). For each assigned aquifer, the diversity of amphipods per 1000 km 2 of surface catchment area (from Booker & Whitehead, 2017) and aquifer area (from Ministry for the Environment, 2015) was calculated.

| RE SULTS
Of the 65 wells sampled, 54 (83%) contained stygofauna with amphipods found in 34 wells (52%) and isopods in 15 wells (23%) ( Table 1). All amphipod and isopod specimens were considered to be stygobionts as they lacked body pigments or pigmented eyes (c.f., Marmonier et al., 1993). Other taxa that were found in more than five wells included cyclopoid and harpacticoid copepods, Syncarida, Ostracoda, Acarina, Annelida, Nematoda, and Gastropoda. From the wells containing amphipods and isopods, 186 amphipods and 42 isopods were collected and processed for their COI sequences, with successful sequences obtained from 154 individuals (68% overall success rate). Morphological identification to family or occasionally genus was possible for most specimens. Mounting and dissection of specimens for morphological assessment were largely precluded as many of the specimens were very small (e.g., adult, brooding amphipods <2 mm long) or their often damaged condition meant that any available tissue was required for DNA extraction. However, subsequent re-examination of morphologies following COI sequencing (based on their clustering within trees; Figures 3 and 4), allowed us to allocate further specimens to established families or genera. TA B L E 1 Numbers of wells within the four New Zealand regions and seven aquifers that were sampled for stygofauna, in which stygofauna (all stygofauna, amphipods, and isopods) were collected, and from which COI gene sequences were successfully obtained (amphipods, isopods, and combined amphipods and isopods) Note: Aquifer names are modified from Ministry for the Environment (2015). * refers to complex of multiple aquifers; ** "Mataura" is the name of a nearby river and not an official aquifer name.
Of the 228 individual crustaceans collected, successful sequences were obtained from 129 amphipods (69% success) and 25 isopods (60% success). These sequences were obtained from a total of 25 wells, with amphipods and isopods successfully sequenced from 23 and 11 wells, respectively ( Table 1). Seventeen of the wells were in Canterbury, with 5-6 wells in each of the aquifers (Table 1). Although sequences from the Canterbury region were over-represented, most aquifers, apart from the Ruataniwha aquifer in Hawke's Bay, contained at least two wells where sequences were obtained (Table 1). Individually, amphipod taxa were successfully sequenced from at least one well in each aquifer. However, isopod sequences were available only for the Canterbury and Tasman aquifers (Table 1).
The 154 COI sequences were assigned to 45 BINs, comprised of nine isopod BINs and 36 amphipod BINs ( Table 2). The number of specimens available for sequencing and success of sequencing differed between regions, catchments, and wells within catchments ( Species accumulation curves (based on BINs) were unsaturated both for the Canterbury region and for all wells combined ( Figure 3).
Thirty-five amphipod BINs (78% of all taxa and 83% of amphipods) were new records on BOLD. Of these, only three could be attributed to known morphologically described genera. The endemic amphipod genus Ringanui was assigned to three BINs ( Figure 3).
The Canterbury/Rangitata Levels BIN (ADL2688, Figure 3) probably comprised one of the two described species (reported range Waimakariri-Ashley to Rangitata Levels aquifers and Temuka; Fenwick, 2006), whereas two other BINs (ADL5144, ADL5178) are probably undescribed species endemic to the Waimea aquifer.
Substantial morphologically cryptic diversity was identified at the family level. The Family Paraleptamphopidae includes three described genera, two of which are hypogean, whereas our analysis found 27 BINs representing several potential genera within a large paraleptamphopid clade (Figure 3).
Cryptic diversity was also apparent within the Isopoda. Eight specimens that were originally morphologically identified as Cruregens fontanus were assigned to five well-separated BINs (>92% support, Figure 5), one each in Motueka and Waimea aquifers (BINs ADP0923, 3149), one shared between aquifers within the Canterbury region (ADL3492) and two appear to be single-aquifer endemics (ADL2602, ADP4594) within the Central Plains aquifer ( Figure 5). Similarly, three BINs of the phreatoicid isopods were distinguished (>91% support) from specimens initially identified as Phreatoicus typicus and P. orarii ( Figure 5).
Each of the 45 genetically distinct isopod and amphipod BINs was restricted to one region; 27 BINs were unique to Canterbury, 16 to Tasman, one to Hawke's Bay, and one to Southland (  Figure 4). Most isopod BINs (89%) were also from single aquifers while specimens assigned to one BIN (ADL3492) were collected in two aquifers (Central Plains, six specimens; Waimakariri-Ashley, two specimens) (Figure 4).
Of the 39 BINs (87%) restricted to individual aquifers, 29 (64% of all BINs) were found at one well within the aquifer, and nine (20%) BINs confined to a single aquifer were found in two wells. Three BINs occurred at three wells within an aquifer.
Increasing spatial separation was commonly correlated with greater genetic variability. For example, specimens of the Family Paracrangonyctidae from the Motueka aquifer (ADL3783) were genetically distinct from Waimea specimens (ADL2540) (Figure 4), and those from Canterbury's Waimakariri-Ashley system (ADL5568) were even more genetically divergent, reflecting the much greater geographic distance of the Waimakariri-Ashley system from the two Tasman aquifers (c. 150 km cf. <1 km between adjacent headwater tributaries).
We found very high amphipod diversity in four aquifers within two of the regions studied. Comparisons based on estimated numbers of BINs per 1000 km 2 of aquifer area and catchment area (

| Physical and chemical parameters
The sixty-five sampled wells varied in physical size (depth range: 2.7-39 m, diameter range: 50-1200 mm), chemical parameters (e.g., conductivity range 1.2-1014 μS), and nutrient status (e.g., NO 3 -N range 1.0-11,000 mg/m 3 , Appendix 1). The water in wells from which specimens with successful genetic sequences were collected was cooler (median spot water temperature 12.9°C) than wells where either taxa were not collected or genetic sequencing failed (median water temperature 13.5°C, ANOVA:
There were no differences in well depth or nutrient concentrations between wells that had specimens resulting in successful COI sequences compared with wells where no successful sequences were obtained, or from which no amphipods or isopods were collected.
Seventy percent (n = 18) of the wells from which successful COI sequences for amphipods and isopods were processed had steel casings, while six wells had PVC casings and two were undetermined.

| D ISCUSS I ON
Of the 45 putative species (BINs) identified from the COI sequences, 78% were previously unrecorded on BOLD. Of these, only three could be attributed to established genera, indicating that current knowledge of New Zealand's stygofaunal diversity is extremely low.
Morphologically cryptic taxa were common, as has been found in other genetic studies of groundwater taxa (e.g., Delić et al., 2017;Eme et al., 2018). For example, one currently recognized isopod species (Cruregens fontanus) was assigned to five well-separated BINs, and over 20 species were found within the Family Paraleptamphopidae, particularly within the genus Paraleptamphopus. Six specimens of phreatoicid isopods also showed cryptic diversity. Specifically, three BINs were found in the vicinity of the Central Plains aquifer, whereas Chilton (1894), and Wilson and Fenwick (1999) previously reported a single species within the Central Plains and Waimakariri-Ashley aquifers.
The few examples of taxa found in more than one aquifer were possibly stygophilic, migrating between aquifers via permanent and/or intermittent surface water connections.
Endemism appeared high, with all species found in only one region, 87% attributed to single aquifers and more than 50% recovered only from single wells. However, we caution that the actual levels of local endemism (e.g., well, aquifer) are likely to be somewhat lower, as our sampling within individual wells and aquifers was not comprehensive (we obtained successful sequences for isopods and amphipods from 38% (n = 25) of the 65 wells sampled). Due to their subterranean nature, the sampling of groundwater ecosystems is inherently challenging (Hancock & Boulton, 2009;Korbel et al., 2017).
Specifically, collection is often restricted to pre-existing wells that are predominantly located in areas of human activities. Likewise, it can be difficult to effectively deploy sampling equipment within a well and samples may not adequately represent invertebrate biodiversity within the surrounding aquifer (Ficetola et al., 2019;Larned, 2012), particularly when sampling is only possible on a single occasion. For example, in the Pilbara region of Australia, multiple sampling methods and multiple visits were required to capture most of the species present within a given well (Eberhard et al., 2009).
In our study, only aquifers that were better sampled (i.e., more wells and more sequenced individuals), such as the Canterbury aquifers, yielded taxa that were not aquifer specific. This implies that further sampling is likely to increase the geographic range of some of the apparent single-aquifer endemic taxa. However, between 42% and 71% of taxa in our most intensively sampled region (Canterbury) were aquifer-specific indicating that additional sampling would also reveal additional species, many of which are likely to be range-restricted. In the United States, sampling over nearly 40 years after an early survey of cave fauna led to a <20% decline in frequency of reported county-specific endemism, while the absolute number of endemic species increased nearly threefold . While we are unable to definitively identify aquifer-specific taxa based on our study, other studies of groundwater stygofauna indicate that aquifer-specific endemism is likely to be common (e.g., Culver & Sket, 2000;Ferreira et al., 2007;Bradford et al., 2010;Murphy et al., 2013 in France, after 200 years of study, stands at 380 species, although this is likely to be an under-estimate due to incomplete sampling (Ferreira et al., 2007) and the presence of several morphologically cryptic species (e.g., Wattier et al., 2020;Westram et al., 2011). In the United States, 300 species of cave-dwelling aquatic groundwater species are known . Extrapolating sampling effort and species caught suggests that the Pilbara region of Western Australia may contain 500-550 species (Eberhard et al., 2009) and 21 species of amphipods were present in Australia's smaller Yilgarn region (Cooper et al., 2007).
We found high richness of stygofaunal amphipods in three of the more intensively sampled aquifers (>6.0 BINs/1000 km 2 ).
We also estimated extremely high stygofaunal amphipod richness (>80 BINs/1000 km 2 ) within the Waimea aquifer, probably resulting from its interconnectedness with two other aquifers and the extremely high hydraulic transmissivity (20,000 m 2 /day) reported for parts of this aquifer, including adjacent to the rivers (White & Rosen, 2001). These amphipod richness values are similar to the highest reported for total stygofaunas elsewhere: 6.6 total stygofaunal species/1000 km 2 for karst in the Balkan Peninsula , although those measures of richness preceded DNA investigations which revealed substantial cryptic stygofaunal diversity. The New Zealand stygofaunal amphipod richness is twice that reported for the total stygofauna in the Pilbara region of Australia (or 3.1 species/1000 km 2 ) and much greater than that reported for amphipods from Australia's Yilgarn region (0.05 species/1000 km 2 ; Cooper et al., 2007). The higher richness found in these New Zealand aquifers is more remarkable because it does not include other taxa such as isopods, copepods, ostracods, syncarids, platyhelminths, and oligochaetes known from these New Zealand aquifers (e.g., Fenwick, 2000;Larned et al., 2014;. As with stygofaunas elsewhere, many of described and new species of amphipods and isopods inhabiting New Zealand's alluvial groundwater likely have restricted geographic distributions. The number, size, and complexity of New Zealand's aquifer systems, hydrologically separated by extensive hills and mountains, are probable reasons for the country's high stygofaunal diversity. This appears true even where headwater tributaries almost join (e.g., Motueka and Waimea catchments). Stygofaunal populations also appear genetically isolated between aquifers within the relatively homogeneous landscape of Canterbury Plains, where there are no obvious geohydrological barriers. However, the complexity of the plains' subsurface geology and hydrogeology (Bradshaw & Soons, 2008;Davey, 2006) may hydrologically separate individual aquifers, leading to genetic isolation and at least some, short-range, endemism.
Climatic events appear to be the main factor in changing hydrological connectivity and genetic isolation for tectonically stable continents like Europe, North America, Australia, and Africa (Collins et al., 2019;King & Leys, 2014;Lefébure et al., 2006;Witt et al., 2006). Both glaciations and aridity are strongly implicated in creating hydrological barriers that isolated populations of stygofauna, leading to genetic divergence and speciation. These two types of climatic events are also likely to be important drivers of hydrological and genetic isolation of aquifers within New Zealand's alluvial plain systems.
Tectonic events may also have a role via lateral and/or vertical displacement creating barriers and/or changing groundwater flow directions (Trontelj et al., 2009;Craw & Waters, 2007). For example, most (if not all) of New Zealand's larger plain systems are fragmented with recent and historic faults, including the Pacific-Australian tectonic plate boundary. The effect of active faulting (e.g., the Greendale Fault responsible for the 2011 Christchurch eartquake), on dispersal and gene flow is unknown, although it may create physical barriers within an aquifer. Shearing and shaking could consolidate alluvium, reduce interstitial pore spaces and hydrological connectivity, or uplift may misalign strata to subdivide an aquifer (Cox et al., 2012;Rutter et al., 2016). Alternatively, tectonic activity may breach existing hydrological barriers between adjacent aquifers (e.g., bedrock fractures through ranges or breaks in confining layers may create new hydrological connections) and facilitate gene flow.

| Summary
There have been repeated calls for accelerated scientific work to identify groundwater biodiversity, which is threatened with extinction before being discovered, identified, and ideally assigned a conservation status and protected (Gladstone et al., 2021;Mammola et al., 2019). Like most countries, knowledge of groundwater fauna is exceptionally poor in New Zealand. Our results support common findings of high biodiversity and short-range endemism in groundwater faunas internationally (e.g., Boulton, 2020;Gladstone et al., 2021) and likewise for the use of genetic data in identifying morphologically cryptic species, which are common in groundwaters (Boulton, 2020;Delic et al., 2017;Eme et al., 2018;Gladstone et al., 2021). By contributing to knowledge of the biodiversity and spatial distributions of groundwater taxa, we hope to help address some of the knowledge gaps inhibiting conservation of groundwater biodiversity (e.g., Boulton, 2020;Mammola et al., 2019Mammola et al., , 2020.

ACK N OWLED G M ENTS
We are extremely grateful to S. Mammola

CO N FLI C T O F I NTE R E S T
The authors declare no conflicting interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Residual cuticular material after DNA extraction for each speci- All data associated with the paper are also available from the NZ Biological Heritage National Science Challenge data repository at https://doi.org/10.34721/ qk3f-6y64.

A PPE N D I X 1
Summary environmental characteristics for the 65 sampled wells. The number of sites in each category is listed in the median column for categorical variables. Some wells were missing environmental data owing to logistical sampling constraints and equipment failure

A PPE N D I X 2
The following keys are available and were used to morphologically identify amphipods and isopods: