Quantitative Trait Locus Analysis of Mating Behavior and Male Sex Pheromones in Nasonia Wasps

A major focus in speciation genetics is to identify the chromosomal regions and genes that reduce hybridization and gene flow. We investigated the genetic architecture of mating behavior in the parasitoid wasp species pair Nasonia giraulti and Nasonia oneida that exhibit strong prezygotic isolation. Behavioral analysis showed that N. oneida females had consistently higher latency times, and broke off the mating sequence more often in the mounting stage when confronted with N. giraulti males compared with males of their own species. N. oneida males produce a lower quantity of the long-range male sex pheromone (4R,5S)-5-hydroxy-4-decanolide (RS-HDL). Crosses between the two species yielded hybrid males with various pheromone quantities, and these males were used in mating trials with females of either species to measure female mate discrimination rates. A quantitative trait locus (QTL) analysis involving 475 recombinant hybrid males (F2), 2148 reciprocally backcrossed females (F3), and a linkage map of 52 equally spaced neutral single nucleotide polymorphism (SNP) markers plus SNPs in 40 candidate mating behavior genes revealed four QTL for male pheromone amount, depending on partner species. Our results demonstrate that the RS-HDL pheromone plays a role in the mating system of N. giraulti and N. oneida, but also that additional communication cues are involved in mate choice. No QTL were found for female mate discrimination, which points at a polygenic architecture of female choice with strong environmental influences.

Restrictions to gene flow between species can be divided into factors that act before fertilization, called prezygotic isolation barriers, and those that cause postzygotic isolation after fertilization. Some "speciation genes" causing postzygotic isolation have been identified (Orr 2005;Presgraves and Glor 2010), but less is known about the genetic basis of traits responsible for prezygotic isolation (Arbuthnott 2009). Differences in mating behavior often form primary reproductive barriers, and evolve rapidly through sexual selection in the early stage of the speciation process (Grant and Grant 1997). Indeed, prezygotic isolation seems to evolve at lower levels of overall genetic divergence than postzygotic isolation, at least in sympatry (Coyne and Orr 1989), suggesting that sexual isolation can evolve quickly. Previous studies investigating the genetic architecture of prezygotic isolation barriers in insects have suggested that a polygenic basis is common. For example, in a quantitative trait loci (QTL) study of two closely related cricket species, Laupala paranigra and L. kohalensis, Shaw et al. (2007) found that male calling song differences are due to many genes of small to moderate effect. Gleason and Ritchie (2004) reported six QTL for the differences in male courtship song interpulse interval between Drosophila simulans and D. sechellia. In a review of the genetic basis of female mate preferences and species isolation in Drosophila, Laturney and Moehring (2012) concluded that, although females appear to use the same traits for both within-and between-species mate choice, to some extent a different genetic basis appears to underlie these choices (see also Arbuthnott 2009). In insects, candidate genes for female choice are usually sought among those involved in auditory or olfactory systems or among receptors in the brain that process these signals. Surprisingly, very few studies have considered both male and female mating signals simultaneously, which is required for a full understanding of the genetic basis of prezygotic isolation.
There are many advantages of Nasonia for genetic study of reproductive isolation barriers. One is its haplodiploid reproduction: males are haploid and develop from unfertilized eggs, whereas females are diploid and develop from fertilized eggs. As dominance effects do not exist in haploids, haplodiploidy greatly facilitates quantitative genetic analysis of traits in males, such as genetic linkage mapping and QTL studies (Gadau et al. 1999(Gadau et al. , 2002Koevoets and Beukeboom 2009;Loehlin et al. 2010aLoehlin et al. , 2010bNiehuis et al. 2011;Gadau et al. 2012). Another advantage is the feasibility of interspecific crosses in the laboratory. In nature, the four Nasonia species are reproductively isolated due to infection with species-specific strains of Wolbachia bacteria that cause cytoplasmic incompatibility and hybrid breakdown in interspecific crosses Werren 1990, 1993;Bordenstein and Werren 1998;Bordenstein et al. 2001). Antibiotic (e.g., tetracycline) curing in the laboratory allows interspecific crosses and genetic analysis of speciesspecific traits Werren 1990, 1993;Werren and Loehlin 2009;Sharon et al. 2011). Other advantages of the Nasonia species complex are the availability of full genome sequences of the four species, and high density marker maps . This makes positional cloning of candidate genes identified with QTL studies feasible, as recently demonstrated for a wing size difference (Loehlin et al. 2010a(Loehlin et al. , 2010b, and a pheromone dimorphism (Niehuis et al. 2013).
In the Nasonia species complex, differences in courtship behavior and sex pheromones appear to be responsible for premating isolation between species. All Nasonia species perform a complex mating ritual that consists of a series of interactions between the male and female and ends with female receptivity and copulation (Whiting 1967;van den Assem 1986;van den Assem and Werren 1994;Beukeboom and van den Assem 2002;Velthuis et al. 2005;Burton-Chellew et al. 2007). We previously found 14 QTL for male courtship behavior in interspecific crosses between N. vitripennis and N. longicornis (J. Gadau, C. Pietsch, J. van den Assem, S. Gerritsma, S. Ferber, L. van de Zande and L. W. Beukeboom, unpublished data). Velthuis et al. (2005) reported three major recessive loci for female mate choice between N. vitripennis and N. longicornis. Nasonia males release a long-range sex pheromone to attract virgin females (Ruther et al. 2007(Ruther et al. , 2008, and a different pheromone to induce receptivity in courted females (van den Assem et al. 1980;Ruther et al. 2010). Niehuis et al. (2011) identified genes for alkene biosynthesis with a high similarity to Drosophila in a QTL study of cuticular hydrocarbon differences between N. giraulti and N. vitripennis. Furthermore, Niehuis et al. (2013) showed that N. vitripennis is the only Nasonia species whose males biosynthesize the sex pheromone component (4R,5R)-5-hydroxy-4-decanolide (RR-HDL), besides (4R,5S)-5-hydroxy-4-decanolide (RS-HDL) and 4-methylquinazoline. By genetic mapping and gene knockdown, they narrowed down the genetic basis of enzymes involved in the synthesis of this pheromone component to three closely linked genes on chromosome 1. Finally, following Steiner et al. (2006), Buellesbach et al. (2013) showed that female cuticular hydrocarbon profiles are used by males as cues for interspecific mate discrimination. However, the genetic basis of female discrimination behavior has not yet been investigated.
In this study, we investigate the species-specific quantitative differences of a component of the male sex pheromone (4R,5S)-5-hydroxy-4decanolide and female preference in reciprocal interspecific crosses between N. giraulti and N. oneida. This is the youngest species pair in the Nasonia complex that exhibits strong assortative mating: N. oneida females discriminate strongly against N. giraulti males, but N. giraulti females are less choosy against heterospecific males (Raychoudhury et al. 2010). We use QTL analysis with a SNP marker map including candidate genes to investigate the genetic architecture of quantitative changes in male sex pheromone and female preference.

Experimental design and strains
We set up reciprocal interspecific crosses between N. giraulti and N. oneida using two inbred Wolbachia-cured iso-female strains RV2X(u) and NONY11/36TET. RV2X(u) is descended from the wild type N. giraulti strain RV2, collected in Rochester, New York (Breeuwer and Werren 1995). NONY11/36TET is derived from N. oneida strain NONY11/36, collected in Brewerton, New York (Raychoudhury et al. 2010). These are the same strains used for the Nasonia Genome Project ).
Wasps were cultured in a climate room with constant temperature at 25°, 16:8 h light/dark cycle, and 45% relative humidity. Reciprocal interspecific crosses were set up by mating males of one species with virgin females of the other species. F 1 virgin females were collected and allowed to oviposit to generate recombinant F 2 hybrid haploid males. These males were subsequently mated to females of either of the two pure species to generate a F 3 generation that is referred to as clonal sibships of hybrid F 3 females ( Figure 1). Individuals within a sibship have identical genotypes, as their hybrid fathers are recombinant haploids and produce identical sperm, whereas their mothers are derived from pure species lines that are iso-female and inbred. This experimental set-up allowed for repeated testing of identical genotypes and a quantitative, rather than binary, estimate of mate discrimination (see below).

Behavioral assays
The typical courtship behavior of Nasonia includes several stages (van den Assem 1986; van den Assem and Beukeboom 2004;Clark et al. 2010). The male produces a long-range sex pheromone as described above to attract the female. After a latency period, the male recognizes the cuticular hydrocarbon profile of the female and mounts the female ). Subsequently, he starts his display by touching the antennae of the female with its own, followed by movements of the head, called "head-nods," together with vibration of the wings, followed with a pause. This pattern of head-nods and pauses is termed a cycle, which is repeated several times. At each first head nod in a cycle, the male deposits an as yet unidentified short-range pheromone (aphrodisiac) on the female's antennae. After a number of cycles, the female will either accept the male by lowering her antennae and raising her abdomen for copulation, or reject the male, who will dismount and terminate courtship. Successfully mated males back up on top of the female after copulation to perform postcopulatory courtship consisting of a few similar cycles of head-nods and pauses with a pheromone released, and subsequently dismount.
Observations were performed in a climate room with constant temperature at 25°, a 16:8 h light/dark cycle and 45% relative humidity. Virgin males and virgin females were placed individually in glass tubes (height 30 mm, diameter 10 mm) 1 d prior to the mating trial to allow plenty time for males to release the long-range sex pheromone, as evidenced by white dots marked by males on the surface of the glass tubes (Steiner and Ruther 2009). Males and females were subsequently paired in no-choice experiments by joining the two glass tubes. Couples were then observed under a stereo binocular microscope until first copulation, or for a maximum of 10 min. All females and males were used only once.
The following sequential courtship components were scored: (a) "interest", males and females approached each other; (b) "cross direction", which sex passed the border of the two joint glass tubes first; (c) "border cross time", the time when the first wasp passed the border of the two joint glass tubes; (d) "latency time", the time from the start of the tube joining until the male mounting on the female; (e) "mounting", the male mounted on top of the female; (f) "arrest", the female becomes immobile after the male's mounting; (g) "display", the male performs repeated cycles of head-nods and wing vibrations interrupted with pauses; (h) "attempted copulation", the male attempts to copulate; (i) "normal copulation", copulation occurs. After observation, males were frozen at -20°for chemical and DNA analysis, and females were given two fly pupae as host to produce offspring. For each reciprocal cross, about 250 individuals were divided in two groups of 125 males to be crossed to either parental species females, resulting in the observation of 494 48-to 72-hr-old F 2 hybrid males. In the next generation, female mate discrimination was measured with the same experimental set-up. A total of 2922 24-to 48-hr-old virgin hybrid F 3 females were tested with a virgin pure male of either parental strain. Female mate discrimination was scored as "mate rejection", i.e., when the male mounted the female but the female did not become receptive (Velthuis et al. 2005). In most cases five females were tested per sibship, and the mate discrimination values could take a value between 0 (all females within a sibship accepted their mate) and 1 (all females rejected their mate). Pairs in which the male did not mount the female within the 10 min of observation were discarded, resulting in a total of 2210 hybrid F 3 females, of which 2148 (578 sibships) with genotypic information were used for QTL mapping of mate discrimination.
Chemical analysis N. giraulti males produce a sex pheromone composed of (4R,5S)-5hydroxy-4-decanolide (RS-HDL) and 4-methylquinazoline in an abdominal gland in large amounts to attract females (Niehuis et al. 2013;Ruther et al. 2014). To quantify the amount of RS-HDL, males were killed after behavioral assays by freezing them at -20°. Pheromones were extracted by immersing the male abdomen of each F 2 individual in 80 ml CH 2 CL 2 [containing 88 ng methyl palmitate (1.1 ng/ml) as an internal standard] in glass vials for 2 hr. After pheromone extraction, the abdomen was removed from the solvent, and extracts were stored at -20°. Next, the solvent of the extracts was evaporated up to 1 ml using a gentle stream of nitrogen. The entire residual volume was used for the chemical analysis. Gas chromatography/mass spectrometry (GC/MS) analysis was performed with a HP 6890 gas chromatograph coupled with a HP 5973 Mass Selective Detector (Hewlett Packard, Waldbronn, Germany). The GC (split/splitless-injector in splitless mode for 1 min, injected volume: 1 ml at 250°) was equipped with a DB-5 Fused Silica capillary column (30 m · 0.25 mm ID, df = 0.25 mm, J&W Scientific, Folsom, CA). Helium served as a carrier gas with a constant flow of 1 ml/min. The following temperature program was used: start temperature 60°, temperature increase by 5°per minute up to 300°, and kept at 300°for 10 min. The electron ionization mass spectra (EI-MS) were acquired at an ionization voltage of 70 eV (source temperature: 230°). The software HP Enhanced ChemStation G1701AA Version A.03.00 was used for recording and analysis of chromatograms and mass spectra. Quantification via integrated peak areas and internal standard was performed using the same software.

SNP markers development
Using the Nasonia genome sequences, a total of 57 single nucleotide polymorphisms (SNP) markers were developed with 10 cM intervals covering all five chromosomes. The available N. giraulti .bam file , and N. oneida short reads FASTA file (Human Genome Sequencing Center at Baylor College of Medicine), were assembled on N. vitripennis scaffolds Nvit_1.0 version. In order to estimate the position of the SNPs on the genetic map, the aligned sequences of N. giraulti and N. oneida were compared using CLC Genomics Workbench (CLC bio A/S) to those containing the SNP markers between N. giraulti and N. vitripennis developed by Niehuis et al. (2010). Following the same procedure, an additional 50 SNPs were identified in candidate genes that are known to be associated with mating behavior, courtship song, circadian rhythm, and sex pheromones (Supplemental Material, Table  S1). Identification of candidate genes used information from Nasonia and Drosophila courtship behavior, courtship song, circadian rhythm, pheromone production, and pheromone detection (Gleason et al. 2005;Kankare et al. 2010;Niehuis et al. 2013).
For confirmation, all developed SNPs on chromosome 1 were checked for their amplification and polymorphism in N. giraulti and N. oneida. After primer design with Primer3, DNA was amplified by PCR with annealing temperature of 58°, and Sanger sequenced on an ABI3130xl or ABI3730. Sequences were checked and validated with Chromas (version 2.33, Technelysium Pty Ltd), and aligned with MEGA5 (Tamura et al. 2011). Sequences were blasted to the N. giraulti and N. oneida scaffolds with CLC Genomics Workbench to verify the position of the SNP. After these initial confirmatory tests, chip design was performed by Illumina (San Diego, CA). Only the SNPs with quality scores of at least 0.6 were placed on the chip, consisting of 55 background markers and 41 candidate genes. Reciprocal intraspecific and interspecific crosses between Nasonia giraulti and N. oneida.
The parental interspecific cross [panels (B) and (C)] generates genetically identical F 1 hybrid females that produce unique recombinant F 2 hybrid male offspring of a 50:50 genetic mixture of the two parental species. Backcrossing of F 2 hybrid males to parental strain females yielded sibships of hybrid F 3 females with a 75:25 ratio of either parental species' genome. The haploid hybrid male genotype is indicated with the paternal species on top and the maternal species below. The diploid hybrid F 3 backcrossed female genotype consists of a recombined hybrid set, and a pure species set, resulting in 75% genome of one species, and 25% of the other species on average. The diploid hybrid female genotype is indicated with the paternal species first followed by the maternal species. The square brackets indicate species cytoplasm, which is maternally inherited. Male mating behavior and male pheromone quantity were investigated of individual F 2 hybrid males (gray shading), and female mate discrimination was investigated on clonal sibships of F 3 females (gray shading). G, N. giraulti; O, N. oneida. SNP genotyping DNA was extracted from the frozen heads and thoraxes of individual F 2 hybrid males, and the two parental couples (whole wasps) with a standard high salt chloroform protocol (Maniatis et al. 1982). SNP genotyping of 480 DNA samples was performed in five 96-semi-skirt-plates, including 476 F 2 hybrid males and four parental wasps, using the Gold-enGate Genotyping assay on the Illumina BeadExpress platform of NERC Bio-molecular Analysis Facility (Project NBAF653) at University of Sheffield, UK. The BeadExpress raw data were processed using Illumina GenomeStudio (version 2011.1) to determine the sample genotypes, based on signal intensity ("R") and allele frequency (the ratio of signal intensities for the two possible bases, "Theta"), relative to the cluster position for a given SNP marker.
The most crucial aspect of SNP analysis is to define boundaries between clusters of genotypes. A cluster is a group of SNPs that fall into either one of the homozygous alleles or the heterozygous alleles (Hoffman et al. 2012). Nasonia males, being haploid, provide a great advantage in this respect because they have only two homozygote (hemizygote) classes and therefore typically have clear clusters. Any SNP with only one cluster or low call frequency in a cluster was discarded (four out of 96 SNPs, three background markers and one in a candidate gene). In total, genotype information of 475 hybrid males with 92 well-typed SNP markers was successfully obtained for further QTL analysis (Table S1).

Statistical analysis
Statistical analysis of all phenotypic data was performed with the R statistical software (version 3.1.2). A generalized linear model (glm) was used to identify the effect of different factors in the F 2 hybrid male crosses, including male genotype, female partner species, their interaction, and cross type, on all investigated behavioral traits and pheromone quantity. The glm for F 3 females included female genotype, male partner species, their interaction, and cross type on female mate discrimination. The best statistical model was built up from a full model with all possible factors, followed by removal of nonsignificant explanatory factors. Chi-square (x 2 ) tests were used to compare the likelihood of the different models. Post hoc tests (e.g., Tukey or multiple comparisons) were subsequently performed to confirm the differences between groups. ANOVA was used to estimate the effect of partner, mother, father, grandmother, and grandfather species and their interactions, on all the investigated traits in hybrid crosses. ANOVA was also used to estimate the between-and within-sibship variance in hybrid female mate discrimination as a proxy for broad-sense heritability of this trait. The between-sibship variance is the total phenotypic variance (V P ) and the within-sibship variance is the environmental variance (V E ) as all genotypes within a sibship are identical.

Genetic mapping and QTL analysis
The R/QTL package in R (version 3.1.2) was used to identify genetic regions contributing to the variation in phenotypic quantitative traits . A linkage map was generated from the genotypes of 475 F 2 hybrid males with the set of 92 SNP markers. F 3 female genotypes were inferred from the genotypes of their F 2 recombinant father and their pure species mother (iso-female line). Recombination fractions between all pairs of markers were estimated using the Lander-Green algorithm (Lander and Green 1987) to obtain precise genetic distances between markers of different linkage groups. A likelihood of odds (LOD) score was calculated for each individual at each marker according to Lincoln and Lander (1992) to correct for genotyping errors in map construction. The map did not deviate from the expected marker order and spacing based on the annotated genome. R/QTL uses a hidden Markov model to calculate QTL genotype probabilities and effects given the observed marker data, with allowance for genotyping errors. All QTL scans of F 2 male traits were performed separately for both backcross types. QTL scans of F 3 female mate discrimination were performed separately for either partner species. First, standard interval mapping was applied (using the "scanone" function) with a single-QTL genome scan with a stepsize of 1 cM, with a normal model for all phenotypic traits based on Haley-Knott regression (Haley and Knott 1992). Permutation tests based on 10,000 permutations yielded 5% genome-wide LOD significance thresholds. Next, a two-dimensional genome scan with a stepsize of 2 cM with Haley-Knott regression was performed (using the "scantwo" function), allowing the estimation of additive and epistatic effects by evaluating several types of models, (1) a full model of additive and epistatic effects, (2) an additive model, (3) an epistatic model, (4) a comparison between the full model and the best single-QTL model, and (5) a comparison between the additive model and the best single-QTL model. Finally, a multiple-QTL model was fit with QTL identified from the one-and two-QTL scan (using function "fitqtl"). The fit of the two-QTL model was compared to the reduced model in which a single-QTL is omitted. Possible interactions between the detected QTL and other potential QTL in the multiple-QTL model were investigated (using the "addint" and "addqtl" functions respectively). Next, a fully automated stepwise algorithm (function "stepwiseqtl") optimized the penalized LOD scores (Manichaikul et al. 2009), which yielded the best final QTL model. Separate QTL analysis for partner species yielded partner-specific QTL for some of the traits, suggesting QTL · partner interactions. To quantify this interaction effect, we performed additional QTL scans, which included partner as a covariate and the QTL · partner as an interaction covariate.
Wasp strains are available upon request. File S1 presents ID numbers and genomic locations of SNPs of Table S1. File S2, File S3, and File S4 present phenotypic data of courtship behavior, male pheromone quantity, and female mate discrimination. File S5 contains genotypic data of recombinant F 2 hybrid males.

Data availability
The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

Phenotypic analysis
Courtship behavior in pure species and hybrid males: Females were the first to walk over to the male in 75-88% of the cases, and this did not differ between N. giraulti and N. oneida in pure and hybrid crosses (glm, effect of females: x 2 = 1.20, P = 0.274; effect of males: x 2 = 5.90, P = 0.116; no interaction effect and no cross type effect). Border cross times also did not differ between species in pure and hybrid crosses (glm, effect of females: x 2 = 0.63, P = 0.427; effect of males: x 2 = 6.06, P = 0.109; no interaction effect and no cross type effect).
Latency times did not differ between the pure species crosses, but differed significantly in interspecific crosses and crosses with hybrid males (Figure 2). N. oneida females always had longer latency times than N. giraulti females when they have the same male partner species (glm, effect of females: x 2 = 38.22, P , 0.001; effect of males including the hybrids: x 2 = 12.22, P = 0.007; no significant interaction effect; cross type effect: x 2 = 48.82, P , 0.001). Surprisingly, latency times of N. giraulti females are shorter in crosses with N. oneida males than in crosses with their own species males (mean 6 SE; 129.7 6 18.3 sec vs. 199.9 6 22.6 sec, t-test, t 87.83 = 2.42, P = 0.018). Latency times of N. oneida females are longer in interspecific than intraspecific crosses (mean 6 SE; 284.8 6 27.2 sec vs. 207.5 6 26.5 sec, t-test, t 88.76 = 2.04, P = 0.045). In crosses with hybrid males, latency times are also longer for N. oneida females than N. giraulti females (mean 6 SE; 238.4 6 13.6 sec vs. 155.8 6 9.7 sec, t-test, t 312.15 = 4.99, P , 0.001) albeit with a strong male partner species effect (Table  S3, ANOVA, F 1,367 = 25.56, P , 0.001). These results indicate that the longer latency times of the interspecific crosses are reduced in the hybrid male crosses, consistent with the intermediate genotype of the hybrid males.
Several discrete stages in the Nasonia mating sequence are distinguished. Figure 3A shows the proportion of pairs that reach each subsequent stage for intraspecific and interspecific N. giraulti and N. oneida crosses. The majority of intraspecific crosses results in copulation, and the proportion is not very different for the two species. In interspecific crosses, N. giraulti females have high rates of transitions from one category to the next, resulting in overall higher proportions of completing the mating ritual than N. oneida, with copulation rates of 90.0% (N. oneida male · N. giraulti female) and 62.5% (N. giraulti male · N. oneida female). The main moment of interruption of the mating sequence is at the mounting stage similar to the intraspecific crosses. In contrast to the pure species crosses that are typically not interrupted after mounting, matings of N. oneida females, with N. giraulti males broke up during all subsequent stages. The same holds for matings with hybrid males which in most combinations result in lower final copulation rates than in the pure species crosses (Figure 3, A and B, see also Table S2). There is a strong female partner species effect (Table S3, ANOVA, F 1,490 = 16.80, P , 0.001) and a grandfather species effect (Table S3, ANOVA, F 1,490 = 5.64, P = 0.018) on copulation success in hybrid crosses. The mounting stage is again the most discriminatory (glm, effect of females: x 2 = 15.57, P , 0.001; effect of males: x 2 = 8.59, P = 0.035; no significant interaction effect; cross type effect: x 2 = 27.65, P , 0.001), with a strong female partner species effect in hybrid crosses (Table S3, ANOVA, F 1,490 = 16.05, P , 0.001), but interruptions of the mating sequence also occur during later steps in contrast to the pure species crosses. Overall, N. giraulti females accept pure and hybrid males more often than do N. oneida females.
The number of courtship cycles of a male depends on female acceptance rate and correlates strongly with the total courtship time (correlation coefficient, r 559 = 0.520, P , 0.001). The number of courtship cycles is consistently higher for crosses that do not result in copulation (mean 6 SE; 3.48 6 0.46 vs. 1.89 6 0.06, Mann-Whitney U-test, W = 12,846, P = 0.049). Among intraspecific successful  Mounting is the most discriminatory stage in both types of crosses, but subsequent stages of mating behavior in the hybrid male crosses are also more often terminated. The genotype labeling is as in Figure 1. matings, the number of cycles is higher for N. giraulti than N. oneida (mean 6 SE; 2.35 6 0.10 vs. 1.02 6 0.02, Mann-Whitney U-test, W = 1750, P , 0.001, Table 1). Interspecific matings with N. oneida females require more cycles than pure crosses (mean 6 SE; 1.65 6 0.18 vs. 1.02 6 0.02, Mann-Whitney U-test, W = 1076.5, P = 0.001), but interspecific matings with N. giraulti females take fewer cycles than pure crosses (mean 6 SE; 1.74 6 0.10 vs. 2.35 6 0.10, Mann-Whitney U-test, W = 634, P , 0.001). In hybrid crosses, a strong female partner species effect (Table S3,    Hybrid male pheromone quantity: N. giraulti males produce the long-range sex pheromone RS-HDL in larger quantities than N. oneida males ( Figure 4A, mean 6 SE; 12.35 6 4.80 ng vs. 1.32 6 0.54 ng, Mann-Whitney U-test, W = 1072, P = 0.001). Hybrid males have on average intermediate pheromone levels with no difference between the two reciprocal crosses (mean 6 SE; giraulti-oneida hybrid 11.00 6 1.20 ng, oneida-giraulti hybrids 11.65 6 1.47 ng, Mann-Whitney U-test, W = 32,602, P = 0.098), but the variation in pheromone quantities is much larger due to many transgressive phenotypes ( Figure 4A). There is, however, an interaction effect between the two types of hybrid males and female partner species on pheromone quantity of hybrid males (glm, x 2 = 7.81, P = 0.005), suggesting that pheromone quantity affects female response.
Hybrid female mate discrimination: Mate discrimination is observed in interspecific crosses between pure N. oneida females and pure N. giraulti males (81.6% acceptance), as well as in interspecific crosses between hybrid F 3 females whose genomic compositions were 75% of one species and 25% of the other species and pure N. giraulti males (65.3% and 60.5% acceptance) (Table 2 and Figure 6). Figure 7 shows the mating steps of hybrid F 3 females (see also Table S2). The mating patterns differ according to male partner species (glm, effect of females: x 2 = 3.63, P = 0.304; effect of males: x 2 = 140.97, P , 0.001; no significant interaction effect; cross type effect: x 2 = 150.55, P , 0.001). Matings with N. giraulti males result in fewer successful copulations (708 out of 1474) than matings with N. oneida males (1012 out of 1448). The most prominent difference with pure species crosses (Figure 3A) is that mating disruption does not occur in the mounting step, but during the display and initial courtship ("interest") stages. Analysis of variance for mate discrimination among F 3 females showed a strong partner species effect (Table S3, ANOVA, F 1,589 = 164.69, P , 0.001). Crosses with a N. giraulti male partner had significantly higher between sibship variances than crosses with a N. oneida male partner, which suggests that among-sibship genetic variation underlying female mate discrimination is expressed mainly in the presence of N. giraulti males but not in the presence of N. oneida males (Table 3 and Figure S1).

QTL analysis
Linkage map: A linkage map corresponding to the five Nasonia chromosomes was constructed from the 92 SNP markers in the N. giraulti-N. oneida hybrid cross ( Figure S2). None of the SNPs showed deviation from Mendelian segregation.
Male pheromone: Three QTL were found at a 5% genome-wide level of significance for male pheromone quantity. One QTL was present on chromosome 4 in both the crosses with N. giraulti and N. oneida females, and this may be the same QTL as their positions overlap ( Figure 4B). Another QTL was present on chromosome 1 in crosses with N. giraulti females only. The two-QTL scan split the QTL on chromosome 1 with N. giraulti partners into two QTL. The multiple-QTL scan yielded one additional QTL on chromosome 5 in crosses with N. oneida females. Total variances explained out of the final model are 10.2% for the QTL on chromosomes 4 and 5 in crosses with N. oneida females, and 14.9% for the three QTL in crosses with N. giraulti females (Table 4 and Figure 4B). For the QTL on chromosomes 4 and n 5, hybrid males with the "G" allele of the associated marker (C4M8, C4M10, and C5M6) have higher pheromone levels than males with the "O" allele in both cross types, consistent with the larger pheromone quantities in pure N. giraulti males. The QTL on chromosome 1 in crosses with N. giraulti females had opposed effects, i.e., higher pher-omone quantity in crosses with N. oneida females (Table 4). There were significant additive effects (LOD = 7.21, P = 0.001), but no epistatic effects (LOD = 0.90, P = 0.997).
Male courtship behavior: Because of the strong female partner species effect on male courtship, the male courtship QTL analysis was performed separately for backcross type. In the F 2 hybrid male crosses with either pure species females, the QTL scans for cross direction yielded one QTL on chromosome 4 for crosses with N. oneida females, explaining 6.1% of the phenotypic variance, but none for N. giraulti as partners (Table 4). No significant QTL were detected at 5% significance level for border cross time, latency, mounting rate and number of cycles as part of the courtship, even after inclusion of the significant covariate partner species. If we take a 20% genome-wide significance level to further obtain suggestive QTL, peaks are visible for latency on chromosome 1 at 108 cM (LOD = 1.32, P = 0.122); for mounting rate on chromosome 1 at 116 cM (LOD = 0.97, P = 0.171); and for number of cycles, on chromosome 4 at 35.9 cM (LOD = 1.08, P = 0.139), and on chromosome 5 at 52 cM (LOD = 1.13, P = 0.126) (data not shown). There were two significant QTL at the 5% level for copulation success in crosses with N. giraulti females, one on chromosome 1, and one on chromosome 3, explaining 4.8% and 2.3% of the phenotypic variance, but none with N. oneida females. Hybrid males with the "G" allele of marker C1M2 have lower copulation success than males with the "O" allele, but the effect is opposite for marker C3M6 (Table 4).
Female mate discrimination: The two reciprocal types of F 2 hybrid males were backcrossed with either a N. giraulti or a N. oneida female, yielding a total of 2148 hybrid F 3 females, for which mate discrimination was scored. Because female mate discrimination appeared to be expressed only in the presence of N. giraulti males, but not in the presence of N. oneida males, we performed QTL analysis separately for the two species of males ( Figure 6B). As expected, we found suggestive QTL for female mate discrimination (chromosome 3, at QTL mapping results for female mate discrimination. Strong mate discrimination occurs in interspecific crosses between pure N. oneida females and pure N. giraulti males, and in interspecific crosses between hybrid F 3 females whose genomic compositions were 75% of one species and 25% of the other species, and pure N. giraulti males but not pure N. giraulti females. Number of males accepted and number of pairs observed are shown in Table 2. The dashed line in (B) shows the 5% genome-wide significance level from permutation tests of the single-QTL genome scan. Upper and lower panel show results for females mated to N. oneida and N. giraulti males, respectively. Sample sizes are shown in the left upper corner, and species names of the male partner in the right upper corner.

Figure 7
Mating behavior progress in hybrid females. Proportions reaching subsequent stage in the mating behavior process are shown for different hybrid F 3 females. Failure of female arrest and male display occur more often in crosses with N. giraulti (dashed lines) than with N. oneida males (solid lines). The genotype labeling is as in Figure 1.
52.2 cM, LOD = 2.08, P = 0.090) when F 3 females were mated with N. giraulti males (explaining 3.3% of the phenotypic variance), but no QTL when mated with N. oneida males. The relatively weak QTL effects detected for this trait are consistent with the fact that broad-sense heritability for this trait is only about 0.2. The inclusion of partner as both additive covariate and interactive covariate makes little difference in the QTL analysis ( Figure S3).
Candidate genes: Forty candidate genes for mating behavior were used as markers in our QTL analysis (Table 4 and Figure 4B). Four candidate genes on chromosome 4, protein-1-like, disco, lateNAy, and Ato, were associated with the QTL for cross direction in crosses with N. oneida females, and for pheromone quantity in both crosses with N. giraulti and N. oneida females. All candidate genes on chromosome 1 (e, fru, fix_nod, XP001602953, per, dNA, cycle, lateNAy, csp, MpK2(2), and Rhodophilin (B)) and chromosome 3 (TipE, nonA, dco, beethoven, headnod, and acyl-CoA) were associated with QTL for copulation success in crosses with N. giraulti females as this QTL spans the complete chromosome. Three candidate genes on chromosome 5, slo, Fmr1, and Dy(2), were associated with the QTL for pheromone quantity in crosses with N. oneida females. Candidate genes, e, fru, csp, Mpk2(2), and Rhodophilin(B), on chromosome 1 and an additional one, So (1), on chromosome 4 were associated with QTL for pheromone quantity in crosses with N. giraulti females.

DISCUSSION
The goal of this study was to investigate the genetic architecture of male courtship traits and female preference in two recently diverged Nasonia species, N. giraulti and N. oneida, that exhibit strong prezygotic isolation. We performed a QTL analysis of male pheromone quantity, male courtship behavior, and female choice by crossing the two species. The differences in mate discrimination between the two species were generally confirmed in our behavioral analysis. Raychoudhury et al. (2010) showed that N. oneida females discriminate strongly against N. giraulti males, whereas N. giraulti females are less discriminatory. We also found a stronger mate discrimination of N. oneida females in interspecific crosses, although N. giraulti males were not accepted at a rate as low as 20% reported by Raychoudhury et al. (2010). Females were found to typically approach the males in our single pair crosses. N. oneida females had longer latency times than N. giraulti females, and these were even longer in interspecific crosses with N. giraulti males.
There was a strong partner species effect in latency, mounting rate, number of courtship cycles, and copulation success in the hybrid male crosses. Stronger discrimination against heterospecific males was further confirmed by more frequent breaking up of mating behavior stages in interspecific crosses with N. oneida females, particularly in the mounting stage, resulting in lower overall copulation rate.
An important factor in the mating process is the long-range male sex pheromone component (4R,5S)-5-hydroxy-4-decanolide (RS-HDL), which is used to attract females from a distance (Ruther et al. 2007). It is produced at an almost 10-fold higher level by N. giraulti than N. oneida males. Although this pheromone functions to attract females, it apparently also affects female behavior at later stages of the mating process, such as their willingness to accept the male for copulation. We used hybrid crosses to obtain males with different pheromone quantities and tested their mating behavior in trials with pure species females. Interestingly, hybrid male pheromone quantities were sometimes higher than pure N. giraulti males, suggesting that pheromone production is partly disrupted (transgressive phenotypes) in hybrid individuals. Hybrid males with more pheromone quantity had significantly higher mating success with pure females of either species. Hybrid females were found to reject N. giraulti males more often than N. oneida males, but surprisingly this did not depend on their relative proportions of either species genotype. Interestingly, male rejections appeared to happen mostly during the display stage after mounting in N. giraulti but not in N. oneida, which suggest that the communication between the two partners is partly disrupted. As the long-range sex pheromone RS-HDL is considered to play a role before the mounting stage to attract females from a distance, a different trait may be responsible for the breaking up of the display stage. Good candidates are the cuticular hydrocarbon pattern , or the short-range volatile pheromone that males produce in their mandibular glands, and deposit on the female's antennae during male head-nods (van den Assem et al. 1980). Another possibility is male courtship song, which differs between species (W. Diao, L. van de Zande, M. G. Ritchie, and L. W. Beukeboom, unpublished data).
The role of the RS-HDL pheromone in the mating process of the various Nasonia species deserves further attention. N. oneida males apparently produce the RS-HDL component at much lower quantities than the other Nasonia species. This does not seem to negatively affect their acceptance by N. giraulti females, who accepted N. oneida males at even higher rates than their own males. In contrast, N. oneida females discriminate strongly against N. giraulti males. This suggests that the RS-HDL pheromone of N. giraulti males functions as a repellent to N. oneida females. A similar effect was found for the sex pheromone component (4R,5R)-5-hydroxy-4-decanolide (RR-HDL), which is produced only by N. vitripennis males, and which makes them less attractive to N. longicornis and N. giraulti females (Niehuis et al. 2013;Ruther et al. 2014). This explanation is consistent with rapid shifts in sender and receiver cues in establishing prezygotic isolation in this genus Niehuis et al. 2013). However, it leaves the question open as to how mate attraction occurs within the N. oneida species without the long-range male sex pheromone.
n Table 3 Between-and within-sibship variance of mate discrimination in hybrid females used in the QTL analysis An important assumption in our pheromone analysis is that the amount of pheromone measured in the male abdomen after mating is representative of the amount of pheromone released before the courtship ritual to attract the female. Ruther et al. (2007) showed that the pheromone titer in the male abdomen depends on the age of the male. Amounts of HDL are close to zero in freshly emerged males, increase within the first 2 d after emergence, and remain at a constant level on d 3. We consistently used virgin males of 2-3 d old to control for this age effect on pheromone synthesis. It is not known how fast pheromone release depletes the amount of pheromone present in the abdominal gland when a female partner is around. As males can attract many females in a row, and a positive correlation between abdominal pheromone level and mate acceptance was found in both Nasonia species, our assumption of a correlation between quantity of the long-range male sex pheromone released and measured in the abdomen is likely valid.
QTL studies on hybridizing species have been informative about the underlying genetic architecture of interspecific mate discrimination.
They have revealed a few QTL with large effect in some organisms though many QTL with small effects is more typically found (Coyne et al. 1994;True et al. 1997;Ritchie and Philips 1998;Doi et al. 2001;Takahashi et al. 2001;Shaw and Parsons 2002;Henry et al. 2002;Gleason and Ritchie 2004;Arbuthnott 2009). The vast majority of examples concern Drosophila species, in which male signal traits are typically components of the courtship song. Many courtship song QTL have been identified and in some cases the underlying gene has been implicated (Gleason and Ritchie 2004;Etges et al. 2007;Schäfer et al. 2010;Ellison et al. 2011;Lagisz et al. 2012;Chung et al. 2014). Much less is known about the genetic basis of female receptivity, but auditory and/or olfactory receptors are likely candidates (Laturney and Moehring 2012). Due to a lack of studies covering a wider array of insect species, a consistent pattern has not yet emerged about the genetic architecture of courtship differences between species. A complicating factor is that the genetic basis is often different for intraspecific and interspecific matings (e.g., Gleason and Ritchie 2004;Arbuthnott 2009). Our study found two QTL for copulation success, and four for n Significant QTL with positions, effect and explained genetic variance were identified with multiple-QTL regression models. a SNP markers corresponding to candidate genes listed in Table S1.
pheromone quantity on chromosomes 1, 3, 4, and 5 in the Nasonia oneida-giraulti species pair. Although this suggests that there may be some genomic regions with moderate effect involved in these traits, the resolution of our QTL analysis may have been limited, and QTL of minor effect may have gone undetected. Pheromone synthesis depends on many enzymatic reactions (Chung and Carroll 2015), and may therefore be considered as a complex trait. Niehuis et al. (2013) detected two genomic regions, one on chromosome 1 and one on chromosome 4, involved in the production of a pheromone component that is specific to N. vitripennis. Whether these coincide with our observed QTL on chromosomes 1 and 4 requires a higher resolution QTL study. In the study of cuticular hydrocarbon differences between Nasonia species, Niehuis et al. (2011) identified at least 42 QTL. They also listed several candidate genes for cuticular hydrocarbon production. One of these genes is protein-1-like (LOC100118619), which corresponds to our C4M8 marker, and that is centrally positioned under the pheromone quantity QTL on chromosome 4. As pheromone quantity was correlated with mating success, this gene is a candidate gene for mate discrimination in Nasonia. There are also a number of other candidate genes on chromosome 4, including So(1), disco, late-NAy, and Ato, that were used as markers in our QTL study, and that correspond to the QTL peaks. These genes were selected for their known function in mating behavior, courtship song, circadian rhythm, or pheromone production and detection (Gleason et al. 2005;Kankare et al. 2010;Niehuis et al. 2013; J. Gadau, C. Pietsch, J. van den Assem, S. Gerritsma, S. Ferber, L. van de Zande, and L. W. Beukeboom, unpublished data). More detailed functional analyses, such as knock down studies, are needed to establish the possible role of these genes in RS-HDL quantity and interspecific mate discrimination.
Despite the large size of the mapping population (578 sibships), no significant QTL were detected for female preference at the 5% significant level, although some peaks were visible at the 20% genome-wide significance level. The broad-sense heritability of this trait is 0.2, so our results are probably consistent with a polygenic basis of female interspecific discrimination and strong environmental effects. However, the result is in contrast with two other studies. Velthuis et al. (2005) found a few major QTL for heterospecific male acceptance in crosses between N. giraulti and N. longicornis. We found major QTL on chromosomes 1, 2, 3, and 4 in interspecific crosses of N. giraulti and N. oneida when confronted with N. vitripennis males (M. C. W. G. Giesbers, B. A. Pannebakker, L. van de Zande, and L. W. Beukeboom, unpublished data). N. giraulti and N. oneida females discriminate strongly against N. vitripennis males, and this is mediated by the additional pheromone component, (4R,5R)-5-hydroxy-4-decanolide (RR-HDL) (Niehuis et al. 2013). Taken together, these results demonstrate that the genetic architecture for mate discrimination in the Nasonia species complex consists of loci with major effects, and loci with minor effects, and differs to some degree based on the species pair considered.
There are several possible explanations for the lack of female mate discrimination QTL in our study. Our experimental design did not allow for directly estimating heritability, as females were not mated with similar males. Comparison of the between-sibship variance and withinsibship variance revealed a V P of 0.2-0.3 in crosses with N. giraulti male partners, consistent with values obtained from a response to selection study (M. C. W. G. Giesbers, B. A. Pannebakker, L. van de Zande, and L. W. Beukeboom, umpublished data). However, no substantial V G was found in crosses with N. oneida male partners, which suggests that there is little genetic variation for mate discrimination in interspecific crosses with N. oneida. A possible explanation for our inability to detect mate discrimination QTL is that mate discrimination is highly polygenic and plastic, i.e., with many small effect loci that went undetected because of insufficient power, and whose effect is strongly modulated by environmental conditions. Although our experimental conditions were kept as constant as possible (including temperature, humidity, time of day, and wasp age), there may have been uncontrolled factors, such as developmental differences due to host quality that affect Nasonia mating behavior.
Hybridization studies can be very informative about the speciation process. One major unresolved question concerns the kind of genetic changes that cause reproductive isolation. The Nasonia species complex is very suitable for genetic studies of speciation as different species can be intercrossed in the laboratory and confronted with heterospecific and hybrid (recombinant) mating partners. We have shown that male courtship behaviors and pheromone quantities differ between two recently diverged species and have a genetic basis. We also demonstrated that female interspecific mate discrimination in Nasonia is partly governed by these male traits, albeit in a complex way probably involving many genes of small effect and strong environmental effects. It is evident that the communication between the two sexes relies on multiple different cues, and that these cues partly differ in intraspecific and interspecific mating interactions. Further progress can be made with more detailed genomic studies and functional knockdowns of candidate genes in this system.