New insights on angiosperm crown age based on Bayesian node dating and skyline fossilized birth-death approaches

Domain: earth_science
Published: 2025-03-07
Editorial outcome: major_revision

Abstract
Despite considerable work in recent years, pinpointing the time when angiosperms originated has been challenging. However, the rapid development of molecular clock methodology has provided new tools to resolve this conundrum. In particular, the fossilized birth-death model establishes a rich interplay between molecules and stratigraphy by incorporating fossils explicitly into dating analyses. In this study, we apply Bayesian node dating and the skyline fossilized birth-death model, which differ in how the calibration is applied, to estimate the crown age of angiosperms. Node dating analyses with different calibration strategies show that the posterior distribution is strongly constrained by the effective prior at the node of crown angiosperms, dominated by the maximum age constraint. Using the skyline fossilized birth-death model, we reveal that assigning different priors for origin time resulted in similar crown ages for angiosperms. Moreover, the oldest fossils play a significant role in time estimates, and the dating results are robust to sampling assumptions of extant taxa. Our dating analyses indicate a largely Triassic crown age (255–202 Ma) for angiosperms, the period when mammals, dinosaurs, and squamate reptiles first appeared, and highlight the potential of morphological data to redefine the timeline of angiosperms.

Introduction
Angiosperms, also known as flowering plants, are the most diverse, widely distributed and adaptable group of all major lineages of plants. The emergence of angiosperms was one of the most important evolutionary events in the history of the Earth, forming the structural and energetic basis of the vast majority of extant terrestrial biomes and reshaping major ecosystems around the globe. In particular, the diversification of angiosperms is an important driver of the decline or expansion of other plant groups such as ferns and conifers1,2. The time of origin of angiosperms also has a direct bearing on the discussion of their geographical origin and on the understanding of their paleoenvironmental conditions3–5. A reliable evolutionary timeline of angiosperms is needed to improve our understanding of how angiosperms came to dominate terrestrial ecosystems, and to explain patterns of diversity in other groups (e.g., insects and amphibians).

It is widely accepted that the earliest unequivocal fossils assignable to angiosperms appear in the early Cretaceous of northern Gondwana, ~133–125 million years ago (Ma)4,6–8. However, many molecular clock studies have pushed the origin of living angiosperms (the crown age) back to the early Jurassic or even earlier9–11, much older than fossil evidence. This age difference between molecular and fossil data was referred as the “Jurassic gap”10, which hinders our understanding of angiosperm evolution. Importantly, current Bayesian dating analyses based on large numbers of nuclear genes are limited and the taxon sampling is relatively sparse, and there has been no consensus on the crown age of angiosperms (e.g., 139–136 Ma12, 199–167 Ma13, 247–154 Ma14, 255–222 Ma11, 267–187 Ma10).

Calibration is the key factor in molecular clock dating, typically based on the fossil record or biogeographic events, converting the distances from molecular data into absolute times and rates15. One mainstream approach, referred to as node dating (ND), uses probability distributions (e.g., uniform or Cauchy distribution) to calibrate one or more internal nodes in the tree to estimate the ages of the other nodes16,17, reflecting our knowledge of the amount of elapsed time. Each calibration density is specified by the user to accommodate uncertainty in the relationship of the node age to the age of the reference fossil, taking into account processes such as fossilization, recovery, and morphological evolution. Although the ND approach has been widely used in estimating angiosperm timescales11,12,18–21, there are still some concerns about calibrating divergence time. Firstly, only the oldest fossil can be useful in constructing the calibration prior for the node in ND analyses. However, the existence of different morphospecies for different parts of the plant makes it difficult to identify the appropriate oldest fossils, with pollen being an extreme example. Secondly, unlike the minimum bound in the calibration distribution, the maximum bound is much less restrictive. When the maximum age of a lineage divergence has not yet reached consensus, specifying different maximum constraints could produce widely varying ages22. Thirdly, the user-specified calibration priors typically deviate from the prior that has effect in dating (called the effective prior or induced prior), due to interactions among calibration densities and the underlying tree and clock priors. Even though the tree prior can be constructed conditional on the calibration densities, the truncation effect resulting from the requirement that ancestral nodes be older than descendant nodes, still contributes to the deviation23.

An attractive alternative approach, named tip dating or total-evidence dating, directly includes fossils as terminal taxa in the tree24, thus obviating the need to specify calibration densities that are often assigned based on arbitrary criteria. Unlike ND analyses, tip dating allows for inclusion of all available fossils from a specific group, rather than only referring to the oldest one, and thus can take advantage of abundant time information. As the tree prior in the tip-dating approach, the fossilized birth-death (FBD) model represents a significant improvement over the uniform tree model, since the uniform tree prior has no birth, death, and sampling rate parameters, and it is usually considered a vague or uninformative prior24–26. The FBD model provides a coherent mechanistic description of species diversification, fossil preservation, and taxa sampling processes that produced the tree27,28. Recent advances extend the FBD model to allow the speciation, extinction, and fossil recovery rates to vary over time in a piecewise-constant manner26,29, referred as skyline FBD (SFBD) model. Furthermore, the SFBD model was extended to accommodate diversified sampling for larger, species-rich clades30, which could correct overestimation of divergence times resulting from violations of the random-sampling assumption26,31. Moreover, the SFBD model with many time epochs has recently proved to be robust to violations of sampling assumptions30. Despite the enthusiastic reception of the FBD process25,32–35, no studies have so far used the FBD model to estimate the origin time of angiosperms.

In this work, we expand the taxon sampling to use more angiosperm fossils to infer the timeline of angiosperms. Specifically, we focus on two questions: (1) the longstanding debate on the crown age of angiosperms; and (2) determining differences in divergence time estimation between the SFBD and ND approaches. The SFBD approach we use here differs from total evidence dating (TED) because no morphological characters are used. Instead, it relies on a fixed tree topology including fossils, and the effective prior on divergence times is induced by the SFBD process. To investigate key factors on divergence time estimation, we explore several calibrating strategies in ND analyses; while in SFBD analyses, we evaluate the impacts of fossil sampling, sampling strategy, density of extant species, and priors on the origin time. Our study not only provides the foundation for a better understanding of the patterns of angiosperm diversification and the evolution of key traits, but also yields valuable insights for estimating divergence times using the FBD model.

Results
We reconstructed a well–supported phylogenetic tree using a concatenated supermatrix including 324 nuclear genes from 670 angiosperm taxa and 78 species of gymnosperms (Supplementary Fig. 1). Along the backbone of the phylogeny, Amborellales was strongly supported as the sister lineage to all other angiosperms, followed successively by Nymphaeales and Austrobaileyales (Supplementary Fig. 1). We found eudicots and magnoliids to be sister groups with full support, and Chloranthales was the closest relative of eudicots-magnoliids (Supplementary Fig. 1). A sister relationship between the Ceratophyllales and monocots also had high support (Supplementary Fig. 1).

As the full molecular dataset is computational intractable for the subsequent dating analyses, we subsampled taxa at the genus level of eudicots, monocots, and gymnosperms, but retained the diversity of all sampled families, resulting in 349 angiosperm species and 15 gymnosperm species as an outgroup.

Impacts of maximum bound on time estimates in the node dating approach
In the ND analyses, we firstly assessed the relative fit of two relaxed-clock models to our data using Bayes factors. The results indicated that the independent-rates model had a higher marginal likelihood than the autocorrelated-rates model (Supplementary Data 1). We then explored whether maximum bounds are responsible for the age disparity of angiosperms using the independent-rates model. Based on the 114 fossil calibrations, we employed five fossil calibration strategies (see Methods for details) to accommodate different maximum bounds of crown angiosperms, mesangiosperms, and eudicots. To have a better understanding of the general impact of priors on posteriors, we compared the user-specified priors, effective (induced) priors, and posteriors for these three nodes.

Strategies 1 and 2 inferred a Carboniferous origin of crown angiosperms (strategy 1: 366–313.2 Ma; strategy 2: 365–304.9 Ma; 95% highest posterior density (HPD) interval), a Carboniferous–Permian diversification of mesangiosperms, and a Permian–Jurassic origin of eudicots (Supplementary Figs. 2 and 3). When reducing the maximum bound to 247.2 Ma for crown angiosperms and mesangiosperms, we obtained younger crown ages of angiosperms (strategy 3: 255.9–219.7 Ma; strategy 4: 292.0–249.5 Ma) and mesangiosperms. We also noted that changing from a narrow uniform distribution (strategy 3) to a more diffuse Cauchy distribution (strategy 4) for crown eudicots not only varied its own age estimate, but also pushed the crown age of angiosperms to be older (Fig. 1, Supplementary Figs. 4 and 5). By further reducing the maximum constraint in strategy 5, we obtained a much younger age of 147.2–140.9 Ma for crown angiosperms (Supplementary Figs. 6 and 7). The HPD intervals were quite narrow for the three nodes, indicating that the calibrations were too restrictive and the maximum constraint might be too young.

Moreover, we found that the effective prior has a much older and narrower age range than the calibration range for crown angiosperms (Fig. 1). Particularly in strategies 1–4, the effective prior assigns less than 5% probability of age within 50 Myr of the oldest fossil. In extreme cases, most prior density was in the tail of the specified prior (e.g., strategy 2). These results indicated that without molecular data, the effective priors do not result in a crown age for angiosperms in the early Cretaceous. In other words, the effective prior structure excludes the user-specified prior assigned for the crown node of angiosperms. Analyses with molecular data revealed posterior distributions that nearly overlap with the effective priors for angiosperms (strategies 1–3), implying that there is a strong dependence of the posterior estimates on the effective priors. Interestingly, the posteriors estimated in strategies 4–5 deviated from their effective priors (Fig. 1 and Supplementary Fig. 7). Additionally, we found that varying the birth and death rates in the five strategies did not affect the time estimations (Supplementary Fig. 8). Overall, our results indicated that the divergence times of angiosperm lineages were highly sensitive to the maximum constraints of the calibrated nodes. There is not yet any consensus on the evolutionary history of angiosperms among these calibrating strategies.

Given that the phylogenetic relationships of the five groups (Ceratophyllales, Chloranthales, magnoliids, monocots and eudicots) of mesangiosperms are still uncertain10,11,20,36, we evaluated the impact of alternative phylogenetic hypotheses based on calibration strategy 4 (Supplementary Fig. 9A). The estimated divergence times of crown angiosperms, mesangiosperms, monocots, and eudicots inferred from different topological relationships were almost identical (Supplementary Fig. 9B). These results indicated that different phylogenetic relationships had little effect on the estimated ages of angiosperms, probably because there is equal support for the various relationships among the genes.

Impacts of extant and fossil sampling schemes on time estimates under the skyline fossilized birth-death approach
To make direct use of the fossil occurrences rather than relying on controversial node calibrations to date the divergence times, we applied the SFBD approach assuming either random or diversified sampling of extant taxa, using 166 fossils and 349 extant angiosperm species. The results showed that the crown age of angiosperms under the diversified sampling (254.9–201.6 Ma) was similar to the random sampling (253.8–206.6 Ma) (Fig. 2A; Supplementary Figs. 10 and 11). The crown ages of major groups of angiosperms produced similar results, for instance, the mesangiosperms were estimated at 208.9–163.3 Ma under diversified sampling and 208.2–175.4 Ma under random sampling (Fig. 2A). Similarly, for shallow phylogenetic levels, the age intervals across most families under both sampling strategies were consistent (Supplementary Fig. 12). The effective priors under the two sampling strategies were generally older and had wider intervals (random sampling: 277.3–146.1 Ma; diversified sampling: 249.4–132.1 Ma) than posteriors for the crown ages of angiosperms.

Meanwhile, we explored the effects of overestimated sampling fraction (at family level) on time estimations under the two sampling schemes, and found the effects were generally limited (Fig. 2B). To investigate whether the observed distribution of fossil ages was consistent with that implied by the SFBD model, we sampled fossil ages given posterior estimates of divergence times and SFBD rates. Overall, the generated fossils had wider age intervals than the observed fossil ages under diversified sampling and random sampling (Fig. 2C). In addition, two sampled stem-group angiosperm fossil ages do not fall within the observed age range, likely because the fossils close to the root are more sensitive to origin time. Despite unavoidable uncertainties, the sampled fossil ages were moderately correlated with the observed data (Fig. 2C), implying that the selected fossils are compatible with the SFBD model in use.

As the SFBD process was conditioned on the origin time (tor), we further considered the impact of the choice of different priors for the origin on divergence time estimation. The dating results showed that the crown age of angiosperms was much less sensitive to the maximum age on the tor than under ND approach (Fig. 3). When the maximum value of tor was specified as 247.2 Ma, the crown age of angiosperms was restricted to this interval (diversified sampling: 242.6–197.1 Ma; random sampling: 242.4–213.2 Ma); while using an offset-exponential prior with mean of 240.63 Ma and offset of 125 Ma, the estimated time was similar with the main analyses (i.e., tor ~U(125, 365.63)) (Fig. 3).

Furthermore, we explored the effects of fossil sampling on estimating divergence times using four sampling strategies. Our results showed that when only including the oldest fossils (strategy A), the crown-group angiosperm time is largely in agreement with that estimated using all fossils (random sampling: 239.2–177 Ma; diversified sampling: 253.7–194.1 Ma) (Table 1 and Supplementary Figs. 13 and 14). In contrast, the corresponding age estimate from strategy B was relatively young (random sampling: 194.9–138.6 Ma; diversified sampling: 191.3–151 Ma) (Table 1, Supplementary Figs. 15 and 16). Removing fossils from stem-group angiosperms (strategy C) led to dramatically younger estimates of the crown age of angiosperms (156.7–151.3 Ma) and the major lineages under random sampling, while slightly older estimates under diversified sampling (Table 1, Supplementary Figs. 17 and 18). When the outgroup taxa with fossils were included (strategy D), we observed a much older crown age of angiosperms (random sampling: 277.3–225.3 Ma; diversified sampling: 281.0–228.5 Ma) (Table 1, Supplementary Figs. 19 and 20). This is possibly due to different diversification and sampling patterns between the outgroup and angiosperms, which are not accommodated in the SFBD model.

Evolutionary timescale of angiosperms
Using the two dating approaches, our results indicate that divergence times of angiosperms are sensitive to calibrations. However, we found that for the major groups of angiosperms, the width of the HPD intervals was similar between SFBD and ND approaches (~50 Ma) (Fig. 3), implying that the two methods are equally capable of estimating precision. There is no doubt that both approaches are undergoing continual refinement and have their own advantages and weaknesses. It is clear that the crown age of angiosperms depends on underlying assumptions about the fossil record and that there is still no consensus on maximum bounds in the ND approach. Given this, and our finding that assignment of different priors for the origin in the SFBD model resulted in relatively consistent crown ages of angiosperms compared to the ND approach (Fig. 3), we considered using the date estimates across SFBD main analyses under two sampling schemes to interpret the evolutionary history of angiosperms.

The SFBD analyses supported that crown-group angiosperms originated at 254.9–201.6 Ma (Fig. 4 and Supplementary Figs. 21 and 22). The crown group mesangiosperms originated in the late Triassic to late Jurassic (208.9–163.3 Ma). Apart from Amborellales, Nymphaeales and Austrobaileyales (namely ANA grade) diverged at 249–192.3 Ma and 227.2–175.8 Ma, respectively. The estimated timeline suggests that crown-group magnoliids arose in the early Jurassic to early Cretaceous (176.6–137 Ma). The Piperales and Canellales diverged at 162.1–111.1 Ma, Magnoliales and Laurales split at 155.9–126.4 Ma, with all four orders of magnoliids originating in the Cretaceous. We inferred that the crown-group monocots also arose during the early Jurassic to early Cretaceous (183.8–142.9 Ma). Within the monocots, the Acorales were the earliest to diverge, followed by the Alismatales splitting at 150.7–111.2 Ma. Furthermore, several radiations have occurred among the commelinids during the Cretaceous (128.8–103.3 Ma: Poales; 122.7–94.7 Ma: Commelinales and Zingiberales; 98.8–76.4 Ma: Zingiberales) (Fig. 4 and Supplementary Figs. 21 and 22).

We inferred an age of 182.6–146.9 Ma for crown eudicots. The orders corresponding to ~88% (31/35) of the core eudicots originated from the Aptian to the K-Pg (Cretaceous-Paleogene) boundary (125–66 Ma) in the Cretaceous (Fig. 4). This age range corresponds closely with the extensive angiosperm fossil record of this period (Supplementary Information), which implies that the Cretaceous was an important period in shaping the morphological diversity of angiosperms. This period has also been reported as a critical time for the co-evolution of angiosperms with insects, birds, and mammals37–40, which are key components of contemporary terrestrial biodiversity. Among the nine species-rich families of beetles (>10,000 species), seven of them originated in the Cretaceous and feed on plants or decomposing plant materials and woody tissues41. The recent discovery of a remarkably well-preserved short-winged flower beetle from mid-Cretaceous amber offers direct evidence of pollen-feeding in a Cretaceous beetle42. This amber was associated with pollen aggregations and coprolites consisting mainly of pollen, confirming that diverse beetle lineages visited early-diverging lineages of angiosperms in the Cretaceous. By examining the evidence for matching dates of origination of crown angiosperms and key insect pollinator lineages, the study highlighted the importance of the K-Pg extinction event in the context of angiosperm and insect co-speciation37. These results not only contribute to the elucidation of the origin and evolution of angiosperms, but also provide a theoretical basis for understanding key questions such as the diversification of other lineages that have co-evolved with angiosperms.

Methods
Extant taxon sampling
We sampled 670 angiosperm taxa (representing 669 angiosperm species) including 504 eudicots, 121 monocots, 29 magnoliids, four Chloranthales (five transcriptomes of Chloranthales representing four species), two Ceratophyllales, and nine early-diverging angiosperms, representing 58 orders and 278 families. Another 78 gymnosperm taxa were used as an outgroup. Our data of the 144 genomes and 579 transcriptomes were obtained from the public databases, and other 25 transcriptomic data were downloaded from the NCBI SRA database and assembled by Trinity v2.4.064 with default settings (min_kmer_cov=2). The complete list of data information is provided in Supplementary Data 2.

Orthology inference
The orthology inference was based on identifying homologs followed by cleaning, aligning, and cutting homolog trees65. First, we selected the complete nuclear genomes of 19 representative species (Ginkgo biloba, Papaver somniferum, Amborella trichopoda, Nymphaea colorata, Ceratophyllum demersum, Vitis vinifera, Actinidia eriantha, Coffea canephora, Helianthus annuus, Arabidopsis thaliana, Eucalyptus grandis, Prunus persica, Musa schizocarpa, Oryza sativa, Asparagus officinalis, Zostera marina, Cinnamomum kanehirae, Liriodendron chinense, Piper nigrum) to generate putative orthologous groups (OGs).

The CD-HIT (-c 0.995 -n 5) was used to reduce amino acid redundancy66. Initial homology search was carried out using all-by-all BLASTP of 19 complete genomes with an E value cutoff of 10 and -max_target_seqs 1000. We used a hit_fraction cutoff of 0.5 and set “IGNORE_INTRASPECIFIC_HITS=True” to filter raw blast output. The filtered hits were used as input for MCL (MCLv14-137)67, with the E value cutoff set to 10−5 and an inflation value of 1.5. Clusters containing fewer than 19 species were further removed. Then, we built the 3503 homolog trees by performing the tree estimation, trimming, masking, and cutting deep paralogs (keep the parameter default). Finally, we used homologs that had no duplicated taxon and are one-to-one orthologs to prune all homologous gene trees and inferred a set of 902 orthologs. The transcriptomic and genomic data of the remaining 728 taxa were then incorporated into 902 OGs through Hmmsearch (HMMERv3.2.1)68.

Phylogenetic analyses
All 902 OGs were aligned using MAFFT v7.31069 with the option “-localpair -maxiterate 1000” and then trimmed using Gblocks 0.91b70, allowing 50% gap positions. To remove the short sequences, trimAL v1.471 was used for pruning with the option resoverlap 0.5 -seqoverlap 50. All alignments were discarded if the length was below 100 amino acids or species coverage was lower than 90%. Finally, 324 OGs from a diverse selection of 748 taxa were obtained for phylogenetic inference.

We constructed phylogenetic relationships using the concatenated supermatrix that was partitioned by gene in IQ-TREE v2.0.572. The best-fit substitution model for each partition was determined and branch support was assessed with 1000 ultrafast bootstrap replicates.

Divergence time estimation
Molecular dataset and fossil occurrences. Our divergence time analyses were conducted using two Bayesian methods: node-dating (ND) and SFBD approaches. We performed preliminary analyses of the concatenated dataset (324 OGs for 748 taxa), but such a large dataset is computationally intractable. Previous study showed that the time estimates were largely consistent in full compared to reduced data sets, and suggested selecting a subset of informative genes for molecular dating analyses9. “Clock-like” genes were defined as those evolving in a clock-like manner73, which can minimize errors associated with model mis-specification and have been used for molecular dating9,22,74. Therefore, we used the SortaDate74 to identify the clock-likeliness of each gene. According to the three criteria sequentially: “1=minimal root-to-tip variance, 3=least topological conflict against the species tree, and 2=tree length (discernible information content)”, the top 100 genes of 324 OGs were selected for molecular-clock dating analyses. In addition, we reduced taxa at the genus-level of eudicots, monocots, and gymnosperms, but retained the diversity of all sampled families. The final dataset included 27,326 amino acids representing 100 nuclear genes from 364 taxa (349 angiosperm species and 15 gymnosperm species as an outgroup) for the dating analyses.

Following best practice principles of fossil relationships75, we selected 186 fossils that could reliably indicate the occurrence of seed plant clades down to family level. These contained 160 angiosperm fossils, 6 stem angiosperm lineage fossils, 11 gymnosperm fossils, and 9 extinct seed lineage fossils. The detailed fossil information including fossil ages, locality, fossil type, and relationships with taxa in the tree, are provided in Supplementary Data 3 and Supplementary Information. Most of (not all) the fossils were cited from the previous study14.

Node dating analyses. For the ND approach, divergence times were estimated using the program MCMCTree in the PAML v4.9j76, based on a fixed tree topology of the 364 extant taxa. The 27,326 amino acids were treated as a single partition, and the substitution model used was LG+F+Γ4.

To evaluate the relative fit of two relaxed-clock models (i.e., autocorrelated-rates model and independent-rates model) to the datasets, we selected 32 b values to estimate the marginal likelihood of a given model using a “stepping-stone” algorithm in MCMCTree. The output of MCMCTree was parsed and the marginal likelihood for the chosen relaxed-clock model was calculated using the R package mcmc3r77.

We consistently used a time unit of 100 MY (million years) across the ND analyses. For the evolutionary rate, we used the better-fit independent-rates model which assumes the branch rates follow the same lognormal distribution. To specify a sensible prior on the mean rate, we calculated the genetic distance between Arabis alpina and Ginkgo biloba based on LG+F+Γ4 model using CODEML76. Fossil evidence suggested that the divergence time of the two species is ~307 Ma78, which indicated that the mean substitution rate was assigned to a gamma hyperprior G(2, 26) with 2/26=0.076 amino acid substitutions per site per 100 MY. The variance (rate-drift) parameter σ2 was assigned G(1, 10).

The tree prior follows a birth-death process. As MCMCTree does not support estimating the speciation and extinction rates, we had to use fixed values for them. Based on the diversification rate estimates from previous sudy79, we set the birth rate (λ) at 9.88 and the death rate (μ) at 2.77. In addition, we assessed the impact of varying birth and death rates (λ=0.988, μ=0.277) on the time estimates. The sampling fraction of extant species (ρ) was set to 0.001 (364 species compared with the 352,000 extant angiosperms species) (https://www.mobot.org/MOBOT/research/APweb/). The root age was constrained to be between 307 and 365.63 Ma, based on the oldest fossil of Acrogymnospermae78 and the first record of seeds in the form of preovules that satisfy the criteria of the seed habit13.

We applied 114 fossil calibrations. In most cases, the uppermost boundary of the stratigraphic interval of the fossil attribution can be used as the minimum age for the calibration distribution. Nodes 1–9 were set with uniform distributions with soft bounds (pL=pU=0.025). Truncated Cauchy distributions with hard minimum bounds (pL=1e−300) were used at nodes 11–12, 14–38 and 40–114 (Supplementary Data 3, Supplementary Figs. 2–6). We implemented five strategies to test the effects of different maximum bounds on nodes 10 (crown angiosperms), 13 (crown mesangiosperms), and 39 (crown eudicots). In strategy 1, we used uniform distribution with 365.63 Ma as the maximum age for three nodes, which is based on the first record of seeds in the form of preovules that satisfy the criteria of the seed habit13. In strategy 2, Cauchy distributions truncated at 125 Ma were used on nodes of crown angiosperms, mesangiosperms, and eudicots. In strategy 3, we changed the maximum age in strategy 1 from 365.63 Ma to 247.2 Ma for nodes of crown angiosperms and mesangiosperms, based on sediments devoid of angiosperm-like pollen below their first report in the middle Triassic80; and applied the maximum age of 128.63 Ma for node of crown eudicots, based on the maximum age of the oldest potential age of tricolpate pollen80. In strategy 4, the nodes of crown angiosperms and mesangiosperms have the same specified priors as in strategy 3, while crown node of eudicots was assigned a Cauchy distribution. In strategy 5, following previous study12, we used we used 139.35 Ma for nodes of crown angiosperms and mesangiosperms, and Cauchy distribution for node of crown eudicots. Furthermore, to explore the effects of different topologies on the estimation of angiosperm crown age, we selected four topologies published in recent years10,11,20,36 and estimated the timeline of angiosperms based on strategy 4.

The MCMC analyses were run for 5 million generations sampled every 500 generations after discarding 1 million generations as burn-in. We conducted each analysis twice to check consistency and confirmed the effective sample sizes (ESS) of all parameters were larger than 200 by Tracer v.1.7.181.

Skyline fossilized birth-death dating analyses. For the SFBD approach, we used the BDSKY v1.5.0 package for BEAST v2.7.529,82 to infer the divergence times.

When morphological characteristics of extant species and fossils are unavailable, setting both accurate and precise topological constraints of fossil placements should be beneficial for analyses on FBD trees48. Following best practice principles75, we scrutinized the placement of fossils used in our study. Instead of integrating over the uncertainties of fossil placement through MCMC, which is computationally intractable in our case, we fixed the tree topology, as well as the position of the fossils in the tree. Detailed information on the placement of the fossils in the tree were available in the Supplementary Information. The fossil ages were given uniform distributions according to their stratigraphic intervals (Supplementary Data 3) and were sampled as part of the divergence time estimation. As the gymnosperms likely have different diversification and sampling patterns from the angiosperms, which would violate the SFBD model assumption, we used 349 angiosperm species and 166 fossils (including 6 stem group of angiosperm fossils) in the main analyses. Time was divided into ten intervals equally spaced from the root to the present. The speciation, extinction, and fossil‑sampling rates were all assumed independent across intervals. We specified diffuse Beta(1, 1) prior for the turnover (μ/λ) and fossil sampling proportion (ψ/(ψ+μ)), and Exp(1) prior for the diversification rate (λ−μ) in each interval. The sampling fraction of extant species (ρ) was fixed to 0.001 as in the ND analyses, and the sampling strategy was assumed either random or diversified26,30. But we also set ρ=0.67 (sampling fraction at family level) to explore the effects of overestimated ρ on time estimation. Considering the relationships of many stem branches of angiosperms remain highly uncertain and should not be relied upon46, we conservatively conditioned the SFBD process on the origin time (tor) using the minimum age constraint at 125 Ma based on the oldest fossil of angiosperms. The maximum age constraint at 365.63 Ma for tor is according to the first records of seeds in the form of preovules that satisfy the criteria of the seed habit, which represents a soft maximum constraint on the divergence of crown seed plants13. Since the true origin of the clade must have taken place before the oldest fossil date from within the clade, but we do not know how much earlier the origin was than this fossil. Therefore, we assigned the tor a uniform distribution (125, 365.63), assuming that the origin time has the same probability in this interval.

To explore the impact of the choice of different priors for the tor on divergence time estimation, we performed two additional analyses: (1) the origin time was assigned to an offset-exponential prior with offset at 125 Ma and mean of 240.63 Ma; (2) the origin time was set to a uniform prior (125, 247.2). Furthermore, we employed four strategies to explore the effects of fossil sampling on estimating divergence times in the SFBD approach. In strategy A, we selected the oldest fossils from the clades containing fossils, forming a subset of 106 fossils. In strategy B, we excluded the oldest fossil from each clade and the remaining 60 fossils were used to perform the dating analyses. In strategy C, we excluded the 6 stem angiosperm fossils and used the remaining 160 fossils. In strategy D, the 11 gymnosperm fossils, 9 extinct seed lineages fossils, and 15 gymnosperm outgroup species (186 fossils and 364 extant species) were all used for dating. The prior for the origin time was changed to U(307, 365.63) to accommodate the inclusion of the outgroup.

Still, the amino acids were treated as a single partition, and the substitution model used was LG+Γ4 (BEAST does not support +F). We also used the independent lognormal relaxed clock model83 for the evolutionary rate, and the implementation has been optimized to improve its computational performance84. We also calculated the amino acid pairwise distance between Amborella trichopoda and Actinidia chinensis with LG+Γ4+F model by the package CODEML in PAML package v4.9j. According to Tricolpate pollen grain fossil placement, the divergence time between the two species was about 125 Ma, which indicated that the mean rate was assigned to G(2, 0.00085) with 2×0.00085=0.0017 substitutions per site per time unit (1 My) (note that the change of timescale and that BEAST used shape-scale parameterization while MCMCTree used shape-rate parameterization in which rate=1/scale). The standard deviation of the rate was given G(2, 0.05) prior with mean 0.1.

For each run, the MCMC was executed for 400 million generations with sampling every 10,000 generations and a discarded burn-in fraction of 0.2. We ran each analysis twice and checked consistency and ESS > 200 for important parameters with Tracer v.1.7.1. Maximum-clade-credibility (MCC) trees were summarized from the combined samples using TreeAnnotator v2.6.0. To investigate the dating results and for better comparison with the ND results, we pruned the fossils from the MCC trees to produce annotated trees containing only extant taxa.

For all the analyses, we explored the effective priors for the divergence times induced by the SFBD process (by fixing the likelihood to 1.0 in the MCMC). Besides, we also performed a simple predictive check for the fossil ages in the main analyses. Specifically, we fixed the divergence times in the MCC tree, the speciation, extinction, and fossil-sampling rates in each time interval were all to their posterior medians. All fossils were given U(0, 500) distributions for their ages, and were all assumed as tips (this treatment is because the MCMC sampler does not support moving fossil ancestors along the tree). The fossil ages were then sampled using MCMC and were compared with the actual ages used in the dating inference.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_3": "Paper summary:\nIn this paper, the authors contrasted two highly used methodological approaches (ND vs Tip dating) to date the phylogeny of\nangiosperms. I think the paper is interesting from both a methodological, as well as an empirical perspective. From a\nmethodological perspective, I think the main contribution of this paper is showing the robustness of tip dating approaches on\nage estimates, in comparison to node dating approaches (Fig 3). Additionally, this paper provides a good example of how\nthe prior distributions in the ND approach can affect the divergence time estimates---a patten that has been previously\ndiscussed in the literature. From an empirical perspective, the authors addressed the longstanding question of when\nangiosperms first appeared and how they diversified using an impressive dataset.\nDear authors,\nThank you so much for this revised version of the manuscript, and your great efforts to address our extensive comments in\nthe previous round. This version is definitively a great improvement on the previous one. I am particularly pleased with the\nincreased level of detail in the methods, and the explanation of the reasoning for methodological choices. I think the addition\nof figure 3, the more extensive discussion of your results, together with a clearer explanation of your goals in the last\nparagraph of the introduction allow the reader to much better understand the impact of this paper.\nI think that your paper is at the stage in which you don’t need to conduct any more analyses. Nonetheless, I have several\nsuggestions (mainly on figures and some text clarifications) that I think will help to clarify a few details of your paper.\nCongratulations on an outstanding improvement!\nMain comments:\n• The abstract needs to be rewritten. In particular, the first sentences (l. 16-21) feel confusing, e.g. what conundrum you are\ntalking about? Please highlight your results from fig 3! They are in my opinion the most novel result of your study.\n• Overall, I would prefer to have more detailed figure captions. I expect your paper to be of general interest, ranging from\nempiricists to more method-oriented public, and with the amount of analyses you conducted, I think the readers would\nbenefit from more details on the figure legends. For example, for each, make sure to clearly state if these are results of ND or\ntip-dating analyses. Fig 2C labels needs much more detail on what “sampled fossil” is.\n• Thank you for showing figure 1 with the same axis in your response! I understand that you want to show the overlap of\nposterior and induced prior in the strategy 5. But I strongly prefer the figure you show in your response over the one you have\nin your main text. I think it is misleading to have different axis for different strategies, and very much prefer to be able to\nvisually contrast the priors and posteriors across strategies. The fifth strategy priors are very different from the first 4\nstrategies, and I think it is necessary to highlight that. That said, I think there are other ways in which you could solve the\nissue of not seeing the detail of the distributions on strategy 5, for example by adding a supplemental figure. Please re-\nconsider using the figure you sent in your response as figure 1 in the main text.\n• Fig 2C. In the las round of review I commented a lot on figure 2C. With your much improved description of the methods, I\nnow understand what you did to obtain this figure. But when I look at the figure now, I find the x axis label “sampled fossil\nage” confusing/not clear enough. I think a better label would be something like “Ages from posterior predictive simulations”\nor “simulated ages”, and I would appreciate if you briefly clarified how you obtained these ages in the label of the figure. I am\nstill unsure about why showing the regression line in this plot vs the 1 to 1 line. It seems more intuitive to show the 1 to 1 line\nto me, since that would be the expectation under a theoretical perfectly fit model.\n• Fig. 3 is a great addition!\n• I think the last paragraph of the introduction improved considerably! I think you can still highlight more the two guiding\nquestions of the main contributions of your paper: (1) node vs sophisticated tip dating approach, and (2) longstanding\nquestion on angiosperms divergence times.\n• Like reviewer #2, I am also confused by the fact that you use the first record of pre-ovules as the lower bound for the origin\ntime prior. I think this is a very old boundary since as far as I understand there are multiple groups of seed plants outside of\nstem branch of angiosperms to which these ovules could have belonged. Could you please add a more detailed justification\nof this choice following your description in line 610?\nIn-line comments:\n• L.90 – I am not sure all readers will understand that you are talking about a time uniform tree prior, maybe briefly explain?\n• L. 108 – “The SFBD approach we use here differs from tip dating and total evidence dating”. I agree your study is not a total\nevidence dating, because you don’t use morphology to infer fossils placement, but I actually consider your study to be tip\ndating, since it uses a FBD tree prior.\n• L.140 – “ we employed five fossil calibration strategies to accommodate different maximum bounds of crown angiosperms,\nmesangiosperms, and eudicots”. Thank you for adding extensive information in the methods section on what prior\ninformation informed your different strategies. I think the detailed information is fine in Methods, but you should mention in\nthis part of the manuscript that your strategies are informed and described later (maybe point to the section)--- otherwise the\nreader is left confused about why you chose those strategies.\n• It was hard to follow the explanations on lines 149 to 151. Can you specify what panels in fig 1 are you referring to?\n• Line 175, why did you choose strategy 4?\n• Line 290. Typo on “When”\n• L. 334 “In the absence of morphological data, fossil placement is first constrained by the researcher, and then refined by the\nprior information from the FBD process” --- I don’t follow, what do you mean by refined? Isn’t fossil placement fixed in your\nanalyses?\n• L. 355. Please clarify this sentence “The impact of fossils in our study is similar to the node dating approach, as the prior\ndivergence times are basically informed by the SFBD process.” --- I don’t follow hoe the impact of fossils in node dating is\nthe same as in tip dating.\n• I think the manuscript requires an additional proof reading for syntax errors/typos. A few examples below, but make sure\nyou check all the text.\no L. 339, get rid of “more”?\no L. rephrase as “the most influence in divergence time estimation”?\no L. 365 – rephrase to: “ Given the absence of methods that account for variation….”?\no L. 369 – rephrase to “potential bias when including outgroup sampling”?", "reviewer_4": "The manuscript by Ma et al. aims to investigate uncertainty around the age of the origin of angiosperms by comparing two\nmethods of dating the phylogeny: node dating, and tip dating through the skyline fossilised birth-death model. This is the first\ntime I have reviewed this paper, although I have seen the previous round of reviews and the authors’ responses to these.\nI see one major issue with the manuscript which is yet to be raised by the reviewers thus far. The implemented node dating\nmodel uses birth and death rates which are fixed (chosen by the authors) and constant throughout the evolutionary process\n(lines 552-556). The skyline FBD model was parameterised instead with diversification and turnover, with these values\nallowed to change across ten equally spaced time intervals (lines 600-604). It is not clear why a skyline model was used for\nthe FBD analysis, as using rates which were constant through time would make the assumptions of the models more\ncomparable. In my eyes, it is highly possible that allowing birth and death rates to change through time in one analysis, but\nnot the other, could be confounding the different results seen here, currently attributed to differences in the dating method. I\ncannot see the inferred diversification and turnover rates (or corresponding birth and death rates) from the skyline FBD\npresented, which means we currently do not know how much these rates are inferred to change through time in this model,\nand therefore how important this difference in assumptions is to the differences we see between the two models. The\nauthors need to present these results and explain how much (a) the difference in rate change assumptions, and (b) the\ndifference in dating method, makes to the results.\nI also find it surprising that there is no section in the manuscript which directly addresses the complex nature of the data\nbeing put into the analyses, particularly as this is a factor which may influence whether node dating or tip dating is chosen as\nthe most appropriate approach for dating a phylogeny. For plants specifically, the existence of different morphospecies for\ndifferent parts of the plant makes it much more difficult to identify appropriate fossils, with pollen being the extreme example\nof this. The paper comments on how node dating only uses the oldest fossil (line 73), but does not discuss how difficult it can\nbe to identify the appropriate oldest fossil, and how much that might contribute to the difficulties described here.\nI agree with (the previous) Reviewer 2 that not enough information is given concerning the origin prior implemented in the\nFBD analysis, and root prior implemented in the node dating analysis. Practically, setting an origin prior in an FBD analysis\nfollows the same logic as a node dating prior, and so direct comparison should be drawn between the age information\nplaced on the oldest nodes between the two analyses. The oldest fossil evidence should serve as the minimum age of the\nprior, while the maximum age should be clearly justified, as it often relies on uncertain assumptions regarding evidence\nexternal to the clade of interest. This is especially important as the age of the oldest nodes is the major result being\ncompared, and is likely highly dependent on the ages used in these priors (evidenced in lines 290-292). The root age of the\nnode dated analysis is currently stated but not explained (line 559), while the origin prior in the FBD mentions the age of\nTricolpate pollen, but does not explain why this pollen is considered a reasonable upper age boundary for the clade (lines\n607-610).\nI appreciate the difficulties involved in inferring the placement of the fossils in the\nphylogeny, but I think it is important for the manuscript to address the limitations of fixing their topological placement on the results. To what extent could potentially\nincorrect fossil placement impact the results of either analysis?\nI think this manuscript addresses an interesting problem but would need the above issues to be addressed before I am\nconvinced that the findings presented here are reliable."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #3 (Remarks to the Author):\nPaper summary:\nIn this paper, the authors contrasted two highly used methodological approaches (ND\nvs Tip dating) to date the phylogeny of angiosperms. I think the paper is interesting\nfrom both a methodological, as well as an empirical perspective. From a\nmethodological perspective, I think the main contribution of this paper is showing the\nrobustness of tip dating approaches on age estimates, in comparison to node dating\napproaches (Fig 3). Additionally, this paper provides a good example of how the prior\ndistributions in the ND approach can affect the divergence time estimates---a pattern\nthat has been previously discussed in the literature. From an empirical perspective, the\nauthors addressed the longstanding question of when angiosperms first appeared and\nhow they diversified using an impressive dataset.\nReply: Thank you very much for your positive comments.\nDear authors,\nThank you so much for this revised version of the manuscript, and your great efforts to\naddress our extensive comments in the previous round. This version is definitively a\ngreat improvement on the previous one. I am particularly pleased with the increased\nlevel of detail in the methods, and the explanation of the reasoning for methodological\nchoices. I think the addition of figure 3, the more extensive discussion of your results,\ntogether with a clearer explanation of your goals in the last paragraph of the\nintroduction allow the reader to much better understand the impact of this paper.\nI think that your paper is at the stage in which you don’t need to conduct any more\nanalyses. Nonetheless, I have several suggestions (mainly on figures and some text\nclarifications) that I think will help to clarify a few details of your paper.\nCongratulations on an outstanding improvement!\nReply: Thank you very much for your constructive suggestions, which are invaluable\nin enhancing our study.\nMain comments:\n• The abstract needs to be rewritten. In particular, the first sentences (l. 16-21) feel\nconfusing, e.g. what conundrum you are talking about? Please highlight your results\nfrom fig 3! They are in my opinion the most novel result of your study.\nReply: We have rewritten the abstract in our revised version. Specially, we have\nrephrased the first sentences to show the conundrum described here, as follows:\n“Despite considerable work in recent years, pinpointing the time when angiosperms\noriginated has been challenging. However, the rapid development of molecular clock\nmethodology has provided new tools to resolve this conundrum”. Furthermore, we\nhave highlighted the results from fig. 3 in the abstract: “Using the skyline fossilized\nbirth-death model, we reveal that assigning different priors for origin time resulted\nin similar crown ages for angiosperms.”\n• Overall, I would prefer to have more detailed figure captions. I expect your paper to\nbe of general interest, ranging from empiricists to more method-oriented public, and\nwith the amount of analyses you conducted, I think the readers would benefit from more\ndetails on the figure legends. For example, for each, make sure to clearly state if these\nare results of ND or tip-dating analyses. Fig 2C labels needs much more detail on what\n“sampled fossil” is.\nReply: We have supplemented more detailed figure captions. In particular, we have\nclearly stated if these are results of ND or tip-dating analyses for each figure.\nAdditionally, we have replaced “sampled fossil” with “simulated ages” as you\nsuggested, and added more detail on what “simulated ages” is as follows: “The age\nwas obtained from posterior predictive simulations by fixing the divergence times in\nthe MCC tree and assuming that all fossils were tip fossils.”\n• Thank you for showing figure 1 with the same axis in your response! I understand that\nyou want to show the overlap of posterior and induced prior in the strategy 5. But I\nstrongly prefer the figure you show in your response over the one you have in your\nmain text. I think it is misleading to have different axis for different strategies, and very\nmuch prefer to be able to visually contrast the priors and posteriors across strategies.\nThe fifth strategy priors are very different from the first 4 strategies, and I think it is\nnecessary to highlight that. That said, I think there are other ways in which you could\nsolve the issue of not seeing the detail of the distributions on strategy 5, for example by\nadding a supplemental figure. Please re-consider using the figure you sent in your\nresponse as figure 1 in the main text.\nReply: As you suggested, we have shown the figure 1 with the same axis in the main\ntext, and illustrated detailed information of strategy 5 in the supplementary figure 7\nin our revised manuscript.\n• Fig 2C. In the las round of review I commented a lot on figure 2C. With your much\nimproved description of the methods, I now understand what you did to obtain this\nfigure. But when I look at the figure now, I find the x axis label “sampled fossil ages”\nconfusing/not clear enough. I think a better label would be something like “Ages from\nposterior predictive simulations” or “simulated ages”, and I would appreciate if you\nbriefly clarified how you obtained these ages in the label of the figure. I am still unsure\nabout why showing the regression line in this plot vs the 1 to 1 line. It seems more\nintuitive to show the 1 to 1 line to me, since that would be the expectation under a\ntheoretical perfectly fit model.\nReply: We have replaced “sampled fossil ages” with “simulated ages”, and used 1 to\n1 line in figure 2C in the revised manuscript.\n• Fig. 3 is a great addition!\nReply: Thank you very much for your valuable suggestions.\n• I think the last paragraph of the introduction improved considerably! I think you can\nstill highlight more the two guiding questions of the main contributions of your paper:\n(1) node vs sophisticated tip dating approach, and (2) longstanding question on\nangiosperms divergence times.\nReply: We have pointed two guiding questions in the last paragraph of the\nintroduction, as follows: “Specifically, we focus on two questions: 1) the\nlongstanding debate on the crown age of angiosperms; and 2) determining\ndifferences in divergence time estimation between the SFBD and ND approaches.”\n• Like reviewer #2, I am also confused by the fact that you use the first record of pre-\novules as the lower bound for the origin time prior. I think this is a very old boundary\nsince as far as I understand there are multiple groups of seed plants outside of stem\nbranch of angiosperms to which these ovules could have belonged. Could you please\nadd a more detailed justification of this choice following your description in line 610?\nReply: We have added more detailed justification of this prior (lines 613-617). In\npractice, the origin time is often highly uncertain. There are indeed multiple groups\nof seed plants outside of stem branch of angiosperms, but the relationships of many\nstem branches of angiosperms remain highly uncertain and should not be relied upon\n(Sauquet et al. 2022). Moreover, the age of the oldest fossil sample (stem group of\nangiosperms: Petriellaea triangulate) in our dataset is included in this age interval\n(U(125, 365.63)). Thus, we conservatively used the currently oldest recognized fossils\nof crown angiosperms as lower bound (125 Ma).\nIn addition, we explored the impact on time estimates for crown angiosperms\nusing the age of the Petriellaea triangulate (242.7 Ma) as the lower bound. The\nresults indicated that the divergence time of crown angiosperms under two sampling\nstrategies (diversified: 253.7–199.6 Ma; random: 256.3-198.3 Ma) are in agreement\nwith the use of conservative prior (diversified: 254.9–201.6 Ma; random: 253.8–206.6\nMa).\nReference:\nSauquet, H., Ramírez-Barahona, S. & Magallón, S. What is the age of flowering\nplants? J. Exp. Bot. erac130 (2022).\nIn-line comments:\n• L.90 – I am not sure all readers will understand that you are talking about a time\nuniform tree prior, maybe briefly explain?\nReply: We have added a brief explanation of the uniform tree prior as follows: “As\nthe tree prior in the tip-dating approach, the fossilized birth-death (FBD) model\nrepresents a significant improvement over the uniform tree model, since the uniform\ntree prior has no birth, death, and sampling rate parameters, and it is usually\nconsidered a vague or uninformative prior.”\n• L. 108 – “The SFBD approach we use here differs from tip dating and total evidence\ndating”. I agree your study is not a total evidence dating, because you don’t use\nmorphology to infer fossils placement, but I actually consider your study to be tip dating,\nsince it uses a FBD tree prior.\nReply: We have rephrased this sentence as “The SFBD approach we use here differs\nfrom total evidence dating because no morphological characters are used. Instead, it\nrelies on a fixed tree topology including fossils, and the effective prior on divergence\ntimes is induced by the SFBD process.”\n• L.140 – “we employed five fossil calibration strategies to accommodate different\nmaximum bounds of crown angiosperms, mesangiosperms, and eudicots”. Thank you\nfor adding extensive information in the methods section on what prior information\ninformed your different strategies. I think the detailed information is fine in Methods,\nbut you should mention in this part of the manuscript that your strategies are informed\nand described later (maybe point to the section)--- otherwise the reader is left confused\nabout why you chose those strategies.\nReply: We have supplemented the related description in the revised version, as follows:\n“Based on the 114 fossil calibrations, we employed five fossil calibration strategies\n(see Methods for details) to accommodate different maximum bounds of crown\nangiosperms, mesangiosperms, and eudicots…”.\n• It was hard to follow the explanations on lines 149 to 151. Can you specify what\npanels in fig 1 are you referring to?\nReply: We have illustrated specific strategies in this sentence, as follows: “We also\nnoted that changing from a narrow uniform distribution (strategy 3) to a more diffuse\nCauchy distribution (strategy 4) for crown eudicots…”.\n• Line 175, why did you choose strategy 4?\nReply: Many studies including our dating results indicated that different\nphylogenetic relationships had little effect on the estimated ages (e.g., Barba-\nMontoya et al. 2018; Nie et al. 2020). Additionally, the 247.2 Ma was extensively used\nas the maximum bound for crown angiosperms (e.g., Barba-Montoya et al. 2018;\nMorris et al. 2018; Yang et al. 2020; Zhang et al. 2020; Su et al. 2021), which was\nbased on sediments devoid of angiosperm-like pollen below the first report in the\nMiddle Triassic. Compared to other maximum limits for crown age of angiosperms,\nthe 247.2 Ma (strategy 4) was considered more reliable.\nReferences:\n1. Barba-Montoya, J., Dos Reis, M., Schneider, H., Donoghue, P.C.J. & Yang, Z.\nConstraining uncertainty in the timescale of angiosperm evolution and the\nveracity of a Cretaceous Terrestrial Revolution. New Phytol. 218, 819–834 (2018).\n2. Nie Y, Foster CSP, Zhu T, Yao R, Duchêne DA, Ho SYW, Zhong B. Accounting\nfor uncertainty in the evolutionary timescale of green plants through clock-\npartitioning and fossil calibration strategies. Sys. Biol. 69, 1–16 (2020).\n3. Morris, J.L. et al. The timescale of early land plant evolution. Proc. Natl. Acad.\nSci. 115, E2274 (2018).\n4. Yang, L. et al. Phylogenomic insights into deep phylogeny of angiosperms based\non broad nuclear gene sampling. Plant Commun. 1, 100027 (2020).\n5. Zhang L, Chen F, Zhang X, et al. The water lily genome and the early evolution\nof flowering plants. Nature 577, 79–84 (2020).\n6. Su, D. et al. Large-scale phylogenomic analyses reveal the monophyly of\nbryophytes and Neoproterozoic origin of land plants. Mol. Biol. Evol. 38, 3332–\n3344 (2021).\n• Line 290. Typo on “When”\nReply: “When” >> “when”\n• L. 334 “In the absence of morphological data, fossil placement is first constrained by\nthe researcher, and then refined by the prior information from the FBD process” ---\nI don’t follow, what do you mean by refined? Isn’t fossil placement fixed in your\nanalyses?\nReply: The fossil placement is fixed in our analyses. To avoid confusion, we have\nrephrased these sentence as follows: “In the absence of morphological data,\nincorporating prior knowledge from previous studies to establish accurate\ntopological constraints of fossil placements has important implications for\ndivergence time estimation.”\n• L. 355. Please clarify this sentence “The impact of fossils in our study is similar to the\nnode dating approach, as the prior divergence times are basically informed by the SFBD\nprocess.” --- I don’t follow how the impact of fossils in node dating is the same as in\ntip dating.\nReply: We previously showed that the oldest fossils have an important role in\nestimating time in either ND or SFBD approaches. Considering that this section\nemphasizes on the performance of the SFBD model, we removed this sentence and\nrephrased this part as follows: “Our SFBD dating analyses show that the oldest\nfossils and the stem fossils exert the influence in estimating time, as the prior times\nare basically informed by the SFBD process.”\n• I think the manuscript requires an additional proof reading for syntax errors/typos. A\nfew examples below, but make sure you check all the text.\nReply: We have checked all the text to ensure that the revised manuscript has no\nsyntax errors/typos.\n• L. 339, get rid of “more”?\nReply: We have removed the word “more”.\n• L. rephrase as “the most influence in divergence time estimation”?\nReply: We have rephrased this sentence as follows: “In the absence of morphological\ndata, incorporating prior knowledge from previous studies to establish accurate\ntopological constraints of fossil placements has important implications for\ndivergence time estimation.”\n• L. 365 – rephrase to: “Given the absence of methods that account for variation….”?\nReply: We have rephrased this sentence to: “Given the absence of methods that\naccount for variation…”.\n• L. 369 – rephrase to “potential bias when including outgroup sampling”?\nReply: We have rephrased this sentence to: “…potential bias when including\noutgroup sampling.”\nReviewer #4 (Remarks to the Author):\nThe manuscript by Ma et al. aims to investigate uncertainty around the age of the origin\nof angiosperms by comparing two methods of dating the phylogeny: node dating, and\ntip dating through the skyline fossilised birth-death model. This is the first time I have\nreviewed this paper, although I have seen the previous round of reviews and the authors’\nresponses to these.\nI see one major issue with the manuscript which is yet to be raised by the reviewers\nthus far. The implemented node dating model uses birth and death rates which are fixed\n(chosen by the authors) and constant throughout the evolutionary process (lines 552-\n556). The skyline FBD model was parameterised instead with diversification and\nturnover, with these values allowed to change across ten equally spaced time intervals\n(lines 600-604). It is not clear why a skyline model was used for the FBD analysis, as\nusing rates which were constant through time would make the assumptions of the\nmodels more comparable. In my eyes, it is highly possible that allowing birth and death\nrates to change through time in one analysis, but not the other, could be confounding\nthe different results seen here, currently attributed to differences in the dating method.\nI cannot see the inferred diversification and turnover rates (or corresponding birth and\ndeath rates) from the skyline FBD presented, which means we currently do not know\nhow much these rates are inferred to change through time in this model, and therefore\nhow important this difference in assumptions is to the differences we see between the\ntwo models. The authors need to present these results and explain how much (a) the\ndifference in rate change assumptions, and (b) the difference in dating method, makes\nto the results.\nReply: Thank you very much for your comments. The two dating methods have\ndifferent model assumptions and parameters, for example, the fossil sampling rate is\ninvolved in the FBD model but not in ND approach. Thus, the time differences\nobserved between the two dating approaches were not solely due to variations in birth\nand death rates. For taxa that have undergone diversification for hundreds of million\nyears, fossils are usually patchy and non‑uniformly sampled, and the speciation and\nextinction rates change over time. In our data, we also observed that the inferred net\ndiversification and turnover rates dramatically varied over time (Fig. 1).\nIn addition, it has been reported that very old age estimates are produced under\nthe constant‑rate FBD model, which are not supported by the fossil record (e.g.,\nO’Reilly et al. 2015; Ronquist et al. 2016; May et al. 2021). One possible reason is\nthat assuming constant rates will impose an unrealistic prior belief on the\ndiversification and sampling processes, which may influence the divergence time\nestimates (Zhang et al. 2023). Allowing for flexible variation in speciation, extinction\nand fossil sampling rates over time could push the dramatically old age estimates\nobserved under the constant‑rate FBD model towards the present (Zhang et al. 2023).\nSimilarly, we also found that crown age of angiosperms estimated using the\nconstant‑rate FBD model (diversified sampling: 344.9–277.3 Ma; random sampling:\n344.7–283.1 Ma) were significantly older than the estimates obtained from the\nskyline FBD model (diversified sampling: 253.7–199.6 Ma; random sampling:\n256.3–198.3 Ma). Therefore, we analyzed our dataset using the skyline FBD model.\nFig. 1 Posterior estimates of the SFBD rates (median and 95% HPD CI) under random sampling\nand diversified sampling.\nReferences:\n1. O’Reilly J.E., dos Reis M., Donoghue P.C.J. Dating tips for divergence‑time\nestimation. Trend. Genet. 31:637–650 (2015).\n2. Ronquist F., Lartillot N., Phillips M.J. Closing the gap between rocks and clocks\nusing total‑evidence dating. Philos. Trans. R. Soc. London B Biol. Sci.\n371:20150136 (2016).\n3. May M.R., Contreras D.L., Sundue M.A., Nagalingum N.S., Looy C.V., Rothfels\nC.J. Inferring the total‑evidence timescale of marat-tialean fern evolution in the\nface of model sensitivity. Syst. Biol. 70:1232–1255 (2021).\n4. Zhang, C., Ronquist, F. & Stadler, T. Skyline fossilized birth–death model is\nrobust to violations of sampling assumptions in total-evidence dating. Syst. Biol.\n72, 1316–1336 (2023).\nI also find it surprising that there is no section in the manuscript which directly\naddresses the complex nature of the data being put into the analyses, particularly as this\nis a factor which may influence whether node dating or tip dating is chosen as the most\nappropriate approach for dating a phylogeny. For plants specifically, the existence of\ndifferent morphospecies for different parts of the plant makes it much more difficult to\nidentify appropriate fossils, with pollen being the extreme example of this. The paper\ncomments on how node dating only uses the oldest fossil (line 73), but does not discuss\nhow difficult it can be to identify the appropriate oldest fossil, and how much that might\ncontribute to the difficulties described here.\nReply: In the revised manuscript, we have rephrased this sentence to express the\ncomplex nature of the data as follows: “However, the existence of different\nmorphospecies for different parts of the plant makes it difficult to identify the\nappropriate oldest fossils, with pollen being an extreme example.”\nI agree with (the previous) Reviewer 2 that not enough information is given concerning\nthe origin prior implemented in the FBD analysis, and root prior implemented in the\nnode dating analysis. Practically, setting an origin prior in an FBD analysis follows the\nsame logic as a node dating prior, and so direct comparison should be drawn between\nthe age information placed on the oldest nodes between the two analyses. The oldest\nfossil evidence should serve as the minimum age of the prior, while the maximum age\nshould be clearly justified, as it often relies on uncertain assumptions regarding\nevidence external to the clade of interest. This is especially important as the age of the\noldest nodes is the major result being compared, and is likely highly dependent on the\nages used in these priors (evidenced in lines 290-292). The root age of the node dated\nanalysis is currently stated but not explained (line 559), while the origin prior in the\nFBD mentions the age of Tricolpate pollen, but does not explain why this pollen is\nconsidered a reasonable upper age boundary for the clade (lines 607-610).\nReply: We have added more detail to explain the origin prior implemented in the FBD\nanalysis (lines 613-617) and the root age of the node dating analysis (lines 558-561).\nIn practice, the origin time is often highly uncertain. There are indeed multiple\ngroups of seed plants outside of stem branch of angiosperms, but the relationships of\nmany stem branches of angiosperms remain highly uncertain and should not be\nrelied upon (Sauquet et al. 2022). Moreover, the age of the oldest fossil sample (stem\ngroup of angiosperms: Petriellaea triangulate) in our dataset is included in this age\ninterval (U(125, 365.63)). Thus, we conservatively used the currently oldest\nrecognized fossils of crown angiosperms as lower bound (125 Ma).\nIn addition, we explored the impact on time estimates for crown angiosperms\nusing the age of the Petriellaea triangulate (242.7 Ma) as the lower bound. The\nresults indicated that the divergence time of crown angiosperms under two sampling\nstrategies (diversified sampling: 253.7–199.6 Ma; random sampling: 256.3-198.3 Ma)\nare in agreement with the use of conservative prior (diversified sampling: 254.9–\n201.6 Ma; random sampling: 253.8–206.6 Ma).\nReference:\nSauquet, H., Ramírez-Barahona, S. & Magallón, S. What is the age of flowering\nplants? J. Exp. Bot. erac130 (2022).\nI appreciate the difficulties involved in inferring the placement of the fossils in the\nphylogeny, but I think it is important for the manuscript to address the limitations of\nfixing their topological placement on the results. To what extent could potentially\nincorrect fossil placement impact the results of either analysis?\nReply: Recently, it has been indicated that incorrect fossil placement can result in\nbiases in FBD model parameters as well as time estimates (Luo et al. 2020, 2023). In\npractice, setting both accurate and precise topological constraints of fossil\nplacements should be beneficial for analyses on FBD trees (Luo et al. 2023).\nFollowing best practice principles (Parham et al. 2012), we scrutinized fossils used\nin our study and were confident with fixed fossil placements. Thus, instead of\nintegrating over the uncertainties of fossil placements through MCMC, which is\ncomputationally intractable in our case, we fixed the tree topology, as well as the\nposition of the fossils in the tree.\nReferences:\n1. Luo, A., Duchêne, D.A., Zhang, C., Zhu, C. & Ho, S.Y.W. A simulation-based\nevaluation of tip-dating under the fossilized birth–death process. Syst. Biol. 69,\n325–344 (2020).\n2. Luo, A., Zhang, C., Zhou, Q., Ho, S.Y.W. & Zhu, C. Impacts of taxon-sampling\nschemes on Bayesian tip dating under the fossilized birth-death process. Syst. Biol.\n72, 781–801 (2023).\n3. Parham, J.F. et al. Best practices for justifying fossil calibrations. Syst. Biol. 61,\n346–359 (2012).\nI think this manuscript addresses an interesting problem but would need the above\nissues to be addressed before I am convinced that the findings presented here are\nreliable.\nReply: Thank you very much for your suggestions. In the revised manuscript, we\nhave added more detail to explain the origin prior implemented in the FBD analysis\nand the root age of the node dating analysis. Furthermore, we have supplemented\nsome sentences to express the complex nature of the data. The detailed response was\nprovided regarding the use of the skyline model in the FBD analyses, as well as\ninferring position of the fossils. We hope that our revisions and responses\ndemonstrate that our findings are reliable."}, {"round": 2, "reviewer_comments": {"reviewer_4": "I appreciate the efforts that the authors have gone to in addressing this round of reviews. However, I do not feel that my concerns have been adequately addressed.\nThe first point I raised was that the two methods implemented make different assumptions about whether speciation and extinction rates change through time. The authors have validly justified their use of the skyline FBD, in which these rates can change - I agree that this is very likely the better approach (as evidenced by the plotted skylines, which it would be great for the authors to add to the paper supplement as well as in response to me). However, their node dating approach assumed that speciation and extinction rates did *not* change through time, which the authors have agreed is a poor assumption. In this case, why should these results also be presented in the paper, if this model is being heavily violated? The authors' response to my comment includes the fact that \"the two dating methods have different model assumptions and parameters\", but the way these two methods are presented in the paper is as equally valid approaches. The abstract also includes the sentence \"Bayesian node dating and the skyline fossilized birth-death model, which differ in how the calibration is applied\", which implies that the treatment of sampling information is the only difference between the two approaches. In my eyes, the tone of the paper aims to present a \"fair\" comparison between the node dating and FBD approaches, but the lack on discussion of this specific difference in model assumptions seems problematic to me.\nThe third point I raised concerns the information given in the methods around the description of the origin and root age priors in the two analyses. I appreciate that the authors have added more explicit information about the specific age data that is being used to inform these priors. But at present, there is not enough explanation of the relationship between this data and the prior shape - the manuscript just reads that the prior is \"based on\" this data. We know that the true origin of the clade must have taken place before the oldest fossil date from within the clade, but we do not know how much earlier the origin was than this fossil. The use of the age of this known oldest fossil as a lower bound on the prior is reasonable, and I acknowledge that the authors investigated some sense of uncertainty in this lower bound, but how was the upper bound decided? The upper bound provides a vital constraint on the analyses and can be very influential on the results. In their response the authors state that \"the origin time is [often] highly uncertain\", which is exactly why it is important to explain to readers the decision which has been made, and explain how this decision might affect the results."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #4 (Remarks to the Author):\nI appreciate the efforts that the authors have gone to in addressing this round of reviews.\nHowever, I do not feel that my concerns have been adequately addressed.\nThe first point I raised was that the two methods implemented make different\nassumptions about whether speciation and extinction rates change through time.\nThe authors have validly justified their use of the skyline FBD, in which these\nrates can change - I agree that this is very likely the better approach (as evidenced by the plotted\nskylines, which it would be great for the authors to add to the paper supplement as well\nas in response to me). However, their node dating approach assumed that speciation and\nextinction rates did *not* change through time, which the authors have agreed is a poor\nassumption. In this case, why should these results also be presented in the paper, if this\nmodel is being heavily violated? The authors’ response to my comment includes the\nfact that “the two dating methods have different model assumptions and parameters”,\nbut the way these two methods are presented in the paper is as equally valid approaches.\nThe abstract also includes the sentence “Bayesian node dating and the skyline fossilized\nbirth-death model, which differ in how the calibration is applied”,\nwhich implies that the treatment of sampling information is the only difference between the two\napproaches. In my eyes, the tone of the paper aims to present a “fair” comparison\nbetween the node dating and FBD approaches, but the lack on discussion of this specific\ndifference in model assumptions seems problematic to me.\nReply: Thank you very much for your valuable suggestions. In our revised\nmanuscript, we have added the skyline plots to the Supplementary Fig. 23. We further\nestimated the timeline of angiosperms using ND approach in BEAST2 that allowing\nbirth and death rates to change across ten equally spaced time intervals. Consistent\nwith the constant birth-death rates analyses in MCMCTREE, different maximum\nconstraints conferred considerable influence on inferring the crown age of\nangiosperms. For instance, when the maximum constraint is set at 365.63 Ma under\nvariable birth and death rates, the crown age of the angiosperms is inferred to be\nPermian-Carboniferous (336.6-272.4 Ma). While reducing the maximum bound to\n247.2 Ma, we obtained younger divergence time of angiosperms (247.2-235.3 Ma).\nThus, the crown age of angiosperms mostly depends on underlying assumptions\nabout the fossil record in the node dating, regardless of whether the birth and death\nrates change through time. Based on these results, we think that the most important\ncontributor for the different time results obtained from FBD and ND is the\ndifferences in the dating method.\nWe have expanded discussion about speciation and extinction rates change\nthrough time in two model assumptions as follows: “Additionally, these two dating\nmethods make different assumptions about whether speciation and extinction rates\nchange through time in our study. We observed that the inferred net diversification\nand turnover rates varied over time in the SFBD analyses, as shown in the plotted\nskylines (Supplementary Fig. 23). While the implemented node dating model assumes\nbirth and death rates which are constant throughout the evolutionary process in the\nMCMCTREE. Although it has been shown that the most important factor affecting\ntime estimation in node dating is fossil calibration, future research could provide a\n‘fair’ comparison (i.e., allowing birth and death rates to change through time)\nbetween ND and FBD approaches.”\nThe third point I raised concerns the information given in the methods around the\ndescription of the origin and root age priors in the two analyses. I appreciate that the\nauthors have added more explicit information about the specific age data that is being\nused to inform these priors. But at present, there is not enough explanation of the\nrelationship between this data and the prior shape - the manuscript just reads that the\nprior is “based on” this data. We know that the true origin of the clade must have taken\nplace before the oldest fossil date from within the clade, but we do not know how much\nearlier the origin was than this fossil. The use of the age of this known oldest fossil as a lower bound on the prior is reasonable, and I acknowledge that the authors\ninvestigated some sense of uncertainty in this lower bound, but how was the upper\nbound decided? The upper bound provides a vital constraint on the analyses and can be\nvery influential on the results. In their response the authors state that “the origin time is\n[often] highly uncertain”, which is exactly why it is important to explain to readers the\ndecision which has been made, and explain how this decision might affect the results.\nReply: Thank you very much for your comments. We have supplemented more\ninformation to explain the relationship between the data and the prior shape. In the\nSFBD analyses, the maximum age constraint at 365.63 Ma is according to the first\nrecords of seeds in the form of preovules that satisfy the criteria of the seed habit,\nrepresenting a soft maximum constraint on the divergence of crown seed plants. The\nage constraint at 247.2 Ma is based on sediments devoid of angiosperm-like pollen\nbelow their first report in the Middle Triassic, representing a soft maximum\nconstraint on the divergence of crown angiosperms. Since some stem lineages of\nangiosperms were included in the FBD analyses, it is reasonable to use 365.63 Ma\nas the upper bound for the ‘origin time’. We further explored the prior distributions\n(uniform and offset-exponential) for the ‘origin time’ on time estimates. As you\nmentioned, the true origin of the clade must have taken place before the oldest fossil\ndate from within the clade, but we do not know how much earlier the origin was than this\nfossil. Therefore, we assigned the ‘origin (t )’ a uniform distribution, assuming\nor\nthat the probability that the time of origin in this interval is the same.\nWe have supplemented more detail to explain the prior for “origin time” as\nfollows: “…we conservatively conditioned the SFBD process on the origin time (t )\nor\nusing the minimum age constraint at 125 Ma based on the oldest fossil of\nangiosperms. The maximum age constraint at 365.63 Ma for t is according to the\nor\nfirst records of seeds in the form of preovules that satisfy the criteria of the seed habit,\nwhich represents a soft maximum constraint on the divergence of crown seed plants13.\nSince the true origin of the clade must have taken place before the oldest fossil date\nfrom within the clade, but we do not know how much earlier the origin was than this\nfossil. Therefore, we assigned the t a uniform distribution (125, 365,63), assuming\nor\nthat the origin time has the same probability in this interval.”."}, {"round": 3, "reviewer_comments": {"reviewer_4": "I thank the authors for responding to my last review, I am now satisfied that my concerns have been addressed."}, "author_response": ""}]

Decision letter