Fresh leaves for genome sequencing were collected from a female E. ulmoides plant ‘Huazhong No. 8’, which was used to assemble the Female V1 genome. To observe the morphological structure of female and male flower bud development, flower buds were collected from female ‘Huazhong No. 6’ and the male ‘Huazhong No. 11’ during the entire development period. Samples were obtained every 5 days from the beginning of bud germination in late April to the end of September, once every 10 days from early October to January of the following year, and once every 3 days from February to early April of the following year, with 20 female and male flower buds being collected each time. The experimental materials for sex differentiation analysis were derived from the collected flower buds. Male and female flower buds at the floral organ induction stage (the inflorescence primordium formation stage), floral organ morphological differentiation initial stage (the pistil and stamen differentiation stage), and flower organ maturity stage were selected for transcriptome sequencing, with three biological repetitions for each sample. To analyze the expression of α-linolenic acid synthesis genes, four different tissues (stem, bark, leaf and fruit) of female ‘Huazhong No. 8’ were extracted. All plant samples were collected from the Yuanyang Experimental Base of Research Insitute of Non-timber Forestry, Chinese Academy of Forestry. Plant materials were immediately immersed in liquid nitrogen, transported to the laboratory, and then stored in a -80°C freezer for sequencing.

After collection, the materials were immediately immersed in FAA fixative solution (70% ethanol: glacial acetic acid: 38% formaldehyde =18:1:1), transported to the laboratory, and then stored in a -80° freezer. Materials preserved in FAA fixative solution were sliced using a conventional paraffin-sectioning method. Then the slices were observed and photographed using OLYMPUS optical microscope.

Total genomic DNA was extracted from fresh leaves using the QIAGEN^® Genomic DNA extraction kit (Cat#13323, QIAGEN) for PacBio sequencing on a PacBio Sequel II instrument. For Hi-C library construction and sequencing, restriction endonucleases (HindIII/MboI) were used for chromatin digestion, and after biotin labeling, flat-end binding and DNA purification. Hi-C samples were prepared and sampled for DNA quality testing. The Hi-C fragment consisted of removed biotin that was fragmented by ultrasound, end-repaired, and after the addition of base A, the biotin-containing fragment was isolated and then sequenced to form the jointed product. Then, Polymerase Chain Reaction (PCR) conditions were tested, and DNA was amplified to obtain the library products (Kong and Zhang, 2019). An Illumina HiSeq X Ten sequencer was used for sequencing after the library quality control was confirmed.

RNA was extracted using a TRIzol kit (Invitrogen) according to the manufacturer’s instructions. The RNA concentration was measured using a NanoDrop spectrophotometer, and RNA purity and integrity were determined by an Agilent 2100 Bioanalyzer and 1% agarose gel electrophoresis, respectively. After passing the quality inspection, a cDNA library was constructed according to the instructions accompanying the cDNA Library Construction Kit (NEB). Paired-end sequencing was performed on the cDNA library using the Illumina HiSeq X-10 platform, and 150 bp paired-end sequences were generated.

For the genome assembly of female E. ulmoides plant ‘Huazhong No. 8’, first, sub-sequences from PacBio sequencing were assembled using Falcon v2.0.5 (https://github.com/PacificBiosciences/FALCON/) with the parameter ‘–max diff 100 –max cov 100 –min cov 2 –min len 5000’, the pure third-generation assembly software officially launched by PacBio. Then, Arrow (Chin et al., 2013) was applied to align the third generation data to the preliminary assembly for correction with default parameters. After the initial genome assembly was completed, Pilon v1.22 (Walker et al., 2014) with the parameter ‘–mindepth 10 –changes –fix bases’ was used for iterative polishing. The Hi-C data were aligned to the genome assembly via Juicer v1.5 software with the following parameter settings: -s DpnII -t 20 (Durand et al., 2016b). Finally, the draft genome of E. ulmoides was assembled onto chromosomes by 3D-DNA pipeline with the parameter ‘-m haploid -s 4 -c 17 -j 15’ (Dudchenko et al., 2017). The contact matrix and heatmap of chromosomes were drawn with Juicebox (Durand et al., 2016a).

As for the male E. ulmoides tree (SNJ) genome reassembly, named Male V2 genome, we first downloaded the previous genome assembly data from the NCBI under accession number SRP095726. Next, the downloaded genome was refined using Hi-C data with the same method as that used for the female genome assembly. Benchmarking Universal Single-Copy Orthologs (BUSCO) v4.1.4 (https://busco.ezlab.org/) was used to assess the integrity and accuracy of the final genome assembly with the Eudicotyledons_odb10 database (Simão et al., 2015).

TRF v4.07 (Benson, 1999) (tandem repeats finder) software (https://tandem.bu.edu/trf/trf.html) was applied to predict tandem repeats of the genome. Transposable elements (TEs) were identified by combining de novo-based and homology-based approaches. In the de novo-based approach, RepeatModeler v1.0.11 (Flynn et al., 2020) (http://www.repeatmasker.org/RepeatModeler/) and LTR-FINDER v1.05 (Xu and Wang, 2007) were used to build the repeat library (LTR length 100 to 5000 nt; length between two LTRs: 1000 to 20,000 nt), and then, RepeatMasker v4.0.7 (Tarailo-Graovac and Chen, 2009) (http://www.repeatmasker.org/) was used to identify and classify repeats under the parameter ‘-a -e ncbi -q -norna -nolow -div 30 -cutoff 225’. For the homology-based approach, the known repetitive sequences database Repbase (Bao et al., 2015) (https://www.girinst.org/repbase/) was searched by RepeatMasker to find homologous sequences in the E. ulmoides genome.

Gene prediction for the E. ulmoides genome was performed by integrating three different methods: de novo prediction, homology-based prediction and transcriptome-based prediction. De novo gene prediction was achieved using Augustus v3.3 (parameter: -strand=both -genemodel=partial -gff3=on -species= arabidopsis) (Stanke and Morgenstern, 2005) (http://bioinf.uni-greifswald.de/augustus/), SNAP (Korf, 2004) and GlimmerHMM v3.52 (Majoros et al., 2004). For the homology approach, BLASTP (Camacho et al., 2009) was used to map protein sequences onto the E. ulmoides genome, and GeneWise v2.4.1 (Birney et al., 2004) was used to align the homologous genome sequences with the matching proteins. RNA-seq sequences were mapped to the genome assembly using PASA v2.0.2 to identify putative exon regions and splice junctions. Finally, EVidenceModeler v1.1.1 (Haas et al., 2008) (https://evidencemodeler.github.io/) was employed to integrate the gene sets predicted by the three methods. For gene annotation, EggNOG-mapper v4.5 (Huerta-Cepas et al., 2017) software was employed with emapper DB 4.5.1 to obtain the Clusters of Orthologous Genes (COG), eggNOG, Gene Ontology (GO) and KEGG pathway information for each gene. HMMER v3.2.1 (Prakash et al., 2017) was applied to search the Pfam database (http://pfam.xfam.org/) and obtain the protein domains. Functional annotation of protein coding genes was performed using the UniProt database (https://www.uniprot.org/).

Four types of non-coding RNAs (ncRNAs), tRNA, rRNA, miRNA, and snRNA, were annotated in the E. ulmoides genome. TRNAscan-SE v1.4 (Lowe and Chan, 2016) (http://lowelab.ucsc.edu/tRNAscan-SE/)software was employed to identify the tRNAs with tRNAscn-SE -i -q and eukaryote parameters. rRNAs were detected using BLASTN alignment of rRNA sequences known related species. Other ncRNAs were predicted using INFERNAL v1.1 (Nawrocki et al., 2009) software with the parameter ‘cmsearch –ga –incE 0.01 -E 10.0’ to search the Rfam v12.0 (Nawrocki et al., 2015) database (http://rfam.xfam.org/).

To identify whole genome duplication (WGD) events in E. ulmoides, BLASTP was used to perform a homology search using the female and male genome, and MCScanX v0.8 (Wang et al., 2012) (https://sourceforge.net/projects/mcscanx/) was employed to detect syntenic blocks. The collinearity analysis of male and female was conducted by JCVI v1.2.1. Ks (synonymous substitution) values for paralogous blocks were calculated by KaKs_Calculator v2.0 with the parameter ‘-i test.axt -m MA -c -o’ (Wang et al., 2010), and a frequency distribution graph of the Ks values was drawn to identify WGD events in E. ulmoides. The WGD event time was estimated by the formula ks/2r, with a r value of 8.25 e-9.

To identify MADS-box genes, HMMER v3.2.1 software with the SRF-TF domain (PF00319) and K-box (PF01486) were used to search against the E. ulmoides genome proteins (–cut_ga, E value<1×10⁻⁵), which were obtained from Pfam. The search results containing MADS-box domain or K-box domain were retained, and incomplete functional domain sequences were removed. Using all Arabidopsis MADS-box genes as queries, the predicted E. ulmoides MADS genes were checked by BLASTP searches (E value<1×10⁻⁵). The phylogenetic tree was then constructed using RAxML v8.2.12 with the GTRGAMMA substitution model and 1000 bootstraps on the CIPRES website (https://www.phylo.org/portal2/home.action). The phylogenetic tree was visualized using FigTree v1.4.4. The combination of genome annotation, homology search and MADS-box gene expression in male and female flower buds, enabled the ABCDE model genes to be identified. Using all female E. ulmoides MADS-box genes as queries, BLASTP was used to map protein sequences onto the male genome. For these protein sequences, the phylogenetic tree was constructed through MEGA-X using the maximum likelihood method and the JTT matrix-based model with 1,000 bootstraps. Multiple sequence alignments were performed by CLUSTALW program and viewed in GeneDoc. The gene structure of these MADS-box genes was demonstrated by TBtools. The data of Nymphaea colorata was downloaded from BIG Data Center with the accession number GWHAAYW00000000. And the phylogenetic trees of genes related to sex differentiation in E. ulmoides and N. colorata were constructed using MEGA7 with the Neighbor Joining method.

The KEGG database (https://www.kegg.jp/kegg/kegg1.html) was queried to identify genes involved in α-linolenic acid synthesis and metabolism. First, known genes related to α-linolenic acid biosynthesis and metabolic pathways in KEGG were selected and their protein sequences were downloaded. The downloaded protein sequences were aligned with the E. ulmoides genome using the NCBI online comparison function BLASTP (https://blast.ncbi.nlm.nih.gov/Blast.cgi), and the annotation information was combined to predict genes involved in α-linolenic acid synthesis and metabolism.

FastQC v0.11.9 was used to estimate the quality of RNA sequences. To obtain clean sequences, Trimmomatic v0.36 was used to remove adapter sequences and low-quality sequences (with a base quality value less than 30 or greater than 5% unknown bases). To estimate the gene expression level, clean sequences from RNA-seq were mapped onto the final Female V1 genome assembly using Tophat2 v2.1.1 with the parameters: –b2-sensitive -N 2 -p 6 -g 10 –read-edit-dist 2 –read-mismatches 3 –read-gap-length 2 (Kim et al., 2013), and then fragments per kilobase per million (FPKM) of each gene was performed by Cufflinks v2.2.1 with the following parameter: -p 6 -g file.gff -u –library-type fr-unstranded (Trapnell et al., 2010) (http://cole-trapnell-lab.github.io/cufflinks/). Finally, the values of log₂(FPKM+1)-that represented the gene expression level and the differentially expressed genes (DEGs) were identified using Cuffdiff v2.2.1, with a p-value (FDR) < 0.05 and |log2 (fold change)| > 1.

With the growth and development of female and male flower buds of E. ulmoides, the cells of the flower bud meristem continued dividing, and the internal and external morphological structure gradually changed. The floral organ development of the female and male flower buds of E. ulmoides can be divided into four stages: inflorescence primordium formation, bract differentiation, pistil and stamen differentiation, and pistil and stamen morphological formation (Figure 2).

In the inflorescence primordium formation stage, the female and male flower buds exhibit similar morphological characteristics. The buds are young and green, approximately triangular, with white fluff on the surface (Figures 2A1, C1). The internal anatomical structure shows that the flower buds are apical and conical, with thick cytoplasm (Figures 2B1, D1).

In the bract differentiation stage, the female and male flower bud morphological characteristics remain relatively similar, and males and females cannot be accurately distinguished. The outer bracts of female and male flower buds gradually lignify into scales, and the color of the flower buds gradually changes from green to light brown (Figures 2A2, C2). The number of bracts increases, and the meristem at the bud apex widens and flattens (Figures 2B2, D2).

In the pistil and stamen differentiation stage, the internal structure of female and male flower buds showed obvious morphological differences. Round small protuberances form between the axils of the primordia of the female buds, namely the pistil primordia, and several stamen primordia that form between the axils of the male buds are clustered between the bracts. The male and female buds can be distinguished according to the sectioning results (Figures 2B3, D3).

In the pistil and stamen morphological formation stage, female flower buds grow faster longitudinally and are approximately conical in shape, while the stamens are approximately spherical (Figures 2A4, C4). A single pistil primordium is differentiated from the base to the center on both sides of the female flower bud, while the male flower bud is differentiated from several stamens to form a stamen cluster inside the bract (Figures 2B4, D4). The pistil shape is pinnate, with a central depression at the top that forms two protuberances, namely the pistil stigma. After the pistil completes its morphogenesis, and stamen pollen sac differentiates into the obvious endothecium, middle layer, and tapetum (Figures 2B5, D5). Then, the internal and external morphology of the flower bud does not change, and the bud enters the dormancy stage.

A combined strategy of PacBio and Hi-C sequencing technologies was used to construct the high-quality female E. ulmoides genome assembly (Female V1). The PacBio long sequences were corrected and assembled into the preliminary genome assembly using FALCON. The obtained genome size was 957 Mb, which consisted of 2,662 contigs with contig N50 that was 1.3 Mb in length. Finally, with the assistance of 184.7 Gb of Hi-C sequences, the assembled contigs were anchored to 17 pseudochromosomes of 1.01 Gb with superscaffold N50 that was 51.89 Mb in length, and the GC content of the genome was 35.14% (Table 1). The interactions of contigs on pseudochromosomes were detected to adjust their order and direction in the Hi-C heat map, which also indicated the high quality of Female V1 genome assembly (Supplementary Figures 1, 2). The Male V2 genome was re-constructed using Hi-C technology. The genome size of the new version (Male V2) was 1.24 Gb, with 23 superscaffolds and the superscaffold N50 that was 48.30 Mb long, and the GC content was 35.19% (Table 1). Genomic alignments of chromosomes and superscaffolds between Female V1 and Male V2 are showed in (Figure 3C).

To assess the quality of the two newly assembled genomes, BUSCO was carried out. Among the 2,121 plant-specific orthologs, 2,024 (95.4%) and 2,026 (95.5%) BUSCOs were identified in the Female V1 and Male V2 assembly respectively. In addition, 1,975 (93.2%, Female V1) and 1,954 (92.1%, Male V2) of the identified orthologs in the two genomes were considered to be complete, which were higher than the 90.0% completed genes identified in the Male V1 genome (Supplementary Figure 3 and Supplementary Table 1). The BUSCO results suggested that both genomes were assembled with high quality.

The Female V1 and Male V2 genome consisted of 68.26% (688.75 Mb) and 62.25% (770.88 Mb) repeats, respectively (Table 1 and Figure 4). In both genomes, LTR (long terminal repeat) was the most abundant type of repetitive element, accounting for 40.15% (405.08 Mb) and 36.6% (453.30 Mb) of the genome size, respectively. Additionally, 59.05 Mb (5.85% of the Female V1 genome) and 64.56 Mb (5.21% of the Male V2 genome) DNA transposons were identified, which were the second highest repeated sequence in both genomes (Supplementary Table 2). The LTR-RT insertion time was estimated to provide a retrospective view of the expanded history of LTR-RTs in the two newly assembly E. ulmoides genomes. LTR-RTs gradually accumulated over 5 million years ago, and there was an insertion explosion approximately 2 million years ago in both genomes, with the Male V2 genome undergoing another peak at 0.15 Mya (Figure 3A).

A combination of ab initio prediction, RNA-seq, and homology-based search was used to predict the protein-coding genes in the E. ulmoides genome. Overall, 31,665 and 37,998 genes were predicted in the Female V1 and Male V2 genome, and the average protein-coding gene size was 6,273 bp and 6,199 bp respectively (Table 1). Among these predicted protein-coding genes, 24,049 (75.95%) genes from the Female V1 genome and 27,284 (71.8%) genes from the Male V2 genome were functionally annotated in the GO, KEGG, EggNOG, Pfam and UniProt databases (Supplementary Table 3). In addition, 2,488 and 2,910 noncoding RNAs were separately annotated in the two newly assembled genomes. The Female V1 genome contained 976 tRNAs, 178 rRNAs, 1,141 snRNAs and 193 miRNAs. Similarly, four types of noncoding RNAs were identified in the Male V2 genome, including 1,216 tRNAs, 214 rRNAs, 1,253 snRNAs and 182 miRNAs (Supplementary Figures 4, 5 and Supplementary Table 4).

WGD analysis was performed by calculating synonymous substitutions per synonymous site (Ks) values in the two newly assembled genomes. The density distribution of Ks values exhibited two peaks both in Female V1 and Male V2, with Ks values of 0.45 and 2.29 (Figure 3B). After the early γ duplication event that affected eudicots, E. ulmoides experienced a more recent WGD event approximately 27.3 Mya. The collinearity analysis of male and female revealed a 1:1 synteny pattern (Supplementary Figure 6). The divergence time between female and male E. ulmoides was defined with a Ks peak value of 0.075 (Figure 3B).

In the Female V1 and Male V2 genome, 76 and 81 MADS-box genes were separately identified, and further phylogenetic analysis of these genes was performed (Supplementary Figure 7). Based on the transcriptome analysis, 30 MADS-box genes of Female V1 were differentially expressed at different developmental stages of male and female flower buds. In addition, nine genes involved in the ABCDE model of floral development were identified in these MADS-box genes, which included three A-class genes (EuAP1, EuFUL1, EuAGL6), one B-class gene (EuAP3), one C-class gene (EuAG), two D-class genes (EuSHP1, EuSHP2), and two genes (EuSEP1, EuSEP2) that were associated with the E function (Supplementary Table 5). Nine MADS-box genes involved in the ABCDE model of floral development were identified in Male V2. We found high homology between male and female lineages (Supplementary Figure 8) and the gene structure and domain sequence of these genes were similar (Supplementary Figures 9, 10). In addition, we compared the obtained sex differentiation genes with the related sex differentiation genes in Nymphaea colorata, and the genes involved in sex differentiation in E. ulmoides and N. colorata showed high homology (Supplementary Figure 11).

The expression levels of the ABCDE genes at the floral organ maturity stage were higher than those in the floral organ induction stage and morphological differentiation initial stage. Except for the EuFUL1 and EuSHP2 genes, the other seven ABCDE genes were differentially expressed between the male and female flower buds. Seven ABCDE genes were differentially expressed during the floral organ maturity stage, with four differentially expressed genes in the floral organ morphological differentiation initial stage, and only one gene differentially expressed in the floral organ induction stage. The EuAP3 gene was differentially expressed in three different developmental stages of male and female flower buds (Figure 5 and Supplementary Figure 12).

According to the functional annotation, 12 genes encoding key enzymes involved in the main synthesis and metabolic pathways of linolenic acid were identified in the Female V1 genome. Among them, six genes of four key enzymes related to the synthesis of α-linolenic acid were identified, including two acyl-ACP desaturase (FAB2) genes—EU0119133 and EU0120166, which are homologues of SAD and perform the same function of converting stearoyl-ACP (C18:0-ACP) to oleoyl-ACP (C18:1-ACP) (Kazaz et al., 2020); one fatty acyl-ACP thioesterase A (FATA) gene—EU0103200; two omega-6 fatty acid desaturase (FAD2, FAD6) genes—EU0105412 and EU0128492; and one omega-3 fatty acid desaturase (FAD7) gene—EU0103017. Six key enzyme genes for α-linolenic acid metabolism were also identified, including four lipoxygenase (LOX) genes—EU0114412, EU0114414, EU0119724, and EU0119725—and two alpha-dioxygenase (DOX) genes, EU0107025, EU0131326 (Supplementary Table 6).

In addition, the expression levels of 12 identified key enzyme genes were compared between four parts of the fruit, stem, bark, and leaf. A heatmap of gene expression levels was also drawn (Figure 6). The expression levels of α-linolenic acid synthesis-related genes were higher in fruits and leaves than that in bark and stems. Overall, the expression levels of α-linolenic acid synthesis-related genes were much higher than those of metabolism-related genes. In particularly, the gene EU0103017—encoding FAD7, a key enzyme for linolenic acid synthesis—was highly expressed.