Introduction
Wasps in the family Cynipidae (Hymenoptera) form intimate associations with their plant hosts through the formation of galls – plant induced structures that provide protection and enhanced nutrition for developing larvae (Ronquist & Liljeblad, Reference Ronquist and Liljeblad2001; Stone et al., Reference Stone, Schonrogge, Atkinson, Bellido and Pujade-Villar2002). Morphological and phylogenetic revisions support the division of the Cynipidae into several distinct tribes, within which the oak gall wasps (Cynipini) form a monophyletic assemblage (Liljeblad & Ronquist, Reference Liljeblad and Ronquist1998; Rokas et al., Reference Rokas, Nylander, Ronquist and Stone2002, Reference Rokas, Melika, Abe, Nieves-Aldrey, Cook and Stone2003; Ronquist et al., Reference Ronquist, Nieves-Aldrey, Buffington, Liu, Liljeblad and Nylander2015). However, within the Cynipini, new genera and species are frequently being described (Pujade-Villar et al., Reference Pujade-Villar, Romero-Rangel, Chagoyan-Garcia, Equihua-Martinez, Estrada-Venegas and Melika2010; Medianero, Nieves-Aldrey & Melika, Reference Medianero, Nieves-Aldrey and Melika2011; Ide et al., Reference Ide, Wachi and Abe2012; Pujade-Villar et al., Reference Pujade-Villar, Hanson, Medina, Torres and Melika2012; Melika et al., Reference Melika, Tang, Sinclair, Yang, Lohse, Hearn, Nicholls and Stone2013; Tang et al., Reference Tang, Sinclair, Hearn, Yang, Stone, Nicholls, Schweger and Melika2016) and existing species and genera are often re-assigned making inferences into the biological characteristics of members of this group difficult (Drown & Brown, Reference Drown and Brown1998; Melika & Abrahamson, Reference Melika and Abrahamson2007). While the vast majority of the oak gall wasps persist at low population level densities, due in part to regulation by parasitoids and inquilines (Washburn & Cornell, Reference Washburn and Cornell1981; Moriya et al., Reference Moriya, Inoue, Otake, Shiga and Mabuchi1989; Stone et al., Reference Stone, Schonrogge, Atkinson, Bellido and Pujade-Villar2002), occasionally species are known to outbreak and cause host plant damage and even mortality (Eliason & Potter, Reference Eliason and Potter2000; Cooper & Rieske, Reference Cooper and Rieske2007; Pujade-Villar et al., Reference Pujade-Villar, Cibrian-Tovar, Barrera-Ruiz and Melika2014).
One such outbreaking species in the northeastern United States is the newly described black oak gall wasp, Zapatella davisae Buffington et al. (Buffington et al., Reference Buffington, Melika, Davis and Elkinton2016). This species has caused extensive tree damage and mortality to black oaks, Quercus velutina Lamarck (Fagales: Fagaceae), primarily in the coastal regions of Connecticut, Massachusetts, New York and Rhode Island, with defoliation greatest in Massachusetts on Cape Cod and Martha's Vineyard (Davis, Reference Davis2017). In these regions efforts to control outbreaking populations include injections of systemic pesticides into black oak trees including imidacloprid or emamectin benzoate (Davis, Reference Davis2017; Davis & Elkinton, Reference Davis and Elkinton2018). Unfortunately, basic information about the life history of Z. davisae, such as its geographic distribution are currently unknown (Davis, Reference Davis2017). It is also unclear whether Z. davisae reproduces sexually, asexually, or has strictly alternating sexual and asexual generations at different times of the year – a form of cyclical parthenogenesis that is unique to the Cynipidae (Atkinson et al., Reference Atkinson, McVean and Stone2002; Stone et al., Reference Stone, Schonrogge, Atkinson, Bellido and Pujade-Villar2002).
While Z. davisae was described as asexual based on the lack of observed males in sampled populations (Buffington et al., Reference Buffington, Melika, Davis and Elkinton2016), definitively identifying the mode of reproduction for a species can be a difficult process for several reasons (as reviewed in Schurko et al., Reference Schurko, Neiman and Logsdon2009). These include, but are not limited to (1) that the presence of males does not guarantee that sexual reproduction occurs (Schaefer et al., Reference Schaefer, Domes, Heethoff, Schneider, Schön, Norton, Scheu and Maraun2006); (2) that field observations of sex can be misleading because females in species that engage in intercourse do not necessarily utilize the sperm (Neiman & Lively, Reference Neiman and Lively2005) and (3) that many species are described as asexual due to the lack of males, but that cryptic sex can be common (Stone et al., Reference Stone, Atkinson, Rokas, Nieves-Aldrey, Melika, Acs, Csoka, Hayward, Bailey, Buckee and McVean2008). Due to these difficulties, molecular surveys can be used to provide additional insights into a species’ mode of reproduction (Schurko et al., Reference Schurko, Neiman and Logsdon2009). These methods include examinations for statistical deviations from the Hardy–Weinberg equilibrium (HWE; i.e., deviations from expected allele-frequencies for a randomly mating sexual population; Hardy Reference Hardy1908), and the detection of strong linkage disequilibrium (LD; i.e., the lack of statistical independence for alleles from different loci; Lewontin & Kojima Reference Lewontin and Kojima1960; Slatkin Reference Slatkin2008). Unfortunately, both measures of HWE and LD are known to show deviations from expectations for reasons other than lack of sexual reproduction (e.g., gene-flow, selection, expanding/contracting population sizes, etc.). Therefore, additional molecular measures have been proposed. One measure, termed the Meselson effect, seeks to quantify allelic differentiation at a locus under the assumption that for sexual lineages recombination will constrain the differentiation of alleles whereas in asexual lineages the lack of recombination will result in greater differentiation of alleles (Birky, Reference Birky1996). A second approach seeks to calculate statistical estimations for clonal population structure based on the genetic distance between individuals while accounting for possible scoring errors during molecular genotyping (Meirmans & Van Tienderen, Reference Meirmans and Van Tienderen2004). However, as with deviations from HWE and the presence of LD, these molecular methods are not always definitive, and combinations of field surveys and molecular surveys are therefore desirable to provide a more accurate assessment of reproductive mode than either approach independently.
Therefore, in an effort to aide future research of this important pest species, we utilized next-generation sequencing to develop genomic resources for Z. davisae with the specific goals of: (1) generating sequence data that can be applied to the studies of Z. davisae and other oak gall wasps, (2) generating microsatellite loci that can be used in future studies of the population genetics of Z. davisae, and (3) utilizing the polymorphic microsatellite loci to examine the mode of reproduction for populations of Z. davisae on Cape Cod and Nantucket.
Materials and methods
Ion-torrent library preparation and sequencing
Whole genomic DNA was extracted from 20 individual wasps collected on Cape Cod using the Qiagen PureGene Tissue kit (Qiagen Corp.) following the manufacturer's recommendations. DNA concentration and purity for each extract were then estimated using a NanoDrop 2000c spectrophotometer (ThermoFisher Scientific, Inc.). The two extracts with the highest concentration of DNA (ng μl−1) and the highest ratios of DNA to contaminants (i.e., 260/280 and 260/230 ratios) were then selected for next-generation sequencing library preparation. These included one individual each from Dennis, MA (C04) and West Harwich, MA (C12). Libraries were individually prepared for sequencing on the Ion Torrent™ PGM (ThermoFisher Scientific, Inc.) platform, using the NEXTflex™ DNA-seq kit (Bioo Scientific®, Inc.) with individual barcoding adaptors for each library (A5 for C04, and A7 for C12) by the Functional Genomics Laboratory at the University of California Berkeley. Library concentrations were then quantified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.), and pooled prior to sequencing at the DNA Analysis Facility at Science Hill, at Yale University on an Ion Torrent™ PGM 400 using the 318 v2 Chip kit.
Genomic assembly
Raw sequence reads were assembled into contigs using the de novo assembly program SOAP de novo 2 (Luo et al., Reference Luo, Liu, Xie, Li, Huang, Yuan, He, Chen, Pan, Liu, Tang, Wu, Zhang, Shi, Liu, Yu, Wang, Lu, Han, Cheung, Yiu, Peng, Zhu, Liu, Liao, Li, Yang, Wang and Lam2012) with the default settings running on the University of Massachusetts Green High Performance Computing Cluster. Assembled contigs were compared with the NCBI GenBank ‘nt’ sequence database (downloaded 11 July 2016) using the ‘blastn’ search algorithm (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990; Benson et al., Reference Benson, Cavanaugh, Clark, Karsch-Mizrachi, Lipman, Ostell and Sayers2013), with one best match sequence retained for each contig.
Microsatellite development
Contigs from the SOAP de novo 2 assembly were then used as reference sequences for microsatellite locus identification using the Phobos v. 3.3.11 (Mayer, Reference Mayer2010) plugin in Geneious v. 8.1.6 (Kearse et al., Reference Kearse, Moir, Wilson, Stones-Havas, Cheung, Sturrock, Buxton, Cooper, Markowitz, Duran, Thierer, Ashton, Mentjies and Drummond2012), iMSAT (Andersen & Mills, Reference Andersen and Mills2014), and MSATCOMMANDER (Faircloth, Reference Faircloth2008). For the analyses using Phobos, candidate microsatellite loci were filtered to only include exact repeats, a minimum unit length of two repeats, a maximum unit length of ten repeats, and to have no more than two successive ‘N’ calls in each contig. For the analyses using iMSAT, candidate microsatellite loci were filtered so that they would not include additional polymorphic sites within 150 bp of the target repeat. For the analyses using MSATCOMMANDER, candidate microsatellite loci were filtered to include only those loci with five or more repeat units. In addition, the raw sequence reads were assembled and used to identify microsatellite loci with the program QDD v. 3.1 (Meglecz et al., Reference Meglecz, Costedoat, Dubut, Gilles, Malausa, Pech and Martin2010, Reference Meglecz, Pech, Gilles, Dubut, Hingamp, Trilles, Grenier and Martin2014) using the default parameters and sequences with successfully designed primers were then compared with the RepBase database of transposable elements (Jurka, Reference Jurka2000) to avoid loci that might be duplicated in the genome.
Primers for candidate microsatellite loci were then modified so that the forward primer would include a 5′ M13 tail (5′-TCCCAGTCACGACGT-3′) and the reverse primer would include a 5′ poly-T tail (5′-GTTT-3′) allowing for the use of a single 6-FAM fluorescently labeled M13 primer during primer testing, following the method presented by Schuelke (Reference Schuelke2000). Candidate loci were then filtered to only include those from which primers could be designed with an optimal T m of 55–60°C, and priority for screening was given to those loci with the largest number of repeat units.
Microsatellite marker testing
To verify the utility of candidate microsatellite loci for Z. davisae population genetics research, primer pairs were tested using field collected-individuals from two populations in Massachusetts (Nantucket and Dennis). Individuals were collected by dissecting Q. velutina twigs with visible galls, and specimens identified as Z. davisae were then placed in 95% ethanol and stored at −20°C prior to molecular analyses. DNA was extracted from 13 field-collected individuals from Nantucket and 19 field-collected individuals from Dennis using the Qiagen DNeasy kit following the manufacture's recommendations. For each candidate locus, polymerase chain reactions (PCR) were performed using the following 25 µl reaction conditions: 18 µl water, 5 µl Promega GoTaq 5× buffer (Promega, Corp.), 0.5 µl dNTP (Promega), 0.1 µl 6-FAM M13 primer (100 µM), 0.05 µl candidate microsatellite forward primer (10 µM), 0.5 µl candidate microsatellite reverse primer (10 µM) and 0.2 µl Promega GoTaq (Promega). For all primer pairs, a touchdown protocol was used that cycled through annealing temperatures of 57–50°C (Andersen & Mills, Reference Andersen and Mills2014). PCR products were visualized via electrophoresis using 1.5% agarose gels, and for reactions where >50% of samples had visible bands, PCR products were genotyped at the DNA Analysis Facility on Science Hill at Yale University on a 3730xl 96-Capillary Genetics Analyzer (Applied Biosystems, Inc.). Genotype scores were determined in comparison to the LIZ-500 size standard through the microsatellite plugin in Geneious v. 8.1.6. Candidate loci from which all genotyped samples had the same fragment length were recorded as ‘monomorphic’, and candidate loci from which more than one fragment length was observed among genotyped samples (though with no more than two alleles for any one individual) were recorded as ‘polymorphic’.
Population genetic surveys for Z. davisae
For each polymorphic locus, the following standard population genetic statistics were recorded: the number of alleles observed among populations (Num), the effective number of alleles (Eff_num), the observed heterozygosity (H o), the expected heterozygosity within populations (H s), the expected total heterozygosity (H t), the total heterozygosity corrected for sampling from a limited number of populations (H′t), and the inbreeding coefficient (G IS) using the program GenoDive V 2.0b27 (Meirmans & Van Tienderen, Reference Meirmans and Van Tienderen2004). These statistics were also calculated for each of the two surveyed populations. For each locus and each population, HWE and evidence LD were tested using GenePop (Raymond & Rousset, Reference Raymond and Rousset1995; Rousset, Reference Rousset2008). To determine whether a specific microsatellite locus showed possible evidence of the Meselson effect (i.e., that two alleles from the same locus show a high level of differentiation), we compared the fragment length of alleles with each locus and calculated ‘coefficients of variation’ (CV) following Corral et al. (Reference Corral, Piwczynski, Sharbel, Schön, Martens and Dijk2009). Briefly, for each locus the standard deviation of the number of repeat motifs (σ) is divided by the mean number of repeats (μ) to estimate the CV. Lastly, to test whether populations of Z. davisae exhibit clonal population structure (i.e., the presence of more clones than would be expected under random-mating), we used the ‘Assign Clones’ algorithm implemented in GenoDive with the stepwise-mutation model. Following the recommendations provided in the GenoDive manual (available at: http://www.bentleydrummer.nl/software/GenoDiveManual.pdf), to account for scoring errors and/or possible mutations between clones, analyses were run using multiple genetic distances using the stepwise mutation model, of 0, 1, 10 and 20. For analyses at each distance, missing data were scored as ‘not counted’, and the presence of clonal diversity was assessed using the ‘Corrected Nei's diversity index’. Significance was determined based on 999 bootstrap permutations, with alleles randomized over individuals within populations (i.e., the default settings for GenoDive).
Results
Ion-torrent library preparation and sequencing
Sequencing of the individual barcoded libraries resulted in the generation of 5,032,726 sequence reads. Only 14 sequence reads could be assigned to sample C12, with the vast majority of sequence reads (4,136,285) being assigned to sample C04, while an additional 896,427 sequence reads were unassigned to either sample. Mean read lengths were 204, 244, and 243 bp to library C12, C04, and unassigned, respectively. After removal of low quality (<Q20) reads, 1,017,672,661 bp of Z. davisae associated genomic data was retained. Raw sequence reads for Z. davisae are available at the NCBI Short Read Archive (accession number SRR5570816) as part of BioProject PRJNA386803.
Genomic assembly
De novo alignment of raw sequence reads resulted in the assembly of 832 unique contigs. The average length of contigs was 380 ± 263 bp, with an average coverage of 11x ± 15x. Sixty-three contigs were of sufficient length (≥250 bp) to conduct ‘blastn’ searches. Of those 63 contigs, the vast majority (n = 52) were matched with published sequences from insects, and the majority of those (n = 44) matched published hymenopteran sequences. The full search results using ‘blastn’ including the taxonomic distributions for each sequence are presented in Supplementary table S1.
Microsatellite development
In total, 1252 candidate microsatellite loci were identified across the four utilized software programs. Phobos identified the most candidate loci (n = 671), followed by QDD (n = 550). Both iMSAT (n = 16) and MSATCOMMANDER (n = 15) identified far fewer candidate loci. The number of candidate loci identified by each program for different repeat motifs, as well as the average, maximum and minimum numbers of repeats are presented in table 1.
Table 1. Summary of candidate microsatellite loci identified by different repeat motifs using Phobos, iMSAT, MSATCOMMANDER, and QDD.
Count, the number of candidate loci for each motif; average repeats, the average length of repeat units; Max repeats and Min repeats, the maximum and minimum lengths of repeat units, respectively.
After filtering candidate microsatellite loci to those that had optimal primer sites, 120 candidate microsatellite loci were retained, from which 36 pairs of candidate loci were tested. These included 11 loci identified by Phobos, 12 loci identified by iMSAT, ten loci identified by MSATCOMMANDER, and four loci from QDD (note: one locus, 971, was identified by both iMSAT and MSATCOMMANDER). Results from QDD indicated the presence of many additional optimal candidate loci, however, because this program is computationally more intensive than the others examined, results were not available until later stages of primer testing. A complete list of all candidate microsatellite loci, including their repeat motifs, repeat lengths, fasta formatted primer sequences (in 5′ to 3′ orientation), fasta formatted sequence for the microsatellite repeat region (including flanking regions), and the program from which they were identified are provided in Supplementary table S2.
Of the 36 loci tested, 18 were monomorphic (five from Phobos, eight from iMSAT, five from MSATCOMMANDER), five were polymorphic (two from Phobos, two from iMSAT, one from QDD), 13 loci (three from Phobos, two from iMSAT, five from MSATCOMMANDER, three from QDD) either failed to amplify or had bands indicating the amplification of PCR products of the target fragment length but resulted in fragment lengths that could not be scored (generally the result of stutter or nonspecific amplification). Based on these results, iMSAT had the highest percentage of candidate loci that resulted in scoreable fragment lengths (83%), while QDD had the highest percentage of candidate loci that were polymorphic (25%).
Population genetic surveys for Z. davisae
Population genetic statistics were generated for the five polymorphic loci. Among individuals collected on Nantucket and in Dennis, three loci (0781, 1333, and 2121) were bi-allelic, one locus (0887) was tri-allelic, and one locus (1643) was hexa-allelic. A summary of the genetic diversity for each of the polymorphic microsatellite loci is presented in table 2. Populations from Nantucket and Dennis had similar average numbers of alleles (n = 2), as well as similar average heterozygosity (H o = 0.177 for Nantucket and H o = 0.198 for Dennis). A summary of genetic diversities for both populations is presented in table 3. Both populations displayed highly significant departures from the HWE (P < 0.0001 for both) while there was no evidence of LD within or across populations (P > 0.05 for all comparisons). Following Corral et al. (Reference Corral, Piwczynski, Sharbel, Schön, Martens and Dijk2009), estimates of CV were calculated for three of the five polymorphic microsatellite loci (table 4). These values ranged from 0 for locus 2121, to ~3 for locus 781 and locus 1643.
Table 2. Summary of allelic diversity for polymorphic microsatellite loci.
n, number of genotyped individuals; Num, number of alleles observed among populations; Eff_num, effective number of alleles; H o, observed heterozygosity, H s, heterozygosity within populations; H t, total heterozygosity; H′t, corrected total heterozygosity; G IS, inbreeding coefficient; HWE, probability of deviation from the Hardy–Weinberg equilibrium.
Table 3. Summary of allelic diversity for populations of Z. davisae.
n, number of genotyped individuals from each population; Num, number of alleles observed among populations; Eff_num, effective number of alleles; H o, observed heterozygosity, H s, heterozygosity within populations; H t, total heterozygosity; H′t, corrected total heterozygosity; G IS, inbreeding coefficient; HWE, probability of deviation from the Hardy–Weinberg equilibrium.
Table 4. CV for polymorphic microsatellite loci.
CV, coefficient of variation. Calculated by dividing the standard deviation of the number of repeat motifs for a locus (σ) by the mean number of repeats (μ).
Significant presence (P < 0.001) of clonal population structure was observed at genetic distances of 0, 1, and 10, but not at a distance of 20. Individuals could be assigned to one of seven unique clones with a distance threshold of 0, to five unique clones with a distance threshold of 4, three unique clones with a distance threshold of 6, two unique clones with a distance threshold of 10, and all individuals were assigned as the same clone with a distance threshold of 28.
Discussion
Oak gall wasps (tribe Cynipini) have fascinated researchers for generations with their capacity to manipulate their hosts into producing a diversity of plant-induced structures (Ronquist & Liljeblad, Reference Ronquist and Liljeblad2001; Cook et al., Reference Cook, Rokas, Pagel and Stone2002; Stone et al., Reference Stone, Schonrogge, Atkinson, Bellido and Pujade-Villar2002; Askew et al., Reference Askew, Melika, Pujade-Villar, Schonrogge, Stone and Nieves-Aldrey2013). Many of these species can even produce morphologically distinct galls on different plant tissues, and/or on the same host tissue but at different times in the year (Cook et al., Reference Cook, Rokas, Pagel and Stone2002). However, the study of these complex plant–insect interactions has been hampered by an overall lack of clarity as to the overall number of species and genera that make up the Cynipini – likely due to the vast number of species yet to be found and the taxonomic validity of existing species (see references in Melika & Buss, Reference Melika and Buss2002). While several recent morphological and phylogenetic studies have clarified some of the higher level relationships among the tribes in the Cynipidae (e.g., Liljeblad & Ronquist, Reference Liljeblad and Ronquist1998; Rokas et al., Reference Rokas, Nylander, Ronquist and Stone2002; Rokas et al., Reference Rokas, Melika, Abe, Nieves-Aldrey, Cook and Stone2003; Ronquist et al., Reference Ronquist, Nieves-Aldrey, Buffington, Liu, Liljeblad and Nylander2015), relationships within the Cynipini remain in flux as evident by the numerous revisions and re-assignments of species and genera, and the frequent identification of novel species (e.g., Buffington & Morita Reference Buffington and Morita2009; Pujade-Villar et al., Reference Pujade-Villar, Romero-Rangel, Chagoyan-Garcia, Equihua-Martinez, Estrada-Venegas and Melika2010; Medianero, Nieves-Aldrey & Melika Reference Medianero, Nieves-Aldrey and Melika2011; Ide, Wachi & Abe Reference Ide, Wachi and Abe2012; Pujade-Villar et al., Reference Pujade-Villar, Hanson, Medina, Torres and Melika2012; Melika et al., Reference Melika, Tang, Sinclair, Yang, Lohse, Hearn, Nicholls and Stone2013; Tang et al., Reference Tang, Sinclair, Hearn, Yang, Stone, Nicholls, Schweger and Melika2016).
Development of microsatellite loci
The lack of taxonomic clarity that currently defines species in the Cynipini presents challenges for the study and control of the recently described black oak gall wasp, Z. davisae, that has been causing extensive damage and tree mortality of black oaks, Q. velutina, on Cape Cod, Long Island, Martha's Vineyard, and Nantucket Island (Buffington et al., Reference Buffington, Melika, Davis and Elkinton2016). For example, little is known about its phylogenetic placement in the Cynipini, and/or its geographic distribution in the northeastern United States (Davis, Reference Davis2017). Recent DNA sequencing by Davis (Reference Davis2017) of the mitochondrial locus cytochrome oxidase I suggests that individuals of Z. davisae in the northeastern United States have multiple close matches (>97% similarity) with individuals from Florida (BOLD: BBHYH546-10), Nova Scotia (GenBank: KR804678.1, GenBank: KR782782.1), Ontario (BOLD: CNROF073-13, BOLD: CNTIF3344-15, BOLD: CNRGD596-15, GenBank: KR409819.1, GenBank: KR414957.1), and Tennessee (BOLD: GMGSM442-13). However, due to the limited number of available sequences in public databases, and the lack of resolution provided by a single mitochondrial locus, it is unclear whether populations in the northeastern United States are recently introduced from the southeastern United States (Davis (Reference Davis2017) reports that the highest % similarity matches were to individuals in Florida and Tennessee), or whether Z. davisae has a broad geographic range that includes much of eastern North America. Therefore, we hope to use the microsatellites developed in this study to sample populations of Z. davisae from throughout the eastern United States to determine whether outbreaks are the result of enemy release (which would be predicted if populations in the northeastern United States were recently introduced), or due to some other factor (if it has a broad geographic range that includes much of the eastern United States). In addition we hope to use these microsatellite markers to determine the relationship between Z. davisae and Callirhytis ceropteroides Bassett, a species previously reported to have caused black oak damage and mortality in Connecticut and New York State (Bassett, Reference Bassett1900; Pike et al., Reference Pike, Robison, Abrahamson, Ozaki, Yukawa, Ohgushi and Price2006), as understanding the biology and natural enemies of C. ceropteroides may help inform the control of Z. davisae as well.
Mode of reproduction
In addition to causing damage to black oaks, Z. davisae has been described as asexual based on the lack of observed males. Recent molecular work for oak gall wasps in the genus Andricus, however, found evidence for widespread cryptic sex among previously assumed asexual species (Stone et al., Reference Stone, Atkinson, Rokas, Nieves-Aldrey, Melika, Acs, Csoka, Hayward, Bailey, Buckee and McVean2008). Therefore, we were curious as to whether or not we could utilize molecular methods to determine the reproductive mode of Z. davisae. While these methods have their advantage of overcoming many of the biases inherent in field observations, definitive determination of a species’ sexual status can still be difficult (Ellegren, Reference Ellegren2004; Selkoe & Toonen, Reference Selkoe and Toonen2006; Bhargava & Fuentes, Reference Bhargava and Fuentes2010). Unfortunately, our efforts for combining field observations with comparisons to multiple genetic signatures (deviations from the HWE, presence of LD, presence of the Meselson effect, and statistical tests for clonal population structure) were inconclusive. Our field surveys collected only females, and both of our sampled populations showed highly significant (P < 0.0001) deviations from HWE and a significant presence of clonal population structure was observed at several genetic distances. However, deviations from HWE and the presence of clonal population structure could be the result of other factors, such as genetic bottlenecks following founder events and are therefore not necessarily indicative of asexual reproduction. In addition, our efforts to follow the approach of Corral et al. (Reference Corral, Piwczynski, Sharbel, Schön, Martens and Dijk2009) to determine the presence of the Meselson effect were inconclusive as we cannot compare measures with our polymorphic loci in reference to a known sexual population, as done by Corral et al. (Reference Corral, Piwczynski, Sharbel, Schön, Martens and Dijk2009). Lastly, we found no evidence for the presence of LD within or among populations of Z. davisae. While the lack of LD among Cynipini asexual populations is similar to that found by Stone et al. (Reference Stone, Atkinson, Rokas, Nieves-Aldrey, Melika, Acs, Csoka, Hayward, Bailey, Buckee and McVean2008), in total our results are inconclusive as to the sexual/asexual status of Z. davisae. However, we hope the developed microsatellites will be used in future studies that include populations of Z. davisae from throughout its geographic distribution and/or include sexually reproducing congeners to further explore the mode of reproduction for Z. davisae.
Conclusions
Here we provide a description of recently developed genomic resources that can be applied to the study of the biology and the control of Z. davisae, an important pest species of Q. velutina (black oaks) in the northeastern United States. These resources include over one billion base pairs of sequence data, over 900 assembled contigs, and 120 candidate microsatellite loci. We tested the utility of 36 of the microsatellite loci for Z. davisae population genetics research and found that half of them resulted in monomorphic genotype scores, while only five were polymorphic. Monomorphic loci are not normally useful for population genetic studies (see Ellegren, Reference Ellegren2004; Selkoe & Toonen, Reference Selkoe and Toonen2006), however; here we provide the primers and thermocycler conditions for these loci as supplemental information because our sampling of Z. davisae comes from two geographically proximal locations that may both be recent introductions and it is likely that when populations from greater geographic distances are observed that some of these microsatellite loci may also be found to be polymorphic. In addition, while the mutation dynamics of microsatellite loci can make their amplification difficult among closely related species (Primmer et al., Reference Primmer, Moller and Ellegren1996; Scribner & Pearce, Reference Scribner, Pearce and Baker2000; MacDonald et al., Reference MacDonald, Brookes, Edwards, Baker, Lockton and Loxdale2003; Selkoe & Toonen, Reference Selkoe and Toonen2006), numerous studies have successfully cross-amplified microsatellite loci, and thus some of the monomorphic loci presented here may be polymorphic in other species of gall wasps. Finally, while our genetic tests to determine the reproductive mode of Z. davisae were inconclusive, we hope that these microsatellite loci, the raw sequence reads, and the assembled contigs can be used for future studies of Z. davisae populations from across its distribution, in addition to studies of the population genetics of Cynipini species more broadly, and can contribute to the study of the evolution of gall formation and of reproductive modes in this fascinating group.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485318000858
Acknowledgements
The authors thank H. Broadley for her valuable feedback and for kindly ordering supplies, J. Boettner for building a functional molecular lab on a shoestring budget and for collecting last-minute samples, E. Mooshain for field assistance and laboratory identification work, D. Newman, M. Labbé, B. Normark, and two anonymous reviewers for their constructive comments and suggestions on an earlier version of this manuscript. Funding was provided by the University of Massachusetts Amherst Department of Environmental Conservation (MJD), the Arborjet Corp, Woburn, MA and the Woodbourne Arboretum, Melville, NY (MJD and JSE), the U.S. Forest Service Northern Research Station to NPH. A portion of this research represents part of Honors Thesis project of CPC, as part of the University of Massachusetts Amherst Commonwealth College Honors Program.
