Two tissue samples of M. justitia were obtained from ethanol preserved specimens in the Western Australian Museum (WAM P.35191-001, 2016, Barrow Island, bore SF-RD3-02-20160229-02, and WAM P.35192-001, 2016, Barrow Island, bore SF-RD7-01-20160302-01) and two tissue samples of M. veritas (WAM P.35198-001, 2012, Cape Range, Bundera area bore S1, and WAM P.35200-001, 2015, Cape Range, Camerons Cave bore CC2; Suppl. material 1: table SS1) for DNA extraction and mitochondrial genome sequencing. Total genomic DNA (gDNA) was extracted using QIAGEN blood and tissue kit (Qiagen; Venlo, Netherlands) following the manufacturer’s protocol. The concentration of gDNA was measured using a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Genomic DNA libraries were prepared using a DNA Prep library kit (Illumina; New York, USA), a bead-based transposome complex to tagment gDNA, targeting the 1-9ng of input gDNA. The genomic libraries were sequenced on an Illumina MiSeq (Illumina, San Diego, CA, USA) with standard 500-cycle V2 chemistry (250 bp paired-end sequencing run) following the manufacture’s protocol.

Mitochondrial genome assembly was completed using the program MitoZ version 3.6 (Meng et al. 2019) with data quality control, filtering, assembly and annotation performed. Software defaults for assembly and annotation were used with ‘vertebrate’ specified as the taxonomic source to guide the selection of mitochondrial genome scaffolds. Gene annotations were confirmed using Galaxy Australia (https://usegalaxy.org.au) with MITOS2, a publicly accessible web-based platform for computational biology. Structural features of annotated mitochondrial genomes, including gene length, nucleotide composition, and start/stop codons, were extracted in Geneious Prime version 25.0.3 (Kearse et al. 2012).

The mitochondrial control region was selected as the targeted qPCR assay for M. justitia and M. veritas given its high variability and ability to provide greater resolution for distinguishing closely related populations, making it useful for assessment of both intraspecific and interspecific variation. Consensus sequences were aligned and conserved primer regions targeting ~230bp region for M. justitia and M. veritas were identified, with the final assay designed using Primer3 in Geneious Prime. The blast function was used in Geneious Prime to perform in-silico PCR across the NCBI GenBank database to retrieve any sequence that would bind with four or fewer mismatches to the forward and reverse primers. gBlocks of 280bp for M. justitia and M. veritas control region were obtained from Integrated DNA Technologies (IDT) which allow for the primer binding site, a 10bp overhang at the 3’ end of the strand, and an 8bp sequence (AGCTAGCT) insert to allow for detection of synthetic DNA verses target DNA when future processing of samples occurs.

A range of annealing temperatures was assessed to optimise assay conditions. To validate the optimal annealing temperature, a gradient qPCR was performed in seven reactions for each species ranging from 55 °C to 61 °C (run in 2-degree increments) using a theoretical input of approximately 4 × 10¹⁰ copies per ul of synthetic DNA (gBlocks; Integrated DNA Technologies) and serial diluted 10-fold across the seven reactions to ~42,450 copies per ul. Those conditions found to generate the lowest cycle threshold (CT) value, optimal amplification curve, and sharp melt peak were considered to represent the optimal annealing conditions and were therefore implemented throughout the remainder of the study.

A standard curve using the optimal temperature, with analysis performed on ten replicates for each of the 7-points in a 10-fold dilution series of the synthetic M. justitia and M. veritas control region DNA (theoretical ~0.04 copies to 42,450 copies; Suppl. material 1: tables S2, S3). The definitions of Klymus et al. (2019) for LOD and LOQ were implemented, where LOD represents the lowest standard concentration of template DNA that produced at least 95% positive replicates and LOQ represents the lowest standard concentration that could be quantified with a coefficient of variation (CV) value below 35%. To determine these values, data from the replicate standard curves in RStudio using the script published by Klymus et al. (2019) with the parameters left as default were used.

For exclusion testing, eDNA extracts were obtained from a metabarcoding study (A. Lawrie pers. comm.) of surface water from the Swan Coastal Plain of Western Australian from which several fish taxa were identified, including Geophagus (cichlids), Pseudogobius (gobies), Cypriniformes (minnows and carps), Cyprinodontiformes (toothcarps), Gobiiformes (gobies), Carassius auratus (Goldfish), Gambusia holbrooki (Eastern Gambusia), Geophagus brasiliensis (Pearl Cichlid) and Leptatherina wallacei (Western Hardyhead) commonly found in surface water in a freshwater system in southern Western Australia (Suppl. material 1: table S4). Criteria for successful exclusion tests were for no amplification to occur ≤ 40 cycles in any of the 7-point (10-fold dilution series) of the fish eDNA from surface water with no expected amplification. Further replicate series of 7-point (10-fold dilution series) of the fish eDNA from surface water was spiked with two concentrations of synthetic DNA, 424 and 42 copies/ul (M. veritas) and 42 and 4 copies/ul (M. justitia) as positive controls with expected amplification. A total of 40 reactions including non-template controls for exclusion testing were conducted. Positive amplified products (£40 cycles) were checked for expected size using QIAxcel (Qiagen; Venlo, Netherlands) and sent for sanger sequencing at AGRF to confirm the presence of the target species.

For inhibition testing, an eDNA extract from a previous metabarcoding study (A. Guthrie pers. comm.) that showed PCR-inhibitory properties where the neat (i.e., no dilution) showed no amplification and the 1/10 and 1/100 dilution failed to show a 3.3 cycle shift when tested with a species-specific probe assay. Two concentrations of synthetic DNA, 424 and 42 copies/ul (M. veritas) and 42 and 4 copies/ul (M. justitia) were spiked in a 7-point replicate series where the inhibition was diluted in a 10-fold series. A total of 32 reactions including non-template controls for inhibition testing were conducted. Positive amplified products were checked for expected size using QIAxcel.

The targeted qPCR assay was in-situ validated on groundwater samples collected from the Cape Range peninsula in Western Australia where M. veritas is known to occur. A total of six groundwater eDNA samples collected during a Blind Cave Eel (O. candidum) study and Cave Gudgeons (M. veritas) were visually observed to co-occur at the time of groundwater collection. These samples were tested in triplicate reactions at two dilutions to account for potential inhibition: neat (i.e., no dilution; n = 18 reactions) and 1/10 (n = 18 reactions), alongside positive and negative controls for a total of 38 reactions. Positive amplified products (≤40 cycles) were checked for expected size using QIAxcel and sent for sanger sequencing at Australian Genome Research Facility, Perth, Western Australia to confirm the presence of the target species.

Sequencing data from two M. justitia (AWFS-F21-191 with 8,190,706 reads and AWFS-F21-192 with 7,368,386 reads) tissue samples were assembled to complete, circular sequences with high read coverage (>30×) across each mitogenome. Sequencing data from two M. veritas (AWFS-F21-189 with 1,197,618 reads and AWFS-F21-190 with 634,158 reads) resulted in a 332 bp fragment of the ND5 and ND6 gene not being assembled with the data retrieved (AWFS-F21-189; M. veritas), although a near circular sequence with > 15× coverage was achieved for the remainder of the mitogenome. Sample AWFS-F21-190 with less than 635,000 reads failed to produce a mitogenome. The final assembly length being 16,588bp, consistent with expected mitogenome sizes, and each contained the full complement of mitochondrial genes for vertebrates. As found in other gudgeons (Eleotris acanthopomus; MH367493) we identified two rRNA, 22 tRNA, and 13 protein-coding genes, plus the putative control region and these were annotated and deposited on GenBank (PX694302, PX694303, PX694304). High-throughput sequencing achieved robust coverage for three out of four tissue samples, ensuring the completeness and reliability of the resulting contigs. Each assembled mitogenome is therefore considered essentially complete and well-validated, with only minor caveats in one problematic region (ND5 and ND6; M. veritas) that are inherent to mitochondrial genome assembly (Ribarska et al. 2022). The closest mitochondrial match (78.6%) currently on GenBank was to Oxyeleotris marmorata (KP663727), the marble gudgeon found in Asia. The aligned M. justitia mitogenomes from tissue collected in 2016 from two bores on Barrow Island, Western Australia, contained 12 single nucleotide polymorphisms (SNPs) corresponding to 0.1% sequence variation. An alignment of the M. justitia mitogenomes with the M. veritas genome, excluding the CR and 322 bp region of the ND5 and ND6 gene showed 2,538 SNPs (11.3% sequence variation) amongst the three mitogenomes. These well-characterised genome assemblies offered the necessary resolution for downstream analyses, including in-silico assay design. The accuracy of the builds was further validated through reference-based mapping and consensus sequence confirmation, establishing a solid foundation for future genetic and conservation monitoring studies.

A successfully designed genus-specific assay, herein called MilyeringaCR, for detection of M. veritas and M. justitia and which amplifies a 236 bp fragment of the control region, excluding primer binding sites has been validated. The 263 bp fragment of the control region exhibited 85 SNPs corresponding to 29.9% sequence variation between the two species within the targeted region. The optimised MilyeringaCR assays were in 25 mL reaction volumes consisting of 2 mM MgCl₂ (Fisher Biotec, Australia), 1 × Taq polymerase buffer (Fisher Biotec, Australia), 0.4 mM dNTPs (Astral Scientific, Australia), 0.1 mg bovine serum albumin (Fisher Biotec, Australia), 0.8 mM of each primer (MilyeringaCR654F and Milyeringa903R), 0.2 mL of AmpliTaq Gold (AmpliTaq Gold, ABI, USA), 0.6 µL SYBR Green and 2 mL of template eDNA or synthetic DNA. The cycling conditions were: initial denaturation at 95 °C for 5 min, followed by 50 cycles of 95 °C for 30 s, 57 °C for 30 s, 72 °C for 45 s, a final extension at 72 °C for 10 min followed by a melt curve from 60 °C to 95 °C with data capture every 1 °C increment. The primer sequence for the MilyeringaCR assay are as follows: MilyeringaCR654F 5’-ACATTTCAGAGTGCAGCGC-3’ and MilyeringaCR903R 5’-GCTTTCCTGTT TCCGGGGG-3’.

During in-vitro testing, the inclusivity criteria (amplification of target DNA at ≤ 40 cycles) and exclusivity criteria (amplification of non-target DNA at ≥ 40 cycles) were met successfully. The assay was tested against numerous fish taxa found in surface water in southern Western Australia (see methods) with no amplification observed using the MilyeringaCR assay.

Analysis of the dilution standard curve for M. justitia synthetic DNA showed good precision (R² = 0.998), a slope of -3.46, and efficiency of 94%. Analysis of the standard curve for M. veritas synthetic DNA showed lower precision (R² = 0.862), a slope of -3.66, and efficiency of 88%. Based on the criteria applied, the discrete LOD for M. justitia was calculated to be 1.36 copies/reaction (where 95% of the replicates showed amplification) and the LOQ was 6 copies/reaction (as determined by decay modelling using a CV of 0.35). The discrete LOD for M. veritas was calculated to be 12.8 copies/reaction and the LOQ was 524 copies/reaction. Using the modelled calculation of Klymus et al. (2019), the calculated LOD for M. justitia was 1.36 copies/reaction and M. veritas to be 12.8 copies/reaction (Fig. 1).

Quantitative PCR performance varied between the two species with the genus-specific eDNA assay. The M. justitia assay demonstrated high sensitivity, with a modelled limit of detection (LOD) of 1.36 copies per reaction and a modelled limit of quantification (LOQ) of 6 copies per reaction. Empirical testing further supported this sensitivity, with consistent amplification observed at 4 copies per reaction. In contrast, the M. veritas assay showed a higher LOD of 12.8 copies per reaction and a modelled LOQ of 524 copies per reaction. However, empirical results indicated reliable amplification at much lower concentrations, with all 10 replicates amplifying successfully at 41 copies per reaction (Table 1). The MilyeringaCR assay targeting M. justitia showed high sensitivity and precision, while the M. veritas assay demonstrated reliable amplification below its modelled LOQ, indicating effective detection at lower concentrations (Table 1).

The MilyeringaCR assay, when tested in the presence of inhibition, failed to amplify M. justitia synthetic DNA, at neat and 1/10 dilution concentrations when four and 42 copies/reaction were spiked in. The remaining five dilution series reactions showed CT values ranging from 30.9 to 34. Similarly for M. veritas, failed amplification was observed at neat and 1/10 dilution concentrations, as expected, when 41 and 411 copies/reaction were spiked in and the remaining five reactions showed CT values ranging from 32.1 to 37.8. The MilyeringaCR assay was able to detect M. justitia and M. veritas DNA in the presence of environmental PCR inhibitors with CT values < 40 cycles.

Testing of field samples demonstrated that the assay could detect M. veritas DNA in groundwater samples from Cape Range, with successful PCR amplification observed in 15 of the 18 reactions (83%) of the neat eDNA groundwater extracts with melt peaks identical to the positive control. Sanger sequencing confirmed the amplicons to be 100% match to the reference mitogenome indicating successful in-situ validation as the MilyeringaCR was able to amplify its intended target DNA. The 1/10 dilution reactions showed 11 of the 18 (61%) positive reactions had CT values > 40 cycles and the criteria for successful inclusion tests were for amplification to occur ≤ 40 cycles. The amplicon size for the 1/10 dilution reactions was not confirmed on QIAxcel, although the melt peaks were identical to the positive controls. No amplification was observed in laboratory or PCR negative controls, while all positive controls amplified successfully.

Conservation monitoring cave gudgeons (Milyeringa spp.) presents several significant challenges due to their biology, habitat, and conservation status. These small, eyeless fish live in lightless subterranean aquifers and caves, with narrow ranges, deep and shallow underground and completely in the dark, but sometimes venturing to the surface (Humphreys 2001; Chakrabarty 2010; Page et al. 2018). However, most Milyeringa habitats are physically inaccessible to humans, limiting the use of visual observation or direct sampling. Compounding this, Milyeringa spp. occur at extremely low densities, making them difficult to detect using traditional methods like visual surveys (cave diving and cameras), netting or baited trapping. Both M. veritas and M. justitia are species of conservation concern, limiting the ability to capture or handle individuals for study (DCCEEW 2025). Monitoring must be non-invasive and ethically sound, which restricts the use of some conventional techniques. In the last decade, eDNA technology has revolutionised the survey of rare aquatic fauna, including subterranean fauna, and provides a viable alternative to conventional monitoring practices (Gorički et al. 2017; White et al. 2020; van der Heyde et al. 2023).

The genus-specific targeted qPCR assay, presented here, will dramatically improve detection sensitivity of Milyeringa spp. in their subterranean aquifer habitats. The assay demonstrates very high specificity and sensitivity for M. justitia with a modelled LOD of 1.36 copies/reaction and LOQ of 6 copies/reaction which should significantly improve conservation monitoring of this endangered species endemic to the groundwater on Barrow Island, Western Australia. Despite Barrow Island being a Class A Nature Reserve since 1910, is the site of the large Gorgon gas processing project and, until recently, an active oil field since the 1960s. The known habitat of M. justitia coincides with the island’s industrial zone (Saccò et al. 2022c). This situation presents unique conservation challenges: the primary threat being habitat disturbance from industrial activities (Moore 2019b). Therefore, this genus-specific targeted qPCR assay will fundamentally improve the ability to detect and routinely monitor this cryptic subterranean fish.

The MilyeringaCR assay developed for M. veritas demonstrates moderate specificity and sensitivity, with a modelled LOD of 12.8 copies per reaction and a LOQ of 524 copies per reaction. Indeed, the sensitivity of the assay was substantially lower for M. veritas (LOD: 12.8 copies/reaction; LOQ: 524 copies/reaction) than for M. justitia. This discrepancy is unlikely to be due to differences in primer binding efficiency, as in-silico validation confirmed identical primer binding sites between the two species. Furthermore, both assays were validated using synthetic DNA standards, ensuring precise quantification of input DNA molecules during PCR reactions. These controls increase confidence that the observed differences in sensitivity reflect true assay performance rather than experimental variability.

While the modelled LOQ for the M. veritas assay was estimated at 524 copies per reaction using the Klymus et al. (2019) R script with a coefficient of variation threshold of 35%, empirical testing demonstrated consistent amplification at much lower concentrations. Specifically, all 10 replicates amplified successfully at 41 copies per reaction, inclusive of the inhibition testing and the exclusion testing, indicating reliable detection below the modelled LOQ. This discrepancy suggests that while the assay may not meet the strict precision criteria for quantification at lower concentrations, it remains effective for detection. Reporting both the modelled LOQ and the empirically supported detection threshold provides a more comprehensive view of assay performance, particularly for field applications where consistent detection is often more critical than precise quantification. This highlights the importance of interpreting modelled thresholds alongside empirical performance, particularly when assessing an assay utility for field-based eDNA applications.

Overall, this genus-specific eDNA assay offers a robust and non-invasive tool for detecting Milyeringa spp. in subterranean environments, overcoming the limitations of traditional survey methods. Its demonstrated sensitivity, particularly for M. justitia, supports its application in routine conservation monitoring and conservation management of these cryptic and threatened cave fishes.

The mitochondrial genome data generated in this study provide preliminary insights into the genetic diversity of M. justitia on Barrow Island. Sequencing of the mitogenomes (n = 2) revealed intraspecific variation, with 12 SNPs detected, corresponding to 0.1% sequence divergence. In contrast, comparative analysis with the mostly complete mitogenome of the closely related M. veritas identified 2,538 SNPs, representing 11.3% sequence divergence. These findings could suggest limited haplotype diversity within the M. justitia population, potentially indicative of a population bottleneck, founder effect or restricted gene flow. However, given the small sample size, this apparent homogeneity should be interpreted cautiously. Broader spatial and temporal groundwater sampling on Barrow Island may uncover additional mitochondrial diversity. Such efforts would enhance our understanding of population structure and evolutionary history for this species of conservation concern.

Comparison of the 236 bp fragment of the control region from the M. veritas mitogenome with sequences derived from groundwater samples collected at two distinct sites in Cape Range revealed a 100% match. This perfect sequence identity supports the utility of eDNA with a targeted qPCR assay for detecting this species in groundwater systems. The control region, known for its high mutation rate relative to other mitochondrial loci, is typically expected to exhibit considerable haplotypic variation. The absence of such variation in the current dataset may reflect limited sampling or a genetically homogeneous population. However, given the small fragment size and the early stage of sampling, it is possible that additional haplotypes were not captured in this study. Future targeted qPCR assay efforts incorporating broader spatial and temporal sampling may reveal greater genetic diversity within M. veritas populations in Cape Range. These findings underscore the importance of continued eDNA surveillance to refine our understanding of population structure and connectivity in subterranean aquatic ecosystems.

This study demonstrates the critical role of eDNA in overcoming the inherent challenges of monitoring cryptic subterranean species such as Milyeringa spp. The development and validation of this genus specific assay provides a sensitive, non-invasive tool for detecting these endangered cave fishes in groundwater habitats. The higher sensitivity of the assay to M. justitia compared to M. veritas underscores the importance of species-specific optimisation for reliable detection. Preliminary mitochondrial genome analyses revealed limited haplotype diversity within M. justitia and sequence identity between M. veritas and Cape Range groundwater samples, suggesting potential genetic homogeneity; however broader sampling is required. Future targeted qPCR assay surveys across spatial and temporal scales will be essential to refine our understanding of population structure, connectivity, and evolutionary history. By integrating molecular tools with conservation management, this work establishes a foundation for routine conservation monitoring in subterranean ecosystems, offering a scalable approach to safeguard species of high conservation concern in environments where traditional survey methods are impractical.