The angiosperm radiation has been linked to sharp declines in gymnosperm diversity and the virtual elimination of conifers from the tropics. The conifer family Podocarpaceae stands as an exception with highest species diversity in wet equatorial forests. It has been hypothesized that efficient light harvesting by the highly flattened leaves of several podocarp genera facilitates persistence with canopy-forming angiosperms, and the angiosperm ecological radiation may have preferentially favoured the diversification of these lineages. To test these ideas, we develop a molecular phylogeny for Podocarpaceae using Bayesian-relaxed clock methods incorporating fossil time constraints. We find several independent origins of flattened foliage types, and that these lineages have diversified predominantly through the Cenozoic and therefore among canopy-forming angiosperms. The onset of sustained foliage flattening podocarp diversification is coincident with a declining diversification rate of scale/needle-leaved lineages and also with ecological and climatic transformations linked to angiosperm foliar evolution. We demonstrate that climatic range evolution is contingent on the underlying state for leaf morphology. Taken together, our findings imply that as angiosperms came to dominate most terrestrial ecosystems, competitive interactions at the foliar level have profoundly shaped podocarp geography and as a consequence, rates of lineage diversification.

1. Introduction

The angiosperm radiation has been well documented although key elements are poorly understood [14]. The concurrent, sharp decline in diversity of several non-flowering land plant groups is consistent with competitive exclusion by angiosperms, although based on available data more complex hypotheses (e.g. climatic perturbations and biogeographic effects) cannot be ruled out [2]. Furthermore, the mechanism(s) facilitating the rise of the angiosperms remain uncertain, in part because ‘key innovations’ appear to have arisen repeatedly through the timescale of their evolution [3,4]. In contrast, gymnosperms appear to have been less adventurous in the exploration of morphological and ecological space [4] and as a consequence, can provide outstanding case studies linking functional morphology, ecology and their evolutionary consequences in the context of long-term environmental change [5,6].

Conifers were a diverse and important component of the vegetation until the Late Mesozoic but in the present tend to dominate only where the physical environment is limiting (such as on colder sites and on oligotrophic soils). The drivers of conifer displacement by angiosperms are still debated [711]. An influential hypothesis [5,8] that has been widely accepted in the literature [9] invokes vegetative growth rate differences and, broadly, predicts that the relatively inefficient stem and leaf hydraulics of conifers lead to lower maximum growth rates, placing them at a competitive disadvantage in the ‘regeneration niche’ relative to angiosperms on productive sites. While Bond's hypothesis [5] seems broadly applicable to Northern Hemisphere conifers (e.g. Pinaceae), Podocarpaceae, a predominantly Southern Hemisphere family, is enigmatic in attaining its highest species richness within angiosperm-dominated equatorial forests [6] (see the electronic supplementary material, table S1) and is arguably the most successful conifer family within these settings. Understanding the origins of this exceptional pattern could shed light, more generally, on the nature and evolutionary dynamics of conifer–angiosperm competitive interactions and how these have shaped contemporary broad phytogeographic patterns.

A remarkable aspect of the Podocarpaceae is the diversity of leaf morphologies ranging from scale-like and needle-like to functionally broad (see the electronic supplementary material, table S1). Several distinct podocarp lineages develop bilaterally or bifacially flattened leaves, cladodes and short shoots comprising distichously arranged broad leaves analogous to the compound leaves of angiosperms [6,12] (figure 1; electronic supplementary material, table S1). Typically, conifers have a single vein leaf anatomy that hydraulically limits maximum leaf width [13] and therefore the foliar efficiency of light harvesting [14]. In contrast, angiosperms have evolved leaves with a high-density network of veins that efficiently supply water to the mesophyll independent of leaf width [13,15]. A number of podocarp lineages appear to have partially overcome the conifer hydraulic limitation, developing flattened photosynthetic surfaces analogous to angiosperm leaves that possess either multiple veins, or an extravenous hydraulic network comprising sclerified tissues in the mesophyll that facilitates radial transport of solutes from the midrib to the leaf margins [15]. Leaf flattening in conifers is associated with high photosynthetic efficiency under low light conditions [12,13,16,17] and a plausible hypothesis links the relative success of Podocarpaceae in productive, angiosperm-dominated environments with a high degree of shade tolerance as a consequence of leaf and shoot flattening [6,12]. Although the precise timing of modern rainforest development remains uncertain, new data show a surge in angiosperm leaf vein density in fossil deposits and phylogenetic reconstructions, providing strong evidence that the climate [18] and productivity of angiosperm-dominated rainforest arose in the period between the Late Cretaceous and the Palaeogene [19]. As angiosperms came to dominate productive terrestrial environments, changing light environments may have favoured the evolution of shade-tolerant leaf and shoot architecture among podocarps [6,12] or alternatively, forced them into marginal habitats where angiosperm competition is reduced [5]. One difficulty in testing these ideas is that the fossil record of Podocarpaceae is relatively poorly known for the Late Cretaceous [6].

Figure 1.

Examples of the diversity of shoot flattening in nine genera of Podocarpaceae: (a) Retrophyllum, (b) Dacrycarpus, (c) Falcatifolium, (d) Acmopyle, (e) Podocarpus, (f) Nageia, (g) Prumnopitys, (h) Phyllocladus and (i) Sundacarpus. As well as obvious planation of short shoots (ad,g,i) podocarps have evolved (f) multi-vein leaves and (h) phylloclades, and internal conducting sclereids (see text). See the electronic supplementary material, table S1 for further details.

Phylogenetic reconstructions of evolutionary relationships have proved valuable in testing hypotheses of vegetation change where the fossil record is limited or ambiguous [20,21]. If the spread of modern rainforests has favoured shade-tolerant podocarps, phylogenetic patterns may indicate shifts in evolutionary tempo among Podocarpaceae that can be related to leaf morphology and are coincident with the timing of the rise of angiosperm ecological dominance. Here, we develop a molecular phylogeny for Podocarpaceae using fossil constraints to estimate the timing of lineage divergences within the family and from this perspective we consider the influence of leaf evolution on the distribution and diversity of podocarp lineages.

2. Methods

(a) Molecular data

We assembled a dataset comprising DNA sequences for three chloroplast fragments (partial matK, rbcL and trnL-trnF intron and spacer, ca 3900 bp) for 118 taxa, including 115 Podocarpaceae, representing 19 genera and an Araucariaceae outgroup (see the electronic supplementary material, table S2). For sequences generated de novo, PCR and sequencing followed standard protocols, and where possible, data for each species were obtained from the same accession. Sequence alignment was performed manually (rbcL and matK) or using CLUSTALX [22] followed by manual adjustment (trnL-trnF intron and spacer). For each of the target regions, PCR and sequencing primers are listed in electronic supplementary material, table S3.

(b) Phylogeny and divergence time estimates

The Bayesian molecular clock implementation BEAST v. 1.5.2 [23] was used to estimate the phylogeny of Podocarpaceae under an uncorrelated lognormal-relaxed clock model of rate variation among lineages and a GTR+I+Γ model of sequence evolution that was unlinked across data partitions. We incorporated five fossil constraints (representing the oldest reliably placed macrofossil associated with each constrained node, electronic supplementary material, table S4) to provide absolute estimates of lineage divergence times for Podocarpaceae. For four of the constraints, we assumed a lognormal prior distribution of node ages with the fossil age used to approximate the minimum age of the associated stem group. Stem group placement provides an objective criterion when a fossils' membership of the crown group is uncertain [24] but we used a broad confidence interval on the prior (i.e. much older node ages are not rejected a priori) in the event that the fossil belongs within the crown. In addition to minimum age, parameters of the lognormal prior are the mean (peak probability) and the standard deviation. We experimented with three different lognormal prior means: 2.5, 3.0 and 3.5, resulting in median prior node age estimates that are 20 (mean = 2.5)—40 per cent (mean = 3.5) older than the oldest reliably placed fossil for that lineage. We used a standard deviation of 1.0 yielding a prior distribution that includes ages up to three times older than the minimum age on the node (see the electronic supplementary material, table S4). For the root constraint, we used a normal prior distribution with its mean approximating the fossil age, and a standard deviation of 15.0. While a Triassic origin for the family seems well-constrained (e.g. [6,25,26]), we avoided the use of a ‘hard’ lower bound on the prior given that the issues of correct fossil assignment are likely to increase with clade age (e.g. [27]). In addition to the above, we constrained a maximum root age of 320 millions of years before present (Ma) approximating the age of stem group conifers [26].

For each of our three constraint scenarios (varying the prior mean, see above), we ran two separate analyses over 1 × 107 generations sampling topology and model parameters every 1000 generations. We used Tracer v. 1.4 [28] to assess convergence between runs, and after excluding an appropriate burn-in fraction, the posterior estimates of topology and parameter values were obtained from the pooled sample of the six separate analyses.

(c) Character states and ancestral state reconstructions

We coded leaf morphology for the sampled podocarp species as either functionally broad (foliage flattening, FF), which includes lineages that have two-dimensionally flattened single-veined leaves; multiple-veined leaves; and modified branchlets (cladodes) that function as the primary photosynthetic organs (figure 1). The alternative state, imbricate leaves (IL), includes needle-like (as in Dacrydium and Dacrycarpus) and scale-like (as in Lagarostrobos and Microcachrys) morphologies [29]. Most Dacrycarpus species have bimorphic foliage (the exception being Dacrycarpus compactus, which develops needle-like leaves only) and can display either needle-like or significantly bilaterally flattened leaves on a single plant [30] and these were coded as polymorphic (see the electronic supplementary material, table S1).

We used a reverse jump Markov chain Monte Carlo (RJ-MCMC) approach implemented in BayesTraits v. 1.0 [31,32] to reconstruct leaf morphology on our phylogenies. As input to these analyses, a sample of 1000 (post burn-in) topologies were drawn from our pooled phylogenetic analyses. The ‘ratedev’ parameter was manually set so that the acceptance rate ranged between 20 and 40 per cent. Priors of the rate parameters were estimated using a hyperprior approach [32] with an exponential distribution, its mean seeded from a uniform distribution on the interval of 0–10. Five separate analyses were run over 4 × 107 generations, and probability values greater than or equal to 0.95 for a given state at that node were taken to indicate positive support for that state. When nodal probabilities fell below 0.95, the Bayes factor (BF) was used to assess the strength of support for alternative states, with a BF of two to five indicating positive, and a BF of greater than five indicating strong support [33] for a given state at that node. BF estimates were based upon five separate runs of 4 × 107 generations for each character state, with analysis settings as outlined above. Outputs were visually inspected to ensure stability of BF estimates.

We were interested in correlations between leaf evolution and environmental productivity, as indexed by bioclimatic parameters. For climatic estimates, georeferenced locality data for 87 podocarp species were sourced from the Global Biodiversity Information Facility (GBIF; accessed through GBIF Data Portal, data.gbif.org, 2010-10). The selected taxa were included in our phylogeny and in most cases are represented by 10 or more GBIF data points, or otherwise have very localized distributions (e.g. Dacrycarpus kinabaluensis). Diva-GIS [34] was used to derive values for bioclimatic indices based upon temperature and precipitation from the WorldClim [35] interpolated climate surfaces (see the electronic supplementary material, table S5). Following the thermal classification of Nix [36], we coded a binary character including a microthermal (temperate; mean annual temperature, MAT < 12°C) and meso-megathermal (paratropical/tropical; mesothermal MAT 14–20°C, meso-megathermal interzone MAT 20–24°C; megathermal MAT > 24°C) climatic response. We also recognized a microthermal–mesothermal interzone (MAT 12–14°C) [36], including predominantly upper montane tropical species, which were coded as polymorphic (0/1).These grouping were supported by a cluster analysis (Gower Metric) based upon three climatic variables (MAT, coolest month mean temperature and mean annual precipitation; the electronic supplementary material, figure S1) that summarize environmental energy and its annual variation, and annual moisture regime. The mycoheterotroph Parasitaxus usta was coded as missing data for climatic range. For taxa included in our phylogeny, but not in the climatic estimates, character coding was extrapolated from related taxa with similar/overlapping distributions (see the electronic supplementary material, table S5).

To test for contingency among climatic range and leaf morphological traits, we compared the fit of our data with an independent versus a dependent model of trait evolution [32]. The independent model assumes that the state of one character has no effect on the rate of evolution of the other while for the dependent evolutionary model transitions between states of one trait are influenced by the underlying state of the other. These models were compared using Bayestraits with RJ-MCMC optimization and settings as outlined above. Support for the independent versus dependent models was assessed using a BF.

For all trait optimizations, Araucariaceae were excluded to avoid conditioning our reconstructions on a sparsely sampled, evolutionarily distant outgroup [37].

(d) Diversification rates and diversification rate shifts

Following Magallón & Sanderson [38], we derived a 95 per cent CI on the expected number of species included within a hypothetical stem group for each interval of time from its origin onwards (time = 0–170 Ma, assuming a relative extinction fraction of 0.95). Using the mean age estimate for the stem group, standing diversities for the genera were compared with these sets of critical values, and those exceeding the upper or lower critical values were considered unexpectedly species-rich or poor, respectively, given the estimated overall diversification rate for the family. Calculations were performed using the R package [39] GEIGER [40]. We also derived a confidence interval (as above) on the expected number of species through time based upon a rate estimate for the node where a significant diversification rate shift was detected (see below).

We used LASER [41] to test if the phylogeny of Podocarpaceae, firstly, could have been generated under a homogeneous rate of diversification across lineages, by contrasting the likelihoods of a model where all lineages have diversified at a constant rate with a variable rate diversification model. Secondly, in rejecting a time-homogeneous rate, we wished to identify the branch where the largest diversification rate shift has occurred. Under the variable rate diversification model, the maximum likelihood (ML) diversification rate shift point is the node with the highest combined likelihood obtained by summing the lineage-specific likelihood estimates from each bipartite sub-tree [42]. For these analyses, we used a sample of 300 topologies from the phylogenetic analyses using BEAST. These were pruned to represent generic-level sampling, with the diversities of the terminals (genera) assigned according to Farjon [43] (see the electronic supplementary material, table S1). Analyses in LASER were performed under a relative extinction rate of 0 (no extinction) and 0.95 (95% of lineages go extinct).

3. Results and discussion

(a) Phylogeny and divergence times

Our phylogenetic analyses provide a well resolved and statistically well-supported hypothesis of evolutionary relationships among Podocarpaceae (figure 2; the electronic supplementary material, figure S2). As with previous molecular phylogenetic studies, the monophyly of conventionally recognized genera is strongly supported with the exception of Prumnopitys, which includes Sundacarpus (Prumnopitys s.l.) [4446]. Most intergeneric relationships are also unambiguously resolved, including a ‘prumnopityoid’ clade (node 4) and a ‘podocarpoid’ clade (node 11), which is sister to a ‘dacrydioid’ clade (node 12; podocarpoid–dacrydioid clade; node 10). The position of Phyllocladus has not been well resolved in previous analyses, although here we find a strong support for its placement within Podocarpaceae (node 2) rather than its treatment as a distinct family (e.g. [47]) (figure 2).

Figure 2.

Phylogeny of Podocarpaceae (generic-level sampling, refer the electronic supplementary material, figure S2 for full sample) derived using Bayesian-relaxed clock analysis of plastid DNA sequences. Branches are proportional to time, and horizontal grey bars represent the 95% HPD of divergence times for the associated node. All branches have a PP = 1.0, unless indicated adjacent to the branch. Pie charts represent posterior probabilities of reconstructed ancestral states for leaf morphology using Bayesian optimization criteria (dark grey, IL; light grey, FF). Numbers adjacent to nodes and highlighted taxon groups are referred to in the text and/or electronic supplementary material, table S6. Most Dacrycarpus (asterisked) are polymorphic for leaf morphology.

We estimate a Triassic–Jurassic origin (205 (177–237) Ma) for the split between Podocarpaceae and Araucariaceae, consistent with the early fossil record of both of these families [6]. The Podocarpaceae crown group (node 1) began to diversify from the Mid-Jurassic to Mid-Cretaceous (145 (99–194) Ma), and most of the extant genera by the Early Cretaceous (e.g. Saxegothaea, Microcachrys, Pherosphaera and Phyllocladus) to the Palaeogene (e.g. Afrocarpus, Dacrydium, Falcatifolium and Nageia). By and large, our age estimates indicate Cenozoic radiations for the crown groups of most podocarp genera (figure 2; the electronic supplementary material, figure S2).

(b) Leaf evolution and the angiosperm ecological radiation

Our reconstructions of leaf morphology (figure 2; the electronic supplementary material, table S6) indicate that extant Podocarpaceae were ancestrally IL, and therefore, that FF types are derived within the family. FF has evolved independently at least five, and probably six, times among extant podocarps: separately within each of Phyllocladus, Prumnopitys s.l., Saxegothaea, Acmopyle, Falcatifolium and is the most probable state for the most recent common ancestor (MRCA) of Podocarpoid clade (node 11, including Podocarpus, Afrocarpus, Nageia and Retrophyllum). FF taxa with putative podocarp affinities are also known from the Triassic and Jurassic, although evidently these have left no close relatives (e.g. [25]). Flexibility in leaf morphology appears to be an important evolutionary response of Podocarpaceae to changing environmental conditions (temperature, rainfall and light climates) [6,48] throughout the history of the family.

As well as being derived, the evolution of FF appears to be a relatively recent phenomenon among podocarps, arising predominantly from the Late Cretaceous through to the Palaeogene (figures 2 and 3). Our phylogenetic reconstructions are broadly consistent with the Australasian macrofossil record of Podocarpaceae where there are several first appearances of extant FF lineages in the Palaeogene including Acmopyle (Late Palaeocene), Falcatifolium (Middle Eocene), Podocarpus (Early Eocene), Phyllocladus (Middle Eocene) Prumnopitys (Early Palaeocene) and Retrophyllum (Middle Eocene), along with several globally extinct FF genera including Smithtonia and Willungia [6]. These fossil dates have been taken to suggest response to changing light environments associated with the expansion of Modern rainforests in Southern Australia during the Early Cenozoic [6,12]. FF conifers have light response characteristics well suited to growth in the shade [12,17] and planated shoot architecture that is associated with efficient foraging for light [16]. A shade-tolerant strategy appears to be an important aspect of conifer coexistence in the understorey and canopy alongside angiosperms [17,49].

Figure 3.

Cumulative number of podocarp lineages (log scale) through time partitioned by leaf morphology (a) FF; (b) IL based upon ancestral state reconstructions for Podocarpaceae (most probably state at each node). In each instance, the solid line represents the mean node height and the dashed lines, the upper and lower limits of the 95% HPD of divergence times for the maximum clade credibility tree inferred from the Bayesian-relaxed clock analyses. The grey bar represents the approximate timing of the rise of fossil angiosperm leaf vein densities to the levels comparable with Modern tropical rainforest species (modified from [19]).

An association between leaf traits and lineage accumulation rates is supported by our analyses (figure. 3). Among FF groups, lineage accumulation through time is exponential (log-linear), commencing from the Late Cretaceous to the Palaeogene (figure 3a). This pattern is consistent with a temporally homogeneous birth–death rate, mediated by a relatively low rate of extinction [38]. In contrast, the IL lineages show an initially high rate of lineage accumulation that subsequently plateaus, and show a further recent upturn (ca 20 Ma; figure 3b). While a biological interpretation is not straightforward, a preponderance of ‘early’ IL lineages is consistent, for example, with a declining diversification rate mediated by temporally declining speciation [50]. Alternatively, an anti-sigmoid curve (figure 3b) could suggest a mass extinction event, which can generate a punctuated phylogenetic pattern [51]. We note a reasonable coincidence with the timing of angiosperm vein density expansion [19] and the radiation of FF podocarps (figures 2 and 3). Vein density evolution in angiosperms has been linked not only to the transition of angiosperms from understorey to canopy [19], but also to the formation of rainforest at low latitudes owing to high transpiration-fed recycling of moisture [18]. Both changing light climates associated with the formation of angiosperm-dominated canopies, and a rapid transition of angiosperm leaf vein densities that dramatically increased humidity [18], could have preferentially favoured leaf flattening among Podocarpaceae. Conifer scale and needle leaves are typically associated with high light environments [12,13] and conceivably, declining diversification rates among IL lineages could reflect competitive interactions with angiosperms.

(c) Biogeography of leaf morphology

We found that the overall podocarp radiation is characterized by a significantly heterogeneous diversification rate (LRT; p < 0.012, based upon 300 sampled topologies and two relative extinction rates). We consistently detected a rate increase (mean likelihood difference between maximum rate shift and alternative nodes, 3.9 ± 1.4 SD) associated with the MRCA of the podocarpoid and dacrydioid clade (node 10) commencing ca 89-58 Ma (figure 2). We used this ‘new’ rate estimate (0.029 spp per million years, assuming a relative extinction rate of 0.95) to derive a 95% confidence envelope around the expected diversity of clades through time (the electronic supplementary material, figure S3). The majority of FF clades fall within the bounds of this envelope and the hypothesis of an elevated diversification rate relative to Podocarpaceae overall cannot be rejected in these instances. However, the same is true of Dacrydium, of which all species develop needle-like leaves, and Dacrycarpus, which is commonly imbricate leaved. Interestingly, the higher diversification rate can be applied to virtually all ‘tropical’ clades (see the electronic supplementary material, table S4 and figure S3) while genera with ‘temperate’ geographical distributions tend to be species poor, suggesting a role for geography in explaining clade-specific diversities.

We contrasted models of trait evolution and found that dependency among leaf morphology and bioclimatic range was favoured over a model where these traits evolve in an independent manner (BF ≈ 5.6). For the independent model, paired rates are expected to be equal, while the distribution of parameter rate estimates from our data strongly suggests this is not the case (figure 4; the electronic supplementary material, figure S4). For the dependent evolutionary model, the posterior probability of a non-zero rate for each rate category can be estimated by the proportion of sampled models that assigned a value of zero to that parameter. Both meso-megathermal transitions among FF lineages, and microthermal transitions among IL lineages, have non-zero rate estimates (posterior probability, PP = 0 and ≈0.01, respectively) while for the remaining parameter estimates, a zero rate cannot be rejected (PP ranging from ca 0.27 to 0.9; the electronic supplementary material, figure S4). These analyses, therefore, imply that the direction of climatic range evolution is contingent upon the underlying state for leaf morphology (figure 4).

Figure 4.

Flow diagram showing the eight possible transitions (qij) between binary traits for leaf morphology and climatic range (1, imbricate leaves, microthermal; 2, imbricate leaves, meso-megathermal; 3, flattened foliage, microthermal; 4, flattened foliage, meso-megathermal). The direction of the arrow indicates the direction of change in a trait, and arrow thickness is commensurate to the mean rate estimate from the posterior sample. Transitions for which a zero rate is rejected are indicated by the solid arrows, while for dashed arrows, a non-zero rate has a posterior probability of less than 0.75 (refer to the electronic supplementary material, figure S4).

Geographical range size and clade diversities are linked and, in particular, both speciation and survivorship rates are positively related to the geographical extent of a taxon [52,53]. In the context of our findings, range effects on speciation and extinction rates could explain contrasting diversities among podocarp lineages (see the electronic supplementary material, figure S3). Podocarps, in general, are strongly limited by drought [6,54] and in the Present, the majority of IL lineages is associated with high rainfall upper montane and alpine conditions, or they grow in high-latitude temperate rainforests as long-lived pioneers or riparian species [55]. Such suitable environments are restricted in the Southern Hemisphere, as mesic temperate climates have contracted since the opening of the Southern Ocean from the Oligocene [56]. In Southern Australia, both increasing rainfall seasonality (a shift from summer to winter rainfall) and increasing fire frequency through the Neogene have been linked to sharp contractions and extinction of podocarp lineages [6,57]. The ability to coexist with productive canopy-forming angiosperms may have favoured the expansion of shade-tolerant (FF) podocarp lineages into geographically extensive and relatively stable humid equatorial climates while in contrast, IL lineages have been ‘cornered’, and clade growth has reached ecological limits. We propose that the angiosperm radiation may have primarily driven a shift in podocarp geography, mediated at the foliar level, consistent with Bond's hypothesis [5] for global phytogeographical patterns.


We thank J. Conran, G. Keppel, staff of the Royal Botanic Gardens, Edinburgh, the Botanic Gardens Trust, Sydney and South Australian Botanic Gardens, Adelaide and Mt. Lofty, for supplying tissues for DNA extraction and sequencing. P. Hollingsworth, S. Graham, P. Crane and five anonymous reviewers are thanked for comments on the manuscript. This research was supported under the Australian Research Council's Discovery Projects funding scheme (project number DP0665859 awarded to A.J.L.).

  • Received March 15, 2011.
  • Accepted May 18, 2011.


View Abstract