RADseq processing, loci discovery and variant calling
DNA was extracted from blood stored in SET buffer with salt extraction protocol (Aljanabi & Martinez, 1997). Quality controls were performed with QUBIT to ensure sufficient starting material. Restriction-associated DNA tags were obtained with double-digesting the DNA of 388 owls with PstI/BamHI restriction enzymes following a customized version of Petersen et al (2012). Library construction, sequencing, demultiplexing were performed by Bioname©, Finland. Quality checks were performed with FASTQC (Andrews, 2014). We utilized the new assembly as a reference for mapping the RADseq reads. Mapping was performed with bwa-mem2 with a minimum mapping quality of 20 incorporating both single and paired-end reads and discarding unmapped reads (Li, 2013). Coordinates-alignment files were fed into Stacks v2.4 pipeline for variant calling (Rochette, Rivera‐Colón, & Catchen, 2019). We constructed 6 different SNP panels in order to explore the impact of filtering and subsequent captured loci in exploring associations with colour phenotypes. Specifically, we tested filters on minor allelic frequency (MAF = 0.01; MAF = 0.05), loci presence (R) within the population (no filter; R = 0.5; R =0.8), while holding a constant maximum heterozygosity of 0.50 and always collecting one random SNP per RADtag. Due to the expected association between colour and genes involved in the melanin-production pathway mentioned above, we tracked the number of reads mapped to those regions across the different filtering steps. We accounted for RADseq reads mapped within contigs where those genes where found, within the respective coding sequence, as well as in a range of +/- 20Kb flanking the coding sequence. We further calculated the percentage (in base pairs) of coding sequences covered by our double-digest RAD-sequencing strategy.