Taking eDNA underground: Factors affecting eDNA detection of subterranean fauna in groundwater

Stygofauna are aquatic fauna that have evolved to live underground. The impacts of anthropogenic climate change, extraction and pollution on groundwater pose major threats to groundwater health, prompting the need for efficient and reliable means to detect and monitor stygofaunal communities. Conventional survey techniques for these species rely on morphological identification and can be biased, labour‐intensive and often indeterminate to lower taxonomic levels. By contrast, environmental DNA (eDNA)‐based methods have the potential to dramatically improve on existing stygofaunal survey methods in a large range of habitats and for all life stages, reducing the need for the destructive manual collection of often critically endangered species or for specialized taxonomic expertise. We compared eDNA and haul‐net samples collected in 2020 and 2021 from 19 groundwater bores and a cave on Barrow Island, northwest Western Australia, and assessed how sampling factors influenced the quality of eDNA detection of stygofauna. The two detection methods were complementary; eDNA metabarcoding was able to detect soft‐bodied taxa and fish often missed by nets, but only detected seven of the nine stygofaunal crustacean orders identified from haul‐net specimens. Our results also indicated that eDNA metabarcoding could detect 54%–100% of stygofauna from shallow‐water samples and 82%–90% from sediment samples. However, there was significant variation in stygofaunal diversity between sample years and sampling types. The findings of this study demonstrate that haul‐net sampling has a tendency to underestimate stygofaunal diversity and that eDNA metabarcoding of groundwater can substantially improve the efficiency of stygofaunal surveys.


| INTRODUC TI ON
Groundwater comprises over 90% of the available freshwater on Earth and is a key resource for surface life (i.e., natural ecosystems, human consumption, agriculture, industry and urbanization; Danielopol et al., 2003). Groundwater also supports diverse, complex ecosystems of subterranean fauna, largely comprising shortrange endemic invertebrate and vertebrate species that represent exceptional levels of species richness and often ancient evolutionary lineages reflective of the past climatic and biogeographical history (Des Châtelliers et al., 2009;Fenwick et al., 2021;Juan et al., 2010;Ribera et al., 2010). Stygofauna (aquatic subterranean fauna) perform critical water purification services; they consume and stimulate microbial processes and biomasses which are essential to the ecological integrity of groundwater systems (Smith et al., 2016), and can facilitate hydraulic connectivity. Concerns over the impact of humanmediated disturbances (i.e., aquifer drawdown, contamination and climate change; Mammola et al., 2019) on groundwater-dependent communities has led to an increased need to characterize, describe and monitor them (Saccò, Guzik, et al., 2022;Stumpp & Hose, 2013).
Except for a few examples-Switzerland (SWPO, 1998) and Australia (EPA (WA), 2003; NSW-SGDEP, 2002)-predicting and monitoring the impacts of proposed developments on subterranean fauna are still rarely required within a regulatory setting (Danielopol & Griebler, 2007;Niemiller, Taylor, & Bichuette, 2018). In part, this is because several major impediments still exist that limit the ability to monitor and assess the diversity and distribution of stygofauna species (Cardoso et al., 2011;Mammola et al., 2019), including (i) specialist taxonomic expertise required to morphologically identify the small and highly adaptive forms of stygofauna (Bork et al., 2008;Dumas & Fontanini, 2001); (ii) DNA barcoding is often used to verify morphological species identification due to the high numbers of unknown and cryptic species ; (iii) subterranean environments are difficult to access and, in most cases, sampling is often limited to low-yield methods that require specialist equipment (Saccò, Guzik, et al., 2022); (iv) visual observation and assessment is extremely challenging at subterranean sampling access points such as caves, springs and boreholes (see Section 2 for details; Sorensen et al., 2013). There is therefore a need for new methods that will increase the accessibility, efficiency and accuracy of subterranean monitoring (Gibson et al., 2019;Mammola et al., 2022;Saccò, Guzik, et al., 2022).
Environmental DNA is genetic material that has been shed into the environment by organisms, and can be collected by extracting DNA from various substrates including water and sediment (Taberlet et al., 2012). Research in marine systems, combining next generation sequencing (NGS) and eDNA, has demonstrated that eDNA metabarcoding is an important new molecular tool for biodiversity assessment (Miya et al., 2015;Stat et al., 2017). In other ecosystems, eDNA-based monitoring can significantly decrease the time, expense and resources required compared to traditional survey methods in other ecosystems (Fediajevaite et al., 2021), and has the potential to improve on existing stygofauna sampling and monitoring approaches. Improvements have already been observed; recent studies have successfully used single-species assays to detect rare and endangered subterranean species (Boyd et al., 2020;Gorički et al., 2017;Niemiller, Porter, et al., 2018;Vörös et al., 2017) and investigate genetic differences (haplotypes) between populations . Additionally, Saccò, Blyth, Humphreys, Karasiewicz, et al. (2020), Saccò, Guzik, et al., 2022) have highlighted the potential value of eDNA for community-level detection of subterranean fauna and their diets (Korbel et al., 2017;Saccò, Campbell, et al., 2022;West et al., 2020) in a large range of habitats and for all life stages, thereby avoiding the need for: (i) specialized taxonomic expertise (Bohmann et al., 2014;Thomsen & Willerslev, 2015); and (ii) destructive sampling for identification. However, before an eDNA metabarcoding approach can be used as a standard method for monitoring and assessment of groundwater communities, it is essential that the appropriate testing and development of the methods are undertaken to assess the efficacy, accuracy and speed of species detection. Here we outline some ecosystem-specific factors that may impact the development of eDNA protocols for subterranean habitats.
To date, the utility of eDNA-based surveys in subterranean habitats has been undertaken for the detection of priority species (Baker et al., 2020;White et al., 2020), stygobitic communities from caves (West et al., 2020), and from boreholes and springs (Korbel et al., 2017;Oberprieler et al., 2021). Common challenges identified from eDNA studies for subterranean ecosystems include a critical lack of reference sequence availability and biases in metabarcoding assays (Korbel et al., 2017;Oberprieler et al., 2021). Further, the most appropriate sampling protocols for underground systems have yet to be formally tested and standardized. For instance, sample type (water, sediment, etc.), volume and depth are known to influence the diversity detected through eDNA (Jeunen et al., 2020;Koziol et al., 2018;West et al., 2020). Conversely, it is unknown if bore characteristics (e.g., diameter) or order of sampling (i.e., performing haul-net sampling first) will affect the species diversity that can be detected from eDNA.
Environmental factors are known to influence the detection and persistence of eDNA (e.g., high temperatures promote DNA degradation through enzymatic activity and denaturation, while anoxic environments limit degradation and stabilize eDNA; Barnes et al., 2014). However, for stygofauna, ecological factors, such as life history strategies, niche occupation and possible seasonal/temporal variability that might affect eDNA shedding rates and longevity, are poorly known. We know that there is a general pattern of decreasing species richness with depth, attributed to lower levels of oxygen and nutrients in deeper water (Brunke & Gonser, 1999;Datry et al., 2005). Further, previous studies have indicated that groundwater recharge can influence the abundance and life-history stages of some stygofauna (Datry et al., 2005;Hyde et al., 2018;Saccò, Blyth, Humphreys, Karasiewicz, et al., 2020) probably through increases in organic carbon leaching into the groundwater (Datry et al., 2005;Herrmann et al., 2020;Saccò, Blyth, Humphreys, Middleton, et al., 2020). However, recharge rates can be quite low (Crosbie et al., 2010), with the groundwater in many parts of Australia being tens of thousands, if not millions, of years old (Meredith, 2009). There is much developmental work required to establish protocols for subterranean environments.
Barrow Island, located ~55 km off the Pilbara coast of Western Australia, is the largest island in the tropical waters of the northwest continental shelf of Western Australia (Moro & Lagdon, 2013). At 202 km 2 the island has been listed as a Class A nature reserve since 1910 for its extreme short-range endemic mammal, herpetological and subterranean fauna, including stygofauna and troglofauna (terrestrial subterranean fauna). The island also supports resource infrastructure including an operating oilfield (since 1960) and natural gas plant (Chevron, 2006). The subterranean biodiversity value of Barrow Island is internationally recognized (Burbidge, 2004) with an unprecedented 63 stygobitic species representing 12 orders .
Intensive sampling of groundwater has occurred since 2001, through Environmental Impact Assessment and monitoring in accordance with Western Australia's EPA guidelines (Eberhard et al., 2005;EPA, 2003;Humphreys et al., 2013). Given the high conservation status, the potential impacts of resource development and the legislated requirement to monitor for environmental impacts on Barrow Island, there is a real need to detect stygofauna as accurately as possible.
Here, we used the groundwater ecosystem at Barrow Island as an exemplar for understanding the factors that may affect eDNA detection of a stygofaunal community from groundwater. In addition to its distinct ecosystem, Barrow Island is also suited for scientific study due to the extensive number of sampling bores, many of which have been purpose-drilled to access the groundwater and stygofauna community for monitoring purposes. These bores have been cased and capped to retain the integrity of the bore and its groundwater profile. Using key parameters that are known to impact eDNA detection in various environments (i.e., substrate; Jeunen et al., 2020;Koziol et al., 2018) we investigated factors important in developing a robust eDNA sampling approach for effective assessment and monitoring of stygofaunal communities. Our aims were as follows: 1. To compare eDNA metabarcoding with a universal eukaryote assay (18S) to conventional stygofaunal survey methods (i.e., whole specimens were collected and identified; see Methods below). This comparison can be used to determine the strengths and limitations of both methods and the degree to which they complement each other. This comparison includes stygofauna and troglofauna which can also be detected using both methods.
2. To investigate whether substrate type (shallow water, deep water, sediment) and sampling protocol (order of sampling) affect the stygofaunal diversity and community composition detected through eDNA metabarcoding. The degree of overlap between substrate types will determine whether it is necessary to sample multiple substrates to detect the most diversity.
3. To examine whether temporal differences affect stygofauna communities detected through haul-nets and eDNA. Dynamism in the system and/or stochasticity would indicate the need for multiple sampling events to monitor the system.

| Study site
Barrow and Boodie islands are located northwest of mainland Western Australia on the edge of the Australian continental shelf.
The climate of the islands is arid and subtropical with hot, humid summers with daily mean maxima reaching 34°C and daily mean minima of 20°C. The winter daily mean maxima are lower at 26°C and daily mean minima of 17°C (Moro & Lagdon, 2013). Rainfall is seasonal, most prevalent in November-April, and dependent on the passage of summer monsoons. The average annual rainfall has historically been 320 mm; however, the first year of sampling for this study (2020) was during a particularly dry period when the average annual rain for the previous 3 years was only 105 mm (www.bom.gov.au). The following year of sampling (2021) was during a wetter year with 377 mm of rain, the majority occurring between March and May.
Groundwater on Barrow Island comprises an exceptional example of an anchialine system (Humphreys, 2001); that is, stratified waters where a freshwater lens that originates from seasonal rainfall overlies seawater with a transitional zone in between. The transitional zone is of comparable thickness to the freshwater layer having been expanded due to tidal forces. As a smaller island just south of Barrow, Boodie Island has a slimmer freshwater lens and is almost entirely marine. Access to the groundwater and in situ stygofaunal sampling is conducted via existing bores that have been installed for specific operational or monitoring purposes. The bores are typically constructed of 50-or 100-mm-diameter, PVC casing that is slotted at discrete intervals depending on their purpose. Slots vary from 0.5 to 5 mm in width. Bores with larger slot sizes can have mesh across them to prevent debris from accumulating inside the bore. Depth of bores varies depending on the location on the island (bores at higher elevation need to be deeper to reach the water column) and the initial purpose of the bore. Bores designed to monitor stygofauna or groundwater have slotted intervals in shallower parts of the aquifer (i.e., water several metres deep below the water table) compared to those constructed for the cathodic protection system (tens of metres below the water table). Many bores are slotted below the water table.

| Sample collection and processing
Stygofaunal specimens and eDNA samples were collected from 17 bores in December 2020 and November 2021 ( Figure 1). An additional two bores were sampled in 2021 (November) from Boodie Island. Groundwater at the base of a cave on the west side of Barrow Island was also sampled in November 2021. The cave and Boodie Island sites are particularly difficult to access and were unavailable for sampling during the 2020 sampling period. Bores chosen for sampling were spatially stratified across the island and high in stygofaunal diversity according to previous monitoring efforts . Two "control" bores (BGMW08, WFS2MW04) with no previous detection of stygofauna were also selected.

| Stygofaunal collection and processing
Whole specimens used to assess a baseline for morphological species diversity were collected using haul-nets (similar to plankton-netssee Saccò, Campbell, et al., 2022;Saccò, Guzik, et al., 2022), which were 45 or 95 mm in diameter, depending on the size of the PVC casing. Five (Allford et al., 2008) to six (EPA guidance) haul-net samples were taken at each bore, half with 150μm mesh and half with 50μm mesh, alternating mesh size between hauls. We favoured haul-net sampling for stygofauna prior to collection of eDNA samples to ensure that live specimens were not dispersed by the eDNA sampling equipment. Haul-net samples were sorted and identified using a dissecting microscope and specimens of each morphotype were stored in glass vials in 100% ethanol and kept frozen at −20°C. Stygofauna experts (N.S. and M.T.G.) provided the morphological identifications.

| eDNA sample collection
Water samples were collected for eDNA analysis after hauling. Four 1-L water samples were collected from the upper 2 m of the water column using a 1-L plastic disposable bailer (shallow water samples).
One 1-L water sample was collected from the bottom of the bore using a 1-L steel bailer (deeper water sample). Sediment samples (one per bore) were collected from the deeper eDNA water samples where there was an abundance of sediment left over that could not be filtered. With six samples per bore (i.e., 4× 1-L eDNA, 1× 1-L deeper, 1× sediment), there were 102 samples collected in 2020, and 120 samples in 2021. For five bores the scope was extended to include collecting the eDNA samples before and after the haul-net samples to determine if the order of sampling was significant. This resulted in 12 eDNA samples per bore for these five bores with a total of 30 samples collected before nets and 30 after. All equipment was sterilized between bores in 10% bleach for 10 min and rinsed in reverse osmosis (RO) water, and disposable gloves were used at each borehole to prevent contamination. Fresh bleach solution and RO water was used for every borehole for equipment that was not decontaminated in the laboratory ahead of time, and samples of the RO water used for rinsing equipment (i.e., field controls) were collected to assess sources of contamination. Additional quality control samples included: two blank filter membranes carried in the field and a sample of all RO water sources for a total of 34 field controls. eDNA and haul-net samples were kept on ice until they could be transferred to a 4°C refrigerator at the end of the collecting day. Water physicochemical properties were measured for each borehole using a YSI ProDSS either 2 days before or the day after sampling, to avoid cross-contamination between boreholes.
Measurements of turbidity (FNU), conductivity (μS cm −1 ), salinity (psu), temperature (°C), dissolved oxygen (%sat, mg L −1 ), pH, and depth (m) were recorded for each borehole. Water chemistry at various depths was only measured for bores with a water column that was greater than several metres.

| DNA barcoding of specimens and custom reference library
A leg was removed from 15 morphologically identified specimens (Table S1) collected via haul-nets and DNA was extracted using the DNeasy Blood and Tissue kit (Qiagen) following the manufacturer's instructions with the following modification: genomic DNA was eluted from the silica column with 50 μL of AE buffer. The DNA was then amplified using PCR with universal eukaryote18S primers (Table S2; Pochon et al., 2013). These primers have previously been used for metabarcoding assays (~400 bp) to detect stygofauna (West et al., 2020) and broad enough to detect a wide range of organisms, which might otherwise require multiple assays (Pochon et al., 2013). Each PCR for barcoding was carried out in duplicate in a total volume of 25 μL containing: 1× AmpliTaq Gold PCR buffer (Life Technologies), 2.5 mm MgCl2, 0.25 mm dNTPs, 0.4 μm each of forward and reverse primers (Integrated DNA Technologies), 0.4 mg mL −1 bovine serum albumin (Fisher Biotec), 0.6 μL of 5× SYBR Green (Life Technologies), 1 U AmpliTaq Gold DNA Polymerase (Life Technologies), 2 μL of genomic DNA template, and made up to volume with Ultrapure Distilled Water (Life Technologies). The cycling conditions were initial denaturation at 95°C for 5 min, followed by 40-45 cycles of 95°C for 30 s, 52°C for 30 s, 72°C for 45 s, a melt curve stage of 95°C for 15 s, 60°C for 1 min and 95°C for 15 s, and a final extension at 72°C for 10 min. The expected PCR product size (~400 bp) was visualized on a 2% agarose gel and sent to the Australian Genome Research Facility (Perth) for Sanger sequencing in both directions. The forward and reverse sequences were then aligned, the primers were removed and the consensus sequences were used to create a custom library for reference analyses (i.e., taxonomic identification) of the 18S eDNA sequences obtained from the water and sediment samples. All barcoding sequences generated were submitted to GenBank (Table S1). The custom DNA sequence library also included sequences of stygofauna from the Pilbara region sourced from BOLD and GenBank (www.ncbi.nlm. nih.gov/genbank), as well as previous consulting work done in the area for a total of 207 18S sequences. DNA extracts were screened for quality and quantity of DNA using qPCR to determine the presence of inhibitors and the quantity of target template molecules present in each DNA extract (Murray et al., 2015). Three dilutions (neat, 1:5 and 1:10) were used for a subset of 30 samples, while the remaining were assessed using only the neat and 1:10 dilution. Each PCR for the 18S universal assay was carried out with the same master mix as described above. The cycling conditions were initial denaturation at 95°C for 5 min, followed by 45 cycles of 95°C for 30 s, annealing at 52°C for 30 s, 72°C for 45 s, a melt curve stage of 95°C for 15 s, 60°C for 1 min and 95°C for 15 s, and a final extension at 72°C for 10 min. All PCR plates included a negative control and a positive control of chicken DNA.

| eDNA laboratory processing
Each sample, including all controls, was assigned a unique combination of multiplex identifier (MID) tags for the 18S universal assay. These MID tags were incorporated into fusion tagged primers, and none of the primer-MID tag combinations had been used previously in the laboratory to prevent cross-contamination. Fusion PCRs were done in duplicate to minimize PCR stochasticity, and the mixes were prepared in a dedicated ultraclean room before DNA was added. The PCRs were done with the same conditions as the standard qPCRs described above. Samples were then pooled into approximately equimolar concentrations to produce a PCR amplicon library that was size-selected to remove any primer-dimer that may have accumulated during fusion PCR. Size selection was performed (250-600 bp) using a PippinPrep 2% ethidium bromide cassette (Sage Science). Libraries were cleaned using a QIAquick PCR Purification Kit (Qiagen) and quantified using Qubit Fluorometric Quantitation (Thermo Fisher Scientific). Sequencing was performed on the Illumina MiSeq platform using the 500-cycle paired-end V2 as per the manufacturer's instructions.

| Bioinformatics
DNA sequences were processed using eDNAFlow, an automated bioinformatics workflow designed for analysis of eDNA metabarcoding data (Mousavi-Derazmahalleh et al., 2021). This pipeline uses adaptorremoval (Schubert et al., 2016) and fastqc (Andrews, 2010) to quality filter and demultiplex sequences using "obitools" (Boyer et al., 2016). The minimum Phred quality score was set at 20, with a minimum alignment of 12 for the paired-end reads, minimum length of 100 bp and maximum two primer mismatches allowed, but no mismatches allowed in the MID-tag sequences.
The demultiplexed sequences were then denoised using usearch unoise3 (Edgar, 2016) to create zero-radius operational taxonomic units (ZOTUs) with a minimum abundance of 8, the recommended value for the denoising algorithm (Edgar, 2016). The ZOTU sequences were then queried against a custom barcode reference library for stygofauna and GenBank (NCBI) using blastn with the following parameters (−outfmt "5"-perc_identity 95-qcov_hsp_ perc 95-query-max_target_seqs 10). Taxonomic assignment was performed using a simple lowest common ancestor algorithm on megan (Huson et al., 2007) with a minimum score of 450.
This allowed us to parse stygofauna and troglofaunal taxa from nonsubterranean taxa. ZOTUs belonging to subterranean taxa were then extracted an aligned along with sequences from the custom reference database using muscle (Edgar et al., 2004) in geneious 2021.0.3 (https://www.genei ous.com) with 10 iterations of the default settings. The resulting tree was used to infer putative species with a Poisson tree processes (PTP) model (Zhang et al., 2013). Some subterranean ZOTU taxon names were modified to reflect the species inferences from the model such that all ZOTUs assigned to a species had the same taxa name, and no two subterranean "species" had the same name. For example, three ZOTUs identified as "Atyidae" were renamed to Stygiocaris stylifera because they belonged to the same putative species. Details on initial and modified taxon names can be found in the taxon table in the dryad digital repository.

| Statistical analysis
Statistical analyses were performed in R 4.0.5 (R Core Team, 2018).
We used the R package "phyloseq" (McMurdie & Holmes, 2013) to filter the data set. This included removing samples with fewer than 2000 reads, where rarefaction curves appeared to asymptote ( Figure S1) and pruning ZOTUs that were present in extraction controls and rinsate blanks (>5 reads). This removed 18 ZOTUs that were detected in extraction controls and 166 ZOTUs from the rinsate controls, although only three of those were stygofauna. Two of the potential stygofaunal ZOTUs removed were nematode ZOTUs identified as Monhysterida found in RO water used to rinse equipment. ZOTUs with fewer than five reads were removed, and abundances were converted to presence or absence. We selected only the ZOTUs from presumed subterranean taxa (both stygofaunal and troglofaunal), and combined the detections from all sample types for the comparison between subterranean taxa detected using morphological identification of haul-net samples, and taxa detected using eDNA. Differences in the taxa detected using each method were visualized using a bubble plot to compare taxa at the order level for both stygofauna and troglofauna.
The effect of substrate types (deeper, shallow, sediment) on stygofauna community composition was tested on a data set containing the one sediment sample, one deeper water sample and one randomly selected shallow water sample per borehole. The effect of sample type on stygofaunal ZOTU richness was tested using a Friedman rank sum test for repeated measures. Pairwise comparisons were performed using paired Wilcoxon sign rank tests with a Holm adjustment for multiple comparisons (Holm, 1979). For the data from 2021, the analysis was performed on the data for all sites, and without the three extra sites (Boodie Island and Ledge cave), to ensure the addition of those sites was not influencing the patterns observed. To examine which stygofaunal taxa could be detected using each sample type, the ZOTUs were collapsed by taxon using the tax_glom function from "phyloseq" (McMurdie & Holmes, 2013).
The taxa detected by each substrate were then illustrated using a Venn diagram function (ps_venn) from the "MicEco" package (Russel, 2021).
For the five bores where we collected eDNA samples both before and after using the haul-nets, we tested whether the order of sampling had an effect on stygofauna ZOTU richness using an analysis of variance (ANOVA) with "bore" as a random factor. Differences between community composition before and after haul-net sampling were likewise tested using a permutational multivariate analysis of variance (PERMANOVA, distance = "jaccard", permutations = 999) with bore as a random factor. Finally, a multipattern analysis was run using the R package "indicspecies" (Cáceres & Legendre, 2009) to determine which taxa were more prevalent in the sampling before or after haul-net samples.
An asymptotic regression curve was fitted to the accumulation curves for each sampling replicate using the "AR.2" function in the R package "drc" (Ritz et al., 2015). The replicates for each borehole were then merged and a Mantel test from the R package "ade4" (Dray & Dufour, 2007) was used to test the correlations between similarity in community composition and the distances between the bores. The clustering of bores was analysed using the SIMPROF (++Similarity Profile Analysis) function in the "CLUSTIG" package (Whitaker & Christman, 2014).

| Comparison of morphological vs. eDNA identification of subterranean fauna
Twenty-four subterranean taxa were morphologically identified from the net-haul samples taken at Barrow Island during 2020 and 2021. Of these, 17 were stygofaunal ( Table 1) and seven were troglofaunal. Using the 18S universal metabarcoding assay we identified 36 different subterranean taxa: 26 stygofaunal and 10 troglofaunal.
Here "taxa" refers to the lowest taxonomic identification available for the specimen, and all the ZOTUs assigned to the same "putative species" using the species delimitation model. Overall, more subterranean taxa were identified using an eDNA metabarcoding approach than with morphology alone (36 vs. 24; Table 1).
The stygofauna orders identified using a morphological approach were three copepod groups ( were detected (e.g., Hymenoptera, Coleoptera) in some bores, but these were incidental and the data cannot provide accurate comparison between methods for surface-dwelling taxa.
Some discordance between methods (eDNA vs. morphology) was observed at a few bore locations. For example, bore WFS2MW4, a "control" bore where no stygofauna have been detected on previous surveys (unpublished data from Chevron), was confirmed as a control bore using the haul-net sampling approach during both sampling periods (2020-2021). However, eDNA metabarcoding of water and sediment detected Cylcopoida (Copepoda) at this bore. Similarly, BGMW08 was previously identified as a "control" bore, although in 2020 a copepod taxon (Harpaticoida) was identified from haul-net samples and in 2021 no specimens were collected. However, the eDNA sampling at BGMW08 detected several nonarthropod invertebrates or "worm" taxa (Catenulida and Haplotaxida), but not copepods. Another bore that revealed inconsistent results was MW34, which historically hosted a relatively rich (seven taxa) stygofauna community (unpublished data from Chevron), but in 2020 no stygofauna were identified using the haul net. With eDNA, however, we detected five stygofaunal taxa, including some previously detected, such as Decapoda. In 2021, after a year of heavy rainfall, Amphipoda, Cyclopoida and Ostracoda specimens were again collected at this site using the haul-net samples, but were not detected using eDNA metabarcoding.

| Effect of substrate on eDNA metabarcodingbased detection of stygofauna
In 2020, we detected higher stygofaunal ZOTU richness in the shallow water and sediment samples than the deeper water samples, although not significantly so (Figure 3). In 2021, our results showed a strong decrease in stygofauna ZOTU richness in shallow water samples (p = .016), a weak increase in the sediment (p = .12) and no change in the deeper water samples (p = .77; Figure S2). As a result, in 2021, the stygofauna ZOTU richness in sediment samples was significantly higher than the other two substrates (p < .05; Figure 3).
In terms of compositional differences between the sampled substrates, 100% of stygofauna taxa identified could be detected from the shallow water samples in 2020 and 82% from sediment samples

| Effect of sampling order on eDNA metabarcoding-based detection of stygofauna
Water samples collected before haul nets were visibly observed to contain whole specimens of stygofauna on two occasions; these individuals were removed from the sample with sterilized forceps before filtering. For some of the samples collected before haul-nets there were one or two stygofauna ZOTUs with very high read abundances ( Figure S3), indicating some specimens may have been caught on the filter paper. However, from presence/absence data the eDNA samples collected prior to haul-net sampling were no richer in stygofauna ZOTUs than the samples collected after haul-nets (p = .115) although "bore" was a significant factor (p = .001; Figure S4). Similarly, the order of sampling was not significant in explaining variation in stygofauna community composition (p = .131), but communities varied significantly between bores (p = .001; Figure S5). Multipattern analysis indicated ZOTUs such as those identified as Stygoridgewayia trispinosa could be associated with eDNA samples taken before or after haul-nets depending on the bore (Table S3). The number of specimens per bore collected through haul-net samples was similar between the 2 years of sampling, ranging from one to three at bore SFRD3-03 to 362-45 at bore X62M (see Table S4).

| Physicochemical properties
Dissolved oxygen and salinity showed an inverse relationship with increasing water depth, with oxygenated freshwater found at the surface and oxygen-depleted salt water in deeper water ( Figure S5). Temperature remained relatively constant at ~30°C, while pH remained neutral (~7) for most bores though some became slightly more acidic (pH = 6.5) with depth (M52, L4A; Figure 1).

F I G U R E 3
Stygofauna richness and composition differences between the eDNA sample types: deeper water samples, regular shallow water samples and sediment samples. (a) Stygofauna ZOTU richness is shown above and below (b) is a Venn diagram depicting the number of stygofaunal taxa detected in each sample type.

| Community composition
The asymptotic regression models showed that 90% of the stygofaunal taxa in each borehole could be detected with five 1-L water replicates (Figure 4). SIMPROF analysis indicated that the communities detected in the different years tended to cluster together based on the bore ( Figure 5). However, there were some exceptions such as MW34 and S62M MW03, where the communities from the different years were found in separate clusters ( Figure 5). Boreholes that were geographically close did not necessarily have more similar stygofauna community composition than those that were geographically distant in either 2020 (r = −.008, p = .514) or 2021 (r = −.116, p = .681). There were more stygofaunal taxa detected in 2021 than in 2020, using eDNA metabarcoding (24 vs. 13, respectively), but no difference was observed with the specimens that were morphologically identified (13 vs. 13; Table 1). The most common stygofaunal taxa included Cyclopoida (14 bores), Catenulida (13 bores), Calanoida (12 bores) and Decapoda (10 bores).

| DISCUSS ION
Here we present the first metabarcoding study which examines the factors affecting community-wide eDNA metabarcoding-based biodiversity assessments in groundwater. As such, it is intended as the first step in developing a robust protocol for effective eDNAbased stygofaunal detection and monitoring. Environmental impact assessment and monitoring for subterranean fauna is currently conducted using traditional approaches such as haul-net sampling, but, as with most sampling protocols, these methods have inherent biases (Saccò, Campbell, et al., 2022;Saccò, Guzik, et al., 2022). We show that broad eukaryote 18S metabarcoding assays were more successful than haul-net sampling in detecting the breadth of major stygofaunal taxon groups known from groundwater on Barrow Island, but that key taxa were notably absent from metabarcoding results compared to those collected using haul nets and identified morphologically. eDNA was also able to detect stygofauna at locations that haul nets did not. In testing the protocols for sampling groundwater fauna for eDNA metabarcoding, we show the shallow F I G U R E 5 Stygofauna community composition from bore samples. Scatterpie plots show the stygofaunal orders detected at each borehole using eDNA metabarcoding. Dendogram depicts SIMPROF cluster analysis of bore stygofauna community composition. Different colours indicate significant clusters (p < .01).
water and sediment samples were the most effective substrate for detecting diversity in the majority of taxa. However, temporal variation in detection of taxonomic diversity in all samples was observed between 2020 and 2021.

| Comparison of eDNA metabarcoding vs. morphological identification detection of stygofauna
eDNA metabarcoding consistently detected more taxa than were identified morphologically from haul-net samples (Figure 2). Some of the taxa detected using eDNA metabarcoding but missed by haul-nets included soft-bodied taxa (catenulid flatworms, Annelida) that may be more difficult to catch in haul-nets or identify because of their small size, but leave sufficient DNA for eDNA detection ( Figure 2). These soft-bodied taxa have been identified as firstand second-level predators common to groundwater habitats. In particular, they feed on protozoan and metazoan prey and bacterial biomass which contribute to groundwater purification and nutrient cycling, important for ecosystem functioning (Cunillera-Montcusí et al., 2022;Herrmann et al., 2020). Platyhelminthes (Catenulida) and copepods (e.g., Cyclopoida) are considered high-level predators in some groundwater systems (Herrmann et al., 2020), although recent studies suggest they may also be feeding on roots (Saccò, Campbell, et al., 2022) and bacteria (Saccò et al., 2021). These organisms are themselves preyed on by stygobiont teleosts, if present (Humphreys & Feinbergz, 1995), and are important for the food web structure and trophic dynamics of groundwater systems. Previous haul-net specimens from Barrow Island have been dominated by crustaceans, with Annelida and Chordata being far less numerous in collections . The differences in abundance between these groups may be an artefact of the sampling method, rather than a reflection of their true abundances in this groundwater ecosystem (Herrmann et al., 2020;Humphreys & Feinbergz, 1995).
While morphological identification detected exclusively the larger organisms, eDNA metabarcoding could detect smaller, lower-level consumers, such as Cercozoa and Ciliophora, as well as fungi and bacteria (Herrmann et al., 2020;Oberprieler et al., 2021). The absence of taxon groups from trophic levels other than those that are free-swimming, in haul-nets, indicates that traditional methods of sampling are frequently incapable of detecting organisms below tertiary consumers, providing a narrow window into the diversity of stygofauna communities inhabiting groundwater systems.
With haul-nets, we were able to collect specimens from all crustacean orders known to occur in groundwater ecosystems at Barrow Island . Most of these taxa were also detected using 18S eDNA metabarcoding. Some crustaceans were reliably detected with both methods (Calanoida, Cyclopoida, Bathynellacea and Decapoda) while Ostacoda and Amphipoda were detected more often from haul-net samples than by eDNA metabarcoding. We know that the 18S eukaryote assay can amplify these taxa (based on in silico studies; our unpublished data), but primer bias against crustaceans appears to have resulted in false negatives as the DNA from these taxa are "drowned out" by the many other taxa this broad assay can detect (Korbel et al., 2017;Pochon et al., 2013).
Crustaceans tend to shed less DNA in their environment than softbodied organisms (Andruszkiewicz Allan et al., 2021), which may further bias results against these taxa.
Thermosbaenacea and Isopoda are important components of the subterranean ecosystem at Barrow Island that probably represent Tethyan relicts (Humphreys, 1993;Page et al., 2018) However, these taxa were not detected using eDNA metabarcoding even though they are known to occur on the island and were collected using haul-nets. Further investigation showed there were no mismatches in the primer binding sites for these taxa, but the amplicon length for thermosbanaeceans and isopods was found to be  Table 2). The lack of specificity in the morphological identifications probably occurred because of the cryptic differences between species, and poorly resolved taxonomy Humphreys et al., 2013). Like Korbel et al. (2017), the inconsistent detection of key stygofaunal taxa using eDNA metabarcoding indicates that an 18S assay alone is not yet good enough as a "stand-alone" method for stygofaunal surveys, even with the use of a custom reference database. While the development of more specific assays may increase the specificity and detection probability of eDNA metabarcodingbased surveys, currently, using a combination of both morphological characteristics and eDNA-based survey methods yields a greater diversity than haul-net sampling and morphological identification alone (as is currently standard protocol). Additionally, eDNAbased methods are notoriously unreliable for abundance estimates (Elbrecht & Leese, 2015), while morphological identification provides some population data on relative abundances of taxa from the specimens collected (Eberhard et al., 2009;Halse et al., 2014;Humphreys et al., 2013). Morphological and eDNA-based assessments of stygofaunal diversity are thus complementary, and combining these survey approaches provides a more complete assessment of the communities present in the aquifer under investigation.

| Troglofauna
Like West et al. (2020) we detected troglofauna (e.g., Schizomida, Pseudoscorpiones) from eDNA metabarcoding of groundwater, and troglofauna were also incidentally collected in haul-nets as they were raised up from the borehole access point from the groundwater. These organisms typically live in geology surrounding the groundwater habitat and can fall in and can be eaten by stygofaunal predators; in this way troglofaunal detections from groundwater are analogous to using water samples from estuaries (Polanco et al., 2021) and ponds (Harper et al., 2019) to detect terrestrial diversity. While bores can act as large pitfall traps, depending on the location and size of the slots in the pipe, such results are not consistent between bores and frequently unknown. While there are methods to assess troglofaunal communities (e.g., litter traps, Javidkar et al., 2018), the application of eDNA-based methods in this area is untested. Scrapes of soil from inside drilled holes may provide an adequate eDNA substrate; however, this is difficult on Barrow Island because the PVC casing of the bores prevents the collection of scrapes.
In addition to troglofauna, other taxa opportunistically detected using eDNA and occasionally haul-nets were surface-dwelling organisms (e.g., reptiles, ants, cockroaches, collembolans), animals that either fell in the boreholes or DNA being leeched through the soil to the groundwater. Other subterranean eDNA studies have also detected terrestrial taxa from groundwater, whether they sampled boreholes (Oberprieler et al., 2021) or caves (West et al., 2020). We avoided contact between sampling equipment and the ground to minimize contamination, but any access to groundwater also risks exposing the groundwater to animals/DNA from the surface. This is why eDNA detection of an organism in groundwater is not necessarily proof of occupancy. However, this potential problem can be at least partly ameliorated by expert knowledge of the taxa concerned and their biology.

| Effect of substrate and sampling order on the diversity of stygofauna detected using eDNA
Methodological protocols for sampling eDNA metabarcoding are known to be critical considerations if these methods are to be used in future environmental impact assessment and monitoring (Díaz-Ferguson & Moyer, 2014;Dickie et al., 2018). In applying this knowledge to groundwater, we investigated the effectiveness of sampling various substrate types and the order of sampling to examine the impact of detecting community-wide diversity and composition using eDNA metabarcoding. For substrate type, we found that the treatments tested identified different components of the community (Figure 3). In shallow water, eDNA metabarcoding detected 100% of the total stygofauna detected in 2020, and 55% in 2021. In contrast, the deeper water samples detected only one additional stygofaunal taxon, and was also the substrate with the lowest total richness (Figure 3). This result agrees with previous traditional monitoring surveys using sampling techniques that were able to sample at various depths (Bork et al., 2008;Hahn, 2005).
The drop in diversity in deeper water has generally been attributed to the decrease in oxygen and food supply with increasing depth ( Figure S5; Brunke & Gonser, 1999;Datry et al., 2005).
Inverse to oxygen, salinity on Barrow Island increases with depth, and variation in salinity is known to influence stygofaunal community composition in coastal aquifers (Shapouri et al., 2016) and island cave systems (West et al., 2020). Salinity is one of the main drivers of adaptation, speciation and community assemblage in aquatic systems (Cunillera-Montcusí et al., 2022). Salinization can also alter bottom-up affects with the mobilization of base cations, changes in pH and the mobilization of nutrients (Cunillera-Montcusí et al., 2022). Several stygofaunal taxa (e.g., Harpacticoida, Decapoda, Ostracoda) have been previously associated with freshwater or brackish environments (West et al., 2020) while others (Thermosbaenacea) are primarily marine and appear to tolerate variable levels of salinity (Jaume, 2008).
The taxa we detected primarily in deeper water samples (more saline, less oxygen) tended to be worms (Polychaeta) that can survive a range of salinities (Glasby, 1999). The stygofauna richness in deeper water samples was consistently low, as was the proportion of unique taxa in this substrate. In conclusion, if we were to exclude the deeper water samples, we would still be able to detect 95%-100% of the total taxonomic diversity detected using eDNA metabarcoding.
Mineral particles form electrostatic bonds that can protect the DNA from damage (Hou et al., 2014) resulting in concentrations of DNA in sediments being up to 1800 times higher than in water (Turner et al., 2015). Consequently, the sediment samples were expected to contain the greatest diversity of stygofauna; however, our results were variable. In 2021, sediment samples were the richest, and detected the most stygofauna taxa, while in 2020 all the taxa detected from sediment samples could also be detected using shallow water samples ( Figure 3). Rainfall events are responsible for both carbon and nutrient infiltration into groundwater systems, which is exploited by microbes and are direct and indirect diet sources of stygofauna (Saccò et al., 2021). Buxton et al. (2018) also found increased detection probability in pond water rather than sediment, and that there were seasonal fluctuations in both substrates. Considering the magnitude of change observed in the sediment samples, it is probable that we detected organisms that were alive in the last year, rather than very old DNA. This is reassuring because it means we were detecting primarily recent diversity rather than legacy eDNA, Order of sampling (nets or eDNA first) did not make a significant difference to assessment of stygofauna richness or community composition, although there appeared to be variation between bores ( Figure S4). For example, at one bore (SFRD3-03), S. trispinosa was associated with eDNA samples collected prior to nets, but for another (X62M) they were associated with samples collected after nets. The samples collected before disturbing the water column are typically less turbid and therefore easier to filter with fewer inhibitors than those collected after. Therefore, for improved field sampling methodologies (eDNA and haul-nets), it is possible to conduct the collection of groundwater first for eDNA work and then proceed with the haul-net collection of fauna for manual identification in the laboratory. However, we have not explicitly tested whether sampling for eDNA might impact the haul-net results. Such an undertaking would be extremely difficult because additional haul-net samples are more likely to be affected by the preceding net hauls than the preceding eDNA samples. In the present study, specimens were occasionally observed in the water samples, and larger individuals could be easily removed at the water filtering stage. Further, even if organisms were caught on the filter membrane, they did not overwhelm the results of stygofaunal diversity ( Figure S4). In future, water samples could be prefiltered to ensure any specimens caught in eDNA samples are removed and identified, if required.

| Bore properties
There was not much overlap in stygofauna communities between sampling locations (bores), and bores closer together did not reveal any greater similarity in community composition ( Figure 5).  (Hahn & Matzke, 2005). However, while the abundance of the taxa may be skewed, the diversity of taxa tends to remain the same in the bore/trap as in the aquifer (Bork et al., 2008;Hahn & Matzke, 2005).
The exception to this may represent the ability of different methods to detect naturally rare taxa (Saccò, Guzik, et al., 2022). In fact, the increased abundance of taxa in bores may increase the detectability of stygofauna using eDNA methods, allowing sampling of an enriched community. Therefore, purging bores before sampling is not recommended for stygofaunal surveys specifically. Furthermore, if pumps are to be used, then the first flow-through will always contain the highest number of specimens (Oberprieler et al., 2021;Sorensen et al., 2013), and, consequently, the most DNA. However, purging a bore may be appropriate when the aim is to characterize the microbial communities of the aquifer (Korbel et al., 2017). Sampling the aquifer directly may require a greater volume of water and/or more replication, possibly using pumps with in-line filtration to sample a larger area.
This study took place at a single location, on an island with freshwater overlying saltwater. The results provide a valuable case study; however, they may not be representative of mainland groundwater systems. For example, Mouser et al. (2021) found that the velocity of groundwater increased the detectability of cavefish and cave crayfish using eDNA-based surveys, while the groundwater on Barrow is relatively still. Flowing water may also reduce the difference between inside and outside the bore, and further render purging the bore unnecessary. Similarly, while 5× 1-L eDNA samples may be sufficient to obtain ~90% of the stygofaunal diversity (Figure 4), our model may be specific to this groundwater system, location and metabarcoding assay. Further studies in diverse ecosystems are needed to determine the across-system applicability of these results.

| Temporal variation
There was considerable variability in communities and substrates between the two sample years, as was expected since sampling occurred during a 3-year drought (2020), and the following year when rainfall increased (2021). Groundwater recharge plays an important role in nutrient availability and oxygenation (Datry et al., 2005;Saccò, Blyth, Humphreys, Middleton, et al., 2020), with recharge being related to rainfall (Hendrickx & Walker, 2017).
Under high rainfall, even with low recharge rates, there are increases in carbon inputs that increase the food resources for stygofauna (Saccò et al., 2021;Saccò, Blyth, Humphreys, Middleton, et al., 2020). Not only were there significant shifts in stygofaunal ZOTU richness found in the different substrates (Figure 2), several bores also showed inconsistent results between years. For example, MW34, historically rich in stygofauna, produced no specimens using haul-nets in 2020, while eDNA of previously caught taxa could still be detected. During the drought period between 2017 and 2020 the stygofauna in MW34 may have still been present but at lower abundance or different life stages such that they were not found in haul-net samples, but were still detected using eDNA metabarcoding. This hypothesis is supported by the fact that stygofauna were again collected from MW34 in 2021. The findings of Herrmann et al. (2020) show that bacteria depend on the availability of organic or inorganic electron donors, while protists and Metazoa often have a broad prey spectrum and the hydrochemistry-driven distribution patterns were stronger for bacteria than for the eukaryotic communities. There may also be an element of stochasticity involved, considering not all bores showed the same trends. BGMW04, for example, produced specimens of Harpacticoida in 2020 but not in 2021. As we only had two sampling periods a year apart, we cannot comment on the rate at which stygofauna communities may change and how they may vary across smaller timescales. Seasonality can affect not only stygofauna communities (Saccò, Blyth, Humphreys, Karasiewicz, et al., 2020), but also eDNA degradation and detection probabilities (Troth et al., 2021). Repeated sampling, preferably in different seasons, is necessary to assess stygofaunal diversity using both morphological identification and eDNA metabarcoding.

| CON CLUS ION
This study demonstrated that eDNA metabarcoding has the potential to improve the assessment of stygofaunal communities by widening the breadth of biodiversity that can be detected. It can act as a noninvasive tool that is complementary to morphological assessments and able to detect additional taxa (i.e., food web complexity and trophic levels) that are often missed (e.g., Platyhelminthes, Annelida). The study has also demonstrated the need for multiple substrates and sampling events to account for natural temporal variability in an arid and subtropical continental island off the coast of Western Australia. To standardize subterranean eDNA metabarcoding for environmental monitoring we must address the bias in taxa detected by improving and validating metabarcoding assays for stygofaunal groups so that all major taxa can be detected. Sequencing specimens for primer design will naturally also improve reference databases inclusive of provenance and intraspecies variability, which is particularly important when working with ancient linages of shortrange subterranean endemic species with the future goal of making reliable molecular, species-level identifications possible.
While further optimization is necessary, we have shown the feasibility of this method for assessing stygofaunal communities in these diverse and largely unstudied subterranean ecosystems.
Groundwater is an important global resource that hosts ancient lineages, with many short-range endemics. If we do not have the tools to accurately assess these communities, there is no way of knowing how they are affected by human activity.

ACK N O WLE D G E M ENTS
We acknowledge the traditional owners of the land on which this research was undertaken and pay our respects to Elders past, present and emerging.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available on the dryad digital repository at https://doi.org/10.5061/dryad.dv41n s22g.