Plant material

We used Triticum aestivum (2n = 6× = 42, AABBDD), var. Chinese Spring; T. timopheevii (2n = 4× = 28, AAGG), K-38555; Ae. speltoides (2n = 2× = SS), K-389; Ae. tauschii (2n = 2× = DD), К-1662; and T. urartu (2n = 2× = AA), IG45298*. The accessions were from the wheat germplasm collection of the N.I. Vavilov All-Russian Institute of Plant Genetic Resources RAN (St Petersburg, Russia) and genbank ICARDA (Syria), and were maintained at the Institute of Cytology and Genetics.

BAC-library screening and DNA isolation from individual BAC-clones

The 5BS-specific BAC-library of T. aestivum (43,776 BAC-clones, mean insert size 122 Kbp, 15-times coverage of 280 Mbp length chromosome arm) was obtained from the Institute of Experimental Botany (Olomouc, Czech Republic), kindly provided by Professor J. Doležel. A copy of the library is maintained at the Institute of Cytology and Genetics SB RAS at -80 °C. For BAC-library screening we ran PCR analysis of 2D BAC-pools with primers specific to 5S rDNA coding sequences (5SrDNA_F: 5′-GAGAGTAGTACTAGGATGGG-3′; 5SrDNA_R: 5′-GGAGTTCTGACGGGATCCGG-3′). PCR was performed in a 20 μl reaction mixture containing 0.5 μl of cell culture as a template, 0.25 pM of specific forward and reverse primers, 2 μl of PCR buffer (65 mM Tris-HCl, pH 8.9; 1.5 mM MgCl₂; 16 mM (NH₄)₂SO₄; 0.05% Tween 20), 0.2 mM of each dNTP and 1 unit of Taq DNA polymerase. After initial denaturation at 94 °C for 4 min, 35 cycles were run at 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s. PCR fragments were separated by electrophoresis on 1% agarose gels.

DNA of selected BAC clones was isolated using the NucleoSpin 96 Flash kit (Macherey-Nagel, Germany).

Sanger BAC-end and IonTorrent sequencing

The BAC-end sequences for selected BAC-clones were obtained with the universal M13 Reverse (5′-CAGGAAACAGCTATGAC -3′) and T7 forward (5′-TAATACGACTCACTATAGGG-3′) primers using the BigDye3.1 Terminator kit (Applied Biosystems, USA). Each 20 μl reaction contained ~200 μg of BAC-DNA, 1.5 μl of BigDye 3.1, 0.25 pM of specific forward, or reverse primer, 4 μl of5X buffer and deionized water. After preliminary denaturation at 95 °C for 5 min, 80 cycles were run at 95 °C for 30s, 55 °C for 15 s, and 60 °C for 4 min. The reaction products were precipitated using ethanol, and separated in a 3730XL DNA Analyzer (Perkin Elmer Cetus, USA).

BAC clones marked by 5S rDNA were included in the pool of 134 BAC-clones belonging to different locations of chromosome 5BS, and were collectively sequenced as one sample on an IonTorrent platform (Thermo Fisher Scientific). The sequencing and assembly were described in Nesterov et al. [23].

Pyrosequencing and assembly of BAC clones

The selected BAC-clones of the 5BS chromosome were sequenced into two pools. BAC pool 52 consisted of 6 clones including 5S rDNA-tagged clones TaaCsp5BS010O13 and TaaCsp5BS025F09. Ten overlapping clones of the 5BS chromosome (http://www.wheatgenome.org/), including BAC TaaCsp5BS096G09 with 5S rDNA sequences, were 89 pool. BAC pools were shotgun sequenced using a GS FLX/454 pyrosequencing platform (Roche). Sequence data of clone TaaCsp5BS096G09, as a part of ctg4 have been made available (https://urgi.versailles.inra.fr/gb2/gbrowse/wheat_phys_5BS_v1/, https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies).

For both pools, a library of random BAC-fragments with sizes ranging from 400 to 1000 bp was created. Also, for both pools a paired-end library of 6–10 kbp was constructed. The DNA sequencing was conducted according to Sequencing Method Manual GS FLX+ series with Titanium L+ Kit (Roche). The reads, belonging to E.coli DH10B and pIndigoBAC-5 vector were removed. The selected shotgun and paired-end reads were de novo assembled into contigs (for shotgun reads) and scaffolds (for paired-end reads) using the GS DeNovo Assembler V 2.9.

Identification and subsequent analysis of 5S rDNA-tagged sequences in selected BAC-clones

First, we searched for 454-contigs longer than 700 bp (shotgun reads) and scaffolds (paired-end reads) tagged by 5S rDNA, using: (1) BLASTn searches with sequences of complete 5S rDNA units (both coding and non-transcribed sequences); and (2) BLASTn searches with obtained BES (BAC-end sequences) [24]. Additionally, for verification and elongation of some contigs we used the IonTorrent BAC-clone sequencing data and PCR product sequences obtained with specific primers (Additional file 1: Table S1). The 6 selected 5S rDNA-containing contigs and scaffolds (hereafter named “5SrDNA_fragments”) were annotated by BLASTn search in TREP (http://botserv2.uzh.ch/kelldata/trep-db/index.html) and NCBI non-redundant nucleotide databases (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The nucleotide sequences of 5S rDNA fragments were deposited in GenBank under the accession numbers: MF467437 for pool_52 (5 unordered pieces) and MF467438 for pool_89.

Primers were developed with the Primer3 program [25]. To test for the degree of BAC-clone overlap, 3 Insertion Size Based Polymorphism (ISBP) primer pairs were designed and tested on the individual BAC-DNA templates (Additional file 2: Figure S1). ISBP exploits knowledge of the sequence flanking a TE to PCR-amplify a fragment spanning the junction between the TE and the flanking sequence [26].

To date the LTR-retrotransposon insertion events for the autonomous transposable elements, we analyzed the nucleotide divergence rate between two LTRs in cases when both LTRs were present in the element’s structure. To determine the LTR boundaries, each element was compared with itself using Blast2seq (https://blast.ncbi.nlm.nih.gov/Blast.cgi). In addition, the presence of the characteristic motifs, 5′-TG-3′ and 5′-CA-3′, at the beginning and end of each LTR, respectively, was taken into account. Each pair of LTRs was aligned using the ClustalW algorithm within the MEGA4 program [27]. Sequence divergence was calculated using the Kimura two-parameter method [28] with complete deletion option. To convert this term into the insertion date, we used the following equation: T = D/2r, where T is the time elapsed since the insertion; D, the estimated LTR divergence; and r, the substitution rate per site per year [29]. We applied a substitution rate of 1.3 × 10⁻⁸ mutations per site per year for the plant LTR retrotransposons [30].

5S rDNA sequence alignment and cluster analysis

A search of 5S rDNA sequences was performed in contigs and scaffolds and also in primary reads employing BLASTn search with 5S rDNA coding sequence as a query. In order to find spacer sequences between 5S rDNA coding regions we identified all reads containing two or more copies of 5S rDNA and extracted the sequences located between 5S copies. The coding and spacer sequences were analyzed separately. The alignment of 5S rDNA spacer sequences was done by MUSCLE program [31] included into MEGA 4 software. The calculation of 5S rDNA coding sequence divergence was done using the pairwise-deletion option and Kimura-2-parameter model. The cluster analysis of 5S rDNA coding sequences and spacers was performed with the CD-HIT program [32]. A neighbor-joining phylogenetic tree was obtained with the MEGA4 program using the pairwise deletion option with bootstrap replicates of 500.

Secondary 5S structure analysis

Secondary structure modelling was carried out using an online tool at the RNAfold web server (http://rna.tbi.univie.ac.at/). The secondary structures were based on minimum free energy (MFE) calculations using a loop-based energy model and the dynamic programming algorithm introduced by Zuker and Stiegler [33]. For both groups of 5S rDNA (pool_52 and pool_89 representing short and long units, respectively) a consensus model of secondary structures was constructed. The program setting was as follows: isolated nucleotides were avoided; vote for dangling energies on both sides of a helix in any case.

Fluorescence in situ hybridization (FISH)

Metaphase chromosome preparation, FISH and chromosome identification were performed according to Salina et al. [34] with minor modifications. The total number of analyzed metaphases from individual plants for each probe was 15–30. For FISH, we used the PCR-amplified sequences Short5S and Long5S (94 bp and 131 bp, respectively). The probes were amplified from BAC-clones TaaCsp5BS010O13 and TaaCsp5BS096G09 by PCR with specific primers Short5S-F (5′-GCGTGCACTGGTGCGGTTGAG-3′) and Short5S–R (5′- GACGATTGCACATTGCTTTGGC-3′); Long5S-F (5′-GGAAAAAACTCGTGTTGCTGC -3′) and Long5S–R (5′-CTCACTACCATTACAACCGTTC -3′) using the following program: 35 cycles at 94 °C for 45 s, at 55 °C for 45 s, at 72 °C for 30 s. The primers were designed to the representative spacer sequences of Short5S and Long5S types (determined by cluster analysis with a threshold of 98% identity), and obtained single band PCR product for each type was verified by sequencing. The mean overall sequence divergence calculated for presumable amplification regions of Short5S and Long5S 5S rDNA sequences (obtained from 454 pyrosequencing) corresponded to 0.8–0.9% for both probes. The divergence level between Short5S and Long5S probes at their complete length was 45%, wherein the small region of 62 bp has 84% identity between probes and there was no identity for the remaining sequence.

The PCR-derived probes were labeled with biotin or digoxigenin. Biotinylated probes were detected with fluorescein avidin D (Vector Laboratories, United States). The hybridization signal was enhanced using fluorescein anti-avidin (Vector Laboratories, United States). The digoxigenin-labeled probes were detected with antibodies to anti-digoxigenin-rhodamine (Fab fragments, Sigma-Aldrich, United States). The preparations were embedded in Vectashield mounting medium (Vector Laboratories), containing 0.5 μg/ml DAPI (4′, 6-diamidino-2-phenylindole, Sigma-Aldrich, United States) for chromosome staining. The chromosomes were examined with an Axioskop 2 Plus (Zeiss) microscope and recorded with a VC-44 (PCO) CCD camera.

To identify chromosomes carrying signals, we used the probes pSc119.2 [35] and pAs1 [36].

The work was performed at the Collective Center for Microscopic Analysis of Biological Objects (SB RAS, Novosibirsk).

Identification, sequencing and assembly of 5S rDNA tagged BAC-clones

The 5S rRNA gene number in BAC pools was approximately calculated. For pool 52 the 5S rDNA sequence coverage was established as 3600, whereas the contig coverage by 454-reads was 42, thereby the 5S rDNA copy number was assessed as 86. The 5S rDNA sequence coverage for pool 89 was 1500, and the contig coverage 24, which gives the 5S rDNA copy number as 63.

In order to increase the assembly quality for pool 52, paired end 454-sequencing was performed (Additional file 3: Table S2). Consequently, the assembled reads were arranged into 11 scaffolds with lengths ranging from 2159 to 164,054 bp. It is noteworthy that the lengths of the individual scaffolds are comparable with the BAC-clone length.

The data of shotgun and paired-end pyrosequencing were used for analysis of 5S rDNA-tagged genomic fragments from 5BS chromosome of the bread wheat (Additional file 1: Table S1).

Structural organization of 5S rDNA-tagged sequences of chromosome 5BS

For pool 52, five 5S rDNA fragments were identified and just one 5S rDNA fragment was found for pool 89. The 5S rDNA fragments were annotated using BLASTn searches against the Triticeae transposable element database (TREP) and NCBI nucleotide database (Table 2). All 5S rDNA fragments consist of a combination of transposable elements (TE) and 5S rDNA. No unique sequences representing potential genes or pseudogenes were found. The TEs in 5S rDNA fragments appear as single elements (or their component parts) interspersed within the 5S rDNA, or as multiple, nested TE insertions (predominantly LTR-retrotransposons). The 5S rDNA sequences, mainly located at the ends of 5S rDNA fragments, represent partial sequences from 92 bp up to complete units. In 5S rDNA fragment 3, the 5S rDNA was introduced within the fragment sequence (Fig. 1). In just one case (5S rDNA Fragment 4) the LTR-retrotransposon Fatima was inserted into the coding sequence of 5S rDNA, while in all other cases the TEs were inserted into spacer sequences.

The overlap between BAC-clones TaaCsp5BS010O13 and TaaCsp5BS025F09 from pool 52 was tested by PCR, with specific primers designed to TE insertions (ISBP-markers) (Additional file 2: Figure S1). We showed that BAC clone TaaCsp5BS025F09 completely overlaps with TaaCsp5BS010O13, which is apparently longer. Since pools 52 and 89 consist of different groups of overlapping BAC clones, we can conclude that pool 52 and pool 89 are attributed to different genomic locations on chromosome 5BS of T. aestivum.

The analysis of pool 52 and pool 89 sequences implies their different structural organization. In pool 52, 5S rDNA sequences were interrupted by multiple TE insertions, whereas in pool 89 only one TE insertion was observed on one of the flanking regions (Table 2). The repetitive 5S rDNA units of pool 89 were apparently merged in the assembly process, suggesting a single tandem cluster of 5S rRNA genes adjoining the Nusif LTR-retrotranposon. This is confirmed by the correspondence of pool 89 to the 5BS Illumina pseudomolecule (IWGSC RefSeq v1.0, http://www.wheatgenome.org/). Simultaneously, no correspondence between BAC-clones TaaCsp5BS010O13 and TaaCsp5BS025F09 from pool 52 and the Illumina data was found.

5S rDNA coding and spacer sequence analysis

The 5S rDNA coding and spacer sequences were extracted from the shotgun 454-reads of pool 52 and pool 89. The 5S rDNA sequences were assembled into two files for NTS, and the coding sequences for each pool analyzed separately.

A total of 1511 and 615 copies of 5S rDNA sequence were identified in reads obtained for BAC pools 52 and 89, respectively. For analysis, we removed 5S rDNA sequences which occurred less than 5 times since these are likely to be a result of sequencing errors. Cluster analysis with 99% sequence identity threshold yielded 30 types of complete 5S rDNA coding sequences for pool_52 and 21 type for pool_89. The sequences were aligned (Additional file 4) and mean overall sequence divergence was calculated: for pool 52, the value was 2.8% and for pool 89–1.1%. Thus, the coding sequences in pool 52 showed more sequence heterogeneity than in pool 89.

To evaluate the functionality of the rRNA genes at the different loci of the 5BS chromosome, we used a special program that predicted a consensus secondary structure for each pool of RNA sequences and gave an accompanying estimate of its thermodynamic stability (see Methods). In fact, both pools produced identical molecular shapes, consisting of three functional domains (Fig. 2). The thermodynamic stability of the 5S RNA for pool 89 was slightly higher than for pool 52, although in both cases it is quite high in comparison with what is typical for functional genes (~ 50 kcal / mol) [37]. Nevertheless, both pools differed significantly in the number of incompatible pairs (15 for pool 52 and 7 for pool 89). In other words, there are approximately twice as many gene sequences in pool 52 producing RNA structures with unpaired nucleotides after folding.

Sequences of pool 52 and pool 89 differ in the structure of NTS types. In pool 89 we found 27 complete NTS sequences. The complete sequences, with mean length of 368 bp and low sequence heterogeneity, as well as partial spacer sequences, were subjected to cluster analysis (Additional file 5). Three clusters were obtained at 98% sequenced identity threshold. For the alignment we took representative sequences of these three clusters (Additional file 6). All pool 89 NTS sequences were attributed to LongS1 type (as designated by Baum and Bailey [10]).

In pool 52 we found 144 complete NTS sequences which were divided into two major clusters of 88 and 10 sequences (under a threshold of 98% identity) with a mean spacer length of 295 bp, and 39 minor clusters containing from one to four sequences of 207 to 390 bp (Additional file 5). The 22 NTS originating from minor clusters had large insertions from 30 to 95 bp, that cannot be sequencing errors.

Construction of the neighbor-joining tree showed that most NTS sequences from pool 52 (133) correspond to the ShortA2 type (in agreement with the classification of Baum and Bailey [10]), whereas 8 minor NTS clusters are closer to ShortA1 and ShortG1 units. All 27 NTS sequences from pool 89 correspond to the LongS1 type (Fig. 3).

The overall mean sequence divergence with the pairwise deletion option between the sequences of tree branches is 21% between ShortA2 and ShortA1, 24% between ShortA1 and LongS1, and 29% between ShortA2 and LongS1. Within each of tree branches the mean overall divergence values were: 8% for ShortA2, 11% for ShortA1 and 12% for LongS1. The alignment of representative pool 52 cluster sequences derived from both ShortA2 and ShortA1 branches of the tree (Additional file 7) showed that the major representative cluster of 88 spacer sequences was completely identical to the ShortA2 type; while the two minor clusters from ShortA1 showed 19% divergence in comparison with ShortA1 and ShortG1.

To check the affinity of pool 52 and pool 89 NTS sequences to the presumed wheat progenitors, T. urartu, Ae speltoides and Ae. tauschii, we first performed the GenBank BLAST search with representative sequences 89_cluster_0 and 52_cluster_29 as queries. Among the 5S rDNA sequences of diploids, the Ae. speltoides sequence AY316211.1 is a single highly homologous sequence (99% of identity) to the ShortA2 type, whereas the remaining 29 Ae. speltoides sequences from GenBank showed only the 75–76% of identity to this type. The similar level of homology to ShortA2 type is characteristic of 5S rDNA from T. urartu, T. monococcum and Ae. tauschii.

The LongS1 type (pool 89) has a greater affinity to the corresponding sequences of 5S rDNA from diploids: for Ae. speltoides 5S rDNA there are 18 GenBank sequences with identity from 91 to 99%, for Ae. tauschii - 19 sequences with 82–87% and 30 sequences with 82–84% homology from T. urartu and T. monococcum. It should be noted that both LongS1 and ShortA2 show high identity to numerous 5S rDNA sequences from tetraploid T. turgidum. We added the best-matching GenBank 5S rDNA sequences of T. urartu, T. monococcum, Ae. speltoides, Ae. tauschii and T. turgidum to the phylogenetic tree (Fig. 3).

FISH of 5S rDNA NTS probes

For FISH, we used two probes representing the major type 5S rDNA spacer sequences for pool 52 and pool 89. Previously, the complete 5S rDNA unit sequences were usually used as a FISH probe, thus Mukai and coauthors [7] used the rye 5S rDNA probe pScT7, and Badaeva et al. [15, 38] used the clone pTa794 [39]. Baum and coauthors [11] assessed whether the different Triticeae 5S rDNA units that were assigned to haplomes could be used as a FISH- probe for tracing the origin of polyploid wheats from their diploid progenitors, or chromosomal remodeling during speciation. But interpretation of results in these cases is rather difficult due to a high level of cross-homology between different units, especially at their coding regions and 5′- and 3′- ends of NTS. Therefore, for greatest specificity in the FISH signal, we designed specific primers allowing the amplification of a short sequence within NTS to be used as a probe. The probe Short5S is 94 bp long from a major cluster (52_cluster_29) and the probe Long5S is - 131 bp, from representative copies of 5S rDNA of pool 89 (Additional file 6, Additional file 7).

Figure 4 represents the distribution of both probes on T. aestivum chromosomes. Probes Short5S and Long5S had different chromosomal locations: 10 hybridization sites on the distal parts of chromosomes 5BS, 1BS, 5DS, 5AS, 1DS were revealed for Long5S and the strongest signal was at 5BS (Table 3). Short5S showed the strongest hybridization signals on 1BS and weak signals on 1DS and 5BS. Tetraploid wheat T. timopheevii also had strong Long5S and weak Short5S signals on chromosome 5GS, and weak Long5S signal on chromosome 5AS. FISH of Short5S and Long5S to chromosomes of diploid species T. urartu, T. monococcum, and Ae. speltoides revealed no sites for Short5S hybridization. Ae. tauschii have both sites of 5S rDNA units, Short5S on 1DS and Long5S on 5DS. It should be noted that in polyploid wheats, Short5S units have less sites of localization than Long5S, but in most cases localization is close to Long5S on chromosomes. Copy number variation or level of 5S rDNA unit divergence may be the cause of differences in signal intensities between Short5S and Long5S in one chromosome loci. As shown in this study for chromosome 5B, the differences in signal level of Short5S and Long5S are connected with a high level of divergence of Short5S compared to Long5S (Fig. 2), despite approximately the same number of copies of these genes in the 5BS loci (84 and 75, Table 1).

Species	Long5S	Short5S
T. aestivum (AABBDD)	5BS > 1BS > 5DS > 5AS > 1DS	1BS > 1DS > 5BS (very weak)
T. urartu (AuAu)	5AS	–
T. monococcum (AmAm)	5AS	–
Ae. speltoides (SS)	5SS	–
Ae. tauschii (DD)	5DS	1DS
T. timophevii (AAGG)	5GS > 5AS	5GS – very weak