Discussion

The ability to rapidly discover novel nucleotide variation in germplasm collections and mutant populations is a fundamental tool for functional genomics and breeding of plants and animals. For example, large germplasm collections have been maintained for many species, but limited resources mean little has been done to characterize gene coding and regulatory sequence variation that may prove important for unlocking their full potential. Further, efficient methods for accurate discovery and cataloguing of nucleotide variation can be important to measure and protect biological diversity in developing countries (Gepts 2004). For functional genomic studies, screening of a TILLING population for the recovery of induced mutations predicted to alter gene function typically requires interrogation of thousands of individuals. While reduced-representation genome methods such as exome capture sequencing have been described to recover rare single nucleotide variants in large genome species such as wheat, the approach remains cost-prohibitive for smaller scale projects and for understudied species where research funds are limited (Krasileva et al. 2017). Amplicon sequencing provides a low-cost alternative that is advantageous in that it is highly flexible with regard to population size and the choice of sequence regions for evaluation.

The flexibility and throughput of amplicon sequencing is increased through fragmentation of PCR products. More bases can be interrogated from a PCR product, and fragmentation into small pieces allows use of a range of different short-read sequencing platforms. Fragmentation of low molecular weight PCR products, however, is potentially more challenging when compared to working with genomic DNA where starting fragment sizes can exceed 20,000 bp when using standard genomic DNA extraction protocols (Huynh et al. 2017). For example, recovery of a 500 bp fragment from a 1000 bp PCR product requires that a single double strand break be induced at the mid-point of the amplicon. Owing to this, optimal conditions for sonication-based fragmentation of a single amplicon were evaluated. Trials showed that sonication for 30 sec with a peak power of 50 W, Duty factor of 40 and 200 cycles per burst produced a smooth distribution of fragment sizes that peaked near 400 base pairs, a suitable size for Illumina library preparation and sequencing. While different sized amplicons may require subtle optimizations, it was reasoned that sonication at this setting may be suitable for a complex pool of PCR products of varying sizes. A pool of amplicons ranging between 670 and 1513 bp was next evaluated because this is the size range of amplicons used in traditional gel-based TILLING approaches that have served as a benchmark for mutation discovery in pools (Till et al. 2003, 2004; Talamè et al. 2008; Gottwald et al. 2009; Bovina et al. 2011). Conditions were sought to allow higher throughput screening of existing TILLING populations where gene targets and oligonucleotide primers have been previously validated.

Evaluation of fragment size distribution of the libraries prepared from pooled PCR products showed a broad distribution of fragments, suggesting suitability for sequencing. Analysis of post-sonication fragment sizes after library preparation, and also of insert sizes post-mapping, showed a broad distribution of fragment sizes with median insert size of 348.23 bp, close 350 bp, one of two optimal insert sizes for Illumina library preparation. Over-accumulation of DNA fragments at specific molecular weights was observed. However, this did not appear to affect the experiment as all pools produced high base and mapping quality with the majority of amplicons in pools producing high read coverage. Over-accumulation of specific molecular weights may indicate that some sequence contexts or fragment sizes are less likely to experience sonication-induced cleavage. The GC content of the amplicon may influence this, as a weak negative correlation between GC content and mean coverage was observed. Internal sequences of DNA fragments longer than the sequencing read length will go unsequenced affecting coverage across an amplicon. However, a weak positive correlation between GC content and amplicon size was also observed, and coverage differences may be explained by the method used for amplicon quantification (see below). Further tests are required to determine if there is any dependence on the concentration of input PCR product, as the phenomenon of accumulation of specific molecular weight fragments was most pronounced in the highest concentration pool number 24. It is interesting to note that while read coverage is generally even across fragmented amplicons, there is a consistent pattern of higher coverage and reduced base quality near but not directly at amplicon termini. This pattern is expected in paired end reads if there is a lower probability for fragmentation of DNA sequences near the termini of the PCR amplicon. Importantly, with the exception of the increases near the amplicon termini, read coverage is consistent across the length of amplicons. A weak correlation was observed between fragment size and median sequencing coverage in the experiment. While the sonication procedure cannot be ruled out as a contributing factor, it is also likely that this observation reflects errors in quantification of PCR products prior to pooling. Intercalating DNA binding dyes will produce brighter band images (with a higher pixel density) in larger DNA fragments compared to smaller fragments of the same molar amount. The quantification used for the barley experiment did not correct for amplicon size. Therefore concentrations were over-estimated for larger fragments and therefore these amplicons were over-diluted compared to smaller fragments. While slightly underrepresented, quantity and coverage from the larger amplicons was suitable for mutation discovery. Based on the consistency of coverage across the length of amplicons, it was concluded that the mean coverage per amplicon in each pool is a useful metric to evaluate the performance of the entire population. The population consisted of 1536 data points (32 amplicons × 48 pools). In total, 2.08% of data points had coverage values less than 5% of the mean for the amplicon in all pools. These were considered assay failures even though some high quality reads were produced. Assay failures measured in this way were slightly lower in cassava and tomato non-fragmented amplicon experiments performed on the same MiSeq sequencer (0.5 and 1.19%, respectively). Assay failures may arise due to a combination of variables including genomic DNA quality, genomic DNA quantification, PCR quality, PCR quantification and library quantification (see below). Interestingly, PCR product evaluation and quantification methods were different in all three experiments. Every tomato PCR product was quantified by capillary electrophoresis using the Fragment Analyzer. This proved costly, and so for the cassava experiment all samples were analyzed using 96 well e-gels and digital image quantification using Image J (Duitama et al. 2017; Huynh et al. 2017). For barley, the procedure used in cassava was streamlined by replacing digital image quantification with a rapid visual estimation to determine average yield in nanograms for each amplicon. The method of PCR product quantification trends with failure rates based on sequencing coverage, and therefore may be related. Nevertheless, assay failure is low in all experiments and from this we conclude that using a fragmented PCR approach does not significantly add to assay failures.

Comparisons between sequencing metrics from fragmented and non-fragmented PCR products allow further evaluation of the sonication-based approach. Data from direct sequencing of non-fragmented PCR products pooled in a similar fashion (Figure S5) shows both coverage and base quality drop toward the center of amplicons that are near 600 base pairs. In contrast, coverage increases in the center of shorter amplicons as expected when applying paired-end sequencing. Based on this data, we conclude that the sonication-based approach produces more consistent data independent of starting amplicon size. This allows for more precise experimental design, taking into consideration genomic DNA sample pooling and amplicon pooling so that required depth of coverage can be obtained for accurate variant calling.

There are many parameters in addition to amplicon fragmentation that are important for a successful assay. These include proper quantification and pooling of genomic DNA, proper quantification and pooling of PCR products, proper quantification and pooling of sequencing libraries, and producing sufficient read coverage to ensure recovery of rare mutations based on how many genomic DNA samples are pooled together. It is difficult to evaluate if a DNA quality issue results in a subset of genomic DNAs in a pool underperforming in PCR. All other parameters can be investigated. Quantification of PCR products is evaluated through coverage variations within the same amplicon between different pools, and also through coverage variations for all pools combined calculated between different amplicons. Ideally, all PCR products in all pools will be represented at the same molarity and thus produce the same sequencing coverage. This may be impossible to achieve. In the barley experiment, maximum coverage within the same amplicon ranged between 2 and 11 fold higher than the mean (Table S12). Coverage between different amplicons ranged between 2092.7 and 13445.8 (6.42 fold). This is in comparison to a 4.94 fold difference between the minimum and maximum coverage between amplicons in the tomato TILLING experiment. Thus, the careful quantification of each PCR product used for tomato TILLING, which utilized robotic liquid handling and capillary analysis of each non-pooled PCR product, may provide only incremental assay improvements. In addition to PCR concentration, variation in sequencing library quantification will affect the mean coverage of each pool in an experiment. Library quantification in the barley experiment was performed using fluorimetry. One noticeable error was observed whereby pool 24 was highly over-represented (mean coverage of 52170.41) compared to the next highest pool (mean coverage 8854.62, pool 31). Library quantification using fluorimetry allowed rapid quantification compared to quantitative PCR and the over-representation of pool 24 is assumed to have resulted from a human error. Further tests are required to determine if library to library mean coverage can be made more similar by adding additional quantification methods such as qPCR. Nevertheless, assay success can be measured as having achieved sufficient quality and coverage for pools and amplicons in order to recover mutations. Expected variations can be built into the experimental design when choosing the number of PCR amplicons to sequence and the level of genomic DNA pooling to employ. With regards to coverage of induced mutations, 88% of the experiment produced sufficient data for recovery of heterozygous mutations with a read coverage of the mutant allele of at least four. This coverage can be increased by increasing the number of sequencing reads in the experiment.

Numerous DNA sequencers can potentially be used for pooled amplicon experiments when amplicons are fragmented. Quality and throughput can vary between different platforms and between different sequencer models of the same platform, as reported for Illumina sequencers (Liu et al. 2012). In theory, any platform is suitable provided sufficient coverage and quality is achieved. MiSeq 2x300PE sequencing was used in this study to provide direct comparison with previous non-fragmented amplicon experiments run on the same machine. Higher throughput and lower per-base costs are possible on other sequencing platforms such as the Illumina HiSeq system. Cost estimations can also be considered to further evaluate the utility of sonication of larger amplicons vs. direct sequencing. For simplicity, quality and coverage are assumed to be identical for a non-fragmented 600 bp amplicon vs. a fragmented 1200 bp amplicon. When considering major consumables and assays employing capillary quantification of all PCR products, it is estimated that it will cost approximately 1.9x more to screen two 600 bp amplicons vs. one 1200 bp amplicon that has been sonicated (Table S16). Cost estimations do not take into account the extra time for PCR amplification and FA machine run times when performing more PCR reactions producing shorter amplicon sizes. Comparison of barley, tomato and cassava data suggests that capillary quantification may only provide small assay improvements. Without this step, sequencing of 600 bp amplicons using 2x300PE sequencing is 1.3x more costly than sequencing 1200 bp amplicons (Table S16). Additional time and cost savings may be achieved by using amplicons larger than those reported here. For example, while gene sizes vary, median gene length in plants such as barley can be thousands of base pairs (Wang et al. 2017; Mascher et al. 2017). An elegant experimental design would be for one amplicon to cover an entire coding sequence. Long PCR, however, can be more challenging to optimize, and careful testing should be performed to avoid unforeseen biases and to ensure that PCR amplification produces and equal representation of molecules from complex pooled genomic DNAs. Errors from PCR must also be considered. Previous studies using traditional gel-based approaches showed using Taq polymerase with amplicon sizes up to 1500 bp resulted in less than a 5% false positive error rate (Till et al. 2006a, 2010). PCR based errors are further mitigated when applying a three-dimensional genomic DNA pooling strategy whereby each genomic DNA is represented in three unique pools and true mutations are found only once in each pooling dimension (Tsai et al. 2011). When applying this criterion to mutation calling the likelihood of false positive errors due to PCR amplification of genomic DNA or through post-ligation library amplification is reduced. Further, library preparation kits that do not employ post-ligation amplification are available if errors from library preparation are a concern. However, more studies are needed before direct comparisons between traditional enzymatic-based TILLING assays and sequencing-based assays can be made. Higher pooling in sequencing assays (256x vs. 8x for enzymatic assays) means that rare errors will be more likely to be detected.

An optimal experimental design for sequencing of fragmented pooled amplicons will balance assay throughput, cost, and accuracy. The ability to maximize recovery of true positives while minimizing false positives is important. The choice of bioinformatics tools is a critical component. Previous studies in tomato TILLING by Sequencing tested six different mutation calling tools and estimated a range of accuracy from 89.33 to 29.33 (Gupta et al. 2017). In that study, the CAMBa tool that was purpose built for TILLING performed the best. In the present study, CAMBa successfully identified all three positive controls, but only the homozygous mutation was identified in all three pools using GATK HaplotypeCaller with sample ploidy set at either 600 or 2000. The homozygous mutation is present in DNA pools at an expected frequency of 0.0039 while the heterozygous mutation is expected at 0.00195. Using the amplicon coverage values for the pools, the homozygous mutation is expected to be sequenced 14x in pool 11, 35.65x in pool 19 and 34.26x in pool 41. The heterozygous mutations are expected at 4.98x, 7.91x, 14.26x, 11.81x, 8.04x, and 20.15x respectively, in pool numbers 15, 20, 37, 14, 26, and 46. Lower representation of heterozygous mutations may have contributed to the failure to call the mutation with HaplotypeCaller. Analysis of novel putative mutations identified in the mutant screen is ongoing. True mutations resulting from this work will be used to further evaluate the performance of mutation discovery software. Variant calling tools are constantly evolving and their improvement will likely lead to improvements in pooled amplicon assays such as TILLING.

We conclude that fragmentation of pooled amplicons by ultrasonication provides a suitable method for producing even coverage sequencing for targeted recovery of nucleotide variation in large populations of samples. Alternative methods for amplicon fragmentation, such as use of dsDNA fragmentase have been applied for TILLING with genomic pooling of 64-fold (Burkart-Waco et al. 2017). This suggests that multiple methods may be considered for experiments employing complex pools of PCR products. In theory, any method that can produce a distribution of fragments with optimal size ranges for the chosen library preparation method should be suitable. Further tests are needed, however, to determine if alternative fragmentation methods will produce similar results in highly pooled samples. Given all of the experimental variables in a pooled amplicon experiment that can affect variant discovery, and the fact that variant calling algorithms are still evolving, careful testing of amplicon fragmentation, amplicon pooling, genomic DNA pooling and bioinformatics tools using positive controls is advised prior to choosing an experimental design. In addition to recovery of induced mutations it is envisioned that combining sample pooling with amplicon fragmentation can improve screening efficiencies for natural allelic variations. For example, the same approach used for barley mutants can be applied for sequencing of germplasm banks where rare alleles may provide an important resource for understanding gene function and allow for the genetic improvement of domesticated and semi-domesticated species.