The timescale of early land plant evolution.

Establishing the timescale of early land plant evolution is essential for testing hypotheses on the coevolution of land plants and Earth's System. The sparseness of early land plant megafossils and stratigraphic controls on their distribution make the fossil record an unreliable guide, leaving only the molecular clock. However, the application of molecular clock methodology is challenged by the current impasse in attempts to resolve the evolutionary relationships among the living bryophytes and tracheophytes. Here, we establish a timescale for early land plant evolution that integrates over topological uncertainty by exploring the impact of competing hypotheses on bryophyte-tracheophyte relationships, among other variables, on divergence time estimation. We codify 37 fossil calibrations for Viridiplantae following best practice. We apply these calibrations in a Bayesian relaxed molecular clock analysis of a phylogenomic dataset encompassing the diversity of Embryophyta and their relatives within Viridiplantae. Topology and dataset sizes have little impact on age estimates, with greater differences among alternative clock models and calibration strategies. For all analyses, a Cambrian origin of Embryophyta is recovered with highest probability. The estimated ages for crown tracheophytes range from Late Ordovician to late Silurian. This timescale implies an early establishment of terrestrial ecosystems by land plants that is in close accord with recent estimates for the origin of terrestrial animal lineages. Biogeochemical models that are constrained by the fossil record of early land plants, or attempt to explain their impact, must consider the implications of a much earlier, middle Cambrian-Early Ordovician, origin.

The establishment of plant life on land is one of the most significant evolutionary episodes in Earth history. Terrestrial colonization has been attributed to a series of major innovations in plant body plans, anatomy, and biochemistry that impacted increasingly upon global biogeochemical cycles through the Paleozoic. In some models, an increase in biomass over the continents, firstly by cryptogamic ground covers followed by larger vascular plants, enhanced rates of silicate weathering and carbon burial that drove major perturbations in the long-term carbon cycle (1, 2), resulting in substantial drops in atmospheric CO2 levels (3–6) (but see ref. 7) and increased oxygenesis (8). It also led to new habitats for animals (9) and fungi (10), major changes to soil types (11), and sediment stability that influenced river systems and landscapes (12). Attempts at testing these hypotheses on the coevolution of land plants (embryophytes) and the Earth System have been curtailed by a lack of consensus on the relationships among living plants, the timescale of their evolution, and the timing of origin of key body plan innovations (13). Although the megafossil record provides unequivocal evidence of plant life on land, the early fossil record is too sparse and biased by the nonuniformity of the rock record (13) to directly inform the timing and sequence of character acquisition in the assembly of plant body plans. Therefore, in attempting to derive a timescale for phytoterrestrialization of the planet, we have no recourse but to molecular clock methodology, employing the known fossil record to calibrate and constrain molecular evolution to time. Unfortunately, the relationships among the four principal lineages of land plants, namely, hornworts, liverworts, mosses, and tracheophytes, are unresolved, with almost every possible solution currently considered viable (14). In attempting to establish a robust timeline of land plant evolution, here we explore the impact of these conflicting phylogenetic hypotheses on divergence time estimates of key embryophyte clades.

Early morphology-based cladistic analyses of extant land plants suggested that the bryophytes are paraphyletic, but yielded conflicting topologies (15–17). Molecular phylogenies have been no more certain, with some analyses supporting liverworts as the sister to all other land plants (18), with either mosses (19–21) (Fig. 1), hornworts (22–27) (Fig. 1), or a moss−hornwort clade (28) (Fig. 1) as the sister group to the vascular plants. Variants on these topologies have been suggested, such as a liverwort−moss clade as the sister group to the remaining land plants (29) (Fig. 1). More recently, the debate has concentrated upon two hypotheses: hornworts as the sister to all other land plants (14, 30–34) (Fig. 1) or monophyletic bryophytes sister to the tracheophytes (14, 35, 36) (Fig. 1). Transcriptome-level datasets support both topologies (14), but sequence heterogeneity makes inferring relationships among these early land plants difficult (36).

Fig. 1.

The seven alternative hypotheses considered in the dating analyses. (A) Monophyletic bryophytes; (B) liverwort–moss sister clade to tracheophytes; (C) mosses, liverworts, and hornworts as successive sister lineages to tracheophytes; (D) a moss–liverwort sister clade to other embryophytes; (E) hornworts, mosses, and liverworts as successive sister lineages to tracheophytes; (F) mosses, hornworts, and liverworts as successive sister lineages to tracheophytes; and (G) a moss–hornwort sister clade to tracheophytes.

Here we attempt to establish a timescale of early land plant evolution that integrates over the contested topological relationships among bryophytes and tracheophytes. To achieve this, we constructed 37 fossil calibrations with minimum and soft maximum constraints, following best practice (37). This requires that calibrations are established on the basis of (i) a specific fossil specimen reposited in a publicly accessible collection, (ii) an apomorphy-based justification of clade assignment, (iii) reconciliation of morphological and molecular phylogenetic context of clade assignment, (iv) geographic and stratigraphic provenance, and (v) justification of geochronological age interpretation. Thus defined, these calibrations were combined with existing genetic data (14) in a Bayesian relaxed molecular clock analysis in which we also explored the impact of genetic dataset size and competing calibration strategies, as well as alternative substitution models, on divergence time estimates (Table 1). We find that topology and dataset size have minimal impact on age estimates, but slightly more variance in clade age estimates occurred when using alternative calibration strategies. We conclude that embryophytes emerged within a middle Cambrian to Early Ordovician interval and, regardless of topology, all four major lineages of land plants had diverged by the late Silurian. These dates are older than those used in the latest biogeochemical models (6, 8), and thus our results have implications for simulations of atmospheric chemistry and climate during the Paleozoic.

Table 1.

Summary of the analyses performed employing the seven alternative hypotheses, removal of the embryophyte constraints, and trimming dataset size

	Node distribution
	Uniform								Skew-t 850,000	Cauchy 850,000
Dataset	850,000	1.7 million	435,000	290,000	19,000	2,000	850,000	1.7 million†
Dataset no.
A Monophyletic	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
B Hornworts−sister	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
C Hornworts−liverworts−mosses	✓	X	X	X	X	X	X	X	X	X
D Liverworts−mosses−sister	✓	X	X	X	X	X	X	X	X	X
E Liverworts−mosses−hornworts	✓	X	X	X	X	X	X	X	X	X
F Liverworts−hornworts−mosses	✓	X	X	X	X	X	X	X	X	X
G Liverworts−sister	✓	X	X	X	X	X	X	X	X	X
Monophyletic (embryophytes only)	X	✓	X	X	X	X	X	X	X	X
Hornworts−sister (embryophytes only)	X	✓	X	X	X	X	X	X	X	X
Monophyletic (Chara–embryophytes)	X	✓	X	X	X	X	X	X	X	X
Hornworts−sister (Chara–embryophytes)	X	✓	X	X	X	X	X	X	X	X

All input topologies are based on the 290,718-nucleotide dataset, except for the Chara-embryophytes topology, which is based on the likelihood phylogeny of 1.7 million nucleotides.

A correlated model was used to estimate substitution rates on branches rather than the uncorrelated model used in all other analyses.

Results

Topology.

The competing hypotheses of relationships among bryophytes and tracheophytes all produce congruent age estimates across the phylogeny (Fig. 2 and Tables 2 and 3). Age estimates of key nodes (Embryophyta, Tracheophyta) are very similar regardless of the underlying topology (Fig. 2 and Tables 2 and 3). At the full range of uncertainty across topologies, the 95% highest posterior density (HPD) of ages for the embryophyte node ranges from the mid-Cambrian (Series 2; 515.2 Ma) to Early Ordovician (473.5 Ma) (Table 2), with the bulk of the distributions in the Cambrian (Fig. 2). There is a slightly higher variance in the estimated age of tracheophytes between the different topologies, but there is overlap in all of the 95% HPD age ranges (Fig. 2 and Tables 2 and 3). Estimates for the age of crown tracheophytes range from Late Ordovician (Katian; 450.8 Ma) to the latest Silurian (419.3 Ma).

Fig. 2.

Age estimates for the seven topologies used in analyses, highlighting the 95% HPD age uncertainty for embryophytes and tracheophytes. Age estimates are shown for (A) monophyletic bryophytes, (B) hornworts−sister, (C) hornworts−liverworts−mosses, (D) liverworts−mosses−sister, (E) liverworts−mosses−hornworts, (F) liverworts−hornworts−mosses, and (G) liverworts−sister.

Table 2.

The 95% HPD age estimates for of embryophytes and tracheophytes from divergence time analyses using the seven alternative topologies

Topology	Embryophytes, Ma	Tracheophytes, Ma
Dataset no.
A Monophyletic	514.8–473.5	450.8–431.2
B Hornworts−sister	515.2–482.1	450.8–430.4
C Hornworts−liverworts−mosses	515.2–483.3	450.7–419.3
D Liverworts−mosses−sister	514.9–477.7	450.8–431.1
E Liverworts−mosses−hornworts	515.1–480.8	450.7–427.9
F Liverworts−hornworts−mosses	515.1–483.2	450.7–428.5
G Liverworts−sister	514.9–478.4	450.8–428.2

Table 3.

The 95% HPD age estimates for named nodes in the analyses using the two main topologies of early land plants (monophyletic, hornworts−sister)

Clade	Monophyletic, Ma	Hornworts−sister, Ma
Viridiplantae	972.4–669.9	968.0–676.7
Streptophyta	890.9–629.1	875.4–637.4
Embryophyta	514.8–473.5	515.2–482.1
Bryophytes	506.4–460.3	N/A
Marchantiophyta	443.6–405.3	442.0–405.3
Marchantiopsida	354.9–228.0	357.9–228.0
Bryophyta	448.6–344.3	443.0–343.4
Tracheophyta	450.8–431.2	450.8–430.4
Lycopodiophyta	432.5–392.8	431.2–392.8
Euphyllophyta	437.6–402.2	435.7–402.2
Monilophyta	411.5–384.9	409.3–384.9
Spermatophyta	365.0–330.9	365.0–329.8
Acrogymnospermae	337.2–308.4	335.9–308.4
Pinopsida	301.3–172.4	302.8–172.1
Angiospermae	246.5–197.5	246.6–195.4
Mesangiospermae	180.4–139.5	177.6–139.2
Magnoliids	149.9–118.9	149.1–119.1
Piperales	103.7–51.4	106.7–50.6
Eudicotyledoneae	125.0–119.7	124.2–119.7
Monocotyledoneae	128.5–114.5	128.5–114.6

N/A, not applicable.

The two main hypotheses of early land plant relationships (monophyletic bryophytes and hornworts-sister) give congruent estimates for all nodes across the tree (Fig. 3 and Table 3). For example, the age estimates based on the two topologies are similar for Viridiplantae (972.4 Ma to 669.9 Ma), Streptophyta (890.9 Ma to 629.1 Ma), and Angiospermae (246.6 Ma to 195.4 Ma).

Fig. 3.

Detailed phylogenies showing the congruent age estimates produced using the monophyletic (A) and hornworts−sister (B) topologies.

Dataset Size.

Infinite site plots describe the relationship between clade age and uncertainty (95% HPD of clade age estimates). As the volume of sequence data increases, it is anticipated that clade age estimates should converge on a straight line, with residual dispersion reflecting uncertainty in calibrations that cannot be overcome by additional sequence data (38). We explored the impact of dataset size based on the monophyletic bryophytes topology, trimming the original dataset (1.7 million nucleotides) based on taxon completeness by 50%, 75%, 99%, and 99.9%. As expected, the resulting infinite sites plots reveal greater uncertainty (

Fig. 4.

Infinite site plots showing the effects of including more sequence data on the precision of age estimates. All ages are plotted using the monophyletic bryophytes topology with (A) datasets including all sites, and datasets trimmed so sequences are complete for (B) 50%, (C) 75%, (D) 95%, and (E) 99.9% of taxa.

Dating Strategies.

Across all alternative dating strategies, the age estimate for crown Embryophyta ranges from 583.1 Ma to 470.0 Ma (Fig. 5 and Table 4), which is larger than the range across the different topologies (515.2 Ma to 473.5 Ma). The greatest variance is seen when the embryophyte constraint is removed, resulting in older age estimates in the hornworts–sister topology, with an age distribution that stretches into the Proterozoic (to the middle Ediacaran), compared with the bulk of the distributions that fall within the Cambrian for all other age estimates (Fig. 5).

Fig. 5.

The estimated ages of embryophyte and tracheophyte divergence is more variable due to differences in modeling compared with differences in dataset size or topology. Using the monophyletic topology, the impact on age estimation was tested by using alternative strategies to model substitution rates, age constraints, and by excluding outgroups. An asterisk (*) denotes analysis performed on hornworts−sister topology.

Table 4.

95% HPD age estimates for embryophytes and tracheophytes in analyses after removing all nonembryophyte lineages, employing a correlated clock model, and applying different strategies for the shape of prior node age constraints (uniform unless stated)

Dating strategies	Embryophytes, Ma	Tracheophytes, Ma
Monophyletic no outgroup	515.0–473.6	450.8–430.1
Hornworts−sister no outgroup	515.1–478.6	450.8–430.1
Monophyletic correlation	514.0–470.0	450.9–440.7
Hornworts−sister correlation	514.4–475.0	450.9–439.8
Monophyletic no embryophyte constraint	535.3–475.7	450.8–431.4
Hornworts−sister no embryophyte constraint	583.1–489.2	450.8–431.7
Monophyletic cauchy	515.3–470.4	472.2–424.2
Hornworts−sister cauchy	534.0–471.4	463.4–423.2
Monophyletic skew-t	493.8–470.7	457.7–422.7
Hornworts−sister skew-t	497.3–471.1	444.8–422.4
Monophyletic (Chara–embryophytes)	514.9–476.6	450.9–436.7
Hornworts−sister (Chara–embryophytes)	515.2–484.1	450.9–434.5

There is greater variance when these uncertainties are used compared with the smaller variance seen on dating analyses using the alternative topologies.

We employed different parametric distributions (uniform, Cauchy, skew-t) to express the prior probability of divergence timing relative to the minimum and soft maximum constraints. This often has a dramatic impact on divergence time estimates (39–41); however, different prior distributions have minimal impact on age estimates for embryophytes. The largest difference is seen with the younger age estimates produced using the skew-t distribution (Fig. 5), but both the skew-t and Cauchy models produce younger mean estimates for embryophytes compared with the uniform distribution (Fig. 5). Similarly, there is a younger estimated age for tracheophytes with the skew-t and Cauchy models compared with the uniform distribution (Fig. 5). The age of the tracheophyte node ranges from 472.2 Ma to 422.4 Ma across all alternative dating strategies.

Discussion

Our results demonstrate that divergence time analyses of early land plant evolution are largely insensitive to tree topology and dataset size; however, they show some sensitivity to calibration strategy and, in particular, the calibration on crown Embryophyta. This clearly demonstrates the informative nature of the calibration on crown Embryophyta, which is comparatively narrow in its temporal range (515.5 Ma to 469.0 Ma). The soft maximum constraint on the age of this clade is based on the maximum age of the oldest-possible nonmarine palynomorphs, encompassing all possible total-group embryophyte records (). Land plant spores are encountered commonly among marine palynomorph assemblages, and they have the same fossilization and sampling potential as acritarchs. However, the oldest-possible embryophyte records are preceded stratigraphically by thick sequences bearing only marine palynomorphs. These marine palynomorphs demonstrate that the conditions required for preserving embryophyte remains obtained and, thus, the absence of land plant spores constitutes evidence that embryophytes were not present at this time (42). Thus, we discount the results of the divergence time analyses in which the embryophyte calibration is not employed. Similarly, the skew-t and Cauchy distributions, which reflect a nonuniform probability of divergence timing between the minimum and maximum constraints, suggest younger clade ages. However, these nonuniform distributions are unduly informative, since we have no insight or additional evidence that might inform the probability of the time of divergence between minimum and maximum constraints. Hence, we reject the ensuing results in favor of those based on a uniform distribution which reflects equal probability of divergence timing between minimum and maximum constraints. Since the remaining sources of uncertainty have little impact, a holistic timescale encompassing all relevant uncertainties is, effectively, that represented in Fig. 2. It is difficult to foresee how higher precision can be achieved while also maintaining accuracy. We have shown that additional sequence data and topological uncertainty have little material impact, both perhaps as a consequence of the short temporal succession of clade divergences among early embryophytes and attendant issues such as incomplete lineage sorting. Improved taxon sampling among liverworts and hornworts (especially) is likely to yield more precise estimates for divergences among bryophytes on some topologies, as would improved sampling of their fossil record—which our analyses predict to extend deep into the Lower Paleozoic.

It is possible that a Total Evidence approach (43), integrating living and fossil species, both morphological and molecular data and evolutionary models, will leverage some increased precision. Perhaps more importantly, such an approach might provide a means of more precisely dating the origin of land plant body plan innovations (e.g., stomata, leaves, rooting systems) that have been considered influential in the evolution of the Earth System (44). In the interim, our evolutionary timescale achieves precision while also integrating all of the principal sources of uncertainty, providing a framework for inferring plant evolutionary history, the veracity of its fossil record, and the impact of phytoterrestrialization on the evolution of global biogeochemical cycles.

The Origin of the Embryophytes and Tracheophytes.

Considering the 95% HPDs of divergence times across all topologies, the origin of crown embryophytes is dated to 515.1 Ma to 470.0 Ma (middle Cambrian–Early Ordovician). However, all of the mean estimated ages are resolved within the Phanerozoic across all alternative topologies and dating strategies, and the majority are dated to around 500 Ma (middle Cambrian Series 2). Only one analysis has a 95% HPD that stretches into the Proterozoic. The full span of age estimates for the crown tracheophyte node is 472.2 Ma to 419.3 Ma (Floian, Early Ordovician to the late Silurian). Only one analysis has a 95% HPD that stretches to the Early Ordovician, with those using a uniform prior resulting in estimated mean ages close to the Ordovician−Silurian boundary (∼444 Ma). The span of the tracheophyte stem lineage ranges across all analyses from 25.1 My to 60.0 My; these intervals are shorter for the paraphyletic topology than the monophyletic bryophytes topology (35.5 My and 51.6 My, respectively) ().

Impacts of Alternative Topologies and Dating Strategies on Divergence Time Estimates.

The impact of analytical uncertainty on the estimated age of Embryophyta is minimized by the use of carefully selected temporal information from the fossil record. Differences in topology had a minimal impact on divergence time estimates for Embryophyta (Fig. 5 and Table 2). For each topology, the posterior age estimates conform largely to the specified calibration constraints on clade age (∼511 Ma to 469 Ma). Potential differences in age estimates for embryophytes only appear when the specified age constraint for this node is removed. On the hornworts–sister topology, age estimates for Embryophyta extend into the Proterozoic without the embryophyte calibration, whereas the monophyletic bryophytes topology yields congruent age estimates with or without the user-applied embryophyte age constraint (Fig. 5). Thus, topology can influence the estimated ages for nodes, but only when we ignore germane evidence from the fossil record. Therefore, the use of well-researched and justified fossil constraints, when incorporated alongside tests of model uncertainty, adds confidence in the conclusion of an Early Phanerozoic origin for embryophytes.

There are only minor differences across topologies for the estimated age of tracheophytes, as all trees produce comparable mean estimates (Table 2). One topology, hornworts−liverworts−mosses, produces a younger age from the 95% HPD interval (419 Ma) compared with all other trees (430 Ma), but this younger age is anomalous (i.e., slightly younger than the minimum derived directly from fossil evidence at 420.7 Ma) and has little overall support; the bulk of the posterior age of tracheophytes for the hornworts−liverworts−mosses tree is above 430 Ma.

Comparisons with the Fossil Record.

The first unequivocal embryophyte body fossil taxon, Cooksonia cf. pertoni, appears in the Wenlock [minimum age of 426.9 Ma (45)]. The first account of crown tracheophyte body fossils is shortly after, in the Ludlow [minimum age of 420.7 Ma (46)], followed by an apparent explosion of diversity in the Early Devonian (13). Our mean age estimates are older for both nodes, by 40 My for the embryophytes and 20 My for crown tracheophytes. However, in both cases, this is a consequence of a dearth of continental lithofacies before the late Silurian−Early Devonian (47). The earliest known fossils of embryophyte affinity are permanently fused tetrahedral tetrad cryptospores [sensu stricto Steemans (48), Wellman (49)] that have a long history of occurrences within marine deposits (13) from the Middle Ordovician [Dapingian; 469 Ma (50)]. Cryptospores of unclear affinity from the Cambrian [sensu stricto Strother (51)], while not considered unequivocally embryophyte, informed our soft maximum constraint (515.5 Ma). Our middle Cambrian−Early Ordovician estimate for the origin of crown embryophytes is compatible with an embryophyte interpretation; however, our results do not suggest that they reflect a protracted cryptic earlier evolutionary history. Likewise, the dispersed record of trilete spores that first appear in the Katian (Late Ordovician) (52), followed by an explosion of diversity in the Silurian (13), indicates an earlier origin for tracheophytes that is congruent with our estimates.

The main challenge in testing our divergence time estimates for the bryophyte lineages is their very poor representation in the rock record (13). Nevertheless, our results establish a predictive temporal framework for the stratigraphic intervals in which to prospect for fossils implied by the ghost lineages in our evolutionary timescale. Regardless of the topology, we date the first and second divergences within the bryophytes between 496.5 Ma and 456.2 Ma (late Cambrian–Late Ordovician) and 478.7 Ma and 438.0 Ma (Early Ordovician–early Silurian), respectively. The oldest credible candidate bryophyte fossil is the Pragian (Early Devonian) Riccardiothallus devonicus (53), although the security of its classification is limited by preservation of only gross morphology. The mismatch between the estimated ages and unequivocal fossil finds is contributed to by their low fossilization potential, principally because bryophytes do not biosynthesize lignin. When body fossils occur, they are often too poorly preserved to allow recognition of synapomorphies. However, some extant bryophytes produce permanent tetrads and dyads (54, 55) similar to the cryptospores. The wall ultrastructure of cryptospores, known from as early as the Middle Ordovician, is similar to the multilaminate walls observed in permanent tetrads produced by extant liverworts, such as Sphaerocarpos (56). The presence of liverwort-like spores in the Middle Ordovician is not incongruent with the estimated dates of divergence of the liverworts across all topologies in our analyses. Sporangia described from the Late Ordovician of Oman are significant fragments of plant anatomy recovered from very rare instances of nonmarine Ordovician rocks (57). The spore masses contain either dyads or tetrads, the former displaying multilaminate walls, and most specimens preserve at least a partial covering, making it very difficult to argue that they are anything but land plant sporangia (57, 58). Unfortunately, our understanding of the parent plants of cryptospores, the cryptophytes, is restricted to much later charcoal Lagerstätten in the Pridoli and Lochkovian (59, 60). These fossil plants possess a combination of both bryophytic and tracheophytic characters, and thus their taxonomic position is currently unclear (60). The confirmation of the main synapomorphy for the tracheophytes, the presence of vascular tissues, is particularly difficult to demonstrate, due to the minute size and fused nature of these fossils.

Theories on the process of terrestrialization have long argued for a close temporal relationship between the emergence of land plants and terrestrial animals, particularly arthropods, substantiated by their approximately concurrent first fossil occurrence in terrestrial facies (61, 62). However, this is likely an instance of pseudocongruence, with lineages of differing antiquity exhibiting coeval stratigraphic first occurrences because of secular variation in the preservation of Lower Paleozoic terrestrial facies (40). Thus, a shift from dominantly marine to terrestrial facies results in a telescoped first stratigraphic appearance of disparate terrestrial lineages (63). The results of our divergence time analyses indicate a much earlier (∼70 My to 80 My) origin of land plants, but, surprisingly, this remains congruent with the latest divergence time estimates for three or four independent transitions to terrestrialization among arthropod lineages (hexapods, arachnids, and, perhaps, twice among myriapods) (64). Thus, although our results corroborate the view that the early fossil records of terrestrial arthropods and land plants are temporally misleading, they also corroborate the hypothesis of a close temporal relationship between the emergence of land plants and terrestrial animals, with plants creating habitats suitable for terrestrial arthropods.

Comparisons with Previous Studies.

Previous analyses indicate either a Proterozoic (mainly Cryogenian) (65–67) or Phanerozoic (68–70) origin of the embryophytes. Of the latter, dates range from the Early Ordovician [∼474 Ma to 477 Ma (69, 70)] to early Silurian [435 Ma to 425 Ma (68)]. The majority of our results are congruent with a Phanerozoic origin, but with older estimated ages (middle Cambrian; Fig. 5), reflecting the use of Cambrian cryptospores as a soft maximum constraint on crown embryophyte divergence. In comparison with the fossil record, a Phanerozoic origin of the embryophytes is more tenable than the Proterozoic, which is effectively precluded by the absence of embryophyte remains in marine sequences that nevertheless preserve sporopollenin acritarchs (42).

The origin of the crown tracheophytes has been fixed as a calibration point in most previous studies. Estimated ages include the Late Ordovician (446 Ma) (67), mid-Silurian (432 Ma to 434 Ma) (69), and late Silurian (423.95 Ma) (70). Our analyses are most congruent with the older ages estimated by Clarke et al. (67), around the Ordovician–Silurian boundary (Fig. 5), as a result of the application of an older taxon for the calibration [e.g., Zosterophyllum instead of Leclercqia as in Smith et al. (69)], and a soft maximum age constraint using the first occurrence of trilete spores in the Katian.

Few molecular clock studies focus on bryophyte divergence and, as such, often have restricted analyses to stomatophytes (mosses vs. vascular plants) (65, 66), including very few taxa. Estimates for the first bryophyte divergence begin as early as the Cryogenian (65, 66), with further studies suggesting the Ediacaran to late Cambrian (632 Ma to 499 Ma) (67), late Cambrian to late Silurian (490 Ma to 425 Ma) (68), Late Ordovician (458 Ma) (70), and mid-Devonian (383 Ma) (69). Our age estimates are most congruent with an Early Paleozoic divergence. Where previous studies have included all bryophyte lineages, the second divergence has been estimated from the Early Cambrian–Middle Ordovician (67) (524 Ma to 460 Ma), early Silurian (70) (440 Ma), and Mississippian (69) (335 Ma). Our age estimates are more congruent with the older estimates.

Implications for Hypotheses on the Coevolution of Land Plants and Climate.

The evolution and geographical spread of the embryophytes across Paleozoic continents undoubtedly had a major impact upon global biogeochemical cycles. To test hypotheses on the coevolution of land plants and Earth’s System, biogeochemical models rely on a well-substantiated phylogeny and timeline of embryophyte divergence and character acquisition. The GEOCARB (3, 4) and COPSE (5) biogeochemical models include parameters for the evolution and geographical spread of tracheophytes and their enhancement of silicate weathering rates, resulting in simulations that show a significant decrease in atmospheric CO2 levels in the Devonian [from ∼16× to ∼3× present atmospheric level (PAL)] and the rise in O2 levels to 1.5× PAL by the end of the Carboniferous. However, these models are undermined by their use of the body fossil record to establish a timescale for plant evolution and innovations. These weaknesses can be overcome by considering the divergence time estimates of key innovations from molecular clock studies.

Our results demonstrate that embryophytes were present on land from the middle Cambrian−Early Ordovician interval, and minimally, by the early Silurian, the four major lineages of land plants had already diverged and were constituents of early cryptogamic ground covers (71). Plants had already evolved key adaptations for survival and proliferation on dry land by the early Silurian (e.g., development of an embryo, alternation of generations, aerial sporophytes, sporophyte branching, cuticle, stomata, vascular tissue, sporopollenin-coated spores), including interactions with early soils and nutrient extraction from minerals (rhizoids, rhizomes, and symbiosis with mycorrhizal fungal partners). The results of our analyses suggest that the majority of these characters had evolved within a middle Cambrian−Early Ordovician interval. Modern cryptogamic covers, that comprise bryophytes, lichens, fungi, algae, and cyanobacteria, are capable of significant mineral weathering (72, 73), in particular via symbiotic mycorrhizal fungal partners accessing phosphorous (7), a limiting nutrient, which results in a positive feedback mechanism with increasing biomass of the host plant. As such, the timing of divergence and weathering capabilities of these early ground covers has been underestimated in these biogeochemical models.

Conclusions

The origin and evolution of land plants has transformed the terrestrial biosphere. Our understanding of the timing and nature of this formative episode is undermined by uncertainties associated with the incompleteness of the plant fossil record and the evolutionary relationships of the living land plant lineages. We establish an evolutionary timescale that integrates over these uncertainties, estimating the living clade of land plants to have emerged in the middle Cambrian−Early Ordovician, and the living clade of vascular plants to have appeared in a Late Ordovician−Silurian interval. These are in close accord with estimates for the timing of terrestrialization of arthropod lineages. These results underscore the importance of taking an integrative approach to the establishment of evolutionary timescales, which can only be derived through application of molecular clock methodology (74). Future attempts to explore the role of plant phylogeny in the evolution of global biogeochemical cycles must integrate this recalibrated timescale for plant evolution, rather than relying on the fossil record alone.

Methods

Dating Analyses.

We conducted all dating analyses in MCMCTree within the software PAML version 4.8 (75), and all analyses were prepared using MCMCTreeR in R (https://github.com/PuttickMacroevolution/MCMCTreeR).

Genetic Data.

We used two datasets from the published nucleotide alignments of Wickett et al. (14) for all analyses. For the first dataset, we used the full 852-gene alignment of 1,701,170 nucleotides. We used a subset of these data that were filtered by Wickett et al. (14) to maximize coverage of sites and genes, remove potential contamination, and exclude the third codon position. These data consist of 290,718 nucleotides. Unless specified, all subsequent analyses were conducted using the dataset of 1,701,170 nucleotides.

We estimated topology using topological constraints to enforce each of the seven hypotheses (Fig. 1) but leaving all other relationships unconstrained, using the focal dataset of Wickett et al. (14) (290,718 nucleotides, trimmed and third codon removed). For each hypothesis, we constrained tracheophytes, each bryophyte group (liverworts, hornworts, mosses), and the non-embryophytes. With each of these constraints, we left all other relationships as polytomies. We estimated these topologies in RAxML 8.2 (76) in a nonpartitioned, nucleotide GTR + Γ model.

We explored the impacts of dataset size (number of nucleotides) and site completeness. Plots of infinite sites were used to gauge the potential increase in precision gained by adding more sequence data. We compared infinite site plots of the original sequence data (852 genes, 1.7 million nucleotides) to data we trimmed by site completeness so that only sites complete for 50%, 75%, 99%, and 99.9% of species were included; this produced datasets with 850,000, 435,000, 19,000, and 2,000 nucleotides, respectively. These initial analyses indicated that there is not much effect in adding more sequence data (Fig. 4), and thus, for comparisons of all seven hypotheses, we employed the dataset trimmed by 50% completeness (850,000 nucleotides).

Rate Priors.

To incorporate deviations of a strict molecular clock, we set the IGR model that treats branch rates as being samples from independent and identically distributed log-normal distributions (77, 78). This distribution is given a prior mean rate for branches (μ), and variance σ2 that models the overall rate variability on branches across the phylogeny. In MCMCTree, the mean rate is given a prior gamma distribution with user-specified shape and scale. To obtain a suitable prior on the substitution rate (μ), we compared the pairwise distance between Arabidopsis thaliana and Rhynchostegium serrulatum using the GTR + Γ + F model in baseml version 4.8 (75). For the smaller dataset, this resulted in a substitution rate of 0.08−10 changes per nucleotide site per year after assuming a divergence time of 469 Ma. In the larger dataset, this value was 0.09−10 nucleotide substitutions per site per year. As in dos Reis et al. (79), we fixed the shape parameter of the gamma distribution prior on rate to 2, and, from this, set the scale parameter to 25. For the larger dataset, these figures were set to shape 2 and scale 22. We estimated these parameters for each of the subsets of the larger dataset. We set the prior on rate variability (σ2) as a gamma distribution with shape 1 and scale 10.

Time Priors.

For the priors on branching times, we set the prior birth−death process with parameters of birth = 1, death = 1, and fraction of sampled species = 0, which produce a uniform kernel for the branching times. The time prior or the prior for divergence times for all nodes in the tree is generated in conjunction with the specified node age densities based on the fossil record. The specified calibration densities and the effective time prior can be very different (41). To ensure our priors on divergence times were appropriate, we ran the model without sequence data to obtain the effective priors.

Fossil Ages and Prior Node Distributions.

In each analysis (unless stated), we applied temporal node constraints to 37 nodes, including the root. The location of the 37 nodes is shown in . We applied node distributions using minimum and maximum constraints following protocols outlined in Parham et al. (37). For full phylogenetic and age justifications of each fossil calibration, see .

Three strategies were applied to specifying the prior distributions on node ages. In strategy i, uniform distributions were applied to all internal prior node ages with a hard minimum age and a soft maximum age, allowing 0.001% probability of an age younger or older than the given minima and maxima. For strategy ii, we applied skew-normal distributions with the mode of the distribution above the minimum age and 0.001% and 97.5% probability tails at the maximum and minimum ages. For strategy iii, we applied Cauchy distributions with a hard minimum and a 97.5% probability at the maximum age. For strategy iii, the root node was set as a uniform distribution. For each strategy, we assessed the shape of prior and posterior distributions on the 37 nodes to which we applied data from the fossil record (shown for the hornworts−sister topology, ). The specific parameters used for input into MCMCTree are shown in .

Analyses of large datasets can be highly time-consuming. Therefore, we implemented the approximate likelihood calculations available in MCMCTree (80, 81). We obtained estimates of branch lengths in baseml (82), and, in the program, these maximum likelihood estimates are then used to obtain the gradient and Hessian matrix of the branch lengths. These estimates were then used to calculate the approximate likelihood (81) in the divergence time analyses.

Two of the key nodes we were primarily interested in dating were crown embryophytes and crown tracheophytes. We conducted several sensitivity analyses to explore any potential variation in the age estimates. We tested the effect of removing the nonembryophyte (algal) species from the analysis so the embryophyte node became the root node. In a separate analysis, we removed the user-applied node constraint for embryophytes. We also explored the impact of applying a correlated clock model to the data (80). Additionally, we explored the effect of using topologies based on the maximum likelihood tree search of the 1.7-million nucleotide dataset; the largest difference in this topology is that Chara vulgaris is sister to embryophytes rather than Zygnematophyceae in the main analyses. Finally, we explored the effects of codon partition by comparing the posterior age estimates of a single partition (all codons in a single alignment) and a partition of each codon (three alignments for positions 1, 2, and 3). These analyses indicated that partition did not have any meaningful influence on posterior age estimates for all nodes ().