Genomes of WT05 and Hale were downloaded from the projects of PRJNA838012 (Lu et al., 2022) and PRJNA16585 (Chan et al., 2010) in the National Center for Biotechnology Information (NCBI) database, respectively. The Rc039 genome was downloaded from the Oil Plants database (http://oilplants.iflora.cn). As available at scaffold ordering, Hale genome was assembled into the chromosomal level by RagTag (Alonge et al., 2019; Alonge et al., 2022) using Rc039 and WT05 as the reference genome, respectively. The Hale genome assembly was further evaluated using QUAST (version 5.2.0) (Gurevich et al., 2013) and BUSCO (version 5.7.1) (Simão et al., 2015) to assess its quality and completeness. Standard metrics such as N50, L50, number of contigs/scaffolds, largest contig/scaffold, GC content, and BUSCO scores were calculated to provide a comprehensive assessment of the assembly. These results are summarized in Supplementary Table S1 . Additionally, the completely assembled Hale genome has been compressed and made available as Supplementary File 2 .

To assess the genomic collinearity and identify paralogous genes between the Rc039 and WT05 assemblies of castor bean, we performed all-against-all best reciprocal BLASTP searches. Paralogous genes were identified, and collinear blocks across the genomes were constructed using MCScanX (https://github.com/wyp1125/MCScanX) (Wang et al., 2012). Each block contains at least five paralogous genes. The original sequences of Chr2–5, 7, 9, and 10 of WT05 were reversely complemented. The chromosome cohort of WT05 was reordered according to the collinearity with Rc039. SubPhaser (version 1.1) (Jia et al., 2022) was used to cluster the 15-mer sequences in the castor bean genome. This analysis focuses on paralogous genes because WT05 and Rc039 are assemblies from the same species, and therefore the comparison aims to reveal gene duplication events, genomic rearrangements, and structural variations.

Orthologous genes among Mercurialis annua, Manihot esculenta, and castor bean were identified using OrthoFinder2 (version 2.2.7) (Emms and Kelly, 2019) with the parameter -S diamond. Subsequently, all single-copy orthologs were subjected to multiple sequence alignment using MAFFT (version 7.407) (Katoh and Standley, 2013), and the poorly conserved blocks were trimmed using trimAl (Capella-Gutiérrez et al., 2009) with default parameters. Finally, the consensus sequence was merged into a supergene. The phylogenetic tree was constructed using RAxML (version 8.1.2) (Stamatakis, 2014) with 100 bootstrap replicates and PROTGAMMAAUTO model and visualized using iTOL (Letunic and Bork, 2021). Protein sequences of Mercurialis annua and Manihot esculenta were downloaded from the projects of PRJEB52246 and PRJNA234389 in the National Center for Biotechnology Information (NCBI) database, respectively.

We built a TE library of the castor bean genome using RepeatModeler (version 2.0.3) (Flynn et al., 2020) in a de novo approach. For this purpose, we first integrated the genome assemblies of Hale, Rc039, and WT05 into a single comprehensive assembly. This integrated assembly provided a robust basis for TE identification across different castor bean accessions. Our initial step involved running RepeatModeler on the integrated genome assembly, which generated 1,742 consensus sequences. We meticulously followed the procedural guidelines outlined in the relevant literature (Goubert et al., 2022) to scrutinize and identify the TE library generated by RepeatModeler. Our initial step involved the compilation of a prioritized list of candidates. The second step is the manual curation process. The detailed steps for the manual curation of a family are described as follows: (A) Selection of query sequence: Begin by selecting a query sequence to initiate the manual curation process; (B) BLAST search: Utilize the putative TE sequence as a query in a BLAST search against the reference genome. Record the obtained hits; (C) Extension of hit coordinates: As the prospective TE families often represent truncated versions of the actual TE, extend the hit coordinates in the genome by a specified number of bases both upstream and downstream. This ensures capturing as much TE sequence as possible; (D) Genomic sequence extraction: Extract the genomic sequences corresponding to the extended coordinates using “bedtools getfasta” (Quinlan and Hall, 2010). Subsequently, generate a multiple sequence alignment (MSA) and save it to a file; (E) Manual curation and consensus generation: Visualize the MSA in an alignment viewer and perform manual curation to generate a consensus sequence. Utilize the “cons” function from the EMBOSS (Rice et al., 2000) package to produce a consensus sequence from the MSA.

We employed TEClass (version 2.0.3) (Bickmann et al., 2023) and DeepTE (Yan et al., 2020) tools to identify TEs in the “unknown” sequences classified by RepeatModeler. Although TEClass is not yet published, we chose this tool for its promising performance in preliminary tests, which showed high accuracy in identifying various TE types through its efficient machine learning algorithm. To ensure the reliability of our results, we also conducted cross-confirmation using DeepTE. For TEClass, we set a threshold of 0.7 to increase the classification confidence. The parameters used for DeepTE were “-m P -sp P” for plant TE models, as described in its respective study. Our analysis of the RepeatModeler results revealed 353 “unknown” sequences. After reclassification using DeepTE, 266 of these sequences were identified as other TE families. Using TEClass2.0, 134 sequences were reclassified. Comparing the results of both tools, we found 103 sequences reclassified by both, with 70 sequences identified as the same TE class by both tools ( Supplementary Table S2 ). The results are tabulated in Supplementary Table S3 . A total of 134 TEs of the unknown family were re-identified, culminating in 462 consensus sequences ( Supplementary Table S4 ).

We followed the naming conventions suggested in the literature (Goubert et al., 2022) for renaming TEs. One method, similar to that used by the TE repository Repbase (Bao et al., 2015), is to use the format “superfamilyX_yYyy”, where X is a unique number and yYyy is a four-letter identifier for the species in question—for example, two Ty3-retrotransposons elements from the genome of the castor bean Ricinus communis could be named Ty3RT-1_rCom and Ty3RT-2_ rCom, respectively. Finally, we annotated TE copies by searching the library using RepeatMasker (version 4.1.2) (Chen, 2004) with the parameters -e rmblast -cutoff 250 -xsmall -s -gff, and then we used RepeatCraft (Wong and Simakov, 2019) for post-processing the RepeatMasker annotations to generate less fragmented copies. Additionally, the curated consensus sequences and the TE annotations for each assembly are included as Supplementary Files . These files are compressed into a single archive named Supplementary File 3.zip , which contains the following: Curated_Consensus_Sequences.fasta, Hale_TE_Annotations.gff, Rc039_TE_Annotations.gff, and WT05_TE_Annotations.gff.

The structure of a full-length LTR-RT with LTR comprised the following, a pair of dinucleotide palindromic motifs flanking each LTR, the internal region including protein-coding sequences for gag, pol, and env, and a 5-bp target site duplication (TSD) flanking the element. LTRharvest (Ellinghaus et al., 2008) and LTR_FINDER (Xu and Wang, 2007) were used to identify flLTR-RTs in the castor bean genome. LTR_retriever (Ou and Jiang, 2018) was used to integrate the results of the two methods, and the flLTR-RT’s insertion time was calculated. The nucleotide mutation rate was set as μ = 6.9 * 10⁻⁹ (Ou and Jiang, 2018; Xu et al., 2021). Parameters were set as “gt ltrharvest -similar 90 -vic 10-seed 20 -seqids yes -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1” and “ltr_finder -D 15000 -d 1000 -L 7000 -l 100 -p 20 -C -M 0.9”.

We further analyzed the flLTR-RT sequences identified in the three genome assemblies. We performed clustering of the flLTR-RT sequences from the Hale, Rc039, and WT05 assemblies using CD-HIT (Fu et al., 2012). Each cluster’s representative sequence was selected, resulting in a total of 1,381 representative sequences. These representative sequences were then compared against the consensus sequences identified by RepeatModeler using BLAST+ (version 2.2.31) (Camacho et al., 2009). The BLAST+ parameters were set to a similarity threshold of 80% and a coverage threshold of 80%. This analysis revealed that 896 out of the 1,381 representative flLTR-RT sequences overlapped with the consensus sequences. The files castor_ltr.cluster, blast_results.txt, analyze_blast_results.py, and overlapping_sequences.txt were compressed into a single archive as Supplementary File 4.zip .

The identified flLTR-RT was clustered using vmatch dbcluster (Kurtz, 2003) with the parameter setting “-dbcluster 80 80 -identity 80 -exdrop 5 -seedlength 15 -d”. The genomic specificity of a cluster was defined by a decision tree algorithm: (1) if ≥90% members in a cluster are from a genome, then this cluster is assigned to this genome; (2) if less than 10% members in a cluster are from a genome, then this cluster is assigned to the shared cluster by another genome; (3) the remaining clusters are assigned to a common cluster of all genomes. Full-length LTR-RT sequences were queried against the CDD database of NCBI (Marchler-Bauer et al., 2015; Lu et al., 2020) to identify the reverse transcriptase domain. Full-length LTR-RTs of the Copia superfamily sequences from each genome were clustered using CD-HIT (Fu et al., 2012) and then were aligned using MAFFT (Katoh and Standley, 2013). Phylogenetic trees were constructed using IQ-TREE2 (Minh et al., 2020) and visualized using iTOL (Letunic and Bork, 2021).

Kimura distances between genome copies and TE consensus from the library were determined using buildSummary.pl, calcDivergenceFromAlign.pl, and createRepeatLandscape.pl (in RepeatMasker utility directory) on alignment files (.align files) after genome masking.

TEs significantly shape the genomic landscape of the three castor bean lineages studied. To comprehend the distribution and composition of transposable elements (TEs) within these lineages, we generated reference libraries using both de novo and homology-based approaches. Subsequently, we conducted a thorough manual inspection of these libraries. Upon detailed annotation of the individual genomes using our comprehensive TE library, we observed that the overall repeat sequence content is relatively consistent across the genomes of Rc039, WT05, and Hale, accounting for 39.17%, 34.10%, and 36.83% of the genomic sequences, respectively ( Figure 1A ). However, it is important to note that this overall consistency does not preclude variations in TE distribution within specific genomic compartments. Such variations, particularly the higher density of TEs in Rc039’s intron regions as shown in Figure 1B , complement the overall genomic TE content pattern, demonstrating the heterogeneous distribution of TEs across different genomic regions. The composition of most TE families is likewise congruent among these genomes—for example, of the 462 identified families, 447 (97%) are present in similar proportions in the three genomes, i.e., we found less than a twofold change of the proportion between genomes ( Supplementary Table S4 ). A total of 10 most abundant families account for 20% of the total transposable elements ( Supplementary Figure S1 ), consistent with other findings in wheat (Wicker et al., 2018) that a few families contribute to the vast majority of the copy number of TEs. The TEs are mainly distributed in intergenic or noncoding regions. Even at gene loci and their immediate surrounding regions, TEs distributed bias to noncoding regions, such as introns and gene flanking regions ( Figure 1B ).

In the TE landscape, a few families notably dominate in terms of copy number. Specifically, 26 families each contain over 1,500 copies, whereas a substantial number of TE families each have fewer than 500 copies. As an example, when evaluating at the superfamily level for Copia, 10 out of the 110 subfamilies hold 20,853 copies, which represents 57% of all the copies in this superfamily. A parallel pattern is found in the hAT-Ac superfamilies, where three subfamilies, making up 17% of the total, encompass 48% of all hAT-Ac copies.

LTRharvest (Ellinghaus et al., 2008) and LTR_FINDER (Xu and Wang, 2007) were used to identify flLTR-RTs in the castor bean genome. LTR_retriever (Ou and Jiang, 2018) was used to integrate the results of the two methods, and the flLTR-RT’s insertion time was calculated. Table 1 provides a comprehensive summary of the genome assemblies and flLTR-RTs identified in the castor bean genomes, including details such as genome size, genome completeness (BUSCO values), scaffold N50, and the number of flLTR-RTs. The data in Table 1 show that the Hale genome has significantly fewer flLTR-RTs compared to the Rc039 and WT05 genomes. This difference strongly correlates with genome assembly quality metrics such as N50 and genome completeness (BUSCO values). A multiple linear regression analysis revealed that N50 is a significant predictor of the number of flLTR-RTs (p = 0.036), indicating that higher N50 values are associated with a greater number of flLTR-RTs. Although the BUSCO values also showed a positive relationship, it was not statistically significant (p = 0.470). Specifically, Hale exhibits the smallest N50 at 0.50 Mb and a BUSCO completeness of 93%, suggesting a lower assembly quality which likely contributes to the reduced identification of flLTR-RTs. In contrast, Rc039, with the highest N50 of 32.06 Mb and a BUSCO completeness of 98%, shows the highest abundance of flLTR-RTs. This correlation suggests that better assembly quality facilitates a more comprehensive and accurate identification of flLTR-RTs, particularly evident in the Ty3-retrotransposons and unknown superfamilies.

We also compared the length distribution of flLTR-RTs in the three genomes ( Supplementary Figure S2 ). In the Ty3-retrotransposons superfamily, the Hale genome had fewer flLTR-RTs ranging from 10,000 bp to 15,000 bp (38) compared to Rc039 (105) and WT05 (87). A Kruskal–Wallis rank sum test was conducted to compare the length distribution of flLTR-RTs among the three genomes for both the Copia and Ty3-retrotransposons superfamilies. For the Copia superfamily, whose lengths range from 4,000 bp to 5,500 bp, the test indicated no significant difference in the length distribution among the three genomes (Kruskal–Wallis chi-squared = 0.267, df = 2, p-value = 0.875). However, for the Ty3-retrotransposons superfamily, the test indicated a significant difference in the length distribution among the three genomes (Kruskal–Wallis chi-squared = 19.54, df = 2, p-value = 5.715e-05), suggesting that the assembly quality of the Hale genome may have led to fewer flLTR-RTs in the Ty3-retrotransposons superfamily.

Additionally, a Kruskal–Wallis rank sum test was conducted to compare the length distribution of TE copies among the three genomes, indicating no significant difference (Kruskal–Wallis chi-squared = 2.1313, df = 2, p-value = 0.3445) ( Supplementary Figure S2C ). Therefore, the assembly quality of the Hale genome may not significantly impact the TE copy annotation results compared to the other two castor genomes.

In the castor bean genome, the distribution of insertion times of flLTR-RTs belonging to the Copia superfamily was similar: the expansion started six million years ago and continued until about two million years ago ( Figures 2A–C ). The distribution of insertion time of flLTR-RTs belonging to the Ty3-retrotransposons superfamily in the three genomes was significantly different. In Rc039 and WT05, the expansion activity has continued to the present. In Hale, however, the peak of the expansion activity occurred two million years ago.

To further explore the evolutionary history of full-length LTR retrotransposons in the three genomes, we performed a phylogenetic analysis ( Figure 2D ). This phylogenetic tree demonstrates that nearly every branch contains full-length LTR retrotransposons of the Copia superfamily from all three genomes, with no genome-specific clustering. This suggests that these retrotransposons have followed similar evolutionary paths across the three genomes. Additionally, a small number of them are clustered at the basal branches of the tree, while the rest of the tree consists of clades with more branches. This evolutionary pattern is consistent with the earlier analyses of insertion times. The presence of numerous closely related branches in the terminal clades further supports the recent amplification events of LTR retrotransposons. Although the assembly level of the Hale genome was inferior to that of the other two genomes, the results of the length distribution and insertion time distribution of the full-length LTR retrotransposons showed that this did not affect the Copia superfamily.

Upon examining the genomes of the Ethiopian accession (Rc039) and the Kenyan accession (WT05) of the castor bean, we observed a pronounced chromosomal collinearity for coding genes ( Figure 3A ). Both of these genomes were sequenced and assembled independently with high precision at the chromosomal level (Xu et al., 2021; Lu et al., 2022). Local synteny gene blocks were observed between the chromosomes of Rc039 and WT05, encompassing 32,152 paralogous genes, constituting 69.4% of all identified high-confidence genes, and totaling 46,325. The conservation level of homologous chromosomes between Rc039 and WT05 surpasses that of non-homologous chromosomes. Moreover, the similarities extend to gene length and intergenic distances, which remain consistent not only between Rc039 and WT05 but also in the inbred cultivar accession, Hale ( Figures 3C, D ).

Unlike gene collinearity, flLTR-RTs exhibit reduced collinearity between homologous chromosomes across varied genomes compared to non-homologous chromosomes ( Figure 3B ). Notably, flLTR-RT’s movements between chromosomes are more prevalent than those within a single chromosome ( Figure 3E —VII and Supplementary Figure S3 ).

We evaluated the dynamics of flLTR-RTs and genes, considering the collinearity between different chromosomes. While our initial analysis suggested a weak correlation between the gene collinearity matrix ( Figure 3A and Supplementary Table S5 ) and the flLTR-RT collinearity matrix ( Figure 3B and Supplementary Table S6 ), as indicated by a Spearman’s rho of 0.06 (Mantel test, p = 0.32), a further examination revealed a more nuanced relationship. Although the main patterns of gene and flLTR-RT collinearity are broadly conserved, with identical inversions observed in both as exemplified by Chr10, the proportion of collinear gene pairs across different chromosomes is 20.3%, which contrasts with a notably higher proportion of 58% for flLTR-RTs. This marked difference underscores the dynamic inter-chromosomal activity of flLTR-RTs, which is likely due to their active transposition and a higher frequency of insertion events compared to inter-chromosomal gene exchanges. Consequently, while the primary patterns of collinearity are conserved, the flLTR-RTs show an enhanced tendency for inter-chromosomal rearrangements, suggesting that their evolutionary trajectories exhibit distinct dynamics due to the influence of active transposable elements.

Furthermore, TEs exhibit an increasing density from the distal regions toward the centromere on the chromosome ( Figure 3E –I–VI). In contrast, genes demonstrate the reverse trend. Within the WT05 genome, the distribution pattern of most TEs is inversely related to that of genes, for instance, retrotransposons show a negative correlation with gene distribution on chromosomes (rho = −0.15, p < 0.01).

The origins of the Hale breed may be traced to the Ethiopian Rc039 lineage. We reassembled Hale’s scaffolds to the chromosome scale using Ethiopian Rc039 and Kenyan WT05, respectively, as reference genomes. The different k-mers with a range from 5 to 51 were used to test their potential effects on the different genomes. With the exception of the 5-mers and the 7-mers, which lacked sufficient differential k-mers for phasing, all other k-mer lengths were able to phase all three genomes perfectly. The number of the identified differential k-mers stabilized at 15-mers is close to the conclusions reached by previous related studies (Jia et al., 2022).

In Figure 4 , four genomes are represented: Rc039, WT05, and the chromosome-level Hale genome assembled using these two genomes as templates. It can be seen from the clustering tree in Figure 4 that the k-mers belonging to WT05 cluster together, while the k-mers belonging to Rc039 and the two Hale genomes cluster together. The k-mers of WT05 exhibit clear distinctions from the other three genomes, while the differential k-mers in the two Hale genomes are notably more similar to those in Rc039. Hale consistently demonstrated stronger genetic ratios with Rc039 than with WT05 ( Figure 4 ). These findings suggest that the domestication of Hale may be related to the ancestor of Rc039, although further phylogenomic analysis based on host genes is required to confirm this evolutionary relationship.

We analyzed 3,571 full-length LTR retrotransposons in castor bean genomes, grouping them based on sequence similarity. Using an 80/80 cluster criterion (80% sequence identity and 80% coverage) (Wicker et al., 2007), we identified 78 families common to Rc039, WT05, and Hale ( Figure 5A ). The phylogenetic trees, constructed using sequences from the three major shared families, reveal that these genomes are closely intertwined, pointing to a shared, undifferentiated ancestor ( Supplementary Figure S4 ). We analyzed the proportions of flLTR-RT families shared among the Rc039, WT05, and Hale genomes. We found that the correlation of the proportions of flLTR-RT families within each genome was high, with a Pearson correlation coefficient of 0.98 between each genome pair (p < 0.05), indicating statistically significant similarities ( Supplementary Table S7 ). Further pairwise comparisons showed a correlation of 0.96 between WT05 and Hale for the proportions of shared flLTR-RT families within these genomes ( Supplementary Table S8 ) and a correlation of 0.97 between Rc039 and Hale for the proportions of shared flLTR-RT families within these genomes ( Supplementary Table S9 , Pearson test, p < 0.05), suggesting strong and statistically significant similarities among these genomes’ flLTR-RT families. Our analysis also revealed that within the same subfamily of flLTR-RTs, not only is there a close evolutionary relationship but also the proportion within the genomes is similar.

To determine the evolutionary relationship among the accessions, we performed a phylogenomic analysis based on orthologous genes ( Figure 5B ). This analysis included Mercurialis annua, Manihot esculenta, and the three castor bean genomes (Hale, Rc039, and WT05). The results show that Rc039 and Hale are the closest in evolutionary terms, indicating a closer evolutionary relationship between these two genomes. This supports our hypothesis that Hale may have been domesticated from an ancestral population related to Rc039.

To infer transposition times, TEs are classified according to their K (Kimura) values. The extent of sequence alignment between a TE copy and its counterpart in the TE reference library sheds light on the timing of its transposition. Copies showing a significant similarity (K < 25) indicate a recent transposition activity, whereas copies with a reduced similarity (K > 25) suggest more distant transposition events ( Figure 6 ) (Chalopin et al., 2015).

The castor bean genome has experienced multiple TE activity bursts. The first notable burst, centered on K = 10, is marked by a significant increase in the unknown superfamily, while the other superfamilies remain comparatively unchanged. Following this, additional activity bursts are observed, including a prominent surge near K = 5, predominantly characterized by the expansion of the Ty3-retrotransposons superfamily. This consistent growth of the Ty3-retrotransposons superfamily aligns with the evolutionary patterns seen in flLTR-RTs within the castor bean genome. As shown in Figure 6 , there appear to be four distinct TE activity bursts in total.