Cryptic species can be phylogenetically old despite strong sex-biased dispersal

The impact of strongly differentiated populations on species delimitation due to limited or sex-biased dispersal remains challenging and under-explored in the framework of integrative taxonomy. The Mediterranean chafer beetle genus Pachypus is remarkable for its extreme female philopatry, with entirely wingless and subterranean females. This makes Pachypus an interesting case study. Based on a dataset of over 900 protein-coding genes (metazoan universal single-copy orthologs; mzl-USCOs), we investigated phylogeny, species delimitation, gene flow, and population differentiation to provide an integrative assessment of species boundaries. Integrative consideration of all results led to the recognition of 14 mostly morphologically cryptic species, including several new taxa. Most inferred speciation events occurred in the time between the end of the Messinian salinity crisis (about 5.3 million years ago) and the early Pleistocene. Phylogenetically old species and lack of recent speciation was unexpected because of the extreme philopatry, the morphological similarity of the species, and the high degree of differentiation observed among populations of the same species. Speciation was partly associated with the disruption of previously more connected ranges after the Messinian salinity crisis (MSC). This also helps clarify the extent to which the Mediterranean dried out during the MSC, since land connections in the circum-Tyrrhenian region must have persisted long enough for flightless Pachypus females to disperse across drifting land areas connecting the Apennine Peninsula and Africa. We found evidence for historical gene flow between species, while more recent gene flow between populations is low, which is potentially the cause of considerable over-splitting found in the Bayesian Phylogenetics & Phylogeography (BPP) species delimitation analysis. We showed that integrating the outcome of the BPP species delimitation with genealogical divergence index (gdi) values proved to be helpful in some cases but was inconclusive in many others. Generalized Mixed Yule Coalescent (GMYC) and Poisson Tree Processes (PTP) analyses were less prone to over-splitting. This illustrates how species delimitation analyses of cases with restricted or sex-biased dispersal and highly differentiated populations can serve as empirical tests of the utility and robustness of delimitation approaches.

Methods
Metazoan single-copy orthologs (mzl-USCOs) were sequenced in 171 Pachypus individuals by targeted enrichment and mapped against the reference sequences with bwa. Phylogenetic analyses were conducted with IQ-TREE and ASTRAL. Time-calibrated phylogenetic trees were constructed with MCMCTree. Species delimitation was performed with BPP, tr2, SODA, GMYC and PTP. Population clustering was performed with ADMIXTURE, and intra-population differentiation and gene flow were analysed with ADMIXTOOLS and Dsuite.