Delimiting species is a challenge, especially in scenarios of diversification with gene flow and when species are now allopatric where reproductive isolation cannot be directly tested. Continental burrowing crayfishes of the genus Parastacus present a disjoint distribution in southern South America. One of the species is P. nicoleti, which lives in underground waters in swampy and wooded areas of southern Chile. A previous assessment based on mitochondrial DNA sequences suggest that the taxon may represent a species complex. Here, using thousands of nuclear genomic single-nucleotide polymorphisms obtained via RADSeq from 81 specimens collected at 27 localities throughout the distributional range of the species, we apply an integrative species delimitation approach to test species boundaries and to investigate some aspects of the speciation process. Our analyses corroborate previous results; a scenario that we favor suggests that the P. nicoleti encompasses seven distinct species. Additionally, demographic analyses show that the distinct species have followed distinct trajectories in size change during the last 17.5 million years and that speciation in this group occurred both in strict isolation as well as in the presence of gene flow.