Whole genome profiling physical map and ancestral annotation of tobacco Hicks Broadleaf

WGP tag generation

A total of 1107 384-well plates from four libraries (425 088 BAC clones) was subjected to WGP™ (Keygene NV, Wageningen, The Netherlands, http://keygene.com), as described by van Oeveren et al. (2011). BAC pooling was performed in a two-dimensional format, with each pool consisting of 48 clones. DNA isolated from these pools was digested using EcoRI and MseI, bar-coded adapters were ligated, and restriction fragments were amplified. Amplicons were subsequently sequenced from the EcoRI restriction site end using the Illumina Genome Analyzer II. In total, this yielded 1718 million reads with a 78 nt read length, a valid sample identification tag and proper restriction site sequence. Deconvolution was performed using the first 31 nt, the first 51 nt and the full 70 nt reads (excluding the bar-code sequences), resulting in half of the total number of reads being assigned as WGP tags to individual BACs. After filtering these WGP tags using several quality criteria, a total of 1.2 million WGP tags were mapped to 361 034 BACs for the 51 nt set (Table 1). The results for all three sets are presented in Table S1. The 51 nt WGP tag analysis was selected as the final dataset because the analysis demonstrated that this read length provided the maximum number of deconvoluted reads.

Physical map construction

After conversion of the WGP tag data to pseudo-mobility numbers, the FPC software (Soderlund et al., 1997) adapted for use with sequence tags as described by van Oeveren et al. (2011) was used to assemble sequence-based physical BAC maps. A cut-off value of 10⁻²⁵ was used with a subsequent DQ step to produce a high-stringency map. Additional end-to-end contig and singleton merging steps at 10⁻¹⁵ were performed to obtain a normal-stringency map with fewer contigs. Using the 51 nt set and normal-stringency settings, the tobacco physical map comprised 9750 contigs containing 330 632 BACs (77.8% of the BACs tested; 91.5% of the tagged BACs). Of the remaining BACs, 30 402 did not link to any other BAC and were classified as singletons, and 64 064 BACs had no deconvoluted tags. The estimated mean and N₅₀ contig sizes were 462 and 689 kbp, respectively, and the calculated genome coverage was 4508 Mbp (Table 1). The results for both stringency settings for all three sets are presented in Table S2.

Ancestral origin of tobacco BACs and WGP contigs

Sequence tags flanking EcoRI restriction sites for the whole genomes of N. sylvestris and N. tomentosiformis, the two closest relatives to the ancestral contributors of the S and T genomes to tobacco, were obtained by digesting genomic DNA of these lines using EcoRI and MseI, followed by amplification and Genome Analyzer II sequencing using the same method as applied for WGP map construction of N. tabacum. After quality filtering, 1 089 317 N. sylvestris and 1 035 343 N. tomentosiformis tags were obtained, compared with the 1 239 733 unique 51 nt tags obtained for N. tabacum during construction of the physical map, thus allowing determination of the ancestral origin of the tobacco tags, BACs and WGP physical map contigs by P value calculations for S or T tag enrichment. This analysis showed that 60.7% of the BACs are of S origin and 37.4% of T origin. The origin of 0.4% of the BACs could not be determined despite their S and/or T tags, and these BACs were classified as having an undefined origin. Finally, 1.5% of the BACs have no S or T tags and are thus of unknown origin. BACs consist of tobacco DNA fragments of approximately 100 kb, and as such should originate from only one of the two ancestors. Being able to assign an ancestral origin to 98.1% of the BACs indicates that S and T tags are not present randomly in BACs, thus validating the origin determination of WGP tags. The P value calculations for S or T tag enrichment further showed that 53.7% of the WGP contigs are of S origin, 45.8% are of T origin, 0.5% are of undefined origin, and only four are of unknown origin (Figure 1 and 2). Thus the S and T tags are rarely found together in the same WGP contig, providing additional confidence in the quality of the constructed physical map.

Figure 1

Log₁₀ frequency of the weighted proportion of S tags. A value of 0 indicates a T origin (red bars) and a value of one indicates an S origin (blue bars), according to the enrichment P values (see Experimental procedures). The gray bars indicate counts of: (a) BACs. (b) WGP contigs of undefined origin, i.e. for which unequivocal assignment to S or T was not possible. Absolute counts are given above each bar.

Comparing the ancestral origin of BACs to the ancestral origin of the WGP contigs in which they were placed showed that the majority of the BACs of S or T origin were found in WGP contigs with the same origin annotation (98.0 and 94.0% for BACs of S or T origin, respectively) (Table 3). BACs of undefined origin were almost evenly placed in WGP contigs of S or T origin (43.8 and 53.1%, respectively), whereas 78.1% of the BACs of unknown origin were placed in WGP contigs of T origin and 19.4% in WGP contigs of S origin.

Taking advantage of the ordering of the BACs in a WGP contig (where no ordering is available for the tags in a BAC) and to complement the enrichment P values, the number of domains composed of BACs of S or T origins was estimated for each WGP contig, as described in Experimental procedures. WGP contigs with a clear origin are ideally composed of only one domain, whereas WGP contigs of undefined origin either possess more than one domain or are artifacts created during the WGP physical map construction. The number of domains counted in WGP contigs of S, T or undefined origin showed that only one domain was identified for 94.1 and 93.2% of the WGP contigs of S or T origins, respectively (Table S3). Enrichment P values, which are only based on counts and do not consider contig units, are therefore well complemented by the domain structure, as some S and T contigs contain several domains, with the majority corresponding to the assigned origin.

Linking the genetic and physical maps

The N. tabacum genetic map was constructed from a Hicks Broadleaf × Red Russian F₂ mapping population (Bindler et al., 2011). It consisted of 24 linkage groups comprising 1776 unique loci determined by 2318 SSR markers. SSR amplification tests on N. sylvestris and N. tomentosiformis were then used to annotate tobacco linkage group regions according to their putative S or T ancestral origin. Of the 2318 SSR markers of the N. tabacum genetic map, 918 were experimentally linked to 910 BACs. In the WGP physical map, 802 of these BACs were contained in 725 WGP contigs.

The number of BACs and WGP contigs of S, T or undefined origin mapped to these S or T linkage group regions is shown in Table 4. Almost 80 and 65% of the WGP contigs of S or T origins were indicated as linked to S or T linkage group regions, respectively 80% S or 65% T.

For each linkage group, Figure S1 shows the location of SSR markers of S or T origin used for the genetic map construction and that of the BACs of S or T origins that are linked to the tobacco genetic map. The linkage groups are colored according to their S or T annotation (based on Bindler et al., 2011). For linkage group 22, a color inversion in the S and T assignment is visible in Figure 1 of Bindler et al. (2011). Recalculation of ancestral origins after correcting this possible inversion increased the fraction of BACs for which the predicted origin corresponds with the genetic map and decreased the number of conflicting cases by approximately 2% (Table S4).

BAC libraries

Four BAC libraries of the allotetraploid genome of N. tabacum cv. Hicks Broadleaf (inbred line PI 552397 or TC 311 in the United States Department of Agriculture, Agricultural Research Service, Germplasm Resources Information Network (USDA-ARS-GRIN database, http://www.ars-grin.gov/), comprising a total of 425 088 BAC clones (1107 384-well plates) and approximately 10.4-fold genome coverage, were used for construction of the WGP map. These libraries were (i) a HindIII library (Microsynth AG, Balgach, Switzerland, http://www.microsynth.ch) consisting of 112 896 clones with an estimated mean insert size of 100 kb, representing approximately 2.5-fold genome coverage, (ii) a BamHI library (Microsynth), consisting of 146 304 clones with an estimated mean insert size of 100 kb, representing approximately 3.4-fold genome coverage, (iii) an EcoRI library consisting of 69 120 clones, with an estimated mean insert size of 125 kb, representing approximately 1.9-fold genome coverage, and (iv) a HindIII library, consisting of 96 768 clones with an estimated mean insert size of 125 kb, representing approximately 2.7-fold genome coverage.

BAC library pooling, DNA isolation and WGP sample preparation

Individual BACs of the four libraries stored in 384-well plates were pooled in a 2D format prior to DNA isolation. Specifically, BACs were 2D-pooled by taking six plates at a time in a 2 × 3 layout (termed a superpool) and pooling each row over two plates (48 BACs) and each column over three plates (48 BACs), yielding 96 pools per superpool. BAC pools were subjected to isolation of high-concentration DNA, and WGP sample preparation was performed essentially as described by van Oeveren et al. (2011). Briefly, AFLP® templates (Vos et al., 1995) were prepared from pooled BAC DNA by digestion using five units of EcoRI and two units of MseI. Next, adapter ligation was performed using either a P5 EcoRI adapter containing a 5 or 6 nt sample identification tag in combination with a unique P7 MseI adaptor (when one or two superpools were pooled per lane for sequencing), or using a P5 EcoRI adaptor containing a 5 nt sample identification tag and P7 MseI adaptor with a 3 nt sample identification tag (when three superpools were pooled for sequencing in one lane). PCR was performed in 20 μl volumes, and contained 5 μl of 10-fold diluted restriction ligation mixture, 30 ng Illumina P5 primer (5′-AATGATACGGCGACCACCG-3′), 30 ng Illumina P7 primer (5′-CAAGCAGAAGACGGCATACGA-3′), 0.2 mm dNTPs, 0.4 units of Amplitaq (Applied Biosystems, http://www.appliedbiosystems.com) and 1× Amplitaq buffer. Next, equal amounts of BAC pool PCR reaction mixtures were pooled per superpool, and purified using a QIAquick PCR purification kit (Qiagen, http://www.qiagen.com). When one or two superpools were sequenced per lane, each BAC pool was bar-coded using one of 192 different 5 or 6 nt sample identification tags in the P5 EcoRI adapter. When three different superpools were combined for sequencing in one lane, each BAC pool was tagged by a combination of one of 96 different P5 EcoRI adapters with a 5 nt sample identification tag, and one of three different P7 MseI adapters with a 3 nt sample identification tag, such that 288 BAC pools could be pooled per lane for sequencing.

Sequencing

A total of 90 lanes of either 76 or 78 cycles of sequencing, divided over 15 runs, were performed using an Illumina (San Diego, CA, USA, http://www.illumina.com) Genome Analyzer II. Sequencing with 76 cycles was used when one or two superpools were pooled per lane, whereas 78-cycle sequencing with a second priming event at the MseI side was used when three superpools were pooled per lane, in which case the 3 nt sample identification tags at the MseI side were determined. Each run used a flow cell with eight lanes of physically separated samples, with the same set of sample tags used for each lane. Genome Analyzer II runs were performed with seven lanes of tobacco WGP samples at 5.5 pm concentration, each covering one, two or three superpools represented in 96, 192 or 288 row and column pools, respectively. The Illumina pipeline software was used to extract sequence reads of 76 or 78 nt length from the images. An additional quality filter was applied to select only those reads for which all base calls have at least a Solexa quality of 0 (equivalent to a Phred quality of 3, http://www.phrap.com/phred) on the Illumina GA scale. All sequence data have been deposited in the European Nucleotide Archive (ENA) under accession number ERP001765.

Deconvolution and filtering

Sequence reads were split into three parts in the case of 76-cycle sequencing or in four parts in the case of 78-cycle sequencing to enable assignment of unique tags to pools and to allow consecutive deconvolution into individual BACs. In the case of 76-cycle sequencing, the first 5 or 6 nt represented the sample (i.e. BAC pool) identification tag at the EcoRI side, the next 6 nt matched the EcoRI restriction site of the adapter, and the remaining nucleotides defined the ‘WGP tag’. In the case of 78-cycle sequencing with combinatorial bar-coding at the EcoRI and MseI side ends, a fourth section representing the 3 nt sample identification tag at the MseI side was obtained.

Three different WGP tag lengths were analyzed to test the effect of read length on WGP map resolution: 31, 51 and 70 nt (all including the EcoRI restriction site). The assignment of unique WGP tags to individual BACs was based on the following criteria. First, a specific WGP tag must occur in two pools to indicate its unique position on the plate – one column and one row pool with both being represented by at least three reads. Second, if WGP tags are inadvertently observed in a third or fourth pool, the number of reads in these other pools must be less than a tenth of those in the smallest true pool. WGP tags not matching these criteria were discarded. A computer script was used to recognize and trim the sample identification tags and the restriction site part of the sequence reads and to perform the deconvolution. Unique WGP tags were defined by grouping them in 100% identical read sets. The output of this procedure consisted of a list of all WGP tags, the corresponding number of reads, and the identification number of the BAC they were assigned to. Finally, a filtering step was applied to remove WGP tags matching vector, Escherichia coli or chloroplast sequences, those containing homopolymer sequences of 5 nt or longer, and WGP tags occurring on just a single BAC.

Contig building

Contig building was performed using the FPC program (Soderlund et al., 1997). This software tool was originally developed for analyzing BAC fingerprint data, i.e. restriction fragments determined by their length. The WGP tags were adapted for use in FPC as described by van Oeveren et al. (2011) by converting them into numbers to yield pseudo restriction fragment sizes for which the software was originally designed. As the WGP tags are uniquely defined by their sequence composition, FPC was used at the highest stringency setting of tolerance (value = 0). Various cut-off values were tested, specifying the threshold on the probability of BAC coincidence, i.e. the likelihood that different BACs containing partly overlapping sets of WGP tags originate from the same genomic region. The output of FPC consisted of a list with contigs and the corresponding order of BACs within each contig. The genome coverage, mean contig size and N₅₀ contig size in million base pairs were estimated by multiplying FPC band units by the mean distance between two WGP tags. The latter was estimated by dividing the mean BAC insert size by the mean number of WGP tags.

Linking of BACs to the genetic map

To identify BACs containing SSR markers, 1185 SSR primer pairs were screened from pools and superpools of 196 BamHI and 193 HindIII BAC library plates representing approximately threefold coverage of the tobacco genome (approximately 149 000 BACs).

Determination of ancestral origin

Nuclear genomic DNA was isolated from seedling leaf material of N. sylvestris TW138 and N. tomentosiformis TW142 in accordance with a protocol for high-molecular-weight nuclear DNA extraction (Liu and Whittier, 1994). As the N. tomentosiformis line did not yield enough nuclear DNA, genomic DNA isolation was performed using the Nucleon plant genomic DNA extraction kit (Hologic Gen-Probe Inc., Bedford, MA, USA, http://www.gen-probe.com) according to the manufacturer's instructions to obtain sufficient material. DNA samples were used for EcoRI/MseI (E/M) AFLP template preparation as described by Vos et al. (1995), using two different bar-coded adaptors suitable for Genome Analyzer II sequencing, in a similar approach to the WGP sample preparation. Seven lanes of a single Genome Analyzer II run were used with a read length of 75 nt, yielding 63.9 million reads for these two samples. Reads for which one or more sequenced nucleotide had a phred quality score lower than 20 (i.e. a probability of error of 0.01) were discarded. For both ancestral species, unique sequences of 51 nt length were grouped into tags, and tags found in both species or with only one read were removed.

The S or T ancestral origin of each BAC included in the physical map construction and each WGP contig was determined by calculating the S or T enrichment P values. The proportion of S tags and T tags, weighted according to the number of reads in N. sylvestris and N. tomentosiformis, respectively, was computed for each BAC or WGP contig. If the BAC or WGP contig was enriched with S or T tags), the estimated proportion was far from a random selection of tags from the original pool. Therefore, an enrichment P value for the S and T tag proportions was estimated using the hypergeometric distribution for the null hypothesis. The P values obtained were subsequently corrected for multiple testing effects using the Benjamini–Hochberg false discovery rate (Benjamini and Hochberg, 1995). The ancestral origin of a BAC or WGP contig was predicted to be S or T if the enrichment P value for S or T was <10⁻⁶, and the enrichment P value for T or S, respectively, was above 10⁻⁶, otherwise the origin was undefined.

To elucidate the ancestral structure of a WGP contig, each WGP contig unit, as defined by the FPC software, was assigned a putative label S or T if the number of S or T BACs covering the unit exceeds the 95% quantile of a binomial distribution with parameter 0.5. Otherwise the unit was labeled ‘undefined’. The advantage of using this quantile over a constant threshold is that it penalizes low-coverage regions of the WGP contig.

Once putative labels were assigned, the WGP contig was split into connected S, T or undefined segments by mimicking a classification and regression tree approach (e.g. Hastie et al., 2001). To control the level of smoothing, a minimum segment size of 40 was chosen. The binary rules obtained were then parsed to extract the structure of the WGP contig. Each segment of a given origin was called a domain of the WGP contig.

Scaffolding

The position on a WGP contig covered by all the BACs containing a given WGP tag, and by no other BACs, was calculated for all WGP tags belonging to only one WGP contig. These WGP tags were mapped to the Sol Genomics Network version 1 assembly of the Tobacco Genome Initiative reads, and perfect matches were retained. Assembly contigs with at least two mapped WGP tags from two parts of the same WGP contig were selected, as they could be mapped with orientation to the WGP contig. Non-overlapping assembly contigs mapping to the same WGP contig were scaffolded. Assembly contigs with mapped WGP tags from two different WGP contigs were used to scaffold these WGP contigs. The positions of the WGP tags and of the BACs in the 51 nt normal-stringency WGP physical map are shown in Tables S6 and S7, respectively.