Abstract
Improvements to massively parallel sequencing have allowed the routine recovery of natural and induced sequence variants. A broad range of biological disciplines have benefited from this, ranging from plant breeding to cancer research. The need for high sequence coverage to accurately recover single nucleotide variants and small insertions and deletions limits the applicability of whole genome approaches. This is especially true in organisms with a large genome size or for applications requiring the screening of thousands of individuals, such as the reverse-genetic technique known as TILLING. Using PCR to target and sequence chosen genomic regions provides an attractive alternative as the vast reduction in interrogated bases means that sample size can be dramatically increased through amplicon multiplexing and multi-dimensional sample pooling while maintaining suitable coverage for recovery of small mutations. Direct sequencing of PCR products is limited, however, due to limitations in read lengths of many next generation sequencers. In the present study we show the optimization and use of ultrasonication for the simultaneous fragmentation of multiplexed PCR amplicons for TILLING highly pooled samples. Sequencing performance was evaluated in a total of 32 pooled PCR products produced from 4096 chemically mutagenized Hordeum vulgare DNAs pooled in three dimensions. Evaluation of read coverage and base quality across amplicons suggests this approach is suitable for high-throughput TILLING and other applications employing highly pooled complex sampling schemes. Induced mutations previously identified in a traditional TILLING screen were recovered in this dataset further supporting the efficacy of the approach.
Next generation sequencing (NGS) techniques have had a profound impact on biological research. While technologies continue to advance, whole genome approaches remain costly for projects involving the analysis of species with large genomes or those involving the interrogation of many individuals. A variety of reduced representation genome sequencing approaches have been described to circumvent this issue (Hirsch et al. 2014). One powerful approach to evaluate sequence variation in targeted regions is by sequencing PCR amplicons. Amplicon sequencing has been applied for the discovery and characterization of both natural and induced mutations in plant and animal populations (Rigola et al. 2009; Tsai et al. 2011; Pan et al. 2015; Duitama et al. 2017). For example, sequencing has been used to increase throughput of mutation discovery for the reverse-genetics technique known as TILLING (Targeting Induced Local Lesions IN Genomes)(McCallum et al. 2000; Till et al. 2018). For TILLING, mutant populations are typically created using mutagens such as ethyl methanesulfonate (EMS) or a combination of sodium azide and N-Nitroso-N-methylurea (MNU) that induce primarily single nucleotide variants (SNV). While the frequency of mutations varies between species and ploidy level, mutation densities in diploid species have been reported to range between 1 mutation per 150,000 and 1 mutation per 700,000 base pairs (Jankowicz-Cieslak and Till 2015). TILLING screens typically aim to recover multiple mutations from a single target to increase the chances of recovering deleterious alleles. Therefore, a typical TILLING assay involves the screening of 3000 or more mutant individuals. To increase throughput, genomic DNAs are pooled together prior to mutation screening. The application of new sequencing technologies to improve TILLING throughput has been termed TILLING by Sequencing (Tsai et al. 2011). In addition to pooling genomic DNA, multiple amplicons are produced from each gDNA pool and combined prior to massively parallel sequencing to further improve throughput. The quantitative nature of NGS methods means that rare induced and natural SNV mutations can be effectively recovered from pools of over 200 individuals (Pan et al. 2015; Duitama et al. 2017). The most common sequencing platform used for TILLING by Sequencing is short read sequencing-by-synthesis from the company Illumina (Tsai et al. 2011; Guo et al. 2015; Pan et al. 2015). Similar approaches have been used in population genetics studies to identify rare alleles in large germplasm collections (Duitama et al. 2017).
Direct sequencing of PCR amplicons is potentially advantageous in that an even read coverage across the target can theoretically be achieved, thus providing consistent recovery of sequence variants. This approach is disadvantageous, however, due to the fact that the length of sequenced amplicons is limited to the maximum read lengths of the sequencer used (e.g., 500-600 bp for the Illumina MiSeq)(Pan et al. 2015; Duitama et al. 2017). Thus, multiple amplicons are required to screen an entire gene, necessitating extra liquid handling steps and also additional work in adjusting PCR amplicon concentrations to a similar level prior to sequencing. An alternative approach is to produce longer PCR products and to fragment them prior to library construction. Several approaches for fragmentation have been described such as nebulization, enzymatic cleavage, and ultrasonication (Head et al. 2014). Most approaches are developed for the fragmentation of genomic DNA samples that contain high molecular weight molecules. We describe here the optimization of ultrasonication for PCR products amplified from the genomes of Coffea arabica (Arabica coffee) and Hordeum vulgare (barley) ranging between 670 and 1513 bp. We further show suitable coverage can be achieved in a pool of 32 distinct PCR products generated from PCR amplification of pools of 256 genomic DNA samples prepared from chemically mutagenized barley. This suggests that the approach can be easily adapted for pooled amplicon sequencing in different species.
Materials and Methods
Genomic DNA and PCR
Genomic DNA from Coffea arabica was provided by Margit Laimer. PCR primers were designed using the Primer3 program with parameters previously developed for TILLING assays(Till et al. 2006b). Primer sequences were TCGATTCGATTCGTTGACACCCCTA and TGGATGATGGATGGGAATGTGGTTC. PCR was performed as previously described (Till et al. 2015). Amplification of a single PCR product was assayed using agarose gel electrophoresis and confirmed via capillary electrophoresis (Figure 1A).
Comparison of different fragmentation parameters of a single PCR amplicon. Panel A shows the non-fragmented control of the single 964 bp PCR product as assayed by capillary electrophoresis. Panel B shows the results of test parameter 1 in Table 1. This is an example of incomplete fragmentation as the original 964 bp PCR product can be observed at a high concentration. Panel C and D are technical repeats showing results of fragmentation using test #13 parameters. The y axis in all graphs shows relative abundance of product in relative fluorescence units (RFU). The x axis shows molecular weight in base pairs. LM and UM are the lower and upper marker, respectively.
Genomic DNA from chemically mutagenized barley was prepared using a modified CTAB protocol (Szarejko et al. 2017). DNA samples were adjusted to a similar concentration and pooled. Two hundred and fifty-six samples were combined in each pool. Forty-eight unique pools were prepared. Primers were designed to target regions of the H. vulgare genome as previously described (Table S1)(Slota et al. 2017). PCR primers amplify part of a coding sequence and thus target names are assigned for each amplicon that is unqiue from gene annotations. PCR was performed on pooled genomic DNA samples using 32 unique primer pairs. Individual PCR reactions were performed for each primer pair. The concentration of each PCR product was estimated by running products on a 96-well eGel system with lambda DNAs standard of known concentrations (Huynh et al. 2017). Average concentration was estimated visually and all samples adjusted to a final concentration of 200 ng in TE buffer. PCR products were then pooled together. DNA sequences of each resulting amplicon were prepared from the barley reference genome and combined into one multifasta file for downstream bioinformatic analysis (Gupta et al. 2017). GC content for each amplicon was extracted from the multifasta file using the Emboss infoseq tool (Carver and Bleasby 2003). Average GC content for barley was calculated from the coding sequences from the International Barley Genome Sequencing Consortium version 2 assembly (IBSC v2, INSDC Assembly database version 95.3) (International Barley Genome Sequencing Consortium et al. 2012).
GC percentage for each sequence was extracted using the Emboss infoseq tool. Minimum and maximum GC percentage was extracted using the following command: awk ′BEGIN{first = 1;} {if (first) {max = min = $3; first = 0; next;} if (max < $3) max=$3; if (min > $3) min=$3;} END {print min, max}′. Average GC was calculated using the following command: awk ′{sum = sum+$3 ; sumX2+=(($3)^2)} END { printf “Average: %f. Standard Deviation: %f \n”, sum/NR, sqrt(sumX2/(NR) - ((sum/NR)^2))}′ .
Sonication of PCR products
Sonication of PCR products was performed using a Covaris M220 ultrasonicator with microTUBE AFA bead split tubes (coffee amplicons) or microTUBE AFA FiberPre-slit (barley amplicons) in 60 μl volumes except where indicated with parameters adjusted according to Tables 1, S2 and S3. Fragmentation of PCR products was assayed using a Fragment Analyzer with the low sensitivity 1kb separation matrix with 30 cm capillaries (cat #DNF935). Analysis of data were performed using the PROSize Data Analysis Software. Performance of sonication of pooled PCR products was independently tested in 48 pools. Pooling of genomic DNA for each of the 48 pools was done such that each individual genomic DNA is represented in triplicate in the experiment (in three different pools). A total of 4096 DNAs from unique mutant lines are represented in the 48 pools.
Sequencing and data analysis
Quantified PCR products were normalized to approximately 200ng and pooled together according to a three-dimensional pooling scheme whereby samples were arrayed in a 16 × 16 × 16 grid for pooling (Gupta et al. 2017). A total of 48 pools were created for library preparation. Sequencing libraries were prepared using TruSeq Nano DNA HT Library Prep Kit (Illumina, Catalogue # FC-121-9010DOC) with 200 ng of starting pooled PCR products according to manufacturer’s recommendations. Libraries were quantified using a Qubit fluorimeter, concentrations normalized to a common molarity and pooled as previously described (Datta et al. 2018). Sequencing was performed on an Illumina MiSeq using 2x300 PE chemistry according to manufacturer’s protocol. Reads were mapped to amplicon sequences during the sequencing run using BWA with the MiSeq BWA default settings (Li and Durbin 2009). Sequencing coverage and mapping quality were calculated using Qualimap version v.2.2.1 with default settings except that the number of windows was set to amplicon length to retrieve per base metrics used in Figure 3 (Okonechnikov et al. 2016). Base qualities were prepared using pysamstats and the baseq feature (https://github.com/alimanfoo/pysamstats). Mean coverage data per amplicon for all 48 pools was extracted from the Qualimap output. Principle component analysis of all 48 pools was produced using the multi-sample BAM QC feature of Qualimap. Correlation coefficients were calculated using the LibreOffice CORREL function. Coverage data used to calculate mean coverage for tomato amplicon sequencing comes from Gupta et al. (2017). Coverage data used to calculate mean coverage for cassava sequencing was extracted using nucleotide number binning rather than amplicon name as data were mapped to a contiguous sequence rather than individual amplicons (Duitama et al. 2017). Data for barley, tomato and cassava were collected on the same MiSeq instrument using 2x300 PE sequencing.
Screening of data for previously identified SNP mutations was performed using the CAMBa pipeline for three dimensionally pooled samples as previously described (Tsai et al. 2011; Gupta et al. 2017). The HaplotypeCaller tool from GATK 4.0.10.0 was used as a second method to call SNPs (McKenna et al. 2010). Processing was carried out in parallel using the following command: ls *.bam | parallel java -jar gatk-package-4.0.10.0-local.jar -Xmx32g HaplotypeCaller -I {} -R reference.fa -O {.}.vcf–sample-ploidy N, where sample ploidy N was set to either 600 or 2000 in different analyses (Tange 2011). Resulting VCF files were first filtered to remove common variants appear above 5% allele frequency in pools (bcfftools view -i ‘AF < 0.05’). True mutations should be present in only one sample in a pool and therefore the expected frequency is 1/256 for homozygous mutations and 0.5/256 for heterozygous mutations. The 16 VCF files from a single dimension pool where then grouped together for further analysis. True induced mutations should appear only one time in each dimension. Text files containing variant positions were created, and variants unique to a single file were extracted using awk ′END {for (R in r) {split(r[R], t, SUBSEP) if (!t[1]) print t[3], t[2]}}{k = $1 SUBSEP $2 SUBSEP $3r[k] = c[k]++ SUBSEP FILENAME SUBSEP $0 }′ *.txt > output. This resulted in three files containing unique variant calls for each pooling dimension. True mutations should be present one time in each dimension and so should be represented in each of the three lists. Matching values fitting these criteria were extracted using grep -wFf. Original VCF files were evaluated in cases where positive control mutations were not identified using this method. This was how the mutation in VDE_2_27 in pool 14 was recovered.
Data availability
The raw datasets analyzed during the current study are available in the Sequence Read Archive, Accession PRJNA422048 at https://www.ncbi.nlm.nih.gov/bioproject/422048. M3 seed from positive control mutants is available upon request while supplies last. Table S1 contains target names, primer sequences, amplicon lengths and percent GC of each amplicon. Table S2 contains pre-set parameters of sonication of genomic DNA. Table S3 contains additional tested sonication parameters. Table S4 contains sequencing statistics for mutant pool 10. Table S5 contains sequencing statistics for mutant pool 11. Table S6 contains sequencing statistics for mutant pool 15. Table S7 contains sequencing statistics for mutant pool 19. Table S8 contains sequencing statistics for mutant pool 20. Table S9 contains sequencing statistics for mutant pool 37. Table S10 contains sequencing statistics for mutant pool 41. Table S11 contains summary sequencing statistics for 48 barley pools. Table S12 contains coverage data and calculations for 32 barley amplicons in 48 pools. Table S13 contains coverage data from a tomato TILLING by Sequencing experiment using non-fragmented amplicons. Table S14 contains coverage data from cassava sequencing using non-fragmented amplicons. Table S15 contains frequencies of mutations called with the CAMBa and GATK HaplotypeCaller tools. Table S16 contains cost estimations for experiments using fragmented and non-fragmented pooled amplicons. Figure S1 shows Fragment Analyzer profiles of a single amplicon subjected to different default sonication settings. Figure S2 shows Fragment Analyzer profiles of a single amplicon subjected to different sonication parameters. Figure S3 shows Fragment Analyzer profiles of a single amplicon subjected to additional sonication parameters. Figure S4 shows Fragment Analyzer profiles of pooled amplicons before and after sonication. Figure S5 shows read coverage, base coverage and mapping quality for tomato amplicons not subjected to sonication. Figure S6 shows insert sizes for 48 barley sequencing pools. Figure S7 shows principle component analysis of sequencing statistics for 48 barley pools. Figure S8 contains a graph of mean coverage and amplicon size in 48 barley pools. Figure S9 contains a graph of mean coverage and percentage of GC in amplicons in 48 barley pools. Figure S10 contains a graph of percentage of GC and amplicon size in base pairs. Supplemental material available at FigShare: https://doi.org/10.25387/g3.8052821.
Results
Evaluation of sonication parameters established for genomic DNA on PCR products
A PCR product of 964 base pairs (bp) amplified from Coffea arabica genomic DNA was created for initial experiments to optimize DNA fragmentation. The aim of the work was to identify conditions resulting in a fragmented PCR product between 350 and 500 bp using a Covaris focused-ultrasonicator. The manufacturer’s parameters established for genomic DNA shearing were initially tested (Table S2). Fragment analysis showed limited fragmentation and retention of a high concentration of 964 bp product at all pre-defined parameters (Figure S1).
Optimization of parameters for sonication of a single PCR amplicon
Based on the results from the default parameters for genomic DNA we hypothesized that halving the power, duty factor and/or cycles/burst may improve fragmentation of lower molecular weight PCR products. Fourteen distinct sonication parameter combinations were chosen and evaluated (Table 1). Test numbers 1 through 3 were designed to assess the effect of the number of cycles/burst on the fragment size and distribution. Tests 4 through 6 were designed to assess the effect of the duty factor on the fragment size and distribution. Tests 7 through 9 were designed to assess the effect of the peak power on the fragment size and distribution. Tests 10 through 12 were designed to assess the effect of the duration of the sonication on the fragment size and distribution, and tests 13 and 14 were designed to test the combination of high duty factor with an average number of cycles/burst and very low duty factor with the highest cycles/burst respectively.
The extent of PCR fragmentation was evaluated using a capillary Fragment Analyzer (FA), (Figure S2). Relative concentration of DNA fragments was evaluated owing to the fact that the absolute concentration varies according to amount of input PCR product used in the assay. Based on this analysis, parameters of test # 13 were chosen to provide the best fragmentation as they produced a broad distribution of fragments with a peak of 390 bp, which was close to 350 bp, one of two sizes recommended for Illumina library preparation (Figure 1). Additional modifications, whereby average power (Peak Power and Duty Factor) were held constant while time, cycles/ burst and volume were modified, showed no substantial improvement (Table S3, Figure S3).
Sonication of complex pools of PCR amplicons
Having established optimal conditions for fragmentation of a single amplicon, the parameters were next tested in a TILLING by Sequencing experiment containing 32 amplicons ranging between 670-1513 base pairs. The GC base percentage in the amplicons ranged from 35.16 to 68.66 with a median of 50.28 (Table S1). For comparison, GC values of the entire coding sequence assembly of the barley genome build 2 was calculated as having a minimum GC percentage of 22.86, maximum of 86.96 and median of 52.48 (see methods). The PCR primers producing amplicons tested in this study had been previously validated and used in traditional TILLING assays employing gel-based cleavage assays using an eightfold genomic DNA pooling strategy (Stolarek et al. 2015a, 2015b; Szarejko et al. 2017). Thus, the amplicon set represents the complexity of a TILLING project, and also allows testing of the feasibility of applying sonication-based fragmentation and next generation sequencing to an already existing TILLING platform. Fragmentation was performed on 48 pools followed by preparation of Illumina libraries (Figure S4). Fragmentation profiles of the prepared libraries showed more complex patterns with average sizes of approximately 650 bp, that includes the size of the target DNA (known as the insert) plus the ligated adapter sequences. The presence of higher concentration peaks at specific molecular weights is also observed (Figure 2).
Fragmentation profiles of Illumina sequencing libraries prepared from 32 pooled amplicons produced by PCR amplification of genomic DNA from 256 pooled barley samples. PCR products from forty-eight genomic DNA pools were subjected to sonication using test parameter 13 in Table 1. Fragmentation profiles of two different libraries are shown (top and bottom panel).
Sequencing libraries prepared from pooled fragmented amplicons
To evaluate the effect of fragmentation on sequencing coverage and quality, the Illumina libraries were subjected to 2x300 PE sequencing on a MiSeq instrument. Sequencing coverage profiles, per base mapping quality, and per base sequence quality of four unique amplicons are shown (Figure 3, Tables S4, S5, S6, S7, S8, S9, S10). Similar data were graphed for two pools and two amplicons from a previously published project where TILLING screening was performed by direct sequencing (without fragmentation) of short amplicons derived from PCR from pools of 64 tomato genomic DNAs (Gupta et al. 2017) (Figure S5). The minimum insert size in all 48 barley pools was 259 bp, maximum was 386 bp, and median was 348.23 bp (Table S11). Pools show a distribution of sizes, with accumulation of higher concentration of DNA at specific molecular weights. This is pronounced in the most highly concentrated pool number 24 (Figure S6). At the pool level, mean coverage ranged from 1492 to 50030 with a median of 6365. PCA analysis showed pool 24 to be the outlier with coverage of 50030 (Figure S7). Percent GC ranged from 47.64 to 54.33 and mapping quality from 59.06 to 59.76 (Table S11).
Read coverage, base quality and mapping quality from 2x300PE sequencing of targets EXPB4 (1350 bp) (A), RTH3 (1028 bp) (B), PRT1 (979 bp) (C) and ALS3_1 (670 bp) (D). The number of aligned reads (coverage) at each position in the amplicon is marked in blue. The average base quality per position is marked in yellow and mapping quality in red. Coverage is even across the majority of each amplicon with slight increases in coverage and concomitant decreases in quality score toward the end of each amplicon. Shorter amplicons have longer regions of increased coverage at the ends.
To compare sonicated PCR products to direct sequencing of PCR products, average coverage values for each amplicon in each of the 48 barley pools was tabulated (Table S12). Data were also tabulated from previously published tomato and cassava experiments that employed direct sequencing of PCR products without fragmentation (Table S13, S14) (Duitama et al. 2017; Gupta et al. 2017). In barley, 1536 coverage values (32 amplicons × 48 pools) were evaluated. The minimum mean coverage was observed in the 1153 bp amplicon named FC1_11 (value = 2093). Maximum coverage was in the 778 bp BAK1_1 (134456), with a mean coverage of 6621 for the total collection of amplicons. A weak negative correlation was observed between amplicon size and mean coverage (r = -0.38), and between percentage of GC in the amplicon and mean coverage (r = -0.41) (Figures S8, S9 and Table S12). A weak positive correlation was observed between amplicon size and percentage of GC (r = 0.25) (Figure S10 and Table S12). To estimate failures, data points were counted where coverage fell below 5% of the mean coverage for the tested amplicon. This detects genomic DNA pools where a single amplicon produces much less coverage than the same amplicon in other DNA pools. Using this criterion, 2.08% of the experiment failed. This is in comparison to 0.5% for tomato (n = 2408) and 1.19% for cassava (n = 2538).
To evaluate the suitability of the experiment to recover induced mutations, the number of times a mutant allele would be sequenced was calculated. In the pooling scheme tested, a heterozygous mutation should be present at a frequency of 0.00195 (1/(256*2)). A coverage of 512 is therefore needed to identify the heterozygous mutation 1 time in a pool. Using this, the percentage of the assay where heterozygous mutations are expected at a specific coverage was calculated. One hundred percent of the experiment produced data for 1x coverage. This reduced to 96.74% for 2x coverage, 88.15% for 4x and 62.30% for 8x coverage (Table S12). Homozygous mutations are twice the concentration in an individual plant and therefore 88.15% of the experiment produced data for recovery of homozygous mutations at 8x coverage.
Mutation discovery in complex pools
To evaluate the use of sonication for TILLING assays involving complex pooled PCR products, the mapped reads from sequencing were subjected to two mutation discovery pipelines. Three plants harboring previously identified induced mutations were used to evaluate mutation discovery (Szarejko et al. 2017). The CAMBa pipeline discovered all positive controls (Table S15). Control 1, a plant harboring a C to T homozygous mutation at position 416 in the PRT1_2 target was identified at allele frequencies of 0.0078 in pool 11 (first dimension of pooling), 0.0103 in pool 19 (second dimension of pooling) and 0.0117 in pool 41 (third dimension of pooling). The expected frequency of a homozygous mutation is 0.0039. Control 2, a plant harboring a C to T heterozygous mutation in target EXPB4_23 was identified in pool 15 (first dimension) at a frequency of 0.00778, in pool 20 (second dimension) at 0.0024 and in pool 37 (third dimension) at 0.00199. The expected frequency of a heterozygous mutation is 0.00195. The heterozygous T to A mutation in control 3, VDE_2_27, was identified in pool 14 (0.0026), 26 (0.0051) and 46 (0.0060). GATK HaplotypeCaller was performed using two different ploidy levels to recover rare alleles. Variant calling at ploidy settings 600 and 2000 both recovered the homozygous mutation in control 1. Allele frequencies in pools 11, 19, and 41 for ploidy 600 were 0.0033, 0.0050 and 0.0067, respectively. For ploidy 2000 the frequencies were 0.0035, 0.0055 and 0.0070. Control 2 was not identified at either ploidy setting. The control 3 mutation was identified only in pool 14 at frequencies of 0.0033 and 0.0030 for ploidy settings 600 and 2000, respectively. No trend was observed between sequencing coverage of the amplicon in the pool and observed vs. expected allele frequencies (Table S15).
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 en