Skip to main content

Genomics of sex allocation in the parasitoid wasp Nasonia vitripennis



Whilst adaptive facultative sex allocation has been widely studied at the phenotypic level across a broad range of organisms, we still know remarkably little about its genetic architecture. Here, we explore the genome-wide basis of sex ratio variation in the parasitoid wasp Nasonia vitripennis, perhaps the best studied organism in terms of sex allocation, and well known for its response to local mate competition.


We performed a genome-wide association study (GWAS) for single foundress sex ratios using iso-female lines derived from the recently developed outbred N. vitripennis laboratory strain HVRx. The iso-female lines capture a sample of the genetic variation in HVRx and we present them as the first iteration of the Nasonia vitripennis Genome Reference Panel (NVGRP 1.0). This panel provides an assessment of the standing genetic variation for sex ratio in the study population. Using the NVGRP, we discovered a cluster of 18 linked SNPs, encompassing 9 annotated loci associated with sex ratio variation. Furthermore, we found evidence that sex ratio has a shared genetic basis with clutch size on three different chromosomes.


Our approach provides a thorough description of the quantitative genetic basis of sex ratio variation in Nasonia at the genome level and reveals a number of inter-related candidate loci underlying sex allocation regulation.


The study of sex allocation is one of the most successful areas in evolutionary biology [1,2,3]. Theoretical predictions of optimal resource allocation to male and female offspring in response to environmental conditions are now supported by a wealth of empirical data [1, 4,5,6,7,8,9] . This is particularly true for local mate competition theory [10, 11]. Briefly, local mate competition theory describes sex allocation dynamics when related males (including full siblings) compete for mates in locally structured populations. When kin competition for mates is high, termed local mate competition by Hamilton (1967), female-biased sex ratios, that limit kin competition for mates and maximise the number of available mates for those competing males, are favoured [12]. Local mate competition is maximal when the male offspring from a single female compete for mates (their sisters) on a discrete mating patch. When local mate competition is reduced, with increasing numbers of unrelated males competing for mates, less female-biased sex ratios are favoured. The logic of Hamilton’s local mate competition theory extends to situations where females experience environments that vary in the expected level of local mate competition amongst sons, for instance when different numbers of females oviposit together on a patch. Facultative sex allocation is then predicted, and has been shown in a diverse range of organisms, from malaria parasites to fig wasps [3, 13]. Perhaps the most extensive exploration of local mate competition theory has occurred in parasitoid wasps, however, especially in Nasonia vitripennis [14, 15].

Despite the success of sex allocation theory at the phenotypic level, we still know rather little about the underlying mechanisms of sex allocation at the molecular level [16,17,18,19]. This is important for at least three reasons. First, deviations in sex ratio from theoretical predictions, showing variation around the predicted optimal sex ratio for a given number of foundresses, still remain [20] and require explanation. In Nasonia vitripennis, extending and testing basic local mate competition theory has revealed that understanding how female parasitoids obtain and use information relating to the expected level of local mate competition their sons will experience can refine theoretical models and explain much of this deviation from the predicted optimal value [15, 21,22,23,24,25]. However, independent from local mate competition, there is also genetic variation in sex allocation in N. vitripennis [19, 26,27,28,29,30], which can result from the genetic architecture constraining adaptation [31]. Second, adaptive sex allocation provides an opportunity to explore the mechanisms underlying adaptation, including the genetic architecture and molecular evolution of a well-characterised trait (for recent discussions see for example [32,33,34,35,36,37];). For instance, to what extent does genetic architecture constrain sex allocation? Third, phenotypic plasticity is at the heart of facultative sex allocation, and there has long been interest in how such plasticity is encoded in the genome [38]. Altogether, sex allocation is a well-characterised plastic trait, offering the opportunity to dissect the genetic basis of plasticity.

To date, studies on the genetic basis of sex ratio have mainly focussed on quantitative genetic analysis across a range of organisms, including Drosophila [39, 40], ants [41], snails [42], fish [43], turtles [44], birds [45], pigs [46] and humans [47]. Several species of parasitoid wasp also show genetic variation for sex ratio [48,49,50,51,52], but most work has focused on Nasonia vitripennis [26,27,28, 30, 53] and studies of the molecular mechanisms underlying sex allocation are still scarce. Recently, mutation accumulation studies have explored how mutation augments additive genetic variation in sex ratio [16], suggesting that genes governing sex ratio are pleiotropic, also influencing other fitness-related traits. A linkage mapping study of sex ratio variation within a natural population of N. vitripennis identified a major quantitative trait locus (QTL) on chromosome 2 and three weaker QTL, one on chromosome 3 and two on chromosome 5. In addition, the QTL data revealed scope for pleiotropy between offspring sex ratio and clutch size [19]. Finally, recent work has explored patterns of gene expression associated with sex allocation in Nasonia. Thus far, patterns of differential gene expression associated with oviposition behaviour have been revealed in females, but there is no suggestion that the facultative response to local mate competition cues made during oviposition is associated with any change in gene expression [17, 18, 54]. This suggests that, whilst oviposition is associated with changes in gene expression, phenotypically plastic sex allocation is not.

Here, in order to better understand the genomics of facultative sex allocation and oviposition, we extend our characterisation of the genetic basis of sex ratio variation in N. vitripennis with a genome-wide association study (GWAS) for single foundress sex ratios using iso-female lines derived from the recently developed outbred laboratory strain HVRx [55]. The iso-female lines capture a sample of the genetic variation in HVRx and we present them as the first iteration of the Nasonia vitripennis Genome Reference Panel (NVGRP 1.0). This panel provides an assessment of the standing genetic variation for sex ratio in the study population. Our approach provides the first genomic visualization of the quantitative genetic basis of sex ratio variation in Nasonia and reveals a number of candidate molecular processes underlying sex allocation.


Molecular variation in NVGRP

The NVGRP was constructed by isolating mated females from the HVRx outbred laboratory population, followed by 9 generations of inbreeding. 34 NVGRP lines were sequenced using Illumina HiSeq 2000, and the reads were mapped to the N. vitripennis reference genome. The mean sequence depth was 5.4X per line, with on average 230.2 Mb (96.47%) of the effective reference sequence covered per line (Supplementary Table 1). SNP calling was performed using JGIL, which uses the relatedness between the iso-female lines in the NVGRP to simultaneously call SNPs based on the aggregate coverage across lines, resulting in a decreased genotyping error rate, even at low coverages per individual iso-female line [56]. By further taking into account the expected allele frequencies and error probabilities after 9 generations of inbreeding from the outbred HVRx laboratory population, JGIL called 205,691 SNPs based on the aggregated coverage across all the NVGRP lines (mean aggregated coverage per SNP 207.88, s.e. = 1.41). The residual heterozygosity in the NVGRP lines was low, on average the proportion of segregating SNPs across lines was 2.68% (s.e. = 0.35%), which is less than is theoretically expected after 9 generations of inbreeding (13% given an inbreeding coefficient of F = 0.87). The majority of the SNPs in the NVGRP (96,920) are found in intergenic regions. 76,142 SNPs occur in the introns, 6250 SNPs in the 5′ UTR and 7086 in the 3′ UTR sequences. In the coding regions, there are 12,344 synonymous SNPs and 6723 nonsynonymous SNPs (6634 missense and 89 nonsense SNPs), resulting in a dN/dS ratio of 0.545.

NVGRP population genetics

We calculated the genome-wide polymorphism (the number of segregating sites (Watterson’s θ) and nucleotide diversity (π)), both per chromosome, and in 400 kb non-overlapping windows (Fig. 1, Supplementary Table 2). Averaged over the entire genome, π = 1.268 *10− 3 (s.e. = 2.504*10− 5), and θ = 1.331*10− 3 (s.e. = 2.330*10− 5). The average polymorphism differs per chromosome for both π (F4,539 = 3.102, P = 0.015) and θ (F4,539 = 3.332, P = 0.010). Nucleotide diversity is higher towards the tips of the chromosomes (Fig. 1). Comparison of the NVGRP to the HVRx source population shows lower levels of both the nucleotide diversity (π) and the number of segregating sites (Watterson’s θ) among the 34 NVGRP lines. The difference between the NVGRP and the HVRx is largest for the number of segregating sites, which suggests loss of low frequency alleles, and is further reflected in the higher genome-wide Tajima’s D values in the NVGRP (Fig. 1). Linkage disequilibrium (LD) decays to half the of the initial value in 17.8 kb (half-decay distance, Fig. 2), with an average pairwise r2 value of 0.08 (s.e. = 0.0005). The extend of LD varies across and within the chromosomes. Some regions on chromosomes 1 and 3 exhibit half-decay distances over 8.9 Mb and 36 Mb respectively (Fig. 3). Pairwise FST values over all variable sites between all the NVGRP lines were on average 0.309 (s.d. = 0.046) with a minimum of 0.145 and a maximum of 0.465 (Supplementary Figure 1). This indicates that inbreeding on the one hand has resulted in highly divergent lines, but, on the other hand, ample overlap in genetic variation still exists between all pairs of lines (i.e., we did not find lines completely divergent for all SNPs). Furthermore, the pairwise FST values did not correlate with the absolute differences in sex ratio between lines selected for the GWAS (Spearman’s rs = − 0.034, P = 0.542) indicating that no correction for population structure in the GWAS is needed.

Fig. 1
figure 1

Mean nucleotide diversity π (upper panel), mean number of segregating sites Watterson’s θ (middle panel) and Tajima’s D (lower panel) over non-overlapping 400 kb windows across the chromosomes in NVGRP (red and blue lines) and the HVRx laboratory outbred population (grey and black lines)

Fig. 2
figure 2

Decay of linkage disequilibrium with physical distance. Dots shows the r2 among pairs of SNPs, red solid line gives the non-linear least squares fit of r2 on the distance between pairs of SNP. Dashed line indicates the half-decay LD distance at 17.8 kb

Fig. 3
figure 3

Linkage disequilibrium half-decay distance across the genome. For each SNP on each of the five chromosomes, the log10(LD half-decay distance) is plotted

Sex ratio GWAS

Nasonia vitripennis is naturally infected with cytoplasmic incompatibility inducing Wolbachia [57]. To avoid possible confounding effects of Wolbachia infection status on the sex ratio, we collected sex ratio data from NVGRP lines with the same infection status. The majority of the 26 NVGRP tested positive, while 8 lines had lost their Wolbachia infection during 36 generations of culturing after initial inbreeding (Supplementary Figure 2). Among the Wolbachia-positive lines, line 62 was a clear outlier, showing very male-biased sex ratios and low clutch sizes, indicative of a transient Wolbachia infection status (Supplementary Figure 3) resulting in within-line cytoplasmic incompatibility. We therefore excluded line 62 from further analysis. Overall, sex ratio data from 25 NVGRP lines was used in our study. Across these 25 NVGRP lines, there was highly significant among-line variation in sex ratio (Binomial GLM: F24,918 = 15.244, P < 0.0001). This among-line variation represents a broad-sense heritability of H2 = 0.106 for sex ratio (Likelihood Ratio = 56.176, P < 0.0001). For clutch size, the 25 NVGRP lines showed a broad-sense heritability of H2 = 0.078 (Likelihood Ratio = 36.61, P < 0.0001). Among the NVGRP lines, our data showed a weak, non-significant negative correlation between sex ratio and clutch size (Supplementary Figure 4, linear regression: b = − 0.0015 (s.e. = 0.0008), t = − 1.823, P = 0.08).

In total, 18 SNPs were significantly associated with sex ratio, according to our empirical false discovery rate (FDR) threshold of 0.1. (Fig. 4, Table 1, Supplementary Table 3). These SNPs represent a single peak on chromosome 1 and show a high degree of linkage (mean r2 = 0.97, Supplementary Table 4). All SNPs represent common variants, with an average minimal allele frequency of 0.416. These 18 SNPs were associated with 9 annotated loci (Table 1). While the majority of the significant SNPs are predicted to either be intron variants (8) or upstream variants (7) of genes, we found two significant SNPs with a potential effect on the expression of the associated loci. NC_015867.2_30,570,243 is causing a synonymous change (p.Ala224Ala) in the splice region of exon 4 of cilia- and flagella-associated protein 52 (CFAP52, LOC100124028), while NC_015867.2_30,604,770 (P = 5.00 *10− 7) is a variant in the 5′ untranslated region (UTR) of ADP-ribosylation factor-like protein 4C (Arl4C, LOC100124030). We will discuss both these SNPs in more detail below.

Fig. 4
figure 4

Manhattan plot for offspring sex ratio (top panel) and brood size (bottom panel) in N. vitripennis in the GWAS experiment, showing –log10(P)-values of the single marker regressions for every polymorphic SNP across the chromosomes of the Nasonia vitripennis Genetic Reference Panel (NVGRP). Red line indicates the empirical threshold at a q-value of 0.1, corresponding to P = 2.00*10–6, or –log10(P) = 5.70 for sex ratio. Green highlighted SNPs show the 400 kb windows in which the P-values for sex ratio and brood size overlap more than expected by chance

Table 1 Description of the SNPs significantly associated to sex allocation in the NVGRP. Table shows the SNP ID, the location (chromosome and position), the reference and alternative nucleotide, the minor allele frequency (MAF), the gene (NCBI Nasonia vitripennis Annotation Release 102 Gene ID, the description) for which SNPEff predicted an effect its variant definition according to the Sequence Ontology terminology and the -log10(P) value of the single marker regression for that SNP

We found no SNPs significantly associated with clutch size variation across the lines. However, despite this absence, we did detect significant enrichment of overlap between the P-values of the sex ratio and clutch size GWAS in 400 kb windows on chromosomes 2, 3 and 5 (Supplementary Figure 5). Most notably, two regions on chromosome 5, between 4000 kb – 4800 kb and between 8400 kb – 8800 kb show enriched overlap in test statistics, indicating a shared genetic background for these two traits (Fig. 4). Although our finding of 18 strongly linked variants in the same genomic region suggests that the number of variants in the NVGRP segregating for sex ratio is limited, this overlapping background of sex ratio with clutch size on three different chromosomes, in addition to the a large number of sub-significant peaks for sex ratio throughout the genome, suggests that sex allocation is a complex trait with polygenic regulation and epistatic interactions.


Sex allocation has been extensively studied at the phenotypic level, especially in species such as parasitoid wasps. However, our understanding of the genetic architecture of the mechanisms associated with sex allocation is much more rudimentary. Here we have shown that at least 18 single nucleotide polymorphisms are associated with variation in sex allocation across experimental lines drawn from a single population. While only a single peak significantly passed the genome-wide threshold, both the presence of a large number of subsignificant peaks (Fig. 4), and evidence of an overlapping genetic background on chromosome 5, confirm the polygenic nature of the genetic variation in sex allocation in Nasonia. While the emerging picture is that variation in most complex traits is typically influenced by a large number of loci with pleiotropic alleles, polygenic traits represent a major challenge in terms of understanding the molecular genetic basis of phenotypes and their adaptive evolution [37]. This is the case in terms of both understanding how interactions among DNA sequences constrain or promote the production of adaptive phenotypes, and also in terms of reconstructing the path evolution has taken at the molecular level, including identifying key substitutions. Put simply, we are likely to have an embarrassment of riches, although even then, we may well only be able to discover a fraction of evolutionarily relevant SNPs [33]. Moreover, unlike candidate gene approaches, we are often faced with genes or putative regulatory sequences about which we know little or nothing, either through a lack of annotation, or through the lack of a clear link between proposed gene ontology and the phenotype of interest. In our case, the 18 SNPs on chromosome 1 that are associated with sex ratio variation among our lines have predicted effects on 9 annotated loci.

Of these, two significant SNPs are predicted to have an effect on the expression of the loci they are associated with. NC_015867.2_30,570,243 (P = 1.12*10− 6) is a variant causing a synonymous change (p.Ala224Ala) in the splice region of exon 4 of cilia- and flagella-associated protein 52 (CFAP52, LOC100124028). CFAP52 is highly conserved and its human ortholog, WDR16, occurs in sperm tails [58]. Furthermore, human low-motility sperm shows low expression of this testis specific transcript WDR16 [59]. While no direct studies have been done into the role of CFAP52 in Nasonia sperm, CFAP52 shows high expression levels in N. vitripennis testis [60, 61]. Splice region changes could affect the functionality of CFAP52, and because of its role in sperm motility, affect male fertility, as was found for splice variants in other cilia- and flagella-associated proteins in humans (CFAP43/44 and CFAP70 [62, 63]). Such an effect could directly impact sex allocation in the haplodiploid N. vitripennis, in which fertilization is the key trigger to female development [64]. Results such as these would help explain the small influence of males in sex allocation observed in N. vitripennis [65]. The observed variation among strains in sex ratio could well be due to differences in male ejaculate quality and fertilization success. Such variation in male fertility might also explain why the observed single foundress sex ratios are slightly higher (i.e. more males) than we might expect from the number of females a single male can inseminate, and would act as a form of ‘fertility insurance’ [3, 66]. Whether the splice region variant CFAP52 observed in the NVGRP can account for such variation requires a more thorough characterization of the potential mis-splicing and is outside the scope of the present study [67].

The other SNP that may well influence gene expression is NC_015867.2_30,604,770 (P = 5.00*10− 7). This is a variant in the 5′ untranslated region (UTR) of ADP-ribosylation factor-like protein 4C (Arl4C, LOC100124030), but does not generate a premature start codon. The 5′ UTR is the mRNA region directly upstream from the start codon and plays a role in the regulation of transcript translation. While it is hard to predict the exact effect of a 5′ UTR variant in a non-model system, 5′ UTR variants can affect the translation efficiency and hence the expression level of proteins [68]. Arl4 are group of proteins that belongs to the family of ADP-ribosylation factors, which are involved in the regulation of vesicular trafficking processes [69]. Arl4C has been implicated in the regulation of vesicular traffic of cholesterol and of transferrin in endosomes [70] and likely plays a similar role as the related Arl4A [71]. Interestingly, Arl4A is strongly expressed in adult testis in mice [72] and targeted disruption of Arl4A resulted in a reduced sperm count [73]. Functional studies are needed to determine if Arl4C also plays a role in Nasonia spermatogenesis, and whether the identified 5′ UTR variant affects Arl4C expression and alters the Nasonia sex ratio phenotype, again via effects on sperm quality or quantity.

While no significant SNPs were found for clutch size, we did detect significant enrichment of overlap between the P-values of the sex ratio and clutch size GWAS in 400 kb windows on chromosomes 2, 3 and 5, with those on chromosome 5 providing the strongest evidence (i.e. P-values were non-randomly associated with each other for the two traits across these regions, whereas there should be no association if there is no shared genetic basis: Fig. 4, Supplementary Figure 5). This overlap indicates a shared genetic basis for the two traits, compatible with theoretical predictions on the genetic basis of sex ratio. Based on the observed natural variation and estimates of mutational parameters for sex ratio, the genetic variation for sex ratio is predicted to be maintained by selection on pleiotropic loci with effects on other fitness related traits [16]. While clutch size is only one of many fitness-related traits that can show pleiotropy with sex ratio, a correlation between these two traits is also expected from theoretical work. In the case of a single foundress per host, local mate competition theory predicts that a mother should only produce enough sons to mate with all of her daughters [10]. For small clutch sizes, this can be a single male, and when clutch size increases, an increasingly male-biased sex ratio is predicted [20, 74, 75]. Interestingly, our previous analysis of the quantitative genetic basis of sex ratio and clutch size in different Dutch N. vitripennis strains, also identified overlapping QTLs for both traits on chromosomes 2 and 5 on a recombination linkage map [19]. Unfortunately, of the eight microsatellite markers associated to the sex ratio in the QTL study, only two could be found on the same chromosomes in the current assembly (Nvit_2.1). The other six are located on unlocalized scaffolds (data not shown), making a direct comparison between these studies difficult and indicating the need for a further improved genome assembly for this species [76].

In addition to revealing some of the molecular genetic variants associated with sex allocation in Nasonia, we have also presented a new genomic resource for the community. Since the publication of the Nasonia genome in 2010 [77], there has been a steadily growing number of studies of the molecular genetics of a range of phenotypes, deploying a number of techniques [76, 78,79,80,81,82,83,84]. In terms of our own work on sex allocation, we have shown for example that facultative sex allocation under local mate competition is not associated with changes in gene expression [17, 18], even though female oviposition of eggs is associated with major changes to the transcriptome (in particular, a down-regulation of metabolic processes: [17, 54]). We have also shown though that disrupting patterns of DNA methylation changes the pattern of facultative sex allocation [85, 86], suggesting that the regulation of gene expression via DNA methylation is important for facultative sex allocation.

Our first iteration of a proposed Nasonia vitripennis Genome Reference Panel (NVGRP 1.0) currently consists of 34 iso-female lines generated from wasps collected from one population, of which 25 lines comprised the GWAS for sex ratio presented here. Whilst the number of lines is currently modest, we have nonetheless captured SNPs associated with sex ratio; moreover, the NVGRP lines exhibit significant broad-sense heritabilities for a range of traits, not just sex ratio, including total lifetime fecundity, longevity, head width, wing length, and starvation resistance (H2 = 0.13–0.58; B.A. Pannebakker, unpublished observations). Importantly, as a species characterised by sib-mating, and thus inbreeding, the creation of iso-female lines in Nasonia does not suffer the problems associated with the inbreeding of natural out-crossers [87], and so we do not see great reductions in fitness in our inbred lines. The effect of sib-mating is reflected in the relatively slow decay of linkage disequilibrium (r2 < 0.1 at 160.2 kb, and r2 < 0.2 at 71.2 kb) as is also observed in selfing plants (e.g. rice Oriza sativa, r2 < 0.1 at 75-150 kb [88] or soy Glycine max r2 < 0.1 at 90-500 kb [89]). Other insects, such as Drosophila melanogaster and Apis mellifera, show much shorter distances at which LD decays (r2 < 0.2 at 10 bp and r2 < 0.2 at 1 kb respectively, [90, 91]).

The extent of LD in the NVGRP lines does show large variation within and across chromosomes, which likely reflects the observed variation in recombination rates in Nasonia [92,93,94]. In direct observations of recombination events from markers segregating in a cross, areas of low recombination were observed near the centre of the chromosomes [92, 93]. The NVGRP shows areas of high LD, but not just at chromosome centres. The LD estimates in the NVGRP, however, are population estimates, integrating historical recombination events that could result in a different pattern than from the present-day recombination estimates observed in the progeny from a cross [95]. To what extent natural inbreeding, or other demographic factors such as the limited population size of NVGRP, is driving the slow decay of linkage disequilibrium and the observed variation in LD in Nasonia, requires further genotyping of additional samples from wild populations.


We present the first iteration of the Nasonia vitripennis Genetic Reference Panel as a community resource for the analysis of complex traits. We found substantial variation for sex allocation in the NVGRP and identified 18 SNPs associated with variation in sex allocation and found evidence for overlapping genetic background of sex ratio with clutch size on different chromosomes. As our data represent only a sample of the standing genetic variation in our study population, it is likely that we have missed other variants that influence sex ratio variation segregating in that population, but we have nonetheless provided the first genomic visualization of the heritability of sex ratio observed in this and other studies on Nasonia [16, 26, 28].

Wild N. vitripennis populations show rather limited population genetic structure across Europe, which perhaps explains our ability to capture significant genetic variation with only a comparatively small sample of lines from one population [96, 97]. Nevertheless, we recognise that the study presented here is as much a proof-of-principle for exploring the molecular genetic basis of sex allocation in Nasonia as anything else, and an expansion of the NVGRP is certainly required. However, we note that as whole-genome sequencing becomes ever cheaper and alternative genotyping methods are developed (such as RAD sequencing and other genotype-by-sequencing techniques [98, 99]), the role of reference panels may change. Shifting from the main focus of a study, they can provide supporting genomic resources, or genetic hypotheses to test, for studies that interrogate genetic variants in the wild more directly [100,101,102,103].


Study organism

Nasonia vitripennis (Hymenoptera, Chalcidoidea) are generalist parasitoid wasps that oviposit in large dipteran pupae, which includes species of Calliphoridae. Females oviposit between 20 and 50 eggs in an individual host, depending on the host species. Development takes approximately 14 days at 25 °C, after which male offspring emerges just before female offspring [104]. Males have short wings and cannot actively fly and remain close to the emergence site where they compete with other males for emerging females, including their own sisters. After mating, only the females disperse to find new hosts.

As with all Hymenoptera, N. vitripennis is haplodiploid, with diploid females developing from fertilised eggs, and haploid males developing from unfertilised eggs. Sex allocation is therefore associated with females controlling the fertilisation of eggs during oviposition, releasing stored sperm from the spermatheca to fertilise eggs as they pass down the oviduct. Unless otherwise specified, wasps were reared on Calliphora spp. hosts at 25 °C, 16 L:8D light conditions.

Base population

We used the HVRx outbred population of N. vitripennis [55] as the base population for our selective breeding experiment and as the source population for our iso-female inbred lines for our genome-wide phenotypic association study. This line was created from wild caught wasps collected from Hoge Veluwe National Park in the Netherlands and is maintained as large outbred population at an effective population size of Ne > 200.

The Nasonia vitripennis Genome Reference Panel 1.0.

NVGRP lines

Iso-female lines were established by randomly collecting 105 virgin females from the HVRx outbred population at HVRx generation 45 [55], which were individually provided with two hosts for 2 days at 25 °C to produce male offspring. Mothers were stored at 4 °C and crossed back to one of her sons (mother-son mating) upon emergence of the male offspring. Mated females were provided with two hosts for 2 days at 25 °C to produce male and female offspring. These lines were inbred by 8 generations of full-sib mating to produce 34 stable inbred lines, followed by random mating. This resulted in an inbreeding coefficient of F = 0.87, equivalent of 10 generations of diploid full-sib matings.

Genome resequencing

DNA was isolated from 60 pooled females at generation 29 for each of the 34 stable inbred lines, using a standard high salt–chloroform protocol [105]. Library construction and genome sequencing were performed by BGI Tech Solutions according to standard Illumina protocols. For each inbred line, a 91 bp paired-end library was constructed and the libraries were run on an Illumina HiSeq 2000 (Illumina, San Diego, CA, USA). Clean reads were aligned to the Nasonia vitripennis genome build Nvit_2.1 (GCF_000002325.3, downloaded from the NCBI FTP site: using BWA v0.6.2-r126 [106]. Duplicated reads were filtered out using samtools v0.1.18 [107] and alignments were generated as BAM files for further analysis. Assembled data have been submitted to the NCBI Short Read Archive in BioProject no. PRJNA387118. SNPs were called for each line simultaneously using the Joint Genotyper for Inbred Lines (JGIL) v1.6 [56] using the following parameters: read mapping quality threshold: 10, number of generations: 10. SNPs with a quality less than 20 were ignored, and the genotypes of any individual for which the SNP quality score was less than 10 were treated as missing genotypes. SNP data were stored as VCF and have been deposited in the EMBL European Variation Archive in Project no. PRJEB33514 with analysis accession no. ERZ1029220.

Population genomics

For the estimation of nucleotide diversity (π), the number of segregating sites (Watterson’s θ), and Tajima’s D, BAM files resulting from the above described analysis were merged (samtools merge) separately for the HVRx outbred population bam files and the inbred line bam files to produce one outbred BAM file and one inbred bam file. BAM files were combined without a priori coverage correction. From these BAM files, pileup files were produced (samtools mpileup) and subsampled to standardize the coverage at 20X, based on the mean coverage depth of the HVRx outbred line, using of PoPoolation [108] and the following settings target-coverage 20, max-coverage 400, min-qual 20, method withreplace. These subsampled pileup files, were used to estimate the nucleotide diversity (π), the number of segregating sites (Watterson’s θ), and Tajima’s D in sliding windows using of PoPoolation [108] and the following settings pool-size 60, min-count 5, min-coverage 10, min-covered-fraction 0.5, fastq-type sanger, window-size 400,000, step-size 200,000.

Linkage disequilibrium (r2) was estimated between all pairs of SNPs within a distance of 1 Mb using the VCF resulting from the JGIL SNP calling pipeline of the inbred lines as input. To estimate LD decay, we fitted the equation r2 = 1/(1 + px) for every focal SNP (using nls method, [109] where x denotes the distance to every other SNP 1 Mb up- and 1 Mb down-stream of the focal SNP. After the value for p was retrieved for every SNP, the distance at which LD equal 0.5 was estimated by (1/p). Then these values were plotted against the position of the SNPs to get an indication of LD variation in the genome. To determine the need for correction based on population structure, we determined the significance of the correlation between pairwise FST values over all SNPs called in the JGIL pipeline and pairwise absolute Euclidian distance in sex ratio for the 26 lines included in the GWAS (see below) using Spearman’s rank correlation. When correlated with the phenotype, population structure can result in false positives in GWAS analyses, requiring a correction based on the relatedness of the samples [110].

Wolbachia detection

Nasonia vitripennis can be infected with the maternally inherited bacterium Wolbachia pipentis, which affects reproduction. To be able to account for the effects of Wolbachia, the NGRP lines were assessed for Wolbachia infection in generation 31 after initial inbreeding using a PCR assay for the Wolbachia specific wsp gene (81F/691R primers [111];. PCR conditions were denaturation at 94 °C for 3 min, then 35 cycles of 94 °C for 1 min, 55 °C for 1 min, 72 °C for 1 min and a final extension at 72 °C for 5 min. PCR products were visualized on a 1% agarose gel stained with ethidium bromide.

Iso-female line GWAS

The focal females used in this experiment were from 26 Wolbachia-positive NVGRP lines. At the point of phenotyping associated with this experiment, these lines had been mass-reared for 36 generations after initial inbreeding.

Hosts can have a large effect on the observed variation in N. vitripennis [112]. To control for such host effects and possible maternal effects, we isolated 40 females (2 day old, mated) from the mass culture of each of the 26 lines into individual glass vials and provided each of the females with three hosts. We then used 40 females per line from the resulting F1 generation in the experiment, resulting in one female per “grandmother”. These experimental females (2 day old, mated) were isolated into individual glass vials and given with a single host for 24 h as a pre-treatment to facilitate egg development, following [22]. Pre-treatment hosts were removed, and each female was given honey solution applied on a piece of filter paper for an additional 24 h. Subsequently, we gave pretreated experimental females access to three hosts for a period of 24 h. One-way escape tubes were fitted to the glass vials after 1-h had passed to allow females to disperse, preventing unnatural levels of superparasitism [113]. Hosts were incubated and we counted and sexed the emerging offspring to calculate sex ratio. In total, we phenotyped 57,465 wasps in 958 broods (mean clutch size = 60.04).

We analysed the variation in sex ratio in two steps. First, we fitted a Generalised Linear Model with binomial error structure and a logit link function to test for significant differences in sex ratio between the iso-female lines. To correct for over-dispersion, common when analysing binomial data, F-tests were used. Second, we determined the between iso-female line variation using linear mixed-effect models on arcsine square root transformed data for sex ratio, and untransformed data for clutch size. Iso-female line was fitted as a random effect and variance components were estimated by REML. This isofemale line analysis estimates the genetic variation as the broad-sense heritability H2 [52, 114,115,116]. All statistical analyses were carried out in R [109], the linear mixed effect models were implemented using the lme function in the nlme package [117].

Genotype-phenotype associations

To perform genotype-phenotype associations, we filtered our NVGRP SNP dataset. All heterozygous sites within NVGRP lines were treated as missing data. We first removed the 9 NVGRP lines that had lost their complete Wolbachia infection and removed all remaining SNPs that were fixed for one allele. Finally, only SNPs that had data for 18 individuals or more were retained. Our final SNP dataset for the analysis consisted of 94,656 SNPs with a minimal allele frequency threshold of 0.04, in which each SNP had on average 1.06 (s.e. = 0,004) missing genotype.

To identify SNPs significantly associated with single-foundress sex ratio, we performed single marker regressions for all SNPs on each chromosome segregating between iso-female lines using a general linear mixed effect model. SNPs segregating within lines were treated as missing data. A Generalized Linear Mixed Model was implemented using the glmer function in the lme4 package [118] in R [109], with sex ratio as the response variable, genotype as the explanatory variable, and line as a random effect. We used a binomial error structure and a logit link function. P-values for the significance of the association between single-foundress sex ratio and each SNP were obtained. Significance thresholds were estimated by permutating iso-female line identity for each SNP. The permutated dataset was then used to estimate an empirical threshold using a q-value of 0.1 as significant, which resulted in an empirical threshold of P = 2.00*10− 6, or –log10(P) = 5.70 for sex ratio.

Similar single marker regressions using a linear mixed model were performed for clutch size, using the lmer function in the lme4 package [118] in R [109], with clutch size as the response variable, genotype as the explanatory variable, and line as a random effect. To determine significance thresholds for clutch size, we permutated iso-female line identity for each SNP. Because the minimum q-value was 0.463, no variant was considered significant for clutch size.

To determine the scope for a shared genetic background, we determined the overlap in test statistics between sex ratio and brood size in 400 kb windows per chromosome. To do this, we calculated the overlap in SNPs for an increasing rank in p-values. We then performed he SuperExact test in the SuperExactTest package [119] in R on this overlap for each cutoff obtained observed p-values. Because genome structure can have an effect on the result, mainly through variation in the number of SNPs in each window, we also performed the same analyses on 100 randomized sites where SNP identity was permutated. Windows were considered significant if the observed p-value was lower than any of the 100 permutated ones. We used an increasing `windows size of 25, 50, 100, 200 and 400 kb windows (Supplementary Figure 5).

Identification of candidate genes

To identify candidate genes, we annotated the VCF file with SnpEff v4.3b [120] based on the NCBI Nasonia vitripennis Annotation Release 102 ( Candidate SNP positions and their effects were filtered from the annotated VCFs using SnpSift v4.3b [120].

Availability of data and materials

Raw DNA sequencing reads deposited in the National Center for Biotechnology Information (NCBI) with the project accession of PRJNA387118. SNP data were stored as VCF and have been deposited in the European Molecular Biology Laboratory (EMBL) European Variation Archive (EVA) in Project no. PRJEB33514 with analysis accession no. ERZ1029220. The NVGRP iso-female lines are available upon request from the laboratory of Bart Pannebakker.

For our analysis, we used the Nasonia vitripennis genome build Nvit_2.1 (GCF_000002325.3, downloaded from the NCBI FTP site:, and the NCBI Nasonia vitripennis Annotation Release 102 (available from:





ADP-ribosylation factor-like protein 4


ADP-ribosylation factor-like protein 4A


ADP-ribosylation factor-like protein 4C


Binary Alignment Map




Cilia- and flagella-associated protein 43/44


Cilia- and flagella-associated protein 52


Cilia- and flagella-associated protein 70


Deoxyribonucleic acid


European molecular biology laboratory


European variation archive


False discovery rate


Generalized linear model


Genome-wide association study


Joint genotyper for inbred lines




Linkage disequilibrium


Minimal allele frequency




National center for biotechnology information


Nasonia vitripennis genome reference panel


Polymerase chain reaction


Quantitative trait locus


Restriction site associated DNA


Restricted maximum likelihood


Standard error


Single nucleotide polymorphism


Untranslated region


Variant call format


WD repeat-containing protein 16


  1. Charnov EL. The theory of sex allocation. Princeton, NJ: Princeton Univ Press; 1982.

    Google Scholar 

  2. Herre EA. Optimality, plasticity and selective regime in fig wasp sex ratios. Nature. 1987;329:627–9.

    Google Scholar 

  3. West SA. Sex allocation. Princeton, NJ: Princeton Univ Press; 2009.

    Google Scholar 

  4. Hamilton WD. Narrow roads of gene land : the collected papers of W.D. Hamilton. Oxford: W.H. Freeman/Spektrum; 1996.

    Google Scholar 

  5. West SA, Shuker DM, Sheldon BC. Sex-ratio adjustment when relatives interact: a test of constraints on adaptation. Evolution (N Y). 2005;59:1211–28.

    Google Scholar 

  6. Frank SA. Foundations of social evolution. Pinceton: Princeton University Press; 1998.

  7. Hardy ICW, editor. Sex ratios. Cambridge: Cambridge University Press; 2002.

    Book  Google Scholar 

  8. Macke E, Magalhães S, Bach F, Olivieri I. Experimental evolution of reduced sex ratio adjustment under local mate competition. Science (80- ). 2011;334:1127–1129.

  9. Macke E, Olivieri I, Magalhães S. Local mate competition mediates sexual conflict over sex ratio in a haplodiploid spider mite. Curr Biol. 2014;24:2850–4.

    Article  CAS  PubMed  Google Scholar 

  10. Hamilton WD. Extraordinary sex ratios. Science. 1967;156(80):477–88.

    Article  CAS  PubMed  Google Scholar 

  11. Hamilton WD. Reproductive competition and sexual selection in insects. In: Blum MS, Blum NA, editors. Sexual selection and reproductive competition in insects. New York: Academic Press; 1979. p. 167–220.

    Google Scholar 

  12. Taylor PD, Bulmer MG. Local mate competition and the sex ratio. J Theor Biol. 1980;86:409–19.

    Article  CAS  PubMed  Google Scholar 

  13. Reece SE, Drew DR, Gardner A. Sex ratio adjustment and kin discrimination in malaria parasites. Nature. 2008;453:609–14.

    Article  CAS  PubMed  Google Scholar 

  14. Godfray HCJ. Parasitoids : behavioral and evolutionary ecology. Princeton, N.J: Princeton University Press; 1994.

    Google Scholar 

  15. Burton-Chellew MN, Koevoets T, Grillenberger BK, Sykes EM, Underwood SL, Bijlsma K, et al. Facultative sex ratio adjustment in natural populations of wasps: cues of local mate competition and the precision of adaptation. Am Nat. 2008;172:393–404.

    PubMed  Google Scholar 

  16. Pannebakker BA, Halligan DL, Reynolds KT, Ballantyne GA, Shuker DM, Barton NH, et al. Effects of spontaneous mutation accumulation on sex ratio traits in a parasitoid wasp. Evolution (N Y). 2008;62:1921–35.

    Google Scholar 

  17. Cook N, Trivedi U, Pannebakker BA, Blaxter M, Ritchie MG, Tauber E, et al. Oviposition but not sex allocation is associated with transcriptomic changes in females of the parasitoid wasp Nasonia vitripennis. G3 genes, genomes. Genet. 2015;5:2885–92.

    CAS  Google Scholar 

  18. Cook N, Boulton RA, Green J, Trivedi U, Tauber E, Pannebakker BA, et al. Differential gene expression is not required for facultative sex allocation: a transcriptome analysis of brain tissue in the parasitoid wasp Nasonia vitripennis. R Soc Open Sci. 2018;5.

  19. Pannebakker BA, Watt R, Knott SA, West SA, Shuker DM. The quantitative genetic basis of sex ratio variation in Nasonia vitripennis: a QTL study. J Evol Biol. 2010;24.

  20. West SA, Herre EA. Stabilizing selection and variance in fig wasp sex ratios. Evolution (N Y). 1998;52:475–85.

    CAS  Google Scholar 

  21. Flanagan KE, West SA, Godfray HCJ. Local mate competition, variable fecundity and information use in a parasitoid. Anim Behav. 1998;56:191–8.

    CAS  PubMed  Google Scholar 

  22. Shuker DM, West SA. Information constraints and the precision of adaptation: sex ratio manipulation in wasps. Proc Natl Acad Sci U S A. 2004;101:10363–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Shuker DM, Pen I, Duncan AB, Reece SE, West SA. Sex ratios under asymmetrical local mate competition: theory and a test with parasitoid wasps. Am Nat. 2005;166:301–16.

    PubMed  Google Scholar 

  24. Shuker DM, Pen I, West SA. Sex ratios under asymmetrical local mate competition in the parasitoid wasp Nasonia vitripennis. Behav Ecol. 2006;17:345–52.

    Google Scholar 

  25. Ivens ABF, Shuker DM, Beukeboom LW, Pen I. Host acceptance and sex allocation of Nasonia wasps in response to conspecifics and heterospecifics. Proc R Soc B Biol Sci. 2009;276:3663–9.

    CAS  Google Scholar 

  26. Parker ED, Orzack SH. Genetic variation for the sex ratio in Nasonia vitripennis. Genetics. 1985;110:93–105.

    PubMed  PubMed Central  Google Scholar 

  27. Orzack SH, Parker ED, Gladstone J. The comparative biology of genetic variation for conditional sex ratio behavior in a parasitic wasp, Nasonia vitripennis. Genetics. 1991;127:583–99.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Orzack SH, Gladstone J. Quantitative genetics of sex ratio traits in the parasitic wasp, Nasonia vitripennis. Genetics. 1994;137:211–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Orzack SH. Sex-ratio control in a parasitic wasp, Nasonia vitripennis. II. Experimental analysis of an optimal sex-ratio model. Evolution (N Y). 1986;40:341.

    Article  Google Scholar 

  30. Orzack SH, Parker ED. Sex-ratio control in a parasitic wasp, Nasonia vitripennis. I. Genetic variation in facultative sex-ratio adjustment. Evolution (N Y). 1986;40:331–40.

    Article  Google Scholar 

  31. Barton N, Partridge L. Limits to natural selection. Bioessays. 2000;22:1075–1084.

  32. Nadeau NJ, Jiggins CD. A golden age for evolutionary genetics? Genomic studies of adaptation in natural populations. Trends Genet. 2010;26:484–92.

    CAS  PubMed  Google Scholar 

  33. Rockman MV. The QTN program and the alleles that matter for evolution: all that’s gold does not glitter. Evolution. 2012;66:1–17.

  34. Jones FC, Grabherr MG, Chan YF, Russell P, Mauceli E, Johnson J, et al. The genomic basis of adaptive evolution in threespine sticklebacks. Nature. 2012;484:55–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Wellenreuther M, Hansson B. Detecting polygenic evolution: problems, pitfalls, and promises. Trends Genet. 2016;32:155–64.

    Article  CAS  PubMed  Google Scholar 

  36. York RA. Assessing the genetic landscape of animal behavior. Genetics. 2018;209:223–32.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Sella G, Barton NH. Thinking about the evolution of complex traits in the era of genome-wide association studies. Annu Rev Genomics Hum Genet. 2019;20:461–93.

    Article  CAS  PubMed  Google Scholar 

  38. Hofmann HA. Functional genomics of neural and behavioral plasticity. J Neurobiol. 2003;54:272–82.

    Article  CAS  PubMed  Google Scholar 

  39. Toro MA, Charlesworth B. An attempt to detect genetic variation in sex ratio in </i>Drosophila melanogaster</i>. Heredity (Edinb). 1982;49:199–209.

    Google Scholar 

  40. Carvalho AB, Sampaio MC, Varandas FR, Klaczko LB. An experimental demonstration of Fisher’s principle: evolution of sexual proportion by natural selection. Genetics. 1998;148:719–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Frohschammer S, Heinze J. A heritable component in sex ratio and caste determination in a Cardiocondyla ant. Front Zool. 2009;6.

  42. Yusa Y. Genetics of sex-ratio variation inferred from parent-offspring regressions and sib correlations in the apple snail Pomacea canaliculata. Heredity (Edinb). 2006;96:100–5.

    CAS  Google Scholar 

  43. Vandeputte M, Dupont-Nivet M, Merdy O, Haffray P, Chavanne H, Chatain B. Quantitative genetic determinism of sex-ratio in the European sea bass (Dicentrarchus labrax L.). Aquaculture. 2007;272:S315.

  44. Janzen FJ. Heritable variation for sex ratio under environmental sex determination in the common snapping turtle (Chelydra serpentina). Genetics. 1992;131:155–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Postma E, Heinrich F, Koller U, Sardell RJ, Reid JM, Arcese P, et al. Disentangling the effect of genes, the environment and chance on sex ratio variation in a wild bird population. Proc R Soc B Biol Sci. 2011;278:2996–3002.

    Article  Google Scholar 

  46. Toro MA, Fernández A, García-Cortés LA, Rodrigáñez J, Silió L. Sex ratio variation in Iberian pigs. Genetics. 2006;173:911–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Gellatly C. Trends in population sex ratios may be explained by changes in the frequencies of polymorphic alleles of a sex ratio gene. Evol Biol. 2009;36:190–200.

    Article  Google Scholar 

  48. Antolin MF. Sex ratio variation in a parasitic wasp I. Reaction norms. Evolution (N Y). 1992;46:1496–510.

    Article  Google Scholar 

  49. Antolin MF. Sex ratio vriation in a parasitic wasp II. Diallel cross. Evolution (N Y). 1992;46:1511–24.

    Google Scholar 

  50. Kobayashi A, Tanaka Y, Shimada M. Genetic variation of sex allocation in the parasitoid wasp Heterospilus prosopidis. Evolution (N Y). 2003;57:2659–64.

    Google Scholar 

  51. Henter HJ. Constrained sex allocation in a parasitoid due to variation in male quality. J Evol Biol. 2004;17:886–96.

    CAS  PubMed  Google Scholar 

  52. Wajnberg E. Intra-population genetic variation in Trichogramma. In: Wajnberg E, Hassan SA, editors. Biological control with egg parasitoids. Wallingford: CAB International; 1994. p. 245–71.

    Google Scholar 

  53. Orzack SH, Parker ED. Genetic variation for sex ratio traits within a natural population of a parasitic wasp, Nasonia vitripennis. Genetics. 1990;124:373–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  54. Pannebakker BA, Trivedi U, Blaxter MA, Watt R, Shuker DM. The transcriptomic basis of oviposition behaviour in the parasitoid wasp Nasonia vitripennis. PLoS One. 2013;8:e68608.

  55. Van De Zande L, Ferber S, De Haan A, Beukeboom LW, Van Heerwaarden J, Pannebakker BA. Development of a Nasonia vitripennis outbred laboratory population for genetic analysis. Mol Ecol Resour. 2014;14:578–87.

    PubMed  Google Scholar 

  56. Stone EA. Joint genotyping on the fly: identifying variation among a sequenced panel of inbred lines. Genome Res. 2012;22:966–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. Raychoudhury R, Grillenberger BK, Gadau J, Bijlsma R, van de Zande L, Werren JH, et al. Phylogeography of Nasonia vitripennis (Hymenoptera) indicates a mitochondrial–Wolbachia sweep in North America. Heredity (Edinb). 2010;104:318–26.

    Article  CAS  Google Scholar 

  58. Hirschner W, Pogoda H-M, Kramer C, Thiess U, Hamprecht B, Wiesmüller K-H, et al. Biosynthesis of Wdr16, a marker protein for kinocilia-bearing cells, starts at the time of kinocilia formation in rat, and wdr16 gene knockdown causes hydrocephalus in zebrafish. J Neurochem. 2007;101:274–88.

    Article  CAS  PubMed  Google Scholar 

  59. Amaral A, Paiva C, Attardo Parrinello C, Estanyol JM, Ballescà JL, Ramalho-Santos J, et al. Identification of proteins involved in human sperm motility using high-throughput differential proteomics. J Proteome Res. 2014;13:5670–84.

    Article  CAS  PubMed  Google Scholar 

  60. Davies NJ, Tauber E. WaspAtlas: a Nasonia vitripennis gene database and analysis platform. Database. 2015;2015:bav103.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Akbari OS, Antoshechkin I, Hay BA, Ferree PM. Transcriptome profiling of Nasonia vitripennis testis reveals novel transcripts expressed from the selfish B chromosome, Paternal Sex Ratio. G3 Genes|Genomes|Genetics. 2013;3:1597–605.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Coutton C, Vargas AS, Amiri-Yekta A, Kherraf Z-E, Mustapha SFB, Le TP, et al. Mutations in CFAP43 and CFAP44 cause male infertility and flagellum defects in Trypanosoma and human. Nat Commun. 2018;9.

  63. Beurois J, Martinez G, Cazin C, Kherraf Z-E, Amiri-Yekta A, Thierry-Mieg N, et al. CFAP70 mutations lead to male infertility due to severe astheno-teratozoospermia. A case report. Hum Reprod. 2019;34:2071–9.

    Article  PubMed  Google Scholar 

  64. Verhulst EC, Beukeboom LW, Van De Zande L. Maternal control of haplodiploid sex determination in the wasp Nasonia. Science. 2010;328(80):620–3.

    CAS  PubMed  Google Scholar 

  65. Shuker DM, Sykes EM, Browning LE, Beukeboom LW, West SA. Male influence on sex allocation in the parasitoid wasp Nasonia vitripennis. Behav Ecol Sociobiol. 2006;59:829–35.

    Article  Google Scholar 

  66. Martel V, Shuker DM, Boulton RA, Damiens D, Boivin G. Sex allocation and the evolution of insemination capacity under local mate competition. Entomol Exp Appl. 2016;159:230–42.

    Article  Google Scholar 

  67. Harvey SE, Cheng C. Methods for characterization of alternative RNA splicing. Methods Mol Biol. 2016;1402:229–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Araujo PR, Yoon K, Ko D, Smith AD, Qiao M, Suresh U, et al. Before it gets started: regulating translation at the 5′; UTR. Comp Funct Genomics. 2012;2012:1–8.

    Article  CAS  Google Scholar 

  69. Wennerberg K, Rossman KL, Der CJ. The Ras superfamily at a glance. J Cell Sci. 2005;118(Pt 5):843–6.

    Article  CAS  PubMed  Google Scholar 

  70. Wei S, Xie C, Abe Y, Cai J. ADP-ribosylation factor like 7 (ARL7) interacts with α-tubulin and modulates intracellular vesicular transport. Biochem Biophys Res Commun. 2009;384:352–6.

    Article  CAS  PubMed  Google Scholar 

  71. Kahn RA, East MP, Francis JW. ARF-Like (ARL) Proteins. In: Ras superfamily small G proteins: biology and mechanisms 2. Cham: Springer International Publishing; 2014. p. 215–51.

    Chapter  Google Scholar 

  72. Jacobs S, Schürmann A, Becker W, Bockers TM, Copeland NG, Jenkins NA, et al. The mouse ADP-ribosylation factor-like 4 gene: two separate promoters direct specific transcription in tissues and testicular germ cell. Biochem J. 1998;335:259–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Schurmann A, Koling S, Jacobs S, Saftig P, Krau S, Wennemuth G, et al. Reduced sperm count and normal fertility in male mice with targeted disruption of the ADP-ribosylation factor-like 4 (Arl4) gene. Mol Cell Biol. 2002;22:2761–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Green RF, Gordh G, Hawkins BA. Precise sex ratios in highly inbred parasitic wasps ( Goniozus- emigratus/gordhi). Am Nat. 1982;120:653–65.

    Google Scholar 

  75. Nagelkerke CJ, Hardy ICW. The influence of developmental mortality on optimal sex allocation under local mate competition. Behav Ecol. 1994;5:401–11.

    Article  Google Scholar 

  76. Lynch JA. The expanding genetic toolbox of the wasp Nasonia vitripennis and its relatives. Genetics. 2015;199:897–904.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Werren JH, Richards S, Desjardins CA, Niehuis O, Gadau J, Colbourne JK, et al. Functional and evolutionary insights from the genomes of three parasitoid Nasonia species. Science (80- ). 2010;327:343–8.

  78. Desjardins CA, Perfectti F, Bartos JD, Enders LS, Werren JH. The genetic basis of interspecies host preference differences in the model parasitoid Nasonia. Heredity (Edinb). 2010;104:270–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Niehuis O, Buellesbach J, Gibson JD, Pothmann D, Hanner C, Mutti NS, et al. Behavioural and genetic analyses of Nasonia shed light on the evolution of sex pheromones. Nature. 2013;494:345–8.

    Article  CAS  PubMed  Google Scholar 

  80. Davies NJ, Tauber E. Deep sequencing analysis of the circadian transcriptome of the jewel wasp Nasonia vitripennis. bioRxiv. 2016:048736.

  81. Paolucci S, Dalla Benetta E, Salis L, Doležel D, van de Zande L, Beukeboom L. Latitudinal variation in circadian rhythmicity in Nasonia vitripennis. Behav Sci (Basel). 2019;9:115.

    Article  Google Scholar 

  82. Dalla Benetta E, Beukeboom LW, van de Zande L. Adaptive differences in circadian clock gene expression patterns and photoperiodic diapause induction in Nasonia vitripennis. Am Nat. 2019;193:881–96.

    Article  PubMed  Google Scholar 

  83. Rago A, Werren JH, Colbourne JK. Sex biased expression and co-expression networks in development, using the hymenopteran Nasonia vitripennis. PLoS Genet. 2020;16:e1008518.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Cohen LB, Edwards R, Moody D, Arsala D, Werren JH, Lynch JA. Genetic, morphometric, and molecular analyses of interspecies differences in head shape and hybrid developmental defects in the wasp genus Nasonia. bioRxiv. 2019:663732.

  85. Cook N, Pannebakker BA, Tauber E, Shuker DM. DNA methylation and sex allocation in the parasitoid wasp Nasonia vitripennis. Am Nat. 2015;186:514–8.

    Google Scholar 

  86. Cook N, Parker DJ, Tauber E, Pannebakker BA, Shuker DM. Validating the demethylating effects of 5-aza-2′-deoxycytidine in insects requires a whole-genome approach: (a reply to Ellers et al.). Am Nat. 2019;194:432–8.

    Article  PubMed  Google Scholar 

  87. Henter HJ. Inbreeding depression and haplodiploidy: experimental measures in a parasitoid and comparisons across diploid and haplodiploid insect taxa. Evolution (N Y). 2003;57:1793–803.

    Google Scholar 

  88. Mather KA, Caicedo AL, Polato NR, Olsen KM, McCouch S, Purugganan MD. The extent of linkage disequilibrium in rice (Oryza sativa L.). Genetics. 2007;177:2223–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Hyten DL, Choi I-Y, Song Q, Shoemaker RC, Nelson RL, Costa JM, et al. Highly variable patterns of linkage disequilibrium in multiple soybean populations. Genetics. 2007;175:1937–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Mackay TFC, Richards S, Stone EA, Barbadilla A, Ayroles JF, Zhu D, et al. The Drosophila melanogaster genetic reference panel. Nature. 2012;482:173–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  91. Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N, et al. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet. 2014;46:1081–8.

    CAS  PubMed  Google Scholar 

  92. Niehuis O, Gibson JD, Rosenberg MS, Pannebakker BA, Koevoets T, Judson AK, et al. Recombination and its impact on the genome of the haplodiploid parasitoid wasp Nasonia. PLoS One. 2010;5.

  93. Desjardins CA, Gadau J, Lopez JA, Niehuis O, Avery AR, Loehlin DW, et al. Fine-scale mapping of the Nasonia genome to chromosomes using a high-density genotyping microarray. G3 genes, genomes. Genet. 2013;3:205–15.

    CAS  Google Scholar 

  94. Beukeboom LW, Niehuis O, Pannebakker BA, Koevoets T, Gibson JD, Shuker DM, et al. A comparison of recombination frequencies in intraspecific versus interspecific mapping populations of Nasonia. Heredity (Edinb). 2010;104:302–9.

    CAS  Google Scholar 

  95. Wallberg A, Glémin S, Webster MT. Extreme recombination frequencies shape genome variation and evolution in the honeybee, Apis mellifera. PLOS Genet. 2015;11:e1005189.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Grillenberger BK, Koevoets T, Burton-Chellew MN, Sykes EM, Shuker DM, Van De Zande L, et al. Genetic structure of natural Nasonia vitripennis populations: validating assumptions of sex-ratio theory. Mol Ecol. 2008;17:2854–64.

    CAS  PubMed  Google Scholar 

  97. Grillenberger BK, Gadau J, Bijlsma R, Van De Zande L, Beukeboom LW. Female dispersal and isolation-by-distance of Nasonia vitripennis populations in a local mate competition context. Entomol Exp Appl. 2009;132:147–54.

    Google Scholar 

  98. Davey JW, Hohenlohe PA, Etter PD, Boone JQ, Catchen JM, Blaxter ML. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat Rev Genet. 2011;12:499–510.

    CAS  PubMed  Google Scholar 

  99. Narum SR, Buerkle CA, Davey JW, Miller MR, Hohenlohe PA. Genotyping-by-sequencing in ecological and conservation genomics. Mol Ecol. 2013;22:2841–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  100. White LC, Fontsere C, Lizano E, Hughes DA, Angedakin S, Arandjelovic M, et al. A roadmap for high-throughput sequencing studies of wild animal populations using noninvasive samples and hybridization capture. Mol Ecol Resour. 2019;19:609–22.

    Article  CAS  PubMed  Google Scholar 

  101. Stapley J, Reger J, Feulner PGD, Smadja C, Galindo J, Ekblom R, et al. Adaptation genomics: the next generation. Trends Ecol Evol. 2010;25:705–12.

    Article  PubMed  Google Scholar 

  102. Gienapp P, Calus MPL, Laine VN, Visser ME. Genomic selection on breeding time in a wild bird population. Evol Lett. 2019;3:142–51.

    Article  PubMed  PubMed Central  Google Scholar 

  103. Bosse M, Spurgin LG, Laine VN, Cole EF, Firth JA, Gienapp P, et al. Recent natural selection causes adaptive evolution of an avian polygenic trait. Science. 2017;358:365–8.

    Article  CAS  PubMed  Google Scholar 

  104. Whiting AR. The biology of the parasitic wasp Mormoniella vitripennis [=Nasonia vitripennis]. Q Rev Biol. 1967;42:333–406.

    Google Scholar 

  105. Maniatis T, Frirsch E, Sambrook J. Molecular: a laboratory manual. New York: Cold Spring Harbor Laboratory; 1982.

    Google Scholar 

  106. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25:1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    PubMed  PubMed Central  Google Scholar 

  108. Kofler R, Orozco-terWengel P, de Maio N, Pandey RV, Nolte V, Futschik A, et al. Popoolation: a toolbox for population genetic analysis of next generation sequencing data from pooled individuals. PLoS One. 2011;6.

  109. R Development Core Team. R: A language and enviroment for statistical computing. 2018.

  110. Sul JH, Martin LS, Eskin E. Population structure in genetic studies: confounding factors and mixed models. PLoS Genet. 2018;14:e1007309.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Zhou W, Rousset F, O’Neill S. Phylogeny and PCR-based classification of Wolbachia strains using wsp gene sequences. Proc R Soc B Biol Sci. 1998;265:509–15.

    CAS  Google Scholar 

  112. Xia S, Pannebakker BA, Groenen MAM, Zwaan BJ, Bijma P. Quantitative genetics of wing morphology in the parasitoid wasp Nasonia vitripennis: hosts increase sibling similarity. Heredity (Edinb). 2020:1–10.

  113. Werren JH. A model for sex ratio selection in parasitic wasps: local mate competition and host quality effects. Netherlands J Zool. 1984;34:81–96.

  114. Wajnberg E. Measuring genetic variation in natural enemies used for biological control: why and how? Genet Evol Biol Control. 2004:19–37.

  115. Lynch M, Walsh B. Genetics and Analysis of Quantitative Traits. Sunderland: Sinauer Assoc; 1998.

  116. Hoffmann AA, Parsons P. The analysis of quantitative variation in natural populations with isofemale strains. Genet Sel Evol. 1988;20:87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  117. Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team. {nlme}: Linear and nonlinear mixed effects models. 2020.

  118. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67:1–48.

    Article  Google Scholar 

  119. Wang M, Zhao Y, Zhang B. Efficient test and visualization of multi-set intersections. Sci Rep. 2015;5:16923.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  120. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80–92.

    Article  CAS  Google Scholar 

Download references


We would like to thank Gabriella Bukovinszkine Kiss, Bertha Koopmanschap (in loving memory) and Jade Green for their help in the laboratory, and Eveline Verhulst for her insightful comments on the molecular biology of sex determination. We also would like to thank two anonymous reviewers for their helpful comments on an earlier version of this paper.


This research was funded by the Netherlands Genomics Initiative (Zenith 93511041) and the by the Natural Environment Research Council (NE/J024481/1). The funding agencies played no part in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations



BAP, LvdZ and DMS conceived the study. NC performed the sex ratio experiments, BAP, NC and JvdH perfomed the statistical and sequence analyses. BAP and DMS wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Bart A. Pannebakker.

Ethics declarations

Ethics approval and consent to participate

The Nasonia vitripennis Genetic Reference Panel strains presented in this article are derived from the HVRx laboratory outbred population, which has been reared in our laboratory since 2012. We obtained our culture from the University of Groningen, The Netherlands, who started it from field lines collected from the Hoge Veluwe National Park in the Netherlands following permission of the Hoge Veluwe management. The authors declare that the collections of specimens comply with institutional, national, or international guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Supplementary Table 1.

Sequencing statistics NVGRP. For each NVGRP iso-female line, the total number of clean reads (million reads), the total number of clean bases (Gb), percentage of nucleotides with a read quality higher than 20 (Q20 percentage), GC content (GC %), total number of reads mapped to Nvit 2.1 assembly (million reads), total number of nucleotides covered (Mb) and the percentage of Nvit 2.1 assembly covered (percentage) is given. Data in SupplementaryTable1.

Additional file 2: Supplementary Table 2

. NVGRP and HVRx Sliding windows analysis (400 kb windows) for pi, theta and Tajima’s D. Data in SupplementaryTable2_PopoolationResults.

Additional file 3: Supplementary Table 3

. Results for sex ratio and clutch size GWAS in NVGRP. Worksheets show: Significantly associated SNPs for sex ratio, All GWAS results for sex ratio and All GWAS results for clutch size. Data in SupplementaryTable3_GWAS_Candidate_and_All_SNPs.

Additional file 4: Supplementary Table 4

. Linkage disequilibrium (r2) between 18 SNPs significantly associated with sex ratio in the 25 NVGRP lines. Data in SupplementaryTable4_r2_Matrix.

Additional file 5: Supplementary Figure 1

. Hierarchical clustering of all 34 NVGRP lines based on the average pairwise FST values over all variable sites between lines. Plot shows clustering based on the complete linkage method in R package hclust. Supplementary Figure 2. Distribution of phenotypic values over all 26 tested NVGRP lines, for (a) sex ratio expressed as proportion males, and (b) clutch size. Supplementary Figure 3. Boxplot of sex ratio (proportion males) and clutch size among all 26 tested NVGRP lines. Plot shows the median, the box indicates the first and the third quartile, the whiskers indicate 1.5 times the interquartile range and values beyond that are labelled outliers (open circles). Supplementary Figure 4. Scatterplot of the mean clutch size and the mean sex ratio for the 25 tested NVGRP lines (excluding outlier line 62). Solid line shows the fit of the linear regression y = 0.224467–0.001477x. Supplementary Figure 5. Multipaneled plot showing the cumulative results of the SuperExactTest for association of the p-values of the sex ratio and clutch size GWAS per chromosome for windows of 25 kb, 50 kb, 100 kb, 200 kb and 400 kb. Black solid line shows observed overlap using GWAS p values for each SNPs, red solid line shows the most significant p-value for each window given from a distribution of 100 permutations. Significant overlap is observed for all window sizes, when the p-values of the association between the sex ratio and clutch size GWAS exceeds the p-values of the permutation.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Pannebakker, B.A., Cook, N., van den Heuvel, J. et al. Genomics of sex allocation in the parasitoid wasp Nasonia vitripennis. BMC Genomics 21, 499 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: