Development of genomic resources for the narrow-leafed lupin (Lupinus angustifolius): construction of a bacterial artificial chromosome (BAC) library and BAC-end sequencing

Construction of a BAC library from NLL cv Tanjil

A NLL cv. Tanjil BAC library was constructed from hydroponically grown seedlings using the Bam H1 restriction enzyme. The library comprised of 111,360 BAC clones which were stored in 292 384-well microtitre plates. The average insert size of the library was approximately 99.7 kb based on the analysis of ca. 250 randomly selected clones. Approximately 2% of the clones contained no insert while the majority (ca. 75%) had insert sizes of between 90 kb and 110 kb, 12.5% above 110 kb, 9.2% between 80 kb and 90 kb, 3% between 50 kb and 80 kb and 0.3% smaller than 50 kb. The coverage of the library was estimated to be around 12 haploid genome equivalents according to the (haploid) genome size estimate of approximately 924 Mbp for NLL [36, 37].

Screening the BAC library for β-conglutin containing clones

To demonstrate the utility of the library, 111,360 BAC clones which were double-spotted onto six nylon membrane filters using a robot were screened for NLL β-conglutin genes using a cDNA fragment from a β-conglutin gene (BETA2) as the probe [38]. This cDNA probe cross-hybridised to the other six members of the β-conglutin gene family due to the high level of homology among the gene sequences [38]. The screening permitted 108 positive clones to be identified. Twelve clones were randomly selected and confirmed to contain β-conglutin gene(s) by PCR with a pair of primers able to recognise all seven β-conglutin gene members.

Characterisation of β-conglutin containing BAC clones

For further analysis, eight of the BAC clones were analysed for the presence of specific β-conglutin genes by PCR using primers and annealing temperatures specific for each β-conglutin gene, followed by sequencing of the PCR product. As shown in Figure 1, each of the eight β-conglutin hybridising BAC clones contained β-conglutin genes, including BETA1, BETA2, BETA4, BETA5, and a new β-conglutin gene, recognised by BETA2 primers, which has been renamed, BETA8. The truncated sequences of the eight β-conglutin genes are shown in Additional file 1. It is not clear why the BETA5 primers recognised BETA1 sequences from BACB but not from BACF. BETA5 primers detected two bands from BACH, however sequence data showed that they both contained BETA5 sequences, indicating a PCR aberration. BETA1 primers identified a product with sequence identity to BETA1 except there was an extra 111 bp, presumably an intron within BETA1. Apart from BACF, which contained BETA1 and BETA8, each BAC clone only contained one β-conglutin gene.

Predicting characteristics of the NLL genome from bioinformatic analysis of a randomly sampled collection of BAC-end sequences

To provide insight into the sequence content and complexity of the NLL genome, we sequenced and analyzed both ends of 9600 BAC clones randomly selected from the library and generated a total of 17856 BAC-end sequences (BESs). Among these sequences, 13895 BESs (77.8%) of BAC-end sequences contained usable sequence data (> 100 bp after trimming) with a total length of ca. 8.89 Mbp which covers 0.96% of the NLL genome. The average length of the BESs was 683 bp with a maximum of 1112 bp and minimum of 101 bp. A total of 3961 of the BESs were discarded, because sequences were too short (< 100 bp), matched to the vector, showed low complexity or were identified as bacterial or organelle sequences.

To estimate the protein-encoding gene content of the NLL genome, the 13895 BESs were compared to the NCBI NR protein database by BLASTx. Protein-matching regions were found in 2667 BESs, 1723 of which also contained matches to known repetitive DNA sequences. Putative gene-encoding regions totalled 483,216 bp, equivalent to 5.4% of the total BES dataset. Based on this value and an estimated average gene length of 956 bp in legume species and 1170 bp over all plant species (source: http://www.phytozome.net) the NLL genome was predicted to contain between 42656 and 52204 genes. Blast2GO assigned 12831 Gene Ontology (GO) terms to 2930 BESs (Figure 2, Additional files 2 and 3). GO annotations were converted to 5448 FunCAT annotations assigned to 2067 BESs (Additional file 4).

Functional analysis of the encoded genes in the 13985 BESs predicted high proportions of the NLL gene content involved in the cell cycle and DNA processing (16.42%), protein binding (15.02%), protein folding and modification (5.75%), organelle proteins (4.76%) and metabolism (3.82%) (Additional files 5 and 6). A summary of gene ontology revealed a relatively higher proportion within the NLL BESs of genes associated with DNA metabolism, proteolysis, auxin biosynthesis, cellular processes, primary metabolism, oxidation and reduction and serine family amino acid metabolism (Figure 2). Conversely, there was a relative scarcity in the BESs of genes involved in organ formation, signal transduction, transmembrane transport, gluconeogenesis, alkaloid biosynthesis, starch and sucrose metabolism and defense/stress responses (Figure 2).

Comparative genomics between NLL and other plant species

To examine the phylogenetic relationship between NLL and sequenced species, the 13895 quality-screened BESs were aligned (via BLASTn) to the NCBI Nucleotide database and the probable phylogenetic distribution of these hits was visualised with MEGAN [39, 40] (Figure 3). MEGAN indicated the level of species similarity and sequence conservation within the randomly sampled BES subset of the NLL genome. About 18.3% of BESs matched sequences available in the NCBI Nucleotide database, while 79.2% of BESs had no hits. This is due in part to the relatively low number of sequences from Lupinus and other closely related genera that are currently available. The majority (84%) of the matched BESs were mapped to eudicotyledon species, 77.9% of which were assigned to the species of the subfamily Papilionoideae (Leguminosae). The three sequenced leguminous species, G. max, L. japonicus and M. truncatula, were highly represented followed by the genera Lupinus and Arachis. In contrast, there were only six NLL BESs that were best aligned to Arabidopsis sequences and no BESs specifically aligned to Oryza sativa. The average percent identity of BLASTn alignments between NLL BESs and G. max, M. truncatula and L. japonicus was 95.8%, 94.4% and 90.6% respectively. However the average identity restricted to predicted gene-encoded regions was 87.0%, 88.5% and 86.7% respectively. A small proportion of BESs (ca. 0.7%) were aligned to proteobacteria, perhaps due to the fact that the tissues for constructing the BAC library were collected from lupin seedlings grown under non-sterile conditions.

Simple sequence repeat (SSR) profiling and SSR markers

However, analysis of the relative frequency of individual groups of SSRs with motif length 1-8 bp revealed some major trends for all six species as well as some distinctions for NLL SSRs (Figure 4). While the di-nucleotide SSRs in NLL were the most abundant of all the class I SSRs (≥ 20 nt), in line with the other 5 species, the mono- and octa-nucleotide SSRs showed higher representation and the tri-, tetra- and penta-nucleotide SSRs lower representation in NLL than those in other species. In contrast, the distribution of the relative abundance of each repeat motif length of the class II SSRs (12-19 nt) in NLL appeared consistent with that in the other species. Penta-nucleotide SSRs were an exception, and were relatively more abundant by about four times in NLL, than those of the other species. In the combination of both class I and class II SSRs, NLL appeared to have higher abundance of mono-nucleotide SSRs and lower abundance of tri- and tetra-nucleotide SSRs. In NLL, most of the penta-nucleotide SSRs had short repeat sequences and grouped into class II (12-19 bp).

To develop SSR markers for NLL, primer-pairs were designed flanking 2023 SSRs using Primer3 [42] with additional criteria described in the Materials and Methods. After taking repetitive sequences into account, there were 1497 non-redundant SSR marker candidates including 24 class I SSRs (≥ 20 bp) and 1455 class II SSRs (12-19 bp). These SSRs and the details of their primers are shown in Additional file 8. Twenty four additional class I SSR markers were designed using less stringent criteria (see Methods). The details of these 24 SSRs are shown in Additional file 9.

As expected the "wild" NLL line, P27255, was the most divergent, with the three domesticated lines being closely conserved (3). Nine out of 24 markers were polymorphic between the parents of the narrow cross (cv Tanjil × cv Unicrop), whereas 14 markers were polymorphic between the parents of the wide cross (P27255 "wild" × 83A:476 "domestic").

Plant material

For the construction of the BAC library, seeds of Lupinus angustifolius L., narrow-leafed lupin (NLL), cv. Tanjil were germinated on moisturised filter paper at room temperature for two days. Seeds were subsequently grown hydroponically in half-strength Hoagland solution in a growth room at 22°C over a 16 h/8 h day/night schedule. After ten days, leaves were collected, frozen in liquid nitrogen and stored at -80°C.

For the Multiplex-Ready PCR assays to determine SSR lengths in four NLL lines, seeds of cv. Tanjil, cv. Unicrop, 83A:476 ("domestic") and P27255 ("wild") were germinated on moisturised filter paper at room temperature for two days. Subsequently, the germinated seeds were transferred into pots and grown in a temperature controlled growth chamber over 16 h/8 h day/night schedule, using fluorescent light at 100 to 120 μE m^-2 s^-1 and a constant temperature of 22°C. Leaf material was harvested from two-week-old plants for DNA isolation.

Construction of the BAC library

BAC library construction was performed at the Australian Centre for Plant Functional Genomics, University of Adelaide, South Australia. The procedure has been described in detail in Shi et al. [54] and high-molecular weight DNA preparations were generated according to Zhang et al. [55]. Leaf tissue was ground in liquid nitrogen and the powder combined with homogenisation buffer. The homogenised sample was mixed with an equal volume of pre-warmed 1% low melting point (LMP) agarose and cast into plugs using plug molds (Bio-Rad). Following lysis, the agarose plugs were cut and digested with the restriction enzyme Bam HI. The DNA slices were size fractionated using a CHEF Mapper XA pulse-Field gel electrophoresis (PFGE) system (Bio-Rad), for 18 h at 11°C and 6 V/cm, using a 1-40 s pulse time and 120° field angle.

The 100-250 kb DNA fractions were excised and the DNA eluted in a Bio-Rad electro-eluter (Model 422), applying 10 mA per tube. Size-fractionated DNA was ligated to Bam HI Cloning-Ready pIndigoBAC-5 vector DNA (Epicenter) and transformed into ElectroMAX E. coli DH10B competent cells (Invitrogen). The transformation was carried out using a Bio-Rad Gene Pulser Xcell and a cuvette with a 1 mm gap at 1800 volts (Bio-Rad). Transformed cells were incubated at 37°C for 1 h in 1 mL LB, plated on to LB agar medium containing 12.5 μg/mL chloramphenicol and grown overnight at 37°C. Colonies were picked into wells of 384-well plates containing 70 μL LB freezing medium using a VersArray Colony Picker and Array System robot (Bio-Rad). Plates were incubated overnight at 37°C and used to make three copies of the library. Libraries were stored at -80°C. To print filters for hybridisation screening, clones from one copy of the library were arrayed in duplicate onto 22 cm × 22 cm positively charged nylon membranes (Amersham Hybond-N+ from GE Healthcare, or Performa II from Genetix) using a Qpix2 robot (Genetix). The three copies of the library were later transferred to and stored at CSIRO, Floreat, Western Australia.

Individual BAC clones were grown overnight at 37°C with vigorous shaking in 5-10 mL LB containing 12.5 μg/mL chloramphenicol, and the cultures were used to prepare BAC DNA by alkaline lysis. DNA of each clone was digested with 0.2 unit Not I restriction enzyme (New England BioLabs) and subjected to PFGE in 1% agarose gels and 1 × TAE buffer, for 18 h at 11°C and 6 V/cm, using a 1-40 s pulse time and 120° field angle. Fragments were photographed under UV light after ethidium bromide staining and insert sizes estimated by comparison to a Lambda Ladder PFG Marker (New England BioLabs).

BAC library screening

A cDNA probe from the NLL β-conglutin gene (Beta2; [38]) was used for screening the BAC library. The probe was derived from the plasmid DNA containing the cDNA fragment of NLL Beta2 by PCR using M13 forward (5'-GTTGTAAAACGACGGCCAGT-3') and reverse (5'-CAGGAAACAGCTATGACC-3') primers. After purification with the Wizard SV 96 PCR Clean-up System (Promega), about 600 ng of the PCR product was labelled with [α-³²P]-dCTP using Ready-To-Go DNA labelling beads (-dCTP) (Amersham Biosciences, Uppsala, Sweden) following the supplier's protocol. After incubating for 30 min at 37°C, the unincorporated nucleotides were removed using Illutra ProbeQuant G-50 Micro Columns (GE Healthcare).

Hybridisation was carried out at 68°C for 2 hr in ExpressHyb Hybridisation Solution (Clontech) following the supplier's protocol. The membranes were washed in 2 × SSC for 30 min and 1 × SSC for 1.5 hours. Subsequently, the membranes were exposed to an Imaging Screen K (Bio-Rad) for 30 min and the images captured using the Molecular Imager FX System and Quantity One software (Bio-Rad).

Characterisation of BAC clones

Twelve clones were randomly selected from those which showed strong hybridisation signals. These clones were used in PCR reactions with a pair of primers designed from the conserved regions of the EST sequences of the seven β-conglutin genes [38]. The primer sequences were: forward 5'-TCCTCGTTGTACTCAATGGT-3' and reverse 5'-GGTTAAGGATATAAGAAGT-3' and the presence of β-conglutin sequences confirmed.

Eight BAC clones were selected for further characterisation. PCR reactions using primers specific for each β-conglutin gene and annealing temperatures (as described in [38]) were carried out with each BAC clone. PCR products were separated by gel electrophoresis and reactions which yielded a band sequenced using the PCR reaction primers. Where two bands were seen from one PCR reaction, each band was excised and purified from the agarose gel using the Qiaquick Kit (Qiagen) following the manufacturer's instructions. Retrieved sequence data was used to identify the best aligned NLL β-conglutin by BLAST analysis.

BAC-end sequencing

The BAC-end sequencing was carried out by the Istituto di Genomica Applicata, Italy. BAC DNA was prepared from randomly selected clones using the MultiScreen Plasmid384 Miniprep Clearing plates (Millipore) coupled with a Biomek FX robot (Beckman). Sequencing of BAC-ends used an ABI3730xl platform (Applied Biosystems) with M13 forward and M13 reverse primers with Big-Dye Terminator chemistry (Applied Biosystems).

Sequence cleaning and trimming

BAC-end sequences (BESs) were trimmed and filtered for quality and sequence contamination using seqclean (http://compbio.dfci.harvard.edu/tgi/software/), which had been modified to use the same BLAST parameters as NCBI VecScreen (http://www.ncbi.nlm.nih.gov/VecScreen/VecScreen.html). BESs were screened against contaminant databases including the UniVec vector database (Build #5.2 - Dec 7, 2009, ftp://ftp.ncbi.nih.gov/pub/UniVec), E. coli BL21 (DE3) complete genome sequence [NCBI NUC: CP001665.1] and plant organelle genomes (all NCBI mitochondrial, chloroplast and plastid genome sequences under the taxon "Eudicotyledons", [NCBI taxid: 71240]).

Estimation of genome statistics, gene content and function

G:C content of BESs was calculated using custom perl scripts. Repetitive DNA content was estimated using RepeatMasker (version 3.2.9, http://www.repeatmasker.org) [56], with the REPBASE repeat database (version 20090604) [57]. RepeatMasker parameters were set to slow/sensitive alignment and restricted REPBASE entries to those of species within the taxon "Eudicotyledons" [NCBI taxid: 71240].

Trimmed NLL BESs were compared to the NCBI NR protein database by BLASTx [58]. BESs regions matching NCBI NR proteins and not matching to RepBase repeats [57] were assumed to contain protein-encoding genes. The proportion of filtered BLASTx-aligned BES regions relative to total BES sequence was used to estimate the percent of the NLL genome containing endogenous protein-encoding genes.

The functional content of the genome was estimated by running the 13985 trimmed BESs through the Blast2GO pipeline using default settings, blastx alignment to NR and ANNEX augmentation [59]. Gene Ontology summary PCA plots were generated via ReviGO [60]. REVIGO analysis focussed on functional annotations within the Lupinus angustifolius BES dataset only and did not involve multiple species comparisons. Assigned GOs were summarised via the 'plant' GOSlim subset (http://www.geneontology.org/GO_slims/goslim_plant.obo). Unfiltered GO terms were also summarised by conversion to MIPS FunCAT terms (http://www.geneontology.org/external2go/mips2go).

Comparative genomics between NLL and other legumes

BESs were aligned via BLASTn [39] to the NCBI Nucleotide database. Phylogenetic distribution BLAST hits was visualised with MEGAN according to the lowest common ancestor of BLASTn BES hits [40]. The phylogenetic relationships in the cladogram presented in Figure 3 were derived from the NCBI taxonomy database.

Screening of BES regions against published markers

Previous studies had published a total of 1104 markers used for the genetic mapping of NLL (summarised in Additional File 10). Of these 892 were PCR primer-based, comprising 74 AFLP, 646 MFLP, 159 RFLP, 6 isosyme and 7 phenotypic markers. The remaining 212 were tested for PCR amplification of BES sequence regions in-silico between a range of 20 to 10000 bp [53].

Identification of SSRs and design of SSR-flanking primers

For the purposes of comparison with previous studies of SSRs in legumes, NLL SSRs were predicted via tandem repeats finder [61] according to the criteria outlined by Mun et al. [41]. SSRs were required to have a total length of ≥ 12 bp, SSR unit length ≥ 1-8 bp and 100% identical repetition of the SSR unit. Predicted SSRs were also distinguished into two classes according to total SSR length as per Mun et al. [41]: class I: SSR length ≥ 20 bp, class II: 12-19 bp. For the purpose of polymorphic marker discovery, additional Class I SSRs were predicted using less stringent criteria: allowing for a SSR unit size of 2-5 bp, a minimum of 75% identity to the SSR unit and a minimum number of five repetitions.

SSR-flanking primers were designed using primer3 [42]. Primer design was optimised for pairs producing 250 bp amplicons, primer lengths of 22 bp, melting temperatures (Tm) of 60°C, maximum difference in forward and reverse Tm of 1°C, maximum allowable primer hairpin lengths of 3 bp and maximum primer mononucleotide repeats of 2 bp. BESs were screened for potential repetitive DNA and organelle genome sequence contamination via sensitive RepeatMasker alignment to the Eudicotyledon subset of RepBase and BLASTn match (e-value < 1e-20) to NCBI Eudicotyledon organelle sequences. BAC-end sequences were also assembled via Cap3 [62] to identify potentially repeated regions (i.e. multi-gene families/uncharacterised repeats). Primer pairs designed within "repetitive" BESs were excluded from consideration as SSR marker candidates.

Multiplex Ready PCR assays to determine SSR lengths in four NLL lines

DNA was isolated using the CTAB method [63] and dissolved in 10 mM Tris HCl (pH 8.0). The forward and reverse primer for 24 class I primer pairs were synthesised with the added nucleotide sequence 5'-ACGACGTTGTAAAA-3' and 5'-CATTAAGTTCCCATTA-3' respectively and called "locus specific primers". A list of the 24 class I SSR primer sequences is presented in Additional file 11. Two generic tag primers, tagF and tagR with the sequences 5'-ACGACGTTGTAAAA-3' and 5' CATTAAGTTCCCATTA 3', respectively, were also synthesised. The tagF primer was labelled at its 5'-end with one of the fluorescent dyes: VIC, FAM, NED and PET (Applied Biosystems). Multiplex-ready PCRs were subsequently carried out as described by Hayden et al. [64]. The multiplexed SSR PCR products were subjected to fragment analysis on an ABI3730 DNA analyser (Applied Biosystems) according to Hayden et al. [64] and SSR allele sizing used the Genemarker software (SoftGenetics LLC).