Ancestral area reconstruction

Question of interest: Can ancestral geographic distribution be inferred in an ancestral-states reconstruction?

Is legitimate use geographic occurrences, altitudinal or bathymetric range like character states?

Can one infer the origin in this way?

The initial answer to the question is, yes.

The ancestral area reconstruction in its beginning was based on the ancestral-states reconstruction method. Throughout time this approach has been changed in order to minimize mistakes when the ancestral range is estimated. Initially Bremer, in 1992, proposed the Camin & Sokal Optimization and some assumptions (e.g. the ancestral range had to be less than the actual range). Thus, the most probable reconstruction for the ancestral node is simply the state that requires the minimum number of steps to explain the distribution of states among terminal taxa given the tree.

However, according to Ronquist 1994, there are problems at using this method to estimate the ancestral area because the areas can follow fused or hierarchical history. In terms of fused history the Camin & Sokal optimization cannot fit because it is based on irreversible parsimony.

In this way, Ronquist proposed to create a character state matrix with costs to different events, assuming that “the cost of an event should be inversely related to the probability of that particular event occurring.” Another problem encountered is polymorphisms which cannot be solved via Camin & Sokal. Nevertheless, Fitch optimization for unordered states appeared to be the solution. But in this way, the dispersion plays a key role. If there is some evidence of it, the areas have to be ordered in that positions. Thus, the important part of develop an optimization method is finding a reasonable transformation series hypothesis specifying the cost of changing from one combination of distribution areas to any other combination and the cost of retaining a particular distribution (Ronquist 1994). But, if the inference is made assuming that dispersion is irreversible then, ancestral area analysis is flawed.

In 1995, Ronquist developed some ideas that are the base of a nonhierarchical approach to historical biogeography, which combines a basic assumption of vicariance with a minimization of dispersal and extinction. According to Ronquist, the optimal ancestral distributions are those that minimize the number of implied dispersal and extinction events. Thus, this new approach applies costs for the events better than a simple character optimization under assumptions as no vicariance, irreversible dispersion and ancestral range limits.

In the original approaches, vicariance was not taken into account. Instead, when Ronquist, 1997 described his method, he affirmed that allopatric speciation, associated with geographical vicariance, may be accepted as the null model in historical biogeography. The dispersion-vicariance analysis derived from character optimization methods, but in contrast to Fitch optimization allows multiple and reticulate relationships among areas.

In addition, changes of the method used for reconstruction can be observed. From parsimony to maximum likelihood (Pagel, 1999; Ree et al., 2005), where the last, compares local and global likelihood, and with two likelihoods at each node. Then the probability is computed from these. The dispersal, extinction and cladogenesis (DEC) model implemented in LAGRANGE (Ree et al., 2005) is based on a likelihood technique and in contrast of DIVA, it requires branch lengths and the time at the base node, linking the amount of change and to the biogeographic history. Also, new methods have recently been proposed, including pseudo-Bayesian versions of DIVA and LAGRANGE (e.g., Wood et al., 2013). Around the Bayesian technique a new method was developed called BayArea (Landis et al., 2013) which reconstructs geographical histories along phylogenetic branches (Matzke, 2013).

Areas of distribution

The biogeographic studies take into account the geographic distribution of organisms. These distributions are assigned to discrete areas (e.g. areas of endemism), that according to Morrone 2009 are a prerequisite for any biographical analysis. However, there have been studies over the time that use distributional ranges and not areas of endemism. So, what is the meaning of using any of these?, or, what do you have to take into account in the delimitation of the areas?

Axelius 1991, observed a problem with the delimitation of the areas and affirmed: “The result of using areas of distribution instead of separated areas may be fatal,” and one can carry on errors in the resulting area cladogram. To avoid it, the areas have to be separate entities that do not overlap. However, the initial analysis did not use areas of endemism, instead, they used isolated areas and in this way escape from erroneous area cladograms.

Functioning as separate areas, the areas of endemism were proposed for this type of analysis and defined as the congruent distribution of at least two species of restricted range (Platnick, 1991). Also, the areas of endemism represent history and primary homology (Morrone, 2001).

Altitudinal and bathymetric data

Respect to using the altitudinal and bathimetric data, there are studies that use these type of records:

Brumfield and Edwards 2007, assessed the evolution into and out of the Andes, analyzing the historical diversification in Thamnophilus. For this study, they reconstructed the ancestral area and inferred the shift from lowlands to highlands based on the elevations at which each species was most commonly observed in the field. The coding used was lowland-restricted (L, 0–1500 m), lowland-to-highland (LH, 0–3050 m), or highland-restricted (H, 500–2500 m). However, as exposed by the authors, “This discrete coding scheme does not account for variance in habitat distributions, but captures what would be considered the ‘optimal’ elevation for each.” (Brumfield and Edwards 2007). They used two approaches to reconstruct the ancestral area, Parsimony and Bayesian. But they found that the Bayesian reconstructions provided some resolution at nodes where the parsimony and DIVA analyzes were equivocal, although the Bayesian results were not statistically significant.

I think that if you can create a discrete coding for reconstructing the ancestral area ancestral using bathymetric records, you can do an analysis similar to Brumfield and Edwards 2007 but the recommendation is to look for or generate your own areas of endemism.

Miranda et al 2013 did a reanalysis of Southern Ocean areas of endemism that can be used for any study like discrete areas. Nevertheless, theoretical and practical frameworks concerning to areas of endemism are complicated in marine biogeography. This is because the nature of marine realm, the oceanic dynamics, the difficulties in establishing thresholds in ecophysiological continuum and the amazingly diverse strategies of dispersal.

But, if you want to know more about regionalization of the marine zones you can see Spalding et al., 2007, they defined the levels of classification and described how the areas are determined. The problem is that ocean boundaries shift continuously with weather patterns, with seasons, and with longer or more random fluctuations in oceanographic conditions. This study was made with the aim of strategically planning and prioritizing new marine conservation measures but taking into account comprehensive biogeographic system to classify the oceans.

Finally, you can reconstruct your ancestral areas via ancestral-states reconstruction method using any of the methods derived from it and follow different approaches such as parsimony, likelihood or bayesian inference. For your reconstruction, remember that you have to use discrete areas that do not overlap. In terms of areas of endemism, you can use those already available or generate new areas. On the other hand, you can do the reconstruction starting from altitudinal or bathymetric data if you discretize the ranges to use.

References

Axelius, B. (1991). Areas of distribution and areas of endemism. Cladistics,7(2), 197-199.

Bremer, K. (1992). Ancestral areas: a cladistic reinterpretation of the center of origin concept. Systematic Biology, 41(4), 436-445.

Brumfield, R. T., & Edwards, S. V. (2007). Evolution into and out of the Andes: a Bayesian analysis of historical diversification in Thamnophilus antshrikes.Evolution, 61(2), 346-367.

Landis, M. J., Matzke, N. J., Moore, B. R., & Huelsenbeck, J. P. (2013). Bayesian analysis of biogeography when the number of areas is large.Systematic biology, syt040.

Matzke, N. J. (2013). Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. University of California, Berkeley.

Miranda, T. P., Cantero, Á. L. P., & Marques, A. C. (2013). Southern Ocean areas of endemism: a reanalysis using benthic hydroids (Cnidaria, Hydrozoa).Latin American Journal of Aquatic Research, 41(5), 1003.

Morrone, J. J. (2001). Homology, biogeography and areas of endemism.Diversity and Distributions, 7(6), 297-300.

Morrone, J. J. (2009). Evolutionary biogeography: An integrative approach with case studies. Columbia University Press.

Pagel, M. (1999). The maximum likelihood approach to reconstructing ancestral character states of discrete characters on phylogenies. Systematic biology, 612-622.

Platnick, N. I. (1991). On areas of endemism. Australian Systematic Botany,4(1), 11-12.

Ree, R. H., Moore, B. R., Webb, C. O., & Donoghue, M. J. (2005). A likelihood framework for inferring the evolution of geographic range on phylogenetic trees.Evolution, 59(11), 2299-2311.

Ree, R. H., & Smith, S. A. (2008). Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Systematic Biology, 57(1), 4-14.

Ronquist, F. (1995). Ancestral areas revisited. Systematic Biology, 44(4), 572-575.

Ronquist, F. (1997). Dispersal-vicariance analysis: a new approach to the quantification of historical biogeography. Systematic Biology, 46(1), 195-203.

Ronquist, F. (1994). Ancestral areas and parsimony. Systematic Biology, 267-274.

Spalding, M. D., Fox, H. E., Allen, G. R., Davidson, N., Ferdaña, Z. A., Finlayson, M. A. X., ... & Robertson, J. (2007). Marine ecoregions of the world: a bioregionalization of coastal and shelf areas. BioScience, 57(7), 573-583.

Wood, H. M., Matzke, N. J., Gillespie, R. G., & Griswold, C. E. (2012). Treating fossils as terminal taxa in divergence time estimation reveals ancient vicariance patterns in the palpimanoid spiders. Systematic Biology, sys092.

A first approach to Jackknife support on Endemism Areas.

Introduction.

The endemism areas are a geographic delimited by the distributional congruence of at less two taxa (Szumik et al. 2006). Those areas are the basis for inferences on historical distributional patterns in diversification process of taxa (Morrone & Crisci. 1995). Even with its importance and its broad use, there is not a statistical metric to allow us estimate the precision of the areas calcualted, like Jacknife or Bootstrap implemented in phylogenetic analysis. Jacknife are a resampling technique used as a variance and bias estimator, and therefore as a precision estimator (Badii et al. 2007). Basically this technique create a new matrix from the original data randomly removing a percentage of observations and re-calculate the corresponding analysis, then the result generated from the resampled matrix is compared with the original data and a average is calculate (Efron, 1982; Badii et al. 2007). The main goal in this work is to attempt the Jacknife method on endemism areas and establish bases for a future functional and replicable technique. 

 ______________________________________________________________________________

Methods.

Using 20.000 occurrences for 426 Mammals species, with distribution on the North Andean Block, endemism areas were calculated using the endemicity index implemented in NDM/VNDM (Szumik et al. 2002; Szumik & Goloboff. 2004). Two types of analysis were made: 1. Using only a single replication, and 2. Using 20 replications. A 0.5° x 0.5° cell gride, 50% of similarity strict consensus and default values were implemented for both analysis type. The R platform was used to create a Jacknife function with 15, 30 and 45 resamplig replics. The Jacknife function resamples the occurrences of the species from a xyd file and for each resampling its consensus areas were calculated and compared with original areas and the recovered percentage of each area was estimated from resampling consensus. 

_______________________________________________________________________________

Results and Discussion.

Consensus areas for the two analysis were different, 4 consensus areas in the single replicas analysis (OCA) and 7 areas in the 20 replicas analysis (OCA20). The OCA overlaps OCA20 on the geographic space (Fig 1.). Only two areas of 4 possible in OCA were recovered by the Jacknife resampling, and those two areas presented a low support value (2% - 6%) for the three resampling replications (Fig2. A), the highest values (0.06 and 0.04) were in the 15 resampling replications. For the 7 possible areas in OCA20, only 3 areas were recovered, also with low support values (2% - 6%) and the highest values were in the 15 replications. Something common for the analysis is that the supported areas in OCA are the same areas supported in OCA20, and one of those areas (Purple for OCA and red for OCA20) presented the highest support in the 15 and 45 replications. The lowest values for the supported areas can be associated with the small consensus area and the few occurrences used (considering all species with the same number of occurrences, each species will have approximately 45 occurrences for North Andean Block), but despite the support values, similar areas were recognized for both analysis by Jacknife, and is clear that those areas couldn't be generated by random. There is no improvement in the support values when the number of resampling replications increase. Otherwise the support values decrease with increasing the number of resampling replications. It appears that small areas are less recovered by an increasing number of resampling replications, and those areas are supported for specific combinations of species occurrences given the low number of occurrences. The Jacknife is a potential method in endemism areas, is necessary repeat this work with bigger data, especially in the number of occurrences and probe an extreme approach, resampling species not occurrences.

_______________________________________________________________________________

Figures.

Figure 1. Geographic ubication of the consensus areas. OCA areas (Blue), overlaping the OCA20 areas (Green). The ubication is in North West Ecuador Region near to Colombia.

Figure 2. Jacknife support values for each resampling replic in both analysis: One single replic (A), and 20 replics (B).

 _______________________________________________________________________________

 Bibliography.

Badii, M. H., J. Castillo, A. Wong & J. Landeros. (2007): Presición en los índices estadísticos: Técnicas de jacknife & Bootstrap. Innovaciones de Negocios 4(1): 63-78, 2007.

Efron, B. (1982):  The Jacknife, the bootstrap, and other resampling plans. Societyof Industrial and Applied Mathmatics. CBMS-NSF Monograph, 38.

Juan J. Morrone & Jorge V. Crisci. (1995): Historical Biogeography: Introduction to Methods. Annu. Rev. Ecol. Syst. 1995. 26: 373-401.

Szumik C., F. Cuezzo, P. Goloboff & A. Chalup. (2002):  An optimality criterion to determine areas of endemism. Syst. Biol. 51: 806-816.

Szumik C. & P. Goloboff. (2006): Areas of endemism. An improved optimality criterion. Syst.Biol. 53: 968-977.

 

Szumik, Claudia, Dolores Casagranada & Sergio Roig Juñent. (2006): Manual de NDM/VNDM: Programas para la identificación de áreas de endemismo. Insituto Argentino de Estudios Filogenéticos, Año V, Vol (3).

________________________________________________________________________________

Supplements

1. R code for resampling and Jacknife Calculation avaible here.

Estimating divergence times under the methods autocorrelated variation rate and independent rate depending on the size of the topology using MULTIDIVTIME and BEAST.

Introduction

Estimating divergence times can be performed by methods based on variation of substitution rates. The AR method establishes that the rates of change between ancestor and descendant are autocorrelated and the rates can follow a log-normal distribution with a mean equal to the parent rate (Thorne, Kishino and Painter, 1998; Kishino, Thorne and Bruno, 2001; Thorne and Kishino, 2002) or can be determined by a non-central χ 2 distribution (Lepage et al., 2007). Another method used in molecular clock analysis is the RR method which assigns random independent rates for each lineage and these rates are drawn from a single underlying parametric distribution such as an exponential or log-normal (Drummond et al., 2006; Rannala and Yang, 2007; Lepage et al., 2007). 

RR method has been implemented in BEAST the AR method has been implemented in MULTIDIVTIME. In terms of the programs, other difference is the priors on node ages, where MULTIDIVTIME uses a dirichlet distribution without an explicit assumption about the biological processes (Kishino, Thorne and Bruno, 2001; Thorne and Kishino, 2002) and in the other hand, BEAST uses a Yule prior and the Birth-Death prior (Yule, 1924; Rannala and Yang, 1996; Yang and Rannala, 1997).

These two methods have been widely used to dating phylogenies and it is recommended to use several calibration points; however, it is not common to find a number of fossil equivalent to the number of nodes in the topology. At the moment, it has not been established what is the minimum number of constraints to achieve correctness divergence times.

Figure 1. Number of citations. Comparison between Bayesian inference and Maximum Likelihood programs.

Objectives

General objective

To determine what are the minimum number of calibration points related to the number of tips of a topology using the Autocorrelated rate and Random rate methods.

Specific objectives

To correlate the divergence time estimated with the divergence time simulated.

To determine the delta of variation when the number of points increase.

Methods

1. Simulations

The trees will be simulated considering four different number of tips: 10, 25, 50 and 100, plus an age of 31Mya at the root, with the package phytools v.0.4-56 (Revell, 2012) in the R language and it will be replicated 3 times. Based on the these trees, the sequences will be simulated using the program Seq-gen v1.3.3 (Rambaut and Grass, 1997), under the model HKY, with base frequencies 0.30[T], 0.25[C], 0.30[A], 0.15[G]; a transition-transversion rate K=5, gamma parameter alpha = 0.5 and a length of 1000 bp.

2. Molecular clock analysis and calibration

The analyzes will be performed using the AR and RR methods respectively in MULTIDIVTIME v9.25.03 (Thorne and Kishino, 2002) (MDT) and BEAST v.1.8.2 (Drummond and Rambaut, 2007; Drummond et al., 2012). In MDT the divergence time will be calculated after 10⁶ generations employing a correlated relaxed lognormal clock. In BEAST the analyzes will be done under an uncorrelated relaxed lognormal clock, Yule speciation process, a normal distribution for the tmrca and different number of generations depending on the authors's recommendation about the relation between it and the number of tips. For both methods the calibration points will be chosen randomly, leaving the base node of the ingroup fixed to 30 Mya. The number of points will be replicated for every tree simulation.

3. Correlation and additional comparisons

The times simulated and the outcomes of these analyzes will be compared and related by a Pearson correlation for each topology of different size and replica. The standard deviation among replicas was calculated and the delta between number of points, with the aim of observe variance in the results.

Results and Discussion

Accurate in dating with AR and RR methods varies according to the number of tips in the phylogeny; by increasing the number of tips the correlation coefficient (R²) decrease, thus the number of calibration points must be in proportion to the number of tips. It is evidenced when for the two methods there is a point from which the correlation decrease (Figure 2, a, b).

A

B

Figura 2. Variation of the correlation coefficient after increase the number of calibrated nodes. A. Autocorrelated rate method. B. Random rate method.

The addition of the second point of calibration does not improve the dating of the phylogeny in approximately 90% of the cases reviewed (Figure 2). Similar results were obtained by Battistuzi et al., 2010, they assess the effect of multiples calibrations into the credibility Intervals (CrIs) size, and they do not found differences in dating with one or two points. On the other hand, the use of exact divergence times can generate this behavior, in which differences between one and two points are not perceived (Benton et al., 2009). Consequently, the minimum number of calibration points is one, because with just one point the value of the correlation coefficient is greater than 0.89, in all cases. However, when I add three points, the coeficient (R²) was high and in almost 70% of the cases was the maximum value.

Standard deviation among replicas was of 10⁻¹, although cases of standard deviation of 10⁻³ were recorded , it demonstrated that there is no changes in using different points and multiples runs of the programs. This is consistent with Hedges et al., 2004, who found that along many methods if the data and times are constant the estimations will be similar. Both methods evaluated in this work are based on Bayesian Inference and the priors play a key role in the estimations of divergence time and in all the analyses, but Brown and Yang in 2010 tested the influence of the priors and conclude that it affect the estimation, however, one can obtain similar results.

A

B

Figure 3. Standard deviation among replicas for each number of tips. A. MultiDivTime. B. BEAST.

The delta of increasing the number of points varied in the order of 10⁻¹ until 10⁻³. Thus, demonstrating, that the differences among number of calibrations vary in a little range. Then, when the position of the unique point was assessed, the position imports, is not the same dating towards the tips that dating the base node of the ingroup, but the value of the coefficient of correlation is not below 0.75 in both cases (Figure 4). According to Battistuzi et al., 2010, there is no significant differences in the relationship of the calibration nodes and the estimated nodes, “this suggests that the relative positions of the two nodes do not affect the time estimations as long as the two nodes are closely related (Battistuzi et al., 2010).” In contrast, with empirical data, the estimations depends on the position, the method, and the dataset (Hug and Roger, 2007), idea discussed by many other authors who assess this problem via simulations and empirical data and found similar results among different methods, or the parameters of the analysis.

Figure 4. Variation of the correlation coefficient according to different positions of the calibrated node.

In addition, to reduce the error by 5% and increase the coefficient of correlation in larger topologies compared to that obtained in small topologies, with a single point calibration, a number of points four time major than the points evaluated are required. However, a solution would be dated by small separate clades within the topology in order to obtain a high coefficient correlation for all phylogeny.

Conclusions

Finally, the number of tips influence the estimation, regardless the method used in estimating divergence times. Although a point is the minimum number, when the number of tips increase with the error occurs exactly the same.

References

Battistuzzi, F. U., Filipski, A., Hedges, S. B., & Kumar, S. (2010). Performance of relaxed-clock methods in estimating evolutionary divergence times and their credibility intervals. Molecular biology and evolution, 27(6), 1289-1300.

Brown, R. P., & Yang, Z. (2010). Bayesian dating of shallow phylogenies with a relaxed clock. Systematic biology, 59(2), 119-131.

Donoghue, P., Benton, M., Ziheng, Y., & Inoue, J. (2009). Calibrating and constraining the molecular clock. Journal of vertebrate paleontology, 29, 89A-89A.

Drummond, A. J., Ho, S. Y., Phillips, M. J., & Rambaut, A. (2006). Relaxed phylogenetics and dating with confidence. PLoS biology, 4(5), 699.

Drummond, A. J., & Rambaut, A. (2007). BEAST: Bayesian evolutionary analysis by sampling trees. BMC evolutionary biology, 7(1), 214.

Drummond, A. J., Suchard, M. A., Xie, D., & Rambaut, A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular biology and evolution,29(8), 1969-1973.

Hedges, S. B., Blair, J. E., Venturi, M. L., & Shoe, J. L. (2004). A molecular timescale of eukaryote evolution and the rise of complex multicellular life. BMC evolutionary biology, 4(1), 2.

Hug, L. A., & Roger, A. J. (2007). The impact of fossils and taxon sampling on ancient molecular dating analyses. Molecular biology and evolution, 24(8), 1889-1897.

Kishino, H., Thorne, J. L., & Bruno, W. J. (2001). Performance of a divergence time estimation method under a probabilistic model of rate evolution.Molecular Biology and Evolution, 18(3), 352-361.

Lepage, T., Bryant, D., Philippe, H., & Lartillot, N. (2007). A general comparison of relaxed molecular clock models. Molecular biology and evolution,24(12), 2669-2680.

Rambaut, A., & Grass, N. C. (1997). Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.Computer applications in the biosciences: CABIOS, 13(3), 235-238.

Rannala, B., & Yang, Z. (1996). Probability distribution of molecular evolutionary trees: a new method of phylogenetic inference. Journal of molecular evolution, 43(3), 304-311.

Rannala, B., & Yang, Z. (2007). Inferring speciation times under an episodic molecular clock. Systematic Biology, 56(3), 453-466.

Revell, L. J. (2012). phytools: an R package for phylogenetic comparative biology (and other things). Methods in Ecology and Evolution, 3(2), 217-223.

Thorne, J. L., Kishino, H., & Painter, I. S. (1998). Estimating the rate of evolution of the rate of molecular evolution. Molecular Biology and Evolution,15(12), 1647-1657.

Thorne, J. L., & Kishino, H. (2002). Divergence time and evolutionary rate estimation with multilocus data.Systematic Biology, 51(5), 689-702.

Yang, Z., & Rannala, B. (1997). Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Molecular biology and evolution, 14(7), 717-724.

Yule, G. U. (1924). An introduction to the theory of statistics. Bull. Amer. Math. Soc. 30 (1924), 465-466