Life history influences rates of climatic niche evolution in flowering plants

Stephen A. Smith , Jeremy M. Beaulieu

Abstract

Across angiosperms, variable rates of molecular substitution are linked with life-history attributes associated with woody and herbaceous growth forms. As the number of generations per unit time is correlated with molecular substitution rates, it is expected that rates of phenotypic evolution would also be influenced by differences in generation times. Here, we make the first broad-scale comparison of growth-form-dependent rates of niche evolution. We examined the climatic niches of species on large time-calibrated phylogenies of five angiosperm clades and found that woody lineages have accumulated fewer changes per million years in climatic niche space than related herbaceous lineages. Also, climate space explored by woody lineages is consistently smaller than sister lineages composed mainly of herbaceous taxa. This pattern is probably linked to differences in the rate of climatic niche evolution. These results have implications for niche conservatism; in particular, the role of niche conservatism in the distribution of plant biodiversity. The consistent differences in the rate of climatic niche evolution also emphasize the need to incorporate models of phenotypic evolution that allow for rate heterogeneity when examining large datasets.

1. Introduction

Heterogeneity in the rate of molecular evolution has been widely documented (Laird et al. 1969; Kohne 1970; Sibley & Ahlquist 1984; Wu & Li 1985; Bousquet et al. 1992; Gaut et al. 1992, 1996, 1997; Andreasen & Baldwin 2001; Bromham et al. 2002; Smith & Donoghue 2008; Ho 2009) and has been attributed to a variety of causes, including differences in body size (Martin & Palumbi 1993; Bromham et al. 1996), metabolic rate (Martin & Palumbi 1993; Gillooly et al. 2005), DNA repair (Sarich & Wilson 1973; Britten 1986) and generation time (Mooers & Harvey 1994; Bromham et al. 1996), with correlations among these variables also being observed (Martin & Palumbi 1993). In plants, there is strong evidence for heterogeneity in the rate of molecular evolution between woody and herbaceous species (Gaut et al. 1992, 1996; Laroche et al. 1997; Kay et al. 2006; Smith & Donoghue 2008), and to a lesser extent between annuals and perennials (Andreasen & Baldwin 2001). Both the differences in rate between woody and herbaceous species, and between annual and perennial plants, probably reflect a difference in generation time (Sturtevant 1965; Wu & Li 1985; Li et al. 1987; Gaut et al. 1992, 1996; Smith & Donoghue 2008). This evidence suggests that plant lineages with longer generation times (i.e. the time to first flower) tend to accumulate nuclear substitutions at a slower rate per unit time than lineages with shorter generation times.

While the number of generations per unit time has been shown to be correlated with molecular substitution rates (e.g. Gaut et al. 1992, 1996; Kay et al. 2006; Smith & Donoghue 2008), it remains unclear whether life-history attributes such as generation time influence the tempo of phenotypic evolution. Early studies suggested that molecular and phenotypic rates were largely decoupled (Kimura 1983; Ohta 1992; Bromham et al. 2002). More recently, Davies & Savolainen (2006) found evidence for such a correlation in plants, with molecular change accounting for 2 to 11 per cent of the variance in the rate of morphological evolution (see also Xiang et al. 2008 for similar results within Cornus). Though evidence for a correlation between rates of molecular and phenotypic evolution is clearly mounting, to date there remains no clear consensus on this issue. A correlation between rates of molecular and phenotypic evolution would be remarkable and would provide evidence for the far-reaching consequences of life history across both genotypic and phenotypic scales.

Before molecular data were available, generation time was thought to exert a profound influence on the rates of phenotypic evolution (e.g. Sinnott 1916; Simpson 1944). The number of generations per unit time can influence the potential for adaptive phenotypic change by modulating the rate at which advantageous mutations arise in a population (Sinnott 1916). However, Simpson (1944, pp. 62–64) postulated that although generation time could cause dramatic differences in phenotypic evolution, this had not been widely observed. Yet phylogenetic information has not specifically been brought to bear on these questions relating generation time, morphological evolution and molecular evolution. Recent advances in large-scale phylogeny construction (Smith & Donoghue 2008; Smith et al. 2009) combined with sophisticated model-based phylogenetic comparative methods provide an ideal framework in which to adequately address this longstanding question.

Here, we demonstrate the correlation between growth form and rates of phenotypic evolution as measured by climate tolerance in flowering plants. We examine this pattern across previously published phylogenies representing a diverse set of flowering plant species. Because generation time is correlated with plant habit, designations such as ‘woody’ and ‘herbaceous’ provide convenient distinctions for testing this pattern of rate heterogeneity. In this study, we focus on climate tolerance, which includes a particularly interesting set of traits to examine. This is not only because abundant spatial data allow for the compilation of very large datasets, but also results pertaining to climate tolerance evolution have broad implications for the evolution of flowering plants and the ability of plants to respond to climate fluctuations.

2. Material and methods

(a) Phylogenetic trees

Smith & Donoghue (2008) previously examined molecular rate heterogeneity and plant habit in five large angiosperm phylogenies: Apiales, Commelinidae (Cantino et al. 2007), Dipsacales, ‘Primulales’ and Moraceae + Urticaceae (Rosales). From this study, we relied on both the time-calibrated phylogenies and the information on plant habit. We also use the uncalibrated branch lengths from these phylogenies to assess the adequacy of molecular branch lengths in properly scaling phenotypic change according to life history (see below).

(b) Geographical and climate data

To estimate a climate tolerance for a given species, we retrieved all GPS coordinates provided by the Global Biodiversity Information Facility (GBIF; http://gbif.org). This resulted in 889 species for Apiales, 3174 species for Commelinidae (commelinids; Cantino et al. 2007), 290 species for Dipsacales, 301 species for Primulales (Ericales, APG II) and 351 species for Moraceae + Urticaceae (Rosales). We then extracted the mean climate data for each of the 19 BIOCLIM variables at 2.5 arc-minutes resolution developed by Hijmans et al. (2005), which describe the major temperature and precipitation dimensions of a given species. The BIOCLIM variables are biologically meaningful layers derived from monthly rainfall and temperature values. Because our datasets involved taxonomically broad samples spanning a wide range in geographical and climatic spaces, we used all 19 BIOCLIM variables in our analyses.

(c) Testing for differences in climatic niche dimensions among growth form

To determine the climatic niche for each species in each clade, we conducted a principal component analysis across all 19 BIOCLIM variables. Because the data represent multicollinear species values, we used a principal component analysis that specifically incorporates phylogenetic information while calculating eigenvectors, eigenvalues, loadings and species scores (Revell in press). While this method takes into account phylogenetic information when estimating the eigenstructure of the data, the resulting PC scores remain in the original species space and therefore require a comparative approach for any subsequent analysis due to shared ancestry. The first two principal axes (PC1 and PC2) were chosen as these represented the major axes of temperature and precipitation dimensions, and accounted for a majority of the shared variation among all species of a given clade. These analyses were performed in R (R Development Core Team 2009) using code generously provided by L. J. Revell.

To assess whether there were mean differences in climatic niche space, we used a phylogenetic MANOVA (multivariate analysis of variance) to test whether significant differences in cross-species trait means among growth form were larger than expected based on a random model of Brownian motion (BM) evolution (sensu Garland et al. 1993). We used the R (R Development Core Team 2009) package GEIGER (Harmon et al. 2008) to generate 1000 Monte Carlo simulations using our input tree topology and time-calibrated branch lengths. All simulations were carried out under a gradual model of BM evolution and we relied on Wilks's λ as our multivariate test statistic. Wilks's λ measures the proportion of variance in multiple dependent variables (i.e. niche axes) that remains unexplained by the independent variable (i.e. growth form). The null distribution of Wilks's λ was compared against the observed Wilks's λ calculated from a conventional MANOVA. If the observed Wilks's λ was less than 95 per cent of the null distribution, then trait differences were greater than expected based on a model of BM evolution. We carried out this analysis separately for each clade.

To determine whether there were differences in the variance of climatic niche space between sister lineages that differ only in predominate growth form, we examined the observed variance in PC1 and PC2, separately. To do so, we calculated a sample variance (mean-squared error; MSE) that is the average squared difference between a tip value and the phylogenetic mean (â), or the maximum-likelihood ancestral trait value at the root node of a given tree (Schluter et al. 1997; Blomberg et al. 2003). In this case, the phylogenetic mean represented the ancestral trait value for a particular PC axis, analysed separately, at the root of a subtree consisting of sister lineages. To estimate a phylogenetic mean, we used a phylogenetic generalized least-squares analysis that incorporated the BM covariance structure of a phylogeny into the calculation of the mean (Martins & Hansen 1997; Garland & Ives 2000). Because we analysed each PC axis separately, the independent variable was treated as fixed to 1 for all values of the dependent variable (the PC axis). In this way, the phylogenetic mean is equivalent to the estimate of the grand mean of the dependent variable, which is also the slope intercept of the regression model (Garland & Ives 2000; Rohlf 2001). We calculated a disparity ratio from the observed disparity in the herb lineage to the disparity observed by the woody lineage. Therefore, a disparity ratio greater than 1 would indicate that the herbaceous lineage exhibits greater tip disparity than the woody lineage. It is important to note that while this is somewhat analogous to the rate analysis of O'Meara et al. (2006), we view this analysis as a more detailed description of the observed climate space occupied by a given lineage.

We also calculated a disparity ratio based on the expected tip variances predicted by the tree structure underlying each growth form lineage comparison. The expected disparity was calculated from eqn (1) in O'Meara et al. (2006) and is a function of three factors: (i) the rate of BM evolution (which, in this case, is irrelevant as it factors out of the numerator and denominator, so we arbitrarily set it to one); (ii) the path length of a tip to the root of the clade of interest; and (iii) the average entry in the covariance matrix (O'Meara et al. 2006). Larger values of factors (i) and (ii) increase disparity while larger values of factor (iii) decrease disparity. Our particular interest was the influence that speciation had on our observed disparity ratio. Speciation events occurring more recently necessarily increase the average entry in the covariance matrix and, consequently, reduce the expected disparity among tip taxa. Thus, our use of the expected disparity was meant to examine whether the observed asymmetry in the ratio of tip variance between herbaceous and woody lineages was greater than the disparity expected based solely on their respective underlying tree topology and branch lengths. Both the observed and expected disparity ratio analyses were carried out using code written by J.M.B. for R (R Development Core Team 2009). It is important to note that while this is somewhat analogous to the evolutionary rate analysis of O'Meara et al. (2006), our method should be viewed as a more detailed description of the observed climate space occupied by a given lineage, whereas the evolutionary rate analysis describes the tempo of filling that space.

(d) Estimating rates of climatic niche evolution

To assess whether there are differences in the rate of climatic niche evolution among woody and herbaceous lineages, we compared the fit of a single-rate model of BM evolution with that of a multiple-rate model. The single-rate model assumes that all lineages accumulated evolutionary changes in climate tolerance at the same rate—σ2, or the variance of phenotypic evolution—while the multiple-rate model assigns a separate rate to lineages that differ in growth form (e.g. σ2HERB and σ2WOODY). We carried out model comparisons using the ‘non-censored’ approach in BROWNIE v. 2.1 (O'Meara et al. 2006). Because the non-censored approach includes information about internal branches when estimating rate parameters, we used the likeliest state reconstructed along each edge, a procedure implemented in BROWNIE to obtain an estimate of the ancestral growth form state (e.g. woody or herbaceous) across all branches in a given tree. The best-fitting model was chosen based on the second-order information criterion, AICc (Akaike 1974; Sugiura 1978), owing to low sample size in some cases.

In phylogenies with branch lengths proportional to substitutions per site (i.e. a non-ultrametric phylogeny), branch lengths may represent estimates of evolutionary rate times time, instead of just time. If phenotypic evolution were proportional to rates of molecular evolution, the phylogenies with branch lengths proportional to substitutions per site would be a better model for variances and covariances of species than the ultrametric trees. Because molecular branch lengths apparently reflect the bias of generation time on accumulated molecular rates of evolution (Smith & Donoghue 2008), we also conducted these rate analyses to assess whether molecular branch lengths better accommodate for generation time by generally supporting a single-rate model.

3. Results

(a) The major axes of climate tolerance

Principal component analysis of the 19 BIOCLIM variables was used to uncover the primary axes of climate variation exhibited among species contained within each of the five angiosperm clades we examined. The first (PC1) and second (PC2) components explained 34.6–43.2 and 23.1–29.7 per cent of the variation, respectively, accounting for a majority of the total variation (electronic supplementary material, tables S1–S5). These axes mainly reflected temperature and precipitation seasonality gradients (figure 1).

Figure 1.

The primary axes (PC1 and PC2) of climate variation exhibited among woody (grey circles) and herbaceous (black crosses) species contained within (a) Dipsacales, (b) Apiales, (c) Moraceae + Urticaceae, (d) Primulales and (e) Commelinidae. Interpretations of the variable loadings for each axis are indicated as axis subheadings.

The variable loadings for PC1 calculated for Dipsacales, Apiales and the commelinids indicated that species separate along a gradient of temperature seasonality (figure 1). In all three clades, annual mean temperature was negatively associated with mean daily temperature range, temperature seasonality, and annual temperature range. For Moraceae + Urticaceae and Primulales, PC1 was associated with both temperature and precipitation patterns in the winter months (figure 1).

For Dipsacales, Apiales and the commelinids, PC2 indicated that species separate along precipitation seasonality. The loadings for these clades showed strong associations with annual precipitation, especially variables associated with precipitation during the warmest months. In contrast, PC2 for Moraceae + Urticaceae and Primulales was associated mainly with temperature conditions during the warmest months. In both cases, PC2 loaded strongest with maximum temperature in the warmest month and mean temperature in the warmest quarter.

(b) Growth form separation in climate tolerance

In all but one clade, there were significant differences in mean climate tolerance among extant woody and herbaceous species, but no more different than would be expected given a model of BM evolution (table 1). That is, mean differences in climate tolerance observed among growth form could have arisen by chance. Primulales was the only clade to vary significantly (Wilks's λ = 0.598, F2,298, p = 0.005; table 1), even after incorporating both chance and phylogeny. However, this is probably due to sampling, as nearly all the woody species sampled for Primulales were contained within Theophrastaceae sensu stricto, a clade of trees and shrubs restricted to the northern regions of South America. The sister families, Myrisinaceae and Primulaceae, are much more widely distributed around the northern Hemisphere. Thus, the appearance of significant niche differentiation among woody and herbaceous species of Primulales may simply reflect geography.

View this table:
Table 1.

Results from a conventional MANOVA and Monte Carlo simulations to test for significant niche differences between herbaceous and woody species relative to those expected based on stochastic BM evolution. The critical value is the 95th percentile obtained from a distribution of 1000 Monte Carlo simulated Wilks's λ-statistics assuming a gradual model of BM evolution. p-values from Monte Carlo simulations are the proportion of simulated Wilks's λ that are greater than the observed test statistic obtained from a conventional MANOVA. Italics indicate p < 0.05.

Although in general we detected no meaningful differences in overall climate tolerance among growth form, we did observe a general asymmetry in the amount of climatic niche space occupied by extant woody and herbaceous species. Lineages composed of predominately herbaceous extant species consistently occupied a greater amount of climatic niche space described by PC1 than sister lineages composed of predominately woody extant species (table 2). Moreover, the observed asymmetry in the herb-to-woody tip disparity ratio was greater than that predicted based on the underlying tree topology and branch lengths. These patterns were similar for PC2, with the exception of sister lineage comparisons within Apiales and commelinids (table 2). In the case of Apiales, species contained within Pittosporaceae + Myodocarpaceae exhibited nearly twice the variation in precipitation tolerance than the strictly herbaceous Apiaceae (table 2).

View this table:
Table 2.

Disparity ratios depicting asymmetry in the observed tip variance (MSE) observed between sister lineages that differ only in predominate growth form exhibited by tip taxa (herb versus woody). The expected ratio is the expected disparity ratio between the sister lineages based solely on the underlying tree topology and branch lengths. Observed tip variance was calculated as the MSE of the tip data measured from the phylogenetically correct mean (â) or the estimated trait value at the root node of the tree. Expected tip variance was calculated from eqn (1) in O'Meara et al. (2006). Capri, Caprifoliaceae + Heptacodium; rest D, rest of the Dipsacales; Pitt, Pittosporaceae; Myodo, Myodocarpaceae; Api, Apiaceae; Mora, Moraceae; Urtica, Urticaceae; Theo, Theophrastaceae; Prim, Primulaceae; Myrs, Myrsinaceae; Areca, Arecaceae; rest C, rest of the commelinids.

(c) Rates of niche evolution

We assessed the fit of two BM models that differed in the number of evolutionary rate parameters. For both PC1 and PC2 in each clade, using time-calibrated branch lengths, a two-rate model that inferred separate rates for woody and herbaceous species was strongly favoured over a single-rate model of BM (table 3). Herbs consistently accumulated changes in climate tolerance at much higher rates than related woody lineages; the estimated rate parameters ranged from approximately 2 to 12 times as high in herbs as in trees/shrubs. In only one case did we estimate tree/shrubs to have a higher rate. Changes in precipitation seasonality (PC2) for herbs within Apiales accumulated at a rate that was nearly three times slower than related woody lineages (table 3). Interestingly, in each clade, results obtained when using molecular branch lengths also supported a two-rate model based on growth form for PC1, with varying results for PC2 (table 4).

View this table:
Table 3.

Parameter estimates and fit of single- and multiple-rate models of BM to the major axes of climate tolerance using phylogenies with branch lengths in units of millions of years. Bold indicates the favoured model based on the sample-size-corrected AICc. Positive ΔAICc denotes significant contribution of the multiple-rate BM. σ2, estimated single-rate BM parameter; σ2HERB, rate parameter for herbaceous lineages; σ2WOODY, rate parameter for woody lineages.

View this table:
Table 4.

Parameter estimates and fit of single- and multiple-rate models of BM to the major axes of climate tolerance. Analyses were performed on phylogenies with branch lengths in units of substitutions per site (rates multiplied by 10−4 to allow for easier comparison). Bold indicates the favoured model based on the sample-size-corrected AICc. Positive ΔAICc denotes significant contribution of the multiple-rate BM. σ2, estimated single-rate BM parameter; σ2HERB, rate parameter for herbaceous lineages; σ2WOODY, rate parameter for woody lineages.

4. Discussion

Our results demonstrate that the rate of climatic niche evolution is growth-form-dependent. We hypothesize that, because growth form in plants is correlated to generation time, these differences in rates of climatic niche evolution are associated with differences in generation time. With the exception of Primulales, these results did not reflect differences in the average climate space occupied by woody and herbaceous species. Instead, the growth form dependency reflected differences in the overall amount of occupied climatic niche space. Across both climate axes (PC1 and PC2), woody species generally exhibited less interspecific variation in occupied climate space than herbaceous species of similar age (table 2). This pattern is clearly linked to differences in the rate of climatic niche evolution. We found that, in all but one of the ten comparisons, the rate of climatic niche evolution in woody lineages ranged from nearly 2 to 12 times slower than related herbaceous lineages. Even molecular branch lengths—which, if phenotypic evolution were proportional to rates of molecular evolution, would be a better model for variances and covariances of species—also strongly favoured a two-rate model for the primary climatic niche axis (table 4). Taken together, our results strongly suggest that life attributes such as generation time can impose limits to climate tolerance evolution in flowering plants.

Recently, there has been renewed interest in the extent to which niches evolve (e.g. Wiens & Donoghue 2004; Wiens & Graham 2005; Donoghue 2008). The theory of niche conservatism suggests that closely related species tend to exhibit very similar habitat requirements, implying that speciation does not necessarily result in evolutionary shifts across niche dimensions (Wiens 2004; Wiens & Graham 2005; Kozak & Wiens 2006). Donoghue (2008) suggests that, when faced with a changing environment, plant movements along corridors may be favoured over evolving key adaptations in place. Our results reveal an additional layer of complexity in understanding plant responses to changing environmental conditions. Because of the slower rates of both phenotypic and molecular evolution in woody lineages, there is an apparent constraint imposed by growth form and, by correlation, possibly generation time. Therefore, woody lineages may disperse more rapidly than accumulated mutation rates can allow for adaptation. Hence, instead of adapting to new climate tolerances, woody species may not diverge far from ancestral climate tolerances before they migrate to available niche space. This is evidenced by the clear reduction in observed climate space occupied by woody species (figure 1; table 2), while herbaceous lineages, lacking molecular and phenotypic constraints imposed by generation time, may be better able to accommodate themselves to new, emerging environments.

The tendency for woody plants to be less labile with regard to niche evolution also has implications for the evolution of the earliest angiosperms. The climatic conditions under which the first angiosperms have evolved have been the subject of considerable study, as these conditions would have implications for the evolution of a suite of morphological traits that characterize angiosperms (e.g. Stebbins 1965; Cronquist 1988; Donoghue & Doyle 1989; Taylor & Hickey 1992; Sun et al. 2002; Field et al. 2004). Hypotheses for the early angiosperm habit have ranged from shrubs (Stebbins 1965; Cronquist 1988) to herbs (Donoghue & Doyle 1989; Taylor & Hickey 1992) to aquatic plants (Sun et al. 2002; but see Friis et al. 2003). More recently, Field et al. (2004), using phylogenetic reconstructions, inferred a woody growth form as the ancestral state for early angiosperms. From our results, a woody growth habit would limit the degree of adaptation of early angiosperms to changing or fluctuating climates. In many studies, the phylogenetic branch leading to angiosperms is long, measured either by molecular substitutions or time, suggesting a low net diversification rate (Magallón & Sanderson 2001). The findings in this paper may help further our understanding of the low net diversification rate characterizing early angiosperm lineages by attributing lower rates of habitat evolution to growth habit.

Although we demonstrate the correlation of growth form and rates of climatic niche evolution, we cannot rule out other factors that would also influence the rates of climate evolution. Specifically, population size may play a significant role, as it can influence the rate at which species adapt (Lynch 2007). For example, if certain lineages exhibit larger population sizes, as some have proposed for both woody (e.g. Petit & Hampe 2006) and herbaceous (e.g. Hamrick & Godt 1996) lineages, strong stabilizing selection may inherently reduce the rate of phenotypic evolution. Alternatively, large population sizes can allow for faster responses to changing selection coefficients. However, large datasets of population size do not yet exist, and current models of phenotypic evolution do not allow for decoupling the selection coefficient (e.g. α; sensu Butler & King 2004) across character states.

Previous studies have suggested correlated rates of evolution between morphology and molecular evolution in plants (Davies & Savolainen 2006; Xiang et al. 2008). The results presented here also demonstrate this pattern, but with climate tolerance, and by additionally examining this rate correlation in the context of growth form evolution. Specifically, the same heterogeneity in the rates of molecular evolution associated with growth form is also found in climate evolution. As it relates to molecular evolution, the generation time hypothesis posits that producing more generations per unit time will lead to a higher observed nucleotide substitution rate per unit time (Wu & Li 1985; Li et al. 1987; Gaut et al. 1992, 1996; Smith & Donoghue 2008). This mechanism could link the observed differences in rates of climate evolution among growth forms. Herbaceous plants that often reproduce in the first or second year of life can exhibit very high rates of nucleotide substitution and niche evolution. Trees and shrubs typically take longer to reach reproductive age (Verdú 2002; Petit & Hampe 2006), and high rates of nucleotide substitution and phenotypic evolution are correspondingly rare.

The consistent detection of two rates of phenotypic evolution, one for each growth form, inherently violates the assumption of a single rate of BM when incorporating phylogeny in a comparative test. If phenotypic rate heterogeneity is consistent with other traits (i.e. plant functional traits), comparative tests consisting of data for both woody and herbaceous species should test for differences in the strength of trait correlations. This would effectively treat branches of differing growth form as separate in a comparative test (Revell & Collar 2009). Growth-form-dependent trait correlations may help to understand additional factors (e.g. diversification rates) also associated with the disparity in phenotypic variation among woody and herbaceous species (e.g. Beaulieu et al. 2008; table 2, this study).

Further study is necessary in order to fully understand not only the extent but also the conditions under which molecular rates and phenotypic rates are correlated. For example, an examination of among-species variation is of particular interest. In this study, we discuss large-scale patterns, and more fine-scaled analyses are needed to determine what differences are apparent in among-species variation between different growth forms. Additionally, the examination of other phenotypic traits (aside from those related to climate tolerance is important), though one limitation is the size of available datasets. For plants, many of these patterns become apparent at very large scales (Davies & Savolainen 2006; Smith & Donoghue 2008) and therefore large-scale morphological datasets are required to fully address these questions.

Acknowledgements

Valuable feedback was received from Brian O'Meara, Jeff Oliver and two anonymous reviewers. S.A.S. has been supported through the National Evolutionary Synthesis Center (NESCent; NSF no. EF-0423641). J.M.B. has been supported through an NSF Assembling the Tree of Life (AToL) award.

Footnotes

    • Received July 7, 2009.
    • Accepted August 26, 2009.

References

View Abstract