Cryptic species can be old despite extreme female philopatry
Abstract
The impact of strong differentiation between populations due to limited or sex-biased dispersal on species delimitation 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 extremely interesting case study. Based on a dataset of over 900 protein-coding genes (metazoan universal single-copy orthologs; mzl-USCOs), we investigated its phylogeny and species boundaries and explored gene flow and population differentiation to allow an integrative assessment of species boundaries. Integrative consideration of all results leads 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. The considerable age of most 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 between populations of the same species. Speciation was partly associated with a disruption of formerly likely more connected ranges after the Messinian salinity crisis (MSC). This also helps to answer to what extent the Mediterranean dried out during the MSC as land connections in the circum-Tyrrhenian region must have existed for longer time intervals in which flightless Pachypus females could disperse between drifting land areas reaching 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. Our study illustrates how species delimitation analyses of study cases with restricted or sex-biased dispersal and highly differentiated populations are useful as empirical test cases for the utility and robustness of species delimitation approaches.
Full text
Figure S1. Collection localities of specimens analyzed in this study.
Figure S2a. Coalescent-based tree calculated with ASTRAL based on individual USCO gene trees of 171 Pachypus individuals and outgroup taxa. Numbers above branches are local posterior probabilities. 0.6 DA3382_L37 836892_L13 DA2567_L2 DA3326_L33 DA4238_L59 DA2606_L15 DA2594_L15 DA2573_L1 DA3383_L37 DA4248_L52 DA4266_L49 836889_L24 DA2564_L23 DA3356_L35 DA2656_L21 DA2810_L18 DA4233_L57 DA3336_L34 DA2506_L25 DA4261_L51 DA2725_L14 DA4251_L53 DA2657_L21 DA2659_L17 DA3422_L39 836893_L13 JE0062_L41 Sparrmannia DA4257_L51 DA3481_L44 DA3373_L36 DA2585_L22 DA2570_L8 DA3327_L33 DA4265_L49 DA3338_L34 DA2568_L8 DA2626_L16 DA2740_L20 DA4254_L50 837890_L10 DA4242_L60 DA4079_L55 834183_L20 DA4263_L49 DA3999_L48 DA4235_L58 DA2794_L19 DA4262_L51 836895_L13 DA4236_L58 DA2730_L6 DA3415_L39 DA3337_L34 DA4239_L59 DA3375_L36 836890_L13 DA2595_L15 DA4232_L57 837884_L6 837885_L6 DA2729_L6 836897_L13 836880_L12 DA4081_L55 836887_L24 836899_L13 DA3331_L33 DA4273_L61 DA2572_L0 DA3973_L23 DA3330_L33 DA2784_L19 DA2577_L11 DA2589_L15 DA3355_L35 DA2597_L15 DA2583_L9 DA3997_L48 DA4237_L58 DA3423_L39 DA4241_L60 DA2719_L19 DA3419_L39 DA3395_L38 837883_L6 836891_L13 DA3421_L39 836888_L24 DA3396_L38 DA4258_L51 DA2779_L18 DA4230_L56 DA2778_L18 DA2505_L26 DA4256_L50 DA2581_L11 DA2562_L27 DA3417_L39 DA2588_L15 DA4267_L49 DA3394_L38 DA3418_L39 836885_L7 DA2586_L15 DA3374_L36 DA2560_L25 DA2632_L16 DA2582_L9 DA4255_L50 836898_L13 DA4078_L55 DA4272_L61 836900_L13 DA2739_L20 DA2726_L14 DA2776_L18 DA2565_L13a JE0061_L41 DA4234_L58 DA4229_L56 836882_L3 DA4076_L55 DA4080_L55 DA2563_L23 DA3974_L23 DA4260_L51 DA4077_L55 837887_L10 DA3393_L38 DA3420_L39 DA4075_L55 DA3416_L39 DA4240_L60 DA3328_L33 836901_L13 DA3998_L48 DA2575_L1 DA4249_L52 DA2580_L11 DA2578_L11 DA2569_L8 DA2571_L0 DA3996_L48 836881_L4 Phyllotocus DA2561_L25 DA2598_L15 DA4253_L50 DA4250_L53 DA2709_L19 836883_L3 DA3329_L33 836894_L13 DA2636_L16 DA3385_L37 DA4252_L50 DA3392_L38 DA4264_L49 836884_L7 DA3384_L37 837886_L6 DA3354_L35 DA2566_L13a DA2742_L20 DA2579_L11 DA3971_L47 DA4231_L56 836886_L5 836896_L13 DA2639_L16 DA2658_L17 DA4259_L51 0.67 1 0.73 0.95 0.99 1 0.58 0.44 1 0.67 1 1 1 1 0.54 0.44 0.99 1 0.94 1 1 1 0.48 0.58 0.65 0.61 1 1 1 0.98 1 0.81 1 0.84 0.56 1 0.75 0.54 1 0.52 1 1 1 1 0.97 0.49 0.99 1 0.48 1 0.78 1 1 0.84 1 0.97 0.5 0.41 1 0.8 0.55 0.51 1 1 0.91 0.85 0.47 0.99 1 1 0.44 1 1 1 1 1 1 1 1 0.92 1 1 1 0.64 0.99 0.78 1 1 0.94 1 0.92 0.83 1 1 0.65 0.37 0.9 0.73 0.96 0.75 1 0.38 1 0.42 0.99 0.72 1 1 0.68 0.47 0.45 0.56 1 0.52 0.95 1 1 1 0.81 0.68 0.85 0.48 0.6 0.78 0.91 1 0.37 0.84 1 1 1 0.48 0.43 1 0.96 0.78 0.45 1 0.42 0.39 0.8 0.52 0.97 0.94 0.92 0.77 1 1 1 1 0.61 1 1 0.8 0.64 0.96 1 0.44 0.75 1 1 1 1 1 1 0.93 0.39 0.51 0.46 0.98 1 sp. 3 sp. 2 (North) sp. 4 sp. 2a excavatus sp. 5 sardiniensis sp. 1 caesus demoflysi melonii1 sp. 6 candidae outgroups melonii2 melonii3
Figure S2b. Coalescent-based tree calculated with ASTRAL based on individual USCO gene trees of 171 Pachypus individuals and outgroup taxa. Numbers above branches are quartet scores. 0.6 DA3393_L38 DA4259_L51 DA2588_L15 DA4242_L60 DA2740_L20 DA2571_L0 DA2726_L14 DA4252_L50 DA4262_L51 DA4255_L50 DA2561_L25 DA4235_L58 DA3419_L39 DA2560_L25 JE0061_L41 DA4075_L55 DA3971_L47 DA4229_L56 DA4237_L58 837887_L10 836883_L3 836897_L13 DA3418_L39 DA4263_L49 DA3385_L37 DA2581_L11 DA2778_L18 836886_L5 Phyllotocus DA3354_L35 836888_L24 DA4267_L49 834183_L20 DA2577_L11 DA4234_L58 837886_L6 DA3998_L48 836891_L13 DA3999_L48 DA2729_L6 DA4232_L57 836896_L13 DA2505_L26 DA3974_L23 DA2579_L11 836889_L24 DA2572_L0 DA4241_L60 DA3355_L35 DA3996_L48 DA2566_L13a DA2565_L13a DA2709_L19 DA2725_L14 DA3416_L39 837883_L6 DA2594_L15 DA3375_L36 DA2779_L18 836893_L13 DA2597_L15 DA3329_L33 DA3997_L48 DA4260_L51 DA2632_L16 DA3383_L37 DA4081_L55 DA2580_L11 DA4238_L59 DA4079_L55 DA3382_L37 DA4078_L55 DA2598_L15 Sparrmannia DA2658_L17 837885_L6 DA3396_L38 DA2563_L23 DA4077_L55 836890_L13 DA2589_L15 836894_L13 DA2586_L15 DA3973_L23 DA2742_L20 DA2575_L1 DA2567_L2 836900_L13 DA2506_L25 DA4253_L50 836882_L3 DA2562_L27 836881_L4 DA2776_L18 DA2569_L8 DA3422_L39 DA4254_L50 DA2657_L21 836898_L13 DA3327_L33 DA4231_L56 DA3423_L39 DA4076_L55 DA2810_L18 DA3338_L34 DA2784_L19 DA4233_L57 DA2636_L16 DA2578_L11 DA2656_L21 DA2739_L20 DA4230_L56 836885_L7 DA2564_L23 DA4273_L61 DA2639_L16 DA4251_L53 DA3330_L33 DA3415_L39 DA3392_L38 DA2582_L9 DA2585_L22 837890_L10 836901_L13 836895_L13 DA2595_L15 DA4261_L51 DA2606_L15 DA4239_L59 DA4266_L49 DA2794_L19 DA4272_L61 DA4257_L51 836880_L12 DA3394_L38 DA3356_L35 DA4080_L55 DA3384_L37 DA2573_L1 DA4240_L60 DA3331_L33 DA2730_L6 DA3420_L39 DA2659_L17 DA4250_L53 DA3328_L33 836884_L7 DA3421_L39 DA4265_L49 DA4236_L58 DA2719_L19 837884_L6 JE0062_L41 DA3481_L44 836887_L24 DA3326_L33 DA2583_L9 836899_L13 DA3374_L36 DA3336_L34 DA4249_L52 DA3395_L38 DA3417_L39 DA4258_L51 DA2570_L8 DA2568_L8 DA4256_L50 DA4248_L52 836892_L13 DA3337_L34 DA2626_L16 DA3373_L36 DA4264_L49 45.03 95.34 34.91 47.04 34.48 51.41 42.31 36.02 35.56 44.91 54.47 35.85 46.03 34.27 35.78 38.51 36.36 38.08 34.51 34.22 41.28 68.73 37.58 34.13 45.26 34.13 38.39 40.53 37.28 35.66 37.76 33.65 37.61 71.8 34.95 59.94 36.08 36.09 67.6 43.81 35.54 34.43 38.38 33.99 39.4 42.42 71.33 57.63 36.12 34.91 38.52 36.19 37.3 38.77 42.71 39.57 34.92 43.69 36.71 56.48 36.82 78.3 36.76 61.37 35.05 34.09 36.93 78.77 33.85 35.31 39.26 41.14 79.3 37.57 34.68 34.36 33.98 41.8 35.96 39.37 37.01 35.44 38.2 37.79 33.53 35.93 39.13 33.91 40.35 33.63 37.1 33.67 69.56 34.39 40.33 86.81 37.27 36.32 35.93 39.74 34.84 36.02 39.27 36 39.19 77.63 37.06 49 37.28 34.98 57.66 34.96 44.3 34.22 46.53 34.56 39.74 54.65 34.65 67.66 89.22 34.31 42.45 42.5 35.99 34.5 34.04 52.32 33.56 70.71 59.35 42.79 37.18 33.89 35.86 35.42 95.34 34.9 58.68 34.33 55.34 53.04 35.28 33.94 36.24 41.87 41.68 41.04 41.18 48.93 36.34 36.39 56.26 38.05 36.71 50.45 35.81 36.77 41.21 35.38 36.54 43.49 34.3 62.94 40.52 78.92 37.74 59.52 33.8 35.35 33.93 outgroups candidae sp. 6 melonii1 melonii2 melo nii3 demoflysi caesus sardiniensis sp. 1 sp. 5 excavatus sp. 2a sp. 4 sp. 2b sp. 3
Figure S3. Maximum-likelihood tree based on concatenated USCO data of all 171 Pachypus individuals and outgroup taxa. Numbers above branches are support values from approximate likelihood ratio tests and ultrafast bootstrapping. 0.004 DA4231_L56 DA2709_L19 836884_L7 DA4259_L51 DA3354_L35 DA3330_L33 DA2567_L2 DA2561_L25 DA2581_L11 DA4254_L50 DA2571_L0 DA3326_L33 DA2563_L23 DA4237_L58 DA3355_L35 DA2739_L20 DA3415_L39 DA3419_L39 DA2778_L18 DA4242_L60 DA4258_L51 836891_L13 DA2566_L13a DA2779_L18 836890_L13 836894_L13 836885_L7 DA3338_L34 836893_L13 DA4263_L49 DA2657_L21 DA2725_L14 DA3999_L48 DA3393_L38 836887_L24 DA2598_L15 DA3997_L48 DA2632_L16 DA4266_L49 DA4240_L60 Sparrmannia DA3373_L36 DA3416_L39 DA2585_L22 DA2730_L6 DA2568_L8 DA4233_L57 DA4076_L55 DA3329_L33 DA4078_L55 DA2656_L21 DA4256_L50 837883_L6 Phyllotocus DA2719_L19 DA4257_L51 DA2626_L16 DA2589_L15 DA3336_L34 DA2594_L15 DA3327_L33 836895_L13 837885_L6 836898_L13 DA2582_L9 DA4248_L52 DA4081_L55 DA4234_L58 DA4232_L57 836900_L13 DA3375_L36 DA4229_L56 DA4238_L59 DA2586_L15 836889_L24 DA3394_L38 DA4273_L61 DA2606_L15 DA3418_L39 DA4264_L49 DA4241_L60 DA3383_L37 DA2569_L8 DA4236_L58 DA3974_L23 DA3374_L36 DA4230_L56 JE0061_L41 DA3417_L39 DA4251_L53 DA2575_L1 DA2740_L20 DA3356_L35 DA4249_L52 836892_L13 DA2595_L15 DA4252_L50 DA3392_L38 DA3422_L39 DA4079_L55 DA2579_L11 DA4261_L51 DA2639_L16 DA4239_L59 DA3337_L34 DA2564_L23 836896_L13 DA2588_L15 DA4260_L51 DA2560_L25 DA4253_L50 DA2729_L6 836883_L3 DA4250_L53 DA3973_L23 DA2726_L14 DA2580_L11 836897_L13 DA3420_L39 834183_L20 DA2810_L18 836888_L24 DA4075_L55 DA2505_L26 DA2573_L1 837886_L6 DA2583_L9 DA4255_L50 DA4265_L49 DA4077_L55 836901_L13 DA2578_L11 DA3328_L33 836882_L3 DA2565_L13a DA4262_L51 DA3423_L39 DA2742_L20 DA2562_L27 DA2794_L19 DA3385_L37 DA4267_L49 DA2570_L8 837890_L10 DA3382_L37 DA3971_L47 DA2659_L17 DA2577_L11 DA3331_L33 836880_L12 DA2572_L0 DA4272_L61 DA2597_L15 DA3395_L38 DA2636_L16 DA2506_L25 DA3481_L44 DA3998_L48 836881_L4 DA2658_L17 DA3996_L48 DA3421_L39 837884_L6 DA3396_L38 JE0062_L41 DA4080_L55 DA4235_L58 836886_L5 DA2776_L18 DA3384_L37 836899_L13 DA2784_L19 837887_L10 100/100 99.1/100 100/100 100/100 100/100 100/100 100/100 100/100 98.8/97 100/100 92.1/65 44/39 95.3/80 92.7/85 100/100 92/53 100/100 99.9/92 100/100 92.8/70 100/97 100/100 100/100 66.4/71 100/100 97.4/81 88.6/67 98.7/86 12.9/46 73/79 40/52 30.6/49 100/100 68.2/75 100/100 100/100 91.7/47 100/100 94.7/90 100/99 100/100 100/100 100/100 100/100 100/100 97.1/91 99.9/67 34.9/51 34.5/57 95.4/92 100/78 95.6/93 95.2/92 97.8/95 83.6/46 100/100 74.7/27 99.2/60 100/100 100/82 100/100 100/100 98/78 99.8/73 100/100 74/75 100/100 96.3/89 82.8/66 99.9/54 92.8/88 100/100 18/49 99.9/99 97.2/82 99.8/100 92.8/75 100/100 100/100 80.3/61 23.7/58 80.3/83 95.7/94 86.5/88 100/100 81.9/80 95/50 100/100 12.8/44 100/100 100/100 100/100 100/100 100/100 100/100 100/100 95.8/93 97.7/66 100/100 92.3/65 98.7/97 100/100 100/100 70/74 13.9/8 98/93 100/100 95.9/83 100/100 77.8/73 100/100 75.7/44 89.5/61 45.8/52 69.9/70 100/100 100/100 99.3/91 81.5/65 97.8/96 69.4/74 100/100 99/57 100/100 79.7/67 99.9/99 100/81 100/100 93.5/89 99.5/94 100/96 100/100 50.9/51 100/97 100/100 100/100 100/100 99.2/36 100/100 100/100 100/100 100/100 9/53 100/100 100/90 43.2/67 99.3/92 100/100 10.5/51 99.9/99 100/78 61.9/69 36.8/43 100/100 100/100 100/100 99.8/98 10.6/42 55.5/72 100/100 99.8/91 100/100 100/100 24/36 99.6/99 100/100 97.9/95 93.2/55 100/100 66/59 90.9/58 outgroups candidae sp. 6 melonii3 demoflysi caesus sardiniensis sp. 1 sp. 5 excavatus sp. 2b sp. 2a sp. 4 sp. 3 melonii2 melonii1
Figure S4. Phylogenetic NeighborNet network of 171 Pachypus individuals based on SNP data.
Figure S5. Maximum-likelihood tree based on concatenated alignment of selected Pachypus individuals and a larger dataset of other scarabaeoid beetles, aligned with hmmalign and including all nucleotide positions. Numbers above branches are ultrafast bootstrap values. 0.2 Sericoides Heterochelus_sp Pachypus_melonii_DA4236 Lucanus_tetradon Scitalia Aphodius_scrutator Euoniticellus_fulvus Pachypus_sp2_DA2578 Trochalus_sp Valgus_hemipterus Liparetus Melinopterus_prodromus Onthophagus_similis Pachypus_caesus_DA3481 Platycerus_caraboides Helophorus_nanus Sisyphus_schaefferi Ablabera_sp Amphimallon_solstitiale Phyllotocus Geotrupes_spiniger Pachypus_sp3_DA2776 Astenopholis_sp Trichius_sp2 Orphnus Trichius_sp1 Sparrmannia_sp Sinodendron_cylindricum Pachypus_sp6_DA4081 Diphucephala Eulasia_bombylius Pleurophorus_caesus Pachypus_sp1_DA4249 Maechidius Thorectes_sp Hoplia_graminea Serica_brunnea Adoretus_sp Pachypus_sp4_DA4229 Pleophylla_fasciatipennis Pachypus_sp2_DA3385 Pachypus_excavatus_DA2582 Bolboceras Lepiserica_sp Hybosorus_sp Triodontella Ocypus_brunnipes Copris_luniaris Anoplotrupes_stercorosus Dorcus_parallelipedus Camenta_innocua1 Telura Sophrops Neoserica_sp Aphodius_haemorhoidalis Pentodon_sp Heteronyx Anomala_sp Cetonia_aurata_pisana Trox_sp Odontotaenius_disjunctus Dynastinae_sp Hoplia_philanthus Pachypus_sardiniensis_DA3972 Aphodius_rufus Pachypus_candidae_DA2571 Camenta_innocua2 Pachypus_demoflysi_DA3998 Allokotarsa_sp Glaresis_sp 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 66 100 100 63 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
Figure S6. Maximum-likelihood tree based on concatenated alignment of selected Pachypus individuals and a larger dataset of other scarabaeoid beetles, aligned with hmmalign and including only first and second nucleotide positions. Numbers above branches are ultrafast bootstrap values. 0.05 Pachypus_sp3_DA2776 Neoserica_sp Ablabera_sp Bolboceras Pachypus_sp2_DA2578 Telura Dorcus_parallelipedus Pachypus_sp2_DA3385 Adoretus_sp Ocypus_brunnipes Pachypus_candidae_DA2571 Geotrupes_spiniger Serica_brunnea Pleophylla_fasciatipennis Trochalus_sp Heterochelus_sp Anoplotrupes_stercorosus Euoniticellus_fulvus Platycerus_caraboides Copris_luniaris Sparrmannia_sp Aphodius_scrutator Anomala_sp Pachypus_sp4_DA4229 Sophrops Triodontella Pachypus_caesus_DA3481 Diphucephala Pentodon_sp Maechidius Melinopterus_prodromus Trox_sp Astenopholis_sp Heteronyx Helophorus_nanus Trichius_sp1 Liparetus Sisyphus_schaefferi Odontotaenius_disjunctus Pachypus_melonii_DA4236 Valgus_hemipterus Trichius_sp2 Aphodius_rufus Eulasia_bombylius Glaresis_sp Thorectes_sp Allokotarsa_sp Camenta_innocua2 Lucanus_tetradon Aphodius_haemorhoidalis Pachypus_sp6_DA4081 Pachypus_demoflysi_DA3998 Pachypus_excavatus_DA2582 Dynastinae_sp Sinodendron_cylindricum Hoplia_philanthus Orphnus Hybosorus_sp Lepiserica_sp Hoplia_graminea Pachypus_sp1_DA4249 Cetonia_aurata_pisana Amphimallon_solstitiale Camenta_innocua1 Sericoides Pachypus_sardiniensis_DA3972 Onthophagus_similis Phyllotocus Scitalia Pleurophorus_caesus 100 100 99 100 100 100 100 100 100 100 100 47 91 100 100 100 100 100 54 100 98 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 96 100 100 100 100 91 100 100 100 99 47 96 100 100 87 100 87 100 100 100 100 100 100 100 100 100
Figure S7. Maximum-likelihood tree based on concatenated alignment of selected Pachypus individuals and a larger dataset of other scarabaeoid beetles, aligned with MAFFT and including all nucleotide positions. Numbers above branches are ultrafast bootstrap values. 0.2 Heterochelus_sp Camenta_innocua2 Sericoides Aphodius_rufus Pachypus_sp1_DA4249 Anoplotrupes_stercorosus Hoplia_philanthus Platycerus_caraboides Pachypus_excavatus_DA2582 Euoniticellus_fulvus Ocypus_brunnipes Valgus_hemipterus Neoserica_sp Pachypus_sardiniensis_DA3972 Sinodendron_cylindricum Anomala_sp Diphucephala Cetonia_aurata_pisana Trichius_sp1 Lepiserica_sp Telura Trochalus_sp Odontotaenius_disjunctus Maechidius Triodontella Amphimallon_solstitiale Pachypus_sp2_DA3385 Hoplia_graminea Astenopholis_sp Pachypus_sp6_DA4081 Trichius_sp2 Sophrops Liparetus Sparrmannia_sp Orphnus Pleurophorus_caesus Dynastinae_sp Sisyphus_schaefferi Phyllotocus Heteronyx Trox_sp Helophorus_nanus Pachypus_sp2_DA2578 Pachypus_sp4_DA4229 Scitalia Melinopterus_prodromus Onthophagus_similis Pachypus_candidae_DA2571 Aphodius_scrutator Eulasia_bombylius Pachypus_sp3_DA2776 Pachypus_demoflysi_DA3998 Thorectes_sp Pentodon_sp Pachypus_caesus_DA3481 Hybosorus_sp Adoretus_sp Allokotarsa_sp Lucanus_tetradon Glaresis_sp Pachypus_melonii_DA4236 Copris_luniaris Geotrupes_spiniger Bolboceras Ablabera_sp Pleophylla_fasciatipennis Camenta_innocua1 Dorcus_parallelipedus Serica_brunnea Aphodius_haemorhoidalis 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 49 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
Figure S8. Maximum-likelihood tree based on concatenated alignment of selected Pachypus individuals and a larger dataset of other scarabaeoid beetles, aligned with MAFFT and including only first and second nucleotide positions. Numbers above branches are ultrafast bootstrap values. 0.05 Pentodon_sp Anoplotrupes_stercorosus Trichius_sp2 Helophorus_nanus Anomala_sp Lucanus_tetradon Dynastinae_sp Orphnus Heterochelus_sp Pleophylla_fasciatipennis Pachypus_sardiniensis_DA3972 Sericoides Bolboceras Neoserica_sp Heteronyx Trichius_sp1 Diphucephala Copris_luniaris Pachypus_melonii_DA4236 Pachypus_excavatus_DA2582 Amphimallon_solstitiale Pachypus_sp1_DA4249 Camenta_innocua1 Pachypus_sp4_DA4229 Liparetus Sisyphus_schaefferi Adoretus_sp Cetonia_aurata_pisana Aphodius_haemorhoidalis Hoplia_graminea Odontotaenius_disjunctus Dorcus_parallelipedus Pachypus_sp6_DA4081 Astenopholis_sp Allokotarsa_sp Melinopterus_prodromus Hybosorus_sp Lepiserica_sp Sparrmannia_sp Sophrops Pachypus_sp2_DA2578 Hoplia_philanthus Camenta_innocua2 Geotrupes_spiniger Euoniticellus_fulvus Ocypus_brunnipes Serica_brunnea Valgus_hemipterus Sinodendron_cylindricum Platycerus_caraboides Eulasia_bombylius Ablabera_sp Thorectes_sp Phyllotocus Onthophagus_similis Pachypus_caesus_DA3481 Pleurophorus_caesus Pachypus_candidae_DA2571 Trox_sp Telura Aphodius_rufus Triodontella Maechidius Pachypus_sp2_DA3385 Trochalus_sp Aphodius_scrutator Pachypus_sp3_DA2776 Scitalia Glaresis_sp Pachypus_demoflysi_DA3998 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 57 100 100 100 100 100 100 100 100 100 100 100 100 84 100 100 100 100 98 100 100 100 100 100 100 100 100 100 100 100 100 100 100 89 100 89 100 100 100 100 100 100 100 100 100 100 95 100 100 100 100
Figure S15. Tree calibrated with MCMCTree based on fixed ML tree topology using the calibration scheme from McKenna et al. (2019). Ado retus sp. Odontotaenius disjunctus Diphucephala sp. Trox sp. Pachypus candidae Pleurophorus caesus Pachypus sp. 3 Eulasia bombylius Ocypus brunnipes Pachypus melonii3 Ablabera sp. Copris lunaris Hoplia philanthus Heteronyx sp. Coprimorphus scrutator Euoniticellus fulvus Trichius sp. 2 Geotrupes sp.iniger Camenta innocua 1 Lepiserica sp. Bolboceras sp. Pachypus sp. 1 Orphnus sp. Sophrops sp. Hybosorus sp. Sparrmannia sp. Bodilopsis rufa Hoplia graminicola Pachypus sp. 5b Sisyphus schaefferi Amphimallon solstitiale Sinodendron cylindricum Cetonia aurata pisana Astenopholis sp. Scitalia sp. Melinopterus prodromus Camenta innocua 2 Sericoides sp. Pleophylla fasciatipennis Pachypus caesus Pachypus sp. 2b Telura sp. Helophorus nanus Valgus hemipterus Phyllotocus sp. Triodontella sp. Dorcus parallelepipedus Pachypus sp. 6 Neoserica sp. Glaresis sp. Otophorus haemorrhoidalis Anoplotrupes stercorosus Liparetrus sp. Allokotarsa sp. Pachypus sp. 4 Pachypus excavatus Onthophagus similis Trichius sp. 1 Pachypus demoflysi Pachypus sardiniensis Heterochelus sp. Pachypus melonii1 Serica brunnea Thorectes sp. Maechidius sp. Anomala sp. Pachypus melonii2 Pentodon sp. Eophileurus sp. Platycerus caraboides Pachypus sp. 5a Pachypus sp. 2a Lucanus tetraodon Trochalus sp. 0.3 A BC D E F G H L R S U 0 10 20 30 40 5060 70 8090 100110120 130140150 160170180190 200210220 230 Upper Lower Middle Upper Lower Upper Paleocene Eocene Oligocene Miocene Pliocene Triassic Jurassic Cretaceous Paleogene Neogene
Figure S16. Results of BPP species delimitation for subclades with the full dataset (left side) and using a dataset where all positions containing missing data and/or gaps were removed (right side) with different theta and tau values. Boxes at tree nodes refer to different combinations of theta and tau, and colors refer to posterior probabilities of species splits at the node (see inset), Subclades: A – A1, B – A2, C – B2, D – excavatus, E – sp. 2/4, F – sp. 3.
Figure S17. Results of various species delimitation analyses mapped on ASTRAL tree. Colored bars (bars 3-10) show entities inferred to be distinct species separated by horizontal lines. Curved lines connect specimens assigned to the same species-level entity. From left to right these are: 1) species based on Eberle et al. (2018); 2) groupings (i.e., subclades) for BPP analysis; results of 3) BPP with full dataset; 4) BPP with dataset where all positions containing missing data and/or gaps were removed; 5) BPP with full dataset, where lineages not supported as full species by gdi were lumped (grey lines are ambiguous cases; numbers are gdi values separating a given lineage from its sisters); 6) bPTP result with maximum likelihood; 7) bPTP result with highest Bayesian support; 8) GMYC; 9) tr2; 10) SODA. 0.4 DA2562_L27 DA3416_L39 DA4272_L61 836888_L24 836890_L13 836884_L7 DA4076_L55 DA3375_L36 DA2659_L17 DA3421_L39 DA2585_L22 DA2709_L19 DA3373_L36 DA4257_L51 DA2810_L18 DA2784_L19 DA2656_L21 DA3418_L39 DA2726_L14 DA3973_L23 DA4253_L50 DA4241_L60 DA3997_L48 DA4239_L59 DA2636_L16 DA4233_L57 DA3419_L39 837884_L6 DA2740_L20 DA4263_L49 DA3998_L48 DA2588_L15 DA3328_L33 DA4251_L53 DA2779_L18 DA2505_L26 DA3422_L39 DA2571_L0 DA2577_L11 DA3385_L37 DA2560_L25 DA3393_L38 836885_L7 DA2567_L2 836887_L24 836897_L13 DA3382_L37 DA2658_L17 DA2580_L11 836893_L13 837890_L10 DA3417_L39 836896_L13 DA4265_L49 DA2572_L0 DA2776_L18 836880_L12 DA4252_L50 DA3338_L34 DA4080_L55 836894_L13 DA2582_L9 DA4260_L51 DA3326_L33 836900_L13 DA4242_L60 DA4248_L52 DA2583_L9 DA4240_L60 DA2579_L11 DA4273_L61 DA4261_L51 DA2778_L18 DA3337_L34 DA2566_L13a DA3384_L37 DA4235_L58 DA3392_L38 DA4237_L58 DA3394_L38 DA4230_L56 DA2730_L6 834183_L20 DA4258_L51 DA3336_L34 DA2563_L23 DA2569_L8 DA2589_L15 DA3383_L37 DA2578_L11 DA2657_L21 836899_L13 DA2725_L14 DA2570_L8 836901_L13 836886_L5 DA4255_L50 836898_L13 DA2626_L16 DA3355_L35 DA4266_L49 836881_L4 DA4250_L53 836895_L13 837887_L10 DA2729_L6 DA3327_L33 DA4236_L58 DA4259_L51 DA2742_L20 837885_L6 DA3374_L36 JE0062_L41 DA2719_L19 DA4232_L57 DA4231_L56 DA2639_L16 DA2632_L16 DA2568_L8 DA4254_L50 DA3331_L33 DA2595_L15 DA3395_L38 DA2565_L13a DA3329_L33 DA3354_L35 DA2586_L15 DA4077_L55 DA2794_L19 DA4078_L55 JE0061_L41 DA4249_L52 DA3420_L39 DA4079_L55 DA2597_L15 DA2594_L15 837886_L6 DA4238_L59 DA4256_L50 DA2739_L20 DA2606_L15 836882_L3 DA4075_L55 DA4081_L55 836889_L24 DA3356_L35 DA3974_L23 DA2598_L15 837883_L6 DA3971_L47 DA2573_L1 836883_L3 DA3396_L38 DA4267_L49 836892_L13 DA4229_L56 DA2575_L1 DA3999_L48 DA3330_L33 DA4234_L58 836891_L13 DA3481_L44 DA2564_L23 DA4264_L49 DA4262_L51 DA2561_L25 DA2506_L25 DA3996_L48 DA3423_L39 DA2581_L11 DA3415_L39 candidae sp. 6 melonii demoflysi caesus sardiniensis sp. 1 sp. 5 excavatus sp. 2a sp. 4 sp. 2b sp. 3 species B2 B1-excavatus B1sp. 2/4 B1 - sp. 3 groups for BPP analysis A1 A2 BPP all data BPP no gaps tr2 BPP+gdi no gaps SODAGMYC bPTP ML bPTP Bayes 0.72 0.87 0.72 0.57 0.72 0.72
Figure S18. Values of the genealogical diversity index (gdi) values mapped onto the respective nodes of the ASTRAL tree of all Pachypus specimens (for BPP priors: beta = 0.04 for theta, beta = 0.02 for tau). Critical gdi values above merging threshold and below splitting threshold are indicated by green and red squares, respectively. Nodes indicated by an empty circle were not subject to gdi analysis due to dataset subdivision (see text). Inset shows the frequency distribution of gdi values. 836892_L13 DA2567_L2 DA3326_L33 DA4238_L59 DA2606_L15 DA2594_L15 DA2573_L1 DA3383_L37 DA4248_L52 DA4266_L49 836889_L24 DA2564_L23 DA3356_L35 DA2656_L21 DA2810_L18 DA4233_L57 DA3336_L34 DA2506_L25 DA4261_L51 DA2725_L14 DA4251_L53 DA2657_L21 DA2659_L17 DA3422_L39 836893_L13 JE0062_L41 Sparrmannia DA4257_L51 DA3481_L44 DA3373_L36 DA2585_L22 DA2570_L8 DA3327_L33 DA4265_L49 DA3338_L34 DA2568_L8 DA2626_L16 DA2740_L20 DA4254_L50 837890_L10 DA4242_L60 DA4079_L55 834183_L20 DA4263_L49 DA3999_L48 DA4235_L58 DA2794_L19 DA4262_L51 836895_L13 DA4236_L58 DA2730_L6 DA3415_L39 DA3337_L34 DA4239_L59 DA3375_L36 836890_L13 DA2595_L15 DA4232_L57 837884_L6 837885_L6 DA2729_L6 836897_L13 836880_L12 DA4081_L55 836887_L24 836899_L13 DA3331_L33 DA4273_L61 DA2572_L0 DA3973_L23 DA3330_L33 DA2784_L19 DA2577_L11 DA2589_L15 DA3355_L35 DA2597_L15 DA2583_L9 DA3997_L48 DA4237_L58 DA3423_L39 DA4241_L60 DA2719_L19 DA3419_L39 DA3395_L38 837883_L6 836891_L13 DA3421_L39 836888_L24 DA3396_L38 DA4258_L51 DA2779_L18 DA4230_L56 DA2778_L18 DA2505_L26 DA4256_L50 DA2581_L11 DA2562_L27 DA3416_L39 DA2588_L15 DA4267_L49 DA3394_L38 DA3418_L39 836885_L7 DA2586_L15 DA3374_L36 DA2560_L25 DA2632_L16 DA2582_L9 DA4255_L50 836898_L13 DA4078_L55 DA4272_L61 836900_L13 DA2739_L20 DA2726_L14 DA2776_L18 DA2565_L13a JE0061_L41 DA4234_L58 DA4229_L56 836882_L3 DA4076_L55 DA4080_L55 DA2563_L23 DA3974_L23 DA4260_L51 DA4077_L55 837887_L10 DA3393_L38 DA3420_L39 DA4075_L55 DA3417_L39 DA4240_L60 DA3328_L33 836901_L13 DA3998_L48 DA2575_L1 DA4249_L52 DA2580_L11 DA2578_L11 DA2569_L8 DA2571_L0 DA3996_L48 836881_L4 Phyllotocus DA2561_L25 DA2598_L15 DA4253_L50 DA4250_L53 DA2709_L19 836883_L3 DA3329_L33 836894_L13 DA2636_L16 DA3385_L37 DA4252_L50 DA3392_L38 DA4264_L49 836884_L7 DA3384_L37 837886_L6 DA3354_L35 DA2566_L13a DA2742_L20 DA2579_L11 DA3971_L47 DA4231_L56 836886_L5 836896_L13 DA2639_L16 DA2658_L17 DA4259_L51 0.65 0.36 0.6 0.22 0.66 0.78 0.89 0.29 0.72 0.99 0.36 0.19 0.79 0.79 0.72 0.91 0.57 0.5 0.28 0.59 0.57 0.02 0.36 0.38 0.6 0.32 0.07 0.94 0.5 0.04 0.49 1 0.1 0.83 0.41 0.19 0.42 0.98 0.4 0.09 0.51 0.04 0.72 0.06 0.66 0.44 0.89 0.39 0.81 0.06 0 0.33 0.6 0.39 0.56 0.72 0.87 1 0.15 0.04 0.51 0.6 DA3382_L37 0 1 2 3 4 5 6 7 8 9 10 1-0.91 Frequency gdi intervals gdi >0.7 gdi <0.2 0.2-0.11 0.1-0.01 0 0.3-0.21 0.4-0.31 0.5-0.41 0.6-0.51 0.7-0.61 0.8-0.71 0.9-0.81 sp. 3 sp. 2b (North) sp. 4 sp. 2a excavatus sp. 5 sardiniensis sp. 1 caesus demoflysi melonii1 sp. 6 candidae outgroups melonii2 melonii3
Figure S19. Partitioned NMDS analyses using SNP data of the major Pachypus clades, A - A1, B - A2, C - B1, D - B2. -0.375 -0.300 -0.225 -0.150 -0.075 0.000 0.075 0.150 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 Coordinate 2 caesus demoflysi sardiniensis candidae melonii1 sp. 6 melonii2 melonii3 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 excavatus sp. 2a sp. 3 sp. 4 sp. 2b -0.300-0.225 -0.150 -0.075 0.000 0.075 0.150 0.225 0.300 0.375 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 sp. 1 sp. 5a (Corsica) sp. 5a (Sardinia) sp. 5b Coordinate 2 Coordinate 1 Coordinate 1 A B D C
Figure S20. Results of ADMIXTURE analysis based on SNP data including all individuals for a number of ancestral populations (K) of 12. 0.0 0.2 0.4 0.6 0.8 1.0 caesus candidae demoflysi excavatus melonii sardiniensis sp. 1 sp. 2b sp. 2a sp. 3 sp. 4 sp. 5 sp. 6 melonii2 melonii3 melonii1
Figure S21. Partitioned ADMIXTURE analyses based on SNP data for the major Pachypus clades, A - A1, B - A2, C - B1, D - B2. Ancestry 0.0 0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0 caesus demoflysi sardiniensis caesus demoflysi sardiniensis K=3 K=4 Ancestry 0.0 0.2 0.4 0.6 0.8 1.0 candidae sp. 6 melonii1 K=5 0.0 0.2 0.4 0.6 0.8 1.0 candidae sp. 6 K=6 Ancestry K=9 K=10 K=11 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 excavatus sp. 2a sp. 3 sp. 4 sp. 2b excavatus sp. 2a sp. 3 sp. 4 sp. 2b excavatus sp. 2a sp. 3 sp. 4 sp. 2b Ancestry 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 sp. 1 sp. 5a sp. 5b sp. 1 sp. 5a sp. 5b sp. 1 sp. 5a sp. 5b K=2 K=3 K=4 0.0 0.2 0.4 0.6 0.8 1.0 A B D C excavatus sp. 4sp. 3sp. 2a sp. 2b melonii3melonii2 melonii1 melonii3melonii2
! Figure S22. Comparison of cross-validation with ADMIXTURE in Pachypus and its major subclades, compared to other case studies using mzl-USCOs, with the x-axis representing K (Kmax = 20 for all cases except the complete dataset of Pachypus, in which we choose Kmax = 50), and the y-axis representing the cross-validation error. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 1357911 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 Pachypus (all) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 12345678910 11 12 13 14 15 16 17 18 19 20 Pachypus A1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 12345678910 11 12 13 14 15 16 17 18 19 20 Pachypus A2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Pachypus B1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 12345678910 11 12 13 14 15 16 17 18 19 20 Pachypus B2 0 0.1 0.2 0.3 0.4 0.5 0.6 12345678910 11 12 13 14 15 16 17 18 19 20 Anopheles 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Chrysis 0 0.1 0.2 0.3 0.4 0.5 0.6 12345678910 11 12 13 14 15 16 17 18 19 20 Darwin's finches 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 12345678910 11 12 13 14 15 16 17 18 19 20 Discoglossus 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Drosophila 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 12345678910 11 12 13 14 15 16 17 18 19 20 Heliconius 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 12345678910 11 12 13 14 15 16 17 18 19 20 Lithobius 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Pleophylla 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 12345678910 11 12 13 14 15 16 17 18 19 20 Pteromalus 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 12345678910 11 12 13 14 15 16 17 18 19 20 Rana 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Sphaerophoria 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Stygopholcus 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 12345678910 11 12 13 14 15 16 17 18 19 20 Taygetis
Figure S23. Plots of IBD testing for cases in which all three H0 hypotheses could be tested (see also Table S5). IBD was rejected in all here shown cases except one (sp4-1 vs. sp4-2+sp4-3). intra-group distances (sp. 4-1) intra-group distances (sp. 4-2+3) inter-group distances intra-group distances (exca1) intra-group distances (exca2) inter-group distances 4.6 4.8 5.0 5.2 5.4 5.6 5.8 0.00 0.05 0.10 0.15 geographical distance intra-group distances (sp. 3) intra-group distances (sp. 4) inter-group distances intra-group distances (sp. 1) intra-group distances (sp. 5) inter-group distances 2.5 3.0 3.5 4.0 4.5 5.0 0.05 0.10 0.15 0.20 0.25 geographical distance sp. 1 vs sp. 5 sp. 3 vs sp. 4 genetic distance exca 1 vs exca 2 sp. 4-1 vs sp. 4-2+4-3 intra-group distances (melo1) intra-group distances (melo3+4) inter-group distances intra-group distances (sp. 2a) intra-group distances (sp. 4) inter-group distances melo1 vs melo 3&4 sp. 2a vs sp. 4 3.4 3.6 3.8 4.0 4.2 0.05 0.10 0.15 0.20 geographical distance 4.6 4.8 5.0 5.2 5.4 5.6 0.00 0.05 0.10 0.15 geographical distance 4.6 4.8 5.0 5.2 5.4 0.01 0.02 0.03 0.04 0.05 0.06 geographical distance genetic distance 4.6 4.8 5.0 5.2 5.4 0.02 0.04 0.06 0.08 geographical distance genetic distance genetic distance