Accurate and testable species delimitation hypotheses are essential for measuring, surveying and managing biodiversity. Today, taxonomists often rely on mitochondrial DNA barcoding to complement morphological species delimitations. Although COI barcoding has largely proven successful in assisting identifications for most animal taxa, there are nevertheless numerous cases where mitochondrial barcodes do not necessarily reflect the species history. For instance, what is regarded as one single species can be associated with two distinct DNA barcodes, which can point either to cryptic diversity or to deep within-species mitochondrial divergences with no reproductive isolation. In contrast, two or more species can share barcodes, for instance due to mitochondrial introgression. These intrinsic limitations of mitochondrial DNA barcoding can only be addressed with nuclear genomic markers, which are expensive, labour intensive, poorly repeatable, and often require high-quality DNA. To overcome these limitations, we examined the use of ultraconserved nuclear genetic elements (UCEs) as a quick and robust genomic approach to address such problematic cases of species delimitation. This genomic method was assessed using six different bee species complexes suspected to harbour cryptic diversity, mitochondrial introgression, or mitochondrial paraphyly. The sequencing of UCEs recovered between 686 and 1860 homologous nuclear loci and provided explicit species delimitation hypotheses in all investigated species complexes. These results provide strong evidence for the suitability of UCEs as a fast method for species delimitation even in recently diverged lineages. Furthermore, this study provided the first conclusive evidence for both mitochondrial introgression among distinct species, and mitochondrial paraphyly within a single bee species.

UCE library preparation

Whole body DNA extractions were performed overnight in a proteinase K buffer at 56°C and purified using a Qiagen Biosprint 96 extraction robot following the manufacturer’s protocol. Extracts were quantified using Qubit v3 (Thermofisher Scientific) and 50 ng DNA per specimen were sonicated to 500 bp fragment length using a Bioruptor ultrasonicator (Diagenode). Two independent dual-indexed libraries each containing 96 specimens were constructed using a Kapa Hyper prep kit (Roche) using one fourth of the manufacturer’s recommended volumes (as described in Branstetter, Longino, Ward, & Faircloth, 2017). PCR amplifications were performed in the recommended volumes. PCR products were quantified using a Qubit v3 and each row of a 96-well PCR plate were pooled equimolarly (i.e. for total of 8 pools). Libraries were UCE enriched using the Hymenopteran v2 hybridization kit (UCE Hymenoptera 2.5Kv2 Principal/Full, myBaits, Arborbiosci). Each enrichment was performed on a single pool of 12 specimens using 500 ng. The enrichment protocol followed the manufacturer’s recommendations with a hybridization step of 24 h at 65°C, followed by a PCR amplification with 14 cycles. Pools were sequenced on a Miseq using the Illumina v3 kits (2 x 300 bp; Illumina).

Bioinformatic processing of UCE data

Fastq reads were demultiplexed on the Miseq and data from all runs were merged and processed mainly using PHYLUCE tools (Faircloth, 2016). Raw data were cleaned with illumprocessor (Faircloth, 2016), a tool wrapped around trimmomatic (Bolger, Lohse, & Usadel, 2014). Clean reads were assembled with SPAdes v3.12.0 (Nurk et al., 2013) using the single-cell flag (“--sc”), careful option (“--careful”) and a coverage cutoff value of five (“--cov-cutoff”). Obtained contigs were mapped against the corresponding UCE reference file using Lastz (Harris, 2007) and matching reads were extracted and aligned by species complex using MAFFT (Katoh & Standley, 2013). Alignments were edge-trimmed using the PHYLUCE “seq-cap” program; a strategy recommended for closely related species (< 30-50 MYA) (Faircloth, 2016). Loci shared by less than 75% of the maximum number of specimens were filtered out. Remaining alignments were concatenated and saved in fasta format. An additional filtering step was applied to remove specimens with more than 90% missing data.