Material of polyploid dogroses was sampled in wild populations in Germany and the Czechia (Supplementary Table S1). Diploid species and tetraploid species with regular meiosis were obtained from various Botanical Gardens (Supplementary Table S1). In addition, we retrieved sequence information from published work stored in the ENA database for bioinformatics analyses (Supplementary Table S1).

Total genomic DNA of Rosa canina was extracted from fresh leaves applying the standard protocol (Rogers and Bendich, 1985). The 5S rDNA repeats were amplified using the primers pr5S-14 and pr5S-15 (Tynkevich and Volkov, 2014b) with the 5’-extensions containing restriction endonuclease NotI recognition site. The PCR products were treated by NotI, ligated into the Eco52I recognition site of the pLitmus 38 plasmid, and used for transformation of Escherichia coli XL_blue line by electroporation method. Selected recombinant clones were sequenced using a BigDye Terminator Cycle Sequencing Kit (Thermofisher Scientific, United States). Clones containing inserts of A and B variants of 5S rDNA were identified by sequence analysis. The sequences were submitted to GenBank under the accession numbers MW349696 and MW349697.

The inserts of cloned 5S rDNA sequences contained genic regions and IGS. In order to increase the specificity of probe hybridization, we amplified the IGS sequences using specific primers annealing to 5S_A and 5S_B variants. The oligonucleotide primers’ sequences for the 5S_A IGS were as follows: A_for: 5′-CCTCTTTTTTCTGTTTCGGT-3′; A_rev: 5′-ATAAACTCCATTCGCTCAG-3′. Primers for the 5S_B variants were: B_for: 5′-ACCCCTCTTTTTGCCTTT-3′; B_rev: 5′-GCTTCGTCTCACTCCTCT-3’. The 25 μl PCR reaction contained 0.1 ng of plasmid DNA as the template, 4 pmol of each primer, 2.4 nmol of each dNTP, and 0.4 units of Kapa Taq DNA polymerase I (Kapa Biosystems). Cycling conditions were as follows: initial denaturation step (94°C, 180 s); 35 cycles (94°C, 20 s; 57°C, 30 s; and 72°C, 30 s). The length of amplified products was 373 nt for the A variant and 394 nt for the B variant. Purified PCR products were labeled by fluorescent dyes (as below) and used in FISH.

For bioinformatic analyses, the whole genomic sequencing data for 19 Rosa accessions were used (Supplementary Table S2). The genome proportion of 5S rDNA families was determined using the total Illumina reads trimmed for quality (Phred score ≥ 30 over ≥ 95% read length). Trimmed reads (typically > 7 million) were mapped to corresponding 5S_A and 5S_B reference sequences (IGS subregion between the primers, Figure 1) using the following parameters: insertion and deletion costs_3, lengths fraction_0.5, similarity fraction_0.8, and deletion cost_2 (Qiagen, Germany). The distribution of SNPs across the 5S rDNA sequences was recorded when the distribution exceeded a threshold of at least 20 identical SNPs over at least 200 reads that covered the variant position and occurred at ≥ 10% frequency. For the alignment, minimal sequence length coverage was 50% and minimum sequence similarity was 90%. For the more distantly related species, Rosa spinosissima and R. persica, similarity threshold parameter was decreased to 80% (for the 5S_B variant). The genome abundance and copy number was calculated from genome proportions according to the formula stated in (Lunerova et al., 2020).

For phylogenetic reconstructions, the consensus 5S_A and 5S_B rDNA sequences were extracted from mapped reads using CLC genomic workbench. Additionally, we added the partial sequence of 5S ribosomal RNA genes from Acaena latebrosa (EU931698.1) and Cliffortia curvifolia (EU931716.1). Paired 250 bp Illumina reads of Geum urbanum (SRA accession ERR2187925) were mapped in a first round to the C. curvifolia (EU931716.1) sequence. Mapping was done with Bowtie2 (Langmead and Salzberg, 2012) implemented in Geneious^® 10.0.9¹ with the lowest sensitivity pre-set. This resulted in two reads out of 22.8 million hitting to a 39 bp conserved region. The two reads were aligned and served in a second round with same parameters as mapping scaffold of 334 bp length. In the second mapping, 50 reads out of 22.8 million were assembled and its consensus was added to the alignment including all rose sequences and the two species from the Acaena clade. Alignments were done with MAFFT v7.450 algorithm (Katoh and Standley, 2013) implemented in Geneious using default parameters. Model test implemented in MEGA X v. 10.1.8 (Kumar et al., 2018) revealed the Tamura–parameter substitution model with invariant sites as most appropriate for the data based on Akaike information criterion (Tamura, 1992). Based on this model, we computed a maximum-likelihood tree in MEGA X whose branch support was evaluated by 1000 bootstrap replicates. Rooted with G. urbanum, this tree was used for the dating approach conducted with MEGA X. Therefore, we used two calibration points along the tree taken from the respective fossils given in Xiang et al. (2017) Rosa germerensis, 48.6 Mya (Edelman, 1975) and Acaena sp., 37.2 Mya (Zetter et al., 1999). A timetree inferred by applying the RelTime method (Tamura et al., 2018) was computed using two fixed calibration constraints. All positions containing gaps and missing data were eliminated (complete deletion option). A neighbor joining tree was constructed with the Seaview software (Gouy et al., 2010).

The fastq reads were initially filtered for quality and trimmed to uniform length using the pre-processing and QC tools in RepeatExplorer2 (Novak et al., 2013). Read length ranged between 100 and 150 bp, depending on sequencing library and the Illumina sequencing platform. After the fastq > fasta conversion and trimming to uniform length, reads were analyzed with the RepeatExplorer2 clustering program using default parameters. We used 1 million paired-end reads, or 1 million single-end reads as inputs for RepeatExplorer2 clustering. This bioinformatic pipeline runs a graph-based clustering algorithm (Novak et al., 2013) that assembles groups of frequently overlapping reads into clusters of reads, representing a repetitive element or part of a repetitive element with a higher order genome structure. The similarity and structure-based repeat identification tools in RepeatExplorer2 aid in identification of the repeats. RepeatExplorer2 uses a BLAST threshold of 90% similarity across 55% of the read to assign reads to clusters (minimum overlap = 55, cluster threshold = 0.01%, minimum overlap for assembly = 40), and the clusters are identified based on the principle of maximum modularity. We also used the SeqGrapheR (Novak et al., 2010) software in virtual space of Ubuntu 18.04 to visualize the specific reads corresponding to the 5S_A and 5S_B variants.

For slide preparations, we used young anthers from flower buds of about 0.5 cm in length, which were harvested during spring 2019. Male meiosis was studied at prophase I (diplotene/diakinesis) where the bivalents and univalents could be easily distinguished from each other. Fresh flower buds were fixed using Carnoy solution (acetone:acetic acid, 2:1 or ethanol:acetic acid, 3:1 in some cases), and stored in 70% ethanol at −20°C. Before slide preparation, anthers were pre-treated by 0.5% PVP and 2% Triton-X100 (Sigma–Aldrich, United States) for 2–5 min followed by enzyme digestion overnight at 10°C in 1% cellulase, 0.2% pectolyase Y23, 0.5% hemi-cellulase, and 0.5%, macerozym R10 (Sigma–Aldrich, United States; Duchefa Biochemie, Netherlands) dissolved in citric buffer (0.04 M citric acid and 0.06 M sodium citrate). FISH followed the procedures described in Herklotz et al. (2018). Anthers were separated and squashed on slides in a drop of 70% acetic acid and fixed in liquid nitrogen.

For FISH, we used two probes derived from the 5S_A and 5S_B clones of the IGS region, respectively, and in addition, an 18S rDNA probe that was a 1.7-kb fragment of the 18S rRNA gene of Solanum lycopersicum (GenBank # X51576.1). The 5S rDNA genic region originated from Artemisia tridentata S4 clone, GenBank # JX101915.1. The probes were labeled by nick translation using Spectrum green dUTPs (Abbott, United States) for 5S_A rDNA, and Cy3-dUTPs (Roche, Switzerland) for 5S_B and 18S rDNA; 5S rDNA was labeled by Atto647N (Jena Bioscience, Germany). Slide preparation and hybridization followed standard protocols (Schwarzacher and Heslop-Harrison, 2000). Chromosomes were counterstained with 1 μg ml^–1 DAPI (4’, 6-diamidine-2’-phenylindole dihydrochloride) diluted in mounting medium for fluorescence (Vectashield, Vector Laboratories, United Kingdom). The slides were scanned using epifluorescence microscopes (Olympus Provis AX70, with cold cube camera, Metasystems, Germany). Imaging software was ISIS (MetaSystems, Germany), and images were optimized for contrast and brightness with Adobe Photoshop CS6 and PS2020.

Two 5S rDNA clones (5S_A and 5S_B) were isolated from R. canina IGS. Sequence analysis revealed some conserved regulatory elements: Box-C within the coding region, the TATA box at –29 (both clones), and T-rich terminators downstream of the coding sequences (Figure 1A). Box-A could not be unambiguously determined due to primer overlap. The 5S_B clone had a long (20 nt) T-tract which appears to be missing or was much shorter in clone 5S_A. By analogy with other 5S rRNA transcripts, the putative transcription started at the first G within the GGG motif following the C at −1 (Tynkevich and Volkov, 2014b). Intragenomic homogeneity was high, and no significant SNPs were revealed in mapping experiments (not shown). Pairwise alignment revealed conserved coding regions, while most of the IGS was dissimilar between both sequences (Figure 1B). We took advantage of considerable sequence divergence between both clones and amplified the locus-specific IGS subregions from plasmids. The resulting 373 bp (5S_A family) and 394 bp (5S_B family) PCR products were subsequently used in FISH.

To determine the abundance of individual 5S rRNA gene families in Rosa genomes, we used available genomic resources (Table 1). The genome proportion of the 5S_A family was in average twice of that of the 5S_B family (Supplementary Tables S2, S3). The copy number per somatic cell (2C) ranged from 80–8000 (5S_A family) and 0–2400 (5S_B family). Both families appeared to be equally homogeneous containing a relatively low number of SNPs consistent with our previous findings obtained by comparison of individual 5S rDNA clones (Tynkevich and Volkov, 2014a, b). The contribution of each family to total 5S rDNA was expressed for each species by pie charts and visualized in a phylogenetic context (Figure 2). The diploid species from sect. Synstylae, all polyploid species, and R. persica (subg. Hulthemia) carried both families. Rosa laevigata contained the 5S_A family in low copy (c. 80 copies/2C), while its 5S rDNA was dominated by the 5S_B family (980 copies/2C) (Supplementary Table S3). Species from sect. Rosa and R. minutifolia (subg. Hesperhodos) carried the 5S_A rDNA family only. Blast searches failed to reveal significant hits of 5S_A and 5S_B sequences in genomic reads from Prunus, Rubus, Fragaria, Cliffortia, Acaena, and Sanguisorba (all Rosaceae) even at relaxed (e = 0.1) stringencies (not shown).

Cloning experiments cannot address the question about the distribution of gene families in the genome. Thus, in order to determine the number and genomic representation of individual 5S rRNA gene families, we applied clustering analysis (Figure 3). The cluster graph shapes provide information about the number and type of 5S gene families revealing potential hybridization and introgression (Garcia et al., 2020). It visualizes divergent IGS families as loops emanating from the bridge region which contains reads derived from a conserved coding region. Any loop can be considered as a separate gene family. The subregions in the graph can be annotated based on the read alignment against the 5S_A and 5S_B reference clones. Visual exploration indicates that there is no other 5S family amplified except of 5S_A and 5S_B types. The cluster graphs obtained from different Rosa genomes were categorized based on their structure into three groups (Figure 3). Group 1 comprised a single species R. laevigata (sect. Laevigatae) with predominant 5S_B type family representation. Group 2 comprised R. persica (subg. Hulthemia) and the majority of diploid species (sect. Synstylae and Indicae) showing a typical two-loop structure representing relatively balanced ratios of both families. Group 3 contained the diploid species R. minutifolia (subg. Hesperhodos) and species of sect. Rosa harboring a single 5S rDNA family (A). All polyploid species (both, those with regular meiosis and those with Canina meiosis) showed a Group 2 profile indicating the presence of both A and B families (Supplementary Figure S1). In sum, quantitative relationships between both 5S rDNA families were confirmed. Moreover, cluster analysis showed that the maximum number of 5S rDNA families in the Rosa genomes is always two, irrespective of the ploidy level.

To determine the phylogenetic relationships between 5S rDNA families, we computed phylogenies based on aligned 5S rDNA consensus sequences (obtained from mapping experiments. Both the maximum-likelihood (Figure 4 and Supplementary Figure S2) and neighbor joining (Supplementary Figure S3) trees separate the A and B 5S rDNA families clearly into two well-supported clades (A and B). Both, clades A and B contained diploid and polyploid species. Except of sect. Rosa, whose members clustered exclusively within the A clade, members of other sections, including Synstylae, Indicae, Laevigatae, Pimpinellifoliae, and Caninae, partitioned their 5S rDNA between both clades. Sequences from R. persica (subg. Hulthemia) were consistently positioned on early diverging nodes at both subclades. The major 5S rDNA family of R. laevigata (sect. Laevigatae) branched off at a rather basal position in clade B. The 5S_B family of 4x R. spinosissima positioned as sister to R. persica. Five 5S rDNAs of 5x species from sect. Caninae clustered together in both clades with negligible resolution between species (Supplementary Figures S2, S4). Out of the diploids, 5S sequences of R. arvensis (sect. Synstylae, B clade) and R. pendulina in (sect. Rosa, A clade) were most closely related to those of the respective Caninae branches. To gauge the length of time these 5S rDNA variants have existed in the Rosa genomes, we used two calibration points (48.6 myrs for R. germerensis and 37.2 myrs for Acaena sp.). We estimated that a common ancestor of both A and B families lived about 32 myrs ago (Figure 4) relatively long before separation of modern clades.

Fluorescence in situ hybridization was conducted to visualize the position and number of 5S rDNA variants on chromosomes in several diploid and polyploid species. The diploids included representatives of subg. Rosa sect. Synstylae (R. arvensis and R. multiflora), sect. Rosa (R. rugosa, R. majalis, and R. nitida), and subg. Hulthemia (R. persica). Meiotic chromosomes from anthers (Figure 5A) were hybridized with rDNA probes derived from the 5S rDNA genic (Figure 5B, shown in white), 5S_B (red), and 5S_A (green) IGS subregions (Figure 5C). Additionally, the same chromosome spreads were re-hybridized with the 18S rDNA probe (shown in cyan, Figure 5D). In R. arvensis and R. multiflora, each 5S_A and 5S_B probe hybridized to one bivalent (one pair of chromosomes). The 5S_B probe was always co-localized on a chromosome bearing also the 18S rDNA signal. In R. rugosa and R. majalis, the 5S_A probe hybridized to a single bivalent, while we did not detect any hybridization signals with the 5S_B probe in accordance to genomic analyses (Figures 2, 3). In both species, the 18S and 5S rDNA loci were separate. However, in R. nitida, the 5S_A probe hybridized to one pair of chromosomes (mitotic metaphase from root tips, Supplementary Figure S5) which carried 18S rDNA signal (NOR). In R. persica, both variant-specific 5S rDNA probes hybridized to a single bivalent each, and these bivalents carried also 18S rDNA sites. The number and position of rDNA loci on chromosomes are summarized in Supplementary Table S4 and are diagrammatically depicted by ideograms (Figure 6).

We further analyzed meiotic (Figure 7) and mitotic (Supplementary Figure S5) chromosomes in polyploid species. As expected meiotic chromosomes of four 5x dogrose species (sect. Caninae) were represented by seven bivalents (pairs of chromosomes) and 21 univalents (Figure 7). In R. canina and R. corymbifera (subsect. Caninae), the 5S_A probe hybridized to one bivalent and three sites on univalent chromosomes. The 5S_B probe hybridized to one bivalent carrying the 18S (NOR) signal and two sites on univalents. The 5S_A and 5S_B signals occurred on different chromosomes except of one univalent chromosome in R. canina where both probes were co-localized. In R. inodora (subsect. Rubigineae), the 5S_A probe hybridized to one bivalent and three univalent chromosomes. The 5S_B probe hybridized to two univalent chromosomes carrying the 18S rDNA signal. Rosa rubiginosa (subsect. Rubigineae) showed a similar distribution of signals like R. inodora except that only one out of two 5S_B univalent chromosomes co-localized with the 18S signal. In addition, there were only two 5S_A sites on univalents. Collectively, these observations indicate that the number of rDNA sites, their chromosome position, and their meiotic behavior differ between subsections Caninae and Rubigineae. FISH on mitotic chromosomes from 4x R. spinosissima is shown in Supplementary Figure S5. In this species, the 5S_B probe hybridized to a chromosome pair which also carried the 18S rDNA signal (NOR) (Supplementary Figure S5). Two other 5S_B signals were colocalized (but did not overlap) with that of 5S_A on non-NOR chromosomes. 18S rDNA and 5S_A signals were localized on two separate chromosomes. Results are summarized in Figure 6 and Supplementary Table S4.

The IGSs of rRNA genes are rapidly evolving sequences, and it is common to find variation even between closely related species. It was therefore striking to observe that the genus Rosa is dominated essentially by only two 5S rDNA families and that no other family was amplified in any of the species analyzed here. Both families occupy different chromosome loci: the 5S_B family was always co-localized with NOR (35S rDNA), while the 5S_A family was mostly but not exclusively (see below) separate (Figure 6). Moreover, R. persica (subg. Hulthemia) amplified both families at similar ratio in its genome (Figures 2, 3). In contrast to all other diploid species (Ma et al., 1997), R. persica is also exceptional in possessing two NORs instead of one per haploid chromosome set. Since R. persica was consistently identified as the earliest divergent lineage in most phylogenies (Fougere-Danezan et al., 2015; Debray et al., 2019), we presume that the configuration with NORs co-localizing with distinct 5S rDNA families (Figure 6) is an ancient condition, while the NOR chromosome without 5S rDNA locus is derived. This assumption is supported by the following observations: First, the Rosa and allies clade contained the 5S_A family which was co-localized with 18S rDNA locus in R. nitida but not in R. majalis and R. rugosa. Second, all members of Synstylae, Indicae, and Pimpinellifoliae contained two 5S rDNA families albeit at differing ratios. For example, R. multiflora (sect. Synstylae) had prevalent 5S_B family, while the 5S_A family dominated in R. chinensis, a member of the closely related sect. Indicae. Similarly, in R. lucieae ([R. wichurana]; sect. Synstylae), both families are likely to be represented by two loci out of which one co-localized with 18S rDNA on the same chromosome (Kirov et al., 2016). Third, R. laevigata (sect. Laevigatae) was dominated by the 5S_B family (Figures 2, 3). This species is sister to the remaining species of the Synstylae and allies clade (Fougere-Danezan et al., 2015; Debray et al., 2019).

Of note, R. nitida differed from R. majalis and R. rugosa in having the 5S_A family co-localized with NOR. Traditionally, R. nitida has been classified into a separate section called Carolinae (Crépin, 1889). However, more recent taxonomies based on molecular markers failed to support this distinction and all three species are now placed within the Rosa and allies clade (Wissemann and Ritz, 2005; Joly and Bruneau, 2006). Interestingly, members of sect. Rosa tend to have much smaller loci of the centromeric satellite repeat CANR4 compared to other species of the genus (Lunerova et al., 2020). However, R. nitida is exceptional in having large abundance of the CANR4 satellite (10 out 14 chromosomes carried strong FISH signals, not shown). These features suggest chromosomal rearrangements accompanying speciation events in sect. Rosa although the basic chromosome number (x = 7) remained unchanged.

Neither A nor B type sequences were found in 5S rDNA of the genera Prunus, Rubus, Fragaria, Cliffortia, Acaena, and Sanguisorba (all Rosaceae). These observations suggest that both 5S rDNA families have their origin in the early evolution of the genus Rosa because the common ancestor of both families was dated at app. 32 myrs ago. Furthermore, molecular dating of diversification within both 5S rDNA clades was dated to app. 5–9 myrs each (Figure 4) matching at least roughly the origin of Synstylae and allies and Rosa and allies (Fougere-Danezan et al., 2015).

Despite considerable interest, the composition of pentaploid (2n = 5x = 35) dogrose genome remains enigmatic (Wissemann, 1999; Nybom et al., 2004; Ritz et al., 2005; Herklotz et al., 2018; Lunerova et al., 2020). Previous studies based on microsatellite markers revealed genetic distinction between bivalent- and univalent-forming chromosomes (Nybom et al., 2006). The analysis of 35S rDNA markers confirmed these assumptions revealing two highly divergent ITS types (named Canina type and Rubiginosa type) present at variable ratios in subsections Caninae and Rubigineae, respectively (Ritz et al., 2005; Kovařik et al., 2008). Here we show that the Canina type ITS is co-localized with 5S_B rDNA locus (NOR), while the Rubiginosa ITS type is not. This Canina type of configuration [equivalent to marker chromosome 1 (Lim et al., 2005) or chromosome 7 (Kirov et al., 2016)] is typical for the bivalent-forming chromosomes in R. canina and R. corymbifera, both from subsection Caninae (Figure 6). In contrast, the Canina-type configuration of both rDNAs is found on univalent chromosomes in R. rubiginosa and R. inodora (subsection Rubigineae). The bivalent-forming chromosomes in this subsection carry the 5S_A family on a non-NOR chromosome. These data support the hypothesis that dogroses from both subsections arose by reciprocal hybridization of closely related species (Bruneau et al., 2007; Herklotz and Ritz, 2017). A very similar composition of dogrose genomes is also supported by the shallow nodes among dogroses in the 5S rDNA phylogeny (Figure 4). However, within the subsections, we detected some variation in number and distribution of loci especially between the univalent sets. For example, R. inodora carried two univalent-forming chromosomes with Canina type configuration, while R. rubiginosa contained only one. Additionally, one of the R. canina univalent chromosomes carried both 5S_A and 5S_B chromosomes co-localized, while this chromosome was not found in the closely related R. corymbifera. The variation between univalent chromosomes is consistent with increased diversity of microsatellite markers on univalent genomes (Nybom et al., 2004) and could be either related to divergence of genome donors and/or to partial degeneration of univalent chromosomes due to their exclusion from recombination in meiosis. Although 5S rDNA pseuodogenes seem to be present in R. rugosa (Tynkevich and Volkov, 2014b), there is no indication for extensive 5S rDNA pseuodogenization in R. canina (Lim et al., 2005) and in other dogroses (this work).

Although it might be preliminary to trace potential genome donors of allopolyploid dogroses, it is notable that the tetraploid species with a regular meiosis, R. gallica (sect. Gallicanae), tends to cluster with dogroses. It is therefore possible that pentaploid dogroses actually arose by pollination of an unreduced non-dogrose tetraploid egg cell by a reduced male gamete from a diploid donor. In support, R. gallica ITS types were occasionally found in dogroses (Herklotz et al., 2018) and R. gallica belonged to a clade together with dogroses in plastid phylogenies (Fougere-Danezan et al., 2015). Moreover, R. gallica shares the distinct morphological feature of partially pinnate sepals with dogroses, which is absent in the remaining species of the genus. However, there are also tetraploid cytotypes within dogroses (e.g., R. villosa L., R. canina) forming triploid egg cells and haploid sperm cells, which would not fit in the proposed scenario so far but might have arisen by another combination of partially reduced gametes. The also occurring higher ploidy (6x) levels which are less frequently found in dogroses rather originated from hybridizations within dogrose species involving unreduced egg cells (Herklotz and Ritz, 2017).

The maintenance of two 5S rDNA families in the Rosa genomes is consistent with the increased stability of 5S loci as compared to 35S loci in allopolyploid genomes documented in several allopolyploid systems (Pedrosa-Harand et al., 2006; Weiss-Schneeweiss et al., 2008; Garcia et al., 2016; Amosova et al., 2019). The reasons for relative stasis of 5S rDNA loci are not well understood, while their position on chromosomes (Garcia et al., 2016) and epigenetic modifications (5S rDNA loci carry mostly heterochromatic landmarks; Simon et al., 2018) have been discussed. One also has to consider the relative scarcity of meiosis driving genetic recombination (and homogenization) in these long-living perennial shrubs. For example, a R. canina individual known as “Rose of Hildesheim” (North Germany) is estimated to be more than 700 years old (Peters and Peters, 2013). Interestingly, Gossypium allopolyploids which represent also perennial shrubs tend to maintain 5S rDNA loci relatively intact over millions of years, while they homogenized their 35S rDNA loci (Cronn et al., 1996).