Reconstruction of genomic repeats by assembling 454 reads

The process of extracting repeat sequences from the set of 454 sequence reads was implemented using TGICL, a program package originally designed for clustering large EST datasets [25]. The procedure consisted of two steps: (i) clustering the reads based on their mutual similarities into groups of overlapping sequences, and (ii) assembling the reads within the clusters to get longer fragments (contigs) representing consensus sequences. Various clustering and assembly parameters were evaluated in order to get optimal performance with our dataset, which compared to ESTs differed in the short reading lengths and in their considerable sequence variability, reflecting the divergence between individual copies of repeated elements. While the clustering parameters allowed for grouping of relatively divergent sequences, the assembly process was more stringent, and thus typically generated multiple contigs from a single cluster.

The clustering resulted in 22,445 clusters, which were composed from two to thousands of reads. The assembly phase then yielded 25,384 contigs ranging in lengths up to 8,214 bp. Individual contigs were assembled from two to 4,327 reads, and the total number of sequence reads assembled into contigs was 233,303 (73 %). The contigs were characterized by their read depth (RD), expressing the average number of reads assembled over individual positions within the contig consensus sequence, and by their genome representation (GR), calculated as RD multiplied by the contig length. These values allowed us to rate the contig sequences based on their genomic copy numbers and proportion in the genome, respectively. Most contigs were short and composed of only a few reads, whereas 90% of assembled reads were assigned into a relatively small subset of 1,578 contigs with the highest GR, thus corresponding to highly repetitive sequences. Most repeats were represented as sets of overlapping contigs, the number of which was proportional to the sequence diversity of the repeat copies. Among the most important in terms of their genomic abundance were contigs that included coding sequences and conserved LTR regions of Ogre retrotransposons and other LTR-retroelements and of the satellite repeat PisTR-B (Additional file 2). Compared to the three previously sequenced pea Ogre elements [21], this set of Ogre-like contigs showed much higher sequence diversity, suggesting they represent several distinct Ogre subfamilies.

The longest contig (CL2Contig6) represented a 8,214 bp fragment of the rDNA gene cluster, including the complete 18-5.8-25S gene sequences (5,820 bp) surrounded by 3' and 5' parts of large intergenic spacer (Fig. 2A). Comparison to previously published partial pea rDNA sequences revealed its identity with a 1,723 bp fragment of the 18S rRNA gene [GenBank: U43011] and 99.5% similarity to a 620 bp sequence region including the 5.8S rRNA gene and both internal transcribed spacers [GenBank: AY839340].

The contig with the highest genome representation was a 5,097 bp fragment of the LTR-retrotransposon Peabody [GenBank: AF083074]. This contig (CL1Contig2066) had a read depth corresponding to 11,000 copies/1C and was estimated to make up 1.3% of the pea genome. It included a complete internal retrotransposon region surrounded by parts of LTR sequences, which could be further extended by finding and aligning overlapping contigs (Fig. 2B). The internal region contained a gag-pol coding sequence (4,368 bp) devoid of stop codons.

The highest read depth (193 reads, corresponding to about 25,000 copies/1C) was found for a 1,213 bp contig CL1Contig751 representing the LTR sequence of a novel Ty1/copia element designated Ps-copia-1/751. The element reconstruction from ten overlapping contigs resulted in identification of the complete LTR and most of the internal region including open reading frame encoding gag-pol polyprotein (Fig. 2C).

In addition to highly repeated sequences it was also possible to at least partially reconstruct less abundant repeats, many of which were novel to the pea genome. For example, an over 7 kb region of a MuDR-like DNA transposon, including 2 putative coding regions, could be reconstructed from 25 overlapping contigs. MuDR elements were estimated to occur in about 2,200 copies in the pea genome, and similar abundance was also found for another DNA transposon family, En/Spm, for which it was possible to reconstruct a 3 kb fragment of the transposase-coding region (not shown).

The clustering and contig building procedure was also found useful for identifying novel tandemly repeated sequences. Assembling overlapping reads into longer contigs facilitated reconstruction of repeats with monomers exceeding the length of single reads and allowed their identification based on the tandem subrepeats present within the contigs. A number of contigs representing potential satellite repeats with monomers from 50 to 867 bp were identified; except for the previously described PisTR-B satellite [20], they all represented novel sequences. Fourteen of the most abundant repeats (Table 1) were used as probes for in situ hybridization on pea mitotic chromosomes in order to test if they have a genomic distribution typical for satellite DNA. Such hybridization patterns, consisting of signals concentrated into limited number of spots corresponding to long arrays of the satellite sequences, were observed for thirteen repeats, whereas only one produced dispersed signals (Fig. 3A–C and Table 1). The signals occurred mostly in (peri-) centromeric and terminal chromosome regions, and each repeat displayed a specific hybridization pattern. No typical centromeric satellite repeats were found, although the repeat TR-11 strongly labeled central centromeric regions in five out of the seven chromosome pairs (Fig. 3B).

Repeat	Monomer [bp]	Cluster	Abundance(a)		Localization on chromosomes(b)							Note
			[%]	CN	1	2	3	4	5	6	7
TR-1	867	CL14	0.14	7,000	PC	-	C, T	PC,C	-	-	-
TR-2	~440	SCL5	0.29	28,000	-	PC	-	-	C	PC	-
TR-3	82	CL19	0.10	51,000	-	-	PC	I	-	-	-	Fig. 3C
TR-4	172	CL16	0.09	22,000	-	-	-	-	-	PC	-
TR-5	54	CL52	0.04	35,000	-	PC	-	-	-	-	-
TR-6	245	CL65	0.04	7,000	-	-	-	-	PC,C	-	-
TR-7	164	CL49	0.05	13,000	PC,C	-	-	-	-	-	-
TR-8	342	SCL58	0.04	5,000	~	~	~	~	~	~	~	Dispersed
TR-9	189	CL68	0.04	8,000	T	T	-	-	-	-	T	Fig. 3A
TR-10	659	CL78	0.03	2,000	-	-	-	-	PC,C	-	-
TR-11	~510	CL9	0.20	17,000	C	C	-	-	C	C	C	Fig. 3B
TR-12	~120	CL92	0.02	7,000	I	PC, I	I	C, I	I	I	I
TR-14	193	CL124	0.02	5,000	-	-	-	T	-	-	-
TR-17	191	CL193	0.01	3,000	-	-	-	-	-	-	I

^(a) Repeat abundance is expressed as its proportion in the genome (in percents of the genome size) and as copy number of monomers per haploid genome (CN).

^(b) Positions of FISH signals on individual chromosomes (n = 7). C, centromeric; PC, pericentromeric; I, intercalary; T, (sub-)telomeric.

Composition of the most repetitive fraction of the pea genome

Although other retroelement families were found in considerably smaller numbers, there were several elements which made up significant proportions of the genome. They included Peabody, which made up 2–3% of the genome and displayed very low sequence variability suggesting its recent amplification. Other important groups of Ty3/gypsy elements were represented by PIGY [24] and Cyclops [27]. Ty1/copia retrotransposons were found to be less frequent, being represented by PDR [28] and a group of SIRE1-like sequences [29]. However, the most abundant was a novel element, Ps-copia-1/751 (Fig. 2C), which made up about 2% of the pea genome. The LTR sequence of this element was estimated to occur in at least 25,000 copies in the pea genome, whereas other regions are less frequent (about 8,000 copies/1C), thus indicating the existence of a large number of solo-LTRs derived from this element (Table 2).

Overall diversity of retrotransposons was also studied by their classification based on the RT domains, which are conserved enough to be reliably identified based on protein similarity searches with RT domains from already known retrotransposons. The searches resulted in finding 452 contigs with similarity to RT domains (E-value <= 0.001). Of them, 222 contigs were related to Ty1/copia, 217 to Ty3/gypsy, and only 13 to LINE retrotransposons. Only four Ty1/copia-like and seven Ty3/gypsy-like contigs spanned full-length RT domains as defined by [30] while the rest contained partial sequences due to the short length and/or termination within the RT domain. Phylogenetic analysis of RT sequences belonging to Ty1/copia and Ty3/gypsy groups revealed that most Ty3/gypsy-like RT domains are related to the previously described pea retrotransposon families Ogre, Peabody, PIGY and Cyclops (Additional file 4). On the other hand, Ty1/copia-like RT domains were related to a greater number of retrotransposon families, most of which have not been previously identified in the pea genome (Additional file 5). RT domains were also used to estimate the copy number of retrotransposons present in pea. These elements were estimated to occur in about 141,000 copies/1C, of which 46,000 (32.6 %) belong to Ty1/copia, 94,000 (67%) to Ty3/gypsy and ~500 (0,4%) to the LINE group (Fig. 4). The highest copy number was estimated for Ogre retrotransposons, which amplified themselves to about 64,000 copies. More detailed analysis of Ogre-like RT domains confirmed the occurrence of distinct Ogre subfamilies sharing their best similarities with elements from different branches of the Ogre clade (Additional file 4). The most abundant were subfamilies from PS and VP/VM branches having about 24,000 and 32,000 copies, respectively (Additional file 4 and Fig. 4). Thus, the total copy number of Ogre elements as well as the abundance of the two major subfamilies were in agreement with estimates based on the abundance of PBS sites described above. Two other clades of Ty3/gypsy retrotransposons, Peabody and envelope-like retrotransposons (including PIGY and Cyclops), were estimated to have 16,000 and 12,000 copies, respectively. Although Ty1/copia retrotransposons were less frequent than Ty3/gypsy, elements from clades 3 and 5 (Additional file 5) reached high copy numbers (10,000 and 26,000 copies, respectively; Fig. 4). Interestingly, clade 3 contained elements from Graminae species such as barley (BARE-1), rice (RIRE-1), and oats (OARE-1) clustered together with a newly identified pea element Ps-copia-1/751. Clade 5 contains elements similar to SIRE1, an envelope-like retrotransposon from soybean [29]. It is not clear yet, however, whether all elements within this clade are genuine envelope-like retrotransposons, i.e. whether they all bear envelope-like gene.

Compared to the sum of retroelement sequences, other classes of repeats represented significantly smaller parts of the pea genome (Table 2). Identified DNA transposons did not exceed 0.5%, and tandem repeats including rDNA gene clusters and various satellite sequences represented about 2.5% of the pea genome. We have also investigated the abundance of micro- and minisatellite repeats with monomers from 2 to 10 bp. Since these repeats usually occur in the genome as short stretches of repeated monomers dispersed within unrelated sequences, we analyzed their frequency in unassembled sequence reads instead of contigs. When considering arrays of at least five consecutive monomers, microsatellites (AAT)_n, (AT)_n, and (AG)_n were found to be most frequent, occurring in about 80,000 genomic loci. Other microsatellite motifs were less abundant and there were differences over two orders of magnitude in the frequency of various microsatellites within the reads (data not shown).

Special attention was paid to detection of telomeric repeats. The Arabidopsis-type telomeric sequence (TTTAGGG)_n [31] was present in 14 reads, with an average number of 18 repetitions per read. This gives a rough estimate of about 229 kb of this sequence in the pea genome, and agrees with experimental observations of this repeat at termini of all chromosomes (data not shown). However, two additional telomere-like repeats were found in the analyzed reads – (TTAGG)_n and (TTTAGG)_n. Although both were less frequent than the Arabidopsis-type repeat, their sequences also spanned the whole reads, suggesting they occur in the genome as longer contiguous arrays. FISH experiments using labeled oligonucleotide probes confirmed this presumption, revealing the presence of both sequences at the termini of all pea chromosomes (Fig. 3D–E). To check for specificity of this assay, control experiments were run using probe sequences differing in a single-base substitution within the repeat monomer ((TTAC G)_n and (TTTAC G)_n). These sequences were not present in the 454 reads and corresponding probes gave no hybridization signals on mitotic chromosomes (data not shown).

Comparative analysis of the repeat composition of pea, soybean and M. truncatula genomes

In addition to the identification of the major components of pea repetitive DNA, we were interested in using the 454 sequence data for comparing the pea genome composition to that of related legume species. We primarily focused on soybean (Glycine max), as it was the only species for which whole genome shotgun reads produced by the same technology were available. The soybean 454 data included 718,589 reads with an average length of 109 nucleotides [19]. To investigate which sequences are shared between these two genomes we identified and analyzed pea 454 reads producing significant similarity hits (E-value <= 1e^-10) in BLAST searches against a database of the soybean 454 reads. A total of 5,482 pea reads (1.7%) matched soybean sequences; in the soybean dataset, 7,209 reads (1.0%) gave significant hits with pea. The composition of these sequences is summarized in Fig. 5A–B. The largest fraction of the matching reads belonged to rDNA sequences, which also displayed the highest sequence similarity between the species (reflected by high BLAST bit scores on Fig. 5B). Ogre elements, representing the most abundant repeats in the pea genome, were also identified in soybean but in much smaller numbers. On the other hand, the two species were found to share relatively large population of SIRE1 sequences. Interestingly, the other abundant families of pea Ty1/copia elements were absent or gave only a few hits in the soybean genome.

In order to compare genome composition of pea and soybean with the legume model species Medicago truncatula, we employed the same strategy as above, except for using a set of 197,570 M. truncatula BAC-end sequences instead of 454 data. We found that pea genome contains considerably more repeats similar to M. truncatula sequences than soybean. A total of 15,510 (4.9%) pea reads gave significant matches with M. truncatula and included all major pea repeat families (Fig. 5D). However, there were differences in genomic proportion and sequence conservation of individual repeats between pea and M. truncatula genomes. For example, abundant pea Ogre sequences were also frequently found in M. truncatula but their sequence similarity was relatively low, whereas Ty1/copia elements (represented mostly by Ps-copia-1/751 family) produced much fewer hits but displayed higher average sequence similarity (Fig. 5D). Contrary to pea, soybean repeats were only poorly represented in M. truncatula. About 2.1% of soybean 454 reads matched M. truncatula sequences, however, except for rDNA (0.25%) no major repeat family was found to be shared by these two species (Fig. 5C). In summary, these results indicate similarity in repeat composition between pea and M. truncatula genomes, and considerable sequence divergence of most repeat families between these two species and soybean.

Genomic DNA isolation and 454 sequencing

Seeds of garden pea (Pisum sativum L. cv. Carrera) were obtained from the Plant Breeding Station at Boršov, Czech Republic. The DNA was extracted from purified nuclei in order to minimize contamination with chloroplast and mitochondrial genomes. The nuclei were isolated from young leaves by grinding 5 g of the tissue in liquid nitrogen, followed by 5 min incubation in 35 ml of ice-cold H buffer [47]. The homogenate was filtered through 48 μm nylon mesh, adjusted to 35 ml volume with 1 × H buffer, and centrifuged at 200 × g for 15 min at 4°C. Pelleted nuclei were resuspended and centrifuged using the same conditions once in 35 ml of H buffer, and once in 15 ml of TC buffer (50 mM TRIS-HCl pH 7.5, 75 mM NaCl, 6 mM MgCl₂, 0.1 mM CaCl₂). The final centrifugation was performed for 5 min only and the nuclei were resuspended in 2 ml of TC. DNA was released by incubating the nuclei with 40 mM EDTA, 0.2% SDS and 0.25 μg/μl of proteinase K for 4 hours at 37°C and purified by phenol extraction and ethanol precipitation. The sequencing was performed by 454 Life Sciences (Branford, CT, USA) using the GS-20 instrument and yielded 322,396 quality filtered sequence reads with the average length of 104 bp. The reads were deposited into NCBI Short Read Archive [48].

Processing and initial analysis of 454 reads

The reads were screened for their similarity to adaptors and primers used in 454 sequencing using BLAST [49] and 595 reads producing blastn hits with expectation (E) value < 1e^-15 were removed from the set. The same procedure was repeated to detect contamination with organellar sequences, removing 2,322 reads with similarity to plant chloroplast genomes and 77 reads similar to mitochondrial DNA. The remaining 319,402 reads were used for all subsequent analyses. These reads included a total of 33.3 Mb, corresponding to 1/129 of the pea haploid genome size (1C = 4,300 Mb [18]). Thus, the calculated repeat copy numbers were multiplied by 129 in order to compensate for the actual genome coverage.

To determine the representation of previously characterized pea repeats in 454 sequences, these repeats were used as queries in blastn searches against a database of the 454 reads, using E-value cutoff of 1e^-10. BLAST outputs were parsed using a BioPerl [50] script determining the number of similarity hits at all positions along the sequences, calculating their average numbers and the corresponding estimates of genomic copy numbers. The filtered set of 454 reads used for this analysis is available for BLAST searches at our web site [51].

Reconstruction of repetitive sequences

Reconstruction of genomic repeats as well as other computer analysis were performed on an IBM xSeries 226 dual processor server with 4 GB RAM running under the Gentoo Linux operating system. In addition to the programs listed below, several scripts written in Perl were used for data processing and analysis. These scripts are available upon request from authors. The repeat reconstruction was done using TIGR Gene Indices clustering tools (TGICL [25]) employing the following parameters, which were optimized for our dataset by evaluating number and size of resulting clusters, the length of the assembled contigs and the presence of chimeric (misassembled) contig sequences. All-versus-all pairwise similarity scores of the reads were produced by mgblast and parsed using tclust, performing transitive-closure clustering of the read pairs sharing a minimum of 90% similarity over at least 70% of the shorter sequence. The reads within individual clusters were assembled into contigs using tgicl run with the -O '-p 80 -o 40' parameters, specifying overlap percent identity and length cutoff for cap3 assembler. Although most clusters included reads derived from only a single family of genomic repeats, the largest cluster (CL1) was composed from a large number of reads corresponding to several unrelated families of mobile elements. This problem was probably caused by the presence of "hybrid" reads including sequences of two different elements (spanning their insertion sites or recombination breakpoints, see also Discussion). To separate these unrelated sequences, the contigs assembled from CL1 were further grouped into sub-clusters using sclust (parameters HEAVY = 3010 SCORE = 150), which performs weighted seeded clustering based on mutual similarities of the contigs. The programs mgblast, tclust, tgicl, cap3 and sclust are included in the TGICL package and the detailed description of parameters can be obtained from the program help.

The assembled contigs were labeled based on the cluster they originated from (CL [number]Contig [number]); in case of CL1 the name also included sub-cluster number (SCL [number]_CL1Contig [number]). The contig sequences were characterized by calculating their average read depth (RD) and genome representation (GR = RD × contig length). The genome representation was also calculated for individual clusters by summing GR values of the corresponding contigs. Contigs from the clusters representing at least 0.01% of the pea genome (GR > 3,300) were subjected to sequence similarity searches (blastn, blastx) against GenBank database and to detection of conserved protein domains using RPS-BLAST [52] in order to identify the type and family of the repeat they originated from. Satellite repeats were identified by detection of tandem subrepeats within the contig sequences using dotter [53]. Total proportion of identified repeats in the pea genome was calculated by summing the GR of the clusters assigned to the same type of repeat. Detailed information about cluster characterization is available from Additional file 3. All contig sequences with their RD and GR values can be downloaded from our website; alternatively, the contigs can be selected by sequence similarity searches against user-provided queries [51].

Sequence analysis

Phylogenetic analysis of retroelements was done using multiple sequence alignment of reverse transcriptase domains and a phylogenetic tree was calculated by neighbor-joining method (bootstrap values calculated from 1000 replicates) using the ClustalX program [54]. Reverse transcriptase sequences used for the analysis were extracted mostly from plant retrotransposon sequences downloaded from RepBase [55, 56]. Graphical representation of the tree was drawn and edited using ATV Forester [57], TreeEdit [58] and Jalview programs [59]. Copy number (CN) estimation based on RT domains was calculated using formula CN = RD_rt × 129 × L1/L2, where RD_rt is the average read depth within the reverse transcriptase domain, L1 is the length of the reverse transcriptase domain in particular contig and L2 is the length of the full size domain. Contigs bearing only partial reverse transcriptase domains and therefore missing in the phylogenetic tree were assigned to appropriate retrotransposon group based on their best Blast hits to elements shown in the tree (Additional files 4 &5).

An alternative approach for copy number estimation of Ogre elements was based on analyzing contigs including the Ogre-specific PBS sequence complementary to tRNA_Arg [26]. Contig assemblies (ACE files produced by cap3 assembler) were visualized using clview [25] and average read depth over the PBS site was calculated for each contig. The contig sequences were compared using dotter and grouped into three subfamilies according to similarity of the sequences surrounding the PBS. The number of PBS-containing reads was summed for each subfamily and used to calculate its genomic copy number.

Microsatellite and telomeric repeats were identified directly in a set of individual 454 reads instead of in the contigs. The reads separated with runs of "N" (to avoid fusing adjacent sequences) were concatenated into a single sequence and analyzed using Tandem Repeats Finder [60]. The results were sorted and evaluated using TRAP [61].

Comparative analysis of the repeat composition in pea, soybean (Glycine max) and barrel medic (Medicago truncatula) was performed using our set of 319,402 pea 454 reads, a set of recently published 718,589 soybean 454 reads [19], and 197,570 M. truncatula BAC-end sequences. The soybean and M. truncatula sequences were downloaded from the NCBI Trace Archive [48] and from TIGR [62], respectively. Sequence similarities were identified by blastn searches between 454 reads from pea and soybean, and between either pea or soybean 454 reads and M. truncatula BAC-end sequences. The 454 reads producing blastn hits with E-value of 1e^-10 or lower to sequences from the other species were considered to include sequences shared between these genomes. For each of these reads, the number of hits to other genome sequences and average blastn bit score of these hits were determined and plotted using Mgraph (created by Louis Gonzalez and Christine Deroo, available from [63]). Grouping the reads according to the repeat type was done by similarity searches to known sequences performed as described above.

Fluorescence in situ hybridization (FISH)

Hybridization probes for potential satellite repeats were prepared by PCR amplification from pea genomic DNA using primer pairs designed according to the contig sequences including tandem subrepeats (Table 1 and Additional file 6). The PCR was performed in 30 μl of reaction mix (1× PCR buffer, 1.5 mM MgCl₂, 0.2 mM dNTPs, 0.2 μM primers, 5 U of Taq polymerase (Promega), 0.2–20 ng of template DNA) for 25 cycles of 1 min at 94°C, 1 min at 50°C and 3 min at 72°C, preceded by initial denaturation (3 min at 94°C) and followed by final extension step (10 min at 72°C). Reaction products were resolved on agarose gel electrophoresis and bands correspoding to the repeat monomers or short multimers (< 800 bp) were excised from ladders of the amplified fragments. The excised DNA was purified and used as a template in PCR labeling by biotin-dUTP as described [20]. 5' biotin-labeled oligonucleotide probes for alternative telomeric repeats (CCTAA)₈, (CCTAAA)₇ and controls (CGTAA)₈, (CGTAAA)₇-CG were purchased from Bioneer Europe. FISH experiments were carried out on chromosome squash preparations [64] obtained from root tip meristems synchronized as described [65]. The chromosome preparations were aged for 5–60 days and post-fixed in 4% formaldehyde/2 × SSC for 10 min at 25°C. The hybridization was performed as described [64] with the following modifications. Probe concentration in the hybridization mix was 0.2 μM for oligonucleotide and 0.5 ng/μl for PCR-amplified fragments. Formamide concentrations in hybridization and washing solutions and incubation temperatures were adjusted according to the probe sequences to achieve 90% stringency for the oligonucleotide probes (Tm calculated using Exiqon Tm prediction 1.1 software [66]) and 80% stringency for the satellite repeat probes. Stringency for both types of probes was calculated according to [64]. Biotin-labeled probes were detected with streptavidin-Alexa Fluor 568 (Invitrogen) and individual chromosome types were distinguished by co-hybridization with Alexa Fluor 488-labeled PisTR-B probe, producing characteristic patterns on each pea chromosome [20]. Chromosomes were counterstained with 4,6-diamino-2-phenylindole (DAPI) and observed using a Nikon Eclipse-600 epifluorescence microscope equipped with a CCD camera. The signals were collected using appropriate filter sets and LUCIA software (Laboratory Imaging).