A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species — TEI-tagged view

Robert J. ElshirepersName, Jeffrey C. GlaubitzpersName, Qi SunpersName, Jesse A. PolandpersName, Ken KawamotopersName, Edward S. BucklerpersName, Sharon E. MitchellpersName (2011date). “A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species.” PLoS ONEtitle·j.

Object of study — a restriction-enzyme, multiplexed genotyping-by-sequencing protocol for genetically diverse, large-genome species

Research question — Can genome complexityterm·dom reduction and optimized barcoding make high-density genotyping simple, reproducible, inexpensive, and usable without a complete reference genome?

Contribution — a simplified GBSterm·mth library protocol, barcodeterm·mth design, filtering/mapping pipeline, and maizeterm·dom/barleyterm·dom validation for high-diversity species

Full text — tagged

A robust, simple genotyping-by-sequencing (GBS)term·mth approach for high diversity species Abstract Advances in next generation technologies have driven the costs of DNA sequencing down to the point that genotyping-by-sequencing (GBS)term·mth is now feasible for high diversity, large genome speciesterm·dom. Here, we report a procedure for constructing GBS librariesterm·mth based on reducing genome complexityterm·dom with restriction enzymes (REs)term·mth. This approach is simple, quick, extremely specific, highly reproducible, and may reach important regions of the genome that are inaccessible to sequence capture approachesterm·mth. By using methylation-sensitive REsterm·mth, repetitive regions of genomesterm·dom can be avoided and lower copy regionsterm·dom targeted with two to three fold higher efficiency. This tremendously simplifies computationally challenging alignment problemsterm·dom in species with high levels of genetic diversity. The GBSterm·mth procedure is demonstrated with maizeterm·dom (IBM) and barleyterm·dom (Oregon Wolfe Barleyterm·dom) recombinant inbred populations where roughly 200,000 and 25,000 sequence tagsterm·mth were mapped, respectively. An advantage in species like barleyterm·dom that lack a complete genome sequence is that a reference mapterm·mth need only be developed around the restriction sites, and this can be done in the process of sample genotyping. In such cases, the consensus of the read clustersterm·mth across the sequence tagged sites becomes the reference. Alternatively, for kinship analyses in the absence of a reference genome, the sequence tagsterm·mth can simply be treated as dominant markersterm·mth. Future application of GBSterm·mth to breeding, conservation, and global species and population surveys may allow plant breeders to conduct genomic selectionterm·dom on a novel germplasm or species without first having to develop any prior molecular tools, or conservation biologists to determine population structureterm·dom without prior knowledge of the genome or diversity in the species.

Introduction During the last decade, extensive public resources were dedicated to genotyping humans, a species with relatively low genetic diversity (about one substitution per thousand nucleotides) [1] – [3] . Many species including maizeterm·dom [4] , [5] , Drosophila [6] , and some bacteria [7] , however, are at least 10 times more diverse than humans (more than one substitution per hundred nucleotides). Besides containing high levels of nucleotide diversity, the maizeterm·dom genome also exhibits frequent transposon-mediated rearrangements that produce extensive presence/absence variation that often encompasses genic regions [8] – [10] . Standard, fixed-sequence approachesterm·mth like single base extension assays or microarrays require invariant primer binding sites in order to obtain consistent results. Such invariant regions are often difficult to find in maizeterm·dom [11] . Furthermore, the large-scale structural variation also complicates DNA sequence alignment, resulting in a maizeterm·dom “reference” genome that contains only 70% or less of the species-wide genome space [12] .

Although abundant diversity is a challenge to assays that rely on scoring fixed positions, it is advantageous to direct sequencing approachesterm·mth because sequencing efficiency for genotyping scales directly with genetic diversity. We have developed a technically simple, highly multiplexed, genotyping-by-sequencing (GBS)term·mth approach that is suitable for population studies, germplasm characterization, breeding, and trait mapping in diverse organisms. This procedure, which can be generalized to any species at a low per-sample cost, is based on high-throughput, next-generation sequencingterm·mth of genomic subsets targeted by restriction enzymes (REs)term·mth.

Next-generation sequencingterm·mth (NGS) technologies have been recently used for whole genome sequencing and for re-sequencing projects where the genomes of several specimens are sequenced to discover large numbers of single nucleotide polymorphisms (SNPs)term·dom for exploring within-species diversity, constructing haplotype maps and performing genome-wide association studies (GWAS)term·dom [13] . Multiplex sequencingterm·mth has also been accomplished by tagging randomly sheared DNA fragments from different samples with unique, short DNA sequencesterm·mth (barcodesterm·mth) and pooling samples into a single sequencing channel [14] . This approach (random DNA shearing followed by barcodeterm·mth tagging) works very well for species with small genomes, including organellar and microbial DNAs, and has been used to rapidly determine the complete chloroplast genome sequences of spruce and several pine species [15] and for discovery and mapping of genomic SNPs in rice [16] , [17] .

Although GBSterm·mth is fairly straightforward for small genomes, target enrichment or reduction of genome complexityterm·dom must be employed to ensure sufficient overlap in sequence coverage for species with large genomes. Enrichment strategies including long range PCR-amplification of specific genomic regions, use of molecular inversion probes, and various DNA hybridization/sequence capture methods [18] are time-consuming, technologically challenging, and can be cost-prohibitive for assaying large numbers of samples. Reducing genome complexityterm·dom with restriction enzymes (REs)term·mth, however, is easy, quick, extremely specific, highly reproducible, and may reach important regions of the genome that are inaccessible to sequence capture approachesterm·mth. By choosing appropriate REs, repetitive regions of genomesterm·dom can be avoided and lower copy regionsterm·dom can be targeted with two to three fold higher efficiency [12] , [19] , which tremendously simplifies computationally challenging alignment problemsterm·dom in species with high levels of genetic diversity.

The value of sequencing restriction site-associated genomic DNA (i.e., RAD tagsterm·mth) for high-density SNP discovery and genotyping was first demonstrated by Baird et al. [20] . Increased efficiency and cost benefits were realized by incorporating a multiplex sequencingterm·mth strategy that uses an inexpensive barcoding system. Because barcodesterm·mth are included in one of the adapter sequences (i.e., they are not added to individual DNA samples by PCR), reagent costs for constructing sequencing libraries are minimized. The location of the barcodeterm·mth, just upstream of the RE cut-site in genomic DNA, also eliminates the need for a second Illumina sequencing (“indexing”) read. The present work describes an even more cost-effective genotyping procedure based on NGS technology (Illumina, Inc.). The barcoding strategy is similar to RAD but modulation of barcodeterm·mth nucleotide composition and length results in fewer sequence phasing errorsterm·dom. Compared to the RAD method, the procedure described here is substantially less complicated; generation of restriction fragments with appropriate adapters is more straightforward, single-well digestion of genomic DNA and adapter ligation results in reduced sample handling, there are fewer DNA purification steps and fragments are not size selected. Costs can be further reduced via shallow genome sampling coupled with imputationterm·mth of missing internal SNPs in haplotype blocks. The following protocol was initially developed for maizeterm·dom, a genetically diverse (see above), large genome species (2.3 Gbp) [21] . We have since used this procedure for genotyping and mapping in several other species. Results for both maizeterm·dom and barleyterm·dom are reported herein.

Methods DNA Samples Samples comprised the parents and 276 recombinant inbred lines (RILs) from a high resolution maizeterm·dom mapping population (IBM [22] ), and the parents and 43 doubled haploid (DH) barleyterm·dom lines from the Oregon Wolfe Barleyterm·dom (OWB) mapping population [23] . The 43 barleyterm·dom lines were selected from the larger set of 83 OWB lines to maximize recombination. High molecular weight DNAs were extracted from leaves of single plants using a standard CTAB protocol [24] .

Choosing REs and Adapter Design Selection of REs that leave 2 to 3 bp overhangs and do not cut frequently in the major repetitive fraction of the genome is of critical importance. A suitable RE for maizeterm·dom and close relatives (teosintes) is Ape KIterm·mth, a type II restriction endonucleaseterm·mth that recognizes a degenerate 5 bp sequence (GCWGC, where W is A or T), creates a 5′ overhang (3 bp), has relatively few recognition sites in the major classes of maizeterm·dom retrotransposons, and is partially methylation sensitive (will not cut if the 3′ base of the recognition sequence on both strands is 5-methylcytosine). Using an RE that leaves an overhang comprising more than one nucleotide is extremely useful in promoting efficient adapter ligation to insert DNA.

Two different types of adapters were used in this protocol. The “barcodeterm·mth” adapter terminates with a 4 to 8 bp barcodeterm·mth on the 3′ end of its top stand and a 3 bp overhang on the 5′ end of its bottom strand that is complementary to the “sticky” end generated by Ape KIterm·mth (CWG). The sequences of the two oligonucleotides comprising the barcode adapterterm·mth are: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCTxxxx and 5′-CWGyyyyAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT and, where “xxxx” and “yyyy” denote the barcodeterm·mth and barcodeterm·mth complement and sequences, respectively ( Figure 1 ). The second, or “common”, adapter has only an Ape KI-compatible sticky end: 5′-CWGAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAG and 5′-CTCGGCATTCCTGCTGAACCGCTCTTCCGATCT ( Figure 1 ). Adapters were designed so that the Ape KIterm·mth recognition site did not occur in any adapter sequence and was not regenerated after ligation to genomic DNA. Adapter design also allows for either single-end or paired-end sequencing on the Illumina, Inc. (San Diego, CA) NGS platforms.

GBSterm·mth adapters, PCR and sequencing primers. (a) Sequences of double-stranded barcodeterm·mth and common adaptersterm·mth. Adapters are shown ligated to Ape KI-cut genomic DNA. Positions of the barcodeterm·mth sequence and Ape KIterm·mth overhangs are shown relative to the insert DNA; (b) Sequences of PCR primer 1 and paired end sequencing primer 1 (PE-1). Binding sites for flowcell oligonucleotide 1 and barcode adapterterm·mth are indicated; (c) Sequences of PCR primer 2 and paired end sequencing primer 2 (PE-2). Binding sites for flowcell oligonucleotide 2 and common adapterterm·mth are indicated.

A compatible set of 96 barcodeterm·mth sequences that have been used for multiplex sequencingterm·mth is provided as supporting information ( Table S1 ). To minimize the possibility of misidentifying samples as a result of sequencing or adapter synthesis error, all pair-wise combinations of barcodesterm·mth differed by a minimum of three mutational steps. Hence, it should be possible to correctly assign samples with single base barcodeterm·mth sequencing errors, or to identify particular adapters with high rates of synthesis error [25] . To avoid the potential loss of sequence quality due to phasing errorsterm·dom caused by reading through a non-variable restriction site prior to the twelfth base, or through an adapter position with a highly skewed base ratio [( http://www.illumina.com/Documents/products/technotes/technote_rta_theory_operations.pdf )], barcodeterm·mth lengths were modulated from 4 to 8 bp and care was taken to maximize the balance of the bases at each position in the overall set. For barcodesterm·mth larger than 5 bases, mononucleotide runs of 3 or more, and barcodesterm·mth that contained sequences of smaller barcodesterm·mth were disallowed.

Preparation of Libraries for Next-Generation Sequencingterm·mth A basic schematic of the protocol used for performing GBSterm·mth is shown in Figure 2 . Oligonucleotides comprising the top and bottom strands of each barcode adapterterm·mth and a common adapterterm·mth were diluted (separately) in TE (50 µM each) and annealed in a thermocycler (95°C, 2 min; ramp down to 25°C by 0.1°C/s; 25°C, 30 min; 4°C hold). Barcodeterm·mth and common adaptersterm·mth were then quantified using an intercalating dye (PicoGreen®; Invitrogen, Carlsbad, CA), diluted in water to 0.6 ng/µL (∼02 pmol/µL), mixed together in a 1∶1 ratio, and 6 µL (∼0.06 pmol each adapter) of the mix was aliquoted into a 96-well PCR plate and dried down. DNA samples (100 ng in a volume of 10 µL) were added to individual adapter-containing wells and plates were, again, dried.

Steps in GBSterm·mth library construction. Note: Up to 96 DNA samples can be processed simultaneously. (1) DNA samples, barcodeterm·mth, and common adapterterm·mth pairs are plated and dried; (2–3) samples are then digested with Ape KIterm·mth and adapters are ligated to the ends of genomic DNA fragments; (4) T4 ligase is inactivated by heating and an aliquot of each sample is pooled and applied to a size exclusion column to remove unreacted adapters; (5) appropriate primers with binding sites on the ligated adapters are added and PCR is performed to increase the fragment pool; (6–7) PCR products are cleaned up and fragment sizes of the resulting library are checked on a DNA analyzer(BioRad Experion® or similar instrument). Libraries without adapter dimersterm·dom are retained for DNA sequencing.

Samples (DNA plus adapters) were digested for 2 h at 75°C with Ape KIterm·mth (New England Biolabs, Ipswitch, MA) in 20 µL volumes containing 1× NEB Buffer 3 and 3.6 U Ape KIterm·mth. Adapters were then ligated to sticky ends by adding 30 µL of a solution containing 1.66× ligase buffer with ATP and T4 ligase (640 cohesive end units) (New England Biolabs) to each well. Samples were incubated at 22°C for 1 h and heated to 65°C for 30 min to inactivate the T4 ligase. Sets of 48 or 96 digested DNA samples, each with a different barcode adapterterm·mth, were combined (5 µL each) and purified using a commercial kit (QIAquick PCR Purification Kit; Qiagen, Valencia, CA) according to the manufacturer's instructions. DNA samples were eluted in a final volume of 50 µL. Restriction fragments from each library were then amplified in 50 µL volumes containing 2 µL pooled DNA fragments, 1× Taq Master Mix (New England Biolabs), and 25 pmol, each, of the following primers: (A) 5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT and (B) 5′-CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT . These primers contained complementary sequences for amplifying restriction fragments with ligated adapters, binding PCR products to oligonucleotides that coat the Illumina sequencing flow cell and priming subsequent DNA sequencing reactions [26] ( Figure 1 ).

Temperature cycling consisted of 72°C for 5 min, 98°C for 30 s followed by 18 cycles of 98°C for 30 s, 65°C for 30 s, 72°C for 30 s with a final Taq extension step at 72°C for 5 min. These amplified sample pools constitute a sequencing “library.” Libraries were purified as above (except that the final elution volume is 30 µL) and 1 µL was loaded onto an Experion® automated electrophoresis station (BioRad, Hercules, CA) for evaluation of fragment sizes. Libraries were considered suitable for sequencing if adapter dimersterm·dom (∼128 bp in length) were minimal or absent and the majority of other DNA fragments were between 170–350 bp. If adapter dimersterm·dom were present in excess of 0.5% (based on the Experion® output), libraries were constructed again using a few DNA samples and decreasing adapter amounts. Guidelines for adapting the protocol to different species including details for performing adapter titrations and are provided in Supporting Information ( Text S1 , Figure S1 and Figure S2 ).

Once the appropriate quantity of adapters was empirically determined for a particular enzyme/species combination, no further adapter titration was necessary. Single-end sequencing (86 bp reads) of one 48- or 96-plex library per flowcell channel, was performed on a Genome Analyzer II (Illumina, Inc., San Diego, CA). See Bentley et al. [26] for details of the sequencing process and chemistry.

Filtering Raw Sequence Data Analyses of the 86 bp sequencing reads were based upon the unfiltered qseq filesterm·mth, since the filtering process that produces fastq filesterm·mth sometimes discarded good reads that aligned perfectly to the reference genome for at least 64 bases. Starting with the qseq filesterm·mth from a flow cell, we first filtered for reads that (1) perfectly matched one of the barcodesterm·mth and the expected four-base remnant of the Ape KIterm·mth cut site (CWGC), (2) were not adapter/adapter dimersterm·dom, and (3) contained no “Ns” in their first 72 bases. These reads were sorted into separate files according to their barcodeterm·mth, with the barcodeterm·mth removed and the remainder of the sequence trimmed to 64 bases (including the initial CWGC). If either the full Ape KIterm·mth site (from partial digestion or chimera formation) or the first 8 bases of common adapterterm·mth (from Ape KIterm·mth fragments less than 64 bases) were detected within 64 bases, the read was truncated appropriately and then filled to 64 bases with polyA.

For maizeterm·dom, subsequent filtering of the reads was then done in two different ways, depending on our purpose. To generate a reference set of 64 base sequence tagsterm·mth to be included in a presence/absence genotype tableterm·mth, only reads with a minimum Q-score of 10 across the first 72 bases) and that occurred at least twice were kept. We opted to use this somewhat low-stringency minimum Q-score cutoff to maximize the number of useful sequence tagsterm·mth. Sequence tagsterm·mth containing random sequencing errors should not occur multiple times in multiple samples and should not map genetically, so they should be filtered out in subsequent steps. To this set of reference tagsterm·mth, the expected 64 base tags from an in silico Ape KIterm·mth digest of the maizeterm·dom reference genome, B73 RefGen v1 [21] , were added (with fragments shorter than 64 bases filled with polyA, as above). To fill in the observed counts in the genotype table, a second pass across the reads for each DNA sample was performed. In this second pass, 64 base reads were counted for each sample (and the count added to the genotype table) if they perfectly matched one of the reference tagsterm·mth, regardless of their minimum Q score. The resulting genotype table was then filtered to remove tags that occurred in 10 or fewer DNA samples; this should remove most of the sequencing errors. For barleyterm·dom, the absence of a reference genome prevented anchoring reads to a physical map. Sequence reads were simply filtered for unique 64 base sequence reads that were present in five or more lines and these were mapped genetically as described below.

All maizeterm·dom and barleyterm·dom sequences were submitted to the National Center for Biotechnology Information (NCBI) Short Read Archive (study SRP004282.1 ).

DNA sequence alignments The filtered sequence reads were first aligned to the maizeterm·dom reference genome (B73 RefGen v1) using the Burrows-Wheeler alignment tool (BWA)term·mth [27] , allowing a maximum of four mismatches and one gap of up to 3 bp. The Basic Local Alignment Search Tool (BLAST)term·mth [28] was used to query reads that were not aligned by BWA, first against the maizeterm·dom reference genome with an e-value cutoff of 1e −2 and then against the National Center for Biotechnology Information (NCBI) nt database using default settings.

Mapping Presence/absence scores for each tag were used in a binomial test of segregationterm·mth versus an independent framework mapterm·mth. For maizeterm·dom, this framework mapterm·mth consisted of 644 SNPs genetically mapped in the maizeterm·dom nested association mapping (NAM) population [29] and then genotyped in the IBM population. The binomial segregation test filtered for sequence tagsterm·mth that co-segregated with only one of the two parental alleles at a given SNP. For each SNP markerterm·dom, the two possible parental sources of a tag were each tested in turn. A “success” was recorded when a tag co-occurred in a RIL with the SNP allele from its presumed parental source, otherwise a “failure” was recorded. The binomial sample size was the number of RILs in which the tag was present and the SNP was not missing or heterozygous. For maizeterm·dom, tests were only performed if the sample size was at least 10. The probability of success was defined as the proportion of the RILs that contained the SNP allele being tested. For maizeterm·dom, a threshold p -value of 0.001 was considered significant for directed tests versus the physically closest SNP, or 0.0001 for elsewhere in the genome.

For barleyterm·dom, mapping was conducted using flanking SNPs and a threshold of p <0.0001 for the binomial test. In practice, a sequence tagterm·mth was mapped in barleyterm·dom only if it always co-occurred with one SNP allele and never the other.

In maizeterm·dom only, biallelic GBS markersterm·mth were identified as follows. Pairs of tags that aligned to exactly the same unique position and strand in the maizeterm·dom reference genome (B73 RefGen v1) and that also co-segregated with the physically closest SNP ( p <0.001) were merged into a single, biallelic markerterm·dom. These markersterm·dom were then re-tested for co-segregation with the physically closest SNP using Fisher's Exact Test ( p <0.001). Biallelic GBS markersterm·mth that passed the latter test were then incorporated into a high density, framework mapterm·mth and ordered according to their positions in the reference genome. To determine how many of the remaining presence/absence GBSterm·mth tags could be genetically mapped in maizeterm·dom, the binomial test of segregationterm·mth was repeated versus this high density framework mapterm·mth, with a threshold of p <0.0001.

Software for the sequence filtering and the mapping analysis was written in Java and is available on SourceForge ( http://sourceforge.net/projects/tassel/ ). This software is part of the TASSEL package but is not currently implemented in the TASSEL GUI.

Results Read quantity and quality Because we are interested in enabling genome wide association studies (GWAS)term·dom in maizeterm·dom, a species where linkage disequilibrium decays within two to three kbp [30] , we need to identify markersterm·dom that cover around one million genomic locations. For this reason we chose to use Ape KIterm·mth, a RE that should cut frequently in the maizeterm·dom genome because it recognizes a degenerate five bp DNA sequence. Of course, if less genome coverage is desired, the protocol can be easily modified to use enzymes that recognize six or more bp.

Out of 1,146,449 high-quality (filtered) reads from IBM parental line B73, 1,125,731 (98%) could be aligned with the maizeterm·dom genomic DNA sequence. BLAST results indicated that the majority of non-aligning reads represented maizeterm·dom sequences that were absent in the reference genome version used for the analysis (B73 RefGen v1). Of the 868,336 GBSterm·mth sequence reads that aligned perfectly to the maizeterm·dom genome (no mismatches), 673,354 (78%) mapped to single genomic locations while 194,982 (22%) mapped to multiple locations, 87,271 (10%) aligned to <5 sites while 107,711 (12%) mapped to ≥5 sites).

Sequences from the maizeterm·dom IBM mapping population (276 RILs) were collected in six lanes of a single flow cell at 48-plex. On average, 2,090 Mbp of DNA sequence data were collected per lane. From a total of 145,836,644 raw reads, 102,505,713 (70%) were “high quality reads” that passed the Illumina filter while 120,438,739 (83%) contained the barcodeterm·mth and the cut site and no “Ns” within the first 72 bases and were not adapter/adapter dimersterm·dom. This observation indicates that, overall, the Illumina filtering parameters seem to be underestimating read quality. Hence, to maximize the amount of useful data, we worked with raw reads from the qseq fileterm·mth. Very few adapter dimersterm·dom were detected (78,375 or <0. 06% of the raw reads). Out of the 25,397,905 rejected raw reads (17% of total reads), only 1,096,513 (0.75% of the total reads) were discarded solely because of “Ns” in the first 72 bases. The remainder of rejected reads was comprised of adapter/adapter dimersterm·dom and sequences that did not contain the barcodeterm·mth and cut site (24,223,017 reads). Of the 1,096,513 reads discarded solely because of “Ns”, only 36,009 contained a single “N” and only 21,005 contained two “Ns”, whereas the majority (1,039,499) contained more than two “Ns”.

From six sequencing lanes, we identified 809,651 sequence tagsterm·mth (at least five times) from one or both flanks of 654,998 of the 2.1 million Ape KIterm·mth cut sites lying within the single copy genomic fraction. These 0.81 million 64 bp sequence tagsterm·mth cover 51.8 Mbp, or 2.3% of the maizeterm·dom genome. We also observed that the Ape KIterm·mth libraries showed a preponderance of smaller fragments ( Figure 3 ), resulting from both a bias toward production of small fragments during the PCR step of library construction, and precise spatial requirements for optimal cluster formation on the sequencing flow cell (i.e., longer fragments produce diffuse clusters that result in low sequence signal intensity). Fragments under 64 base pairs result in the presence of either the common adapterterm·mth or an internal Ape KIterm·mth recognition sequence (from partial digestion or chimera formation) within 64 bases of the end of the barcodeterm·mth. These were fairly common; out of the 120,438,739 reads that passed our initial filtering criteria (possessing a bar code and cut site, etc), 20,585,840 (17%) were from fragments less than 64 bases in length. As noted in the Methods , these were truncated accordingly and filled to 64 bases with polyA.

Fragment size distributions of a virtual Ape KIterm·mth digest of the maizeterm·dom genome and unique (single-copy) Ape KIterm·mth sequence tagsterm·mth from the maizeterm·dom IBM mapping population. Note that for size bins on the x-axis “50” denotes a bin of size 1–50 bp, “100” denotes a bin of size 51–100 bp, etc. The reference genome employed for the maizeterm·dom virtual digest was B73 RefGen v1.

Barcodeterm·mth optimization Our preliminary studies using RE-digested DNA samples and a small number of same-length (8 bp) barcodesterm·mth showed a substantial decline in read quality in multiplexed sequencing reactions compared to control DNA or other barcoded DNA samples that did not include restriction sites (data not shown). This finding suggests that presence of the invariant restriction site recognition sequence at the beginning of each read (i.e., low 5′ sequence variation) caused base calling errors in subsequent cycles, probably because proper sequence phasing on the Illumina Genetic Analyzer is dependent on detecting 12 random nucleotides at the beginning of each sequence. The presence of the invariant RE cut-siteterm·dom at bases nine to 12, therefore, violates the phasing model assumptions ( http://www.illumina.com/Documents/products/technotes/technote_rta_theory_operations.pdf ). Incorporation of variable length barcodesterm·mth substantially improved base calling accuracyterm·mth, although it still appears that the Illumina algorithm sometimes underestimates read quality. Reads that did not pass the Illumina filter sometimes perfectly matched a 64 base tag that was segregating in our mapping populations.

Sample representationterm·dom The six lanes of the maizeterm·dom IBM population sequencing run yielded 120,438,739 GBSterm·mth reads that contained the barcodeterm·mth and the Ape KIterm·mth cut-site (or 20,073,123 reads per lane). On average, 436,372 reads were produced per DNA sample and 95% of samples generated at least 125,000 reads. Evenness of sample representationterm·dom among the maizeterm·dom IBM RILs was acceptable but not optimal. In our best lane from the IBM flow cell, the coefficient of variationterm·mth (cv = standard deviation/mean) for the number of reads containing the appropriate barcodeterm·mth and the cut site was roughly 43% among samples and, among the six lanes, 39.8% of the variance was attributed to DNA sample. Subsequent adjustments to our robotic liquid handling protocols, however, have resulted in greater evenness among samples ( Figure 4 ). Regardless of the disproportionate sample representationterm·dom, we were still able to map a minimum of 90,000 sequence tagsterm·mth in the poorest performing IBM samples.

Coefficient of variationterm·mth of GBSterm·mth reads per sequencing channel for sequential sequencing runs. Each flow cell comprised 6 or 7 sequencing channels. Large boxes represent the standard deviation of the number of reads per sample; whiskers denote minimum and maximum values; small squares are the median values; and lines extending across the boxes are the means for each run. Flow cells are ordered sequentially by run date; number 1 is the first sequencing run and number 11 denotes the last run. The GBSterm·mth read data from the maizeterm·dom IBM population is contained in flow cell 1. The large variation in reads per sample from this flowcell was due to inconsistent pipetting during robotic liquid handling. Subsequent adjustments to our robotic protocols improved evenness among samples (see flowcells 2–11).

Preliminary results for barleyterm·dom were slightly better with respect to uniform sample representationterm·dom ( Figure 5 ). The one channel of the sequencing run produced 27.5 million reads. On average, 427,130 reads were produced per DNA sample (minimum = 145,648; maximum = 643,631) with a coefficient of variationterm·mth (cv) of 23% ( Figure 5 ).

Distribution of reads across 43 barcoded samples in a single flow cell lane for the Oregon Wolfe Barleyterm·dom population.

Mapping and SNP validation Analysis of the maizeterm·dom IBM population provided a preliminary evaluation of the genetic value of multiplex GBSterm·mth skimming. Overall 25,185 biallelic 64 base tags were genetically mapped to their physically closest anchor SNP. No corresponding alternate allele was found for an additional 584,119 tags. By treating these as dominant data (i.e., either present or absent in each RIL), 167,494 could be placed upon the framework mapterm·mth of 25,185 biallelic sequence tagsterm·mth based upon segregation. Alignment to the reference genome detected unique physical positions for 133,129 of the dominant markersterm·mth, 90.8% of which agreed with the genetic positionsterm·dom.

After filtering for tags present in at least 20% of the lines, 2.1 million unique barleyterm·dom tags were retained. These tags were mapped to the OWB framework mapterm·mth of 2,382 markersterm·dom [31] by considering tags as dominant markersterm·mth and anchoring the tags using the reference mapterm·mth. Prior to mapping, the genetic map was collapsed to retain only markersterm·dom that contained unique linkage information in the subset of 43 lines (i.e., SNPs at the same map position were removed) leaving 436 biallelic markersterm·dom. I n all, we mapped 24,186 sequence tagsterm·mth onto the barleyterm·dom genetic map. To determine the utility of using the sequence tagsterm·mth as genetic markersterm·dom, cross validationterm·mth was conducted for one of the OWB lines (OWB003). Tags were mapped without OWB003 and coded according to whether they were present in either the dominant or recessive OWB parent. A graphical genotype of the excluded, control line, OWB003, showed almost perfect agreement between the reference markersterm·dom and GBSterm·mth regarding chromosome segment parent of origin ( Figure 6 ). A cross validationterm·mth error was scored if a previously mapped SNP and a GBSterm·mth tag disagreed on parent of origin. GBSterm·mth markersterm·dom occurring near the OWB003 recombination break points cannot be unambiguously assigned and were excluded when determining genotyping accuracy. Of the 4,596 mapped GBSterm·mth reads present in OWBOO3, 4,533 (99%) identified the correct parent of origin.

Barleyterm·dom GBSterm·mth markerterm·dom validation using a single DH line (OWB003). Upright triangles denote positions of markersterm·dom on the reference genetic map and downward triangles indicate GBSterm·mth reads mapped in this study. Multiple sequence reads are stacked and colors indicate chromosomal segments in OWBOO3 originating from dominant (blue) or recessive (red) parental lines.

Discussion GBSterm·mth offers an alternative to complex, expensive protocols The value of reducing genome complexityterm·dom with REs coupled with multiplex NGS for high-density SNP discovery and genotyping was originally demonstrated with restriction site-associated DNA (RAD) tags [20] . In the RAD protocol, genomic DNA is digested with a six to eight base-cutter RE and a barcoded adapter is ligated to compatible sticky ends. For multiplex sequencingterm·mth, DNA samples, each with a different barcodeterm·mth, are pooled, randomly sheared, size selected (300–700 bp), and a second adapter is ligated after polishing and filling ends. The RAD library preparation procedureterm·mth is substantially more complicated than the one presented here. In addition to its simplicity (no fragment size selection and few enzymatic and purification steps), our protocolterm·mth is time and cost efficient through its use of a single well for genomic DNA digestion and adapter ligation. These processes can be done in the same buffers so that no additional transfers are needed. Currently, the favored enzymes ( Ape KIterm·mth and T4 ligase) do not have complementary temperature regimes so simultaneous digestion and ligation is not possible unless we substitute an expensive, thermostable ligase.

Recently, a multiplex NGS protocol appropriate for Drosophila and other small-bodied species has been published [32] . Although this protocol is similar to the one reported here, it still includes a fragment size selection step and variation in the number of reads between samples was considerably higher (cv = 89%) [32] than what we observed (cv = 23–43%). Clearly for any multiplex sequencingterm·mth protocol, accurate quantification of high molecular weight DNAterm·mth remains a procedural bottleneck, and is the most likely source of sample-to-sample variation in sequence coverage. DNA quantification using intercalating dyes and spectrophotometry give correlated but not very consistent results. Developing a more precise, cheap, high-throughput DNA quantification protocol, therefore, remains an area where this method can be improved. The performance of the barcodesterm·mth, themselves, is less problematic because, over time, poor performing barcodesterm·mth can be identified and removed from the protocol.

GBSterm·mth does not use 5′-phosphorylated adapters Unlike other protocolsterm·mth, GBSterm·mth employs unmodified adaptersterm·mth (i.e., without 5′-phosphate groups). As a result, only one adapter strand is covalently bound to the ends of restriction fragments. As long as DNA samples are not denatured prior to the pooling and PCR steps, however, Taq polymerase rapidly fills the 3′ recessed ends in the presence of dNTPs. End filling occurs either by immediate displacement of the non-ligated adapter strands at low temperatures (during the assembly of PCR reactions) or following the early dissociation of short, non-ligated strands during the initial heating step of the PCR [33] . Use of unphosphorylated adapters has the added benefits of destabilizing formation of adapter dimersterm·dom during library preparation and reducing reagent costs.

GBSterm·mth does not employ divergent “Y” adapters Standard libraries for Illumina sequencing are prepared by ligating a single “Y” or “forked” adapter to both ends of genomic DNA fragments [26] . These adapters, made by annealing oligonucleotides with both complementary and non-complementary sequences, have, at one end, a region of double stranded DNA that is required for T4 DNA ligase to join adapters to genomic DNA. The other end of the adapter is comprised of single stranded, divergent sequences that serve as binding sites for a pair of primers that, after PCR, generate DNA fragments that have different adapter sequences on each end.

The GBSterm·mth protocol employs two different double stranded adapters (barcodeterm·mth and common) that are ligated simultaneously to restriction fragments with “sticky” ends. This means that any combination of adapters (barcodeterm·mth/common, barcodeterm·mth/barcodeterm·mth, or common/common) may be joined to genomic DNA fragments. Because same-ended DNA strands bind to the flowcell but do not produce DNA sequence on the Illumina platform (those having only barcode adapterterm·mth sequences are cleaved from the surface prior to reverse-terminator sequencing and those with two common adapterterm·mth ends lack a binding site for the PE1 sequencing primer), one might predict a reduction in the number of reads per lane with the GBSterm·mth protocol. Read numbers, however, have equaled or exceeded the specifications of the DNA sequencing instruments (both the Illumina GAII and, more recently, the Illumina Hi-Seq). The lack of impact on read number is most likely because same-ended fragments occupy little physical space on the flowcell. This reasoning seems counterintuitive in light of the fact that half the fragments applied to the flowcell are predicted to have the same adapter ends. Because same-ended strands are able to utilize only one of the two oligonucleotides that coat the flowcell surface during bridge amplification [34] , however, formation of DNA clusters (colonies) is inefficient. As a result, these slow-growing, sparse clusters are rapidly overrun by “normal” colonies derived from DNA strands with different ends that prime off both oligonucleotides.

GBSterm·mth accesses regulatory regionsterm·dom and sequence tagterm·mth mapping requires no reference genome As more information is gathered, it is becoming apparent that regulatory regionsterm·dom controlling the expression of plant genes responsible for agronomically important phenotypes are often located in non-coding DNA. For example, regulatory regionsterm·dom of maizeterm·dom genes vgt1 [35] , tb1 [36] and b1 [37] are located 60 to 150 kb from the structural gene. Therefore, systematic discovery and mapping of genetic diversity should not be limited to coding regionsterm·dom. In this sense, the GBSterm·mth procedure allows access to any sequence within “low-copy” genomic regionsterm·dom, including transposable elements and repeat regions that have not proliferated extensively.

Another advantage to the GBSterm·mth approach is that a reference genome need only be developed neighboring the restriction sites, and this can be done in the process of sample genotyping. In such cases, the consensus of the read clustersterm·mth across the sequence tagged sites becomes the reference. Alternatively, for kinship analyses and genomic selectionterm·dom in the absence of a reference genome, the tags can simply be treated as dominant markersterm·mth. While not addressed here, there has also been tremendous progress in the imputationterm·mth of missing data. In biparental mapping populations of species with a reference genome, this can be done with extremely high accuracy [5] , and even in more diverse material, imputationterm·mth accuracies over 99% permit low coverage. In the case of maizeterm·dom, we envision performing whole genome sequencing on a few thousand lines and then projecting their polymorphisms onto hundreds of thousands of additional lines via GBSterm·mth and imputationterm·mth.

We have shown that for an expenditure of $8,000 (USD), approximately 200,000 maizeterm·dom markersterm·dom can be identified and mapped in a very short time. With the methods outlined here, we can process 336–672 samples (48 or 96 samples per channel ×7 channels per flow cell) simultaneously. Multiplexes up to 384 per lane (2,688 samples per sequencing run) are becoming possible as read density improves on sequencers. The economy of scale associated with these improvements is rapidly pushing genotyping below $20 per sampleterm·dom. Projected gains in the near future could result in a further four to five fold reduction to $5 or less per sample. Soon, plant breeders may conduct genomic selectionterm·dom on a novel germplasm or species without first having to develop any prior molecular tools, or conservation biologists may determine population structureterm·dom without prior knowledge of the genome or diversity in the species. These exciting new avenues for applying GBSterm·mth to breeding, conservation, and global species and population surveys are now poised to become an indispensable component of future biology.

Supporting Information Optimizing GBSterm·mth for Other Species.

(DOCX) Click here for additional data file.

Experion® output showing fragment size distribution of an “unoptimized” maizeterm·dom GBSterm·mth library. Note that the x-axis denotes seconds (elution time) and not fragment size (bp). Two discrete peaks are observed, the primer dimer peak at around 45 seconds (∼70 bp), the adapter dimerterm·dom peak between 51 and 52 seconds (∼128 bp). These are followed by a broad peak, the GBSterm·mth library, occurring between 55 and 65 seconds.

(TIFF) Click here for additional data file.

Experion® output showing fragment size distribution of a GBSterm·mth library after adapter amounts were optimized. Note that the adapter dimerterm·dom peak has disappeared.

(TIFF) Click here for additional data file.

GBSterm·mth barcodeterm·mth sequences for Ape KIterm·mth adapters.

(DOCX) Click here for additional data file.

Acknowledgments Kate Guill and Mike McMullen provided the IBM DNAs and the genotypes of previously mapped 644 SNPs in the IBM population.

Footnotes References Associated Data Supplementary Materials Optimizing GBSterm·mth for Other Species.

(DOCX) Click here for additional data file.

Experion® output showing fragment size distribution of an “unoptimized” maizeterm·dom GBSterm·mth library. Note that the x-axis denotes seconds (elution time) and not fragment size (bp). Two discrete peaks are observed, the primer dimer peak at around 45 seconds (∼70 bp), the adapter dimerterm·dom peak between 51 and 52 seconds (∼128 bp). These are followed by a broad peak, the GBSterm·mth library, occurring between 55 and 65 seconds.

(TIFF) Click here for additional data file.

Experion® output showing fragment size distribution of a GBSterm·mth library after adapter amounts were optimized. Note that the adapter dimerterm·dom peak has disappeared.

(TIFF) Click here for additional data file.

GBSterm·mth barcodeterm·mth sequences for Ape KIterm·mth adapters.

(DOCX) Click here for additional data file.