The crustacean marine isopod species Haploniscus bicuspis (G.O. Sars, 1877) shows circum-Icelandic distribution in a wide range of environmental conditions and along well-known geographic barriers, such as the Greenland-Iceland-Faroe (GIF) Ridge. We wanted to explore population genetics, phylogeography and cryptic speciation as well as to investigate whether previously described, but unaccepted subspecies have any merit. Using the same set of specimens, we combined mitochondrial COI sequences, thousands of nuclear loci (ddRAD), and proteomic profiles, plus selected morphological characters using Confocal Laser Scanning Microscopy (CLSM). Five divergent genetic lineages were identified by COI and ddRAD, two south and three north of the GIF Ridge. Assignment of populations to the three northern lineages varied and detailed analyses revealed hybridization and gene flow between them, suggesting a single northern species with a complex phylogeographic history. No apparent hybridization was observed among lineages south of the Ridge, inferring the existence of two more species. Differences in proteomic profiles between the three putative species were minimal, implying an ongoing or recent speciation process. Population differentiation was high, even among closely associated populations, and higher in mitochondrial COI than nuclear ddRAD loci. Gene flow is apparently male-biased, leading to hybrid zones and instances of complete exchange of the local nuclear genome through immigrating males. This study did not confirm the existence of subspecies defined by male characters, which probably characterize different male developmental stages.

All H. bicuspis specimens examined were sampled aboard the research vessels RV Meteor (M85/3), RV Poseidon (POS456) and RV Maria S. Merian (MSM75) during the IceAGE (2011), IceAGE2 (2013) and IceAGE_RR (2018) expeditions. The specimens were sorted at the laboratories of the DZMB (German Centre for Marine Biodiversity Research) and deposited in the collections of the Center of Natural History (Hamburg, Germany; see Supplementary Table S1). Specimens were collected using an epibenthic sledge or VanVeen grab. DNA extractions were performed using the mid-sections of the animals, leaving the cephalothorax and the pleon intact for subsequent morphological analyses and as partial vouchers in the collections. The animals were dissected carefully using a micro scissor and the gut was removed to avoid contaminations. If possible, the same animal was used for proteomics, COI barcoding and ddRAD, transferring half of the dissected tissue in a vial for genetics and the other half into a vial for proteomics. DNA was extracted using the Marine Animal Tissue Genomic DNA Extraction Kit (Neo Biotech) or the Genomic DNA from tissue kit (Macherey-Nagel). DNA was eluted in 70 µL elution buffer. Chelex (Chelating Ion Exchange Resin, Bio-Rad) was used for some of the animals, for which only the COI fragment was sequenced (for protocol see Jennings, Golovan, & Brix, 2020). Based on DNA yield, a subset of samples was selected for ddRAD sequencing. Usually selected samples had >150 ng DNA, to include certain relevant populations a few samples with > 30 ng DNA were selected as well. DNA concentration was measured with a Qubit 3.0 (Invitrogen). If samples had less than 150 ng, they were concentrated via drying. All samples were brought to 24 µl. The protocol described by Peterson et al. (2012) was mostly adhered to, with a few modifications (see also Schwentner & Lörz, 2021; Franchini, Monné Parera, Kautt, & Meyer, 2017).

The files MetaDataxlsx.xlsx and Hellinger_Matrix.csv refer to the proteomic data. All other data belong to the ddRAD.