Faster‐haplodiploid evolution under divergence‐with‐gene‐flow: Simulations and empirical data from pine‐feeding hymenopterans

Although haplodiploidy is widespread in nature, the evolutionary consequences of this mode of reproduction are not well characterized. Here, we examine how genome-wide hemizygosity and a lack of recombination in haploid males affects genomic differentiation in populations that diverge via natural selection while experiencing gene flow. First, we simulated diploid and haplodiploid “genomes” (500-kb loci) evolving under an isolation-with-migration model with mutation, drift, selection, migration, and recombination; and examined differentiation at neutral sites both tightly and loosely linked to a divergently selected site. So long as there is divergent selection and migration, sex-limited hemizygosity and recombination cause elevated differentiation (i.e., produce a “faster-haplodiploid effect”) in haplodiploid populations relative to otherwise equivalent diploid populations, for both recessive and codominant mutations. Second, we used genome-wide SNP data to model divergence history and describe patterns of genomic differentiation between sympatric populations of Neodiprion lecontei and N. pinetum, a pair of pine sawfly species (order: Hymenoptera; family: Diprionidae) that are specialized on different pine hosts. These analyses support a history of continuous gene exchange throughout divergence and reveal a pattern of heterogeneous genomic differentiation that is consistent with divergent selection on many unlinked loci. Third, using simulations of haplodiploid and diploid populations evolving according to the estimated divergence history of N. lecontei and N. pinetum, we found that divergent selection would lead to higher differentiation in haplodiploids. Based on these results, we hypothesize that haplodiploids undergo divergence-with-gene-flow and sympatric speciation more readily than diploids.

Methods
DNA was extracted and ddRAD libraries were prepared for 23 Neodiprion pinetum larvae and 44 N. lecontei larvae collected from multiple locations in Kentucky, as well as 4 interspecific hybrids and an additional 18 N. lecontei samples from an allopatric population in Michigan. These libraries were sequenced using 150-bp paired end reads on an Illumina HiSeq 4000. We also collected whole-genome resequncing data from  a N. virginiana sample to be used as an outgroup in demographic analyses. Neodiprion sequencing reads are available via the NCBI SRA, accession numbers: SAMN23893940-SAMN23893944, SAMN23893948, SAMN23893960-SAMN23893963, SAMN23893965, and SAMN25157024-SAMN25157101.

We aligned demultiplexed ddRAD reads to the N. lecontei reference genome (Nlec1.1 GenBank assembly accession number- GCA_001263575.2) using the very sensitive setting in bowtie2. We only retained reads that aligned to one locus in the reference genome and had a Phred score greater than 30. For the ddRAD dataset, we removed PCR duplicates using a custom script. We called SNPs in samtools. We required all sites to have a minimum of 7x coverage and 50% missing data or less. We also removed SNPs with significantly more heterozygotes than expected under Hardy-Weinberg equilibrium (an indicator of genotyping/mapping error). We removed any individual that was missing more than 70% of the data. We performed all filtering in VCFtools v0.1.13.

We created several datasets with subsets of individuals and additional filtering for each of the population genetic analyses. We generated three data sets with minor allele filtering (MAF, SNPs <0.01 removed): 1) sympatric N. pinetum and N. lecontei for genome-wide patterns of divergence (36,935 SNPs), 2) sympatric N. pinetum, N. lecontei, and hybrids for admixture analysis (35,649 SNPs), and 3) sympatric N. pinetum, N. lecontei, allopatric N. lecontei, and outgroup N. virginiana for ABBA- BABA tests (12,905 SNPs). We also generated a down-sampled dataset (described below) without a MAF filter for estimating site-frequency spectra (SFS) that included sympatric N. pinetum, N. lecontei, and N. virginiana for demographic analyses.