The Southern Ocean benthic sea slug, Doris ‘kerguelenensis’, was long thought to be a single, widespread species. This slug is a direct-developing, simultaneous hermaphrodite that can only crawl along the seafloor as an adult. More recently however, a combination of metabolomics, nuclear and mitochondrial DNA sequencing, has revealed a multitude of divergent lineages within that ‘species’. Through Sanger sequencing of Cytochrome Oxidase I (COI), there is now at least 59 putative species. This radiation is thought to be caused by a unique combination of selection and allopatry facilitated through millions of years of episodic glacial cycles. To gain phylogenetic resolution among these species, generating and interpreting large amounts of sequence data from across the genome was necessary. We employed an exon capture approach, targeting thousands of exons from hundreds of genes from across the 59 putative species. The Maximum Likelihood and ASTRAL tree topologies were congruent with our previous species hypotheses; however, we did uncover some mitonuclear discordance. We then further tested species hypotheses using exon capture data and investigated correspondence with secondary metabolites. By increasing the amount of data, we were able to strengthen the resolution of the phylogeny, and examine the distribution of metabolomic traits across the topology.