BACKGROUND
western corn rootworm  is a beetle  that is native to north america and the dominant maize pest in the us corn belt.
it is estimated that each year western corn rootworm  costs us farmers at least  <dig> billion dollars through yield losses and treatment costs  <cit> .
wcr is native to north america, being first described in the western kansas in  <dig>  <cit> .
throughout the 20th century it spread eastward, colonizing the major maize production areas in the north central us   <cit> .
more recently it invaded europe  <cit> , making it now a global pest of maize production.
the impact of wcr on maize production has been exacerbated by the fact that wcr has proved difficult to control.
over the last few decades wcr has evolved resistance to crop rotation  <cit> , a range of chemical insecticides  and to transgenic maize expressing the cry3bb <dig> protein from bacillus thuringiensis   <cit> .
the genetics behind some wcr resistance mechanisms have been investigated .
for example, a single nucleotide polymorphism  in the gamma-aminobutyric acid  receptor is correlated with cylcodiene resistance  <cit> .
however, for many of the most pressing resistance traits the genetic basis remains unknown.
developing a genetic toolbox for wcr is an important step toward better understanding the genetic basis of insecticide resistance in wcr.
to this end, siegfried et al.
<cit>  produced  <dig> expressed sequence tags  from wcr midgut tissues and screened them for candidate receptors for bt.
among these ests, they reported one cadherin  and  <dig> cathespin-like fragments  that are potential candidates for resistance to bt and other toxins.
studies in other insect species have suggested that cytochrome p450s  are involved in resistance to several insecticides including methyl parathion, carbaryl, carbamate and organophosphates .
in wcr only  <dig> cyps have been characterized, and all belong to the cyp <dig> clade, which is not typically associated with insecticide resistance  <cit> .
furthermore, a variety of other candidate genes for insecticide or bt resistance are lacking for wcr, including esterases and other receptors  <cit> .
a more complete understanding of the diversity of these important resistance associated gene families can lead to a more rapid characterization of current and future resistance mechanisms.
in addition to resistance, questions remain regarding wcr population structure and patterns of gene flow.
previous studies for north american wcr populations have used several types of molecular markers , with all studies indicating large amounts of genetic diversity and a general lack of population structure across the corn belt .
these findings are consistent with the hypothesis that wcr has retained large amounts of genetic diversity during the species’ rapid eastward range expansion.
whether or not this pattern is shared across the entire wcr genome or if certain areas contain less genetic diversity is unknown.
as the amount gene flow among populations largely dictates the spread of resistance alleles, a better understanding of the genomic variation across populations will lead to a more informed assessment of resistance monitoring and mitigation, and prolong the life of various control options.
modern dna sequencing technologies have increased our ability to undertake large-scale sequencing projects for almost any organism, including those with complicated genomes such as wcr  <cit> .
therefore, our first objective for this study was to characterize the wcr transcriptome using publically available data as well as new data generated by  <dig> roche pyrosequencing.
we then mined this data to identify candidate genes in pathways predicted to be involved in insecticide resistance.
our second objective was to develop an expansive set of snps that can be used for more robust inferences of population structure, construction of linkage maps, and genomic scans of selection.
to produce these snps we generated genomic sequences from multiple laboratory and field populations, including five populations that showed reduced response to cry3bb <dig> in lab-based single plant bioassays.
furthermore, our ability to combine both transcriptomic and genomic data provide important information regarding genetic variation and level of selection among candidate resistance genes.
methods
sample collection, est sequencing and assembly
for  <dig> est sequencing wcr third-instar midgut tissues were isolated from  <dig> larvae from the waterman population  and stored in rnalater  prior to extraction.
total rna extraction was performed on four sets of  <dig> midguts using the rneasy mini kit .
approximately 500 μg of total rna were used for two rounds of polya selection using the micropolypurist kit .
following polya selection, cdnas were synthesized using the superscript ii cdna synthesis kit , using random hexamers for priming during first-strand synthesis.
finally, cdnas were purified using a qiaquick pcr purification column .table  <dig>
population details for all  <dig> populations sampled in this study
for each population the state of origin, year collected, source, diapause status, response to cry3bb <dig>  and intra-population genetic diversity  are given.
diapause status indicates whether or not the population undergoes diapause.
all non-diapause  populations derive this phenotype from the brookings, sd source.
all field sourced populations diapause.
response to cry3bb <dig> was determined by larval recovery data in single plant bioassays under laboratory conditions.
each population was reared on cry3bb <dig> expressing  and isoline  maize.
populations with at least 80% relative survival on treatment when compared to control were defined as low response, while medium response populations ranged between 20% and 80%, and high response populations were less than 20% relative survival.
missing data is indicated by na.
before 454-flx sequencing, emulsion pcr reactions were conducted following the manufacturers specifications .
in total,  <dig> , <dig>  <dig> est reads were generated with a mean length of  <dig>  bp, together totaling  <dig> ,184 bp of total sequence.
these  <dig> sequences were supplemented with  <dig>  sanger sequences from ncbi  adding an additional  <dig> ,393 bp.
the sanger sequences derive from a variety of tissues and life cycle stages, including many previously published by siegfried et al.
<cit> .
prior to assembly all  <dig> sequences were trimmed of their adaptors, while all sanger reads were trimmed by searching against the ncbi univec database of common cloning vectors.
<dig> and sanger sequences were merged and assembled with newbler  <dig> .
all  <dig> ests can be found on the ncbi sequence read archive .
transcriptome annotation and characterization
the wcr transcriptome was annotated using blast2go  <cit> .
blastx searches  were performed between wcr unigene sequences and the ncbi non-redundant  database.
following the mapping step, gene ontology  terms with e value < 10- <dig>  annotation cut-off >  <dig>  and go weight >  <dig> were used for annotation.
the “aqua” classification in categorizer  <cit>  was used to categorize the go terms into different go categories.
the go categories for wcr were compared to those from tribolium casteneum .
to find the pathways in which putative proteins of the wcr transcriptome are involved, analysis of kyoto encyclopedia of genes and genomes  was performed using blast2go  <cit> .
for comparative genomics, pairwise blastx searches  between wcr contig sequences and genomes of model insect species  were performed.
results of these blastx searches were also used to calculate the ortholog hit ratio .
the ortholog hit ratio was calculated by dividing the number of non-gap characters in the query  by the length of the subject   <cit> .
pooled population resequencing and snp detection
population genomic diversity was estimated by sampling  <dig> wcr populations ,  <dig> of which originated from field collections made in  <dig> and  <dig> .
where available, we also list each population’s phenotypic status for diapause and response to the cry3bb <dig> bt toxin based on single plant bioassays under laboratory conditions  <cit>  .
the remaining  <dig> populations originated from both diapause and non-diapause laboratory colonies.
among these populations the hopkinton and seneca research colonies are resistant to cry3bb <dig>  <cit> .figure  <dig>
map of wcr field locations.
collection locations are given for the  <dig> field collected populations used in this study.
<dig> collections are marked in blue, while  <dig> collections are marked in red.
for each population we randomly selected  <dig> individuals  from mixed sex populations for illumina sequencing .
all individuals within a population were combined prior to dna extraction .
dna was extracted using the e.z.n.a.
insect dna isolation kit .
illumina sequencing libraries were prepared using the nextera kit  and sequenced on an illumina hiseq  <dig> genome sequencer in a  <dig> × 100 bp paired-end configuration.
after sequencing, barcode adapters were removed and all reads were quality trimmed from both ends to a minimum phred score of  <dig>  all raw illumina sequences used in this study can be found on the ncbi sequence read archive .
paired-end illumina genomic shotgun sequences from all populations were aligned to the combined transcriptome assembly described above using bwa , with default settings.
aligning genomic sequences to a reference transcriptome creates mismatches at exon/intron boundaries because the introns are sampled by the genomic sequences but typically spliced out of the transcripts.
when bwa encounters these mismatches it issues a “soft-clip” flag  <cit> .
to avoid snp calling errors that might be associated with these mismatches, all reads were first filtered to remove any reads that had been soft-clipped.
we also removed all reads with a “mapq” alignment score <  <dig>  finally, snps were extracted from these alignments using the mpileup program from the samtools package , at all sites with at least 5x coverage from nucleotides with a phred sequence quality ≥  <dig>  all snp calls can be found in additional file  <dig>
to validate snps, we used the cleaved amplified polymorphic sequences  technique  <cit> .
potential caps markers were identified in silico by screening alternate snp haplotypes against the commercially available restriction enzyme database found in biopython .
we selected  <dig> snps that created a polymorphism in a recognition site for the enzyme rsai.
a region around each snp was amplified and digested in up to  <dig> wcr individuals .
primer sequences can be found in additional file  <dig>  pcr reactions were performed in a 20 μl mixture containing 1 μl  dna template, 10 μl 2x pcr promega master mix  and  <dig>  mm each primer.
amplified products for each sample were digested with 1 μl  rsai at 37°c for 2 hrs.
the resulting digest was visualized after electrophoresing in a 2% agarose gel stained with ethidium bromide.
estimating population genetic parameters
we estimated nucleotide diversity within populations for each unigene using a bias-corrected estimator  derived for pooled populations of individuals  <cit> .
this estimator requires a minimum count  for the inclusion of a particular snp variant in the calculation.
here we used b =  <dig> as our minimum.
we also disregarded any sites with greater than > 40x coverage to minimize the inclusion of paralogous sequences or undetected copy number variants.
within each population mean π was approximately  <dig> , however in each population we observed  <dig> to  <dig> outlier unigenes .
these outliers were removed prior to the calculation of summary statistics, as these loci are likely impacted by alignment issues given their extreme departure from the mean.
these calculations were done using a custom python script .
we computed pairwise fst between all populations for each unigene using a bias-corrected estimator for pooled population samples as implemented in popoolation <dig> , again using a minimum coverage of  <dig> and a maximum coverage of  <dig> for each snp locus.
window and step sizes were both set to  <dig>  forcing the program to estimate fst once for each unigene.
fst is a measure of genetic differentiation relative to the net genetic diversity of the populations sampled, which makes fst sensitive to fluctuations in net genetic diversity.
to complement fst, several authors have concluded that absolute metrics of genetic differentiation should also be employed  <cit> .
for this study we created a simple statistic  to estimate absolute genetic differentiation between populations using pooled sequence data.
this statistic is an estimator of nei’s dxy statistic  <cit> , though unlike dxy it uses estimated allele frequencies rather than allele counts, making it suitable for pooled sequence data.
to estimate πxy we first recorded a four element sequence profile b <dig>  b <dig>  b <dig>  b <dig>  for every nucleotide site within each pooled population, that contains the frequency of a, c, g, and t bases in the sample.
then for each unigene we calculated the average per-site expected heterozygosity among gametic unions between populations x and y , using the equation below.
where n is the number of nucleotide sites available for a unigene,  and  are the frequencies of base bj at nucleotide site i in populations x and y, respectively.
πxy requires only estimates of within population allele frequencies and not estimates of within population nucleotide diversity, unlike π and fst, and thus πxy does not require bias-correction for pooled populations.
we applied the πxy metric to calculate genetic distances between populations for principal coordinates analysis .
we also applied the πxy metric to characterize the genetic divergence between combined populations with low or high responsiveness to cry3bb <dig> or diapause and non-diapause populations .
to estimate πxy from combined populations we used the average of the nucleotide frequencies within each population to avoid biases that could arise due to uneven sequence coverage between populations.
as with π, we required a minimum coverage of  <dig> and a maximum coverage of  <dig> within each population for a site to be considered.
we also required a minimum count of  <dig> for the inclusion of a particular snp variant.
for each unigene we identified πxy and fst outliers by permuting population membership and calculating both metrics for permuted pairs.
in the case of the cry3bb <dig> response phenotype there are  possible permutations between high and low response populations, and all permutations were used to create the null distribution.
for diapause versus non-diapause populations there are  potential permutations, among which we randomly selected  <dig>  to create the null distribution.
in either case we selected the 99th percentile πxy and fst values from among these permuted populations as our threshold for significance, after removing all unigenes with fewer than  <dig> available nucleotide sites.
this threshold is arbitrary and lacks multiple test correction, though our goal is to produce a list of potential candidates that is inclusive, though it may include some false positives.
also, given the permutation tests outlined above, the most confidence one can have regarding significance is 1/ <dig> or 1/ <dig> , respectively.
our sampling procedure involved estimating population parameters from five individuals per population across approx.
<dig>  loci.
theory suggests that multi-locus estimates of population genetic parameters can be robust even when estimated from just five individuals  <cit> .
though this theory applies specifically to direct sequencing of individuals and not pooled population sequencing, we have conducted a power analysis  of our pooled sampling design, and it confirms that adequate estimates can be made from just five individuals, especially for statistics utilizing the grand mean across all loci.
in light of this, our grand mean genomic estimates  should be quite accurate, despite the small within population sample size.
due to the small within population sample size , we expect our single locus estimates of divergence  to have relatively modest statistical power.
these scans should be viewed as exploratory.
RESULTS
de novoassembly and annotation
the de novo assembly of the wcr transcriptome yielded  <dig>  high quality unigenes totaling  <dig> ,715 bp, in addition to  <dig>  singleton reads.
the length of unigenes varied from 99- <dig>  bp with an average of  <dig>  bp .
to determine the completeness of wcr transcriptome assembly, each unigene was compared to its putative ortholog from t. casteneum, d. melanogaster, and b. mori.
overall, 44-46% of the unigenes  had an ortholog hit ratio >  <dig>  and 54-58% were >  <dig>  .
on the assumption of conserved gene length between species, an ortholog hit ratio near zero indicates a poor assembly while values near one indicate a fully assembled transcriptome  <cit> .
following these criteria the current assembly for wcr appears to be moderately complete.figure  <dig>
transcriptome assembly of wcr.
a) length distribution of all  <dig>  transcripts; each transcript is placed along x-axis in the ascending order based on its length, and b) ortholog hit ratio for assembled contigs calculated after blastx searches to genomes of t. casteneum, d. melanogaster, and b. mori.
about 64%  of the wcr unigenes had one or more hits to protein sequences in the non-redundant  database at genbank while the remaining  had no match .
the majority of top blast hits for wcr unigenes were to insects , whereas a small proportion showed top hits to non-insect arthropods , other eukaryotes , and bacteria .
as expected, more than two-thirds of the top hits for the wcr unigenes were to coleopteran sequences in the nr database such as those of t. casteneum, dendroctonus ponderosae, chrysomela tremula, leptinotarsa decemlineata, phaedon cochleariae, tenebrio molitor, and batocera horsfieldi.
the wcr unigenes with no match to the nr database  may represent novel genes with functions not yet been identified, or they could be a product of assembly errors.
interestingly, the percentage of unmatched unigenes obtained in our study is considerably lower than those obtained for other insect transcriptomes .
an interproscan search  <cit>  of these unmatched wcr unigenes against the nr database revealed hits to the protein signature databases for  <dig>  out of  <dig>  , suggesting that many have homologs in other species that were undetected using blastx alone.
further,  <dig>  of the unmatched unigenes had an open reading frame  of at least 50% of the contig length and  <dig> had an orf of at least 90% of the contig length.
a pairwise comparison to the recently sequenced genome of mountain pine beetle  <cit>  did not improve the annotation significantly as only  <dig> out of  <dig>  unmatched wcr unigenes showed hits .
comparative genomics
using pairwise blastx searches to protein databases for four model insects, significant matches for a majority of wcr unigenes  were obtained.
the blastx search to the t. casteneum database showed the highest number of matches for wcr unigenes , followed by searches to n. vitripennis , b. mori , and d. melanogaster  databases .
several wcr unigenes  showed matches to all four databases.
however, there were a substantial number of wcr unigenes that matched uniquely to t. casteneum , n. vitripennis  and b. mori .
analyzing the pairwise comparisons, the highest level of agreement was observed between wcr and t. casteneum, likely reflecting the evolutionary relatedness between these two insect species and the completeness of the tribolium genome  <cit> .
functional annotation
using blast2go, about 46% of  wcr unigenes were annotated.
the observed gene ontology  terms for each domain  were widely distributed into different categories .
a comparison of percent mappings to each go category in wcr and t. casteneum revealed a highly similar distribution for both insect species.
the majority of contigs assigned to the ‘biological process’ domain were involved in cellular, regulatory, and developmental activities, while most of the contigs under the ‘molecular function’ domain were predicted to have catalytic, binding and transporter functions.
through the kegg-based pathway analysis using blast2go, a total of  <dig>  wcr unigenes were assigned to one or more of  <dig> pathways .
the majority of contig sequences were assigned to pathways for metabolism of organic compounds such as purine, pyrimidine, and glucose.
the information on predicted pathways coupled with the functional annotation of the wcr transcriptome provides useful information for future genetic studies in this insect.
we further characterized potential genes that may be involved in wcr bt resistance, specifically to cry3bb <dig>  <cit> , focusing on detoxification and degradation , esterases, and cathepsins).
we also identified several unigenes identified as putative receptors for bt toxins  such as cadherins, abc transporters, aminopeptidases and glycosyltransferases.
additional file  <dig> lists the wcr unigenes that significantly matched these sequences in other insects, mainly from tribolium.
the number of putative receptors totaled  <dig>  whereas the number of cathepsins, cyps, and esterases were much larger .
these totals mark a considerable expansion beyond what was previously known for wcr candidate resistance gene targets.
we further characterized two of the more important groups, the cyps and cathepsins, because of their known relevance to resistance.
cytochrome p450s  are a superfamily of enzymes which are ubiquitous among organisms and play an important role in metabolizing endogenous and exogenous substances.
in insects, cyp genes are classified into four clades: cyp <dig>  cyp <dig>  cyp <dig>  and mitochondrial p450s  <cit> .
despite their importance for insecticide resistance in wcr and other insects  <cit> , little information is known regarding the diversity or genetic variation of wcr-cyps.
in addition to blastx searches against the ncbi nr database, we identified the putative clades and families of wcr-cyps through a blastp comparison to a recently published tribolium cyp database  <cit> .
the majority of wcr-cyps had highest similarity to clade  <dig> , with a moderate number  showing similarity to clade  <dig> .
the expanded  <dig> sequencing combined with previous resources show a considerable diversity in the number of cathepsins in wcr.
we found a total of  <dig> unigenes that matched known cathepsins, including some with high similarity to those already described in wcr  <cit>  .
a previous est-based sequencing project  <cit> , revealed  <dig> cathepsin-like fragments , corresponding to two different classes, l and b. in our transcriptome data,  <dig> unigenes belonged to class b and  <dig> matched class l; we also identified  <dig> unigenes in the lesser known d and f classes .
cathepsins play a pivotal role in wcr digestion, and have been implicated in the evolution of rotation resistance.
normally, wcr feed and oviposit in corn, but rotation resistant wcr feed on soybean before oviposition.
when these soybean fields are rotated to corn the following year, the rotation resistant wcr eggs hatch and larvae feed on emerging corn plants.
curzi et al.
<cit>  found that wcr feeding on soybean leaves had higher levels of cathepsin-l like activity compared to wcr fed on corn.
this expanded set of cathepsin-l like sequences and targeted resequencing of wcr populations could help enable the development of a molecular diagnostic for rotation resistant wcr.
for all unigenes, including those that may be involved in wcr resistance to cry3bb <dig>  we pooled snp data from different populations with high or low response to cry3bb <dig>  and screened for genes that were πxy and fst outliers .
we identified  <dig> πxy outliers and  <dig> fst outliers .
only one unigene was found on both lists .
this unigene shares sequence similarity with a gene of unknown function from t. castaneum .
among all  <dig> πxy and fst outlier unigenes we found two on our lists of candidates for genes potentially involved in reduced response to cry3bb <dig>  these include one esterase  and one cyp , both of which were fst outliers.
all outliers should be considered preliminary and exploratory as our tests are low-powered given the relatively small within population sample sizes, and because our strategy assumes that populations with a low response to cry3bb <dig> share resistance alleles, though there is no evidence supporting this assumption.
nonetheless, these loci may represent a promising pool of potential candidates for cry3bb <dig> resistance.figure  <dig>
calculated
π
xy
and
f
st
values between populations with high and low response to cry3bb <dig> versus permuted 99
th
percentile values.
each point represents a unigene’s calculated π
xy
or f
st
value  between pooled high and low response to cry3bb <dig> populations and the 99th percentile value from all possible permutations of these populations .
red dots indicate calculated π
xy or f
st values > 99th percentile permutation value , while black dots represent calculated π
xy or f
st values < 99th percentile permutation value.
the blue line represents the boundary between outliers and non-outliers.
wcr snp development and validation
the field collected populations  included a wide geography across the us corn belt, while the laboratory populations  included several common research strains.
each population consisted of  <dig> individuals  that were pooled and sequenced to approximately 20x-30x coverage using whole genome shotgun sequencing.
these sequences were aligned to the reference transcriptome for snp discovery and estimation of population genetic parameters.
in total we found  <dig>  snps segregating among the  <dig> haploid chromosome sets sampled.
snp discovery saturates quickly among the wcr populations we sampled .
after sampling  <dig> populations, few new snps are detected among the remaining populations.
those new snps that are found by adding more than  <dig> populations tend to be rare, often occurring in just one population.
thus, within our unigene collection, it is possible that our data has revealed most of the common snps in corn belt wcr.figure  <dig>
rarefaction curve of snp discovery among the  <dig> wcr populations sampled.
starting with the first population sampled and continuing until the 26th we record the number of new snps each population contributes.
the slope of the rarefaction curve is dependent on population input order, which in this case is arbitrary.
to account for this, we selected  <dig> random input orders and computed the mean number of new snps discovered  and the upper and lower 95% confidence intervals .
to address the validity of our snp predictions and their utility for genotyping we identified  <dig> that were predicted to create a polymorphic restriction enzyme site and used this to genotype individuals from the inbred and randomly mated brookings population using caps .
for  <dig> of the  <dig> snps the caps genotypes showed evidence of both allelic versions of the predicted snp .
the remaining  <dig> snps were homozygous for one allele among all individuals tested, which could reflect inaccurate snp predictions or simply the absence of one allele among the populations tested.
in any case, these results demonstrate that the predicted snps are reliable and can be readily utilized in genotyping assays.
population genomics of western corn rootworm
using the snps extracted from pooled population sequencing we estimated genetic diversity within and between wcr populations.
to estimate within population nucleotide diversity  we used a bias-corrected estimator  <cit> , suitable for pooled population sequences.
for between population comparisons we created a simple statistic , which is an estimator of the dxy statistic  <cit> .
we focused our analyses on the  <dig>  wcr unigenes with a hit to a protein in the ncbi nr database .
for these unigenes, the mean π within each population can be found in table  <dig>  nucleotide diversity within populations ranged from  <dig> - <dig> , with a mean of  <dig> .
as expected, the inbred brookings colony  showed the lowest mean π, and the highest mean π values were found among field collected populations.
our extensive snp dataset reveals several interesting aspects of genetic variation in wcr.
first, the inbred population had experienced  <dig> generations of full-sib mating, which is theoretically expected to fix approximately 41% of the sites segregating among the gametes of the founding sib-mated pair , yet we observed only a slight reduction in π within the inbred brookings population  compared to the randomly mated brookings population  from which it was derived.
this suggests that full-sib mating has been fairly ineffective at removing variation.
second, there appears to be slightly less nucleotide diversity within the  <dig> collections, when compared to the  <dig> collections .
this could reflect annual variation in population sizes across the corn belt, or perhaps some unnoticed deviation in our sampling procedures between these two years.
finally, there is little evidence for a significant reduction of genetic diversity among lab populations, suggesting that maintenance in artificial environments has not caused an erosion of genetic diversity.
this finding mirrors observations by kim et al.
<cit> , who compared a laboratory population derived from the non-diapause brookings population to field populations using microsatellite markers.
among the  <dig> pairwise comparisons between populations, πxy ranged between  <dig>  and  <dig>  with a mean of  <dig> .
from this we observed that the levels of within population nucleotide diversity  are of a similar magnitude as between population divergence .
this indicates that each population maintains considerable genetic variation, and that there is little evidence for substantial genetic differentiation between populations .
pairwise πxy values were also used to generate a principal coordinates analysis among all populations.
diapause status appears to be the most important factor in population differentiation .
the natural state for wcr is diapause, and the trait shared by all  <dig> non-diapause populations was inherited from a single source  <cit> .
this non-diapause trait was developed after several generations of mass selection, and because of its complex inheritance, it has been suggested that it is likely a complex polygenic trait  <cit> .
consistent with this hypothesized mode of inheritance and strong artificial selection for its maintenance, we find  <dig> unigenes that are πxy outliers when comparing all diapause versus all non-diapause populations .
these outlier unigenes may include causal alleles for the non-diapause phenotype, but likely also include linked loci that have been differentiated between diapause and non-diapause populations by genetic hitchhiking.figure  <dig>
principal coordinate analysis using genetic distances between  <dig> wcr populations.
diapause populations are colored in black and non-diapause population are red.
populations with a reduced response to cry3bb <dig> are denoted with open squares, while populations with unknown or low to intermediate resistance are denoted with filled circles.
to more formally address population differentiation we estimated fst pairwise between all populations for all  <dig>  unigenes with a hit to the ncbi nr database.
mean fst between all populations across all unigenes was  <dig>  , and only  <dig> % of all unigenes had fst values >  <dig>  .
the overall pattern of low fst indicates little genetic differentiation among wcr populations, even across distances of approximately  <dig>  km.
these findings are consistent with earlier work in wcr, which also found low levels of genetic differentiation between field collected populations .
they are also consistent with our observation that new snp discovery saturates fairly quickly as new corn belt populations are sampled .
among all pairwise population contrasts, the highest fst values were found between field and lab populations, with nearly all of the most extreme values between a field population and the inbred brookings population .
this observation is expected, as inbreeding impacts allele frequencies and will tend to increase fst.
interestingly, the smallest mean pairwise fst estimate was between the onslow, ia and trent, sd populations .
despite being 508 km apart, these two field collected populations have both demonstrated reduced response to cry3b <dig> in a lab-based single plant bioassays .
it is premature to conclude that these populations derive from a single cry3bb <dig> resistant source, though this intriguing possibility warrants further research.
in the last 75 years wcr have expanded from their native range in the us southwest into the corn belt and further eastward, recently arriving at the eastern seaboard  <cit> .
given this recent and rapid range expansion, it is not surprising that our results and previous efforts  reveal limited genetic differentiation among corn belt wcr populations, even across distances of approximately  <dig>  km.
theoretical treatments of the genetic properties of range expansion typically find that populations near the edges of the range have reduced genetic diversity and increased fst when compared to populations near the center of the range, and that both of these features reset after colonization at a rate that is proportional to the migration rate between populations .
among the corn belt wcr populations, we find high genetic variation and low fst values among the most recently colonized easterly populations sampled .
these populations are now over  <dig>  km west of the eastern edge of the species range, and our genetic observations are consistent with theoretical predictions, in that these populations are no longer in their colonizing phase and are likely gaining new genetic diversity through gene flow with more westerly populations.
this pattern may also suggest that migration within the corn belt is fairly high; enough to spread considerable genetic diversity into populations founded in the last 40 years.
another interpretation could be that these two eastern populations were formed by a large and diverse set of founders, and that this diversity has been simply maintained since founding.
without historical samples it may be difficult to differentiate between these two hypotheses, though future studies could sample recently colonized populations in the eastern us to estimate typical founding population sizes.
we do observe a modest – albeit significant – trend of increasing genetic isolation as a function of geographic distance among field collected populations .
thus, there may be a subtle signature of the range expansion in the form of genetic isolation by distance.
moreover, the significance of the isolation by distance is strongly influenced by the inclusion of the most easterly population .
if we remove this population, the correlation between genetic and physical distance falls to  <dig>  and is no longer significant .
in contrast, if we remove all  <dig> populations from eastern colorado and western kansas  the correlation between genetic and physical distance is still significant .
the importance of this eastern population in detecting isolation by distance is consistent with incomplete homogenization following range expansion, however, because this result is driven by a single population we cannot exclude the hypothesis this population is in some way anomalous.
moreover, because recently formed corn belt wcr populations may not be at migration-drift equilibrium we cannot determine if the weak pattern of genetic isolation by distance is a product of genetic drift overcoming gene flow or if it is a hold-over from population founder effects during the species eastward invasion.
based on theory and these preliminary results we predict that sampling populations further east will tend to increase the extent of detectable of isolation by distance in wcr, and may reveal important estimates of genetic diversity in recently colonized areas.figure  <dig>
genetic isolation among wcr populations as a function of geographic distance.
following rousset  <cit> , pairwise estimates of f
st were linearized ; y-axis) and regressed onto the natural log of the physical distance between populations .
the linear regression line is shown in red along with the pearson’s correlation coefficient.
in a review article, mcdonald and linde  <cit>  posit that knowing information about the population genetics of a pest organism can lead to a better understanding of its evolutionary potential to evade control strategies.
they proposed that species with high levels of genetic diversity and high gene flow between populations have the greatest evolutionary potential to overcome control options through darwinian evolution.
in this study, and in earlier studies , we find evidence to support the hypothesis that wcr has both high levels of genetic diversity and high gene flow between populations.
for example, assuming wcr has a mutation rate  between 10- <dig> or 10-9 bp- <dig> generation- <dig> , the observed values of π imply that wcr field populations maintains approximate effective populations sizes  on the order of  <dig>  to  <dig> million individuals .
this level of ne gives little indication of recurrent population bottlenecks in wcr, which gives wcr populations the ability to generate and maintain many individuals with novel mutations.
this creates a large pool of standing genetic variation that may in turn increase the probability of evolving resistance to control strategies.
moreover, we show very little genetic differentiation between field populations, and that genetic variation has spread quickly eastward following range expansion, both of which are consistent with high rates of gene flow between populations.
thus, if resistance does evolve, there are few barriers preventing it from spreading between populations.
together these two features of wcr population genetics are hypothesized to lead to a high evolutionary potential to evade controls  <cit> , and may help to explain why wcr has a history of evolving resistance to insecticides and agricultural practices  <cit> .
CONCLUSIONS
wcr is one of the most significant agricultural pests in north american, yet at the time of publication wcr has a modest number of nucleotide sequences are available in genbank.
improved genetic and genomic resources for wcr are crucial tools in the ongoing effort to control this pest.
here we produce a large transcriptome assembly for wcr and complement these unigenes with diversity data obtained by whole genome population resequencing.
our wcr transcriptome assembly contains approximately  <dig>  unigenes, including several gene families that have been implicated in insecticide resistance in other species.
we identify a large pool of snps among  <dig> wcr populations, and show that these snps can be readily adapted for use as genotyping markers.
we also use these snps to scan for outliers among candidate bt resistance genes and to understand how population processes have shaped genetic variation in this species.
we find two potential candidate bt resistance genes that show strong patterns of molecular differentiation between high and low response to cry3bb <dig> populations.
we also expand on past population genetic studies in wcr and show that this species has significant genetic diversity, little population structure, and a high evolutionary potential to resist control strategies.
this study will provide a solid foundation for future research on the molecular genetics of wcr.
availability of supporting data
<dig> ests and illumina pooled population sequences used in this study can be found on the ncbi sequence read archive under the accessions srs <dig> and srp <dig>  respectively.
electronic supplementary material
additional file 1: genbank ids for public wcr est sequences.
this spreadsheet contains the genbank ids for the public wcr sequences incorporated our transcriptome assembly.
additional file 2: population level wcr snp genotypes.
a comma separated values file  containing the genotypic status of all  <dig> populations sampled at  <dig>  snp sites.
missing genotypes are denoted by a dash.
the sequence of each unigenes is listed on first usage.
this file has been compressed with the bzip <dig> compression software.
additional file 3: snp genotyping validation study.
this spreadsheet contains the caps snp validation information, including primer sequences, and observed restriction patterns and inferred genotypes.
additional file 4:
power analysis of population genetic sampling procedure.
additional file 5: blast2go annotations.
this spreadsheet contains annotations for  <dig>  wcr unigenes as identified by blast2go.
additional file 6: comparative genomics of wcr.
a venn diagram showing the number of wcr contigs with significant matches  to genomes of t. casteneum, d. melanogaster, n. vitripennis, and b. mori.
the significant matches  were calculated after pairwise comparisons  to each individual genome.
additional file 7: distribution and comparison of go categories.
vertical bars indicate the distribution of wcr  and red flower beetle  go term mappings that belong to each of the three top-level go categories .
additional file 8: kegg pathway annotations.
this file contains two spreadsheets.
the “detailed” spreadsheet contains putative kegg pathways assignments to unigenes in the wcr transcriptome, while the “summary” spreadsheet lists a summary of the counts of all pathways identified.
additional file 9: identification of putative candidate genes for cry3bb <dig> resistance.
this file contains a spreadsheet for each of several putative gene families that may be candidates for resistance to bt toxins, including the unigenes name and descriptions, and results from π
xyand f
st outlier scans between cry3bb <dig> resistant and susceptible populations.
the last two spreadsheets list the results of π
xyand f
st outlier scans between populations with high or low response to cry3bb <dig> for all available unigenes.
additional file 10: characterization of wcr-cyp genes.
blastp result comparing wcr-cyp unigenes to a tribolium cyp data set  <cit> .
additional file 11: identification of candidate genes for non-diapause trait.
this spreadsheet contains the results of π
xyoutlier scans between diapause and non-diapause populations for all available unigenes.
along with the scan results, unigenes name and descriptions are also given.
additional file 12: plot of calculated
π
xy
values between diapause and non-diapause populations and permuted 99
th
percentile values.
each point represents a unigenes calculated π
xy value and the 99th percentile value from  <dig>  permutations of the diapause and non-diapause population labels.
red dots indicate calculated π
xy > 99th percentile permutation value , while black dots represent calculated π
xy < 99th percentile permutation value.
the blue line represents the boundary between outliers and non-outliers.
additional file 13: a histogram of pairwise
f
st
between wcr populations.
this histogram represents f
st estimates for all genes among all pairwise comparisons of the  <dig> wcr populations.
competing interests
lef, rak, m chen, m carroll, bsg, rf and tc are employed by monsanto company and some of the authors are also shareholders in monsanto company, though individually or collectively in no way represent any controlling interest in the affairs of monsanto company.
results described in the manuscript are the product of work financed by monsanto company, and the article processing charge will be paid by monsanto company.
the authors have an obligation under an employment contract with monsanto company to assign the rights of any invention made in the course of their employment with monsanto company to the company.
monsanto company has filed one or more patent applications disclosing and claiming all or some portion of the results of the work described in the manuscript.
the authors are not known to have any financial competing interests in addition to those listed above.
authors’ contributions
lef, rb, rak, m chen, rf, tc, bsg, and apm conceived the study and contributed to experimental design.
lef, rb and apm analyzed the data and drafted the manuscript.
m chen and m carroll contributed materials.
all authors read and approved the final manuscript.
results described in the manuscript are the product of work financed by monsanto company.
rb and apm received additional support from the ohio agricultural research and development center, the ohio state university.
we thank aaron gassmann at iowa state university, bruce hibbard at usda ars, columbia, mo, and b. wade french at the usda north central agricultural research laboratory for sharing wcr colonies.
we are grateful to john cantwell, douglas jones, and leanna trail for their assistance in wcr field collections, oliver ilagan, kaylee kuban, cara vazquez for colony maintenance, and nay thane for producing the transcriptome assembly.
finally, we thank john ledeaux and peter morrell for helpful comments and discussions.

 western corn rootworm is a beetle that is native to north america and the dominant maize pest in the us corn belt. it is estimated that each year western corn rootworm costs us farmers at least billion dollars through yield losses and treatment costs. wcr is native to north america, being first described in the western kansas in. throughout the 20th century it spread eastward, colonizing the major maize production areas in the north central us. more recently it invaded europe, making it now a global pest of maize production. the impact of wcr on maize production has been exacerbated by the fact that wcr has proved difficult to control. over the last few decades wcr has evolved resistance to crop rotation, a range of chemical insecticides and to transgenic maize expressing the cry3bb protein from bacillus thuringiensis. the genetics behind some wcr resistance mechanisms have been investigated. for example, a single nucleotide polymorphism in the gamma-aminobutyric acid receptor is correlated with cylcodiene resistance. however, for many of the most pressing resistance traits the genetic basis remains unknown. developing a genetic toolbox for wcr is an important step toward better understanding the genetic basis of insecticide resistance in wcr. to this end, siegfried et al. produced expressed sequence tags from wcr midgut tissues and screened them for candidate receptors for bt. among these ests, they reported one cadherin and cathespin-like fragments that are potential candidates for resistance to bt and other toxins. studies in other insect species have suggested that cytochrome p450s are involved in resistance to several insecticides including methyl parathion, carbaryl, carbamate and organophosphates. in wcr only cyps have been characterized, and all belong to the cyp clade, which is not typically associated with insecticide resistance. furthermore, a variety of other candidate genes for insecticide or bt resistance are lacking for wcr, including esterases and other receptors. a more complete understanding of the diversity of these important resistance associated gene families can lead to a more rapid characterization of current and future resistance mechanisms. in addition to resistance, questions remain regarding wcr population structure and patterns of gene flow. previous studies for north american wcr populations have used several types of molecular markers, with all studies indicating large amounts of genetic diversity and a general lack of population structure across the corn belt. these findings are consistent with the hypothesis that wcr has retained large amounts of genetic diversity during the species’ rapid eastward range expansion. whether or not this pattern is shared across the entire wcr genome or if certain areas contain less genetic diversity is unknown. as the amount gene flow among populations largely dictates the spread of resistance alleles, a better understanding of the genomic variation across populations will lead to a more informed assessment of resistance monitoring and mitigation, and prolong the life of various control options. modern dna sequencing technologies have increased our ability to undertake large-scale sequencing projects for almost any organism, including those with complicated genomes such as wcr. therefore, our first objective for this study was to characterize the wcr transcriptome using publically available data as well as new data generated by roche pyrosequencing. we then mined this data to identify candidate genes in pathways predicted to be involved in insecticide resistance. our second objective was to develop an expansive set of snps that can be used for more robust inferences of population structure, construction of linkage maps, and genomic scans of selection. to produce these snps we generated genomic sequences from multiple laboratory and field populations, including five populations that showed reduced response to cry3bb in lab-based single plant bioassays. furthermore, our ability to combine both transcriptomic and genomic data provide important information regarding genetic variation and level of selection among candidate resistance genes. methods sample collection, est sequencing and assembly for est sequencing wcr third-instar midgut tissues were isolated from larvae from the waterman population and stored in rnalater prior to extraction. total rna extraction was performed on four sets of midguts using the rneasy mini kit. approximately 500 μg of total rna were used for two rounds of polya selection using the micropolypurist kit. following polya selection, cdnas were synthesized using the superscript ii cdna synthesis kit, using random hexamers for priming during first-strand synthesis. finally, cdnas were purified using a qiaquick pcr purification column.table population details for all populations sampled in this study for each population the state of origin, year collected, source, diapause status, response to cry3bb and intra-population genetic diversity are given. diapause status indicates whether or not the population undergoes diapause. all non-diapause populations derive this phenotype from the brookings, sd source. all field sourced populations diapause. response to cry3bb was determined by larval recovery data in single plant bioassays under laboratory conditions. each population was reared on cry3bb expressing and isoline maize. populations with at least 80% relative survival on treatment when compared to control were defined as low response, while medium response populations ranged between 20% and 80%, and high response populations were less than 20% relative survival. missing data is indicated by na. before 454-flx sequencing, emulsion pcr reactions were conducted following the manufacturers specifications. in total,, est reads were generated with a mean length of  bp, together totaling,184 bp of total sequence. these sequences were supplemented with sanger sequences from ncbi adding an additional,393 bp. the sanger sequences derive from a variety of tissues and life cycle stages, including many previously published by siegfried et al.. prior to assembly all sequences were trimmed of their adaptors, while all sanger reads were trimmed by searching against the ncbi univec database of common cloning vectors. and sanger sequences were merged and assembled with newbler. all ests can be found on the ncbi sequence read archive. transcriptome annotation and characterization the wcr transcriptome was annotated using blast2go. blastx searches were performed between wcr unigene sequences and the ncbi non-redundant database. following the mapping step, gene ontology terms with e value < 10- annotation cut-off >  and go weight >  were used for annotation. the “aqua” classification in categorizer was used to categorize the go terms into different go categories. the go categories for wcr were compared to those from tribolium casteneum. to find the pathways in which putative proteins of the wcr transcriptome are involved, analysis of kyoto encyclopedia of genes and genomes was performed using blast2go. for comparative genomics, pairwise blastx searches between wcr contig sequences and genomes of model insect species were performed. results of these blastx searches were also used to calculate the ortholog hit ratio. the ortholog hit ratio was calculated by dividing the number of non-gap characters in the query by the length of the subject. pooled population resequencing and snp detection population genomic diversity was estimated by sampling wcr populations, of which originated from field collections made in and. where available, we also list each population’s phenotypic status for diapause and response to the cry3bb bt toxin based on single plant bioassays under laboratory conditions. the remaining populations originated from both diapause and non-diapause laboratory colonies. among these populations the hopkinton and seneca research colonies are resistant to cry3bb.figure map of wcr field locations. collection locations are given for the field collected populations used in this study. collections are marked in blue, while collections are marked in red. for each population we randomly selected individuals from mixed sex populations for illumina sequencing. all individuals within a population were combined prior to dna extraction. dna was extracted using the e.z.n.a. insect dna isolation kit. illumina sequencing libraries were prepared using the nextera kit and sequenced on an illumina hiseq genome sequencer in a × 100 bp paired-end configuration. after sequencing, barcode adapters were removed and all reads were quality trimmed from both ends to a minimum phred score of all raw illumina sequences used in this study can be found on the ncbi sequence read archive. paired-end illumina genomic shotgun sequences from all populations were aligned to the combined transcriptome assembly described above using bwa, with default settings. aligning genomic sequences to a reference transcriptome creates mismatches at exon/intron boundaries because the introns are sampled by the genomic sequences but typically spliced out of the transcripts. when bwa encounters these mismatches it issues a “soft-clip” flag. to avoid snp calling errors that might be associated with these mismatches, all reads were first filtered to remove any reads that had been soft-clipped. we also removed all reads with a “mapq” alignment score <  finally, snps were extracted from these alignments using the mpileup program from the samtools package, at all sites with at least 5x coverage from nucleotides with a phred sequence quality ≥  all snp calls can be found in additional file to validate snps, we used the cleaved amplified polymorphic sequences technique. potential caps markers were identified in silico by screening alternate snp haplotypes against the commercially available restriction enzyme database found in biopython. we selected snps that created a polymorphism in a recognition site for the enzyme rsai. a region around each snp was amplified and digested in up to wcr individuals. primer sequences can be found in additional file pcr reactions were performed in a 20 μl mixture containing 1 μl dna template, 10 μl 2x pcr promega master mix and  mm each primer. amplified products for each sample were digested with 1 μl rsai at 37°c for 2 hrs. the resulting digest was visualized after electrophoresing in a 2% agarose gel stained with ethidium bromide. estimating population genetic parameters we estimated nucleotide diversity within populations for each unigene using a bias-corrected estimator derived for pooled populations of individuals. this estimator requires a minimum count for the inclusion of a particular snp variant in the calculation. here we used b = as our minimum. we also disregarded any sites with greater than > 40x coverage to minimize the inclusion of paralogous sequences or undetected copy number variants. within each population mean π was approximately, however in each population we observed to outlier unigenes. these outliers were removed prior to the calculation of summary statistics, as these loci are likely impacted by alignment issues given their extreme departure from the mean. these calculations were done using a custom python script. we computed pairwise fst between all populations for each unigene using a bias-corrected estimator for pooled population samples as implemented in popoolation, again using a minimum coverage of and a maximum coverage of for each snp locus. window and step sizes were both set to forcing the program to estimate fst once for each unigene. fst is a measure of genetic differentiation relative to the net genetic diversity of the populations sampled, which makes fst sensitive to fluctuations in net genetic diversity. to complement fst, several authors have concluded that absolute metrics of genetic differentiation should also be employed. for this study we created a simple statistic to estimate absolute genetic differentiation between populations using pooled sequence data. this statistic is an estimator of nei’s dxy statistic, though unlike dxy it uses estimated allele frequencies rather than allele counts, making it suitable for pooled sequence data. to estimate πxy we first recorded a four element sequence profile b b b b for every nucleotide site within each pooled population, that contains the frequency of a, c, g, and t bases in the sample. then for each unigene we calculated the average per-site expected heterozygosity among gametic unions between populations x and y, using the equation below. where n is the number of nucleotide sites available for a unigene, and are the frequencies of base bj at nucleotide site i in populations x and y, respectively. πxy requires only estimates of within population allele frequencies and not estimates of within population nucleotide diversity, unlike π and fst, and thus πxy does not require bias-correction for pooled populations. we applied the πxy metric to calculate genetic distances between populations for principal coordinates analysis. we also applied the πxy metric to characterize the genetic divergence between combined populations with low or high responsiveness to cry3bb or diapause and non-diapause populations. to estimate πxy from combined populations we used the average of the nucleotide frequencies within each population to avoid biases that could arise due to uneven sequence coverage between populations. as with π, we required a minimum coverage of and a maximum coverage of within each population for a site to be considered. we also required a minimum count of for the inclusion of a particular snp variant. for each unigene we identified πxy and fst outliers by permuting population membership and calculating both metrics for permuted pairs. in the case of the cry3bb response phenotype there are possible permutations between high and low response populations, and all permutations were used to create the null distribution. for diapause versus non-diapause populations there are potential permutations, among which we randomly selected to create the null distribution. in either case we selected the 99th percentile πxy and fst values from among these permuted populations as our threshold for significance, after removing all unigenes with fewer than available nucleotide sites. this threshold is arbitrary and lacks multiple test correction, though our goal is to produce a list of potential candidates that is inclusive, though it may include some false positives. also, given the permutation tests outlined above, the most confidence one can have regarding significance is 1/ or 1/, respectively. our sampling procedure involved estimating population parameters from five individuals per population across approx. loci. theory suggests that multi-locus estimates of population genetic parameters can be robust even when estimated from just five individuals. though this theory applies specifically to direct sequencing of individuals and not pooled population sequencing, we have conducted a power analysis of our pooled sampling design, and it confirms that adequate estimates can be made from just five individuals, especially for statistics utilizing the grand mean across all loci. in light of this, our grand mean genomic estimates should be quite accurate, despite the small within population sample size. due to the small within population sample size, we expect our single locus estimates of divergence to have relatively modest statistical power. these scans should be viewed as exploratory.  de novoassembly and annotation the de novo assembly of the wcr transcriptome yielded high quality unigenes totaling,715 bp, in addition to singleton reads. the length of unigenes varied from 99-  bp with an average of  bp. to determine the completeness of wcr transcriptome assembly, each unigene was compared to its putative ortholog from t. casteneum, d. melanogaster, and b. mori. overall, 44-46% of the unigenes had an ortholog hit ratio >  and 54-58% were > . on the assumption of conserved gene length between species, an ortholog hit ratio near zero indicates a poor assembly while values near one indicate a fully assembled transcriptome. following these criteria the current assembly for wcr appears to be moderately complete.figure transcriptome assembly of wcr. a) length distribution of all transcripts; each transcript is placed along x-axis in the ascending order based on its length, and b) ortholog hit ratio for assembled contigs calculated after blastx searches to genomes of t. casteneum, d. melanogaster, and b. mori. about 64% of the wcr unigenes had one or more hits to protein sequences in the non-redundant database at genbank while the remaining had no match. the majority of top blast hits for wcr unigenes were to insects, whereas a small proportion showed top hits to non-insect arthropods, other eukaryotes, and bacteria. as expected, more than two-thirds of the top hits for the wcr unigenes were to coleopteran sequences in the nr database such as those of t. casteneum, dendroctonus ponderosae, chrysomela tremula, leptinotarsa decemlineata, phaedon cochleariae, tenebrio molitor, and batocera horsfieldi. the wcr unigenes with no match to the nr database may represent novel genes with functions not yet been identified, or they could be a product of assembly errors. interestingly, the percentage of unmatched unigenes obtained in our study is considerably lower than those obtained for other insect transcriptomes. an interproscan search of these unmatched wcr unigenes against the nr database revealed hits to the protein signature databases for out of, suggesting that many have homologs in other species that were undetected using blastx alone. further, of the unmatched unigenes had an open reading frame of at least 50% of the contig length and had an orf of at least 90% of the contig length. a pairwise comparison to the recently sequenced genome of mountain pine beetle did not improve the annotation significantly as only out of unmatched wcr unigenes showed hits. comparative genomics using pairwise blastx searches to protein databases for four model insects, significant matches for a majority of wcr unigenes were obtained. the blastx search to the t. casteneum database showed the highest number of matches for wcr unigenes, followed by searches to n. vitripennis, b. mori, and d. melanogaster databases. several wcr unigenes showed matches to all four databases. however, there were a substantial number of wcr unigenes that matched uniquely to t. casteneum, n. vitripennis and b. mori. analyzing the pairwise comparisons, the highest level of agreement was observed between wcr and t. casteneum, likely reflecting the evolutionary relatedness between these two insect species and the completeness of the tribolium genome. functional annotation using blast2go, about 46% of wcr unigenes were annotated. the observed gene ontology terms for each domain were widely distributed into different categories. a comparison of percent mappings to each go category in wcr and t. casteneum revealed a highly similar distribution for both insect species. the majority of contigs assigned to the ‘biological process’ domain were involved in cellular, regulatory, and developmental activities, while most of the contigs under the ‘molecular function’ domain were predicted to have catalytic, binding and transporter functions. through the kegg-based pathway analysis using blast2go, a total of wcr unigenes were assigned to one or more of pathways. the majority of contig sequences were assigned to pathways for metabolism of organic compounds such as purine, pyrimidine, and glucose. the information on predicted pathways coupled with the functional annotation of the wcr transcriptome provides useful information for future genetic studies in this insect. we further characterized potential genes that may be involved in wcr bt resistance, specifically to cry3bb, focusing on detoxification and degradation, esterases, and cathepsins). we also identified several unigenes identified as putative receptors for bt toxins such as cadherins, abc transporters, aminopeptidases and glycosyltransferases. additional file lists the wcr unigenes that significantly matched these sequences in other insects, mainly from tribolium. the number of putative receptors totaled whereas the number of cathepsins, cyps, and esterases were much larger. these totals mark a considerable expansion beyond what was previously known for wcr candidate resistance gene targets. we further characterized two of the more important groups, the cyps and cathepsins, because of their known relevance to resistance. cytochrome p450s are a superfamily of enzymes which are ubiquitous among organisms and play an important role in metabolizing endogenous and exogenous substances. in insects, cyp genes are classified into four clades: cyp cyp cyp and mitochondrial p450s. despite their importance for insecticide resistance in wcr and other insects, little information is known regarding the diversity or genetic variation of wcr-cyps. in addition to blastx searches against the ncbi nr database, we identified the putative clades and families of wcr-cyps through a blastp comparison to a recently published tribolium cyp database. the majority of wcr-cyps had highest similarity to clade, with a moderate number showing similarity to clade. the expanded sequencing combined with previous resources show a considerable diversity in the number of cathepsins in wcr. we found a total of unigenes that matched known cathepsins, including some with high similarity to those already described in wcr. a previous est-based sequencing project, revealed cathepsin-like fragments, corresponding to two different classes, l and b. in our transcriptome data, unigenes belonged to class b and matched class l; we also identified unigenes in the lesser known d and f classes. cathepsins play a pivotal role in wcr digestion, and have been implicated in the evolution of rotation resistance. normally, wcr feed and oviposit in corn, but rotation resistant wcr feed on soybean before oviposition. when these soybean fields are rotated to corn the following year, the rotation resistant wcr eggs hatch and larvae feed on emerging corn plants. curzi et al. found that wcr feeding on soybean leaves had higher levels of cathepsin-l like activity compared to wcr fed on corn. this expanded set of cathepsin-l like sequences and targeted resequencing of wcr populations could help enable the development of a molecular diagnostic for rotation resistant wcr. for all unigenes, including those that may be involved in wcr resistance to cry3bb we pooled snp data from different populations with high or low response to cry3bb and screened for genes that were πxy and fst outliers. we identified πxy outliers and fst outliers. only one unigene was found on both lists. this unigene shares sequence similarity with a gene of unknown function from t. castaneum. among all πxy and fst outlier unigenes we found two on our lists of candidates for genes potentially involved in reduced response to cry3bb these include one esterase and one cyp, both of which were fst outliers. all outliers should be considered preliminary and exploratory as our tests are low-powered given the relatively small within population sample sizes, and because our strategy assumes that populations with a low response to cry3bb share resistance alleles, though there is no evidence supporting this assumption. nonetheless, these loci may represent a promising pool of potential candidates for cry3bb resistance.figure calculated π xy and f st values between populations with high and low response to cry3bb versus permuted 99 th percentile values. each point represents a unigene’s calculated π xy or f st value between pooled high and low response to cry3bb populations and the 99th percentile value from all possible permutations of these populations. red dots indicate calculated π xy or f st values > 99th percentile permutation value, while black dots represent calculated π xy or f st values < 99th percentile permutation value. the blue line represents the boundary between outliers and non-outliers. wcr snp development and validation the field collected populations included a wide geography across the us corn belt, while the laboratory populations included several common research strains. each population consisted of individuals that were pooled and sequenced to approximately 20x-30x coverage using whole genome shotgun sequencing. these sequences were aligned to the reference transcriptome for snp discovery and estimation of population genetic parameters. in total we found snps segregating among the haploid chromosome sets sampled. snp discovery saturates quickly among the wcr populations we sampled. after sampling populations, few new snps are detected among the remaining populations. those new snps that are found by adding more than populations tend to be rare, often occurring in just one population. thus, within our unigene collection, it is possible that our data has revealed most of the common snps in corn belt wcr.figure rarefaction curve of snp discovery among the wcr populations sampled. starting with the first population sampled and continuing until the 26th we record the number of new snps each population contributes. the slope of the rarefaction curve is dependent on population input order, which in this case is arbitrary. to account for this, we selected random input orders and computed the mean number of new snps discovered and the upper and lower 95% confidence intervals. to address the validity of our snp predictions and their utility for genotyping we identified that were predicted to create a polymorphic restriction enzyme site and used this to genotype individuals from the inbred and randomly mated brookings population using caps. for of the snps the caps genotypes showed evidence of both allelic versions of the predicted snp. the remaining snps were homozygous for one allele among all individuals tested, which could reflect inaccurate snp predictions or simply the absence of one allele among the populations tested. in any case, these results demonstrate that the predicted snps are reliable and can be readily utilized in genotyping assays. population genomics of western corn rootworm using the snps extracted from pooled population sequencing we estimated genetic diversity within and between wcr populations. to estimate within population nucleotide diversity we used a bias-corrected estimator, suitable for pooled population sequences. for between population comparisons we created a simple statistic, which is an estimator of the dxy statistic. we focused our analyses on the wcr unigenes with a hit to a protein in the ncbi nr database. for these unigenes, the mean π within each population can be found in table  nucleotide diversity within populations ranged from -, with a mean of. as expected, the inbred brookings colony showed the lowest mean π, and the highest mean π values were found among field collected populations. our extensive snp dataset reveals several interesting aspects of genetic variation in wcr. first, the inbred population had experienced generations of full-sib mating, which is theoretically expected to fix approximately 41% of the sites segregating among the gametes of the founding sib-mated pair, yet we observed only a slight reduction in π within the inbred brookings population compared to the randomly mated brookings population from which it was derived. this suggests that full-sib mating has been fairly ineffective at removing variation. second, there appears to be slightly less nucleotide diversity within the collections, when compared to the collections. this could reflect annual variation in population sizes across the corn belt, or perhaps some unnoticed deviation in our sampling procedures between these two years. finally, there is little evidence for a significant reduction of genetic diversity among lab populations, suggesting that maintenance in artificial environments has not caused an erosion of genetic diversity. this finding mirrors observations by kim et al., who compared a laboratory population derived from the non-diapause brookings population to field populations using microsatellite markers. among the pairwise comparisons between populations, πxy ranged between and with a mean of. from this we observed that the levels of within population nucleotide diversity are of a similar magnitude as between population divergence. this indicates that each population maintains considerable genetic variation, and that there is little evidence for substantial genetic differentiation between populations. pairwise πxy values were also used to generate a principal coordinates analysis among all populations. diapause status appears to be the most important factor in population differentiation. the natural state for wcr is diapause, and the trait shared by all non-diapause populations was inherited from a single source. this non-diapause trait was developed after several generations of mass selection, and because of its complex inheritance, it has been suggested that it is likely a complex polygenic trait. consistent with this hypothesized mode of inheritance and strong artificial selection for its maintenance, we find unigenes that are πxy outliers when comparing all diapause versus all non-diapause populations. these outlier unigenes may include causal alleles for the non-diapause phenotype, but likely also include linked loci that have been differentiated between diapause and non-diapause populations by genetic hitchhiking.figure principal coordinate analysis using genetic distances between wcr populations. diapause populations are colored in black and non-diapause population are red. populations with a reduced response to cry3bb are denoted with open squares, while populations with unknown or low to intermediate resistance are denoted with filled circles. to more formally address population differentiation we estimated fst pairwise between all populations for all unigenes with a hit to the ncbi nr database. mean fst between all populations across all unigenes was, and only % of all unigenes had fst values > . the overall pattern of low fst indicates little genetic differentiation among wcr populations, even across distances of approximately  km. these findings are consistent with earlier work in wcr, which also found low levels of genetic differentiation between field collected populations. they are also consistent with our observation that new snp discovery saturates fairly quickly as new corn belt populations are sampled. among all pairwise population contrasts, the highest fst values were found between field and lab populations, with nearly all of the most extreme values between a field population and the inbred brookings population. this observation is expected, as inbreeding impacts allele frequencies and will tend to increase fst. interestingly, the smallest mean pairwise fst estimate was between the onslow, ia and trent, sd populations. despite being 508 km apart, these two field collected populations have both demonstrated reduced response to cry3b in a lab-based single plant bioassays. it is premature to conclude that these populations derive from a single cry3bb resistant source, though this intriguing possibility warrants further research. in the last 75 years wcr have expanded from their native range in the us southwest into the corn belt and further eastward, recently arriving at the eastern seaboard. given this recent and rapid range expansion, it is not surprising that our results and previous efforts reveal limited genetic differentiation among corn belt wcr populations, even across distances of approximately  km. theoretical treatments of the genetic properties of range expansion typically find that populations near the edges of the range have reduced genetic diversity and increased fst when compared to populations near the center of the range, and that both of these features reset after colonization at a rate that is proportional to the migration rate between populations. among the corn belt wcr populations, we find high genetic variation and low fst values among the most recently colonized easterly populations sampled. these populations are now over  km west of the eastern edge of the species range, and our genetic observations are consistent with theoretical predictions, in that these populations are no longer in their colonizing phase and are likely gaining new genetic diversity through gene flow with more westerly populations. this pattern may also suggest that migration within the corn belt is fairly high; enough to spread considerable genetic diversity into populations founded in the last 40 years. another interpretation could be that these two eastern populations were formed by a large and diverse set of founders, and that this diversity has been simply maintained since founding. without historical samples it may be difficult to differentiate between these two hypotheses, though future studies could sample recently colonized populations in the eastern us to estimate typical founding population sizes. we do observe a modest – albeit significant – trend of increasing genetic isolation as a function of geographic distance among field collected populations. thus, there may be a subtle signature of the range expansion in the form of genetic isolation by distance. moreover, the significance of the isolation by distance is strongly influenced by the inclusion of the most easterly population. if we remove this population, the correlation between genetic and physical distance falls to and is no longer significant. in contrast, if we remove all populations from eastern colorado and western kansas the correlation between genetic and physical distance is still significant. the importance of this eastern population in detecting isolation by distance is consistent with incomplete homogenization following range expansion, however, because this result is driven by a single population we cannot exclude the hypothesis this population is in some way anomalous. moreover, because recently formed corn belt wcr populations may not be at migration-drift equilibrium we cannot determine if the weak pattern of genetic isolation by distance is a product of genetic drift overcoming gene flow or if it is a hold-over from population founder effects during the species eastward invasion. based on theory and these preliminary results we predict that sampling populations further east will tend to increase the extent of detectable of isolation by distance in wcr, and may reveal important estimates of genetic diversity in recently colonized areas.figure genetic isolation among wcr populations as a function of geographic distance. following rousset, pairwise estimates of f st were linearized; y-axis) and regressed onto the natural log of the physical distance between populations. the linear regression line is shown in red along with the pearson’s correlation coefficient. in a review article, mcdonald and linde posit that knowing information about the population genetics of a pest organism can lead to a better understanding of its evolutionary potential to evade control strategies. they proposed that species with high levels of genetic diversity and high gene flow between populations have the greatest evolutionary potential to overcome control options through darwinian evolution. in this study, and in earlier studies, we find evidence to support the hypothesis that wcr has both high levels of genetic diversity and high gene flow between populations. for example, assuming wcr has a mutation rate between 10- or 10-9 bp- generation-, the observed values of π imply that wcr field populations maintains approximate effective populations sizes on the order of to million individuals. this level of ne gives little indication of recurrent population bottlenecks in wcr, which gives wcr populations the ability to generate and maintain many individuals with novel mutations. this creates a large pool of standing genetic variation that may in turn increase the probability of evolving resistance to control strategies. moreover, we show very little genetic differentiation between field populations, and that genetic variation has spread quickly eastward following range expansion, both of which are consistent with high rates of gene flow between populations. thus, if resistance does evolve, there are few barriers preventing it from spreading between populations. together these two features of wcr population genetics are hypothesized to lead to a high evolutionary potential to evade controls, and may help to explain why wcr has a history of evolving resistance to insecticides and agricultural practices.  wcr is one of the most significant agricultural pests in north american, yet at the time of publication wcr has a modest number of nucleotide sequences are available in genbank. improved genetic and genomic resources for wcr are crucial tools in the ongoing effort to control this pest. here we produce a large transcriptome assembly for wcr and complement these unigenes with diversity data obtained by whole genome population resequencing. our wcr transcriptome assembly contains approximately unigenes, including several gene families that have been implicated in insecticide resistance in other species. we identify a large pool of snps among wcr populations, and show that these snps can be readily adapted for use as genotyping markers. we also use these snps to scan for outliers among candidate bt resistance genes and to understand how population processes have shaped genetic variation in this species. we find two potential candidate bt resistance genes that show strong patterns of molecular differentiation between high and low response to cry3bb populations. we also expand on past population genetic studies in wcr and show that this species has significant genetic diversity, little population structure, and a high evolutionary potential to resist control strategies. this study will provide a solid foundation for future research on the molecular genetics of wcr. availability of supporting data ests and illumina pooled population sequences used in this study can be found on the ncbi sequence read archive under the accessions srs and srp respectively. electronic supplementary material additional file 1: genbank ids for public wcr est sequences. this spreadsheet contains the genbank ids for the public wcr sequences incorporated our transcriptome assembly. additional file 2: population level wcr snp genotypes. a comma separated values file containing the genotypic status of all populations sampled at snp sites. missing genotypes are denoted by a dash. the sequence of each unigenes is listed on first usage. this file has been compressed with the bzip compression software. additional file 3: snp genotyping validation study. this spreadsheet contains the caps snp validation information, including primer sequences, and observed restriction patterns and inferred genotypes. additional file 4: power analysis of population genetic sampling procedure. additional file 5: blast2go annotations. this spreadsheet contains annotations for wcr unigenes as identified by blast2go. additional file 6: comparative genomics of wcr. a venn diagram showing the number of wcr contigs with significant matches to genomes of t. casteneum, d. melanogaster, n. vitripennis, and b. mori. the significant matches were calculated after pairwise comparisons to each individual genome. additional file 7: distribution and comparison of go categories. vertical bars indicate the distribution of wcr and red flower beetle go term mappings that belong to each of the three top-level go categories. additional file 8: kegg pathway annotations. this file contains two spreadsheets. the “detailed” spreadsheet contains putative kegg pathways assignments to unigenes in the wcr transcriptome, while the “summary” spreadsheet lists a summary of the counts of all pathways identified. additional file 9: identification of putative candidate genes for cry3bb resistance. this file contains a spreadsheet for each of several putative gene families that may be candidates for resistance to bt toxins, including the unigenes name and descriptions, and results from π xyand f st outlier scans between cry3bb resistant and susceptible populations. the last two spreadsheets list the results of π xyand f st outlier scans between populations with high or low response to cry3bb for all available unigenes. additional file 10: characterization of wcr-cyp genes. blastp result comparing wcr-cyp unigenes to a tribolium cyp data set. additional file 11: identification of candidate genes for non-diapause trait. this spreadsheet contains the results of π xyoutlier scans between diapause and non-diapause populations for all available unigenes. along with the scan results, unigenes name and descriptions are also given. additional file 12: plot of calculated π xy values between diapause and non-diapause populations and permuted 99 th percentile values. each point represents a unigenes calculated π xy value and the 99th percentile value from permutations of the diapause and non-diapause population labels. red dots indicate calculated π xy > 99th percentile permutation value, while black dots represent calculated π xy < 99th percentile permutation value. the blue line represents the boundary between outliers and non-outliers. additional file 13: a histogram of pairwise f st between wcr populations. this histogram represents f st estimates for all genes among all pairwise comparisons of the wcr populations. competing interests lef, rak, m chen, m carroll, bsg, rf and tc are employed by monsanto company and some of the authors are also shareholders in monsanto company, though individually or collectively in no way represent any controlling interest in the affairs of monsanto company. results described in the manuscript are the product of work financed by monsanto company, and the article processing charge will be paid by monsanto company. the authors have an obligation under an employment contract with monsanto company to assign the rights of any invention made in the course of their employment with monsanto company to the company. monsanto company has filed one or more patent applications disclosing and claiming all or some portion of the results of the work described in the manuscript. the authors are not known to have any financial competing interests in addition to those listed above. authors’ contributions lef, rb, rak, m chen, rf, tc, bsg, and apm conceived the study and contributed to experimental design. lef, rb and apm analyzed the data and drafted the manuscript. m chen and m carroll contributed materials. all authors read and approved the final manuscript. results described in the manuscript are the product of work financed by monsanto company. rb and apm received additional support from the ohio agricultural research and development center, the ohio state university. we thank aaron gassmann at iowa state university, bruce hibbard at usda ars, columbia, mo, and b. wade french at the usda north central agricultural research laboratory for sharing wcr colonies. we are grateful to john cantwell, douglas jones, and leanna trail for their assistance in wcr field collections, oliver ilagan, kaylee kuban, cara vazquez for colony maintenance, and nay thane for producing the transcriptome assembly. finally, we thank john ledeaux and peter morrell for helpful comments and discussions.
among these unigenes we found several gene families that have been implicated in insecticide resistance in other species. using these tools we identified candidate gene families for insecticide resistance and reveal aspects of wcr population history in light of the species’ recent range expansion. here we produce and analyze a transcriptome assembly for wcr. the wcr transcriptome assembly presented here contains approximately unigenes, many of which have high similarity to other insect species. improved genetic and genomic resources are crucial tools for understanding population history and the genetic basis of trait evolution. we also perform whole genome population resequencing, and combine these resources to better understand the evolutionary history of wcr. we used these snps to scan for outliers among the candidate genes, and to understand how population processes have shaped genetic variation in this species. this study highlights the utility of transcriptomic and genomic resources as foundational tools for dealing with highly adaptive pest species. in addition, this species has a history of evolving traits that allow it to escape various control options. western corn rootworm is one of the most significant insect pests of maize in north america. diabrotica virgifera virgiferainsect resistancetransgenic croptranscriptomewestern corn rootwormissue-copyright-statement© the author 2014
 western corn rootworm is one of the most significant insect pests of maize in north america. wcr has dramatically increased its range in the last century, invading key maize production areas in the us and abroad. in addition, this species has a history of evolving traits that allow it to escape various control options. improved genetic and genomic resources are crucial tools for understanding population history and the genetic basis of trait evolution. here we produce and analyze a transcriptome assembly for wcr. we also perform whole genome population resequencing, and combine these resources to better understand the evolutionary history of wcr.  the wcr transcriptome assembly presented here contains approximately unigenes, many of which have high similarity to other insect species. among these unigenes we found several gene families that have been implicated in insecticide resistance in other species. we also identified over unigene based snps among wcr populations. we used these snps to scan for outliers among the candidate genes, and to understand how population processes have shaped genetic variation in this species.  this study highlights the utility of transcriptomic and genomic resources as foundational tools for dealing with highly adaptive pest species. using these tools we identified candidate gene families for insecticide resistance and reveal aspects of wcr population history in light of the species’ recent range expansion. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords diabrotica virgifera virgiferainsect resistancetransgenic croptranscriptomewestern corn rootwormissue-copyright-statement© the author 2014 