The Crocus panrepeatome reveals the links between whole-genome duplications, repeat bursts, and descending dysploidy

doi:10.21203/rs.3.rs-5147620/v1

Download PDF

Research Article

The Crocus panrepeatome reveals the links between whole-genome duplications, repeat bursts, and descending dysploidy

https://doi.org/10.21203/rs.3.rs-5147620/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Dysploidy is a crucial driver for species diversification by changing karyotypes through diploidization after a whole-genome duplication (WGD). A WGD could trigger repeat bursts but our knowledge about the evolutionary links between WGD, repeat bursts, and descending dysploidy is limited. Using Crocus as a model, we performed a panrepeatomic analysis to gain insights into the influence of WGD on repeat bursts and of repeat dynamics in descending dysploidy. We first established a phylogenetic backbone of the genus using chloroplast (cp) and 35S rDNA sequences to which we could anchor our panrepeatome data. We identified a WGD event before the initial divergence of Crocus, and nested WGD events before the divergence of some series coinciding with repeat bursts. We demonstrate repeat-linked chromosome fusions of C. longiflorus chromosomes in its dysploid relative C. vernus. This work demonstrates the links between WGD, repeat bursts, and descending dysploidy.

Evolutionary Genetics

Crocus

descending dysploidy

panrepeatome

repeat burst

repeatome

whole-genome duplication

Dysploidy, a change in base chromosome number resulting from chromosome fusion or fission, is an essential player in species diversification [1–4]. It is a crucial part of post-polyploid diploidization to counter the potential adverse effects of WGD [5–8]. Reproductive isolation resulting from WGD and dysploidy provides a platform for speciation and cladogenesis [2, 9]; hence, indicating a direct link between WGD, dysploidy, and species diversification [2, 10, 11].

WGD, with (allopolyploidy) or without hybridization (autopolyploidy) [12], causes genomic disturbance in the cell that requires genetic and epigenetic regulation to restore the balance of cellular function [5, 6]. This so-called ‘genomic shock’ is often accompanied by epigenetic activation of repeats, spurring bursts of repetitive elements, such as transposable elements (TEs) and tandem repeats (satellite DNA [satDNA]), ribosomal DNA [rDNA], and microsatellites) [6, 13, 14]. These bursts often amplify only a few repeat families among many in a “library” of low-copy repeats shared by related taxa, called the “library hypothesis” [15–18]. Taking an individual taxon’s total repeat content (repeatome) and comparing it with several other repeatomes (panrepeatome) in a phylogenetic context would reveal patterns of the library hypothesis and repeat bursts. Related taxa forming a clade share a pattern of abundance of specific repeat families compared to those from other clades, making shared repeat bursts a diagnostic feature of cladogenesis [2, 15].

It is yet poorly understood which repeat family amplifies in a burst event and whether this process is random or not [17]. However, repeats are frequent sites of chromosomal aberrations (inversions, deletions, duplication, translocations, and fusions) when (double-strand breaks) DSBs occurring within homologous repeats from non-allelic sites are ligated via homologous recombination [19–21]. Indeed, telomeric repeats, satDNA, and rDNA have been described at chromosome fusion points in some species, alleging their involvement in descending dysploidy [22–24].

Descending and ascending dysploidy have been extensively described in holocentric species [22, 25, 26]. Holocentric chromosomes are more tolerant to chromosome fusion or fissions than monocentric chromosomes because kinetochore assembly and mitotic spindle attachment occur at the entire chromosome length, thus ensuring the heritability of chromosome fragments [27]. However, the gain of functional telomeres is required in cases of chromosome fission. Monocentric chromosomes are less tolerant to chromosome fission because acentric chromosome fragments can lead to chromatin loss and gene dosage imbalance, which might be deleterious. Consequently, descending dysploidy is more prevalent than ascending dysploidy in monocentric species, not because fusion occurs more frequently than fission but because fusion allows for easier retention of genetic elements, ensuring survival [1, 2]. Indeed, chromosome numbers of flowering plants are tightly regulated to n = 7–12 by post-polyploid diploidization mechanisms despite several rounds of WGD, making descending dysploidy a vital process of re-diploidization [4, 28–32].

Chromosomal rearrangements such as fusions and inversions have been demonstrated to be involved in local adaptation and, eventually, speciation [33–38]. These events may capture locally adaptive genetic loci and reduce recombination between linked loci in a rearranged locus, thus favoring positive selection and fixation in a population, especially when heterozygotes for such rearrangements are underdominant compared to their homozygote counterparts. These rearrangements, when reinforced by natural selection, can ultimately lead to speciation [37, 39].

Hence, several reports have presented the theoretical association and mechanisms linking WGD, repeat bursts, dysploidy, reproductive isolation, and speciation (Supplementary Fig. 1A). Nevertheless, direct empirical evidence on these mechanisms is still very limited and focused either on a few related taxa, on one or a few aspect(s) of this chain of events, or on holocentric chromosomes [4, 8, 34, 37, 40]. To better understand the connection between these components in monocentric species, it is necessary to comprehensively analyze a phylogenetic group that comprises a broad number of dysploid taxa.

The monocot genus Crocus L. (Iridaceae) comprises 257 taxa (World Plants: https://www.worldplants.de, Accessed: 17.09.2024) exhibiting high chromosome number diversity (2n = 6, 8, 10, 12, 14, 16, 18, 20, 22, 23, 24, 26, 28, 30, 32, 34, 36, 44, 48, 64; Supplementary Fig. 1B, C) [41–45]. In some cases, closely related species with double the chromosome number have similar genome sizes, complicating the correlation between chromosome numbers and ploidy levels [45]. Moreover, closely related sympatric and parapatric species have been observed and are often characterized by chromosome count differences, suggesting that dysploidy contributes to reproductive isolation and might drive Crocus diversification.

Here, we explored the Crocus panrepeatome in a phylogenomic context to find patterns that connect WGD with descending dysploidy. First, we established the phylogenetic relationships within Crocus, inferred the occurrence of WGD events, and anchored them in the phylogenetic tree. Next, we analyzed the correlation between inferred WGD, repeat bursts, and cladogenesis. Then, we investigated the potential role of satDNAs in descending dysploidy using panrepeatome and cytogenetic data. Finally, we provided cytogenetic evidence of chromosome fusion leading to descending dysploidy in ser. Verni using chromosome painting.

Phylogenetic relationships in Crocus

Until now, the phylogenetic backbone of Crocus was not yet resolved. Therefore, we first established a molecular classification of Crocus to have a robust foundation for interpreting the relationships between WGD, repeat dynamics, and dysploidy in a phylogenomic context. We generated phylogenetic trees based on chloroplast genomes and 35S rDNA sequences of 175–178 taxa (Supplementary Table 1). Four Romulea species, as well as Syringodea bifurcata and Afrocrocus unifolius, served as outgroups. We included only those with complete 35S rDNA and plastid genome assemblies for Crocus samples.

In the 35S rDNA phylogeny (Supplementary Fig. S2), S. bifurcata and A. unifolius group, with 93% bootstrap support, as sister to the maximally supported genus Crocus. Both species were lacking in the plastid dataset. Therefore, we refer for divergence time estimations to the 35S rDNA dataset.

The genus Crocus originated about 17.5 Mya ago (16.59–19.74 Mya) and is divided into four subgenera (Orientalis, Carpetanus, Crocus, and Nudiscapus) in the chloroplast as well as in the 35S rDNA dataset (Supplementary Fig S2). In the 35S rDNA phylogeny, subg. Orientalis split first (15.56 Mya [14.9–16.7 Mya]), followed by Carpetanus (15.2 Mya [14.2–16.7 Mya]) and Crocus (13.8 Mya [13.2–15.2 Mya]), with Nudiscapus as the sister clade (bootstrap supports of 97–100%) (Supplementary Fig. S2). The cp phylogeny differs with Carpetanus splitting first, followed by Orientalis (Supplementary Figs. S3, S4). However, while the support of the backbone of the plastid phylogeny is very high (Supplementary Fig. S3), this split of Carpetanus and Orientalis is only weakly supported (51%) in the 35S rDNA phylogeny. Still, support for the four main clades and most of the subclades is strong in both datasets. The only exception with more incongruences and lower support values concerns the subclades of Nudiscapus (Supplementary Fig. S4). In this subgenus, we could detect hybridization between different series (C. sieheanus 2n = 16, series Abantensis or ser. Reticulati x ser. Punctati, C. mazziaricus 2n = 16, 2C = 11.77 pg originated by a cross of species from ser. Punctati species and ser. Cancellati, Supplementary Fig. S4).

In summary, although we could resolve the phylogeny, we found several incongruences that are, at least in some cases, likely caused by hybridization. We used the 35S rDNA tree to anchor our WGD, repeat dynamics, and dysploidy data, considering the biparental evolution of the nuclear genome-derived 35S rDNA, which is more appropriate for analyzing these nuclear genomic events.

Two periods of WGD events in Crocus

We identified WGD events using paralogous genes from genomic and transcriptomic data of five species (subgenus Crocus: C. hadriaticus and C. sativus from ser. Crocus, C. longiflorus and C. vernus from ser. Verni, and subgenus Nudiscapus: C. weldenii from ser. Punctati) to understand whether WGD events occurred before the diversification of Crocus and of its intrageneric lineages.

With the duplicated BUSCO gene sets from the five species, we could detect a WGD event between 14.89–19.74 Mya (Ks = 0.118–0.168) before Crocus diversification, henceforth referred to as Cr-β WGD event (Fig. 1A, B). Part of the vacuolar protein sorting 11 gene (vps11, Ks = 0.122 in C. longiflorus (evolutionary rate [ER; [46]]: 0.84, Ks normalized by the ER: 0.145), used as a phylogenetic marker for samples without omics dataset, also showed duplication in species representing all four Crocus subgenera, supporting the shared Cr-β WGD event by all crocuses (Fig. 1C). Similar results were obtained using other Angiosperm 353 marker genes (Ks ≈ 0.145), corroborating the shared Cr-β WGD event among crocuses (Supplementary Fig. S5). The gene trees using these marker genes resulted in grouping one Crocus paralog with one or both outgroups S. bifucata and A. unifolius and the other Crocus paralog as a sister (Supplementary Fig. S5). This observation indicates that the divergence of Crocus paralogs is older than the actual Cr-β WGD and that the Cr-β WGD event was allopolyploidization involving an ancestor of S. bifucata and A. unifolius. Because allopolyploidization suggests that the paralogous genes have already accumulated some divergence before the merger of the two parental genomes, the actual Cr-β WGD event should be younger than the dated stem age of 16.59–19.74 Mya but before crocuses started to diversify around 14.89–16.74 Mya.

In addition to the Crocus-specific Cr-β WGD, we also detected more recent independent and intrageneric WGD duplication events. We collectively refer to this period of recent intra-Crocus WGD events as the Cr-α WGDs (Fig. 1A, B). Using the Ks values of duplicated BUSCO genes from C. sativus (density peak at Ks = 0.03) and C. hadriaticus (density peak at Ks = 0.02), we inferred the ser. Crocus Cr-α WGD event at 2.0–3.9 Mya at the root of ser. Crocus. Similarly, for ser. Verni, we estimated the Ks values from C. vernus (density peak Ks = 0.02) and inferred the ser. Verni Cr-α WGD event at 2.1–2.7 Mya after divergence from C. longiflorus (Fig. 1A, B).

With the C. weldenii dataset, we could clearly detect the Cr-β WGD, but the Cr- α WGD-related Ks peak is not distinctly visible, just showing a “bump.” The number of duplicated BUSCO genes with Ks values of 0.015–0.035 around this “bump” is much lesser in C. weldenii (110 out of 3,055 duplicated genes) than in other species that experienced the Cr-α WGD (981 out of 3020 in C. vernus, 772 out of 2,995 in C. sativus, and 163 out of 1,314 in C. hadriaticus transcriptome, Supplementary Table 2). These 110 duplicated genes in C. weldenii are more than those in C. longiflorus (16 out of 2,969), which did not experience the Cr-α WGD. This observation suggests that C. weldenii may have experienced introgressed genomic segments that experienced the Cr-α WGD.

In summary, we have shown that at least two periods of WGD events have predated the initial (Cr-β) and several recent intrageneric (Cr-α) diversifications in Crocus.

Repeat bursts after WGD events

To understand whether WGDs spurred repeat bursts in Crocus, we analyzed the Crocus panrepeatome, which includes TEs, satDNAs, rDNA, and microsatellites. First, we quantified the different repeat classes by subgenus, sub-lineages, and samples (Fig. 2A–C). Then, we calculated statistical support for repeat bursts within each sub-lineage (Fig. 2D–F) and anchored the evolutionary interpretation on the 35S rDNA tree.

We analyzed the repeatomes of 219 samples, including 215 Crocus and four Romulea outgroups (Fig. 2A, B; Supplementary Table 3). RepeatExplorer2 and TAREAN generated 142,149 repeat clusters represented by contigs with the highest genome representations (Supplementary Table 4), which grouped into 47 repeat categories based on RepeatExplorer2 annotation and into six major panrepeatome categories: Class I TE, Class II TE, satDNA, rDNA, Unknown, and Organelle (Supplementary Table 5). The sequences clustered into 40,599 uclusters, excluding organellar repeats (see Methods for the difference between “cluster” and “ucluster”). Of these uclusters, 37,878 are specific to Crocus, 2,552 to Romulea, and 169 are shared (Supplementary Fig. S6 and Supplementary Table 6). For comparative purposes, we included Romulea in the panrepeatome dataset, which comprised 7,354 and 33,245 core and accessory uclusters, respectively (Supplementary Table 6).

Sequences within all uclusters ranged from 13 to 20,697 bp (Fig. 2C, Supplementary Fig. S7, Supplementary Table 4). Transposable elements, rDNA, and Unknown repeats are in the longer range, and satDNAs are in the shorter range. Most Unknown repeats fall into the 100–600 bp range, with 126-bp the predominant sequence representation, and satDNAs were dominated by 178-bp repeats.

The genome proportions (GP) of total repeats are significantly higher in Crocus (34.9–73.2%) than in Romulea (18.5–34.1%) (Fig. 2A), and the repeat content significantly contributed to the genus-level genome size (GS) increase (p < 0.001) (Supplementary Fig. S8). However, GS positively correlated with the repeat content only in subg. Crocus and not in other subgenera (Supplementary Fig. S9). In the cytogenetic group, ser. Verni, total repeat content did not affect GS when samples are strictly limited to ser. Verni. However, this relationship is significantly positively correlated when the immediate sister group, ser. Malyi, is included (Supplementary Fig. S10). However, this relationship between genome size and total repeats should be interpreted cautiously, as many samples do not have genome size data (Supplementary Table 1). Chromosome numbers generally did not correlate with genome size at the genus level but showed a strong negative correlation with GS in combined ser. Verni and ser. Malyi (Supplementary Fig. S10).

Regarding overall repeat content, the subgenera Crocus and Nudiscapus have significantly more repeats than Carpetanus and Orientalis combined, and Nudiscapus has significantly more repeats than Crocus. Moreover, repeat content varied between different sub-lineages (Supplementary Fig. S11). Class I TEs (9.6–57.7% GP) significantly contributed to genome size more than Class II TEs (0.3– 3.4% GP) (Fig. 2A, B, Supplementary Figs. S12, 13). Among Class I TE, LTR retrotransposons are the predominant repeats. Ty1/Copia is significantly higher than the Ty3/Gypsy superfamily (Fig. 2B). Angela (Ty1/Copia), Athila, and Retand (Ty3/Gypsy) LTR families have significant contribution to genome size (Supplementary Figs S12). In contrast, Class II TEs, with EnSpm/CACTA having the highest GP (Fig. 2B), have negative correlations with genome size, although not as substantial as those of Class I TEs (Supplementary Fig. S13). In ser. Verni, Angela, Bianca, Retand, SIRE, Tekay, EnSpm/CACTA, and uncharacterized LTR elements significantly increased relative to genome size (Supplementary Fig. S14).

Both satDNAs (0.2–9.7% GP) and rDNAs (0.4–11.9% GP) also contributed considerably to the genome sizes of some samples. Series Adanensis, ser. Crocus, and ser. Scardici have the most abundant satDNA, 35S rDNA, and 5S rDNA content, respectively (Supplementary Figs. S11 and Fig. S15), although C. wattiorum (currently unplaced in sect. Heterogenei, provisional named ser. Wattiorum) has the highest 5S rDNA content (Fig. 2A). Nevertheless, a considerable proportion of repeats (3.6–46.8 GP) remain uncharacterized (Fig. 2A, C). These unknown repeats may include short, non-autonomous mobile elements like miniature inverted-repeat transposable elements (MITE) [47], terminal-repeat retrotransposons in miniature (TRIM) [48, 49], or short interspersed nuclear elements (SINES) [50].

Of the 40,599 uclusters, 8,918 were significantly associated with each group (p < 0.05), and 2,308 were strictly specific (p ≥ 0.05 in other groups) to each group (Fig. 2D and Supplementary Table 6, 7). These group-specific repeats could resolve series within the subg. Crocus, but not within subg. Nudiscapus, except for section Adamii (Fig. 2E).

The 2,308 group-specific repeats represent 2–70% of the total repeats in each sample and a mean of 5–60% per series (Supplementary Table 8), indicating that the group-specific repeats varied in the magnitude of amplification in each sample and sub-lineage (Fig. 2F). These repeats represent 18–70% (mean 60%) of the total repeats in ser. Crocus, 38–70% (mean 60%) in sect. Adamii, and 16–38% (mean 31%) in ser. Verni, indicating that bursts of specific repeats contribute differently to the overall repeat content in different clades. Moreover, both ser. Crocus and ser. Verni, where WGD events were involved in cladogenesis, showed bursts of these repeats, indicating a link between WGD and repeat bursts.

In summary, we have shown that different types of repeats contributed differently to each sample, series or section, and subgenus. We have also demonstrated a direct correlation between WGD and repeat bursts in ser. Crocus and ser. Verni where we have observed the Cr-α WGD. The panrepeatome profile (heatmap) can also enable one to extrapolate potential WGD events in clades where no genomic data is available but shows a similar profile and amount of repeat burst with those in the two series (e.g., ser. Sieberi, sect. Adamii, and ser. Cancellati). Next, we wondered if TEs co-evolved with some tandem repeats like microsatellites.

Microsatellite repeats co-amplified with transposable elements

To comprehensively survey the panrepeatome, we included microsatellites, but we extended the length range to cover 2–30 bp repeats because genomic repeats are not limited to only the longer satDNAs and TEs. We extracted 357,682 microsatellite sequences Supplementary Table 9) and saw a higher representation of hexanucleotides (105 Mb cumulative length), which is second only to dinucleotide microsatellites (189 Mb) (Supplementary Fig. S16). Moreover, microsatellites that are 6-bp, 12-bp, and 18-bp long, all multiples of six, are the most represented sequences in the dataset, indicating a biological function of hexanucleotide microsatellites in Crocus.

All microsatellites clustered into 35,985 uclusters with the 6-bp human-type telomeric repeat (TTAGGG, micro10581, htel) representing the fourth longest cumulative length, only next to AG, AT, and AAC microsatellites as the three highest, respectively (Fig. 3A and Supplementary Tables 9 and 10). Whereas the Arabidopsis-type telomere (TTTAGGG, micro03525, aratel) ranked only 44th. Intriguingly, both telomeric repeat variants only had 53-bp and 50-bp mean array sizes, respectively, which did not cover the input read length (103 bp) (Fig. 3A and Supplementary Table 10), indicating an abundance of short-arrayed loci over the long arrays found at chromosome termini.

FISH analysis of the htel and aratel microsatellites confirmed the distribution of htel at chromosome termini but not of aratel (Fig. 3B). In addition to the chromosome termini, a chromosome painting-like pattern was also observed for htel in all species we analyzed (Fig. 3B and Supplementary Figs. S17-S37). We suspected this chromosome painting-like pattern of htel is caused by the insertion of short htel arrays in dispersed repeats like retrotransposons. So, we analyzed the co-amplification of microsatellites and TEs between Romulea and the Crocus subgenera and found that several microsatellites were significantly co-amplified with TEs (p < 0.05, Fig. 3C and Supplementary Table 11). For example, micro03527 and micro18034 co-amplified with Ucluster00803, which is a TAR element (Ty1/Copia), and micro02204 co-amplified with Ucluster00104 (SIRE) and Ucluster00323 (Bianca) (Fig. 3D). Moreover, micro10581 (htel) and micro03525 (aratel) co-amplified with Athila (Ucluster00339), Bianca (Ucluster00323, Ucluster00395, and Ucluste00398), and Angela elements (Ucluster00303).

This co-amplification is explained by the insertion of short co-evolved microsatellite arrays into their corresponding co-amplified TEs (Fig. 4). For example, micro10581 forms short arrays between the polyprotein domain and the 3’ LTR of the Athila Ucluster00339. Conversely, the insertion of micro10581 in Bianca and Angela elements is between the 5’ LTR and the gag gene. These short arrays within these TEs explain the overall short mean array size of htel microsatellites (Fig. 3A), contrary to what we expected if htel was exclusively found at chromosome ends. In such a case, we expected htel to cover all 103-bp reads, not just the mean 53 bp.

Micro18034 and micro03527 are mainly localized in short arrays at the 3’ LTR of TAR elements (Ucluster00803). Microsat02204 is localized as part of a longer repeat monomer in short arrays towards the 3’ LTR in SIRE elements (Ucluster00104) and after the 5’ LTR in Bianca elements (Ucluster00323). We did not observe a high amplification of micro02204 in our C. vernus LTR dataset because this microsat is amplified more in subg. Nudiscapus than in subg. Crocus (Fig. 3D).

In summary, we have shown that human-type telomeric repeat makes up the chromosome ends in Crocus, and that this repeat, along with the much fewer aratel, are inserted in TEs. Some other microsatellites, like micro02204, also show co-evolution with TEs and lineage-driven amplification. Next, we wondered if the pattern of repeat amplification and distribution in the genome could tell us about their potential involvement in descending dysploidy.

Amplification of CroSat042 suggests a role in descending dysploidy

To explore whether the patterns of satDNA abundance and chromosomal distribution among related taxa in a clade can tell us about the potential role of repeats in descending dysploidy, we analyzed the karyotypes of species in ser. Verni. We focused our analysis on ser. Verni because of this group's high rate of dysploidy (Supplementary Fig. S10) [45] and the readily available materials for chromosome analysis. To do this, we focused on the satDNA fraction of the panrepeatome (pansatellitome) (Fig. 5A, and Supplementary Table 6), as satDNAs are the repeat types found mainly at chromosome fusion points [22].

The 2,146 satDNA sequences ranging from 13–7948 bp were clustered into 564 CroSats with 178 bp as the most abundant monomer length in terms of sequence count and cumulative length (Fig. 5B). Of the 564 CroSats, some show sub-lineage enrichment, while others are shared by several groups (Fig. 5A). Three CroSats (CroSat042, CroSat144, and CroSat137) are relatively abundant in ser. Verni but are also shared by other groups. Most groups in subg. Crocus share the 178-bp CroSat042, but it can also be detected in very high copies in ser. Adanensis within subg. Nudiscapus. The 186-bp CroSat144 is more abundant in ser. Sieberi than in ser. Verni. Finally, the 159-bp CroSat137 is shared in ser. Verni and sister outgroup (ser. Malyi) and in ser. Versicolores.

These three satDNA markers and the 5S and 35S rDNA markers enabled karyotyping of eight C. vernus populations; three C. heuffelianus populations; one population each for C. bertiscensis, C. kosaninii, C. tommasinianus, and C. longiflorus; two populations of a tetraploid C. cf. heuffelianus; and the cultivar C. etruscus ‘Zwanenburg’ (which is a hybrid of C. tommasinianus and C. etruscus) (Supplementary Fig. S17). However, the homologous pairing was not always straightforward for some species as these cytogenetic markers also revealed the highly heteromorphic karyotypes of many species from different populations, like those of C. heuffelianus and C. vernus. No two populations showed the same karyotype (Supplementary Figs S17-S37).

CroSat042, CroSat144, and CroSat137 varied in the number of chromosome loci, distribution, and cumulative lengths within ser. Verni. Species with more chromosomes have fewer CroSat042 loci and, thus, shorter cumulative lengths, such that the number of loci and the cumulative lengths are strongly negatively correlated with chromosome number (Fig. 5C). Whereas CroSat144 has no significant correlation, and CroSat137 is strongly positively correlated.

Regarding chromosomal distribution, CroSat144 was detected on only one chromosome in C. longiflorus, at the interstitial region of the long arm of chromosome 8 but varied in presence or absence and distribution in its dysploid relatives (Supplementary Figs. S17-S37, Supplementary Table 12). CroSat137 was mainly distributed at the primary constriction, indicating centromeric distribution, but also at intercalary loci in some dysploid species. CroSat042 was localized on only one chromosome in C. longiflorus, a small locus in the subtelomeric regions of both chromosome arms (Fig. 5D, Supplementary Fig. S19, Supplementary Table 12). However, dysploids, particularly C. vernus, which has the lowest chromosome count, have more loci, and FISH signals are more intense at subtelomeric, pericentromeric, and intercalary regions (Supplementary Figs. S20-S37).

The increase in loci number and cumulative length of CroSat137 is reasonable because this satDNA is localized in the centromere; hence, the more chromosomes, the more loci and the longer the cumulative length. On the other hand, the increase in loci number and cumulative lengths of CroSat042 as chromosome numbers decrease suggests a role in chromosome fusion, especially since many loci are present at intercalary regions, potentially mediating chromosomal rearrangements at these sites.

In summary, we have shown that comparative pansatellitome and cytogenetics analyses revealed a putative involvement of CroSat042 in satDNA-linked chromosome fusion and descending dysploidy. Nevertheless, relying solely on this correlation does not unambiguously demonstrate the fusion of chromosomes in dysploids.

Crocus longiflorus chromosomes formed fused chromosome blocks in C. vernus

To demonstrate the fusion of C. longiflorus chromosomes in a dysploid relative, we performed chromosome painting using the genome assembly of C. longiflorus to design the paint probes. Then, we compared the FISH painting signals between C. longiflorus and C. vernus. We developed probes for six chromosome arms: CL_2S for the short arm of chromosome 2, CL_2L for the long arm of chromosome 2, CL_3L for the long arm of chromosome 3, CL_11S for the short arm of chromosome 11, CL_13L for the long arm of chromosome 13, and CL_14L for the long arm of chromosome 14 (Fig. 6A).

FISH signals of all six C. longiflorus chromosome arms showed patterns precisely as predicted (Fig. 6B). Each pool showed two signals, indicating one homologous pair and confirming the diploid genome of C. longiflorus. FISH analysis of the same six targets confirmed the fusion of C. longiflorus chromosome arms in C. vernus. Moreover, the six targets showed a more complex distribution in C. vernus (Fig. 6C, D). Four targets (CL_2L, CL_11S, CL_13L, and CL_14L) showed two signals indicating conserved synteny with that in C. longiflorus; one target (CL_3L) showed four signals but on the same chromosome arm, which could indicate chromosome inversion; and finally, one target (CL_2S) showed six signals, indicating either fragmentation of the short arm of C. longiflorus chromosome 2 before fusion in C. vernus or duplicated segments after WGD that rearranged to different chromosomal regions.

Because five of our probes targeted only one of the two chromosome arms for each chromosome, we could not tell whether the entire chromosomes were fused as a unit or whether chromosome arms fragmented before rearrangements. Nevertheless, we could tell from the C. longiflorus chromosome 2 (CL_2S and CL_2L) that one chromosome was fused as a unit as indicated by adjacent CL_2S and CL_2L in C. vernus chromosome 3a. Whereas in C. vernus chromosome 3b, CL_2S and CL_2L are separated into the different chromosome arms.

This comparative chromosome painting approach confirmed three C. vernus karyotype features: 1) The fusions of C. longiflorus chromosomes forming chromosome blocks in C. vernus. 2) These chromosome blocks are often bordered by CroSat042, reiterating the putative role of CroSat042 in a satDNA-linked chromosome fusion. 3) The absence of exact homologous chromosome pair but rather complex heterologous chromosomes, reiterating the highly heteromorphic C. vernus karyotypes. While some chromosome arms show similar painting patterns indicating pairing in meiosis, most do not (Fig. 6C, D).

In summary, we have shown through chromosome painting that the dysploid C. vernus genome is composed of C. longiflorus chromosomes that “fused” to form syntenic blocks, and that cytogenetically visible loci of CroSat042 are in between these blocks.

Here, we presented the panrepeatomics of the monocentric genus Crocus to demonstrate the links between WGD, repeat bursts, and descending dysploidy within the framework of the karyotype evolution cycle (Supplementary Fig. S1). We emphasize the importance of using multiple datasets such as phylogenetics, genomics, repeatomics, and cytogenetics to perform a comprehensive comparative analysis to gain a clear insight into the complex relationships between these evolutionary events that, otherwise, would not be detectable when using only a few of these data sources.

Establishing a robust phylogenetic backbone was crucial to interpreting the panrepeatome data in an evolutionary context. Our more extensive cp and 35S rDNA datasets allowed us to update the phylogenetic relationships in the genus, which have not been achieved in past works [51, 52]. The division of the genus into four subgenera is supported in both datasets. The two first-splitting subgenera, Orientalis (15.6 Mya) and Carpetanus (15.2 Mya), are species-poor. The common ancestor of the two slightly younger subgenera, subg. Nudiscapus and subg. Crocus, split ca. 13.8 Mya and diversified into the species-rich crown group. Both the cp and 35S rDNA tree topologies also mirror the typical Whole Genome Duplication-Radiation Lag Time (WGD-RLT) Model tree topology [53], where a species-rich crown group radiates from a species-poor sister group, indicating a preceding WGD before genus diversification, the Cr-β WGD event in Crocus. This pattern has been observed in many higher-level groups in Angiosperms and in Brassicaceae [53–55].

Most diversification in the genus took place within a relatively short timeframe (Fig. Supplementary Fig. S4), which increases the likelihood for incomplete lineage sorting (ILS). ILS can cause incongruences between cp and rDNA-inferred trees [56] as observed particularly within the subgenus Nudiscapus. Whether these incongruences are a result of ILS or hybridization is difficult to distinguish particularly when both processes occur around the same time in the evolutionary history [57]. These incongruences may not be easily resolved using conventional molecular markers, but previous works have shown the potential of repeat abundance in adding a layer of information in resolving phylogenetic relationships [58–60].

Because repeat bursts could be induced by WGDs [6, 7, 13, 14], recently derived taxa that share a WGD event share a similar repeatome profile [59, 61, 62], which can be best visualized in a heatmap as a panrepeatome profile. With this profile, we could correlate the repeat bursts in ser. Crocus and ser. Verni with the Cr-α WGD events. Both clades also show considerable but varied magnitude of repeat bursts (mean: 60% in ser. Crocus and 31% in ser. Verni). Moreover, using the panrepeatome profile, we could also extrapolate the occurrence of independent Cr-α WGD events in clades where we do not have duplicated gene sets to infer WGD events, such as in ser. Sieberi and ser. Scardici in subg. Crocus, which show well-defined panrepeatome profiles with 45% and 28% repeat burst proportions, respectively (Fig. 2E, F).

Conversely, the panrepeatome profile of many sub-lineages within subg. Nudiscapus, such as most series of sect. Heterogenei, suggest they did not experience the Cr-α WGD. The inability of the group-specific repeats to resolve these groups (Fig. 2E), as well as their low repeat burst proportions (Fig. 2F), favor this hypothesis. This observation suggests that these repeats have remained shared since the Cr-β WGD event. Cladogenesis in these sub-lineages may have been caused by other mechanisms besides WGD, such as homoploid hybridization or introgression [63, 64], the latter being supported by high incongruences between the cp and 35S rDNA trees (Supplementary Fig S4). Moreover, allopolyploids between different series in subg. Nudiscapus exist (e.g., C. sieheanus 2n = 16 and C. mazziaricus, Supplementary Fig S4) and might facilitate introgression through backcrosses to their diploid parental species.

The Ks peaks of C. weldenii support the absence of the Cr-α WGD in ser. Punctati. Potential Cr-α gene introgression donors might come from ser. Cancellati, considering the panrepeatome profile pointing to Cr-α WGD occurrence in this series. Series Cancellati is also characterized by incongruent gene trees (Supplementary Fig. S4), which could indicate introgression and or even an allopolyploid origin for this group. Section Adamii might also have experienced a Cr-α WGD, considering its well-defined panrepeatome profile with the highest repeat burst proportion (mean 60%, Fig. 2F). While sect. Adamii is always monophyletic, its position varies between different gene trees (Supplementary Fig. S4). One possible explanation could be that the Adamii ancestor was affected by hybridization or allopolyploidization.

Some lineages do not show well-defined panrepeatome profiles but instead have similar patterns with a sister clade. For example, ser. Serotini and ser. Versicolores show a high degree of shared repeats. This similar pattern is also seen between ser. Malyi and ser. Verni in subg. Crocus. Moreover, we have shown that C. longiflorus, which did not share the Cr-α WGD, shows a similar repeatome profile to those in ser. Malyi. This suggests that the ser. Malyi and C. longiflorus have diverged from the initial Cr-β WGD, and that the dysploids in ser. Verni diverged after a nested Cr-α WGD following split from C. longiflorus. Evidently, dysploidy has become an important pathway for species diversification in this group after the Cr-α WGD.

Indeed, the dynamics of CroSat042 in ser. Verni have indicated the active involvement of satDNA in descending dysploidy in this group. Although only one chromosome (chr. 7) of C. longiflorus have long CroSat042 arrays at the subtelomeres of both arms to be detectable by FISH, fewer copies are likely present in other chromosomes. Occurrence of DSB at non-allelic sites may have caused amplification of CroSat042 via the break-replication/fusion mechanism [65] before fusion, causing expansion of CroSat042 arrays in fusion points observed in dysploid species in ser. Verni (Supplementary Fig. 17–37). The fact that CroSat042 is also present in other Crocus sub-lineages indicates its role in dysploidy in a wider range of species.

An abundance of microsatellites that are inserted in TEs may also increase the chances of microhomology-mediated break-induced replication repair of DSB [66]. For example, the insertion of htel (which has also been found in other species in Asparagales [67] where Crocus belong) into a diverse family of repeats may allow a higher chance of misrepair of DSBs during replication stress such as during an ecological upheaval [68]. Whether these repeats play a direct role in the high karyotype heteromorphy in ser. Verni, especially in C. vernus, remains to be analyzed. Nevertheless, the complex interplay of sequence-based propensity to replication-stress-induced DSB [69] and the increasing chance of DSB mis-repair with increasing repeat copy number may have contributed to the high dysploidy rate in Crocus.

Here, we have shown the usefulness of panrepeatomics in revealing the links between WGD, repeat bursts, and descending dysploidy in the monocentric genus Crocus by integrating phylogenetics, repeatomics, cytogenetics, and genomics datasets. An older Cr-β WGD event promoted the initial divergence of the genus, and several younger and independent Cr-α WGD events are nested within the initial diversification, promoting cladogenesis. Nevertheless, clades did not seem to always follow a recent WGD but may have formed through other mechanisms such as introgression. This relationship between WGD and cladogenesis was well supported by the panrepeatome profile, corroborating the repeat amplifying effect of WGD and the shared pattern of repeatome profiles among species in a clade. We also showed that the amplification of repeats, such as the satDNA CroSat042, could, during the process, induce chromosomal rearrangements, resulting in dysploid lineages.

Whether CroSat042 is indeed frequently found in chromosome fusion points require a high-quality genome assembly, or whether other repeats have the propensity to be involved in dysploidy needs further investigation. To understand the evolutionary effects of dysploidy in Crocus diversification, it is necessary to test the viability and fertility of F1 plants from crosses of different species with different chromosome numbers. The high karyotype heteromorphy in ser. Verni, especially in C. vernus, may indicate homoploidy. However, the fertility gradient observed in C. vernus individuals may suggest an intrapopulation “homoploidy” that results in different heterozygote combinations of chromosome blocks. Thus, producing fertile or sterile individuals depending on the degree of chromosome block arrangement that favors proper meiosis. Understanding this would require analysis of meiotic chromosome pairing. Whether the observed dysploidy pathway in ser. Verni is universal across Crocus, replicating the cytogenetics work in other series and completing the genome size and chromosome count of unrepresented samples would be necessary. All of this information will allow us to fully understand the karyotype evolution cycle in Crocus.

Plant materials

Plants from previous collections [45] were grown in the experimental garden at IPK, Gatersleben, Germany, and used for either genome size measurements, DNA extraction, or cytogenetic analyses. A summary of

DNA extraction and sequencing

For short read sequencing, genomic DNA was extracted using the DNeasy Plant DNA Extraction Kit (Qiagen) from about 10 mg of silica-dried leaf material according to the manufacturer's protocol. We generated short Illumina reads to assemble the chloroplast genomes, 35S rDNA, and panrepeatome. DNA concentration and quality were afterward checked on 0.8% agarose gels. About 1–2 µg DNA was used to create the Illumina library according to the Illumina TruSeq DNA preparation protocol following the manufacturer’s recommendations or as described by Meyer and Kircher [70]. DNA was fragmented to 400–500 bp using Covaris S220 Focused-ultrasonicator (Covaris, MA, USA), followed by adaptor and barcode ligation. The library was size selected with a BluePippin with 0.75% agarose cassettes (Sage Science; Beverly, MA, USA). Fragment size distribution and DNA concentration were evaluated on an Agilent BioAnalyzer High Sensitivity DNA Chip (Agilent Technologies Inc, CA, USA) using the Qubit DNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, Thermo Fisher Scientific, MA, USA). Finally, the DNA concentration of the library was checked by a quantitative PCR run. Cluster generation was performed on Illumina cBot and paired-end sequencing (2 x 250bp) on either the Illumina HiSeq 2000/2500 or NovaSeq6000 platforms, following Illumina’s recommendation, which included 1% Illumina PhiX library as internal control. The targeted coverage was 0.5×–5× for each sample. Barcoded reads were de-multiplexed using the CASAVA pipeline 1.8 (Illumina, Inc.).

Fourteen nuclear single-copy markers were amplified [45] using the Phusion High Fidelity DNA polymerase (ThermoFisher Scientific, MA, USA). PCR products of each amplicon of one sample were pooled together regardless of their concentration, purified using a NucleoFast plate (Marcherey-Nagel, Düren, Germany), and diluted in 34 µL triple-distilled water. Sixteen microliters of the purified amplicon pool were digested using the NEBNext™ dsDNA Fragmentase kit (New England Biolabs, MA, USA) for an incubation time of 1.5 min, and another 16 µL of the amplicon pool was digested for 4.5 min following the manufacturer’s instructions. Both were pooled, 100 µL triple distilled water was added, and a NucleoFast plate (Marcherey-Nagel) was used to remove too-small fragments and contaminants.

Fifty microliters were subjected to size-selection targeting fragment sizes of 400–600 bp using BluePippin (Sage Science), blunt-end repaired, and used for sequencing library preparation according to the Illumina TruSeq DNA library protocol. For 169 samples, 2% of the amplicon library was added to short-read WGS libraries.

For high-quality genome assemblies of C. longiflorus, C. vernus and C. weldenii, we generated PacBio HiFi and Hi-C reads. High-molecular weight genomic DNA was obtained from 3–8 g of leaves, ground with liquid nitrogen to a fine powder. Nuclei were isolated, digested with proteinase K and lysed with SDS. For purification of HMW DNA phenol-chloroform extraction and precipitation with ethanol [71] was used. The HMW DNA was then dissolved in 50 ml TE (pH 8,0) and precipitated by the addition of 5 ml 3 M sodium acetate (pH 5,2) and 100 ml ice-cold ethanol. The suspension was mixed by slow circular movements resulting in the formation of a white precipitate (HMW DNA), which was collected using a wide-bore 5 ml pipette tip and transferred for 30 sec into a tube containing 5 ml 75% ethanol. The washing was repeated twice. The HMW DNA was transferred into a 2 ml tube using a wide-bore tip, collected with a polystyrene spatula, air-dried in a fresh 2 ml tube and dissolved in 500 µl 10 mM Tris-Cl (pH 8.0). For quantification the Qubit dsDNA High Sensitivity assay kit (Thermo Fisher Scientific, MA, USA) was used. The DNA size-profile was recorded using the Femto Pulse system and the Genomic DNA 165 kb kit (Agilent Technologies Inc, CA, USA).

As a final step, Hi-C libraries were prepared following the standard protocol of the Dovetail Omni-C kit for plants. Briefly, nuclei were isolated from fresh tissue and crosslinked with formaldehyde to preserve chromatin interactions. The crosslinked chromatin was then digested using a restriction enzyme (DpnII) to generate sticky ends. The overhanging ends were filled in with biotinylated nucleotides, followed by ligation under dilute conditions to promote the formation of ligation products between spatially proximal DNA fragments. Following ligation, crosslinks were reversed, and the DNA was purified. The biotinylated DNA fragments, representing ligation junctions, were captured using streptavidin beads and then subjected to further processing, including shearing to a size suitable for Illumina sequencing (ca. 300–500 bp). The fragments were then end-repaired, A-tailed, and ligated to Illumina sequencing adapters. The final Hi-C libraries were amplified by PCR and quantified using qPCR or a bioanalyzer before sequencing. Sequencing was performed on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA), generating paired-end reads of 2×150 bp.

Genome assembly

Genomes (C. longiflorus, C. vernus, and C. weldenii) were assembled using the generated PacBio HiFi and Hi-C data. First, a phased primary assembly was obtained by running Hifiasm [72, 73] using 50 Gb of PacBio HiFi reads as inputs in combination with Dovetail Omni-C reads with the following command: hifiasm -o out.phased.asm.hic --h1 hic.R1.fastq.gz --h2 hic.R2.fastq.gz hifi.reads.fastq.gz. The phased assemblies of each individual haplotype were further scaffolded to chromosome-scale using YaHS [74].

Chloroplast genome and 35S rDNA assembly

The obtained raw sequencing reads were quality-checked and over-represented reads were detected with FASTQC [75]. Adapter trimming of sequence reads was performed with CUTADAPT [76], and reads shorter than 120 bp after adapter removal were discarded. Error correction and removal of clonal reads were carried out using BBtools [77] (Bushnell B. – BBMap – sourceforge.net/projects/bbmap/).

We used GetOrganelle v1.7.7.0 [78] to de-novo assemble a first draft of the plastid genomes and the 35S rDNA sequences using the 2×250bp Illumina short read data as well as taking 2×150bp reads of Afrocrocus unifolius (ERR7617986), Romulea monadelpha (ERR7617928) and Syringodea bifucata (ERR7617929) from Kew Tree of Life Explorer [79]. This toolkit implements Bowtie [80] to initially find reads mapping to a plant chloroplast or rDNA database and SPAdes [81] for de-novo assembly and iterative extension. During the assembly and iteration process Blast+ [82] is used to identify off-target contigs, which are then removed or trimmed. The resulting plastid genome was then used as a reference for mapping the original reads back using Genieous Prime v2023.2.1 (Biomatters Ltd., Auckland, New Zealand) allowing only mapping of paired reads mapped nearby with a minimum overlap of 75 bp and a minimum overlap identity of 98%. The result was manually examined and corrected where necessary.

The initial annotation of the chloroplast genome was conducted using GeSeq v2.03 [83]. This process included the annotation of chloroplast inverted repeats (IRs), the rps12 interspersed gene, protein-coding sequences, transfer RNAs (tRNAs), and ribosomal RNAs (rRNAs), employing thresholds of 55% identity for protein annotation and 90% for DNA and RNA sequences. Additionally, tRNAscan-SE v2.0.7 [84] and Chloë v0.1.0 (https://github.com/ian-small/chloe) were integrated into GeSeq as supplementary annotation tools. Annotations were then manually refined using Geneious Prime v2023.2.1. To assess the presence of chloroplast genomes with varying orientations of the single-copy regions—large single copy (LSC) and small single copy (SSC)—specific motifs from the border regions of the inverted repeat (IR) and single-copy units were selected. The relative orientation of these motifs was determined on the same reads using seqkit v2.6.1 [85]. The start and end positions of the LSC, SSC, and the start positions of the IRs were taken from the gff3 files for each accession to generate a bed file. The first and last 10–35 bp for each of the regions for each accession were extracted by BEDTools [86]. The orientation of these motifs was inferred from the forward and reverse read files for each accession using seqkit v2.6.1 [85], which locates a given motif and its reversed complement. An R script was used to count the reads having the same and differing orientations of motifs of adjacent regions.

Both inverted repeats were compared to identify differences between them. The length of LSC, SSC, IRs, and sizes of genes flanking the IRs and the two single copy regions (rpS19, ycf1, ndhF) were taken from the annotation and compared. The longest and shortest SSC, LSC, IRs and complete plastid genomes were visualized using the R package genoPlotR [87]. The 35S rDNA was annotated using Blast+ [82]. Sequences were deposited in the European Nucleotide Archive under PRJEB79487.

Sequence alignment and phylogenetic analysis

We aligned the 35S rDNA and the chloroplast datasets separately using mafft v 7.490 [88] and manually edited the alignments. The 35S rDNA sequences of Crocus contain fast evolving several small repeat regions in the two internal transcribed spacers (ITS1 and ITS2) that make it difficult to align those regions among the whole genus unambiguously [51]. These parts were excluded.

We performed phylogenetic analyses for (a) the plastid dataset, (b) the 35S rDNA dataset, and (c) a combined dataset. We used IQ-tree v 2.2.2.6 [89, 90] to perform species relationship estimations in a maximum-likelihood (ML) framework. We used MODELFINDER [91] to find the best-fit model for phylogenetic inference as determined by the Bayesian information criterion for the ML analysis. For the chloroplast, the TVM + F + I + G4 (transversion model with empirical base frequencies, allowing for a proportion of invariable sites discrete Gamma model with 4 rate categories) was inferred as best fitting model. We used 1000 bootstrap replications [92] to infer clade support and a burn-in of 250. Additionally, maximum parsimony (MP) analysis was performed in PAUP* 4.0a169 [93] with tree-bisection-reconnection (TBR) branch swapping, steepest descent not in effect, and initial 5000 random addition sequences (RAS) restricting the search to 50 trees per replicate. The resulting trees were afterward used as starting trees in a search with maxtree set to 10,000 if the limit of 50 trees was reached during a replicate in step one. Bootstrap support values were obtained by 500 bootstrap re-samples using the same settings as before, except that we excluded RAS.

Additional phylogenetic analyses were performed for nuclear single-copy genes with MrBayes 3.2.6 [94]. For Bayesian inference (BI), two times four chains were run for one million generations under the GTR + Gamma + I model of sequence evolution for the ITS data set and 1.5 million generations for the chloroplast marker data set, sampling a tree every 1000 generations. Converging log-likelihoods, potential scale reduction factors for each parameter, and inspection of tabulated model parameters in MrBayes suggested that stationarity had been reached. The first 25% of trees of each run were discarded as burn-in. Two independent runs of BI analysis were performed to confirm that separate analyses converged on the same result. The same topology and similar posterior probabilities (pp) of nodal support resulted in each of the two analyses.

Divergence time estimation

We chose the IQtree consensus tree for a prime analysis in TreePL [95]. The branch lengths were multiplied by 100 using R and the R package ape to avoid too short branches, as suggested by the authors. Initially, a dataset included chloroplast coding regions of Crocus, Romulea and Iris (Iris_domestica_MT001880, Iris_rossii_NC056180, Iris_tectorum_NC056093: [96]; Iris_sanguinea_KT626943: [97] was analyzed with the calibration point for the stem age of Iridaceae 60 +/- 10 Mya stem [98]. In a second calibration, the (a) complete chloroplast alignment as well as the (b) 35S rDNA alignment was used for Crocus and Romulea, setting the split of Crocus from Romulea to 21–25 Mya (Crocoideae). 35S rDNA alignment furthermore included Afrocrocus unifolius and Syringodea bifucata.

Identification of whole genome duplications

Genes identified as duplicates by BUSCO 5.4.6 [99] for our three HiFi assemblies, the C. sativus assembly [100](BIG Data Center accession ID: CRA007742) as well as for RNAseq assemblies of C. hardriaticus (Nemati et al. unpublished) were aligned by macse v2.07 [101]. The alignments were converted into axt format using Unix commands. Synonymous substitutions per synonymous site (Ks) Ks rates were calculated with the KaKs_calculator2.0 [102], accounting for variable mutation rates across sites with a maximum likelihood method.

Autopolyploidy arises from the duplication of a single species' genome. As a result, the paralogs formed within a species are derived from identical genomic backgrounds. In this case, the divergence between the paralogs and the Ks value is directly connected with the WGD. In allopolyploidization, the paralogs in each species originated from two different parental genomes, which were already diverged at the time of polyploidization. Therefore, the divergence between the paralogs and their Ks is related to the time of the last common ancestor of the parental genomes and not to the actual WGD. The Ks from an allopolyploidization can result in an older WGD estimation. To avoid an overestimation, we used the maximum stem age as well as the minimum crown age of Crocus and the mean Ks of all genes probably duplicated at this time, considering the standard deviation of 10% for the Ks. We calculated the general substitution rate (r) for Ks per million years (\(\:r\:=\frac{Ks}{2T\times\:{10}^{-9}}\times\:{10}^{-6};T\:\text{w}\text{a}\text{s}\:\text{t}\text{h}\text{e}\:\text{t}\text{i}\text{m}\text{e}\:\text{d}\text{e}\text{r}\text{i}\text{v}\text{e}\text{d}\:\text{f}\text{r}\text{o}\text{m}\:\text{t}\text{h}\text{e}\:\text{d}\text{a}\text{t}\text{e}\text{d}\:\text{p}\text{h}\text{y}\text{l}\text{o}\text{g}\text{e}\text{n}\text{y}\)). To get evidence, if WGDs were of allo- or autoploid origin, we haplophased and analyzed supposed single-copy genes regions. For that we used the BUSCO results of our three HiFi assemblies (C. longiflorus, C. vernus, C. weldenii) and the assembly of C. sativus as well as read data of Afrocrocus unifolius (ERR7617986), Gladiolus decorates (ERR7617923), Romulea monadelpha (ERR7617928) and Syringodea bifucata (ERR7617929) from Kew Tree of Life Explorer [79]. First, short reads of the latter four species were mapped to the BUSCO genes of crocuses that were found duplicated in at least three out of four species (C. longiflorus, C. sativus, C. vernus, C. weldenii) in Geneious with the Geneious mapper. Mapped reads were de-novo assembled with the Geneious assembler. To generate contigs, reads needed an overlap identity of 100%, allowing only paired reads to be considered. This approach was aimed to provide haplophased contigs in cases where more than one allele or copy of the target gene is present. Contigs with < 10× coverage and < 500 bp length were excluded. Contigs fulfilling these criteria were aligned with the respective genes of C. longiflorus, C. sativus, C. vernus and C. weldenii and phylogenetic relationships were inferred as described above.

Quantification of TEs and satDNA

We obtained repeat consensus sequences by graph-based clustering of Illumina WGS reads [103] from 219 samples (species or accession of uncertain taxonomic rank) representing the whole genus (Supplementary Table 1) using RepeatExplorer2 (RE2) [104] implemented in the Galaxy server [105]. We also employed TAREAN [106] to get a better representation of satDNA. Because there are too many accessions for the comparative analysis feature of RE2 to handle, we only ran TAREAN and RE2 individually for each accession and then performed comparative quantification subsequently. Repeat classification was based on the automatically generated output of RE2 and TAREAN. No further characterization was performed for the repeats classified as ‘unknown’ in this present work but will be a subject in future work.

First, we generated a non-redundant (nr) repeat database from the outputs of TAREAN and RE2 to perform a comparative quantification. To do this, we collected all TAREAN and RE2 output with rank annotation (ranks 1–4) and the contig with the highest genome representation for each cluster, generating 142,149 sequences. We then clustered sequences that shared at least 80% identity to a centroid sequence after sorting the sequence lengths by decreasing length as defined in the UCLUST program [107]. To differentiate cluster output nomenclature from RE2/TAREAN and UCLUST, we refer to “cluster” for RE2/TAREAN output and “ucluster” for UCLUST output. The UCLUST clustering resulted in 41,184 uclusters and reduced to 40,599 after removing uclusters annotated as “organelle,” “mitochondria,” and “plastid” from the RE2/TAREAN cluster annotation. Low-copy clusters (RE2 or TAREAN output above CL300) were filtered out, generating a final list of 16,083 nr uclusters. The centroid sequence of each ucluster was used as a ucluster consensus for downstream analyses. To compare the relative abundance of each repeat in each accession, we mapped back eight million reads from each accession to the 16,083 reference sequences using the Burrows-Wheeler Alignment tool [108]. Accessions with fewer reads for mapping were normalized to eight million to have a comparable quantification of repeat abundance.

To quantify the abundance of each repeat category as defined by RE2 and TAREAN output annotation, we summed up all normalized read counts mapped to the reference repeats with corresponding annotation. We then divided this value by the normalized input reads (eight million) and multiplied by 100 to get the relative genome proportions (GP%). Because not all samples have genome size information, we did not normalize the relative GP to the genome sizes for comparison. For the total repeat content of each sample, we summed up the GP of all annotated repeats.

Satellitome characterization

To isolate satellite repeats, we collected all RE2 and TAREAN output (2,146 repeats) that belonged to rank1 (putative high-confidence satellites, 548 repeats) and rank2 (putative low-confidence satellites, 1598 sequences). These repeats belonged to 999 uclusters (Supplementary Table 4). Because of the rapid dynamics of repeats, some repeats match with other repeats that belong to different uclusters, which share only sections of the repeat monomers. We manually curated the satellite repeats and grouped those from different uclusters when they shared substantial conservation after visual inspection using CLC Main Workbench v24.0 (QIAGEN, Aarhus, Denmark). This analysis resulted in a final list of 564 CroSats. These 564 CroSats belong to 537 uclusters because we separated some uclusters more than one CroSats based on noticeable length differences (Supplementary Table 6).

Quantification of microsatellites

Because genomic repeats are not limited to long monomeric repeats, we also analyzed each accession's total microsatellite content. For this analysis, we extended the definition of “microsatellites” to include 2–30 bp monomers covering an array of ca. 100 bp. We used the same input reads as those in RE2 and TAREAN analyses (Supplementary Table 3). However, for the microsatellite analysis, we concatenated all reads by adding 50 ambiguous ‘N’ sequences and used the resulting single ‘super’ sequence as an input for the Tandem Repeats Finder (TRF) [109]. We scanned for a maximum of 500-bp repeats with TRF and filtered the output using TRAP [110] to include only 30-bp microsatellites. The TRAP analysis generated 357,682 microsatellite sequences.

Because the same microsatellite sequence could be represented multiple times within the same sample in different strand orientations or base order as well as between different samples, we generated a non-redundant list of microsatellites by clustering sequences that share at least 90% identity using the UCLUST program [107]. This clustering generated 35,985 microsatellite uclusters (Supplementary Table 10). Because UCLUST could cluster repeats with different lengths if the shorter repeat matches the longer centroid sequence by a given parameter, the resulting uclusters may include microsatellites with different period (or monomer) sizes. We used the period size mode to represent the period size of each microsatellite ucluster and counted the different period size variants in each microsatellite ucluster. All base counts of each microsatellite from all samples were added to make up the total length of each microsatellite ucluster. The mean array size is the product of the mean repeats per locus and the period size.

Statistical analysis for lineage-specific repeat amplification and co-abundance of TE and microsatellites

To identify repeats that were significantly amplified in each sub-lineage (based on the 35S rDNA phylogenetic tree classification), we performed a correlation analysis using the DESeq2 [111] and WGCNA packages [112] used routinely for weighted gene co-expression network analysis. We analyzed the correlation between the repeat abundance and its amplification in each sub-lineage (as a proxy to phenotype). Repeats present in all accessions were categorized as “core” repeats; otherwise, they were classified as “accessory” repeats. Repeats that showed significant amplification (p < 0.05) in each sub-lineage were isolated, resulting in 8,918 repeats. To get the list of strictly sub-lineage-specific amplified repeats we further filtered out from these isolated repeats those that are not significantly amplified (p ≥ 0.05) in others, resulting in 2,308 repeats. We henceforth refer to the repeat grouping that amplifies in sub-lineages as “Repeat Group.”

To test for co-abundance of TE/satDNA and microsatellites between Romulea and the Crocus subgenera (Carpetanus and Orientalis were combined because they have only a few species each), we performed a weighted repeat co-abundance analysis using DESeq2 [111]. Since DESeq2 was originally designed to analyze differential gene expression between tissues or phenotypes, we used the repeats as a proxy for genes and the genus/main lineage as a proxy for tissues/phenotypes. We combined the 41,184 TE/satellite uclusters from RE2 and TAREAN and the 35,985 microsatellite uclusters as the input dataset. We filtered the DESeq2 output for uclusters included in the 16,083 nr TE/sat database and microsatellites with total bases > 10 kb (Supplementary Tables 6, 10). We considered repeats with greater than two-fold change with p-value < 0.05 between Romulea and the Crocus subgenera significantly differentially amplified. We collected all repeats that were significantly amplified in every three comparisons and categorized them as “subgenus-specific,” “shared,” or “genus-specific” (Supplementary Table 11). However, being subgenus-specific does not mean a repeat is absent in other subgenera. This only means that it is not supported statistically compared to Romulea.

Annotation of microsatellite insertions in LTR retrotransposons

We scanned the genome assemblies of C. longiflorus and C. vernus for LTR retrotransposons using DANTE and DANTE-LTR [113] and extracted full-length LTR retrotransposons with known protein domains, flanking LTRs, target duplication sites, and primer-binding sites (DLTP in DANTE-LTR nomenclature). Because the RE2/TAREAN output are just fragments of LTR retrotransposons, we performed BLAST against the full-length LTR database using the 2,209 microsatellite-co-amplified TE/satDNA sequences as query sequence (Supplementary Table 10) to identify the full-length transposon complements of the query sequences. We then identified microsatellite arrays within the LTR retrotransposons using TRF [109] and performed blastn-short [82] against the co-amplified microsatellites using the LTR-derived microsatellite monomers as the query to match co-amplified LTR retrotransposons and microsatellites. We annotated the full-length LTR retrotransposons that co-amplified with microsatellites using CLC Main Workbench using gff files from DANTE LTR and TRF.

Visualization

Circular plot visualization of repeat abundance and bubble plots of the satellitome and TE-microsatellite co-abundance were generated using iTOL [114]. We used the UpsetR package [115] to generate upset plots to show the amplification of repeats in each sub-lineage. We used the ‘geom_flow’ function in ggplot2 v.3.5.1 [116] to generate the alluvial plot that shows the different categories of the panrepeatome. We used ‘prcomp’ function in base R v.4.2.1 for the principal component analysis and ComplexHeatmap R package v.2.12.1 [117] to visualize the repeat abundance. Volcano plots to show the co-abundance of TE/satDNA and microsatellites were generated using ggplot2 [116] with the DESeq2 output [111]. We generated composite figures using Adobe Illustrator CC v28.7 implemented in Adobe Creative Cloud v5.5.0.617 (Adobe, USA).

Chromosome counting and genome size measurement

Chromosome counts were derived from several sources: existing literature, direct observations made during this study, and extrapolations for some individuals based on genome size data combined with published chromosome counts from nearby populations of the same species. For direct observations, roots were cut about 2 cm from the tips, pretreated with 2mM 8-hydroxyquinoline for 5h at room temperature, and then kept in cold water overnight in a refrigerator. Roots were fixed in Carnoy’s solution (3:1 ethanol acid) for 24 hours and stored in 70% ethanol until needed. Slide preparation followed the protocol of [40] and [118]. For slide processing, samples were fixed in a 2% formaldehyde solution (Sigma-Aldrich, catalog number 47608) for 3 minutes, followed by dehydration through an ethanol series (70%, 90%, and 99%). Chromosomes were stained with 1µg/µL DAPI in 2× SSC. Images were captured using a 100× objective of an Olympus BX61 fluorescence microscope (Olympus).

Genome sizes were measured for 139 individuals. Genome size was determined using propidium iodide (PI) as a stain in flow cytometry with a Cyflow Space (Sysmex Partec) flow cytometer, following essentially the procedure described in [119]. We mainly used rye (Secale cereale; 16.01 pg/2C) or pea (Pisum sativum; 9.09 pg/2C) as internal size standards and the buffer CyStain PI Absolute P (Sysmex Partec). To link the genome sizes with the molecular data, we used silica-gel-dried leaves from the same individual used for DNA extraction whenever possible.

Probe design and FISH analysis

To design oligonucleotide probes (oligoprobes) targeting satellite repeats (CroSat042, CroSat137, and CroSat144), we aligned contigs from the same cluster output from RepeatExplorer [104] and we selected at least two 30-bp regions with the highest conservation. Fluorochrome modifications were added at the 5’ and 3’ ends of the 30-bp sequences (Supplementary Table 13). For the 5S and 35S rDNA, we designed oligoprobes targeting the coding regions (18S coding region for the 35S rDNA) following the methods of Waminal et al. [120]. For the telomeric repeats, we concatenated the monomers to generate 30-bp oligoprobes.

We used the current version of the C. longiflorus genome assembly, which is still at the contig level without Hi-C scaffolding, to design oligoprobes for chromosome painting. We selected the top 20 contigs and designed 45-mer single-copy sequences using Chorus2 [121]. We could already see some contigs representing whole chromosomes or chromosome arms as defined by the location of the centromeric repeats. We selected six targets representing individual chromosome arms. We ordered one pre-labeled oligoprobe library with 9,000 single-copy oligos for each chromosome target from Daicel Arbor Biosciences (Michigan, USA).

For fluorescence in situ hybridization (FISH), the FISH hybridization mix contained 32 µl of master mix (50% formamide, 10% dextran sulfate, and 2 × SSC), ~ 25 ng of each probe, and the total volume was brought to 40 µl with nuclease-free water. For the satellite repeat probes, we incubated the FISH hybridization mix on the slide at 25°C for at least one hour, following the method of Waminal et al. [120]. For the chromosome painting probes, we incubated the FISH hybridization mix at 37°C overnight.

Fluorescence microscopy and image analysis

FISH images were captured using an Olympus BX61 microscope equipped with an ORCA-R2 C10600 CCD camera (Hamamatsu, Japan). FISH signals were acquired in gray scale using appropriate filters for fluorochromes used through the software CellSens Dimension v1.11 (Olympus Soft Imaging Solutions, Germany). Signal pseudo-coloring, color level adjustment, and chromosome cutting to generate the karyograms were performed using Adobe Photoshop CC v25.11 implemented in Adobe Creative Cloud v5.5.0.617 (Adobe, USA). Karyotype ideograms were drawn using Adobe Illustrator CC v28.7.

Author Contributions

Conceptualization: Doerte Harpke, Frank Blattner, and Nomar Espinosa Waminal

Data analysis (chloroplast and rDNA trees, WGD estimation, and genome assembly): Doerte Harpke

Data analysis (panrepeatome, chromosome painting probe design, FISH, and karyotyping: Nomar Espinosa Waminal

Data visualization: Nomar Espinosa Waminal and Doerte Harpke

Writing – original draft: Nomar Espinosa Waminal and Doerte Harpke

Writing – review and editing: Nomar Espinosa Waminal, Doerte Harpke, and Frank Blattner

Funding acquisition: Doerte Harpke and Frank Blattner

All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by grants of the Deutsche Forschungsgemeinschaft (DFG grant HA 7550/2 to Doerte Harpke, and HA 7550/4 to Doerte Harpke and Frank Blattner).

Data Availability Statement

The analyzed DNA sequences and raw reads are available through ENA study accession number PRJEB79487.

Acknowledgments

We would like to thank Ina Faustmann, Christina Koch, Birgit Kraenzlin, for help with plant cultivation and lab work; Axel Himmelbach, Ines Walde and Manuala Knauft for performing PacBio Hifi and Hi-C sequencing; and Almiila Ciftci, Osman Erol, Helmut Kerndorff, Thomas Huber, Lorenzo Peruzzi, Janis Rukšāns, Dirk Schnabel and Lulezim Shuka for providing Crocus materials for molecular analysis. We thank the authorities of the respective countries for supporting our work with collection permits. Computational resources for RepeatExplorer analysis were provided by the ELIXIR-CZ project (LM2023055), part of the international ELIXIR infrastructure.

Additional information

Correspondence should be addressed to Doerte Harpke ([email protected]) and Nomar Espinosa Waminal ([email protected])

Escudero, M., S. Martín-Bravo, I. Mayrose, et al., Karyotypic changes through dysploidy persist longer over evolutionary time than polyploid changes. PLoS ONE, 2014. 9(1): p. e85266. DOI: https://doi.org/10.1371/journal.pone.0085266
Mandakova, T. and M.A. Lysak, Post-polyploid diploidization and diversification through dysploid changes. Current Opinion in Plant Biology, 2018. 42: p. 55-65. DOI: https://doi.org/10.1016/j.pbi.2018.03.001
Roalson, E., A. McCubbin, and R. Whitkus, Chromosome evolution in Cyperales. Aliso, 2007. 23(1): p. 62-71. DOI: https://doi.org/10.5642/aliso.20072301.08
Zhan, S.H., S.P. Otto, and M.S. Barker, Broad variation in rates of polyploidy and dysploidy across flowering plants is correlated with lineage diversification. 2021, Cold Spring Harbor Laboratory. https://doi.org/10.1101/2021.03.30.436382
Comai, L., The advantages and disadvantages of being polyploid. Nature Reviews Genetics, 2005. 6(11): p. 836-846. DOI: https://doi.org/10.1038/nrg1711
Fedoroff, N.V., Transposable elements, epigenetics, and genome evolution. Science, 2012. 338(6108): p. 758-767. DOI: https://doi.org/10.1126/science.338.6108.758
Glombik, M., V. Bačovský, R. Hobza, et al., Competition of parental genomes in plant hybrids. Frontiers in Plant Science, 2020. 11(200). DOI: https://doi.org/10.3389/fpls.2020.00200
Dodsworth, S., M.W. Chase, and A.R. Leitch, Is post-polyploidization diploidization the key to the evolutionary success of angiosperms? Botanical Journal of the Linnean Society, 2016. 180(1): p. 1-5. DOI: https://doi.org/10.1111/boj.12357
Soltis, P.S. and D.E. Soltis, Ancient WGD events as drivers of key innovations in angiosperms. Current Opinion in Plant Biology, 2016. 30(Supplement C): p. 159-165. DOI: https://doi.org/10.1016/j.pbi.2016.03.015
Ramsey, J., Polyploidy and ecological adaptation in wild yarrow. Proceedings of the National Academy of Sciences USA, 2011. 108(17): p. 7096-7101. DOI: https://doi.org/10.1073/pnas.1016631108
Winterfeld, G., A. Ley, M.H. Hoffmann, et al., Dysploidy and polyploidy trigger strong variation of chromosome numbers in the prayer-plant family (Marantaceae). Plant Systematics and Evolution, 2020. 306(2). DOI: https://doi.org/10.1007/s00606-020-01663-x
Ramsey, J. and D.W. Schemske, Pathways, mechanisms, and rates of polyploid formation. Annual Review of Ecology, Evolution, and Systematics, 1998. 29(Volume 29, 1998): p. 467-501. DOI: https://doi.org/10.1146/annurev.ecolsys.29.1.467
Josefsson, C., B. Dilkes, and L. Comai, Parent-dependent loss of gene silencing during interspecies hybridization. Current Biology, 2006. 16(13): p. 1322-1328. DOI: https://doi.org/10.1016/j.cub.2006.05.045
Bird, K.A., R. Vanburen, J.R. Puzey, et al., The causes and consequences of subgenome dominance in hybrids and recent polyploids. New Phytologist, 2018. 220(1): p. 87-93. DOI: https://doi.org/10.1111/nph.15256
Belyayev, A., Bursts of transposable elements as an evolutionary driving force. Journal of Evolutionary Biology, 2014. 27(12): p. 2573-2584. DOI: https://doi.org/10.1111/jeb.12513
Ruiz-Ruano, F.J., M.D. López-León, J. Cabrero, et al., High-throughput analysis of the satellitome illuminates satellite DNA evolution. Scientific Reports, 2016. 6(1): p. 28333. DOI: https://doi.org/10.1038/srep28333
Mestrović, N., M. Plohl, B. Mravinac, et al., Evolution of satellite DNAs from the genus Palorus--experimental evidence for the "library" hypothesis. Molecular Biology and Evolution, 1998. 15(8): p. 1062-1068. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a026005
Fry, K. and W. Salser, Nucleotide sequences of HS-α satellite DNA from kangaroo rat dipodomys ordii and characterization of similar sequences in other rodents. Cell, 1977. 12(4): p. 1069-1084. DOI: https://doi.org/https://doi.org/10.1016/0092-8674(77)90170-2
Argueso, J.L., J. Westmoreland, P.A. Mieczkowski, et al., Double-strand breaks associated with repetitive DNA can reshape the genome. Proceedings of the National Academy of Sciences USA, 2008. 105(33): p. 11845-11850. DOI: https://doi.org/10.1073/pnas.0804529105
Schubert, I. and M.A. Lysak, Interpretation of karyotype evolution should consider chromosome structural constraints. Trends in Genetics, 2011. 27(6): p. 207-216. DOI: https://doi.org/10.1016/j.tig.2011.03.004
Huang, K. and L.H. Rieseberg, Frequency, origins, and evolutionary role of chromosomal inversions in plants. Frontiers in Plant Science, 2020. 11. DOI: https://doi.org/10.3389/fpls.2020.00296
Hofstatter, P.G., G. Thangavel, T. Lux, et al., Repeat-based holocentromeres influence genome architecture and karyotype evolution. Cell, 2022. DOI: https://doi.org/10.1016/j.cell.2022.06.045
Murat, F., A. Armero, C. Pont, et al., Reconstructing the genome of the most recent common ancestor of flowering plants. Nature Genetics, 2017. 49(4): p. 490-496. DOI: https://doi.org/10.1038/ng.3813
Murat, F., J.-H. Xu, E. Tannier, et al., Ancestral grass karyotype reconstruction unravels new mechanisms of genome shuffling as a source of plant evolution. Genome Research, 2010. 20(11): p. 1545-1557. DOI: https://doi.org/10.1101/gr.109744.110
Heckmann, S., E. Schroeder-Reiter, K. Kumke, et al., Holocentric chromosomes of Luzula elegans are characterized by a longitudinal centromere groove, chromosome bending, and a terminal nucleolus organizer region. Cytogenetic and Genome Research, 2011. 134(3): p. 220-228. DOI: https://doi.org/10.1159/000327713
Ribeiro, T., C.E. Buddenhagen, W.W. Thomas, et al., Are holocentrics doomed to change? Limited chromosome number variation in Rhynchospora Vahl (Cyperaceae). Protoplasma, 2018. 255(1): p. 263-272. DOI: https://doi.org/10.1007/s00709-017-1154-4
Zedek, F. and P. Bureš, Holocentric chromosomes: from tolerance to fragmentation to colonization of the land. Annals of Botany, 2017. 121(1): p. 9-16. DOI: https://doi.org/10.1093/aob/mcx118
Rice, A., L. Glick, S. Abadi, et al., The Chromosome Counts Database ( CCDB ) – a community resource of plant chromosome numbers. New Phytologist, 2015. 206(1): p. 19-26. DOI: https://doi.org/10.1111/nph.13191
Román‐Palacios, C., C.A. Medina, S.H. Zhan, et al., Animal chromosome counts reveal a similar range of chromosome numbers but with less polyploidy in animals compared to flowering plants. Journal of Evolutionary Biology, 2021. 34(8): p. 1333-1339. DOI: https://doi.org/10.1111/jeb.13884
Cui, L., P.K. Wall, J.H. Leebens-Mack, et al., Widespread genome duplications throughout the history of flowering plants. Genome Research, 2006. 16. DOI: https://doi.org/10.1101/gr.4825606
Soltis, D.E., V.A. Albert, J. Leebens‐Mack, et al., Polyploidy and angiosperm diversification. American Journal of Botany, 2009. 96(1): p. 336-348. DOI: https://doi.org/10.3732/ajb.0800079
Carta, A., G. Bedini, and L. Peruzzi, A deep dive into the ancestral chromosome number and genome size of flowering plants. New Phytologist, 2020. 228(3): p. 1097-1106. DOI: https://doi.org/10.1111/nph.16668
Guerrero, R.F. and M. Kirkpatrick, Local adaptation and the evolution of chromosome fusions. Evolution, 2014. 68(10): p. 2747-2756. DOI: https://doi.org/10.1111/evo.12481
Durmaz, E., E. Kerdaffrec, G. Katsianis, et al., How selection acts on chromosomal inversions, in eLS. 2020. p. 307-315.
Schaal, S.M., B.C. Haller, and K.E. Lotterhos, Inversion invasions: when the genetic basis of local adaptation is concentrated within inversions in the face of gene flow. Philosophical Transactions of the Royal Society B: Biological Sciences, 2022. 377(1856). DOI: https://doi.org/10.1098/rstb.2021.0200
De Vos, J.M., H. Augustijnen, L. Bätscher, et al., Speciation through chromosomal fusion and fission in Lepidoptera. Philosophical Transactions of the Royal Society B: Biological Sciences, 2020. 375(1806): p. 20190539. DOI: https://doi.org/10.1098/rstb.2019.0539
Yoshida, K., C. Rödelsperger, W. Röseler, et al., Chromosome fusions repatterned recombination rate and facilitated reproductive isolation during Pristionchus nematode speciation. Nature Ecology & Evolution, 2023. 7(3): p. 424-439. DOI: https://doi.org/10.1038/s41559-022-01980-z
Jayakodi, M., S. Padmarasu, G. Haberer, et al., The barley pan-genome reveals the hidden legacy of mutation breeding. Nature, 2020. 588(7837): p. 284-289. DOI: https://doi.org/10.1038/s41586-020-2947-8
Noor, M.A.F., K.L. Grams, L.A. Bertucci, et al., Chromosomal inversions and the reproductive isolation of species. Proceedings of the National Academy of Sciences USA, 2001. 98(21): p. 12084-12088. DOI: https://doi.org/doi:10.1073/pnas.221274498
Waminal, N.E., R.J. Pellerin, S.-H. Kang, et al., Chromosomal mapping of tandem repeats revealed massive chromosomal rearrangements and insights into Senna tora dysploidy. Frontiers in Plant Science, 2021. 12(154): p. 629898. DOI: https://doi.org/10.3389/fpls.2021.629898
Brighton, C.A., Cytological problems in the genus Crocus (Iridaceae): I. Crocus vernus aggregate. Kew Bulletin, 1976. 31(1): p. 33. DOI: https://doi.org/10.2307/4108994
Brighton, C.A., Cytology of Crocus sativus and its allies (Iridaceae). Plant Systematics and Evolution, 1977. 128(3-4): p. 137-157. DOI: https://doi.org/10.1007/bf00984551
Brighton, C.A., Cytological problems in the genus Crocus (Iridaceae): II. Crocus cancellatus aggregate. Kew Bulletin, 1977. 32(1): p. 33. DOI: https://doi.org/10.2307/4117257
Brighton, C.A., B. Mathew, and P. Rudall, A detailed study of Crocus speciosus and its ally C. pulchellus (Iridaceae). Plant Systematics and Evolution, 1983. 142(3-4): p. 187-206. DOI: https://doi.org/10.1007/bf00985898
Raca, I., F.R. Blattner, N.E. Waminal, et al., Disentangling Crocus series Verni and its polyploids. Biology, 2023. 12(2): p. 303. DOI: https://doi.org/10.3390/biology12020303
Kuznetsov, D., F. Tegenfeldt, M. Manni, et al., OrthoDB v11: annotation of orthologs in the widest sampling of organismal diversity. Nucleic Acids Research, 2022. 51(D1): p. D445-D451. DOI: https://doi.org/10.1093/nar/gkac998
Hu, J., Y. Zheng, and X. Shang, MiteFinderII: a novel tool to identify miniature inverted-repeat transposable elements hidden in eukaryotic genomes. BMC Medical Genomics, 2018. 11(S5). DOI: https://doi.org/10.1186/s12920-018-0418-y
Witte, C.-P., Q.H. Le, T. Bureau, et al., Terminal-repeat retrotransposons in miniature (TRIM) are involved in restructuring plant genomes. Proceedings of the National Academy of Sciences, 2001. 98(24): p. 13778-13783. DOI: https://doi.org/10.1073/pnas.241341898
Kalendar, R., J. Tanskanen, W. Chang, et al., Cassandra retrotransposons carry independently transcribed 5S RNA. Proceedings of the National Academy of Sciences, 2008. 105(15): p. 5833-5838. DOI: https://doi.org/10.1073/pnas.0709698105
Deragon, J.-M. and X. Zhang, Short Interspersed Elements (SINEs) in plants: Origin, classification, and use as phylogenetic markers. Systematic Biology, 2006. 55(6): p. 949-956. DOI: https://doi.org/10.1080/10635150601047843
Harpke, D., S. Meng, T. Rutten, et al., Phylogeny of Crocus (Iridaceae) based on one chloroplast and two nuclear loci: ancient hybridization and chromosome number evolution. Molecular Phylogenetics and Evolution, 2013. 66(3): p. 617-627. DOI: https://doi.org/10.1016/j.ympev.2012.10.007
Petersen, G., O. Seberg, S. Thorsøe, et al., A phylogeny of the genus Crocus (Iridaceae) based on sequence data from five plastid regions. TAXON, 2008. 57(2): p. 487-499. DOI: https://doi.org/https://doi.org/10.2307/25066017
Schranz, M.E., S. Mohammadin, and P.P. Edger, Ancient whole genome duplications, novelty and diversification: the WGD Radiation Lag-Time Model. Current Opinion in Plant Biology, 2012. 15(2): p. 147-153. DOI: https://doi.org/10.1016/j.pbi.2012.03.011
Walden, N., D.A. German, E.M. Wolf, et al., Nested whole-genome duplications coincide with diversification and high morphological disparity in Brassicaceae. Nature Communications, 2020. 11(1). DOI: https://doi.org/10.1038/s41467-020-17605-7
Tank, D.C., J.M. Eastman, M.W. Pennell, et al., Nested radiations and the pulse of angiosperm diversification: increased diversification rates often follow whole genome duplications. New Phytologist, 2015. 207(2): p. 454-467. DOI: https://doi.org/https://doi.org/10.1111/nph.13491
Jakob, S.S. and F.R. Blattner, A chloroplast genealogy of Hordeum (Poaceae): Long-term persisting haplotypes, incomplete lineage sorting, regional extinction, and the consequences for phylogenetic inference. Molecular Biology and Evolution, 2006. 23(8): p. 1602-1612. DOI: https://doi.org/10.1093/molbev/msl018
Durand, E.Y., N. Patterson, D. Reich, et al., Testing for ancient admixture between closely related populations. Molecular Biology and Evolution, 2011. 28(8): p. 2239-2252. DOI: https://doi.org/10.1093/molbev/msr048
McCann, J., J. Macas, P. Novák, et al., Differential genome size and repetitive DNA evolution in diploid species of Melampodium sect. Melampodium (Asteraceae). Frontiers in Plant Science, 2020. 11. DOI: https://doi.org/10.3389/fpls.2020.00362
Breman, F.C., G. Chen, R.C. Snijder, et al., Repeatome-based phylogenetics in Pelargonium section Ciconium (sweet) Harvey. Genome Biology and Evolution, 2021. 13(12). DOI: https://doi.org/10.1093/gbe/evab269
Dodsworth, S., M.W. Chase, L.J. Kelly, et al., Genomic repeat abundances contain phylogenetic signal. Systematic Biology, 2015. 64(1): p. 112-126. DOI: https://doi.org/10.1093/sysbio/syu080
Schmidt, N., K. Sielemann, S. Breitenbach, et al., Repeat turnover meets stable chromosomes: repetitive DNA sequences mark speciation and gene pool boundaries in sugar beet and wild beets. The Plant Journal, 2023. DOI: https://doi.org/10.1111/tpj.16599
Mata-Sucre, Y., W. Matzenauer, N. Castro, et al., Repeat-based phylogenomics shed light on unclear relationships in the monocentric genus Juncus L. (Juncaceae). Molecular Phylogenetics and Evolution, 2023. 189: p. 107930. DOI: https://doi.org/https://doi.org/10.1016/j.ympev.2023.107930
Zhou, B.-F., S. Yuan, A.A. Crowl, et al., Phylogenomic analyses highlight innovation and introgression in the continental radiations of Fagaceae across the Northern Hemisphere. Nature Communications, 2022. 13(1). DOI: https://doi.org/10.1038/s41467-022-28917-1
Buerkle, C.A., R.J. Morris, M.A. Asmussen, et al., The likelihood of homoploid hybrid speciation. Heredity, 2000. 84(4): p. 441-451. DOI: https://doi.org/10.1046/j.1365-2540.2000.00680.x
Zhang, C.-Z. and D. Pellman, Chromosome breakage-replication/fusion enables rapid DNA amplification. bioRxiv, 2024: p. 2024.08.17.608415. DOI: https://doi.org/10.1101/2024.08.17.608415
Kockler, Z.W., B. Osia, R. Lee, et al., Repair of DNA breaks by break-induced replication. Annual Review of Biochemistry, 2021. 90(Volume 90, 2021): p. 165-191. DOI: https://doi.org/https://doi.org/10.1146/annurev-biochem-081420-095551
Peska, V. and S. Garcia, Origin, diversity, and evolution of telomere sequences in plants. Frontiers in Plant Science, 2020. 11(117). DOI: https://doi.org/10.3389/fpls.2020.00117
Van de Peer, Y., T.-L. Ashman, P.S. Soltis, et al., Polyploidy: an evolutionary and ecological force in stressful times. The Plant Cell, 2020. 33(1): p. 11-26. DOI: https://doi.org/10.1093/plcell/koaa015
Tubbs, A., S. Sridharan, N. Van Wietmarschen, et al., Dual roles of Poly(dA:dT) tracts in replication initiation and fork collapse. Cell, 2018. 174(5): p. 1127-1142.e19. DOI: https://doi.org/10.1016/j.cell.2018.07.011
Meyer, M. and M. Kircher, Illumina sequencing library preparation for highly multiplexed target capture and sequencing. Cold Spring Harbor Protocols, 2010. 2010(6): p. pdb.prot5448. DOI: https://doi.org/10.1101/pdb.prot5448
Dvorak, J., P.E. McGuire, and B. Cassidy, Apparent sources of the A genomes of wheats inferred from polymorphism in abundance and restriction fragment length of repeated nucleotide sequences. Genome, 1988. 30(5): p. 680-689. DOI: https://doi.org/10.1139/g88-115
Cheng, H., G.T. Concepcion, X. Feng, et al., Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nature Methods, 2021. 18(2): p. 170-175. DOI: https://doi.org/10.1038/s41592-020-01056-5
Cheng, H., E.D. Jarvis, O. Fedrigo, et al., Haplotype-resolved assembly of diploid genomes without parental data. Nature Biotechnology, 2022. 40(9): p. 1332-1335. DOI: https://doi.org/10.1038/s41587-022-01261-x
Zhou, C., S.A. McCarthy, and R. Durbin, YaHS: yet another Hi-C scaffolding tool. Bioinformatics, 2022. 39(1). DOI: https://doi.org/10.1093/bioinformatics/btac808
Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. 2010; Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Martin, M., Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 2011. 17(1): p. 10. DOI: https://doi.org/10.14806/ej.17.1.200
BBMap, B.B. BBMap. Available from: http://sourceforge.net/projects/bbmap/.
Jin, J.-J., W.-B. Yu, J.-B. Yang, et al., GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes. Genome Biology, 2020. 21(1). DOI: https://doi.org/10.1186/s13059-020-02154-5
Baker, W.J., P. Bailey, V. Barber, et al., A comprehensive phylogenomic platform for exploring the angiosperm tree of life. Systematic Biology, 2021. 71(2): p. 301-319. DOI: https://doi.org/10.1093/sysbio/syab035
Langmead, B. and S.L. Salzberg, Fast gapped-read alignment with Bowtie 2. Nature Methods, 2012. 9(4): p. 357-359. DOI: https://doi.org/10.1038/nmeth.1923
Bankevich, A., S. Nurk, D. Antipov, et al., SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. Journal of Computational Biology, 2012. 19(5): p. 455-477. DOI: https://doi.org/10.1089/cmb.2012.0021
Camacho, C., G. Coulouris, V. Avagyan, et al., BLAST+: architecture and applications. BMC Bioinformatics, 2009. 10(1): p. 421. DOI: https://doi.org/10.1186/1471-2105-10-421
Tillich, M., P. Lehwark, T. Pellizzer, et al., GeSeq – versatile and accurate annotation of organelle genomes. Nucleic Acids Research, 2017. 45(W1): p. W6-W11. DOI: https://doi.org/10.1093/nar/gkx391
Chan, Patricia P., Brian Y. Lin, Allysia J. Mak, et al., tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes. Nucleic Acids Research, 2021. 49(16): p. 9077-9096. DOI: https://doi.org/10.1093/nar/gkab688
Shen, W., S. Le, Y. Li, et al., Seqkit: A cross-platform and ultrafast toolkit for fasta/q file manipulation. PLoS ONE, 2016. 11(10): p. e0163962. DOI: https://doi.org/10.1371/journal.pone.0163962
Quinlan, A.R. and I.M. Hall, BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 2010. 26(6): p. 841-842. DOI: https://doi.org/10.1093/bioinformatics/btq033
Guy, L., J. Roat Kultima, and S.G.E. Andersson, genoPlotR: comparative gene and genome visualization in R. Bioinformatics, 2010. 26(18): p. 2334-2335. DOI: https://doi.org/10.1093/bioinformatics/btq413
Katoh, K. and D.M. Standley, MAFFT Multiple Sequence Alignment Software Version 7: Improvements in performance and usability. Molecular Biology and Evolution, 2013. 30(4): p. 772-780. DOI: https://doi.org/10.1093/molbev/mst010
Minh, B.Q., H.A. Schmidt, O. Chernomor, et al., IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Molecular Biology and Evolution, 2020. 37(5): p. 1530-1534. DOI: https://doi.org/10.1093/molbev/msaa015
Nguyen, L.-T., H.A. Schmidt, A. von Haeseler, et al., IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution, 2014. 32(1): p. 268-274. DOI: https://doi.org/10.1093/molbev/msu300
Kalyaanamoorthy, S., B.Q. Minh, T.K.F. Wong, et al., ModelFinder: fast model selection for accurate phylogenetic estimates. Nature Methods, 2017. 14(6): p. 587-589. DOI: https://doi.org/10.1038/nmeth.4285
Hoang, D.T., O. Chernomor, A. von Haeseler, et al., UFBoot2: Improving the ultrafast bootstrap approximation. Molecular Biology and Evolution, 2017. 35(2): p. 518-522. DOI: https://doi.org/10.1093/molbev/msx281
Swofford, D.L., Phylogenetic analysis using parsimony (* and other methods). 2002, Version.
Ronquist, F., M. Teslenko, P. van der Mark, et al., MrBayes 3.2: Efficient bayesian phylogenetic inference and model choice across a large model space. Systematic Biology, 2012. 61(3): p. 539-542. DOI: https://doi.org/10.1093/sysbio/sys029
Smith, S.A. and B.C. O’Meara, treePL: divergence time estimation using penalized likelihood for large phylogenies. Bioinformatics, 2012. 28(20): p. 2689-2690. DOI: https://doi.org/10.1093/bioinformatics/bts492
Kang, Y.J., S. Kim, J. Lee, et al., Identification of plastid genomic regions inferring species identity from de novo plastid genome assembly of 14 Korean-native Iris species (Iridaceae). PLoS ONE, 2020. 15(10): p. e0241178. DOI: https://doi.org/10.1371/journal.pone.0241178
Lee, H.-j., G.-H. Nam, K. Kim, et al., The complete chloroplast genome sequences of Iris sanguinea Donn ex Hornem. Mitochondrial DNA Part A, 2017. 28(1): p. 15-16. DOI: https://doi.org/10.3109/19401736.2015.1106521
Li, H.-T., T.-S. Yi, L.-M. Gao, et al., Origin of angiosperms and the puzzle of the Jurassic gap. Nature Plants, 2019. 5(5): p. 461-470. DOI: https://doi.org/10.1038/s41477-019-0421-0
Simão, F.A., R.M. Waterhouse, P. Ioannidis, et al., BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics, 2015. 31(19): p. 3210-3212. DOI: https://doi.org/10.1093/bioinformatics/btv351 %J Bioinformatics
Xu, Z., S. Chen, Y. Wang, et al., Crocus genome reveals the evolutionary origin of crocin biosynthesis. Acta Pharmaceutica Sinica B, 2024. 14(4): p. 1878-1891. DOI: https://doi.org/https://doi.org/10.1016/j.apsb.2023.12.013
Ranwez, V., S. Harispe, F. Delsuc, et al., MACSE: Multiple Alignment of Coding SEquences accounting for frameshifts and stop codons. PLoS ONE, 2011. 6(9): p. e22594. DOI: https://doi.org/10.1371/journal.pone.0022594
Wang, D., Y. Zhang, Z. Zhang, et al., KaKs_Calculator 2.0: A toolkit incorporating gamma-series methods and sliding window strategies. Genomics, Proteomics & Bioinformatics, 2010. 8(1): p. 77-80. DOI: https://doi.org/10.1016/S1672-0229(10)60008-3
Novák, P., P. Neumann, and J. Macas, Graph-based clustering and characterization of repetitive sequences in next-generation sequencing data. BMC Bioinformatics, 2010. 11(1): p. 378. DOI: https://doi.org/10.1186/1471-2105-11-378
Novák, P., P. Neumann, and J. Macas, Global analysis of repetitive DNA from unassembled sequence reads using RepeatExplorer2. Nature Protocols, 2020. 15(11): p. 3745-3776. DOI: https://doi.org/10.1038/s41596-020-0400-y
Community, T.G., The Galaxy platform for accessible, reproducible, and collaborative data analyses: 2024 update. Nucleic Acids Research, 2024. 52(W1): p. W83-W94. DOI: https://doi.org/10.1093/nar/gkae410
Novák, P., L.A. Robledillo, A. Koblızková, et al., TAREAN: a computational tool for identification and characterization of satellite DNA from unassembled short reads. Nucleic Acids Research, 2017. 45(12): p. e111.
Edgar, R.C., Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 2010. 26(19): p. 2460-2461. DOI: https://doi.org/10.1093/bioinformatics/btq461
Li, H. and R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics (Oxford, England), 2009. 25(14): p. 1754-1760. DOI: https://doi.org/10.1093/bioinformatics/btp324
Benson, G., Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res, 1999. 27(2): p. 573-80. DOI: https://doi.org/10.1093/nar/27.2.573
Sobreira, T.J., A.M. Durham, and A. Gruber, TRAP: automated classification, quantification and annotation of tandemly repeated sequences. Bioinformatics, 2006. 22(3): p. 361-2. DOI: https://doi.org/10.1093/bioinformatics/bti809
Love, M.I., W. Huber, and S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 2014. 15(12). DOI: https://doi.org/10.1186/s13059-014-0550-8
Langfelder, P. and S. Horvath, WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008. 9(1): p. 559. DOI: https://doi.org/10.1186/1471-2105-9-559
Novák, P., N. Hoštáková, P. Neumann, et al., DANTE and DANTE_LTR: Lineage-centric annotation pipelines for long terminal repeat retrotransposons in plant genomes. 2024, Cold Spring Harbor Laboratory. https://doi.org/10.1101/2024.04.17.589915
Letunic, I. and P. Bork, Interactive Tree Of Life (iTOL): an online tool for phylogenetic tree display and annotation. Bioinformatics, 2007. 23. DOI: https://doi.org/10.1093/bioinformatics/btl529
Conway, J.R., A. Lex, and N. Gehlenborg, UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics, 2017. 33(18): p. 2938-2940. DOI: https://doi.org/10.1093/bioinformatics/btx364
Wickham, H., ggplot2: Elegant Graphics for Data Analysis. Create elegant data visualisations using the grammar of graphics. Version. 2016, New York: Springer-Verlag. 1-189.
Gu, Z., Complex heatmap visualization. iMeta, 2022. 1(3). DOI: https://doi.org/10.1002/imt2.43
Rodríguez-Domínguez, J.M., L.L. Ríos-Lara, E. Tapia-Campos, et al., An improved technique for obtaining well-spread metaphases from plants with numerous large chromosomes. Biotechnic & Histochemistry, 2017. 92(3): p. 159-166. DOI: https://doi.org/10.1080/10520295.2017.1288927
Jakob, S.S., A. Meister, and F.R. Blattner, The considerable genome size variation of Hordeum species (Poaceae) is linked to phylogeny, life form, ecology, and speciation rates. Molecular Biology and Evolution, 2004. 21(5): p. 860-869. DOI: https://doi.org/10.1093/molbev/msh092
Waminal, N.E., R.J. Pellerin, N.-S. Kim, et al., Rapid and efficient FISH using pre-labeled oligomer probes. Scientific Reports, 2018. 8(1): p. 8224. DOI: https://doi.org/10.1038/s41598-018-26667-z
Zhang, T., G. Liu, H. Zhao, et al., Chorus2: design of genome-scale oligonucleotide-based probes for fluorescence in situ hybridization. Plant Biotechnology Journal, 2021. 19(10): p. 1967-1978. DOI: https://doi.org/https://doi.org/10.1111/pbi.13610

The authors declare no competing interests.

3CrocusdysploidyExtendedData.pdf
Supplementary Figures
4Crocusdysploidysuptables.xlsx
Supplementary Tables

Download PDF

Version 1

posted

You are reading this latest preprint version

The Crocus panrepeatome reveals the links between whole-genome duplications, repeat bursts, and descending dysploidy

Status:

Version 1

Abstract

Figures

Introduction

Results

Repeat bursts after WGD events

Microsatellite repeats co-amplified with transposable elements

Amplification of CroSat042 suggests a role in descending dysploidy

Discussion

Summary

Materials and Methods

Plant materials

DNA extraction and sequencing

Genome assembly

Chloroplast genome and 35S rDNA assembly

Sequence alignment and phylogenetic analysis

Divergence time estimation

Identification of whole genome duplications

Quantification of TEs and satDNA

Satellitome characterization

Quantification of microsatellites

Statistical analysis for lineage-specific repeat amplification and co-abundance of TE and microsatellites

Annotation of microsatellite insertions in LTR retrotransposons

Visualization

Chromosome counting and genome size measurement

Probe design and FISH analysis

Fluorescence microscopy and image analysis

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1