Most subterranean animals are assumed to have evolved from surface ancestors following colonisation of a cave system, however very few studies have raised the possibility of ‘subterranean speciation’ in underground habitats (i.e. obligate cave-dwelling organisms (troglobionts) descended from troglobiotic ancestors). Numerous endemic subterranean diving beetle species from spatially-discrete calcrete aquifers in Western Australia (stygobionts) have evolved independently from surface ancestors; however, several cases of sympatric sister species raises the possibility of subterranean speciation. We tested this hypothesis using vision (phototransduction) genes that are evolving under neutral processes in subterranean species and purifying selection in surface species. Using sequence data from 32 subterranean and five surface species in the genus Paroster (Dytiscidae), we identified deleterious mutations in: long wavelength opsin (lwop), arrestin 1 (arr1), and arrestin 2 (arr2) shared by a sympatric sister-species triplet, arr1 shared by a sympatric sister-species pair, and lwop and arr2 shared among closely related species in adjacent calcrete aquifers. In all cases, a common ancestor possessed the function-altering mutations, implying they were already adapted to aphotic environments. Our study represents one of the first confirmed cases of subterranean speciation in cave insects. The assessment of genes undergoing pseudogenisation provides a novel way of testing modes of speciation and the history of diversification in blind cave animals.

Sequence data of three phototransduction genes (arr1, arr2, and lwop) were obtained from targeted sequencing of exons from nine phototransduction genes derived from parallel studies of 10 stygobiotic and two surface Paroster species (Tierney et al. 2015, 2018; Langille 2020). The study used RNA baits, designed from transcriptome data from the surface species Paroster nigroadumbratus (Tierney et al. 2015), to enrich for arr1 (1128 bp), arr2 (1194 bp) and lwop (1235 bp) exons, targeting the entire genes. The enriched fragment libraries (followed standard library preparation from Meyer and Kircher 2010, and MYbaits targeted protocol following the MYbaits user manual v2 and Langille 2020) of each exon were sequenced using an Illumina MiSeq and 150 bp paired-end reads. Sequence data from lwop were also obtained from a Geneious blast search of whole-genome shotgun libraries, sequenced using an Illumina MiSeq, for a sympatric sister species triplet from the Sturt Meadows calcrete (Hyde et al. 2018):  P. macrosturtensis (14.6 million sequences), P. mesosturtensis (1.0 million sequences), and P. microsturtensis (4.6 million sequences). The dytiscid sequence data obtained via targeted exon capture were assembled, using a combination of bowtie2 (Langmead and Salberg 2012) and MIRA4 (Chevreux et al. 1999), and aligned to the reference sequence using Geneious v.10.2.6 (Kearse et al. 2012), and the plugin ClustalW (Larkin et al. 2007).

Comparative sequence analyses of the lwop, arr1 and arr2 exon data revealed shared deleterious mutations (frameshift and/or stop codons) among taxa (see results) and we further targeted the exons containing these mutations from additional Paroster species using PCR-amplification and Sanger sequencing analyses (see Table S2 for details of primers; Supplementary Information 1 for detailed PCR-amplification protocol; Table S3 for list of species with successfully/not successfully sequenced data). Sequencing was performed using a Prism BigDye Terminator Cycle sequencing kit (PE Applied Biosystems) with 10 mL reaction volumes according to the manufacturer’s protocol, with reaction products purified using a Multiscreen 384 vacuum well SEQ plate (Millipore Sigma). Reactions were capillary sequenced by the Australian Genome Research Facility (Adelaide, Australia).

Missing bases are represented by 'N'.