Population-level consequences of complementary sex determination in a solitary parasitoid
© de boer et al. 2015
Received: 12 March 2015
Accepted: 26 March 2015
Published: 30 May 2015
Sex determination mechanisms are known to be evolutionarily labile but the factors driving transitions in sex determination mechanisms are poorly understood. All insects of the Hymenoptera are haplodiploid, with males normally developing from unfertilized haploid eggs. Under complementary sex determination (CSD), diploid males can be produced from fertilized eggs that are homozygous at the sex locus. Diploid males have near-zero fitness and thus represent a genetic load, which is especially severe under inbreeding. Here, we study mating structure and sex determination in the parasitoid Cotesia vestalis to investigate what may have driven the evolution of two complementary sex determination loci in this species.
We genotyped Cotesia vestalis females collected from eight fields in four townships in Western Taiwan. 98 SNP markers were developed by aligning Illumina sequence reads of pooled DNA of eight different females against a de novo assembled genome of C. vestalis. This proved to be an efficient method for this non-model species and provides a resource for future use in related species. We found significant genetic differentiation within the sampled population but variation could not be attributed to sampling locations by AMOVA. Non-random mating was detected, with 8.1% of matings between siblings. Diploid males, detected by flow cytometry, were produced at a rate of 1.4% among diploids.
We think that the low rate of diploid male production is best explained by a CSD system with two independent sex loci, supporting laboratory findings on the same species. Fitness costs of diploid males in C. vestalis are high because diploid males can mate with females and produce infertile triploid offspring. This severe fitness cost of diploid males combined with non-random mating may have resulted in evolution from single locus CSD to CSD with two independent loci.
Sex determination systems that initiate differentiation between males and females are highly diverse across the animal kingdom, with examples of temperature dependent sex determination (e.g. in reptiles), heterogametic males in mammals and heterogametic females in birds, and a range of other (genetic) mechanisms in amphibians, reptiles and insects [1-4]. Transitions among sex-determining mechanisms have occurred repeatedly across the animal kingdom [1,5], and genetic experiments on the nematode Caenorhabditis elegans confirm that sex determining mechanisms are indeed extremely labile . Despite growing knowledge on how sex determination mechanisms can evolve, the evolutionary drivers of these mechanisms are understood less well [4,7,8].
The insect order Hymenoptera is interesting in this respect because of the intricate relationship between mating structure and sex determination [9,10]. All sexually reproducing Hymenoptera are haplodiploid, with haploid males developing from unfertilized haploid eggs, while females develop from fertilized diploid eggs. Sex is not determined by fertilization alone, however, and diploid males can develop from fertilized eggs that are homozygous at a sex locus under a system called complementary sex determination (CSD) . Heterozygosity at the sex locus leads to female development. The presence of diploid males, suggesting CSD, has been demonstrated in a variety of parasitoid wasps, ants, solitary and social bees and wasps [9,12]. Genetic support for the sex locus is lacking in most of these species due to the limited availability of genetic markers and the difficulty to make controlled inbred crosses in the laboratory, particularly for social species . An exception is formed by bumblebees, in which the sex locus has been mapped [13-15], and honeybees, in which the csd gene has been mapped, cloned and sequenced [13,14,16,17].
In terms of genetics, CSD resembles the self-incompatibility (SI) breeding system of plants, and the major histocompatibility (MHC) locus in vertebrates because heterozygotes in all three systems have a fitness benefit , and high allelic diversity can be maintained through frequency dependent selection. The fitness costs of CSD are potentially even more severe than those associated with SI and MHC because the fusion of gametes with incompatible sex alleles effectively leads to post-zygotic mortality. Diploid males resulting from sex locus homozygotes have near-zero fitness because they are commonly sterile or unviable, and they are produced at the expense of fertile daughters. It is thus expected that populations with CSD experience strong selection on reducing the costs associated with the production of diploid males . In species in which CSD is maintained, selection may act to reduce levels of inbreeding or the cost of diploid males. For example, in Bracon hebetor and Bombus terrestris, behavioural adaptations reduce the level of mating between siblings [19-21], and in Euodynerus foraminatus and Cotesia glomerata, diploid males are able to sire diploid daughters, presumably by producing haploid sperm [22,23]. Alternatively, selection may lead to the evolution from CSD to mechanisms of sex determination with fewer or no diploid males. Interestingly, within the Hymenoptera only one other mechanism of sex determination is currently supported experimentally: maternal imprinting in the parasitoid Nasonia vitripennis [24,25].
We previously showed that CSD in Cotesia vestalis (Braconidae) is likely caused by two unlinked sex loci, based on the rate of diploid male production in laboratory crosses . The genus Cotesia is polymorphic for CSD variants as well as for other life history parameters, which makes it a valuable system to study evolutionary drivers of sex determination . Cotesia vestalis is a solitary Eurasian parasitoid [27,28], which has been introduced for biological control of diamondback moth worldwide . Solitary parasitoids are expected to have lower rates of inbreeding than gregarious parasitoids because siblings emerge from the same host in the latter . This makes it surprising that C. vestalis has evolved a second sex locus whereas the closely related C. glomerata has maintained single locus CSD (sl-CSD) despite high inbreeding rates in nature [31-33]. This raises the question what has driven the evolution of two-locus CSD in C. vestalis. Here, we investigate the hypothesis that non-random mating in nature may have selected for the evolution of two-locus CSD. We expect that under natural conditions, C. vestalis (1) produces diploid males at a low rate; (2) shows non-random mating; and (3) has little population structure. Because very few molecular markers were available for C. vestalis, we generated a superior genomic resource by sequencing the entire genome of C. vestalis by Illumina technology. After de novo assembly we developed a set of single nucleotide polymorphism (SNP) assays for population-wide studies in this species. This approach also provides a backbone for similar studies in closely related species that may lead to further comparative studies of CSD.
Collection of field samples
Cotesia vestalis (Haliday) was collected in Western Taiwan in February of 2008. Two small non-commercial vegetable gardens with apparent diamondback moth infestations were selected in each of four townships: Luhzu (Kaohshiung county), Sihu (Changhua county), Shanhua (Tainan county) and Shingang (Chiayi county). A variety of crops was grown in these fields, including different types of host plants for Plutella xylostella L. (diamondback moth larvae), such as Chinese cabbage, broccoli, stem cabbage, kohlrabi and cauliflower. In most fields, other, non-host crops were grown as well, or were bordering the cabbage field, such as tomato, corn, sweet pepper, papaya, beans and onions.
In each field, 15 to 24 host plants were randomly selected from which we collected all C. vestalis cocoons and diamondback moth larvae (first to third instar) (Additional file 1: Table S1). Cotesia vestalis cocoons were directly transferred to small Petri dishes (5 cm diameter) with some droplets of honey and kept in the refrigerator until transportation to the Netherlands (a maximum of 2 weeks). Diamondback moth larvae were transferred to plastic containers (18x13x6cm) with a square hole in the lid covered with fine mesh for ventilation. Larvae were provided with fresh cabbage leaves every other day for one week at 20-30°C. By then, most larvae had either pupated or a C. vestalis cocoon had emerged. After one week, C. vestalis cocoons were collected into small Petri dishes and kept in the refrigerator as described above; remaining material was discarded. All C. vestalis cocoons collected and/or reared from a single plant were kept separately.
After transportation to the Netherlands, C. vestalis cocoons were kept at room temperature in a laboratory at the University of Groningen. Petri dishes were checked daily for adult emergence and emerged wasps sexed: Females were identified by the presence of an ovipositor with a stereomicroscope. Males and females were then frozen at −20°C for analysis of ploidy.
Analysis of ploidy
To assess the consequences of CSD in a natural field population, we analysed ploidy levels of approximately 2/3rd of all C. vestalis males, including all males of fields 1, 3, 4, 6 and 7 and a random subset of males from fields 2, 5 and 8 (Additional file 1: Table S2). We followed methods described previously for Cotesia parasitoids [34,35]. In short, the head of an individual wasp was ground in 0.5 ml ice-cold Galbraith buffer  with a Dounce tissue grinder. The cell suspension was sieved (40 μm) and stained with propidium iodide (25 μg per sample). DNA content of 2500 nuclei from head tissue was measured per wasp on an Epics® XL™ flow cytometer (Beckman Coulter, Brea, CA). DNA histograms were compared to those of known haploid males and diploid females to classify ploidy levels as haploid, diploid or unknown. After the head of a wasp was cut off, the rest of the body was placed in 96% EtOH and kept at 4°C for genotyping.
SNP detection, assay development and genotyping
De novo assembly and the Cotesia vestalis reference genome
No genome information was available for C. vestalis. In order to build a draft reference genome and to develop SNP assays, we sequenced the entire genome of C. vestalis on a single lane of paired-end sequences (2x100 bp) on an Illumina HiSeq 2000 (Illumina Inc., U.S.A.) instrument. The SNP discovery panel consisted of eight C. vestalis females, one from each field, from a randomly selected plant. DNA from the entire bodies was extracted as a pool using Qiagen’s DNeasy kit following the protocol for tissue DNA with an extended incubation time in buffer ATL (3 h at 56°C). A single genomic library with an insert size of approximately 300 bp was prepared using the Illumina Sample Preparation Kit and sequenced for 100 cycles with the Illumina HiSeq 2000 (paired-end) at the UMCG Sequencing Facility in Groningen, The Netherlands. The raw data files from the sequencing instrument are deposited in the NCBI short read archive under accession number SRP058413 .
Before assembly, Illumina reads were trimmed using an in-house Perl script that trims the sequence as soon as two consecutive bases have a quality score lower than 20. Reads that after trimming had a length smaller than 50 bp were removed from the analysis. To obtain C. vestalis sequence contigs to be used as a pseudo-reference genome, we performed a de novo assembly on the 133 million 100 bp reads using SOAPdenovo version 1.05 . The assembly was done using a k-mer size of 45 and k-mers that were seen only once were removed (option –d). After contig construction, scaffolding was performed using intra-scaffold closure (option –F) and a minimum length for scaffolding of 50 bp. The total size of the assembly was 152 Mb with a contig N50 size of 761 bp and a scaffold N50 size of 2400 bp. With an estimated size of the Cotesia vestalis genome of 190 Mb (estimated with flow cytometry by J.G. de Boer, unpublished data) this suggests that our assembly covers around 80% of the Cotesia vestalis genome. We subsequently concatenated all scaffolds and unassembled contigs into a single artificial chromosome to be used as a reference genome for SNP identification. Assembled scaffolds and contigs have been deposited in Dryad . Contigs and scaffolds were concatenated at random using a spacer sequence of 150 N’s.
Alignments of reads and SNP identification
Individual paired-end reads were aligned against the artificial Cotesia vestalis reference genome obtained from the de novo genome assembly using BWA [37,38]. The resulting BAM file was then used for the identification of putative SNPs using SAMtools and varFilter from the samtools.pl utility . We only considered nucleotide substitutions and ignored small indels. SNPs were filtered that had a mapping quality higher than 20, a minimum read depth of 3 and a maximum read depth of 90 (3x the average read depth, a strategy to avoid orthologous SNPs, e.g. in multi copy genes [40,41]).
SNP selection and assay development
From our list of putative SNPs across the C. vestalis genome, we selected 100 SNPs for genotyping assay development. We first selected the 200 largest scaffolds; they varied in length from 17-58Kb and contained a total of 7,878 SNPs. We then removed SNPs with a minor allele frequency (MAF) <0.2, SNPs that had another SNP within 50 bp up- or downstream, and SNPs with more than 2 alleles. The remaining SNPs were binned in MAF bins of 0.2-0.3 (1,908 SNPs), 0.3-0.4 (1,605 SNPs) and 0.4-0.5 (1,156 SNPs). Per MAF bin, SNPs were ranked by SNP quality score. We then selected the SNPs with the highest quality scores, picking 20 SNPs with a MAF between 0.2-0.3 and 40 each with a MAF between 0.3-0.4 and 0.4-0.5, all on different scaffolds. All selected SNPs had a quality score of more than 200 (based on SAMtools), and an average read depth of 61. High-throughput genotyping assays based on allele-specific forward primers were developed for these 100 SNP sequences at KBioscience (now LGC Genomics, Hoddesdon, U.K.). Kompetitive allele specific PCR (KASP) technology has been used successfully on a wide range of organisms, including plants, vertebrate and invertebrate animals, and, for plants, compares favourably to chip-based methods such as GoldenGate in terms of accuracy and cost .
Because the entire bodies of C. vestalis females were used in the DNA pool of the SNP discovery panel, we selected another female from the same plant in each field to validate the selected SNPs. DNA was extracted from the entire body using Qiagen’s DNeasy kit as described above. These eight females were subsequently genotyped with the KASP assays at KBioscience to establish whether the selected SNP loci were polymorphic.
To assess C. vestalis mating structure in the Taiwanese field population, we genotyped one female per plant of each field at 98 SNPs that were successfully developed into KASP assays. In addition, all males that were identified as diploid with flow cytometry were genotyped, males identified as haploids were not genotyped to save costs. DNA extractions were performed on entire parasitoid wasp bodies at KBioscience with the Kleargene protocol.
All wasps and all loci were included in the analyses because very few genotypes were missing (on average 2.5% of loci were missing per wasp). On average 3.5% of genotypes was missing per locus with a maximum of 11 out of 139 wasps failing to genotype at any one locus.
We first used Structure 2.3.3 [43,44] to investigate genetic clustering within the complete data set without prior information on sampling locations. We performed 10 independent runs each for K = 1 to K = 8 (where K is the number of clusters), using the admixture ancestry model with correlated alleles, a burn-in length of 50,000 steps and 100,000 MCMC steps. Evanno’s delta K method  was used to assess K.
The largest cluster (76 individuals) inferred by Structure was subsequently used to examine the SNP loci in terms of deviations from linkage and Hardy-Weinberg equilibria because we generated the markers de novo and had no information on their suitability as population genetic markers. The other clusters were considered too small for robust statistical analysis. Statistical comparison of our data to expectations under linkage equilibrium were done with likelihood ratio tests with 16,000 permutations and setting the number of initial conditions for expectation maximization to 5  in Arlequin 188.8.131.52 . To account for the large number of statistical tests in this pairwise procedure, p-levels for statistical significance were adjusted with a sequential Bonferroni correction for multiple comparisons . Deviations from Hardy-Weinberg equilibrium (HWE) in the largest cluster were also tested in Arlequin with a Markov chain algorithm with 1,000,000 steps and 100,000 dememorization steps.
Following these steps, 17 SNP loci were excluded (Additional file 1: Table S3) before calculating Wright’s F-statistics  per field and globally for the entire sample, using Arlequin. Pairwise F ST values  were computed to quantify genetic differentiation between fields. We also did a locus-by-locus analysis of molecular variance (AMOVA) with fields grouped into townships to hierarchically partition genetic variance into regions, fields and individuals.
The degree of sibmating α was then used to predict the proportion of diploids that was male (DMP). However, the occurrence of diploid males is the result of matings that are matched in terms of sex alleles, and not of sibmatings per se. In order to predict DMP, we thus need to predict the frequency of matched matings θ first . In a randomly mating population, θ depends on sex allele diversity, with a higher sex allele diversity resulting in a lower θ. Under inbreeding, θ is increased because brothers and sisters share sex alleles at a higher rate, independent of sex allele diversity in the population. Ultimately, θ and DMP depend on: (1) the number of complementary sex loci; (2) sex allele diversity at each locus; (3) degree of sibmating α; and (4) diploid male developmental survival relative to that of females. For C. vestalis, α is derived from our field study using formula (1), and developmental survival of diploid males appeared to be high in laboratory studies , so we will consider it equal to that of females. Laboratory crosses suggest the presence of CSD with two independent loci in C. vestalis, but we will also consider sl-CSD for comparison.
where α = degree of sibmating and k = number of sex alleles.
where α = degree of sibmating, f HET is the frequency of females heterozygous at both sex loci under random mating, and p HET and p HOM are the frequency of matched matings involving females heterozygous at both sex loci or homozygous at one sex locus respectively.
C. vestalis parasitism and diploid males in Taiwan
Overview of data collected from eight fields in Western Taiwan
Sex ratio of Cv
73.7 ± 7.7
6.3 ± 1.0
0.09 ± 0.02
0.57 ± 0.08
639.3 ± 92.9
97.4 ± 12.1
0.15 ± 0.01
0.47 ± 0.02
140.7 ± 14.8
24.1 ± 3.2
0.17 ± 0.02
0.48 ± 0.03
67.0 ± 14.3
9.1 ± 1.5
0.14 ± 0.03
0.41 ± 0.13
115.0 ± 15.7
49.1 ± 6.4
0.43 ± 0.04
0.44 ± 0.04
40.8 ± 6.7
8.7 ± 1.0
0.21 ± 0.05
0.50 ± 0.06
165.0 ± 19.2
7.3 ± 1.3
0.04 ± 0.01
0.51 ± 0.09
257.4 ± 30.0
130.3 ± 16.4
0.51 ± 0.03
0.53 ± 0.02
0.23 ± 0.08
0.49 ± 0.02
Analysis of male ploidy levels revealed the presence of 21 diploid males among 1,405 analysed males (Additional file 1: Table S2). Triploids were not found. The proportion of diploid males among all diploids (calculated per field as: diploid males/(diploid males + females), DMP) was thus low in all fields, varying between 0 and 0.03 (Table 1), with an overall DMP of 0.014. Diploidy was confirmed by genotyping at 98 SNPs for all 21 males by the presence of one or more heterozygous SNP loci.
SNP detection and validation
By aligning individual reads against the de novo assembled reference genome of C. vestalis, we detected 270,476 putative SNPs across the genome (on average one SNP per 590 bp)(sequence of scaffolds and contigs of reference genome, and information about putative SNPs are available in Dryad . Average depth of sequencing at putative SNP sites was 43X. As a first indication of the validity of our SNP calling procedure, we calculated the ratio of transitions:transversions, which was 2.18:1 for our dataset, well in the range of published ratios [40,41,55,56] and clearly different from the ratio that would be produced by chance (1:2) [57-59].
Of the 100 putative SNPs selected for assay development , 98 were polymorphic in the validation panel of eight females. This translates into a SNP conversion rate of 98%, comparing favourably with published studies, because conversion can be as low as around 50%  and references therein.
The set of 98 polymorphic SNPs was used to study genetic structure of 139 C. vestalis females collected in Western Taiwan (Additional file 1: Table S3 and dataset in Dryad . Bayesian clustering without prior information on sampling location suggested that the most likely number of populations in our sample was three (K = 3 using Structure and Evanno’s delta K method ).
We then investigated linkage disequilibrium and deviations from Hardy-Weinberg equilibrium of our SNP marker set to evaluate their suitability as population genetic markers. Individuals of the largest cluster identified by Structure were used for this purpose, using a conservative threshold of 0.9 probability to belong to a specific genetic cluster (76 individuals, Additional file 1: Table S4). Six out of 98 loci deviated significantly from HWE (Additional file 1: Table S3), one of these loci showed heterozygote excess and five were heterozygote deficient. Analysis of linkage disequilibrium indicated significant linkage between 24 (out of 4,753) pairs of loci (after sequential Bonferroni correction) involving 25 loci (Additional file 1: Table S3). Subsequent analyses presented in this paper were done with 81 loci after excluding the 6 loci that deviated from HWE and randomly excluding 11 additional loci so that no pairs of loci were in significant LD anymore (Additional file 1: Table S3). Analyses using the complete set of 98 loci are reported in Additional file 1: (Tables S5, S6 and S7) and revealed no qualitative differences.
Pairwise F ST for C. vestalis parasitoids collected from eight fields in Western Taiwan
Hierarchical partitioning of genetic variance
Sum of squares
Between fields within townships
Mating structure and sex determination load
Genetic variation within individuals of C. vestalis in Taiwan
Obs. Het. ± se
Exp. Het. ± se
0.348 ± 0.018
0.390 ± 0.014
0.394 ± 0.019
0.395 ± 0.013
0.402 ± 0.018
0.401 ± 0.014
0.412 ± 0.020
0.417 ± 0.014
0.382 ± 0.015
0.393 ± 0.013
0.373 ± 0.017
0.387 ± 0.014
0.420 ± 0.020
0.402 ± 0.013
0.359 ± 0.021
0.370 ± 0.016
Most studies on complementary sex determination have been performed in the laboratory. Yet, studying the genetic load associated with CSD (i.e. diploid males) in relation to mating structure of natural populations allows placing sex determination in an evolutionary context. Mating structure and sex determination are tightly intertwined in hymenopteran insects because the production of diploid males increases dramatically under conditions of inbreeding. Hence, mating systems characterized by inbreeding are thought to select for systems of sex determination with fewer or no diploid males . To date, the best example is found in Nasonia parasitoids, which naturally have high levels of sibmating and have no CSD . Instead, sex is determined through maternal imprinting in N. vitripennis, which does not lead to the production of diploid males under inbreeding [24,25]. Conversely, in species with CSD, selection is expected to act on factors that minimize the genetic load associated with CSD, i.e. to minimize the production or cost of diploid males. Here, we made a population-level analysis of CSD in a natural population of the solitary parasitoid wasp Cotesia vestalis by assessing DMP, and evaluating mating and population structure.
Diploid males in natural populations
We detected 21 diploid males in our field study in Western Taiwan, and the overall rate of males among diploids was 0.014, which is low compared to other Cotesia species. A Swiss population of C. glomerata contained diploid males at a rate of 0.06 , while DMP was 0.08-0.13 in a North American population of C. rubecula . Low DMP was expected in C. rubecula based on its two-locus CSD phenotype. However, this species was introduced in North America from its native range in Europe for biological control purposes and a genetic bottleneck might have reduced the sex determination system to effectively sl-CSD, hence resulting in higher rates of diploid males. Similar effects of genetic bottlenecks on DMP have been shown in Soenopsis invicta ants and Polistes chinensis antennalis wasps [63,64]. Conversely, C. glomerata was expected to have higher DMP because it has sl-CSD and a mating system characterized by a relatively high level of inbreeding [32,33]. Low DMP in C. glomerata may be explained by post-copulatory processes, such as biased fertilization of eggs with compatible sperm (in terms of sex alleles) .
Diploid males arise from matings that are matched in terms of sex alleles, and thus DMP depends on the frequency of matched matings, which in turn depends on the degree of sibmating and sex allele diversity [11,54]. We used the estimated rate of sibmating found in the field to predict DMP for C. vestalis for two CSD scenarios - single locus and two independent loci - and a wide range of sex allele numbers (Figure 2). These predictions illustrate that sl-CSD only results in very low DMP when sex allele diversity is very high; even with 100 different sex alleles, 8.1% sibmating still results in DMP of 0.03. Although most empirical estimates of sex allele diversity are in the range of 9–20 different alleles , recent csd sequence analyses in honeybees dramatically increased estimates of sex allele diversity from 11 to 20 alleles to almost 100 alleles locally ([54,65,66]. Two native populations of the fire ant S. invicta also harbour 66–86 sex alleles . These findings suggest that high numbers of sex alleles can be maintained in finite populations despite the effects of random genetic drift [67,68]. Mate choice could also result in lower than expected DMP, for example based on kinship or on sex allele identity . However, lab studies suggest that Cotesia females are not particularly choosy when it comes to males and they readily mate with brothers or with diploid males [34,35,70]. Moreover, controlled laboratory crosses, without an effect of mate choice, also resulted in low DMP in C. vestalis . When we combine these laboratory studies with the field data collected in the present study, we consider CSD with two independent loci and moderate sex allele diversity a more plausible explanation of low DMP in C. vestalis than sl-CSD with high sex allele diversity (Figure 2). Ultimately, genetic mapping of the sex loci and/or sequencing of the sex determination genes in C. vestalis is necessary to provide hard evidence for the existence of two-locus CSD.
Distinct genotypes were found among wasps collected in Luhzu township (fields 1 and 2), and in field 8 in Shanhua township, as indicated by clustering in Structure and supported by significant F ST values between these fields and all other fields (Figure 1 and Table 2). However, the largest cluster of wasps included individuals from all eight fields (Additional file 1: Table S4), and hierarchical partitioning of genetic variance indeed suggests that township and field are no good predictors of population structure in C. vestalis. Overall, we found a low level of genetic differentiation between C. vestalis from different fields and a high level of variation within fields (AMOVA, Table 3). Similar to our findings in C. vestalis, low levels of genetic differentiation between populations and a high level of genetic variation within populations were found for other parasitoids, including Campoletis sonorensis and Chelonus insularis , Cotesia glomerata , and Aphidius ervi . These studies suggest substantial gene flow over long distances, which counters genetic drift within populations. With respect to CSD, while few sex alleles may get lost locally due to genetic drift in such populations, introgression can maintain sex allele diversity at a larger scale. Interestingly, these studies were performed in an agricultural landscape with high frequencies of hosts. In comparison, different results were obtained in other studies on more natural systems of hosts and parasitoids. Global population differentiation of the parasitoid Neotypus melanocephalus, which parasitizes caterpillars of the endangered butterfly Maculinea nausithous that occur at low densities, was high , suggesting limited migration of parasitoids and an effect of small population sizes. Differences in population structure between parasitoid species are probably best explained by a combination of differences in migration rates (dispersal ability)  and population sizes, which in turn depend on host densities .
SNP marker development
We showed that SNP mining by de novo Illumina sequencing of the entire genome is an efficient way to develop a set of molecular markers for a non-model species. Due to its small genome size, a single lane of paired-end Illumina sequencing was sufficient to enable a near-complete (but fragmented) de novo assembly of the C. vestalis genome. Because we used pooled DNA of eight C. vestalis females, mapping reads against our reference genome allowed for detection of a large number of putative SNPs. Our method is an economical, fast, and qualitatively robust way of developing a genome-wide set of molecular markers for non-model species; for example compared to the rather laborious methods of microsatellite development . A drawback may be the specificity of SNPs for the population they were developed on [76,77], but this may be compensated by the large number of putative SNPs available for assay development. Preliminary tests on a small panel of C. vestalis wasps from different populations look promising: the large majority of the set of 98 markers amplified, and 10-45% of the markers was polymorphic in individuals collected from the field in the Netherlands or from laboratory crosses (J. G. de Boer; unpublished data).
In conclusion, we found a low level of sibling mating in a native population of C. vestalis in Taiwan, confirming our expectation of non-random mating. Even though C. vestalis is a solitary parasitoid, sibling matings may occur when siblings emerge from different host individuals that are spatially close to each other, for example on the same plant. This seems a likely explanation in C. vestalis, because its host, diamondback moth, may reach high densities, especially in agricultural landscapes. Comparisons of mating structure of closely related species parasitizing hosts with different spatial distribution patterns could be informative in this respect. Cotesia vestalis indeed produced diploid males at a low rate under natural conditions as we predicted from our laboratory finding of two-locus CSD. Diploid males in C. vestalis are developmentally viable and produce infertile triploid offspring . This is considered the most costly scenario of CSD in terms of effects on population growth rate and extinction risk in small populations in comparison to species in which diploid males die during development or evolved fertile diploid males [9,68,78]. We think that the relatively high fitness cost of diploid males, combined with non-random mating may have selected for evolution from sl-CSD to CSD with two independent loci in C. vestalis.
Availability of supporting data
The datasets supporting the results of this article are available in the Short Read Archive of NCBI (raw Illumina reads Accession SRP058413), and Dryad (the assembled scaffolds and contigs, SNP information and sequences, list of putative SNPs and SNP genotypes) http://doi:10.5061/dryad.jn14s, and is partly available within this article as the following additional tables and figure:
Table S1: Overview of Pl. xylostella and C. vestalis collected per field in Western Taiwan.
Table S2: Flow cytometric analysis of C. vestalis male ploidy levels per field.
Table S3: Observed and expected heterozygosity per locus for parasitoid wasps, and linkage between loci.
Table S4: Assignment of individual wasps to three genetic clusters inferred by Structure.
Table S5: Hierarchical partitioning of genetic variance (AMOVA) of C. vestalis genotyped at 98 SNP loci.
Table S6: Pairwise F ST for C. vestalis parasitoids based on genotypes at 98 SNP loci.
Table S7: Observed and expected heterozygosity and F IS of C. vestalis based on genotypes at 98 SNP loci.
Figure S1: Evanno’s DeltaK plot based on Structure suggesting 3 genetic clusters.
The authors acknowledge the World Vegetable Center in Shanhua, Taiwan, specifically Dr. Ramasamy for hosting us, and Mrs. Lin and Mrs. Yang of the insectary for their help during the fieldwork in Taiwan. Sander Bot is also thanked for his help in the field. We are also indebted to Remco van Poecke for his help in selecting SNPs, Willem van de Poll of Groningen University for his help with flow cytometry, Pieter van der Vlies of the UMCG for assistance in Illumina sequencing and KBioscience for setting up the KASPar assays. We much appreciate feedback on this manuscript from Axios Review. Funding was provided by the Netherlands Science Organisation (NWO grant 863.07.010 to JGdB), the Netherlands Genomics Initiative (NGI Zenith 935.11.04 to BAP) and the U.K. Genetics Society (Heredity field grant to LWB).
- Gamble T, Zarkower D. Sex determination. Curr Biol. 2012;22(8):R257–62.PubMedView ArticleGoogle Scholar
- Cutting A, Chue J, Smith CA. Just how conserved is vertebrate sex determination? Dev Dyn. 2013;242(4):380–7.PubMedView ArticleGoogle Scholar
- Trukhina AV, Lukina NA, Wackerow-Kouzova ND, Smirnov AF. The variety of vertebrate mechanisms of sex determination. Biomed Res Int. 2013:Article Number: 587460.Google Scholar
- Beukeboom LW, Perrin N. The evolution of sex determination. Oxford: Oxford University Press; 2014.View ArticleGoogle Scholar
- Bull JJ. Evolution of sex determining mechanisms. Menlo Park, CA: Benjamin/Cummings Publishing Company; 1983.Google Scholar
- Hodgkin J. Exploring the envelope: Systematic alteration in the sex-determination system of the nematode Caeraorhabditis elegans. Genetics. 2002;162(2):767–80.PubMed CentralPubMedGoogle Scholar
- Gempe T, Beye M. Function and evolution of sex determination mechanisms, genes and pathways in insects. BioEssays. 2011;33(1):52–60.PubMed CentralPubMedView ArticleGoogle Scholar
- Verhulst EC, van de Zande L, Beukeboom LW. Insect sex determination: it all evolves around transformer. Curr Opin Genet Dev. 2010;20(4):376–83.PubMedView ArticleGoogle Scholar
- Heimpel GE, de Boer JG. Sex determination in the Hymenoptera. Annu Rev Entomol. 2008;53:209–30.PubMedView ArticleGoogle Scholar
- Asplen MK, Whitfield JB, de Boer JG, Heimpel GE. Ancestral state reconstruction analysis of hymenopteran sex determination mechanisms. J Evol Biol. 2009;22(8):1762–9.PubMedView ArticleGoogle Scholar
- Cook JM, Crozier RH. Sex determination and population biology in the Hymenoptera. Trends Ecol Evol. 1995;10(7):281–6.PubMedView ArticleGoogle Scholar
- van Wilgenburg E, Driessen G, Beukeboom LW. Single locus complementary sex determination in Hymenoptera: an “unintelligent” design? Front Zool. 2006;3:1–15.PubMed CentralPubMedView ArticleGoogle Scholar
- Hunt GJ, Page RE. Linkage analysis of sex determination in the honey bee (Apis mellifera). Mol Gen Genet. 1994;244(5):512–8.PubMedView ArticleGoogle Scholar
- Hunt GJ, Page RE. Linkage map of the honey bee, Apis mellifera, based on RAPD markers. Genetics. 1995;139(3):1371–82.PubMed CentralPubMedGoogle Scholar
- Gadau J, Gerloff CU, Kruger N, Chan H, Schmid-Hempel P, Wille A, et al. A linkage analysis of sex determination in Bombus terrestris (L.) (Hymenoptera : Apidae). Heredity. 2001;87:234–42.PubMedView ArticleGoogle Scholar
- Beye M, Hasselmann M, Fondrk MK, Page RE, Omholt SW. The gene csd is the primary signal for sexual development in the honeybee and encodes an SR-type protein. Cell. 2003;114(4):419–29.PubMedView ArticleGoogle Scholar
- Beye M, Seelmann C, Gempe T, Hasselmann M, Vekemans X, Fondrk MK, et al. Gradual molecular evolution of a sex determination switch through incomplete penetrance of femaleness. Curr Biol. 2013;23(24):2559–64.PubMedView ArticleGoogle Scholar
- Hedrick PW, Gadau J, Page RE. Genetic sex determination and extinction. Trends Ecol Evol. 2006;21(2):55–7.PubMedView ArticleGoogle Scholar
- Ode PJ, Antolin MF, Strand MR. Brood-mate avoidance in the parasitic wasp Bracon hebetor Say. Anim Behav. 1995;49(5):1239–48.View ArticleGoogle Scholar
- Antolin MF, Strand MR. Mating system of Bracon hebetor (Hymenoptera, Braconidae). Ecol Entomol. 1992;17(1):1–7.View ArticleGoogle Scholar
- Whitehorn PR, Tinsley MC, Goulson D. Kin recognition and inbreeding reluctance in bumblebees. Apidologie. 2009;40(6):627–33.View ArticleGoogle Scholar
- Cowan DP, Stahlhut JK. Functionally reproductive diploid and haploid males in an inbreeding hymenopteran with complementary sex determination. Proc Natl Acad Sci USA. 2004;101(28):10374–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Elias J, Mazzi D, Dorn S. No need to discriminate? Reproductive diploid males in a parasitoid with complementary sex determination. Plos One. 2009;4(6):e6024.PubMed CentralPubMedView ArticleGoogle Scholar
- Verhulst EC, Beukeboom LW, van de Zande L. Maternal control of haplodiploid sex determination in the wasp Nasonia. Science. 2010;328(5978):620–3.PubMedView ArticleGoogle Scholar
- Verhulst EC, Lynch JA, Bopp D, Beukeboom LW, van de Zande L. A new component of the Nasonia sex determining cascade is maternally silenced and regulates transformer expression. Plos One. 2013;8(5):e63618.PubMed CentralPubMedView ArticleGoogle Scholar
- de Boer JG, Ode PJ, Rendahl AK, Vet LEM, Whitfield JB, Heimpel GE. Experimental support for multiple-locus complementary sex determination in the parasitoid Cotesia vestalis. Genetics. 2008;180(3):1525–35.PubMed CentralPubMedView ArticleGoogle Scholar
- Chou L. Notes on Apanteles (Hym: Bracon.) of Taiwan. J Agric Res China. 1979;28(4):299–310.Google Scholar
- Oatman ER. Introduced parasites and predators of arthropod pest and weeds: a world review. Washington DC: USDA; 1978.Google Scholar
- Talekar NS, Shelton AM. Biology, ecology, and management of the diamondback moth. Annu Rev Entomol. 1993;38:275–301.View ArticleGoogle Scholar
- Godfray HCJ. Parasitoids. Behavioral and evolutionary ecology. Princeton, NJ, USA: Princeton University Press; 1994.Google Scholar
- Antolin MF, Ode PJ, Heimpel GE, O'Hara RB, Strand MR. Population structure, mating system, and sex-determining allele diversity of the parasitoid wasp Habrobracon hebetor. Heredity. 2003;91(4):373–81.PubMedView ArticleGoogle Scholar
- Elias J, Dorn S, Mazzi D. Inbreeding in a natural population of the gregarious parasitoid wasp Cotesia glomerata. Mol Ecol. 2010;19(11):2336–45.PubMedView ArticleGoogle Scholar
- Zhou Y, Gu H, Dorn S. Single-locus sex determination in the parasitoid wasp Cotesia glomerata (Hymenoptera : Braconidae). Heredity. 2006;96(6):487–92.PubMedView ArticleGoogle Scholar
- de Boer JG, Ode PJ, Vet LEM, Whitfield JB, Heimpel GE. Diploid males sire triploid daughters and sons in the parasitoid Cotesia vestalis. Heredity. 2007;99(3):288–94.PubMedView ArticleGoogle Scholar
- de Boer JG, Kuijper B, Heimpel GE, Beukeboom LW. Sex determination meltdown upon biological control introduction of the parasitoid Cotesia rubecula? Evol Appl. 2012;5(5):444–54.PubMed CentralPubMedView ArticleGoogle Scholar
- Galbraith DW, Harkins KR, Maddox JM, Ayres NM, Sharma DP, Firoozabady E. Rapid flow cytometric analysis of the cell cycle in intact plant tissues. Science. 1983;220:1049–51.PubMedView ArticleGoogle Scholar
- Li RQ, Zhu HM, Ruan J, Qian WB, Fang XD, Shi ZB, et al. De novo assembly of human genomes with massively parallel short read sequencing. Genome Res. 2010;20(2):265–72.PubMed CentralPubMedView ArticleGoogle Scholar
- Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics. 2009;25:1754–60.PubMed CentralPubMedView ArticleGoogle Scholar
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics. 2009;25:2078–9.PubMed CentralPubMedView ArticleGoogle Scholar
- Kerstens HHD, Crooijmans R, Veenendaal A, Dibbits BW, Chin-A-Woeng TFC, den Dunnen JT, et al. Large scale single nucleotide polymorphism discovery in unsequenced genomes using second generation high throughput sequencing technology: applied to turkey. BMC Genomics. 2009;10:479.PubMed CentralPubMedView ArticleGoogle Scholar
- Kraus RHS, Kerstens HHD, Van Hooft P, Crooijmans R, Van der Poel JJ, Elmberg J, et al. Genome wide SNP discovery, analysis and evaluation in mallard (Anas platyrhynchos). BMC Genomics. 2011;12:150.PubMed CentralPubMedView ArticleGoogle Scholar
- Semagn K, Babu R, Hearne S, Olsen M. Single nucleotide polymorphism genotyping using Kompetitive Allele Specific PCR (KASP): overview of the technology and its application in crop improvement. Mol Breed. 2014;33(1):1–14.View ArticleGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.PubMed CentralPubMedGoogle Scholar
- Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics. 2003;164(4):1567–87.PubMed CentralPubMedGoogle Scholar
- Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8):2611–20.PubMedView ArticleGoogle Scholar
- Guo SW, Thompson EA. A Monte-Carlo method for combined segregation and linkage analysis. Am J Hum Genet. 1992;51(5):1111–26.PubMed CentralPubMedGoogle Scholar
- Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10(3):564–7.PubMedView ArticleGoogle Scholar
- Holm S. A simple sequentially rejective multiple test procedure. Scan J Stat. 1979;6:65–70.Google Scholar
- Wright S. The genetical structure of populations. Ann Eugen. 1951;15(4):323–54.PubMedGoogle Scholar
- Weir B, Cockerham C. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.View ArticleGoogle Scholar
- Pamilo P. Effect of inbreeding on genetic relatedness. Hereditas. 1985;103(2):195–200.PubMedView ArticleGoogle Scholar
- Paxton RJ, Thoren PA, Gyllenstrand N, Tengo J. Microsatellite DNA analysis reveals low diploid male production in a communal bee with inbreeding. Biol J Linn Soc. 2000;69(4):483–502.View ArticleGoogle Scholar
- Stahlhut JK, Cowan DP. Inbreeding in a natural population of Euodynerus foraminatus (Hymenoptera : Vespidae), a solitary wasp with single-locus complementary sex determination. Mol Ecol. 2004;13(3):631–8.PubMedView ArticleGoogle Scholar
- Adams J, Rothman ED, Kerr WE, Paulino ZL. Estimation of the number of sex alleles and queen matings from diploid male frequencies in a population of Apis mellifera. Genetics. 1977;86:583–96.PubMed CentralPubMedGoogle Scholar
- Van Bers N, Van Oers K, Kerstens H, Dibbits B, Crooijmans R, Visser M, et al. Genome-wide SNP detection in the great tit Parus major using high throughput sequencing. Mol Ecol. 2010;19:89–99.PubMedView ArticleGoogle Scholar
- Jonker RM, Zhang Q, Van Hooft P, Loonen M, Van der Jeugd HP, Crooijmans R, et al. The development of a genome wide SNP set for the Barnacle goose Branta leucopsis. Plos One. 2012;7(7):e38412.PubMed CentralPubMedView ArticleGoogle Scholar
- Vignal A, Milan D, SanCristobal M, Eggen A. A review on SNP and other types of molecular markers and their use in animal genetics. Genet Sel Evol. 2002;34(3):275–305.PubMed CentralPubMedView ArticleGoogle Scholar
- Cooper DN, Mort M, Stenson PD, Ball EV, Chuzhanova NA. Methylation-mediated deamination of 5-methylcytosine appears to give rise to mutations causing human inherited disease in CpNpG trinucleotides, as well as in CpG dinucleotides. Hum Genomics. 2010;4:406–10.PubMed CentralPubMedView ArticleGoogle Scholar
- Scarano E, Iaccarin M, Grippo P, Parisi E. Heterogeneity of thymine methyl group origin in DNA pyrimidine isostichs of developing sea urchin embryos. Proc Natl Acad Sci USA. 1967;57(5):1394–400.PubMed CentralPubMedView ArticleGoogle Scholar
- Milano I, Babbucci M, Panitz F, Ogden R, Nielsen RO, Taylor MI, et al. Novel tools for conservation genomics: Comparing two high-throughput approaches for SNP discovery in the transcriptome of the European Hake. Plos One. 2011;6(11):e28008.PubMed CentralPubMedView ArticleGoogle Scholar
- 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. 2010;327(5963):343–8.PubMedView ArticleGoogle Scholar
- Ruf D, Dorn S, Mazzi D. Unexpectedly low frequencies of diploid males in an inbreeding parasitoid with complementary sex determination. Biol J Linn Soc. 2013;108(1):79–86.View ArticleGoogle Scholar
- Ross KG, Vargo EL, Keller L, Trager JC. Effect of a founder event on variation in the genetic sex-determining system of the fire ant Solenopsis invicta. Genetics. 1993;135:843–54.PubMed CentralPubMedGoogle Scholar
- Tsuchida K, Kudo K, Ishiguro N. Genetic structure of an introduced paper wasp, Polistes chinensis antennalis (Hymenoptera, Vespidae) in New Zealand. Mol Ecol. 2014;23(16):4018–34.PubMedView ArticleGoogle Scholar
- Lechner S, Ferretti L, Schoning C, Kinuthia W, Willemsen D, Hasselmann M. Nucleotide variability at its limit? Insights into the number and evolutionary dynamics of the sex-determining specificities of the Honey Bee Apis mellifera. Mol Biol Evol. 2014;31(2):272–87.PubMed CentralPubMedView ArticleGoogle Scholar
- Hasselmann M, Gempe T, Schiott M, Nunes-Silva CG, Otte M, Beye M. Evidence for the evolutionary nascence of a novel sex determination pathway in honeybees. Nature. 2008;454(7203):519–U517.PubMedView ArticleGoogle Scholar
- Yokoyama S, Nei M. Population dynamics of sex-determining alleles in honey bees and self-incompatibility alleles in plants. Genetics. 1979;91:609–26.PubMed CentralPubMedGoogle Scholar
- Hein S, Poethke HJ, Dorn S. What stops the ‘diploid male vortex’? A simulation study for species with single locus complementary sex determination. Ecol Model. 2009;220:1663–9.View ArticleGoogle Scholar
- Thiel A, Weeda AC, de Boer JG, Hoffmeister TS. Genetic incompatibility drives mate choice in a parasitic wasp. Front Zool. 2013;10:43.PubMed CentralPubMedView ArticleGoogle Scholar
- Gu HN, Dorn S. Mating system and sex allocation in the gregarious parasitoid Cotesia glomerata. Anim Behav. 2003;66:259–64.View ArticleGoogle Scholar
- Jourdie V, Alvarez N, Molina-Ochoa J, Williams T, Bergvinson D, Benrey B, et al. Population genetic structure of two primary parasitoids of Spodoptera frugiperda (Lepidoptera), Chelonus insularis and Campoletis sonorensis (Hymenoptera): to what extent is the host plant important? Mol Ecol. 2010;19(10):2168–79.PubMedView ArticleGoogle Scholar
- Hufbauer RA, Bogdanowicz SM, Harrison RG. The population genetics of a biological control introduction: mitochondrial DNA and microsatellie variation in native and introduced populations ofAphidus ervi, a parisitoid wasp. Mol Ecol. 2004;13(2):337–48.PubMedView ArticleGoogle Scholar
- Anton C, Zeisset I, Musche M, Durka W, Boomsma JJ, Settele J. Population structure of a large blue butterfly and its specialist parasitoid in a fragmented landscape. Mol Ecol. 2007;16(18):3828–38.PubMedView ArticleGoogle Scholar
- Kankare M, van Nouhuys S, Gaggiotti O, Hanski I. Metapopulation genetic structure of two coexisting parasitoids of the Glanville fritillary butterfly. Oecologia. 2005;143(1):77–84.PubMedView ArticleGoogle Scholar
- Kraus RHS, von Holdt B, Cocchiararo B, Harms V, Bayerl H, Kühn R, et al. A SNP-based approach for rapid and cost-effective genetic wolf monitoring in Europe based on non-invasively collected samples. Mol Ecol Resour. 2015;15:295-305 doi:10.1111/1755-0998.12307.
- Bradbury IR, Hubert S, Higgins B, Bowman S, Paterson IG, Snelgrove PVR, et al. Evaluating SNP ascertainment bias and its impact on population assignment in Atlantic cod, Gadus morhua. Mol Ecol Resour. 2011;11:218–25.PubMedView ArticleGoogle Scholar
- Rosenblum EB, Novembre J. Ascertainment bias in spatially structured populations: a case study in the eastern fence lizard. J Hered. 2007;98(4):331–6.PubMedView ArticleGoogle Scholar
- Zayed A, Packer L. Complementary sex determination substantially increases extinction proneness of haplodiploid populations. Proc Natl Acad Sci USA. 2005;102(30):10742–6.Google Scholar
- de Boer, JG. Cotesiavestalis_females_Taiwan (study SRP058413). NCBI SRA. 2015. http://www.ncbi.nlm.nih.gov/sra
- de Boer JG, Groenen MAM, Pannebakker BA, Beukeboom LW, Kraus RHS. Data from: Population-level consequences of complementary sex determination in a solitary parasitoid. Dryad. 2015. http://dx.doi.org/10.5061/dryad.jn14s.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.