The final dataset of macrofossils included 22 415 occurrences assigned to 443 genera, of which 263 are extinct and 180 are extant. The fossil occurrences were evenly distributed across the three plant groups described above: spore‐bearing plants (34%), nonflowering seed plants (34%) and angiosperms (32%). Angiosperms encompassed 58% of the diversity (with 257 genera) whereas the remaining genera were evenly distributed among the other two plant groups. The large majority of fossil occurrences included temporal ranges (minimum and maximum ages), which typically derive from the upper and lower temporal boundaries of the stratigraphic units to which fossils were assigned. In our dataset, temporal ranges spanned on average 12.7 Myr with a standard deviation of 12.6 Myr. In order to account for these dating uncertainties in the analyses of origination and extinction rates we considered the temporal ranges as uniform distributions and, from these, we randomly resampled the ages of each fossil occurrence following Silvestro et al . ( 2014b ). We generated 100 datasets using this randomization approach and ran the analyses on all replicates.

Analyses of diversification should ideally be performed on monophyletic groups, that is, clades in a phylogenetic sense. However, although monophyly assessment is straightforward in phylogenetic analyses based on molecular data, many questions remain on the phylogenetic placement of several fossil taxa. In addition, restricting the analyses to clades would require splitting the already fragmentary fossil data into smaller subsets, reducing analytical power and precision. We therefore categorized vascular plant fossil‐genera into three major groups, distinguished by their functional ecology, morphology and relationships (following Kenrick & Crane, 1997 ; Tomescu, 2009 ; Friis et al ., 2011 ; Hao & Xue, 2013 ): (1) spore‐bearing plants (i.e. early tracheophytes, rhyniophytes, lycophytes, sphenophytes, early euphyllophytes, progymnosperms, ferns and related plants, but excluding bryophytes); (2) nonflowering seed plants (i.e. ‘gymnosperms’ including Cycadidae, Ginkgoidae, Gnetidae, Pinidae, Bennettitales, Caytonia and seed ferns); and (3) flowering seed plants (i.e. angiosperms). Fossil occurrences were assigned to the different categories based on the supra‐generic taxonomic information provided by the PBDB and, when available, using the Index Nominum Genericorum (Smithsonian Institution, http://botany.si.edu/ing/ ). We matched the names of fossil‐genera in our dataset with those listed at The Plant List ( http://www.theplantlist.org/ ) to identify which taxa are extant.

Plant macrofossils (e.g. whole flowers, fruits, seeds, leaves) are commonly assumed to be the most reliable source of data to study plant evolution, because they usually contain many taxonomically useful characters and can therefore be more easily and reliably assigned to different plant groups (Wang et al ., 2010 ). By contrast, microfossils (pollen and spores) often lack a sufficient amount of diagnostic synapomorphies that facilitate taxonomic identification at the generic level. We therefore subsampled our original dataset to include only macrofossils (i.e. > 60% of the original dataset) in the analyses.

There are several practical reasons for avoiding species‐level analyses. First, excluding all fossils without assignation to the species level would reduce the full dataset by 47%. Second, > 50% of the species we compiled are known from a single occurrence, which makes it difficult to estimate accurate preservation rates regardless of the methodological approach used. Lastly, species‐level data should arguably be more sensitive to floristic variation due to regional endemism and fossil preservation biases. The interpretation and analytical incorporation of regional variations of species diversity is often problematic, as it could either derive from true biodiversity differences or reflect the fact that regional morphological differences are often described as new species in the absence of more thorough evidence. Thus, we consider that genus‐level analyses constitute a robust way of exploring the diversification dynamics of the plant fossil record.

It is well known that the plant fossil record is highly fragmentary, not least for plants (Cleal & Thomas, 2010 ). This is further complicated by the incomplete nature of fossilized plant organs, varying quality of preservation, and our incomplete knowledge on morphological evolution. These factors lead to difficulties in confidently assessing taxonomic affinity, and the existence of many fossils of uncertain taxonomic placement or ‘ incertae sedis’ . To minimize such problems we analyzed occurrence data at the genus level assuming that this level best captures the long‐term and global‐level dynamics of plant diversification (Rees, 2002 ; Wang et al ., 2010 ; Xiong & Wang, 2011 ; Xing et al ., 2014 ). An additional problem in the plant fossil record is that a single genus may be represented in several ‘fossil‐genera’ described from different plant remains or organs (e.g. in arborescent lycophytes; Cleal et al ., 2012 ). This might have the effect of overestimating plant diversity from some time intervals (e.g. the Carboniferous). Although our dataset includes several different types of plant remains (e.g. leaves, fertile and sterile axes, seed, fruits and roots), leaf samples represent > 70% of all the occurrences. Therefore, we considered the discrepancies between genera and fossil‐genera as unlikely to alter the major fluctuations of plant diversity.

This study is based on plant fossil occurrences downloaded from the Paleobiology Database (PBDB; http://paleobiodb.org/cgi-bin/bridge.pl ) on 6 March 2014. We performed multiple search queries in order to obtain the most comprehensive dataset of the plant fossils, on a global level. The approach that captured the maximum number of records involved querying the database for all described plant classes and orders individually. With this approach we obtained 35 666 entries, which spanned most of the Phanerozoic (the oldest occurrences dating back to c . 424 million yr ago (Ma)). We included all occurrences that were identified to the species level as well as those identified only at the genus level but without specific assignation (e.g. ‘ Psilophyton sp.’). We did not revise fossil identifications provided by the PBDB. Several culling procedures are normally used to avoid including incompletely identified taxa, such as searching and excluding entries identified with aff., cf., sp. or other modifications to binomial species names (Janevski & Baumiller, 2009 ). We initially explored the effect of excluding such records from the analyses, but because we found no differences in genus‐level diversity patterns we opted for using all available records. The dataset thus includes records from all lithologies, environments and continents that were freely available from the PBDB, including their estimated ages in millions of years (Myr). Raw data are available as supplementary materials (Supporting Information Tables S1, S2).

Diversification rate analyses

We carried out the analyses of plant diversification using a hierarchical Bayesian model to infer the temporal dynamics of origination and extinction (Silvestro et al., 2014b. This approach, implemented in the program PyRate (http://sourceforge.net/projects/pyrate/; Silvestro et al., 2014a), uses as input data all fossil occurrences that can be assigned to a taxon, in this case fossil‐genera, to jointly model the preservation and diversification processes. The preservation process infers the individual origination and extinction times of each taxon based on all fossil occurrences and on an estimated preservation rate (expressed as expected occurrences per taxon per Myr). This process is important because the first and last appearances of a taxon are likely to underestimate the true extent of its lifespan (Liow & Stenseth, 2007). The origination and extinction rates are inferred from the estimated lifespans of taxa based on a birth–death process, which models the diversification dynamics underlying diversity changes and constitutes the prior over the origination and extinction times (Silvestro et al., 2014b). In summary, the main parameters jointly estimated by this approach are: preservation rate and its degree of heterogeneity; origination and extinction times; and origination and extinction rates, which can vary through time. All of these parameters are sampled from their posterior distributions using a Markov Chain Monte Carlo (MCMC) algorithm, yielding parameter estimates and credible intervals, here calculated as 95% highest posterior density (HPD) intervals.

Compared to the original model described by Silvestro et al. (2014a,b), we introduced some modifications that better comply with the type of data and the aims of this study, which is mainly focused on variation in origination and extinction at global scale and large temporal ranges. We used here an homogeneous Poisson process (HPP) of preservation instead of the original ‘hat’‐shaped nonhomogeneous process (NHPP; Liow et al., 2010; Silvestro et al., 2014b). Although we believe the NHPP robustly describes fossil abundance in well‐preserved marine or mammalian species, we consider the HPP as a more appropriate assumption in the case of genus‐level plant fossil data due to potentially sparse sampling. We also accounted for varying preservation rates across taxa using the Gamma model, that is, with gamma‐distributed rate heterogeneity (Silvestro et al., 2014b). Because of the large number of occurrences analyzed and of the vast timescale considered, we used here eight rate categories to discretize the gamma distribution, instead of the default four, to accommodate the potential for more variability of preservation rates across taxa.

Furthermore, we broke down the birth–death process into segments of time, defined by the epochs of the stratigraphic geological timescale (Cohen et al., 2013) and estimated origination and extinction rates within these intervals. We adopted this solution as an alternative to the algorithms implemented in the original PyRate release that jointly estimate the number of rate shifts and the times at which origination and extinction undergo a shift. The estimation of origination and extinction rates within fixed time intervals improved the mixing of the MCMC and allowed us to obtain an overview of the general trends of rate variation throughout a long timescale, that is, almost the entire the Phanerozoic, within reasonable computational effort (a few days on a computer cluster). We emphasize that both the preservation and the birth–death processes are still modeled in continuous time and are not based on boundary crossings. Thus, the origination and extinction rates are measured as the expected number of origination and extinction events per lineage per Myr. One potential problem in fixing a priori the number of rate shifts is overparameterization. We overcame this issue by assuming that the rates of origination and extinction are part of two families of parameters following a common prior distribution, with parameters estimated from the data using hyper‐priors (Gelman, 2004). Thus, we used Cauchy prior distributions centered at 0 and with scale (hyper‐)parameters estimated from the data in the MCMC. The use of hyper‐priors is a convenient way to reduce the subjectivity in defining the priors and to reduce the effective parameter space, thus limiting the risks of overparameterization (Gelman, 2004).

We ran PyRate for 5000 000 MCMC generations on each of the 100 randomly replicated datasets. After excluding the first 20% of the samples as burnin phase, we combined the posterior estimates of the origination and extinction rates across all replicates and used them to generate rates‐through‐time plots. Rates of two adjacent intervals were considered significantly different when the mean of one lay outside of the 95% HPD of the other, and vice versa. We looked at the marginal posterior distributions of origination and extinction rates through the largest extinction events documented in geological history, the so‐called ‘Big Five’ mass extinction events in life's history (Jablonski, 2005; McElwain & Punyasena, 2010). In particular, we examined the diversification dynamics at four of the ‘Big Five’ events (the Frasnian–Famennian, Permian–Triassic, Triassic–Jurassic and Cretaceous–Paleogene boundaries; see the 1 section; Figs 1, 2). The Ordovician–Silurian extinction event occurred before the appearance of vascular plants and was therefore not considered. We focused on the magnitude of rate changes, their statistical significance and the uncertainty around those estimates. We measured extinction as the rate of disappearance of fossil‐genera, which are described based on morphology and might not always represent phylogenetic lineages. Thus, it was not possible with our data to discern potential events of pseudo‐extinction (i.e. when two fossil‐genera describe a single lineage evolving to a modified form).

Figure 1 Open in figure viewer PowerPoint et al., 2014a b x‐axis) is given in million years ago (Ma). Solid lines indicate mean posterior rates, whereas the shaded areas show 95% HPD intervals. The diversification dynamics were estimated for vascular plants as a whole (a–c), spore‐bearing plants (d–f), nonflowering seed plants (g–i) and flowering seed plants (angiosperms; j–l). Net diversification rates (gray) are defined as origination minus extinction. The dashed lines indicate the major mass extinction events widely recognized: at the Frasnian–Famennian (FFB), Permian–Triassic (PTB), Triassic–Jurassic (TJB) and Cretaceous–Paleogene (KPB) boundaries. Plant silhouettes were obtained from Generic‐level diversification analysis of all vascular plants. Origination (blue) and extinction (green) rates were estimated using the Bayesian approach implemented in PyRate (Silvestro.,) within time bins defined as epochs of the geologic timescale (shown at the bottom). The timescale (axis) is given in million years ago (Ma). Solid lines indicate mean posterior rates, whereas the shaded areas show 95% HPD intervals. The diversification dynamics were estimated for vascular plants as a whole (a–c), spore‐bearing plants (d–f), nonflowering seed plants (g–i) and flowering seed plants (angiosperms; j–l). Net diversification rates (gray) are defined as origination minus extinction. The dashed lines indicate the major mass extinction events widely recognized: at the Frasnian–Famennian (FFB), Permian–Triassic (PTB), Triassic–Jurassic (TJB) and Cretaceous–Paleogene (KPB) boundaries. Plant silhouettes were obtained from http://phylopic.org

Figure 2 Open in figure viewer PowerPoint Origination and extinction rates of plants across major events of global biodiversity turnover, the Frasnian–Famennian (a), Permian–Triassic (b), Triassic–Jurassic (c) and Cretaceous–Paleogene (d) boundaries. The time of the events is provided in million years ago (Ma). These plots show the sampled posterior distributions of origination (blue) and extinction (green) rates that were estimated before and after each event. Posterior distributions were standardized to relative densities for graphical purposes. Arrows on the top‐right corner of panels indicate significant rate increases (pointing up) or decreases (pointing down) based on the 95% HPD of the two distributions.

We estimated the origination times of vascular plants, seed plants and flowering plants from the posterior distributions of the origination times of the groups analyzed, following Silvestro et al. (2014b). Based on the posterior samples of the ages of origin derived from all replicates, we identified the best‐fitting distribution by testing several distributions (uniform, normal, log‐normal, gamma and Laplace) under maximum likelihood. We emphasize that these estimates are derived after explicitly modeling the fossil preservation process, and accounting for the uncertainties around the estimated ages of the fossil occurrences. After choosing the best‐fitting distribution based on the Akaike Information Criterion (AIC; Akaike, 1974), we computed the estimated parameter values, which may be directly used for other studies as calibrations to date molecular phylogenetic trees using relaxed molecular clock models (e.g. Lepage et al., 2007). The estimated parameters and offset values are given following the parameterization used in BEAST (Drummond et al., 2012), one of the most popular programs for molecular dating. We implemented these fitting functions in the latest version of PyRate.