RAD-Seq of B. terrestris
Thoracic tissue from 90 haploid males (including 44 wildtype and 46 yellow mutant individuals), were extracted using a Qiagen DNeasy Blood and Tissue Kit with RNaseA treatment. Samples were extracted into a 250 l  extraction buffer, dried using a SpeedVac, and re-eluted in 50 l  water. Extracted samples were quantified using the Qubit BR dsRNA kit to ensure at least 250 ng of DNA, and were assessed on a 1% agarose gel to ensure high-quality genomic DNA. Samples were prepared for RAD-Seq using protocols following ref.77,78. This involved digesting 250 ng of DNA per sample in a 50 l reaction with PstI-HF (NEB) for 1 hour at 37C, followed by 20 min. at 80C. Samples were divided into six sets (libraries) of 16, in each of which 16 different P1 barcoded adaptors were ligated using T4 DNA ligase with incubation for 22 C for 1 hour, 65C for 10 mi., and 30 min at room temperature. Each library was pooled and sheared to 300-700 bp using a Covaris S2 sonicator. Gel-based size selection of 300-700 bp fragments was then performed and samples purified using a Nucleospin Gel and PCR clean-up kit. Samples were end-repaired with the Quick Blunting Kit (NEB) and purified with Agencourt AMPure XP (Beckman Coulter) magnetic beads. This was followed by addition of dATP overhangs added with the Klenow exo (NEB), purification with AMpure XP beads, ligation of P2 adaptors, and additional purification with AMpure XP beads. Samples were then PCR amplified (98 C for 30 sec, 16 cycles of 98C for 10 sec and 72C for 1 min, 72C for 5 min) using Phusion High-Fidelity Master Mix with specific P2 Sanger indexing primers for each of the six pools. This process enabled the unique barcoding of all 9 samples. RAD-tag short read (2*150 bp, paired-end) libraries were sequenced using an Illumina HiSeq 2500 in Pennsylvania State University Genomics Core Facility (University Park, Pennsylvania, USA). RAD-Seq analysis was completed using Stacks 2 software79 following the recommendations of a published RAD-Seq analysis protocol80. First, the raw paired-end sequencing reads for six paired-end libraries were filtered and demultiplexed to generate individual-sample level sequence dataset using process_radtags unit of Stacks 2 software. Reads from individual samples (n=90, 44 wildtype and 46 yellow individuals) were aligned to the published B. terrestris genome assembly