Development of genomic resources for Citrus clementina: Characterization of three deep-coverage BAC libraries and analysis of 46,000 BAC end sequences

Background Citrus species constitute one of the major tree fruit crops of the subtropical regions with great economic importance. However, their peculiar reproductive characteristics, low genetic diversity and the long-term nature of tree breeding mostly impair citrus variety improvement. In woody plants, genomic science holds promise of improvements and in the Citrus genera the development of genomic tools may be crucial for further crop improvements. In this work we report the characterization of three BAC libraries from Clementine (Citrus clementina), one of the most relevant citrus fresh fruit market cultivars, and the analyses of 46.000 BAC end sequences. Clementine is a diploid plant with an estimated haploid genome size of 367 Mb and 2n = 18 chromosomes, which makes feasible the use of genomics tools to boost genetic improvement. Results Three genomic BAC libraries of Citrus clementina were constructed through EcoRI, MboI and HindIII digestions and 56,000 clones, representing an estimated genomic coverage of 19.5 haploid genome-equivalents, were picked. BAC end sequencing (BES) of 28,000 clones produced 28.1 Mb of genomic sequence that allowed the identification of the repetitive fraction (12.5% of the genome) and estimation of gene content (31,000 genes) of this species. BES analyses identified 3,800 SSRs and 6,617 putative SNPs. Comparative genomic studies showed that citrus gene homology and microsyntheny with Populus trichocarpa was rather higher than with Arabidopsis thaliana, a species phylogenetically closer to citrus. Conclusion In this work, we report the characterization of three BAC libraries from C. clementina, and a new set of genomic resources that may be useful for isolation of genes underlying economically important traits, physical mapping and eventually crop improvement in Citrus species. In addition, BAC end sequencing has provided a first insight on the basic structure and organization of the citrus genome and has yielded valuable molecular markers for genetic mapping and cloning of genes of agricultural interest. Paired end sequences also may be very helpful for whole-genome sequencing programs.


Background
Citrus, one of the major fruit tree crops is widely cultivated throughout the globe and therefore has a tremendous economical, social and cultural impact in our society. Citrus improvement through traditional techniques, however, is highly impaired due to the unusual combination of biological characteristics of Citrus species, their low genetic diversity and the long-term nature of tree breeding. Citrus are diploid plants with an estimated haploid genome size of about 367 Mb and 2n = 18 chromosomes, which may facilitate the use of genomics tools for crop improvement. Expressed sequence tag (EST) analyses and molecular marker studies strongly suggest that the main commercial citrus cultivars (oranges, lemons and grapefruits) are mostly interspecific hybrids and therefore are heterozygous "species" [1]. In addition, most of the cultivars in these groups, including Clementine varieties, may represent accumulated somatic mutations identified over centuries [2].
The development of citrus genomic resources is in its infancy although in recent years major efforts and goals mostly on functional genomics have certainly been undertaken [3]. Critical functional and expression analyses through microarrays with several platforms have also been published and analyses of ESTs in public databases have been initiated [4,5]. For instance, 401,692 citrus ESTs have been deposited at GenBank and are currently available. This collection constitutes a valuable source for the direct access to the genes of interest and for the development of molecular markers for map-based cloning purposes or marker-assisted selection programs [6][7][8][9]. Moreover, genetic linkage maps have been produced with increasing value and resolution, following the evolution of new marker systems [10,11]. Genetic transformation in citrus is also available [12] and strategies based on genome-wide mutagenesis are being explored. Other innovative resources such as viral-induced gene silencing (VIGS) are being developed and work in citrus proteomics is in progress [13,14]. Thus, current advances in citrus research include the rapid development of functional genomics and molecular biology resources [15] although, on the other hand, basic information on the organization and structure of citrus genome is lacking. The main challenge for a comprehensive and meaningful description of genomes is the integration of the DNA marker-based genetic maps with physical maps, and eventually with DNA sequence of the whole genome, the ultimate physical map. For the generation of high-resolution physical maps, the construction of BAC libraries containing clones with large DNA fragments appears to be indispensable. BAC end sequencing is indeed an important component of physical map development and can be considered a form of low coverage sequencing [16]. Paired end sequences of BACs form an important part of scaffolding whole-genome shotgun programs [17] as well as in BAC based genome projects [18,19]. In addition, BAC clone collections and BAC-based contig maps are powerful tools having multiple applications in genomics including positional cloning. The BAC end sequence provides a random survey of the information contents (genes, transposons, repeats) of unsequenced genomes [20][21][22], and yields molecular markers useful for genetic mapping [23][24][25], and cloning of genes of agricultural interest [26][27][28]. Furthermore, in many agriculturally important species BAC clones and physical maps are being rapidly developed since they are essential components in linking phenotypic traits to the responsible genetic variation, to integrate the genetic data, for the comparative analysis of genomes, and to speed up marker-assisted selection (MAS) for breeding. Thus, BAC libraries have become central for physical mapping, genome analysis, clone based sequencing and sequencing of complex genomes, for both model [18,19] and main crop plants [29][30][31][32].
In citrus, two BAC libraries from Poncirus trifoliata [11] and a hybrid of Citrus × Poncirus [10] have been described in detail. These libraries were constructed as part of a mapbased cloning strategy of genes conferring resistance to citrus tristeza virus that causes significant economic damage and losses to citrus worldwide. Poncirus is a non domesticated genus related to citrus species that produces inedible fruit. However, other efforts to generate BAC citrus resources, for instance in Satsuma or sweet orange have also been accomplished [3].
In this work, we report the characterization of three genomic BAC libraries from Clementine (Citrus clementina) mandarin, a cultivar of great [3]economic importance that has been a main target of recent studies [15]. The development of these new genomic tools also complements the functional genomics platform generated for this species including an extensive EST collection and a 20 k cDNA chip [4,5], expanding further possibilities for isolation of genes of agronomical interest [33,34]. This study provides the most comprehensive, large insert clone resource of any Citrus species and reports the analysis of 46,339 BAC end sequences offering a first detailed insight into the sequence composition of the Clementine genome. The analyses focused on protein coding regions, repeat element composition, microsatellite and single nucleotide polymorphism contents. Additionally, data on gene homology based comparative genomics with poplar and Arabidopsis are also presented. The annotated BACend sequences may well serve as useful resources for physical mapping, positional cloning, genetic marker development and genome sequencing of C. clementina. Three genomic BAC libraries (CCL1, CCER1 and CCH3) were constructed as described in Material and Methods, with DNA from Citrus clementina (var. Clemenules), pEBAC1 as the cloning vector and three different restriction enzymes ( Table 1). The CCL1 library composed of 19,200 clones was generated with MboI partial digestions. EcoRI and HindIII were used for the construction of CCER1 and CCH3, respectively, and 18,432 clones were picked from each one. The three libraries contained 56,064 BAC clones that were arrayed in 146 384-well microtiter plates. It has been reported that the use of three different restriction enzymes resulted in a more accurate coverage of the genomes, since the different GC contents of their recognition sites increases the representation of a higher number of genomic regions [35]. A single library constructed with one restriction enzyme usually cannot provide a full coverage of the genome, as the restriction sites are not uniformly distributed along the genome, and therefore genomic regions having too many or too few restriction sites are not equally represented [36]. Generally, two or more complementary large insert libraries, constructed with different restriction enzymes, have been successfully used for physical mapping of several plant genomes including those of Arabidopsis [37], japonica rice [38], or soybean [30].

BAC library characterization
In order to evaluate the average BAC insert size, 362 BAC clones (about 120 clones from each library) were randomly chosen and the corresponding DNA was extracted, digested with the rare cutter NotI enzyme and analyzed by PFGE. All fragments generated by NotI digestion contained the 8.7 kb vector band and various insert fragments ( Figure 1). The estimated insert sizes ranged from 10 to 330 kb, with an average of 124 kb for CCL1, 132 kb for CCH3 and 127 kb for CCER1. Since the haploid genome size of C. clementina is in the order of 367 Mb, the libraries coverage is predicted to be 19.5 haploid genome equivalents while the probability of finding any specific sequence is greater than 99.999%. It has been estimated that the number of clones representing 10× haploid genomes is adequate for most genome research purposes, including physical mapping [39].

BAC end sequencing
A total of 28,032 BAC clones from the three genomic libraries were selected for end sequencing as described in material and methods and out of this number, 24,221 clones (86% success rate) rendered 46,339 BAC end sequences (BESs). The three libraries contributed approximately with similar number of reads. The average read length was 652.8 bp and the genomic raw sequence produced was 28.6 Mb, which corresponds to almost 8% of the Clementine genome. Table 2 shows a summary of the BAC end sequencing features. The 46,339 BESs were deposited at GenBank with accession numbers from ET068227 to ET114565.

Sequence Annotation
Chloroplast and mitochondrion DNA analysis In order to identify extranuclear sequences, BES were first compared against the Citrus sinensis chloroplastic [40] and the Arabidopsis thaliana [41] mitochondrial genome sequences with an stringent threshold of 1e-15. The comparison indicated that 736 (1.68%) and 46 (0.1%) sequences produced significant matches with the chloroplastic and mitochondrial genomes, respectively. Chloroplastic BESs were assembled and the consensus sequences obtained spanned 101,470 bp, approximately 70% of the chloroplast genome. Chloroplastic and mitochondrial DNA summed up to 480 kb, and therefore, the total genomic DNA obtained was 28,1 Mb. Considering the average size of the BAC clones, the coverage provided by the BAC end sequencing was higher than 8.4 genomic equivalents.

Repetitive DNA Analysis
The repetitive DNA fraction present in Clementine BESs revealed with RepeatMasker, included 9,618 interspersed repeats which extended over 2.55 Mb, 8.95% of the total raw sequence (Table 3). BLASTX search performed to identify coding regions (see next section), showed that 2,173 additional sequences not detected by RepeatMasker also presented high significant similarity to transposable elements (TEs) and, therefore, are also part of the repetitive DNA fraction of the genome. Thus, the number of BESs with interspersed repeats rose to 11,791, approxi-mately a 25% of the total reads, a percentage between those found for Carica papaya (16%) [20] or Musa acuminata (36%) [21]. No significant differences between the 3 BAC libraries were found when the number of BESs carrying repetitive elements was compared (see Additional File 1). The sequence length occupied by transposon elements (TEs) including the additional reads identified through homology search was in this way increased to 3.58 Mb, a fraction corresponding to 12.6% of the total BAC end sequences. This fraction was rather similar to the ratio reported for Brassica rapa (13.8%) in an estimation also based on partial sequencing [22]. Comparisons with the percentages found in fully sequenced genomes showed that the relative occurrence of TEs in Clementine was also  similar to the fractions found in Arabidopsis (10%) [18] and black cottonwood (12.6%) [17] and lower than in rice (35%) [19] and grapevine (38.8%) [42]. Although the proportion of the different TEs largely varied, in all these species, as well as in C. clementina, class I elements (including LINE, Gypsy-like and Copia-like elements) were predominant over class II (including CACTA, MULE and hAT elements). In comparison with the above sequenced plants, the occurrence of Mutator-like, LINEs, Cacta, and Gypsy elements, represented as the percentage of occupied sequence, was in general lower in the genome of C. clementina. In contrast, copia-like elements were relatively more abundant in Clementine than in Arabidopsis, rice and poplar ( Figure 2). The abundance of these elements in citrus has previously been estimated to be relatively high (13%) [43], while the data obtained through partial sequencing in this study suggested a lower preponderance (3.9%; Table III). Copia and gypsy like elements, however, are transcriptionally active in Clementine [43,44] and therefore, could be an important source of genetic variability in this species.
In order to identify low complexity repeats, all BESs were searched against themselves with BLASTN and then classified based on the number of significant hits produced ( Figure 3). After filtration, 17,585 BES producing at least one hit different from themselves were clustered in this way. The results showed that while a high proportion of BESs (82%) displayed a low number of hits, 3,221 reads produced more that 10 hits, suggesting that these BESs may carry non-coding repetitive sequences, i.e. they may be interspersed repeats of lower complexity. The amount of sequence occupied by these repeats was estimated to be 1.12 Mb, 3.94% of the analyzed sequence. This proportion appears to be moderate in comparison with the 23.5% figure reported for poplar [17].  [19], P. trichocarpa [17], and A. thaliana [18] are included.

Analysis of coding regions
Putative low complexity repetitive sequences identified through BLASTN search of BESs against themselves Figure 3 Putative low complexity repetitive sequences identified through BLASTN search of BESs against themselves. The figure represents the frequency distribution of BESs as related to the hit number obtained after BLASTN search of BESs against themselves. BESs producing more than 50 copies were considered to be putative interspersed repeats.
the estimated genome size (367 Mb), and the number of genomic equivalents analyzed (8.4), the gene contents of the genome of Citrus clementina was assessed as 31,000. This gene number is comparable to the estimate reported for three species of similar genome size that have been sequenced to completion to date: rice (Oryza sativa), with a genome size of 430 Mb and 41,042 genes identified [45], Black cottonwood with 55,000 genes in 485 Mb [17], and grapevine (Vitis vinifera), with 30,434 genes in 487 Mb [42]. Other estimates based on BES analysis obtained for plant species were also analogous. For instance, in papaya, the estimated gene contents was 35,526 with a 372 Mb genome size [20], and in Chinese cabbage, with a genome of 529 Mb, 43,000 genes were calculated [22].
In the Clementine SSR set, there were 758 class I (more than 10 repeats) and 3,056 class II (less than 10 repeats) SSRs, with di and trinucloetides accounting for almost 70% of the SSRs, while tetra and pentanucleotides were less represented. In general, those motifs containing A/T nucleotides were far more abundant than G/C rich repeats, specially ATT/TAA and AT/TA tri-and dinucleotides ( Figure 5). A similar distribution was found by Jiang et al [8] and Chen et al. [51] in the analysis of 8,218 and 3,278 citrus SSRs derived from ESTs, respectively. In the Clementine genome, microsatellites were more numerous in non coding sequences (56% of the SSRs) than in putative coding regions (44%), as previously reported in papaya [20], Chinese cabbage [22], and Arabidopsis [18].

GC content in BES of Citrus clementina BAC libraries
Microsatellites are co-dominant, highly polymorphic, and simple to use markers that have been successfully used in studies of genetic diversity [7] and genetic mapping in citrus [8,52] and in many other plants [25,[53][54][55]. The additional markers reported here will certainly contribute to improve the coverage of the Citrus genome for many purposes, including the development of accurate linkage and genetic maps.

Contig Assembly and SNP analysis
Assembly of BESs that did not contain repetitive sequences was performed with CAP3 [56], and a total of 6,461 contigs including 19,057 reads and covering 6.14 Mb of sequence were produced. It has been suggested that C. clementina is an offspring of a C. sinensis × C. reticulata cross and therefore has a heterozygous genome [1]. In order to identify possible polymorphisms affecting the assembly of the readings and the construction of the physical map of this species, a BLASTN search of all the contigs and singlets was performed against themselves. One hundred thirty sequence pairs that presented a single reciprocal BLASTN hit and the same protein hit from the previous BLASTX search (identical accession number) were selected for further analysis. A total of 81 pairs displayed sequence identities higher than 90% and were considered as originated from the same genomic region. These sequences were not assembled in the same contig due to the presence of polymorphisms, similarly to what was reported in the sequencing of the highly heterozygous grapevine variety, Pinot Noir [57].
The contigs generated in the assembly were analyzed in order to identify single nucleotide polymorphisms (SNPs) in the diploid genome of C. Clementina. SNPs are the most abundant and powerful polymorphic markers, since they provide gene-based markers that can be used in the creation of dense genetic linkage maps [58] and, more important, in the identification of genes associated with specific trait loci. BAC end sequences of heterozygous genotypes have also been successfully utilized in SNP discovery and construction of linkage maps [59][60][61]. PolyBayes [62], a software designed to use genomic sequence as a template and base quality values to discern true allelic variations from sequencing errors, was used to reveal SNP polymorphisms. The number of putative SNPs identified in the Clementine sequences that showed P ≥ 0.9 and SNP depth lower than 10 was 6,617, corresponding to 1.08 SNPs per kb. PolyBayes has been successfully utilized in automated high throughput identification of SNPs from EST collections, and it has been shown that when using P ≥ 0.95, 85% of the predicted SNPs are generally validated experimentally [63].
In order to test the accuracy of the SNP prediction carried out, a total of 30 polymorphisms were randomly chosen for experimental validation. Primers were designed on the consensus sequence of the contigs, in order to amplify the region containing the SNP. Out of 30 genomic regions, 3 produced two or more bands in the PCR amplification or yielded sequence reads with a mixture of templates. This observation constitutes a first indication of the level of heterozygosity of the Clementine genome. Furthermore, the sequence analysis of the remaining SNP candidates showed that 24 out of the 27 hypothetical polymorphisms, an 88.9% success rate, were certainly validated (see Additional File 3). Considering the reliability of the prediction method, more than 5000 putative polymorphisms with P ≥ 0.95 were found in this work. Transitions were the most abundant changes (3,546; 53.6%), followed by transversions (2,162; 32.7%) and indels (909; 13.7%) ( Table 4). The transition fraction found in poplar (70%) was substantially higher [17]. The predominant transversion was A/T at a frequency that doubled that of C/G changes (Table 4), an unexpected observation that remains to be explained.
It should be noted that the SNP listing reported in this work constitutes the first set of putative SNPs identified in any Citrus species, and hence provides a completely new resource for genome analysis in this genera. It is also worth mentioning that 4,500 SNPs were located in or close to putative coding regions, and therefore these 'functional SNPs' may provide an inestimable resource for the identification of genes associated with specific trait loci in addition to their utility as molecular markers for genetic and comparative mapping, nucleotide diversity analysis and association studies.
Two heterozygous genomes have been sequenced to completion, the Nisqually-1 poplar and the Pinot Noir grapevine strains. In both cases the frequency of polymorphisms found within these heterozygous genomes, was 2.6 [17] and 4 [57] SNPs per kb, respectively. These rates are in contrast with the frequency of 1.08 SNPs per kb found in Clementine. The reproductive biology of the different species could explain these differences. Gametophytic self-and cross-incompatibility, and apomixy would produce low variability within Citrus species [7], while outcrossing by means of insect and wind pollination, which is the norm for poplar and vitis, would result in highly heterozygous cultivars [17,57].

Comparative Genomics
The non-repetitive fraction of the BESs was also used in a BLASTN search against the complete nucleotide sequence of the genomes of A. thaliana, P. trichocarpa and O. sativa, with 1e-14 as cut off value. The genomic sequences were displayed with chromosomes as single searchable fasta sequences. In order to map the BESs unambiguously on the heterologous complete genomes, only those sequences producing single significant hits were taken into account. Table 5 shows that the Populus genome not only yielded the largest number of significant hits (3-fold more than Arabidopsis and almost 5-fold than rice), but also spanned more than twice the length of the sequence displaying similarity.
The 1567 BESs that produced significant hits with poplar were mapped on the chromosomes of this species. The representation drawn in Figure 6 showed that citrus sequences were rather uniformly widespread on the 19 poplar chromosomes, an observation that can also be deduced from data in Table 6 that in addition shows that the average number of tags per chromosome was 83 while the distance between tags was 213 kb. Considering that both species have similar genome size, the uniform distribution of the Clementine tags on the poplar chromosomes may suggest that the citrus genome is conveniently represented in the BAC clone set.
Following the approach of Lai et al. [20], we used forward and reverse BES read pairs separated by the approximate length of BAC clone inserts (~120 kb), to analyze the microsynteny between C. clementina and Arabidopsis, rice and poplar. To be considered as potentially collinear with the target genome, the citrus mate pairs had to map in the heterologous genome into a region comprised between 10 and 300 kb and be also oriented properly. The analyses of the sequences identified 108 Clementine BAC end pairs that met these criteria in poplar, while no one was found in Arabidopsis or rice. Furthermore, the majority of these BES pairs mapped on the Populus genome at a distance similar to the insert size of the Clementine libraries, suggesting the microsynteny between citrus and poplar is higher than between citrus and Arabidopsis. These results are striking since C. clementina and A. thaliana belong to Sapindales and Brassicales orders (eurosids II clade) that probably split approximately 87 MYA, while P. trichocarpa belongs to Malpighiales order, (eurosids I clade) that diverged from eurosids II around 109 MYA [64]. Moreover, similar results were obtained by Lai et al [20] with papaya, that also exhibited higher level of colinearity with the poplar than with the Arabidopsis genome despite that C. papaya is a basal member of the Brassicales. Although a definitive explanation has not been provided yet, it is currently believed that the genome of A. thaliana has undergone a recent whole genome duplication, followed by subsequent gene losses and extensive local gene duplications [18], which might be responsible of the lack of colinearity with other eurosid II species. Comparative genomics with the recently sequenced genome of the grapevine, provides additional evidence that the genome of Arabidopsis has been thoroughly rearranged as related to an ancient angiosperm genome [42]. The fact that papaya, Clementine, grapevine and poplar are long lived, clonally propagated, woody plants, might apparently cause a deceleration of their molecular clocks, resulting in  Representation of the mapping of C. clementina BESs hits on poplar chromosomes and a collection of the first 6,617 putative SNPs described in citrus that may be very useful for positional cloning, genetic and comparative mapping, nucleotide diversity analysis, and association studies. Finally, comparative genomics through gene homology searches has shown that, in spite of their taxonomic classification, microsynteny between Citrus and Populus is higher than with Arabidopsis that is a phylogenetically closer species.

Clementine genotype
Citrus clementina (Clementine mandarin, var. clemenules) developing leaves were used for BAC library construction.

BAC library construction
BAC libraries were constructed from high molecular weight (HMW) genomic DNA processed at Amplicon Express, (Pullman, Washington) using the method described in [65]. DNA digestion was performed with varying amounts of MboI, EcoRI, and HindIII to identify appropriate partial digestion conditions. pECBAC1 and contained two FRT and one oriV elements, thus resulting in the pBAC(FRT-oriV) vector [66]. Ligations were transformed into DH10B Escherichia coli cells (Invitrogen) and plated on LB agar with chloramphenicol (30 μg/ml), X-gal (20 mg/ml) and IPTG (0.1 M). Clones were robotically picked with a Genomic Solution G3 into 384 well plates containing LB freezing media. Plates were incubated for 18 h, replicated and then frozen at -80°C. The replicated copy was used for BAC end sequencing.

Insert size estimation
To estimate insert sizes, 10 μl aliquots of BAC miniprep DNA were digested with 5 U of NotI enzyme for three h at 37°C. The digestion products were separated by pulsedfield Weld gel electrophoresis (CHEF-DRIII system, Bio-Rad) in a 1% agarose gel in TBE buffer 0,5×. Insert sizes were compared to those of the Lambda Ladder PFG Marker (New England Biolabs). Electrophoresis was carried out for 18 h at 14°C with an initial switch time of 5 s, a Wnal switch time of 15 s, in a voltage gradient of 6 V/ cm.

BAC End Sequencing
BAC clones were inoculated into 96-deep well macroplates and grown for 20 hs at 37°C. Cells were harvested by centrifugation and BACs were purified in 96-well plates by a standard alkaline lysis protocol developed by Genoscope (Paris, France). BAC DNA was precipitated with isopropanol and washed with 70% ethanol. Sequencing was carried out on ABI3730 equipment with "Dye Terminator" process using ABI kit version 3.1. in the Genoscope facility.

Bioinformatics
The software phred [67]was used for base calling, and Crossmatch for vector masking. Repetitive DNA was identified with the RepeatMasker software [68], using the viridiplantae section of the RepBase Update [69] as database. Assembly was performed with CAP3 [56], using read quality and default parameters. Similarity searches were performed with the standalone version of BLAST [51,70], against the NCBI non redundant protein, nucleotide and EST databases available on November 2007 [71]. Parsing of the BLAST results was performed with the Bio::Sear-chIO module from the Bioperl package [72]. Coding sequences were annotated with GO terms using Blast2GO [48]. SPUTNIK [50] was used to identify simple sequence repeats (SSRs), and POLYBAYES [62] to search for SNPs. PCR product purification was done directly or after cutting single bands on agarose gel, using respectively QIAquick ® PCR Purification Kit and QIAquick ® Gel Extraction Kit (Qiagen) Sequencing was carried out on ABI3730 equipment with "Dye Terminator" process using ABI kit version 3.1. SNPs were identified as double peaks by manual inspection of chromatograms, and subsequently validated