Phylogenomics and the rise of the angiosperms

<jats:title>Abstract</jats:title><jats:p>Angiosperms are the cornerstone of most terrestrial ecosystems and human livelihoods<jats:sup>1,2</jats:sup>. A robust understanding of angiosperm evolution is required to explain their rise to ecological dominance. So far, the angiosperm tree of life has been determined primarily by means of analyses of the plastid genome<jats:sup>3,4</jats:sup>. Many studies have drawn on this foundational work, such as classification and first insights into angiosperm diversification since their Mesozoic origins<jats:sup>5–7</jats:sup>. However, the limited and biased sampling of both taxa and genomes undermines confidence in the tree and its implications. Here, we build the tree of life for almost 8,000 (about 60%) angiosperm genera using a standardized set of 353 nuclear genes<jats:sup>8</jats:sup>. This 15-fold increase in genus-level sampling relative to comparable nuclear studies<jats:sup>9</jats:sup> provides a critical test of earlier results and brings notable change to key groups, especially in rosids, while substantiating many previously predicted relationships. Scaling this tree to time using 200 fossils, we discovered that early angiosperm evolution was characterized by high gene tree conflict and explosive diversification, giving rise to more than 80% of extant angiosperm orders. Steady diversification ensued through the remaining Mesozoic Era until rates resurged in the Cenozoic Era, concurrent with decreasing global temperatures and tightly linked with gene tree conflict. Taken together, our extensive sampling combined with advanced phylogenomic methods shows the deep history and full complexity in the evolution of a megadiverse clade.</jats:p>

section_content: Angiosperms are the cornerstone of most terrestrial ecosystems and human livelihoods 1, 2. A robust understanding of angiosperm evolution is required to explain their rise to ecological dominance. So far, the angiosperm tree of life has been determined primarily by means of analyses of the plastid genome 3, 4. Many studies have drawn on this foundational work, such as classification and first insights into angiosperm diversification since their Mesozoic origins [5] [6] [7]. However, the limited and biased sampling of both taxa and genomes undermines confidence in the tree and its implications. Here, we build the tree of life for almost 8,000 (about 60%) angiosperm genera using a standardized set of 353 nuclear genes 8. This 15-fold increase in genus-level sampling relative to comparable nuclear studies 9 provides a critical test of earlier results and brings notable change to key groups, especially in rosids, while substantiating many previously predicted relationships. Scaling this tree to time using 200 fossils, we discovered that early angiosperm evolution was characterized by high gene tree conflict and explosive diversification, giving rise to more than 80% of extant angiosperm orders. Steady diversification ensued through the remaining Mesozoic Era until rates resurged in the Cenozoic Era, concurrent with decreasing global temperatures and tightly linked with gene tree conflict. Taken together, our extensive sampling combined with advanced phylogenomic methods shows the deep history and full complexity in the evolution of a megadiverse clade. 

Flowering plants (angiosperms) represent about 90% of all terrestrial plant species 2 but, despite their remarkable diversity and ecological importance underpinning almost all main terrestrial ecosystems, their evolutionary history remains incompletely known. Since their Mesozoic origins 5, 10, 11, angiosperms have had a pervasive influence on the biosphere of Earth, shaping climatic changes at global and local scales 12, supporting the structure and assembly of biomes 13 and influencing the diversification of other organisms, such as insects, fungi and birds 14. 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. 

High-throughput DNA sequencing methods now enable us to reconstruct phylogenetic trees that broadly represent the evolutionary signal across entire genomes. Target sequence capture 15 has revolutionized plant phylogenetics by unlocking herbarium specimens as a source of sequenceable DNA 16, 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 3, 4, 7, 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. 8), have made target sequence capture consistently applicable across broad taxonomic scales, opening doors to collaboration and data integration 17. 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. 

Here, we present a nuclear phylogenomic tree that includes all 64 orders and 416 families of angiosperms recognized by the prevailing classification 18, using the Angiosperms353 (ref. 8) gene panel. Our sampling of 7,923 angiosperm genera (represented by 9,506 species) amounts to a 15-fold increase compared to previous work 9. 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.
section_name: Phylogenomics and the rise of the angiosperms

section_content: Our phylogenetic tree includes 58% of the approximately 13,600 currently accepted genera of angiosperms (Fig. 1 and Supplementary Table 1 ; ref. 2). 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 Article 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.
section_name: The angiosperm tree of life

section_content: Large 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
section_name: The phylogenomic challenge

section_content: 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 ). A 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 8, 9. Although gene duplication cannot be ruled out 19, the methods we used have been shown to be robust to the presence of paralogues 20. 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 21.
section_name: Article

section_content: Our 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 3, 4, 18, 22. 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 4, 9, 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 18 (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. 4 using plastid genomes (Supplementary Fig. 4 ). 

The 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 2. The relationships among main orders of asterids are stable 9, 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 23. 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 4, 24 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. 

By contrast to asterids, our findings in rosids conflict markedly with plastid-based evidence. First, we resolve Saxifragales, rather than Vitales 4, as sister to the remainder of rosids. In rosids, the fabid and malvid subclades, recovered as reciprocally monophyletic by plastid data 4, 22, 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 25. 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 18 ). 

Notwithstanding 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 9, 26. 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 27. These examples highlight The results illustrated are based on the young tree (maximum constraint at the root node of 154 Ma). See Extended Data Fig. 5 for results based on the old tree. a, Estimated net diversification rate through time (yellow, left y axis) and the level of gene tree conflict through time (blue, right y axis). Net diversification rates are estimated with a model that enables speciation rates to vary between time intervals; the line is the posterior mean and the yellow shaded area is the 95% highest posterior density. Gene tree conflict is calculated from the percentage of gene trees that do not share a congruent bipartition with each species tree branch, with the plotted value being the mean across all species tree branches that cross each 2. 5 Myr time slice. b, Cumulative percentage of extant orders and families that have originated through time. In both a and b, the background grey-scale gradient is the estimated percentage of extant lineages represented in the species tree through time (sampling fraction). 

the importance of areas of poor resolution as waymarkers to biological events meriting further study.
section_name: Phylogenetic insights from nuclear data

section_content: Our tree was analysed in combination with a dataset of 200 fossil calibrations (originally described in ref. 5, with modifications) to estimate divergence times and rates of diversification. Because the age of angiosperms is uncertain 28, 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 5, 28. 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 5 (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. 

With 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.
section_name: Time frame for angiosperm macroevolution

section_name: The diversification of angiosperms

section_content: We 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, r 2 = 0. 51). 

To 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 29, 30. Overall, however, gene tree conflict seems to be reliable corroborating evidence for investigating temporal patterns of angiosperm diversification.
section_name: Diversification linked to gene conflict

section_content: Our 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' 31, 32, is expected given the sudden emergence of diverse angiosperm fossils during the Early Cretaceous 11, [33] [34] [35]. Phylogenetic studies based on single or few genes have also implied that angiosperms diversified rapidly in the Early Cretaceous 7, [36] [37] [38]. 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. 

More 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 14, 39, 40. 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 33, 34,
section_name: Early burst of angiosperm diversification

section_content: highlighting that current molecular dating methods are unable to resolve the age of angiosperms 28.
section_name: Article

section_content: Following 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 41, 42 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. 43, they do not lead to observable global rate shifts across angiosperms.
section_name: A tapestry of lineage-specific histories

section_content: A 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 44. Instead, this finding is consistent with previously proposed positive feedbacks between increased diversity and increased rates of diversification in angiosperms 14, alongside more positive feedbacks, for example, between angiosperm and insect diversification 45, 46. Alternatively, global climatic cooling during the Cenozoic acting as a driver of angiosperm diversification could explain this finding 7, [47] [48] [49]. 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.
section_name: A Cenozoic diversification surge

section_content: The nuclear phylogenomic framework presented here is the result of an ongoing initiative to complete the tree of life for all angiosperm genera 50, 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 50. The standardized panel of nuclear genes in our dataset paves the way for more collaborations and data integration 17, 51, while the open availability of universal tools to sequence them (that is, Angiosperms353 probes 8 ) has made nuclear genomic data more accessible at relatively low cost. The accelerating uptake of this approach [52] [53] [54], which is readily applicable to herbarium collections 16, indicates that large volumes of data will soon become available for a wide range of applications in plant diversity, systematic and macroevolutionary research. Our 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.
section_name: Synthesis

section_content: Any 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 https://doi. org/10. 1038/s41586-024-07324-0.
section_name: Online content

section_content: As part of the Plant and Fungal Trees of Life (PAFTOL) Project at the Royal Botanic Gardens, Kew 50, 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. 

We produced target sequence capture data for 7,561 samples using the universal Angiosperms353 probe set 8 following established laboratory protocols 50, 56. We complemented our dataset with publicly available data for 2,054 species, sourced from the One Thousand Plant Transcriptomes Initiative 9 (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. 57,58) and target capture data (for example, see refs. 59,60). To standardize taxonomy and nomenclature, all species names and families were harmonized with the World Checklist of Vascular Plants 2 and orders with APG IV if possible 18.
section_name: Methods

section_content: Sequence recovery was carried out in two ways, depending on the type of input data. For recovery on the basis of raw reads, that is, Angio-sperms353 data or data mined from the SRA, we used HybPiper v. 1. 31 (ref. 61), embedded in a bespoke pipeline (https://github. com/baileyp1/ PhylogenomicsPipelines). Raw reads were trimmed using Trimmomatic 62 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 63 and refined by joining and trimming gene contigs to match coding regions using Exonerate 64. 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 20. 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. 50. 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. 50, using less stringent recovery thresholds. The average recovery per order is presented in Supplementary Fig. 9.
section_name: Sequence recovery

section_content: To 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.
section_name: Phylogenetic inference

section_content: The 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 ). 

The 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. 65) (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 66, gene trees were estimated using IQ-TREE v. 2. 2. 0-beta 67, 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 68. Before the second iteration, we identified long branch outliers using TreeShrink 69 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 70, enabling the full annotation option (-t 2), having first collapsed poorly supported nodes (ultrafast bootstrap ≤ 30%) in the input gene trees using Newick utilities 71. 

Order-level subalignments. For 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.
section_name: Backbone tree inference.

section_content: We produced global gene alignments by merging the order-level subalignments (before removal of gappy sites) and adding the orphan non-aligned sequences using MAFFT 65, 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 72 (Supplementary Table 5 ) and the average occupancy per gene per order is presented in Supplementary Fig. 11. 

We then estimated gene trees in Fasttree v. 2. 1. 10 (ref. 73), 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. 74), after collapsing branches with poor support (that is, those with support lower than 10% in the Shimodaira-Hasegawa test).
section_name: Global gene alignments and trees.

section_content: Divergence times were estimated by penalized likelihood in treePL 75, 76. 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 77. 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 78, 79. 

Fossil calibrations were based on the AngioCal fossil calibration dataset described in ref. 5. 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 5, 28. Both values are nonetheless considerably older than the oldest known crown group angiosperm fossils of around 127. 2 Ma (ref. 80). 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.
section_name: Divergence time estimation

section_content: At 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 2. 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.
section_name: Sampling extant lineages through time

section_content: Dated 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. 

An 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. 

Time-dependent diversification parameters (speciation and extinction rates) were also explicitly estimated across all angiosperms. These analyses were performed in RevBayes with the dnEpisodicBirthDeath function 81. 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 82. 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). 

Lineage-specific diversification rate estimation was performed in BAMM 83 and RevBayes. For analyses in BAMM, the setBammPriors function from the R package BAMMtools 84 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.
section_name: Diversification rate estimation

section_content: Simulations 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. 

Species 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. 

One-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. 

More methods, results and discussion are available (Supplementary Information; Supplementary Figs. 13-24 and Supplementary Table 8 ).
section_name: Simulations on gene tree conflict

section_content: The 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). . We also thank the Millennium Seed Bank Partnership for supporting access to samples. We acknowledge all national, state and regional authorities who authorized and facilitated the sourcing of these specimens. See also extended acknowledgements in the Supplementary Information. Extended Data Fig. 5 | Angiosperm-wide diversification and gene tree conflict through time. This is equivalent to Fig. 3 but for the old tree (maximum constraint at the root node of 247 Ma). a, Estimated net diversification rate through time (yellow, left y-axis) and the level of gene tree conflict through time (blue, right y-axis). Net diversification rates are estimated with a model that enables speciation rates to vary between time intervals; the line is the posterior mean and the yellow shaded area is the 95% highest posterior density.
section_name: Inclusion and ethics statement

section_name: Reporting summary Further information on research design is available in the

section_content: Gene tree conflict is calculated from the percentage of gene trees that do not share a congruent bipartition with each species tree branch, with the plotted value being the mean across all species tree branches that cross each 2. 5 Myr time slice. b, Cumulative percentage of extant orders and families that have originated through time. In both a and b, the background grey-scale gradient is the estimated percentage of extant lineages represented in the species tree through time ("sampling fraction").
section_name: Extended Data

section_content: April 2023

Data exclusions also excluded samples if more than three accessions were available for the same genus.
section_name: nature portfolio | reporting summary

section_content: All raw data and intermediate files are provided, making every step reproducible. All accessions used are listed in the Supplementary Table 1.
section_name: Reproducibility

section_content: Not applicable -randomization was not required.
section_name: Randomization

section_content: Not applicable -this work does not involve trials or controlled experiments. 

Did the study involve field work?

Yes No

Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.
section_name: Blinding

section_name: Acknowledgements

section_content: The PAFTOL project was funded by grants from the Calleva Foundation
section_name: Acknowledgements

section_name: Data availability

section_content: All raw DNA sequence data generated for this study are deposited in the European Nucleotide Archive under the following bioprojects PRJNA478314, PRJEB35285, PRJEB49212 and PRJNA678873. All analysed data and metadata are available in Zenodo at https://doi. org/10. 5281/ zenodo. 10778206 (ref. 55). The resulting trees and metadata are also available in GBIF (https://doi. org/10. 15468/4njn8b) and Open Tree of Life (https://tree. opentreeoflife. org/curator/study/view/ot_2304). The names used in this work match the World Checklist of Vascular Plants (https://doi. org/10. 34885/jdh2-dr22).
section_name: Data availability

section_content: The code used and developed to perform analyses is available in GitHub at https://github. com/RBGKew/AngiospermPhylogenomics and Zenodo at https://doi. org/10. 5281/zenodo. 10778206 (ref. 55).
section_name: Code availability

section_content: All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

-Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy All raw DNA sequence data generated for this study are deposited in the European Nucleotide Archive under the following bioprojects PRJNA478314, PRJEB35285, PRJEB49212 and PRJNA678873. All analysed data and metadata are available in Zenodo (https://doi. org/10. 5281/zenodo. 10778206). The resulting trees and metadata are also available in GBIF (https://doi. org/10. 15468/4njn8b) and Open Tree of Life (https://tree. opentreeoflife. org/curator/study/view/ot_2304). The names used in this work match the World Checklist of Vascular Plants (https://doi. org/10. 34885/jdh2-dr22).
section_name: Data Policy information about availability of data

section_name: Competing interests

section_content: Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 

Author contributions A. R. Z., T. C., A. A., S. Bellot, D. M. C., O. M. G., P. J. K., I. J. L., H. Sauquet, S. A. S., W. L. E., F. Forest and W. J. B were involved in conceptualization of this work. A. R. Z., T. C., A. A., S. Bellot, D. M. C., O. M. G., H. Sauquet, S. A. S., W. L. E., F. Forest and W. J. B. contributed to the methodology. A. R. Z., O. M., E. F., C. McGinnie, S. R. R., L. Simpson, J. J. C., R. S. C., S. D., L. Pokorny, G. M. A., K. G., K. P. H., A. Hoewener, A. -Q. H., E. M. J., I. A. B. S. K., I. L., D. A. L., E. J. Lírio, J. -X. L., P. Malakasi, N. A. S. P., T. S., J. V., G. K. A., R. L. A., M. S. A., M. A., M. J. A., J. A., W. A., J. -A. C., T. L. P. C., U. E., R. H. J. E., M. Estrella, M. F. F., P. L. F., M. G., L. J. G., S. W. G., J. H., T. Haevermans, J. C. H., O. H., T. R. H., K. K., E. A. K., M. A. K., F. L., C. J. L., D. -Z. L., S. L. -S., T. D. M., V. S. F. T. M., F. A. M., A. K. M., M. J. M., D. J. M., F. J. N., L. Palazzesi, J. P., D. S. P., O. A. P. -E., Y. P., G. M. P., R. K. R., R. R., H. Schaefer, A. S. -L., M. P. S., A. R. G. S., N. S., D. E. S., P. S. S., R. J. Soreng, P. F. S., S. C. K. S., A. H. T., T. M. A. U., J. M. V., J. J. W., T. C. W., Y. Y., S. Z. and I. J. L. reviewed the final manuscript. S. A. S., W. L. E., F. Forest and W. J. B. undertook supervision. P. J. K., I. J. L., F. Forest and W. J. B. acquired funding. V. M. A. B., P. J. K., I. J. L., F. Forest and W. J. B. were responsible for project administration.

section_content: The authors declare no competing interests.
section_name: Competing interests

section_content: The online version contains supplementary material available at https://doi. org/10. 1038/s41586-024-07324-0. Correspondence and requests for materials should be addressed to William J. Baker c, Equivalent to a, but focusing on the largest 25% of diversification rate increases. 

In a, b and c, the number of shifts is extracted from the maximum a posteriori shift configuration, the prior for the number of shifts is set to 10 and the background grey-scale gradient is the estimated percentage of extant lineages represented in the species tree through time ("sampling fraction").
section_name: Additional information Supplementary information

section_content: Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.
section_name: Reporting Summary

section_content: For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.
section_name: Statistics

section_content: The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

The statistical test(s) used AND whether they are one-or two-sided

Only common tests should be described solely by name; describe more complex techniques in the Methods section.
section_name: n/a Confirmed

section_content: A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e. g. means) or other basic estimates (e. g. regression coefficient) AND variation (e. g. standard deviation) or associated estimates of uncertainty (e. g. confidence intervals)

For null hypothesis testing, the test statistic (e. g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Give P values as exact values whenever suitable. 

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e. g. Cohen's d, Pearson's r), indicating how they were calculated

Our web collection on statistics for biologists contains articles on many of the points above.
section_name: A description of all covariates tested

section_content: Policy information about availability of computer code Data collection No software was used.
section_name: Software and code

section_content: Raw reads were trimmed using trimmomatic 0. 39. Data recovery was performed in HybPiper 1. 3. 1, which includes blast 2. 5. 0, exonerate 2. 4. 0, parallel 20200922 and SPADES 3. 13. 0. Phylogenetic inference was performed using ASTRAL-III 5. 7.
section_name: Data analysis

section_content: Not applicable

Note that full information on the approval of the study protocol must also be provided in the manuscript.
section_name: Ethics oversight

section_content: Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.
section_name: Field-specific reporting

section_content: For a reference copy of the document with all sections, see nature. com/documents/nr-reporting-summary-flat. pdf
section_name: Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences

section_content: All studies must disclose on these points even when the disclosure is negative.
section_name: Ecological, evolutionary & environmental sciences study design

section_content: We sampled nearly 8,000 genera and 353 genes to infer a phylogenetic tree for angiosperms, dated and calibrated it with a dataset of 200 fossils, and used this evolutionary time frame to study the diversification of the group.
section_name: Study description

section_content: We aimed to produce a dataset with at least one species per genus for >50% of the ca. 13,600 currently angiosperms genera. For genera with multiple species available, we retained up to three species, selecting primarily by phylogenetic representation followed by amount of data (number of genes and total length of recovery). Twelve species of gymnosperms were used as outgroups, one from each family. The final dataset contains 9,506 species, 7,923 genera, 416 families and 64 orders of angiosperms. It comprises 7,561 samples with data produced by this project and 1,963 samples mined from public repositories.
section_name: Research sample

section_content: We sourced samples from herbarium, living, DNA bank and and tissue bank collections.
section_name: Sampling strategy

section_content: We produced target sequence capture data using the Angiosperms353 probe kit. Total DNA was extracted using CTAB protocol and quantified by fluorometry. Average fragment size was assessed with 1% agarose gel. DNA extracts were diluted to 4 ng/ul. For extracts with high molecular weight, total DNA was fragmented using Covaris M220 Focused ultrasonicator. Genomic DNA libraries were prepared using NEBNext Ultra II DNA Library Prep Kit, following manufacturer's protocol at half volume, and with dual indexing. Libraries were normalised to 10 nM and pooled in equimolar amounts according to the average fragment size and taxonomic groups. Pools included, on average, 20-24 libraries. The pools were hybridised with Angiosperms353 probe kit. Hybridised pools were normalised to 7nM and combined for sequencing. Sequencing was carried out in Illumina platforms MiSeq and HiSeq. 

Timing and spatial scale Data production started in May 2016 and ended in December 2021.
section_name: Data collection

section_content: We excluded samples that failed both phylogenetic and barcode validation, as described in Baker (2022). Syst. Biol. 71: 301-319. We
section_name: Data exclusions