The S. malaccense, S. aqueum, S. jambos, and S. syzygioides genome assemblies were generated from trees growing in the Masoala Hall of the Zurich Zoo in Switzerland. Voucher specimens were deposited in the Zürich herbarium (S. malaccense (ZT-00170996), S. aqueum (ZT-00170994), S. jambos (ZT-00170999), and S. syzygioides (ZT-00170991)). Samples collected from the trees were stored at -80°C until nucleic acid extraction.

High-molecular-weight genomic DNA was isolated from frozen leaves using the “ONT high-molecular-weight gDNA extraction from plant leaves” protocol (Oxford Nanopore Technologies, Oxford, UK). Following the extraction, we performed a size selection step using the Circulomics Nanobind Plant Nuclei Big DNA Kit from PacBio (Menlo Park, CA, USA). (NB-900-801-001).

Total RNA from S. malaccense, S. aqueum, S. jambos, and S. syzygioides were isolated in triplicate from whole leaves (young and mature), lamina of mature leaves, and stems. Total RNA was also isolated in triplicate from S. syzygioides’ buds (in the fruiting stage) and S. jambos’ buds (before and after flowering) and flowers.

Total RNA was extracted from frozen powder using Ambion PureLink Plant RNA Reagent (Ambion by Life Technologies, Carlsbad, CA, USA). The concentration and quality of the total RNA were assessed with an Agilent Bioanalyzer using the Agilent RNA 6000 Nano Kit (Agilent, Santa Clara, CA, USA).

DNAseq libraries were prepared from total gDNA using the Celero PCR workflow with an enzymatic fragmentation kit from Tecan (Männedorf, Switzerland). DNAseq libraries were loaded on an Illumina S2 flow cell and sequenced on the Illumina Novaseq 6000 instrument (Illumina, San Diego, CA, USA) as 2 x 151 bp paired-end reads.

Hi-C libraries were prepared from 0.2 g of frozen leaves using the Proximo Hi-C Kit following the manufacturer’s instructions (Phase Genomics, Seattle, WA, USA) and sequenced on an Illumina HiSeq 4000 instrument (Illumina) as 2 x 151 bp paired-end reads.

mRNA stranded libraries were prepared from 500 ng of total RNA using the Tecan Universal Plus mRNA-Seq library preparation kit with NuQuant^® and sequenced on an Illumina HiSeq 4000 instrument as 2 x 151 bp paired-end reads.

Illumina raw reads generated from DNAseq libraries and Hi-C libraries were cleaned using fastp 0.23.2 (--length_required 75 --low_complexity_filter) (Chen et al., 2018).10.1038/s41597-021-00968-x

Sequencing libraries were generated from high-molecular-weight gDNA and prepared for sequencing on PromethION flow cells (FLO-R0002) by using the ligation sequencing (SQK-LSK109) and flow cell priming (EXP-FLP002) kits (Oxford Nanopore Technologies, Oxford, UK). The base calling was performed by using Guppy 6.1.1 and the super accuracy plant model. Raw ONT reads were cleaned using seqkit 2.2.0 (--min-qual 9 --min-len 5000) (Shen et al., 2016) to discard reads shorter than 5,000 bp or with quality scores lower than 9.

Cleaned Illumina paired-end reads from DNAseq libraries were analyzed by GenomeScope 2.0 and smudgeplot 0.2.4 to estimate the genome size, percentage of heterozygosity, and the ploidy level using a k-mer size equal to 21 (Ranallo-Benavidez et al., 2020).

ONT cleaned reads were corrected with fmlrc2 0.1.7 (--cache_size 13 –K 21 59 79) (Mak et al., 2023) using cleaned Illumina paired-end short-reads from DNAseq libraries. The corrected ONT reads were then assembled using flye 2.9 (--read-error 0.01 --nano-hq) (Kolmogorov et al., 2019) and iteratively polished with ntedit 1.3.5 (-m 2 -i 3 -d 3 -X 0.5 -Y 0.5) using kmer profiles created with nthits 0.0.1 (--solid --outbloom -b 36) for kmers of lengths 60, 50, 40 and 30 (Warren et al., 2019) using Illumina paired-end short reads from DNAseq libraries. Haplotigs were detected and removed from the polished contigs using purge_dups 1.2.5 (Guan et al., 2020) using cutoff of 10, 315 and 645 for S. malaccense, 70, 440 and 960 for S. aqueum, and 60, 410, 960 for S. jambos, 10, 410 and 960 for S. syzygioides.

Cleaned Illumina read pairs generated from Hi-C libraries were mapped to the genomes to remove reads with low mapping scores, duplicated reads, and paired-end reads. Illumina Hi-C read pairs were mapped to the haplotig-purged contigs using minimap2 2.24 (Li, 2018) rather than bwa (Li, 2013) since we noticed that it results in assemblies of equivalent qualities in a shorter time. The scaffolding to a chromosome-scale assembly was performed using yahs 1.1a2 (-r 1000,2000,5000,10000,20000,50000,100000,200000,500000,1000000,2000000,5000000) (Zhou et al., 2022). Hi-C map files were generated with PretextMap 0.1.9 (https://github.com/wtsi-hpag/PretextMap) and used to manually curate the assemblies using PretextView 0.2.5 (https://github.com/wtsi-hpag/PretextView).

The curated genome assemblies were mapped to the S. aromaticum genome (Ouadi et al., 2022) using minimap2 2.24, visualized using a custom R script, and the orientation and names of the chromosomes were set in accordance with those of S. aromaticum. Chromosome-scale assembly completeness was assessed by using the genome evaluation mode of BUSCO 5.4.4 and the eudicots_odb10 lineage dataset (Simão et al., 2015). The QVs of the final assemblies were estimated using yak 0.1 (qv -K 2000000000) with kmer profiles created using yak 0.1 (count -k 31 -K 2000000000 -b37) (Cheng et al., 2021).

The Illumina RNAseq reads from S. malaccense, S. aqueum, S. jambos and S. syzygioides as well as those used for the clove genome annotation were cleaned, and overlapping paired-reads were merged using fastp 0.23.2 (--length_required 75 --low_complexity_filter --merge) (Chen et al., 2018) before being mapped as single cDNA reads to the assemblies using minimap2 2.24 (-ax splice:hq -G5K -N50) (Li, 2018). Gene models were then created for each RNASeq sample using scallop 0.10.5 (--min_transcript_coverage 5 --min_single_exon_coverage 50 --min_splice_bundary_hits 5 --min_mapping_quality 0) (Shao and Kingsford, 2017).This approach was used for the annotation of the clove genome, where it was observed to produce better gene models than by directly mapping paired-reads with a dedicated mapper.

To obtain models for genes that are not expressed in the RNAseq samples, the transcripts from S. aromaticum and E. grandis gene annotations were mapped to the assemblies using minimap2 2.24 (-ax splice:hq -I5G -G5K -N50 -uf) (Li, 2018), and gene models created using bedtools 2.30.0 (bamtobed -bed12) (Quinlan and Hall, 2010) and custom gawk scripts to convert the obtained bed file into a gtf file.

The final gene models were obtained by merging the RNAseq, S. aromaticum, and E. grandis gene models using taco 0.7.3 (--gtf-expr-attr TPM --filter-min-expr 10) (Niknafs et al., 2017) and adding coding sequences using Transdecoder 5.5.0 (LongOrfs -S -m 64; Predict --single_best_only --retain_blastp_hits dmd.tsv) (https://github.com/TransDecoder/TransDecoder/wiki), diamond 2.0.15 (blastp --query longest_orfs.pep --db uniref-malvids.dmnd --max-target-seqs 1 --outfmt 6 --evalue 1e-6) (Buchfink et al., 2015) and gffread 0.12.7 (Pertea and Pertea, 2020).

The eudicotyledons portion of UniProt filtered to remove proteins with poor descriptions was used to annotate the gene models with their best hit using diamond 2.0.15 (blastx --query tx.fa --db eudicotyledons.filtered.dmnd --top 10 --min-score 200 --ultra-sensitive --iterate). The illustration of the regions where genes encoding for putative eugenol synthase were predicted was generated using gggenes 0.4.0 (https://github.com/wilkox/gggenes).

Annotation of transposable elements was carried out using TE-greedy-nester 1.0.0 (--discovery_tool LTRharvest) (Lexa et al., 2020), genometools LTRharvest 1.6.2 (Ellinghaus et al., 2008) and TEsorter 1.3.0 (-db rexdb-plant --min-coverage 10 --max-evalue 0.01 --pass2-rule 70-30-80) (Zhang et al., 2022) with REXdb (Neumann et al., 2019). The insertion age of the predicted transposable elements was then calculated as previously reported (Marcon et al., 2015). In addition, Red 2.0 (Girgis, 2015), GRF 1.0 (Shi and Liang, 2019) and cd-hit 4.8.1 (grf-main -i genome.fa -c 1 -o genome.MITE --min_tr 10; cd-hit-est -i genome.MITE/candidate.fasta -o genome.MITE/clusteredCandidate.fasta -c 0.90 -n 5 -d 0 -aL 0.99 -s 0.8 -M 0; grf-mite-cluster -i genome.MITE/clusteredCandidate.fasta.clstr -g genome.fa -o genome.MITE) (Fu et al., 2012), EAHelitron (Hu et al., 2019), and tantan 39 (-f4) (Frith, 2011) were used to predict repeats, Miniature Inverted-repeat Transposable Elements (MITEs), helitron, and tandem repeats, respectively.

Synteny between the Syzygium species was done by pairwise mapping whole genomes using minimap2 2.24 (Li, 2018), identifying structural variants using syri 1.6 (Goel et al., 2019), and plotting syntenic blocks larger than 20 kb using plotsr 0.5.4 (Goel and Schneeberger, 2022).

Orthologous genes were clustered into HOGs with OrthoFinder 2.5.4 (Emms and Kelly, 2019) using the set of predicted protein sequences from the five species assemblies.

Smudgeplot and GenomeScope 2.0 were used to perform a genome profiling step using Illumina PE short-reads from DNAseq libraries as input and a K-mer length of 21 bp (Ranallo-Benavidez et al., 2020) ( Table 1 ; Supplementary Table 1 ; Supplementary Figures 1 , 2 ). The ploidy level predicted by Smudgeplot was in accordance with previous karyotype studies for the studied S. malaccense and S. jambos specimens (Oginuma et al., 1993; Pedrosa et al., 1999). S. malaccense was predicted to be a diploid specimen (2n = 2x = 22) like S. aromaticum and S. grande. The S. aqueum, S. jambos, and S. syzygioides specimens were predicted as being autotetraploid (2n = 4x = 44). The estimated monoploid genome sizes were similar among the four Syzygium species (343–372 Mb), a size range consistent with the small genome assembly sizes of S. aromaticum (370 Mb) and S. grande (405 Mb) (Low et al., 2022; Ouadi et al., 2022). The heterozygosity rate estimated by the GenomeScope 2.0 ranged from 2.3% for the diploid specimen S. malaccense to 4.3% for the autotetraploid specimen S. aqueum. These heterozygosity rates appeared to be higher than for S. aromaticum (0.18%) (Ouadi et al., 2022) and the average reported by Ellestad et al., who performed a literature review of the genome-wide heterozygosity values estimated using the software GenomeScope and GenomeScope 2.0 (Ellestad et al., 2022). They found that the average value inferred for all plant species assessed was 1.59% (1.10% for diploid plants only) noting that over half of the plant species considered were cultivated for human usage, which could affect the average value accuracy.

The four de novo chromosome-scale assemblies were constructed using long-reads from Oxford Nanopore Technologies (ONT), short paired-end reads from Illumina DNAseq libraries, and Hi-C libraries generated for each Syzygium species ( Supplementary Tables 1, 2 ).

To prevent assembly artifacts possibly caused by heterozygosity and polyploidy of the Syzygium specimens, haplotigs were detected and removed from the polished contigs. The effect of the haplotig removal step was assessed using BUSCO (Benchmarking Universal Single-Copy Orthologs) in genome mode (Simão et al., 2015). After the haplotig removal step, the number of complete and duplicated BUSCOs genes was considerably reduced in the haplotig-purged contigs (3.3% to 6.1%) when compared to the polished contigs (93.6% to 97.1%) ( Figure 1A ). Hi-C data enabled the scaffolding of contigs into 11 chromosomes. On the Hi-C contact matrices, a strong intra-chromosomal signal indicates efficient scaffolding, with the 11 chromosomes of each Syzygium assembly supported by a high number of their respective Hi-C reads ( Figure 1B ).

The final chromosome-scale assemblies for S. malaccense (430 Mb), S. aqueum (392 Mb), S. jambos (426 Mb), and S. syzygioides (431 Mb) consisted of monoploid consensus (11 chromosomes and unplaced sequences) with comparable quality metrics. A high level of quality at the base-scale (quality value [QV] between 44.006 and 45.114), of contiguity (97.5% to 99.8% of the assemblies length anchored on 11 chromosomes) and completeness (BUSCO complete genes scores of 98%) was reached for the four new assembled Syzygium genomes ( Table 2 ; Figure 2 ; Supplementary Tables 3, 4 ).

Despite their high heterozygosity rate, the quality metrics for the genome assemblies of the diploid specimen S. malaccense and the autotetraploids S. aqueum, S. jambos, and S. syzygioides were comparable to those reported for S. aromaticum assembly (370 Mb) (Ouadi et al., 2022). Nevertheless, BUSCO scores revealed a higher percentage of complete and duplicated BUSCOs in the four new assemblies compared to S. aromaticum (2.2%), principally in the genome assembly of the three autotetraploid specimens (3.3% to 5.5%) ( Figure 2 ).

The average number of protein-coding genes predicted for the four newly assembled genomes is 31,119, representing 26.52% of the genome assemblies’ size ( Table 2 ).

The annotation completeness was assessed using the BUSCO method in transcriptome and protein modes and by selecting the whole set of predicted transcripts and proteins for each gene as inputs, respectively ( Figure 2 ; Supplementary Table 3 ). BUSCO results indicated that the annotation completeness is comparable among the four newly assembled Syzygium species, with complete BUSCO scores ranging from 91.9% in S. aqueum assembly to 93.5% in S. malaccense assembly in transcript mode and from 89.3% in S. aqueum assembly to 90.9% in S. malaccense assembly in protein mode. BUSCO scores obtained for S. aromaticum by using the same assessment methods (95% in transcriptome mode and 93.7% in protein mode) were slightly superior to those of newly assembled genomes but still comparable. The loss of complete BUSCOs between the genome and protein mode assessments ranged from 7.2% in S. malaccense assembly to 8.7% in S. aqueum assembly, indicating acceptable quality of the predicted gene models and protein sets.

The genome assembly of S. aromaticum comprised multiple copies of a gene encoding for putative eugenol synthase (EGS), the enzyme that catalyzes the synthesis of eugenol from coniferyl acetate. In total, 15 copies split into 2 loci were reported: a first locus on chromosome 10 comprising 14 copies and a second locus on chromosome 11 with 1 copy (Ouadi et al., 2022). The functional annotation of the four newly assembled Syzygium species genomes revealed fewer genes encoding for putative EGS. One gene encoding for putative EGS was identified in the genome assembly of S. malaccense, two in the genome assembly of S. aqueum, and three copies were found in the genome assemblies of S. jambos and S. syzygioides. All putative EGS genes were located on chromosome 10 except for one of the three copies of S. syzygioides located on chromosome 11 ( Figure 3 ; Supplementary Table 5 ).

Effective lengths of repeat elements, which are different from their genomic length, were calculated by removing the length of the nested elements they contained. The proportions of genome assembly length occupied by predicted genes (25.97% to 27.37%) and repeat sequences (41.34% to 43.02%) appear to be conserved among the four newly sequenced Syzygium genomes ( Table 2 ). Using the same method, repeat elements in Syzygium aromaticum genome assembly represents 39.98%. The most abundant repeat elements identified in the four newly sequenced Syzygium genomes were the LTR-RTs spanning 16.97% of the assembly length for S. syzygioides to 22.35% for S. malaccense. As reported for S. aromaticum and S. grande, LTR-RTs belonging to the Gypsy superfamily were more abundant than elements belonging to the Copia superfamily in the four newly sequenced genomes ( Table 2 ; Supplementary Tables 6–9 ) (Low et al., 2022; Ouadi et al., 2022).

To identify evolutionary structural changes among the Syzygium species chromosomes, we performed a synteny analysis on the four newly assembled genomes, S. aromaticum and S. grande. The alignment of the 11 chromosomes’ DNA sequences of the 6 Syzygium species revealed a high conservation of the chromosomal organization ( Figure 4A ).

No large interchromosomal rearrangements were detected between the chromosomes of the six Syzygium species. A high percentage of the five species’ chromosome lengths were syntenic with S. aromaticum, ranging from 68.45% between S. aromaticum and S. jambos to 73.02% between S. aromaticum and S. aqueum. Intrachromosomal rearrangements such as inversions, translocations, and duplications between the chromosomes of S. aromaticum and those of the other five Syzygium species represented 5% of their 11 chromosomes length on average. In terms of number, the most frequent rearrangements observed between S. aromaticum and the five other species were duplications and translocations with average numbers of 1348 and 1325, respectively, spanning an average of 0.85% to 1.43% of the 11 chromosome lengths. Inversions were found less frequently for all species but occupied a larger fraction of the genome assemblies’ length than duplications and translocations except for S. syzygioides. The percentage of assembly lengths comprising inversions between S. aromaticum and the five other Syzygium species ranged from 0.68% between S. aromaticum and S. syzygioides to 4.83% between S. aromaticum and S. grande. Overall, the size of the inversions was relatively small. For instance, 11 inversions were detected, between chromosome 5 of S. aromaticum and S. grande, representing 17.32% of the chromosome length of S. grande (41,797,999 bp) and 1.87% of its 11 chromosomes (387,620,547 bp) ( Figure 4B ; Supplementary Table 4 ). In contrast, the synteny analysis performed between S. aromaticum and E. grandis revealed 10 intrachromosomal rearrangements on chromosomes 2, 4, 6, 8, 9, and 10 that included large terminal inversions representing up to 40% of the chromosome length of S. aromaticum. The other four chromosomes (1, 3, 5, and 7) of the two Myrtaceae species were highly syntenic (Ouadi et al., 2022). To further investigate the chromosomal architecture evolution of the Syzygium species and verify that these rearrangements were due to biological events rather than assembly artifacts, we also performed DNA alignment of the chromosome sequences of E. grandis with the those of S. malaccense, S. aqueum, S. jambos, and S. syzygioides. Chromosomes 1, 3, 5, and 7 of E. grandis and those of the four newly assembled species were also highly syntenic, and we observed the same 10 rearrangements on chromosomes 2, 4, 6, 8, 9, 10 and 11 ( Figure 4A ; Supplementary Figure 3 ).

To investigate the phylogenetic relationships among gene sequences of S. aromaticum, S. malaccense, S. aqueum, S. jambos, and S. syzygioides, the sets of predicted protein sequences from the five species assemblies were analyzed using OrthoFinder (Emms and Kelly, 2019).

A total of 49,269 hierarchical orthogroups (HOGs) were identified, including 93.7 to 95.2% of each species gene set ( Figure 5A ). Of these, 18,963 (38.5%) HOGs contained genes from all five species, and 4,928 (10%) were species specific. In more detail, 789 were specific to S. aromaticum, 950 were specific to S. malaccense, 940 HOGs were specific to S. aqueum, 1009 HOGs were specific to S. jambos, and 1240 HOGs were specific to S. syzygioides. Pairwise, S. aromaticum and S. aqueum appear to share the lowest number of orthogroups (625). The highest number of shared HOGs inferred between each pair of studied species was found between S. aqueum and S. syzygioides (1218), followed by S. aqueum and S. malaccense (1152), and S. jambos and S. malaccense (1027). The species tree resulting from the analysis of the HOGs divided the Syzygium species studied into two groups based on closer relationships: the first group comprising S. aromaticum and S. aqueum and a second group comprising S. jambos, S. malaccense, and S. syzygioides ( Figure 5B ).

To clarify the dynamic activity of full-length LTR-RTs belonging to the superfamilies Gypsy and Copia within the Syzygium genus, we identified the lineages belonging to each superfamily located on the chromosomes of S. malaccense (429 Mbp), S. aqueum (387 Mbp), S. jambos (416 Mbp), and S. syzygioides (425 Mbp) and estimated their insertion time. Then, we compared the repertoires’ compositions and repeat element insertion times of the four species with those of S. aromaticum (368 Mbp).

We found that S. malaccense and S. aromaticum, the largest and smallest chromosome-scale assemblies of this study, contained the highest (8427) and lowest number (6167) of LTR-RTs in Gypsy and Copia, respectively ( Figure 6A ; Supplementary Tables 6–9 ). In the five Syzygium species’ chromosomes, we identified a higher number of LTR-RTs for Gypsy than Copia, with a ratio of Gypsy to Copia content ranging from 1.09 for S. syzygioides to 1.45 for S. malaccense. The Gypsy superfamily comprised a higher proportion of nested elements (17.37% to 24.47%) compared to the Copia superfamily (7.01% to 9.44%), suggesting distinct accumulation and mobile activity of both superfamilies in all five species. Our results revealed little variation in the number of Copia elements on the chromosomes of S. aqueum (2705 elements) and S. aromaticum (2809), the two smallest chromosome-scale assemblies, and on the chromosomes of S. syzygioides (3290), S. jambos (3324), and S. malaccense (3433). In contrast, we found a notably higher accumulation of Gypsy elements (4994) in the chromosomes of S. malaccense compared to the four other species. The ratio of Gypsy content varied from 1.35 when comparing S. malaccense with S. jambos to 1.49 when comparing S. malaccense with S. aromaticum. It represented a difference in Gypsy effective length of 19,402,234 bp to 21,766,176bp, respectively. In the five Syzygium chromosome-scale assemblies, the most abundant lineage belonged to the Gypsy superfamily, but it varied according to the species. The Gypsy lineage Tekay was the most represented for S. aromaticum (1534 elements), S. jambos (1674 elements), and S. syzygioides (2090 elements). At the same time, for S. malaccense and S. aqueum, we found a higher abundance of the gypsy lineage Ogre (2382 and 1799 elements, respectively). Among the Gypsy superfamily, the most abundant lineages, Tekay and Ogre, were those with the highest proportion of nested elements (19.10% to 28.55% and 16.69% to 27.92%, respectively) in all five species. For S. aromaticum, S. malaccense, and S. syzygioides, the proportion of nested elements belonging to the Athila lineage was also among the highest identified ( Figure 6B ; Supplementary Figures 4, 5 ).

Regarding the Copia superfamily, the most represented lineages on the chromosomes of the five Syzygium species were Ale (608 to 873 elements), followed by the lineage Tork (456 to 762 elements) for S. malaccense, S. aqueum, S. jambos and S. Syzygoides, and the lineage SIRE (502 elements) for S. aromaticum.

The insertion times of 97.13% of the full-length Gypsy and Copia elements identified in the five Syzygium species were estimated (33,861elements). Nearly all elements (97.33%) were inserted in the last 5 million years (32,958 elements) ( Figure 7 ). During this time period, distinct insertion activities of the two superfamilies occurred in the five Syzygium species.

Compared to the other four Syzygium species, the chromosomes of S. aromaticum underwent a more ancient wave of Gypsy insertions (peak at ~2.5 million years ago [Mya]), principally attributed to the Tekay elements, the most abundant lineage in this species ( Figure 7A ). We also found that a few recent insertions (18.02% of insertions) occurred in S. aromaticum chromosomes within the last one million years. In contrast, a recent burst of Gypsy insertions (~0–1 Mya) occurred in four other species chromosomes: most insertions of Gypsy in S. malaccense (44.53%), S. aqueum (44.43%), S. jambos (52.55%), and S. syzygioides (36.45%) were less than one million years old. We inferred that the high number of Gypsy LTR-RTs found in S. malaccense may be attributable to two successive waves of insertions: a peak of Tekay insertions at ~2 Mya and a more recent peak of Ogre at ~1 Mya.

Similar to what we observed for the Gypsy superfamily, the insertion of Copia elements occurred earlier in S. aromaticum compared to the four other species, with fewer recent insertions ( Figure 7B ). Compared to the Gypsy elements, a smaller proportion of recent Copia insertions (less than one million years old) were detected in S. aromaticum (10.66%), S. malaccense (24.16%), S. aqueum (26.49%), and S. jambos (36.90%) suggesting a distinct recent insertion pattern of the two superfamilies in the four species. However, we found a comparable proportion of Gypsy (36.45%) and Copia (32.22%) elements that were less than one million years old in S. syzygioides, the species for which we found the lowest ratio of Gypsy to Copia content (1.09).