Estimating Evolutionary Rates and Times to the last Common Ancestor DEV-2, using partitions vs complete coding region

Introduction

Dengue Virus serotype 2(DENV-2) is one of the four serotypes of the arthropod borne viruses (arbovirus) that belongs to the Flaviridae family. As the other serotypes, it is responsible for causing Dengue fever and severe Dengue. The genome is formed by a single stranded positive-sense RNA of approximately 11kb long. This viral RNA serves as an mRNA for the translation of 10 viral proteins: 3 structural(Capside C, Membrane precursor preM, Envelope E), that assemble the nucleocapside and envelope of mature virions; and 7 non-structural(NS1, NS2a, NS2b, NS3, NS4a, NS4B, NS5), which are mainly responsible for viral replication and other cell functions.

Due to its considerable diversity, DENV-2, has been grouped into 6 lineages, designated as genotypes (Rico-hesse et al., 1997, 1998; S Susanna Twiddy et al., 2002)⁠. 1.) Asian I genotype, represents strains from Thailand, Malaysia, Cambodia, Myanmar, Vietnam and Australia. 2.) Asian II, represents strains from China, Indonesia, The Philippines, Taiwan, Sri Lanka, India, Honduras and Mexico. 3.) Southeast (SE) Asian/American genotype comprises all strains from Southeast Asia, and strains collected in Central and South America and the Caribbean over the last 30 years. 4.) The Cosmopolitan genotype, represents strains distributed in a wide geographic area, including East and West Africa, the Middle East, the Indian subcontinent, Indian and Pacific Ocean Islands and Australia. 5.) American genotype, representing strains from Central and South America, the Caribbean and older strains collected in the Indian subcontinent and the Pacific Islands. Lastly, there is the Sylvatic genotype, representing strains from humans, canopy-dwelling arboreal mosquitoes, and non-human primates collected in West Africa and Southeast Asia. Moreover, the great genetic diversity within the DENV-2 genotypes, reflect their continual divergence and diverse geographic distribution (Chen & Vasilakis, 2011)⁠.

Several studies have been aimed at understanding the epidemiology of DENV, rates and dates of evolution, and selection pressure in its different genes and genotypes (Rico-hesse & Rico-hesse, 1990; Rico-Hesse et al., 1997; S. S. Twiddy, Holmes, & Rambaut, 2003; Weaver & Vasilakis, 2009; Zanotto, Gould, Gao, Harvey, & Holmes, 1996)⁠. Further, a relatively recent review (Chen & Vasilakis, 2011)⁠ reports large scale phylogenetic studies of all dengue serotypes and discusses how factors such as rates of evolution, selection pressures, population sizes and evolutionary constraints influences DENV transmission, pathogenesis and emergence. New approaches reveal the patterns in the spread and epidemiology of viral diseases by phylogeographic studies. On a global scale, DENV-2 has been reported as the serotype with the highest effective population size (Costa, Voloch, & Schrago, 2012)⁠, and the time to the last common ancestor as 338.1 (227.1–462.9) (Patil et al., 2011)⁠. Nonetheless, all these studies have been based only in estimations made from the E gene sequences.

The objectives of this study were to infer the time to the Most recent common ancestor MRCA of the DENV2 genotypes, and estimate the evolutionary rates of the data set, using the complete genome and individual genes.

Materials and Methods

Sequence Data

Complete genome sequences of DENV-1 were retrieved from GenBank(1093 sequences with collection dates (1944 to 2014), from (48) distinct countries, database up to June 2015). In order to delete duplicates, recombinants, and chimeras, I used uclust v.4.0.38 (Edgar, 2010)⁠, with a cut-off id value of 90%. An statistical test for detection of significant recombination events in the dataset was run using SplitsTree v.4.2 (Bruen, Philippe, & Bryant, 2006)⁠. Sequences were aligned using MUSCLE v3.8.31(Edgar, 2004)⁠. From this alignment, I obtained 11 partitions corresponding to the individual genes(10), and the complete coding sequence(1). Substitution model was selected under the Akaike information criterion, using PhyML (Guindon, Dufayard, Lefort, & Anisimova, 2010)⁠⁠, called with the function phymltest(package ape, R v.3.2.1.(R Core Team 2014)).

Phylogenetic Signal

In order to test the phylogenetic signal, ergo, the topological resolution yielded by the data set of each partition, a Likelihood mapping analysis was run for 10 million random quartets using TreePuzzle v.5.3.(Schmidt, Strimmer, Vingron, & von Haeseler, 2002)⁠. The quartets were reconstructed following a ML approach under the selected substitution models for each partition. Thus, an independent analysis was run for each dataset. The results were plotted as a triangle with fully resolved quartets shown as dots concentrated near the vertices, and the unresolved ones are located in the center. (See anexes)

Genotypes Assignation

a Maximum likelihood (ML) tree was inferred using RaxML, using the complete coding sequence (CDS), with the model GTR+G+I.

Estimation of Evolutionary rates and TMCRA

Tree topologies, model parameters, evolutionary rates for each partition dataset, and times of the MCRA for each genotype were co-estimated, using BEAST v.8.1.2.(A. J. Drummond, Suchard, Xie, & Rambaut, 2012)⁠. Inferrences were performed using a relaxed molecular clock (uncorrelated lognormal[rate: 7.5*10^-4, variation: 0.4 (Costa et al., 2012)⁠]. A Bayesian skyline plot was used as a coalescent prior during the inference (a. J. Drummond, Rambaut, Shapiro, & Pybus, 2005; Ho & Shapiro, 2011)⁠, for it is an approach with no dependence on a pre-specified parametric model of demographic history, such as logistic, exponential. Analysis was run for twice for 50 million generations, sampling every 10000. Topology was rooted with a sequence of Japanese Encephalitis Virus(EF107523) from 1988.

Results and Discussion

The mean evolutionary rate for each gene ranged 5.3 10^-4 –1.5 × 10^-4 nucleotide substitutions per site per year, being the NS5, NS2B, and E, the ones with the highest substitution rates. This estimations can be compared to the ones obtained by(S. S. Twiddy et al., 2003)⁠. The complete coding sequence was not account for comparisons, because converge of this parameter is not guaranteed (ESS<100) for the partition. In general, evolutionary rates for the dataset obtained from the E gene sequence were similar compared to the non-structural genes NS2B, NS4A, NS4B, and NS5.

Table 1. Estimated Mean Nucleotide Substitution Rates per Site per Year for each partition of Dengue Virus Serotype 2. 

	Rate(subs/site/year)	Rate 95%HPD
C	4.10E-004	[9.5E-4, 2.4E-4]
preM	4.20E-004	[8.6E-4, 2.5E-4]
E	2.30E-004	[4.9E-4, 1.4E-4]
NS1	3.80E-004	[6.4E-4, 1.1E-4]
NS2A	5.30E-004	[8.1E-4, 1.8E-4]
NS2B	1.40E-004	[4.1E-4, 7.3E-3]
NS3	2.80E-004	[5.2E-4, 1.1E-4]
NS4A	2.50E-004	[6.5E-4, 1.3E-4]
NS4B	2.60E-004	[6.5E-4, 2.4E-4]
NS5	1.50E-004	[7.5E-4, 3.1E-3]
CDS	4.30E-004	[9.4E-4, 1.2E-3]

Regarding the time the most recent common ancestor (MRCA), results are consistent to the times found by (Walimbe, Lotankar, Cecilia, & Cherian, 2014), and are comparable to the times of outbreaks and introductions summarized in (Chen & Vasilakis, 2011)

In average, the partitions whose TMCRA for Asian/American genotype were closest to the proposed time in other studies (using the E gene sequence: 42.757-*44.4), were NS1:46.4, NS4A:45.6, and NS4B:43.7. For genotype Asian1/Asian2, the reference value (using the E gene sequence: 64.189-*70.45), was nearly recovered by NS1: 70.7, and preM: 67.8. For the Cosmopolitan genotype, (using the E gene sequence: 59.8-*67.9), the TMCRA was roughly recovered by NS2B: 63.885, NS4B: 60.929, and the CDS: 77.779. The TMCRA of the sylvatic2 genotype (using the E gene sequence: 121.1-*158.4), was recovered by NS5:157.1, CDS:157.2, NS1:120.1, NS2A:135.3, and NS4B: 126.9. The only exception for this case, would be the American 2 genotype, whose TMCRA estimations were different (103.3*), from the obtained in this work (56.2). * Values obtained by (Patil et al., 2011; Walimbe et al., 2014).

Genes pre-M, NS1, NS2A, NS4B, showed the closest values to the references. Thence, they could be an alternative to use when estimating evolutionary processes in DENV2. Besides, regarding the potential resolution of phylogenetic relationships between clades, these genes have shown good phylogenetic signal (>90%). Also, they have relatively smaller sizes compared to the E gene (with exception of NS1), which may be even practical when running long analyses.

Table 2. Estimated time to the most recent common ancestor (MRCA) for every partition, estimated under a Bayesian skyline coalescent tree prior.

	root	95%HDP	AA	95%HDP	A1/A2	95%HDP	Cosmo	95%HDP	A2	95%HDP	Sylvatic	95%HDP
C	2174.6	[1100.8, 3381.1]	55.1	[38.5, 74.4]	57.2	[50, 72.1]	64.9	[46.2, 87.7]	54.9	[40.1, 77.1]	100.8	[56.7, 160.1]
preM	3268.7	[1866.7, 4821.1]	55.1	[39.9, 72.5]	67.8	[54.3, 83.4]	54.9	[43.9, 68.4]	46.8	[38.3, 59.2]	107.4	[64.3, 165.1]
E	6175.5	[3780, 7200.1]	42.8	[34.5, 52.1]	64.2	[53.8, 77.5]	59.8	[45.5, 76.4]	56.7	[42.3, 78.9]	121.1	[70.2, 179.5]
NS1	5452.1	[3371.3, 7790.7]	46.4	[34.9, 58.7]	70.7	[56.2, 87.3]	56.5	[44.9, 70.7]	59.5	[40.1, 71.5]	120.1	[69.9, 181.3]
NS2A	3362.9	[1944.8, 4966.4]	43.5	[34.9, 54.5]	67.1	[54.1, 83.7]	59.1	[46.6, 74.7]	53.3	[39.5, 72.4]	135.3	[72.4, 206.5]
NS2B	2712.5	[1497.1, 4033.4]	49.9	[35.9, 66.4]	62.9	[51.4, 78.7]	63.9	[45.3, 89.1]	57.2	[43.2, 74.8]	101.4	[61.2, 156.5]
NS3	2914.1	[1502.1, 4102.2]	56.9	[38.3, 73.5]	63.7	[54.1, 76.1]	59.1	[42.8, 83.9]	58.8	[46.4, 71.6]	119.1	[58.2, 215.6]
NS4A	3989.9	[2283.7, 5881.3]	45.6	[35.3, 57.8]	65.3	[52.4, 82.1]	55.1	[43.9, 66.5]	45.1	[38.1, 53.3]	118	[69.6, 187.7]
NS4B	2454.2	[1290.3, 3688.6]	43.7	[35.3, 54.2]	60.9	[51.6, 72.4]	60.9	[46.1, 79.8]	46.1	[37.9, 58.4]	126.9	[73.1, 187.9]
NS5	5675.6	[3371.4, 7790.7]	58.6	[39.9,72.5]	77.9	[56.2, 95.1]	67.3	[52.5, 83.5]	63.5	[49.3, 76.9]	157.1	[72.4, 206.5]
CDS	10944.7	[3459.8, 25368.5]	105	[36.6, 550.8]	117.6	[53.6, 337.7]	77.8	[43.9, 172.7]	76.6	[39.5, 145]	157.2	[75, 290.6]

AA:American/Asian. A1A2:Asian1/Asian2. Cosmo: Cosmopolitan. A2: American2.

One visible problem with the analysis was that the convergence of parameter for the inference based on complete coding region sequence was not reached. Although this problem was expectable, increasing the number of generations to the double of the used in this study, would probably yield better estimations in terms of ESS, and I consider that this is a worth-doing comparison.

Conclusion

In terms of the estimation of evolutionary rates and TMRCA, NS1, NS2A, NS4B, would be good candidates to infer evolutionary processes in DENV-2.

Anexes

A1.) List of accessions and details of the sequences used in the study. (Link to Dropbox)

https://www.dropbox.com/sh/svjbqupk3ydzeom/AABx2P67bR2q4hgd-9-l2P1ea?dl=0

A2.) Substitution models for each partition

Partition	Subs Model
C	GTR+G+I
preM	GTR+G+I
E	GTR+G+I
NS1	GTR+G+I
NS2A	GTR+G+I
NS2B	TN93+G+I
NS3	GTR+G+I
NS4A	TN93+G
NS4B	GTR+G+I
NS5	GTR+G+I
CDS	*

A3.) Phylogenetic signal for each partition.

A. Gene C; B. Gene preM; C. Gene E; D. Gene NS1; E. Gene NS2a; F. Gene NS2b; G. Gene NS3; H. Gene NS4a; I. Gene NS4b; J. Gen NS5; K. Complete CDS.
RQ= Resolved quartets.

References

Bruen, T. C., Philippe, H., & Bryant, D. (2006). A simple and robust statistical test for detecting the presence of recombination. Genetics, 172(4), 2665–2681. doi:10.1534/genetics.105.048975.

Chen, R., & Vasilakis, N. (2011). Dengue--quo tu et quo vadis? Viruses, 3(9), 1562–608. doi:10.3390/v3091562.

Costa, R. L., Voloch, C. M., & Schrago, C. G. (2012). Comparative evolutionary epidemiology of dengue virus serotypes. Infection, Genetics and Evolution, 12(2), 309–314. doi:10.1016/j.meegid.2011.12.011.

Drummond, a. J., Rambaut, a., Shapiro, B., & Pybus, O. G. (2005). Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution, 22(5), 1185–1192. doi:10.1093/molbev/msi103.

Drummond, A. J., Suchard, M. a, Xie, D., & Rambaut, A. (2012). Bayesian P hylogenetics with BEAUti and the BEAST 1 . 7. Molecular Biology and Evolution, 29(8), 1969–1973. doi:10.1093/molbev/mss075.

Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research, 32(5), 1792–1797. doi:10.1093/nar/gkh340.

Edgar, R. C. (2010). Search and clustering orders of magnitude faster than BLAST. Bioinformatics, 26(19), 2460–2461. doi:10.1093/bioinformatics/btq461.

Guindon, S., Dufayard, J.-F., Lefort, V., & Anisimova, M. (2010). New Alogrithms and Methods to Estimate Maximum- Likelihoods Phylogenies: Assessing the performance of PhyML 3.0. Systematic Biology, 59(3), 307–321.

Ho, S. Y. W., & Shapiro, B. (2011). Skyline-plot methods for estimating demographic history from nucleotide sequences. Molecular Ecology Resources, 11(3), 423–434. doi:10.1111/j.1755-0998.2011.02988.x.

Patil, J. a., Cherian, S., Walimbe, a. M., Patil, B. R., Sathe, P. S., Shah, P. S., & Cecilia, D. (2011). Evolutionary dynamics of the American African genotype of dengue type 1 virus in India (1962-2005). Infection, Genetics and Evolution, 11(6), 1443–1448. doi:10.1016/j.meegid.2011.05.011.

Rico-hesse, R., Harrison, L. M., Nisalak, A., Vaughn, D. W., Kalayanarooj, S., Green, S., … Ennis, F. A. (1998). MOLECULAR EVOLUTION OF DENGUE TYPE 2 VIRUS IN THAILAND, 58(1), 96–101.

Rico-Hesse, R., Harrison, L. M., Salas, R. a, Tovar, D., Nisalak, a, Ramos, C., … da Rosa, a T. (1997). Origins of dengue type 2 viruses associated with increased pathogenicity in the Americas. Virology,230(2), 244–251. doi:10.1128/CVI.00363-08.

Rico-hesse, R., Harrison, L. M., Salas, R. A., Tovar, D., Nisalak, A., Ramos, C., … Nogueira, R. M. R. (1997). Origins of Dengue Type 2 Viruses Associated with Increased Pathogenicity in the Americas, 251(230), 244–251.

Rico-hesse, R., & Rico-hesse, R. (1990). Molecular Evolution and Distribution of Dengue Viruses Type 1 and 2 in Nature, 493, 479–493.

Schmidt, H. a, Strimmer, K., Vingron, M., & von Haeseler, A. (2002). TREE-PUZZLE: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics (Oxford, England), 18(3), 502–504. doi:10.1093/bioinformatics/18.3.502.

Twiddy, S. S., Farrar, J. J., Chau, N. V., Wills, B., Gould, E. A., Gritsun, T., … Holmes, E. C. (2002). Phylogenetic Relationships and Differential Selection Pressures among Genotypes of Dengue-2 Virus 1, 72, 63–72. doi:10.1006/viro.2002.1447.

Twiddy, S. S., Holmes, E. C., & Rambaut, a. (2003). Inferring the Rate and Time-Scale of Dengue Virus Evolution. Molecular Biology and Evolution, 20(1), 122–129. doi:10.1093/molbev/msg010.

Walimbe, A. M., Lotankar, M., Cecilia, D., & Cherian, S. S. (2014). Global phylogeography of Dengue type 1 and 2 viruses reveals the role of India. Infection, Genetics and Evolution, 22, 30–39. doi:10.1016/j.meegid.2014.01.001.

Weaver, S. C., & Vasilakis, N. (2009). Molecular evolution of dengue viruses: contributions of phylogenetics to understanding the history and epidemiology of the preeminent arboviral disease. Infection, Genetics and Evolution : Journal of Molecular Epidemiology and Evolutionary Genetics in Infectious Diseases, 9(4), 523–40. doi:10.1016/j.meegid.2009.02.003.

Zanotto, P. M., Gould, E. a, Gao, G. F., Harvey, P. H., & Holmes, E. C. (1996). Population dynamics of flaviviruses revealed by molecular phylogenies. Proceedings of the National Academy of Sciences of the United States of America, 93(2), 548–553.

Zhang, Z., Li, J., & Yu, J. (2006). Computing Ka and Ks with a consideration of unequal transitional substitutions. BMC Evolutionary Biology, 6, 44. doi:10.1186/1471-2148-6-44.

Inferring Ancestral Areas

Can they be inferred in an ancestral-states reconstruction?

Even though this question has been addressed in several ways, the answer is yes, they can! However, there are a few concepts that need to be accounted first.

Ancestral area has been defined as the center of origin of the diversification of a clade (Bremer, 1992)⁠⁠. In other words, it constitutes the ancestral ranges of distribution of a monophyletic group. Establishing this area is one of the central questions in historical biogeography, especially if the objective of a particular study is to assess the contribution of vicariance and dispersal to the speciation and distribution of a group of organisms (F Ronquist & Ronquist, 1995)⁠.

In order to achieve such purpose, two problems have to be considered when doing a biogeographical analysis: Earth history, which seeks to establish area relationships based on the phylogenies of at least two taxa inhabiting the areas of interest, thus using areas as taxa, and taxa as characters; and Taxon history, who aims to infer the biogeographic patterns and events that shaped the history of taxa (Hovenkamp, 1997)⁠. The main forethought with the latter approach, is that inferences are restricted to general patterns, and they can not be assumed arbitrarily.

Areas of Endemism

Regarding the first problem of biogeographical analysis, if one is using the taxon-as-area analogy, it presupposes the existence of discrete areas, as an example, I will call areas of endemism (Hovenkamp, 1997)⁠⁠, who are the first subject of investigation you should include in your analysis. The term “area of endemism” refers to a particular pattern of distribution delimited by the distributional congruence of at least two taxa (Platknick, 1991).

There are several methodological proposals to identify these areas, I am listing two of them, which have different theoretical backgrounds but as any other quantitative method use the distribution of species as data: Parsimony Analysis of Endemicity PAE (Morrone, 1994)⁠⁠, and Endemicity Analysis EA (C. a Szumik, Cuezzo, Goloboff, & Chalup, 2002; C. Szumik & Goloboff, 2004)⁠. PAE groups hierarchically groups area units based on their shared species, using the maximum-parsimony criterion. EA, on the other hand, identifies areas of endemism by assessing the congruence among species distributions, following an optimality criterion. The congruence between a species distribution and a given area is measured by an Endemicity Index EI ranging from 0 to 1. This proposal is implemented in NDM/VNDM (Goloboff, 2005)⁠, and currently it manages to run analyses for higher taxa (C. A. Szumik & Goloboff, 2015)⁠.

Ancestral Areas: Concepts and Models

Once you have established your areas of endemism, you can proceed with the reconstruction of ancestral areas, treating the areas you obtained as discrete characters. For this matter, there has been a long discussion about what definition, approach, and method to use. It formally started with Hennig (1950), who proposed the chorological progression rule, which assumes progression in the areas parallel to the progression in the characters in the cladogram, so that the areas inhabited by primitive species are deemed to be ancestral, whereas the areas inhabited by apomorphic species are situated far away from the center of origin. This rule is based on the assumption that peripatric speciation model is common in nature. Nonetheless, there are many exceptions to this rule.

Bremer(1992), assumed that areas or regions could be treated as binary irreversible characters we could analyze separately (each area as a character), optimizing it to a tree, using Camin-Sokal Parsimony to see which areas were most parsimoniously explained as being part of the ancestral area. Therefore, the selection criterion was that the areas that required the fewest independent losses relative to gains on the cladogram would be the ones most likely to be the central area of the clade.

This approach was ratified by Ronquist(1994). Nonetheless (Ronquist, 1995)⁠ discussed the notion of areas treated as irreversible characters, for it would be only valid if dispersal was irreversible and a region could not be subsequently invaded. Thence, he considered that allowing dispersal events to occur as unordered, reversible events would be more realistic. Therefore, he proposed Fitch's parsimony to optimize characters to trees. These approaches had many conceptual problems because they search for replicated areas in basal clades, ignoring homology and using paralogy to weigh areas and locate centers of origin (Ebach, 1999)⁠, thus tend to overestimate the areas that had less extinction processes, and fail dentifying an area as ancestral.

Dispersal Vicariance Analysis DIVA

This method was presented by (Fredrik Ronquist, 1997)⁠. It uses optimizations with reversible parsimony for estimating ancestral areas. DIVA searches ancestral areas using a three-dimensional cost matrix that gives different costs to events, minimizing the dispersal events needed for explaining the distributions. Unlike previous models, it focuses on mapping area distributions onto the phylogeny, and vicariance events have no cost, whereas dispersals and extinctions have a cost of one per area unit added to the distribution. The optimal reconstruction(s) are those requiring the minimal number of dispersal events. Since this method does not take dispersion into account, and always assumes that speciation is due to vicariance, it represents a problem if species have not followed this sort of event. Due to the aforementioned, it does not model extinction and range expansions, reason what it has been criticized (Kodandaramaiah, 2010)⁠⁠.

Despite these problems, it has been a powerful approach for inferring reticulate biogeographic scenarios that include different combination of events over time, such as the diversification in the Holartic (Sanmartín, 2001)⁠. Two statistical extensions of this model have been proposed: S-DIVA (Yu, Harris, & He, 2010)⁠, that evaluates the alternative ancestral ranges at each node in a tree accounting for phylogenetic uncertainty and uncertainty in DIVA optimization using an statistical framework, and Bayes-DIVA (Nylander, Olsson, Alström, & Sanmartín, 2008)⁠ which uses DIVA to perform reconstructions at all nodes that occur in a summary topology. Bayes-DIVA has been implemented in S-DIVA.

DEC Dispersal Extinction Cladogenesis

This is a continuous-time model for geographic range evolution that enables the inference of ancestral ranges in a likelihood framework (Ree, Moore, Webb, & Donoghue, 2005; Ree & Smith, 2008)⁠⁠. Range contractions and expansions are caused by dispersal to an unoccupied area and local extinction within an area. Given a phylogeny, the distribution of the taxa involved, and an explicit model of Dispersal-extinction and cladogenesis, dispersal and extinction rates are calculated using maximum likelihood. With this model, probabilities of range transitions are computed as a function of time, enabling free parameters in the model, rates of dispersal, and local extinction to be estimated by maximum likelihood. This model can be extended by incorporating fossil and geological information into the rate matrix, which is allowed to vary over time. Also, dating uncertainty can be accommodated by integrating DEC reconstructions over a Bayesian Inference posterior sample of dated trees.

To cite a few examples: (Smith, 2009)⁠⁠ examined uncertainty of divergence-time in a parametric biogeographical analysis of the Northern Hemisphere plant clade Caprifoliea; (Smedmark, Eriksson, & Bremer, 2010)⁠⁠, explored how uncertainty in estimated divergence times affects conclusions in biogeographical analysis, using the group Urophylleae, which has a disjunct pantropical distribution. DEC model has been implemented in Lagrange(http://www.reelab.net/home/software/lagrange/). Although it is considered a merely realistic model, it does not work efficiently when using more than 7 areas (Fredrik Ronquist & Sanmartín, 2011)⁠.

BayArea

It is a Bayesian approach for inferring biogeographic history that extends the application of biogeographic models to the analysis of problems that involve a large number of areas (Landis, Matzke, Moore, & Huelsenbeck, 2013)⁠.

S-DIVA, DEC, and BayArea are implemented in the software RASP, which offers a graphical user interface (GUI) to specify a phylogenetic tree or set of trees and geographic distribution constraints, draws pie charts on the nodes of a phylogenetic tree to indicate levels of uncertainty, and generates exportable graphical results (Yu, Harris, Blair, & He, 2015)⁠.

⁠

BioGeoBEARS

is an R package, authored by Nicholas J. Matzke, that was designed to perform inference of biogeographic history on phylogenies, and also model testing, which includes dispersal, vicariance, founder-event speciation (free parameter j), DEC, DIVA, and BAYAREA, inter alia (Matzke, 2013)⁠. The advantage of using this package is that you can compare the probabilities of each model, and measure the effects of the parameters you use for each model. (Matzke NJ, 2014)⁠, encourages to test the founder-event parameter for the speciation of Island Clades.

Using Altitudinal and bathymetric data: Could that be an alternative?

Regarding the possibility of using other sort of data, such as altitudinal and bathymetric, there are studies that use these type of data. (Yesson, Yesson, & Culham, 2007)⁠ used distribution data and inferred climate preferences to determine the potential distribution of species in the past, present and future, which they called: Phyloclimatic Modelling (Yesson & Culham, 2011)⁠. This proposal was applied in the study of the biogeography of the garden plant Cyclamen. 

This approach is similar to that of (Vasconcelos, Rodríguez, & Hawkins, 2011)⁠, who used a cluster analysis of richness, topography and climate to determine the variable that most affects the distribution pattern of Amphibians in South America, thus, delimiting a new scheme of regionalization. (Brumfield & Edwards, 2007)⁠, on the other hand, reconstructed the ancestral area and inferred the shift from lowlands to highlands based on the elevations at which each species of Thamnophilus was most commonly observed in the field. Nonetheless, as exposed by the authors, the discrete coding scheme they used did not account for variance in habitat distributions, but found the ‘optimal’ elevation for each.

In summary, you can infer the ancestral area of your clade of study via ancestral-states reconstruction method. Also, you can do the reconstruction accounting for altitudinal or bathymetryc data if you discretize the ranges to use. What you can do is first, searching dated phylogenies or date them in case you do not feel positively sure about them. Then, get the distributional data of the species. Once you have done this, establish the biogeographic areas(preferably areas of endemism), or ecoregions in which you will test your hypothesis(es). Use these areas to discretize the distribution of your taxa, and make sure these areas do not overlap. Run an evaluation of patterns and events (vicariance, dispersion, extinction), and with your results, you can infer the ancestral area(s). I personally prefer using DEC model, for it accounts on probabilities of different events. Yet, the main problem would be if you need to include a large number of areas into your analysis.

References

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

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. doi:10.1111/j.1558-5646.2007.00039.x

Ebach, M. C. (1999). Paralogy and the Centre of Origin Concept. Cladistics, 15(4), 387–391. doi:10.1006/clad.1999.0118

Goloboff, P. (2005). NDM/VNDM v. 2.5. Programs for identification of areas of endemism. Retrieved from http://www.lillo.org.ar/phylogeny/endemism/

Hovenkamp, P. (1997). Vicariance Events, not Areas, Should be Used in Biogeographical Analysis. Cladistics, 13(1–2), 67–79. doi:10.1006/clad.1997.0032

Kodandaramaiah, U. (2010). Use of dispersal-vicariance analysis in biogeography - A critique. Journal of Biogeography, 37(1), 3–11. doi:10.1111/j.1365-2699.2009.02221.x

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, 62(6), 789–804. doi:10.1093/sysbio/syt040

Matzke, N. J. (2013). Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Frontiers of Biogeography, 5(4), 242–248. doi:10.5811/westjem.2011.5.6700

Matzke NJ. (2014). Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. Syst Biol, (793117552). doi:10.1080/10635150490522232

Morrone, J. J. (1994). On the identifications of areas of endemism. Systematic Biology. doi:10.1093/sysbio/43.3.438

Nylander, J. a a, Olsson, U., Alström, P., & Sanmartín, I. (2008). Accounting for phylogenetic uncertainty in biogeography: a Bayesian approach to dispersal-vicariance analysis of the thrushes (Aves: Turdus). Systematic Biology, 57(2), 257–268. doi:10.1080/10635150802044003

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; International Journal of Organic Evolution, 59(11), 2299–2311. doi:10.1554/05-172.1

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. doi:10.1080/10635150701883881

Ronquist, F. (1997). Dispersal-vicariance analysis: A new approach to the quantification of historical biogeography. Systematic Biology, 46(1), 195–203. doi:10.2307/2413643

Ronquist, F., & Ronquist, F. (1995). Ancestral Areas Revisited. Systematic Biology, 44(4), 572–575.

Ronquist, F., & Sanmartín, I. (2011). Phylogenetic Methods in Biogeography. Annual Review of Ecology, Evolution, and Systematics, 42(1), 441–464. doi:10.1146/annurev-ecolsys-102209-144710

Sanmartín, I. (2001). Patterns of animal dispersal, vicariance and diversification in the Holarctic. Biological Journal of the Linnean Society, 73(4), 345–390. doi:10.1006/bij1.2001.0542

Smedmark, J. E. E., Eriksson, T., & Bremer, B. (2010). Divergence time uncertainty and historical biogeography reconstruction - an example from Urophylleae (Rubiaceae). Journal of Biogeography, 37(12), 2260–2274. doi:10.1111/j.1365-2699.2010.02366.x

Smith, S. a. (2009). Taking into account phylogenetic and divergence-time uncertainty in a parametric biogeographical analysis of the Northern Hemisphere plant clade Caprifolieae. Journal of Biogeography, 36(12), 2324–2337. doi:10.1111/j.1365-2699.2009.02160.x

Szumik, C. a, Cuezzo, F., Goloboff, P. a, & Chalup, A. E. (2002). An optimality criterion to determine areas of endemism. Systematic Biology, 51(5), 806–816. doi:10.1080/10635150290102483

Szumik, C. A., & Goloboff, P. A. (2015). Cladistics Higher taxa and the identification of areas of endemism. Cladistics, 1–5.

Szumik, C., & Goloboff, P. (2004). Areas of endemism: an improved optimality criterion. Systematic Biology, 53(6), 968–77. doi:10.1080/10635150490888859

Vasconcelos, T. D. S., Rodríguez, M. Á., & Hawkins, B. A. (2011). Biogeographic Distribution Patterns of South American Amphibians: A Regionalization Based on Cluster Analysis. Natureza & Conservação, 9(1), 67–72. doi:10.4322/natcon.2011.008

Yesson, C., & Culham, A. (2011). Biogeography of cyclamen: an application of phyloclimatic modelling. Retrieved from http://centaur.reading.ac.uk/20418/

Yesson, C., Yesson, C., & Culham, A. (2007). Phyloclimatic Modelling Can Estimate Ancestral Areas. Nature Precedings, (1997), 2006. doi:10.1038/npre.2007.478.1

Yu, Y., Harris, a. J., & He, X. (2010). S-DIVA (Statistical Dispersal-Vicariance Analysis): A tool for inferring biogeographic histories. Molecular Phylogenetics and Evolution, 56(2), 848–850. doi:10.1016/j.ympev.2010.04.011

Yu, Y., Harris, A. J., Blair, C., & He, X. (2015). RASP (Reconstruct Ancestral State in Phylogenies): A tool for historical biogeography. Molecular Phylogenetics and Evolution, 87, 46–49. doi:10.1016/j.ympev.2015.03.008