2.2. Reproductive isolation
Mating incompatibility between I. graellsii males and I. elegans and hybrid females is known in populations from the north-west hybrid region (Sánchez-Guillén, Wellenreuther, & Cordero-Rivera, 2011), but not its potential temporal and spatial variation. Thus, the strength of the reproductive barriers involved in preventing copulation in north-central was used to test temporal and/or spatial variation. Similar data from north-west populations are available from an unpublished study (Arce-Valdés, Ballén-Guapacha, Chávez-Rios & Sánchez-Guillén, under review).
For reproductive isolation analysis, last instar larvae of approximately 200 I. elegans and 200 I. graellsii individuals, both from north-central Spain (Mateo, Valbornedo and Villar), were sampled in June-July of 2016-2017. Larvae were transported to the laboratory and maintained until adulthood (for details about larval rearing methodology see Sánchez-Guillén et al., 2005; Van Gossum, Sánchez-Guillén & Cordero-Rivera, 2003). Mechanical isolation to reach the copulation was estimated by measuring incompatibility between the male cerci and female prothorax (males attempt to grasp females) and the incompatibility between male-female genitalia (both genitalia come into contact), in heterospecific crosses of I. elegans with I. graellsiifemales, and I. graellsii males with I. elegansfemales. Additionally, hybrids in the laboratory (from crosses betweenI. elegans and I. graellsii ) were used to produce backcrosses to measure mechanical isolation in both directions (hybrid males with I. elegans females and I. elegans females with hybrid females; and hybrid males with I. graellsii females andI. graellsii males with hybrid females).

2.3 Genome-wide SNP markers

2.3.1 DNA extraction and RAD library preparation

Genomic DNA from head and thorax tissues of 187 individuals from the 20 populations (6-10 individuals per site; Table 1) was extracted with Qiagen DNeasy Blood & Tissue Kit. Extracted DNA was quantified using Nanodrop and visually controlled for DNA degradation using a 1% agarose gel. Five single-digest RAD DNA libraries were processed according to the protocol implemented in Etter et al. (2011) and modified in Dudaniec et al. (2018). All samples, plus five sample replicates, were distributed across the five separately prepared RAD libraries with 40 unique barcodes used per library (sourced from Metabion). Each library was paired-end sequenced (2*100 bp) on a separate lane of an Illumina HiSeq 2500 at SNP&SEQ Technology Platform at Uppsala University, yielding 180 million read pairs per lane (i.e. per library).
2.3.2. Quality checking and SNP calling
Libraries were processed using pipelines within STACKS v2.2 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013; Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011). Raw reads were demultiplexed with process_radtags and PCR clones were identified and discarded with clone_filter using default parameters. Sequence reads were aligned to the I. elegans draft genome assembly (Chauhan et al., 2021) using BOWTIE2 v.2.3 (mismatch allowance per seed alignment of 1, maximum mismatch penalty of 6 and minimum of 2, maximum fragment length of 1000 bp and minimum of 100 bp, Langmead & Salzberg, 2012). We used the ref_map pipeline to detect SNPs using default parameters. Only SNPs with a minor allele frequency of >0.05 and a maximum observed heterozygosity of 0.7 were retained. Moreover, the locus had to occur in 80% of the individuals in a population and in 18 of the 20 populations to be included in the final SNP set. SNP markers were filtered to include only a single random SNP on each RAD tag to create a data set without closely linked loci (using the write_single_snp option in STACKS). Finally, by using the I. elegans reference genome (Chauhan et al., 2021) SNPs were filtered to include only those located on autosomal scaffolds.
Exploratory analyses of population structure revealed possible hybridization in two of the I. graellsii samples from Seyhouse (Algeria); probably with a third Ischnura species (see Supplementary Figure S1). These two samples were removed from further analyses leaving the final total sample size of 185 (Table 1).
2.3.3. Identification of species-specific loci
Diagnostic species-specific loci are powerful markers for assigning later generation hybrids and detecting introgressed alleles in population genetic studies (Hohenlohe, Amish, Catchen, Allendorf, & Luikart, 2011). To provide a list of species-specific markers, alternatively fixed SNPs between the parental species from the allopatric distribution (n=43 I. elegans and n=25 I. graellsii ) were identified using VCFtools v0.1.16 (Danecek et al., 2011). SNPs for each of the two allopatric regions that had only one allele (–max-maf 0) were selected, and then, shared loci between the two allopatric regions were found using the intersect( ) function of R (R Core Team, 2016). Following that, we applied to those loci the Hardy-Weinberg test implemented in VCFtools (–hardy ), and removed loci fixed for the same allele in the two species (HE=0). The remaining 381 SNPs (out of the 5,702 SNPs) were considered as the species-specific markers set.