Genome and transcriptome sequencing

Total DNA for genome sequencing was extracted from young leaves. Leaf RNA was extracted from 18 water lily species: N. colorata, Euryale ferox, Brasenia schreberi, Victoria cruziana, Nymphaea mexicana, Nymphaea prolifera, Nymphaea tetragona, Nymphaea potamophila, Nymphaea caerulea, Nymphaea rubra, Nymphaea ‘midnight’, Nymphaea ‘Choolarp’, Nymphaea ‘Paramee’, Nymphaea ‘Woods blue goddess’, Nymphaea gigantea ‘Albert de Lestang’, N. gigantea ‘Hybrid I’, Nymphaea ‘Thong Garnjana’ and Nuphar lutea. In addition, for transcriptome sequencing we sampled several organs and tissues from N. colorata including mature leaf, mature leafstalk, juvenile flower, juvenile leaf, juvenile leafstalk, carpel, stamen, sepal, petal and root.

For PacBio sequencing, we prepared approximately 20-kb SMRTbell libraries. A total of 34 SMRT cells and 49.8 Gb data composed of 5.5 million reads were sequenced on PacBio RSII system with P6-C4 chemistry. All transcriptome libraries were sequenced using the Illumina platform, generating paired-end reads. For the Hi-C sequencing and scaffolding, a Hi-C library was created from tender leaves of N. colorata. In brief, the leaves were fixed with formaldehyde and lysed, and the cross-linked DNA was then digested with MboI overnight. Sticky ends were biotinylated and proximity-ligated to form chimeric junctions, which were physically sheared to and enriched for sizes of 500–700 bp. Chimeric fragments representing the original cross-linked long-distance physical interactions were then processed into paired-end sequencing libraries and 346 million 150-bp paired-end reads, which were sequenced on the Illumina platform.

Sequence assembly and gene annotation

To assemble the 49.8 Gb data composed of 5.5 million reads, we filtered the reads to remove organellar DNA, reads of poor quality or short length, and chimaeras. The contig-level assembly was performed on full PacBio long reads using the Canu package22. Canu v.1.3 was used for self-correction and assembly. We then polished the draft assembly using Arrow (https://github.com/PacificBiosciences/GenomicConsensus). To increase the accuracy of the assembly, Illumina short reads were recruited for further polishing with the Pilon program (https://github.com/broadinstitute/pilon). The genome assembly quality was measured using BUSCO (Benchmarking Universal Single-Copy Orthologues)23 v.3.0. The paired-end reads from Hi-C were uniquely mapped onto the draft assembly contigs, which were grouped into chromosomes and scaffolded using the software Lachesis (https://github.com/shendurelab/LACHESIS).

Genscan (http://genes.mit.edu/GENSCAN.html) and Augustus24 were used to carry out de novo predictions with gene model parameters trained from Arabidopsis thaliana. Furthermore, gene models were de novo predicted using MAKER25. We then evaluated the genes by comparing MAKER results with the corresponding transcript evidence to select gene models that were the most consistent on the basis of an AED metric.

The evolutionary position of water lily and divergence-time estimation

LCN genes were identified based on OrthoFinder26 results. The orthologues were obtained from six monocots (Spirodela polyrhiza, Zostera marina, Musa acuminata, Ananas comosus, Sorghum bicolor and Oryza sativa) and six eudicots (Nelumbo nucifera, Vitis vinifera, Populus trichocarpa, A. thaliana, Solanum lycopersicum and Beta vulgaris), N. colorata, Amborella, and the gymnosperms G. biloba, P. abies and P. taeda. LCN genes needed to meet the following requirements: strictly single-copy in N. colorata, Amborella, G. biloba, P. abies or P. taeda, and single-copy in at least five of the 12 eudicots or monocots. With G. biloba, P. abies or P. taeda as the outgroup, we identified 2,169, 1,535 and 1,515 orthologous LCN genes, respectively. Furthermore, we trimmed the sites with less than 90% coverage. LCN gene trees were estimated from the remaining sites using RAxML v.7.7.8 using the GTR+G+I model for nucleotide sequences (Fig. 1c) and the JTT+G+I model for amino acid sequences (Supplementary Note 4.1). To account for incomplete lineage sorting and different substitution rates, we applied the multispecies coalescent model and a supermatrix method, respectively, to the LCN genes and found further support for the sister relationship between Amborella and all other extant flowering plants (Supplementary Note 4.2).

We further carefully selected five LCN gene sets (1,167, 834, 683, 602 and 445) from 115 species and applied both a supermatrix method27,28,29 and the multi-species coalescent model to infer the phylogeny of angiosperms (Supplementary Note 4.2). The phylogeny inferred from 1,167 LCN genes is shown in Fig. 1d, with different support values from the multi-species coalescent analyses of the other four LCN gene sets.

To estimate the evolutionary timescale of angiosperms, we calibrated a relaxed molecular clock using 21 fossil-based age constraints7 throughout the tree, including the earliest fossil tricoplate pollen (approximately 125 Ma) associated with eudicots30. We concatenated 101 selected genes (205,185 sites) and fixed the tree topology to that inferred from our coalescent-based analysis of 1,167 genes from 115 taxa. We performed a Bayesian phylogenomic dating analysis of the 101 selected genes in MCMCtree, part of the PAML package31,32, and used approximate likelihood calculation for the branch lengths33. Molecular dating was performed using an auto-correlated model of among-lineage rate variation, the GTR substitution model, and a uniform prior on the relative node times. Posterior distributions of node ages were estimated using Markov chain Monte Carlo sampling, with samples drawn every 250 steps over 10 million steps following a burn-in of 500,000 steps. We checked for convergence by running the analysis in duplicate and checked for sufficient sampling.

We also implemented the penalized likelihood method under a variable substitution rate using TreePL34 and r8s35, as a constant substitution rate across the phylogenetic tree was rejected (P < 0.01) for all cases by likelihood-ratio tests in PAUP36. Three fossil calibrations, corresponding to the crown groups of Lamiales, Cornales and Laurales, were implemented as minimum age constraints in our penalized likelihood dating analysis, except that the earliest appearance of tricolpate pollen grains (about 125 Ma)30 was used to fix the age of crown eudicots. We determined the best smoothing parameter value of the concatenated 101 LCN genes as 0.32 by performing cross-validations of a range of smooth parameters from 0.01 to 10,000 (algorithm = TN; crossv = yes; cvstart = −2; cvinc = 0.5; cvnum = 15). We used 100 bootstrap trees with branch lengths generated by RAxML37 to infer the 95% confidence intervals of age estimates (Supplementary Note 4.2).

Identification of WGD

The N. colorata genome was compared with each of the other genomes by pairwise alignment using Large-Scale Genome Alignment Tool (LAST; http://last.cbrc.jp/). We defined syntenic blocks using LAST hits with a distance cut-off of 20 genes apart from the two retained homologous pairs, in which at least four consecutive retained homologous pairs were required. We then obtained the one-to-one blocks to exclude ancient duplication blocks with QUOTA-ALIGN38.

K S -based paralogue age distributions were constructed as previously described39. In brief, the paranome was constructed by performing an all-against-all protein sequence similarity search using BLASTP with an E-value cut-off of 10−10, after which gene families were built with the mclblastline pipeline (v.10-201) (micans.org/mcl). Each gene family was aligned using MUSCLE (v.3.8.31)40, and K S estimates for all pairwise comparisons within a gene family were obtained using maximum likelihood in the CODEML program41 of the PAML package (v.4.4c)31. We then subdivided gene families into subfamilies for which K S estimates between members did not exceed a value of 5.

To correct for the redundancy of K S values (a gene family of n members produces n(n − 1)/2 pairwise K S estimates for n − 1 retained duplication events), we inferred a phylogenetic tree for each subfamily using PhyML42 with the default settings. For each duplication node in the resulting phylogenetic tree, all m K S estimates between the two child clades were added to the K S distribution with a weight of 1/m (in which m is the number of K S estimates for a duplication event), so that the weights of all K S estimates for a single duplication event summed to one. Paralogous gene pairs found in duplicated collinear segments (anchor pairs) from N. colorata were detected using i-ADHoRe (v.3.0) with ‘level_2_only = TRUE’43,44. The identified anchor pairs are assumed to correspond to the most recent WGD event.

The K S -based orthologue age distributions were constructed by identifying one-to-one orthologues between species using InParanoid45 with default settings, followed by K S estimation using the CODEML program as above. K S distributions for one-to-one orthologues between N. colorata and each of V. cruziana, N. advena, C. caroliniana, I. henryi and Amborella were used to compare the relative timing of the WGD in N. colorata with speciation events within Nymphaeales. K S distributions for one-to-one orthologues between the outgroup species I. henryi and each of N. lutea, N. advena, N. mexicana, Nymphaea ‘Woods blue goddess’, N. colorata, and C. caroliniana were used to estimate and compare relative substitution rates among these Nymphaealean species. Additional comparisons using V. vinifera and Amborella as outgroup species instead of I. henryi gave similar results (data not shown).

Absolute dating of the identified WGD event in N. colorata was performed as previously described46. Briefly, paralogous gene pairs located in duplicated segments (anchor pairs) and duplicated pairs lying under the WGD peak (peak-based duplicates) were collected for phylogenetic dating. We selected anchor pairs and peak-based duplicates present under the N. colorata WGD peak and with K S values between 0.7 and 1.2 (grey-shaded area in Extended Data Fig. 2b) for absolute dating. For each WGD paralogous pair, an orthogroup was created that included the two paralogues plus several orthologues from other plant species as identified by InParanoid45 using a broad taxonomic sampling: one representative orthologue from the order Cucurbitales, two from Rosales, two from Fabales, two from Malpighiales, two from Brassicales, one from Malvales, one from Solanales, two from Poaceae (Poales), one from A. comosus47 (Bromeliaceae, Poales), one from either M. acuminata48 (Zingiberales) or Phoenix dactylifera49 (Arecales), one from the Asparagales (from Asparagus officinalis50, Apostasia shenzhenica46, or Phalaenopsis equestris51), one from the Alismatales (either from S. polyrhiza52 or Z. marina53), one from Amborella, and one from G. biloba54. In total, 217 orthogroups based on anchor pairs and 142 orthogroups based on peak-based duplicates were collected.

The node joining the two WGD paralogues of N. colorata was then dated using the BEAST v1.7 package55 under an uncorrelated relaxed-clock model and an LG+G model with four site-rate categories. A starting tree with branch lengths satisfying all fossil prior constraints was created according to the consensus APG IV phylogeny1. Fossil calibrations were implemented using log-normal calibration priors on the following nodes: the node uniting the Malvidae based on the fossil Dressiantha bicarpellata56 with prior offset = 82.8, mean = 3.8528, and s.d. = 0.557; the node uniting the Fabidae based on the fossil Paleoclusia chevalieri58 with prior offset = 82.8, mean = 3.9314, and s.d. = 0.559; the node uniting the non-Alismatalean monocots based on fossil Liliacidites60 with prior offset = 93.0, mean = 3.5458, and s.d. = 0.561; the node uniting the N. colorata WGD paralogues with the eudicots and monocots based on the sudden abundant appearance of eudicot tricolpate pollen in the fossil record with prior offset = 124, mean = 4.8143 and s.d. = 0.562; and the root uniting the above clades with Amborella and then G. biloba with prior offset = 307, mean = 3.8876, and s.d. = 0.563. The offsets of these calibrations represent hard minimum boundaries, and their means represent locations for their respective peak mass probabilities in accordance with previous dating studies of these specific clades63 (see Supplementary Note 5.3 for an alternative setting of orthogroups).

A run without data was performed to ensure proper placements of the marginal calibration priors, which do not necessarily correspond to the calibration priors specified above, because they interact with each other and the tree prior64. Indeed, a run without data indicated that the distribution of the marginal calibration prior for the root did not correspond to the specified calibration density, so we reduced the mean in the calibration prior of the node combining the N. colorata WGD paralogues with the eudicots and monocots with offset = 124, mean = 4.4397, s.d. = 0.5 to locate the marginal calibration prior at 220 Ma62.

Markov chain Monte Carlo sampling for each orthogroup was run for 10 million steps, with sampling every 1,000 steps to produce a sample size of 10,000. The resulting trace files were inspected using Tracer v.1.555, with a burn-in of 1,000 samples, to check for convergence and sufficient sampling (minimum effective sample size of 200 for all parameters). In total, 263 orthogroups were accepted, and absolute age estimates of the node uniting the WGD paralogous pairs based on both anchor pairs and peak-based duplicates were grouped into one absolute age distribution, for which kernel density estimation and a bootstrapping procedure were used to find the peak consensus WGD age estimate and its 90% confidence interval boundaries, respectively. More detailed methods have been previously described39.

To identify the duplication events that resulted in the 2,648 anchor pairs detected in the genome of N. colorata, we performed phylogenomic analyses to determine the timing of the duplication events relative to the lineage divergences in Nymphaeales as described previously46. Protein-coding genes from 12 species were used, including eight species from Nymphaeaceae and one species from Cabombaceae in Nymphaeales, one species (I. henryi) from Austrobaileyales, plus Amborella and G. biloba. The phylogeny of the 12 species was obtained from Fig. 1d, and the branch lengths in K S units were estimated from 23 LCN genes (selected from the 101 LCN genes used in Fig. 1d, because only 23 are shared across all of the species studied) using PAML31 under the free-ratio model. OrthoMCL (v.2.0.9)65 was used with default parameters to identify gene families. Then, we removed 907 of the 2,648 anchor pairs with K S values greater than five. If the remaining anchor pairs fell into different gene families, thus indicating incorrect assignment of gene families by OrthoMCL, we merged the corresponding gene families and finally obtained 53,243 multi-gene gene families. Next, phylogenetic trees were constructed for a subset of 881 gene families with no more than 200 genes that had at least one pair of anchors and one gene from G. biloba. Multiple sequence alignments were produced by MUSCLE (v3.8.31)40 and were trimmed by trimAl (v.1.4)66 to remove low-quality regions based on a heuristic approach (-automated1).

We then used RAxML (v.8.2.0)67 with the GTR+G model to estimate a maximum-likelihood tree, starting with 200 rapid bootstraps followed by maximum-likelihood optimizations on every fifth bootstrap tree. Gene trees were rooted based on genes from G. biloba if these formed a monophyletic group in the tree; otherwise, mid-point rooting was applied. The timing of the duplication event for each anchor pair relative to the lineage divergence events was then inferred. In brief, internodes from a gene tree were first mapped to the species phylogeny according to the common ancestor of the genes in the gene tree. Each internode was then classified as a duplication node, a speciation node, or a node that has no paralogues and is inconsistent with divergence in the species phylogeny. The parental node(s) of a duplication node supported by an anchor pair were traced towards the root until reaching a speciation node in the gene tree. The duplication event that resulted in the anchor pair was hence circumscribed between the duplication node as the lower bound and the speciation node as the upper bound on the species tree. If the two nodes were directly connected by a single branch on the species tree, the duplication was thus considered to have occurred on the branch. To reduce biased estimations, we used the bootstrap value on the branch leading to a duplication node as support for a duplication event. In total, 497 anchor pairs in 473 gene families coalesced as duplication events on the species phylogeny, and duplication events from 254 anchor pairs in 246 gene families (or from 380 anchor pairs in 364 gene families) had bootstrap values greater than or equal to 80% (or 50%).

Floral scent measurement, gene identification, and functional characterization

We collected floral volatiles of N. colorata using a dynamic headspace sampling system and analysed them using gas chromatography–mass spectrometry (GC–MS) as previously described68. After 2 h of collection from the headspace of detached open flowers of N. colorata in a glass chamber (10 cm diameter, 30 cm height), volatiles were eluted from the SuperQ volatile collection trap using 100 µl of methylene chloride containing nonyl acetate as an internal standard. We then analysed samples using an Agilent Intuvo 9000 GC system coupled with an Agilent 7000D Triple Quadrupole mass detector. Separation was performed on an Agilent HP 5 MS capillary column (30 m × 0.25 mm) with helium as carrier gas (flow rate of 1 ml min−1). We applied splitless injections of 1 µl samples, injection temperature of 250 °C, an initial oven temperature of 40 °C (3-min hold) and a temperature gradient of 5 °C per min increase from 40 °C to 250 °C. Products were identified using the National Institute of Standards and Technology mass spectral database (https://chemdata.nist.gov).

A full-length cDNA of NC11G0120830 was amplified from the open flowers of N. colorata using reverse transcription PCR (RT–PCR), and cloned into pET-32a (MilliporeSigma). After confirmation by sequencing, NC11G0120830 was expressed in E. coli strain BL21 (DE3) (Stratagene) and the recombinant protein produced was purified using a modified nickel-nitrilotriacetic acid agarose (Invitrogen) protocol as previously reported69. For methyltransferase enzyme assays, we used both radiochemical and non-radiochemical reaction systems. The radiochemical reaction system (50 µl) was composed of 50 mM Tris-HCl, pH 7.8, 1 mM substrate, 1 µl 14C-S-adenosyl-l-methionine, and 1 µl of purified NC11G0120830. After 30 min of incubation at room temperature, 150 µl of ethyl acetate was added to extract the 14C-labelled reaction products. The extracts were counted using a scintillation counter (Beckman Coulter) to measure the activity of NC11G0120830. To determine the chemical identity of the reaction product, we performed non-radiochemical assays in which nonradioactive S-adenosyl-l-methionine was used as the methyl donor. The reaction product was collected by headspace solid-phase microextraction and analysed by GC–MS as previously described70.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.