Phylogenomics and the rise of the angiosperms

1Flowering plants (angiosperms) represent about 90% of all terrestrial plant species but, despite their remarkable diversity and ecological importance underpinning almost all main terrestrial ecosystems, their evolutionary history remains incompletely known. Since their Mesozoic origins,,, angiosperms have had a pervasive influence on the biosphere of Earth, shaping climatic changes at global and local scales, supporting the structure and assembly of biomes and influencing the diversification of other organisms, such as insects, fungi and birds. The evolution of terrestrial biodiversity is thus inextricably linked with the macroevolution of angiosperms, which can only be shown using a robust and comprehensive tree of life. Reconstructing such a tree, however, is challenging because of the sheer diversity of angiosperms and the complex phylogenetic signal in their genomes.

2High-throughput DNA sequencing methods now enable us to reconstruct phylogenetic trees that broadly represent the evolutionary signal across entire genomes. Target sequence capture has revolutionized plant phylogenetics by unlocking herbarium specimens as a source of sequenceable DNA, thus removing the chief sampling bottleneck that has obstructed the completion of the tree of life. Although previous work on plants has relied primarily on the widely sequenced plastid genome,,, these technologies now allow us to tap into the evolutionary signal of the much larger and more complex nuclear genome. Universal nuclear probe sets, such as Angiosperms353 (ref. ), have made target sequence capture consistently applicable across broad taxonomic scales, opening doors to collaboration and data integration. As a result, opportunities now present themselves to address fundamental questions in plant evolutionary biology, such as the origin of angiosperms, the tempo and mode of their diversification and the classification of main lineages.

3Here, we present a nuclear phylogenomic tree that includes all 64 orders and 416 families of angiosperms recognized by the prevailing classification, using the Angiosperms353 (ref. ) gene panel. Our sampling of 7,923 angiosperm genera (represented by 9,506 species) amounts to a 15-fold increase compared to previous work. Leveraging a dataset of 200 fossil calibrations, we scale the tree to time, effectively capturing evolutionary divergences for all but the most recent 15% of angiosperm history. Although our tree broadly supports relationships predicted by previous studies primarily based on plastid data, it also shows previously unknown relationships and highlights some that remain intractable despite a vast increase in data. Gene tree conflict is tightly linked to diversification across the tree. We find evidence for high levels of conflict associated with an early burst of diversification, which is followed by an extended period of constant diversification rates underpinned by a tapestry of varied lineage-specific patterns. Diversification then increases in the Cenozoic Era, potentially driven by global climatic cooling. Our results highlight the fundamental role of botanical collections in reconstructing the tree of life to illuminate long-standing questions in angiosperm macroevolution.

4Our phylogenetic tree includes 58% of the approximately 13,600 currently accepted genera of angiosperms (Fig. 1 and Supplementary Table 1; ref. ). Together, the 7,923 genera encompass 85.7% of total known angiosperm species diversity. We produced data for 6,777 of these genera; before this study, 3,154 of these lacked publicly available genomic data, of which 393 lacked any form of DNA sequence data. For the remaining genera, data were obtained from public repositories. Sampling for this project was possible thanks to the collaborative effort of many biodiversity institutions from around the world, including 163 herbaria in 48 countries. More than one-third of species were sourced directly from herbarium specimens, some dating back nearly 200 years. Many phylogenetically problematic lineages with unconventional genome evolution were sampled, such as holoparasites, mycoheterotrophs and aquatics. Many of the species included are threatened and four are extinct (or extinct in the wild). The resulting tree of life presented here is one of the largest genomic trees generated yet for angiosperms as a whole.

5Large genomic datasets present challenges to phylogenetic inference. One issue is accurate homology assessment, which proved intractable across the full span of our dataset, even with the most advanced multiple sequence alignment methods. Another challenge is the efficient search of tree space based on gene matrices that have many more taxa than characters. We overcame both challenges with a divide-and-conquer approach (Supplementary Fig. 1). First, we computed a backbone species tree with sampling limited to five species per family (1,336 (15%) samples in total) and targeted to represent their deepest nodes (Supplementary Fig. 2). We used the backbone species tree to delimit taxon subsets for the construction of order-level gene alignments, which were then merged into global alignments. We then computed global gene trees from the global alignments, using backbone gene trees (inferred during the estimation of the backbone species tree) as topological constraints to reduce tree space while still letting gene trees differ from each other. The smaller number of samples in the backbone dataset permits a more thorough search of tree space, resulting in greater confidence at deeper nodes than could be achieved in an unconstrained global analysis. This approach allows a trade-off between comprehensive sampling and tree search robustness while accommodating putative discordance among gene trees. Finally, we used the global gene trees to generate a global species tree in a multispecies coalescent framework (Supplementary Fig. 3).

6A widespread concern in phylogenomic analysis is the presence of undetected gene copies. Our findings are unlikely to be affected by this because we used genes that have been selected to be mostly single-copy across green plants,. Although gene duplication cannot be ruled out, the methods we used have been shown to be robust to the presence of paralogues. In addition, a full assessment of orthologues was not computationally tractable but should be undertaken when methods become available to fully unravel the complexity of genome evolution at this scale.

7Our results broadly corroborate the prevailing understanding of angiosperm phylogenetic relationships, which rests on three decades of molecular systematic research largely built on data from the plastid genome,,,. We recover all main lineages of angiosperms, namely Amborellales, Nymphaeales, Austrobaileyales, Ceratophyllales and the three larger clades, monocots, magnoliids (including Chloranthales) and eudicots (Figs. 1 and 2). Although some of the relationships among those groups, such as the placement of Amborellales as sister group to all other angiosperms, are well-established and confirmed here, others, such as the placement of Ceratophyllales, which have been unstable in previous work,, remain inconclusive in our results. Despite the contrasting biological properties of the nuclear and plastid genomes (for example, size, copy number, mode of inheritance, recombination and evolutionary rate), which can lead to conflicting phylogenetic results, our findings largely support the mostly plastid-based phylogenetic classification of the Angiosperm Phylogeny Group (Extended Data Fig. 1). For example, 58 of the 64 now accepted orders and 406 of the 416 families are recovered as monophyletic (excluding artefacts; Supplementary Table 1). The most striking exception is the non-monophyly of Asteraceae, the largest angiosperm family comprising the sunflowers and their relatives. Our tree also confirms 85% of the relationships among families recovered by ref. using plastid genomes (Supplementary Fig. 4).

8The overall stability of established relationships is unevenly distributed across the tree, as observed in contrasting patterns in the main eudicot clades, the asterids and rosids, which account for 35% and 29% of angiosperm diversity, respectively. The relationships among main orders of asterids are stable, with a clade comprising Ericales and Cornales sister to all other asterids and the remaining 15 orders divided in two main clades (campanulids and lamiids), both long characterized by their contrasting floral ontogeny. Relationships contrasting with the status quo are mostly restricted to small orders, such as the paraphyly of Aquifoliales, Bruniales and Icacinales. These DNA-defined orders were consistently recovered as highly supported clades in plastome analyses, but they lack morphological cohesion. Given their placement in our phylogenetic tree and unique morphologies, these changes, although small, will alter our understanding of the evolution of asterids.

9By contrast to asterids, our findings in rosids conflict markedly with plastid-based evidence. First, we resolve Saxifragales, rather than Vitales, as sister to the remainder of rosids. In rosids, the fabid and malvid subclades, recovered as reciprocally monophyletic by plastid data,, are substantially rearranged into a grade of four orders subtending two well-supported sister clades, which we designate here as the recircumscribed fabids and malvids. The new fabid clade (Cucurbitales, Fabales, Fagales and Rosales) has long been characterized by symbiotic nitrogen fixation. In the new malvids (Brassicales, Celastrales, Huerteales, Malpighiales, Malvales, Oxalidales, Picramniales and Sapindales), Oxalidales is resolved as two independent lineages, the core emerging closer to Brassicales, Malvales and Sapindales, whereas Huaceae emerges in the position conventionally occupied by Oxalidales, that is, closer to Malpighiales and Celastrales (the former Celastrales–Oxalidales–Malpighiales (COM) clade).

10Notwithstanding the many well-supported confirmatory and new findings, some key relationships remain contentious and cannot be resolved by our data. These areas of high gene tree conflict often coincide with biological processes that confound phylogenetic inference. For example, the uncertain placements of eudicot orders Caryophyllales, Dilleniales and Gunnerales are probably impacted by key whole genome duplications,. The poor support for relationships among magnoliids, monocots, eudicots and Ceratophyllales might be explained by ancient hybridization events, such as that recently proposed for the origin of the monocots. These examples highlight the importance of areas of poor resolution as waymarkers to biological events meriting further study.

11Our tree was analysed in combination with a dataset of 200 fossil calibrations (originally described in ref., with modifications) to estimate divergence times and rates of diversification. Because the age of angiosperms is uncertain, we dated the tree with two different maximum constraints at the angiosperm crown node (154 and 247 million years ago (Ma), termed the young tree and old tree, respectively), which reflect realistic upper and lower bounds for the maximum age of this node,. These different constraints affected age estimates across angiosperms (Extended Data Fig. 2, Supplementary Fig. 5 and Supplementary Table 3). For example, in the young tree, stem node age estimates for Nymphaeales, Austrobaileyales and Ceratophyllales were 153, 152 and 152 Ma, respectively, whereas in the old tree the equivalent age estimates were 245, 244 and 243 Ma. Likewise, for larger clades such as magnoliids, monocots and eudicots, crown node age estimates were 151, 149 and 151 Ma in the young tree and 238, 237 and 241 Ma in the old tree. This range in age estimates is consistent with the most comprehensive comparable study (Extended Data Fig. 3) but our trees provide age estimates for a further 7,000 nodes. In subsequent analyses, we indicate if differing age estimates between the young tree and old tree cause substantially different interpretations of angiosperm diversification.

12With our sampling across angiosperms, we ensured that deeper branching events leading to extant lineages are comprehensively represented, while recognizing that extinct lineages are inaccessible to genomic methods. However, our dated trees are sparsely sampled at the species-level, meaning that branching events are incompletely represented towards the present, limiting diversification inferences in that time window. To address this, we developed a simulation-based approach to quantify the sampling fraction through time. For both dated trees, the lineage representation begins to drop substantially (below 75%) around 50 Ma (Supplementary Fig. 6). However, the most dramatic fall in lineage representation occurs in the most recent 20 Myr, in which it falls from around 50% to slightly more than 1% at present. Our investigation of angiosperm diversification should be interpreted with this broader context in mind. In particular, inferences in the most recent 20 Myr may be updated in the future with denser species sampling.

13We used our dated trees to reconstruct both diversification and gene tree conflict across a broad range of temporal and phylogenetic scales and investigate the relationship between them. We show that throughout angiosperm macroevolution, elevated gene tree conflict was tightly associated with elevated diversification. At a general level, this relationship is visible by simply comparing estimated diversification rates with gene tree conflict across all angiosperms through time (Fig. 3a). Meanwhile, in a branch-specific analysis using the temporal duration of branches as a proxy for the rate at which branches are diversifying, we also show that conflict and diversification rate are positively correlated (Extended Data Fig. 4) (P < 0.001, r2 = 0.51).

14To characterize the theoretical basis of this relationship, we simulated species trees with corresponding gene trees under different diversification scenarios in a multispecies coalescent framework. These simulations showed that gene tree conflict is positively correlated with diversification when caused by incomplete lineage sorting, assuming that effective population size is constant (Supplementary Fig. 7). Our empirical results are largely consistent with such a scenario. Other potential causes of gene tree conflict such as whole genome duplication and hybridization may also be associated with rapid diversification and have been recorded extensively throughout angiosperms,. Overall, however, gene tree conflict seems to be reliable corroborating evidence for investigating temporal patterns of angiosperm diversification.

15Our lineage-through-time (LTT) heatmap and diversification rate estimates through time both indicate an explosive early phase of diversification of extant lineages during the Late Jurassic and Early Cretaceous Periods (Fig. 2b and Fig. 3a). An early burst of angiosperm diversification, popularized as ‘Darwin’s abominable mystery’,, is expected given the sudden emergence of diverse angiosperm fossils during the Early Cretaceous,–. Phylogenetic studies based on single or few genes have also implied that angiosperms diversified rapidly in the Early Cretaceous,–. Our dated tree corroborates the existence of a distinct early burst of diversification, associated with high levels of gene tree conflict (Fig. 3a and Supplementary Fig. 8), further increasing our confidence in this finding.

16More than 80% of extant angiosperm orders originated during the early burst of diversification (Fig. 3b). Although not strictly comparable because of their subjective delimitation, orders represent the main components of angiosperm feature diversity, which have arisen rapidly after the crown node of angiosperms. In the young tree (Fig. 3), the early burst occurs during the Cretaceous, consistent with the hypothesis that a Cretaceous terrestrial revolution was triggered by the establishment of main angiosperm lineages,,. More controversially, the old tree places the early burst in the Triassic Period (Extended Data Fig. 5), which is dramatically at variance with the palaeobotanical record,, highlighting that current molecular dating methods are unable to resolve the age of angiosperms.

17Following the early burst, overall rates of diversification across angiosperms continued at a lower, constant pace for at least 80 Myr (Fig. 3a), during which time around three-quarters of all families originated (Fig. 3b). As expected, this phase of slower diversification was associated with lower levels of gene tree conflict. Despite the constancy of overall rates, diversification during this period was underpinned by a complex tapestry of lineage-specific patterns. This is illustrated by the LTT heatmap, which shows profound differences in diversification trajectories among orders (Fig. 2) and by the estimation of around 160 lineage-specific diversification rate shifts in angiosperms, most of which occur during this period. These rate shifts have a widespread phylogenetic distribution, with most orders containing at least one rate shift and many containing several nested shifts (Supplementary Table 4). The importance of nested rate shifts is highlighted extensively in discussions of evolutionary radiation, and underpins the continual response of diversification to dynamic extrinsic and intrinsic conditions. However, because these rate shifts are temporally scattered, as also shown by ref., they do not lead to observable global rate shifts across angiosperms.

18A second surge in angiosperm diversification occurred during the Cenozoic Era (Fig. 3a). The occurrence of this surge, despite the already high standing diversity of angiosperms at the time, suggests that diversification was unaffected by diversity-dependent processes, that is, the filling of available niche space as clades diversify. Instead, this finding is consistent with previously proposed positive feedbacks between increased diversity and increased rates of diversification in angiosperms, alongside more positive feedbacks, for example, between angiosperm and insect diversification,. Alternatively, global climatic cooling during the Cenozoic acting as a driver of angiosperm diversification could explain this finding,–. Importantly, an even larger Cenozoic surge would probably be inferred with increased sampling that addresses the under-representation of branching events in the recent time window. The temporal distribution of lineage-specific diversification rate shifts may offer some insight into the cause of the Cenozoic surge. Many of the largest diversification rate increases occur during the Cenozoic, whereas the number of diversification rate decreases declines markedly during this period (Fig. 4). These large rate increases may underpin the Cenozoic surge. The expansion of taxon sampling should be given priority to confirm these patterns.

19The nuclear phylogenomic framework presented here is the result of an ongoing initiative to complete the tree of life for all angiosperm genera, a milestone in our understanding of angiosperm evolutionary relationships. This study not only sheds light on much of the deep diversification history of the angiosperms but also lays foundations for future work towards a species-level tree. The standardized panel of nuclear genes in our dataset paves the way for more collaborations and data integration,, while the open availability of universal tools to sequence them (that is, Angiosperms353 probes) has made nuclear genomic data more accessible at relatively low cost. The accelerating uptake of this approach–, which is readily applicable to herbarium collections, indicates that large volumes of data will soon become available for a wide range of applications in plant diversity, systematic and macroevolutionary research.

20Our fossil-calibrated, phylogenomic tree enables a range of unique insights into broad-scale diversification dynamics of angiosperms, substantiating the early burst of diversification anticipated by Darwin while illuminating the complexity and conflict in the lineage histories underlying it. This sets the scene for future research, extending these investigations to shallower phylogenetic scales or digging more deeply into the data to discover the processes driving angiosperm diversification, such as genomic conflict, polyploidy, selection, trait evolution and adaptation. The challenges brought by the scale of this dataset and its ongoing expansion may also catalyse the development of methods which take full advantage of the global proliferation of genomic data.

21As part of the Plant and Fungal Trees of Life (PAFTOL) Project at the Royal Botanic Gardens, Kew, we assembled a nuclear genomic dataset consisting of newly generated data and data mined from public repositories. Our objective was to sample at least 50% of all angiosperm genera, with genera selected in a phylogenetically representative manner on the basis of published research. To avoid excessive imbalance in the tree, we included only one sample per species and a maximum of three species per genus. When several samples were available for the same species, we selected those with the largest amount of data, that is, more genes and a higher sum of gene length. For genera with several species available, the criterion for selection was primarily phylogenetic representation followed by amount of data. One species of each gymnosperm family was selected to form the outgroup, totalling 12 samples.

22We produced target sequence capture data for 7,561 samples using the universal Angiosperms353 probe set following established laboratory protocols,. We complemented our dataset with publicly available data for 2,054 species, sourced from the One Thousand Plant Transcriptomes Initiative (OneKP; 564 samples), annotated and unannotated genomes (151 samples) and the sequence read archive (SRA; 1,339 samples), the last including transcriptomes (for example, see refs.,) and target capture data (for example, see refs.,). To standardize taxonomy and nomenclature, all species names and families were harmonized with the World Checklist of Vascular Plants and orders with APG IV if possible.

23Sequence recovery was carried out in two ways, depending on the type of input data. For recovery on the basis of raw reads, that is, Angiosperms353 data or data mined from the SRA, we used HybPiper v.1.31 (ref. ), embedded in a bespoke pipeline (https://github.com/baileyp1/PhylogenomicsPipelines). Raw reads were trimmed using Trimmomatic to remove low-quality bases and short sequences. In HybPiper, reads were initially binned into genes using BLASTN and an amino acid target file as reference (Supplementary File 1). Individual genes were assembled de novo using SPADES and refined by joining and trimming gene contigs to match coding regions using Exonerate. For genes with paralogue warnings, only the putative orthologue as identified by HybPiper was used. Exclusion of genes with several copies per species has been shown to have negligible impact on species tree inference when it is performed under a multispecies coalescent framework, as described below. Conversely, the inclusion of several copies per species would have rendered our study computationally intractable. Gene sequences from assembled genomes and OneKP transcriptomes were recovered using custom scripts described in ref.. Briefly, the assembled sequences were searched against the target file mentioned above using BLASTN, selecting the best match for each gene and trimming it to the BLAST hit. For a few Angiosperms353 samples that represented the sole accession of their respective families (Ixonanthes reticulata, Mitrastemon matudae and Tetracarpaea tasmannica) and had poor recovery from HybPiper (that is, below 5 kilobase pairs (kb) in total sum of contig length), recovery was undertaken following ref., using less stringent recovery thresholds. The average recovery per order is presented in Supplementary Fig. 9.

24To analyse the dataset, we devised a divide-and-conquer approach. First, we computed a backbone tree, sampling up to five species per family, to test the monophyly of orders and to rigorously explore deep relationships. We used the backbone tree to identify groups (orders or groups of orders) for multiple sequence alignment, with the aim of producing refined subalignments among closely related taxa. Subsequently, the subalignments were merged into global gene alignments and global gene trees were inferred from these using the respective gene trees from the backbone analysis as constraints. Finally, we inferred a multispecies coalescent tree using the estimated gene trees. The inference pipeline is summarized in Supplementary Fig. 1.

25The samples for the backbone were selected so as to represent the crown node and deepest divergences in each family. For families with five or fewer samples (279 families), all samples were included. For those with more than five samples (156), we selected the best sample (most genes and longest sequence) of each consecutively diverging clade (based on published phylogenetic evidence and preliminary analyses of our own data), until five samples were included. To evaluate the extent to which sample selection might affect the backbone tree topology, we inferred 20 backbone replicates, randomly selecting five samples for each family with more than five samples (among the 50% best samples in terms of gene number and gene length recovered). We then summarized the trees to family level and computed Robinson–Foulds distances between the backbone and the 20 replicates (Supplementary Fig. 10).

26The phylogenetic reconstruction of the backbone involved up to two iterations of gene alignment and gene tree estimation, with an intermediate step of outlier removal. This was followed by species tree inference in a multispecies coalescent framework. In the first iteration, all sequences for a given gene were aligned using MAFFT v.7.480 (ref. ) (with ffnsi method, that is, --retree 2 --maxiterate 1000) and with the direction of the sequence adjusted (--adjustdirection). After removing sites with more than 90% missing data with Phyutility, gene trees were estimated using IQ-TREE v.2.2.0-beta, keeping identical sequences in the analysis (--keep-ident), setting the substitution model to GTR + G and estimating branch support with 1,000 ultrafast bootstrap replicates. Before the second iteration, we identified long branch outliers using TreeShrink in ‘all-genes’ mode and rerooting at the centroids of the trees. A second iteration of gene alignment, removal of gappy sites and gene tree estimation was performed for genes with outliers after the removal of outlier sequences. Subsequently, the resulting gene trees were summarized into a species tree using ASTRAL III v.5.7.3, a quartet-based species tree estimation method statistically consistent with the multispecies coalescent model, enabling the full annotation option (-t 2), having first collapsed poorly supported nodes (ultrafast bootstrap ≤ 30%) in the input gene trees using Newick utilities.

27For the order-level subalignments, most orders were analysed individually, following the same method described for the backbone. In some cases, smaller orders (fewer than 50 samples) were analysed together with larger ones if they formed monophyletic groups in the backbone. These groups are: (1) Commelinales with Zingiberales, (2) Dioscoreales with Pandanales, (3) Fagales with Fabales, (4) Columelliales, Dipsacales, Escalloniales and Paracryphiales with Apiales, (5) all magnoliids (Canellales, Laurales, Magnoliales and Piperales) and (6) all gymnosperms together (Cycadales, Ephedrales, Gnetales, Ginkgoales and Pinales). Conversely, orders emerging as non-monophyletic in the backbone were split into monophyletic subgroups as follows: (1) Cardiopteridaceae and Stemonuraceae separate from the rest of Aquifoliales, (2) Dasypogonaceae separate from the rest of Arecales, (3) Collumelliaceae separate from the rest of Bruniales, (4) Oncothecaceae separate from the rest of Icacinales and (5) Huaceae separate from the rest of Oxalidales. The groupings of samples used in the order-level subalignments are provided in Supplementary Table 1. Very small groups, comprising one or two samples (termed orphan sequences), were not included in subalignments and were incorporated directly in global analyses.

28We produced global gene alignments by merging the order-level subalignments (before removal of gappy sites) and adding the orphan non-aligned sequences using MAFFT, with up to 100 refinement iterations. This approach yields alignment across the order-level subalignments without disrupting the structure in the subalignments. The final gene alignments were cleaned by removing gappy sites. A summary of the alignments was produced with AMAS (Supplementary Table 5) and the average occupancy per gene per order is presented in Supplementary Fig. 11.

29We then estimated gene trees in Fasttree v.2.1.10 (ref. ), setting the model to GTR + G, using pseudocounts to avoid biases from fragmentary sequences and increasing search thoroughness (-spr 4, -mlacc 2 and -slownni). We used the gene trees from the backbone analysis to constrain the topology of each respective global gene tree. To avoid propagating error from the backbone analysis to the global analysis, we removed potentially misleading signal from the backbone gene trees before applying them as constraints. First, branches with bootstrap values below 80% were collapsed to avoid enforcing poorly supported relationships. Second, tips placed far from the rest of their order were algorithmically removed (but retained in global gene alignments). Once global gene trees were estimated, outlier long branches were removed using TreeShrink and the set of pruned gene trees was used to compute the global species tree using ASTRAL-MP v.5.15.5 (ref. ), after collapsing branches with poor support (that is, those with support lower than 10% in the Shimodaira–Hasegawa test).

30Divergence times were estimated by penalized likelihood in treePL,. This method is computationally efficient for datasets of this scale and typically estimates similar divergence times to more computationally intensive Bayesian analyses. The coalescent species tree topology was used as the input tree with molecular branch lengths estimated in IQ-TREE, on the basis of a concatenated alignment of the top 25 genes selected by SortaDate. Genes were selected by ranking their corresponding gene trees according to the number of congruent bipartitions with the species tree. We selected genes on this basis because high gene tree conflict leads to error in divergence time estimates,.

31Fossil calibrations were based on the AngioCal fossil calibration dataset described in ref.. We used an updated version of this dataset, referred to as AngioCal v.1.1 (Supplementary Table 6 and Supplementary File 2). Assigning fossil calibrations in this dataset to our tree topology led to 200 unique minimum age calibrations at internal nodes (Supplementary Table 7 and Supplementary Fig. 12). A maximum constraint of 154 or 247 Ma was used at the angiosperm crown node. These two values, respectively, represent a young and old constraint for the maximum age of the angiosperm crown node,. Both values are nonetheless considerably older than the oldest known crown group angiosperm fossils of around 127.2 Ma (ref. ). Both maximum constraints, in combination with all the minimum age constraints, were used to time-calibrate the species tree. Depending on the maximum constraint at the root node, these dated phylogenetic trees are referred to as young tree and old tree, respectively. For both the young tree and old tree, four analyses were performed in treePL, using smoothing values of 0.1, 1, 10 or 100. These different smoothing values assume high to low levels of among-branch substitution rate variation.

32At 1 Myr intervals from the root age of the dated phylogenetic trees to the present, we calculated how many angiosperm lineages would have been present in a hypothetical tree that sampled 100% of extant angiosperm species diversity. We used this to quantify the proportion of extant lineages incorporated by our phylogenetic trees through time (Supplementary Methods). To do this we simulated unsampled diversity on the dated trees: the species diversity of unsampled genera was simulated as a constant-rate birth–death branching process originating in the crown group of its respective family, whilst unsampled species diversity in sampled genera was simulated as a constant-rate birth–death branching process originating at the stem node of the relevant genus. The extant diversity of each simulated branching process was determined using the World Checklist of Vascular Plants. At each time interval, we then calculated the proportional difference between the number of lineages in our dated phylogenetic tree and the hypothetical fully sampled tree.

33Dated trees estimated with alternative smoothing values were very similar (Extended Data Fig. 2 and Supplementary Fig. 5), so diversification rate estimates were only performed with the dated trees estimated with a smoothing value of 10. By contrast, age estimates in the young and old trees differed markedly. Diversification rate estimates were therefore performed for both these dated trees. In each case, the dated trees were pruned such that there was a maximum of one tip for each genus.

34An initial analysis of diversification rates was performed by generating LTT plots as heatmaps for angiosperms as a whole, as well as for each order, with colours representing the steepness of each LTT curve at 5 Myr intervals. To calculate the steepness of the curve, we calculated the running difference between logarithmic corrected cumulative sums of lineages and applied Tukey’s running median smoothing to avoid excessive noise. For order plots, the cumulative sum starts at the first branching point, that is, order crown nodes.

35Time-dependent diversification parameters (speciation and extinction rates) were also explicitly estimated across all angiosperms. These analyses were performed in RevBayes with the dnEpisodicBirthDeath function. The smallest time windows in which rates were estimated were 5 Ma. However, larger windows were used toward the root of the tree such that there were at least 50 branching events in each time window. Three different models were used: equal rates of speciation and extinction across all windows; variable rates of speciation across windows but equal rates of extinction; and equal rates of speciation across windows but variable rates of extinction. Bayes factor comparison was used to compare models and offered strong support for the variable rate models but could not distinguish between the two variable rate models (Supplementary Information), indicating that they are probably from the same congruent set of models for the species tree. In subsequent discussion we primarily refer to results from the variable speciation rate model (for justification see Supplementary Information), although both variable rate models estimate similar patterns of net diversification rates through time (Supplementary Information).

36Lineage-specific diversification rate estimation was performed in BAMM and RevBayes. For analyses in BAMM, the setBammPriors function from the R package BAMMtools was used to define appropriate priors. Different sets of analyses were performed with the prior for the expected number of shifts set to either 10 or 100. These different prior settings had minimal effect on parameter estimates. Clade-specific sampling fractions were specified for each sampled family and a backbone sampling fraction of 1 was used. We therefore accounted for incomplete sampling within families alongside comprehensive sampling of the backbone of the tree. For analyses in RevBayes, the dnCDBDP function was used and the prior for the total number of rate shifts was set to either 10 or 100. Clade-specific sampling fractions cannot be specified with this function. Therefore, the sampling fraction was set to 1 meaning that estimates will become inaccurate toward the present because of unsampled within-family diversity.

37Simulations were based on a multispecies coalescent process. Each species tree contained 100 tips and was simulated as a birth–death branching process with time-dependent rates of speciation and extinction. In experiment 1, the extinction rate was always 0. The speciation rate was 0.75 species Myr−1 at times over 6 Ma, between 6 and 2 Ma the speciation rate was 0.075 species Myr−1 and less than 2 Ma the speciation rate was 0.75 species Myr−1. In experiment 2, the net diversification rates were the same as in experiment 1; however, in this case changes to the extinction rate led to the net diversification rate shifts. Therefore, for all time intervals, the speciation rate was 0.75 species Myr−1. At times over 6 Ma the extinction rate was 0 species Myr−1, between 6 and 2 Ma the extinction rate was 0.675 species Myr−1 and at times less than 2 Ma the extinction rate was 0 species Myr−1.

38Species trees with extinct lineages have extra complexities: first, changes in the extinction rate have a less direct impact on the duration of extant lineages in the species tree compared to changes in the speciation rate (Supplementary Information); and second, the effect of extinction is reduced at times close to the present. This causes shorter branches in the species tree, leading to the so-called ‘pull of the present’. We therefore performed a further analysis that was similar to experiment 2 but with no decrease in the extinction rate at the present. This offered insight into the effect of the ‘pull of the present’ on inferences of gene tree conflict and diversification rates and the relationship between these variables and the timing of rate shifts.

39One-hundred gene trees were simulated along the branches of the birth–death branching processes according to a multispecies coalescent process. For most experiments, the effective population size was 5,000. In one further experiment, which was otherwise the same as experiment 1, the effective population size was 50,000. For each simulated dataset, the degree to which the simulated gene trees exhibited conflicting topologies with the species tree was plotted through time (Supplementary Information). This enabled characterization of the relationship between gene tree conflict caused by incomplete lineage sorting and shifts in speciation and extinction rates in the species tree.

40More methods, results and discussion are available (Supplementary Information; Supplementary Figs. 13–24 and Supplementary Table 8).

41The research described here results from a highly inclusive, large-scale, international collaboration, which has actively encouraged the participation of many individuals from around the world. The authorship comprises many nationalities and is representative in terms of gender, career stage and career path. A total of 163 herbaria from 48 countries provided samples and/or house herbarium vouchers related to samples used in the study (see Acknowledgements). These samples originated from many countries, including Indigenous lands. We recognize the complex histories underlying all natural history collections and the global challenge we face in acknowledging them. We gave priority to recently collected samples and, as a result, most (85%) date from the postcolonial era (estimated here as 1970 onward). To share the benefits of our research, all data generated through this collaboration have been made publicly available before the submission of this work in several data releases, starting in 2019 (see Data availability).

42Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

43Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-024-07324-0.