Recurrent selection shapes the genomic landscape of differentiation between a pair of host-specialized haplodiploids that diverged with gene flow

Understanding the genetics of adaptation and speciation is critical for a complete picture of how biodiversity is generated and maintained. Heterogeneous genomic differentiation between diverging taxa is commonly documented, with genomic regions of high differentiation interpreted as resulting from differential gene flow, linked selection, and reduced recombination rates. Disentangling the roles of each of these non-exclusive processes in shaping genome-wide patterns of divergence is challenging but will enhance our knowledge of the repeatability of genomic landscapes across taxa. Here, we combine whole-genome resequencing and genome feature data to investigate the processes shaping the genomic landscape of differentiation for a sister-species pair of haplodiploid pine sawflies, Neodiprion lecontei and Neodiprion pinetum. We find genome-wide correlations between genome features and summary statistics are consistent with pervasive linked selection, with patterns of diversity and divergence more consistently predicted by exon density and recombination rate than the neutral mutation rate (approximated by dS). We also find that both global and local patterns of FST, dXY, and π provide strong support for recurrent selection as the primary selective process shaping variation across pine sawfly genomes, with some contribution from balancing selection and lineage-specific linked selection. Because inheritance patterns for haplodiploid genomes are analogous to those of sex chromosomes, we hypothesize that haplodiploids may be especially prone to recurrent selection, even if gene flow occurred throughout divergence. Overall, our study helps fill an important taxonomic gap in the genomic landscape literature and contributes to our understanding of the processes that shape genome-wide patterns of genetic variation. Methods We collected Neodiprion lecontei and N. pinetum mid- to late-instar larval colonies from Lexington, Kentucky and surrounding areas. To confirm sex (and therefore ploidy), we reared the larvae to adults in the lab using standard lab protocols. The adults were either preserved in 100% ethanol (stored at -20 °C) or flash frozen and stored at -80 °C. To avoid sampling close relatives, we selected one individual from each larval colony (each colony typically represents a group of siblings). In total, we sampled 20 N. lecontei females and 18 N. pinetum females. We extracted DNA from head and thorax tissue with a Qiagen DNeasy Blood & Tissue Kit. We followed the standard manufacturer protocol for insects, including an optional RNase A step. DNA concentration was measured with a Quant-iT dsDNA High-Sensitivity fluorescence assay (Invitrogen). A single library was prepared using a Tn5 tagmentation protocol following Bendall (2020). Whole-genome resequencing was performed with 150bp paired-end sequencing technology in an Illumina Hi-Seq X sequencer at Admera Health (Plainfield, NJ). The library was first run in a single lane of the sequencer. A subset of samples that had low read counts were then re-pooled and run on a second lane. Demultiplexed reads were cleaned using trimmomatic v0.39 with the following criteria: (a) remove adapters, (b) perform sliding window trimming where a sequence is cut when a window (4bp) drops below a quality threshold (15), (c) remove low quality bases (quality score < 3) from the beginning and end of each read. These trimmed reads are published in the NCBI SRA database (BioProject accession number PRJNA1107580). Then, for each lane separately, the cleaned reads were mapped to the high-quality chromosome-level N. lecontei reference genome (iyNeoLeco1.1, GCA_021901455.1) using the BWA-MEM algorithm in bwa v0.7.17. We used samtools v1.13 to mark PCR duplicates (‘markdup’) and remove ambiguously mapped reads/secondary alignments (‘view -F 1284 -f 0x02’). We then used samtools to merge the filtered reads from the two lanes (‘merge’) and index the resulting bam files (‘index’). All genetic summary statistics were calculated in ANGSD v0.933. Recombination rates (cM/Mb) were obtained from a previously published high density linkage map generated from a cross between divergent N. lecontei populations (Linnen et al. 2018; Herrig et al. 2024). We used the Neodiprion lecontei genome annotation (GCF_021901455.1_iyNeoLeco1.1_genomic.gtf; Herrig et al. 2024) to calculate exon density in 50kbp windows. To estimate mutation rate, we calculated the synonymous substitution rate (dS) in 50kbp windows. To estimate distance from the centromere, we first used Juicebox v1.11.08 to visualize the HiC contacts for each chromosome and estimated the midpoint of each centromere by identifying the local maximum delta in the number of contacts between adjacent loci within each chromosome. We then calculated the distance of each window from the centromere by taking the absolute value of the midpoint of each window subtracted from the midpoint of the centromere.