(a) Step 1: estimating the bioclimatic envelope of bison through time
Estimating a species's bioclimatic envelope is an important first step for reconstructing and understanding its past distribution [13,17]. We included climate data from multiple time periods, which allowed us to compare estimates of bison BEMs calculated independently ‘within’ each time period against ‘pooled’ estimates obtained using the full fossil record. The ‘within period’ analysis refers to climatic niches calibrated using only data from within each time period. The ‘pooled period’ used all the data to create a conglomerate of the climatic niche conditions experienced by the species throughout the Late Quaternary, which was then projected to each time slice. The advantage of this approach is that it limits possible spatial and environmental bias, but it requires accepting that niches are conserved over time (but see ).
Two palaeoclimatic simulations were performed to represent the climatic conditions during the Marine Isotope Stage 3 (MIS 3): the warmer middle part, around 42 thousand years ago (ka) and the colder later part, around 30 ka. We also used one simulation for the Last Glacial Maximum (LGM; approx. 21 ka) and one for the Mid-Holocene (approx. 6 ka). Carbon dioxide levels were specified at 200 ppm for the MIS 3 and LGM simulations , and 280 ppm for the Mid-Holocene simulation . The 0 ka simulation is pre-industrial and built using the same general circulation model, to ensure temporal and spatial comparability of our climatic surfaces and BEM outputs. Sea surface temperatures (SSTs) for the MIS 3 and LGM simulations were taken primarily from CLIMAP , with modifications from GLAMAP-2000 and other sources . SSTs for the Mid-Holocene simulation were prescribed at present-day values . In all cases, insolation was calculated using orbital parameters [24,25]. All simulations were spun up to equilibrium; results are 10-year averages. Areas known to be under ice sheets given palaeoclimate models were masked as unsuitable habitat. Although palaeoclimatic simulations based on different AOGCMs might differ, previous studies comparing the effect of AOGCM (GENESIS v. 2 versus HadCM3) on modelled ranges show that trends in range size are highly correlated for six different megafauna species, including the bison . Together, these climate data provided five time periods (42, 30, 21, 6 and 0 ka) for estimating BEMs.
A comprehensive set of fossil and historic bison localities was assembled from multiple sources [7,26]. Fossils were calibrated using the IntCal09 calibration curve , available through OxCal online (http://c14.arch.ox.ac.uk/). Georeferencing was accomplished using Google Earth and according to best practices as defined in Chapman & Wieczorek  (see the electronic supplementary material, table S1).
Fossils were considered contemporaneous with a climate layer if the calibrated estimate of the radiocarbon-dated fossil was within ±3000 years, following Nogues-Bravo et al. . We note that the bin of ±3000 years may not be appropriate for all climate layers and should not be considered as a hard and fast rule of our suggested framework. Three climate predictors were used: average minimum temperature of the coldest month (tmin), average maximum temperature of the warmest month (tmax) and mean annual precipitation sum (pre). Following ensemble forecasting methodologies , we fitted models with a fully factorial combination of predictor variables allowing an exploration of the resulting range of uncertainties .
All models were fitted with BIOENSEMBLES, a platform for computer-intensive ensemble forecasting of bioclimatic models (e.g. ). Models included three presence-only methods (BIOCLIM, DOMAIN and Mahalanobis ), two presence-background methods (MaxEnt  and GARP ) and four presence–absence methods (GLM, GAM, MARS and GBM). As we do not have true absences in our data, we generated randomly selected pseudo-absences across the cells in the region of interest without records of bison while keeping prevalence constant at 0.13 (the value found in the pooled dataset).
We randomly split the fossil and contemporaneous distributions data into 75% for calibration and 25% for evaluation, repeating the procedure 10 times. Every model run yielded a projection; ‘True Skill Statistics’ (TSS) measured the matching between predictions and observations in the 25% evaluation data. TSS-weights, indicating model performance, were obtained for every model run and eventually used to weight the different models for their ability to predict the data. For each dataset (five periods and one ‘pooled’ set), bison data were modelled using nine model types × seven variable combinations × 10 cross-validated samples for a total of 630 model runs per dataset (i.e. 3780 model runs in total).
To generate a consensus across all individual model projections, first we removed all poorly performing projections (i.e. with TSS < 0.4 in the evaluation data) . Then, we overlaid the remaining projections and considered a site suitable if models agreed at least 40% of the time. This is an arbitrary measure of agreement among models that is less conservative than the 50% consensus threshold used in several forecasting studies (e.g. [30,35]), and which provided reasonable results in tests using artificial data (F. G. Guilhaumon & M. B. Araújo 2012, unpublished data; see also ). With presence–absence methods and MaxEnt, we used the 0.13 prevalence value as the cut-off to convert probabilities or continuous suitability scores (from 0 to 1) into estimates of presence and absence [37,38]. With the distance-based presence-only methods, we used thresholds usually fixed in the literature at 0.95 for BIOCLIM and DOMAIN, and 0.75 for Mahalanobis, while for GARP we used default options to set the threshold internally. Specific details on the parametrization of each model are provided (see the electronic supplementary material, table S2). We used the binary thresholded results (figure 2) from both the ‘within’ and pooled’ probability of occurrences to help calibrate coalescence models.
BEMs predicting suitable habitat for B. bison over five time periods: 42, 30, 21, 6 and 0 ka. The areas for which over 40% of models resulted in suitable habitat are shown in dark shading on the left side of each panel and the overall probability of an area containing suitable bison habitat is shown on the right of each panel. (a) The ‘within period’ refers to climatic envelopes calibrated and projected within each period of time. (b) The ‘pooled period’ uses all the data for each period to create a summary environmental niche space that is projected to each time slice. Fossil localities are shown as black dots. (Online version in colour.)
(b) Step 2: demographic model set-up for bison
The first step to create demographic models from the BEMs involved examining the number of separate populations at each time slice and the relative geographic extents of those populations through time. Criteria for defining a population were: (i) a continuous set of cells separated from other such continuous groups by modelled unsuitable habitats; (ii) clear evidence that bison had existed at locations at some point in the past; and (iii) that enough fossil evidence and ancient DNA were available to generate population parameters usable for demographic model testing.
The approach so far only considers climate drivers in the creation of demographic models. However, biotic interactions (e.g. predation) can dramatically impact demography, and this is especially true with regard to bison. Archaeological evidence supports increased human occupation , decreased numbers of bison fossils  and bison hunting from Alaska to New Mexico [40–42] starting around 10 ka. Most evidence of large-scale bison hunting, including communal bison hunting with use of corrals and jumps, occurred in the Mid-Holocene, beginning at approximately 5–6 ka (reviewed by Bamforth ). For example, the oldest corral site is in Scoggin, Wyoming and dates to approximately 5.2 ka . Furthermore, jump sites, which date back to a similar time period, became increasing frequent after approximately 3.2 ka . Together, evidence of bison hunting suggests that most large-scale hunting occurred after approximately 5 ka and into the Late Holocene. The arrival of European settlers on bison populations during historic times is perhaps an even greater impact on Bison demography. The introduction of firearms and the European horse in approximately AD 1700 resulted in large-scale slaughter for private and commercial interests, which ultimately almost drove the species to extinction by the end of the nineteenth century [43,44]. We used all the information above to create alternative demographic models that represent either ‘within’ or ‘pooled’ BEM results, along with different biotic interactions. The models are described in more detail in the Results section, as well as in figure 3 and table 1.
Demographic models based on BEMs (M1 and M2) and a combination of BEMs and additional data (M3 variants). Each model illustrates the number of populations (number of discs) and relative size (width of disc) through time. Model 3 includes six variants that include large-scale bison hunting by Native Americans starting around 5 ka (M3a), large-scale bison hunting by European settlers (M3b) and a combination of both hunting events (M3c). Each hunting scenario is included in a model based on M1 (e.g. M3a1, M3b1 and M3c1) or M2 (e.g. M3a2, M3b2 and M3c2). (Online version in colour.)
Description of demographic models tested with modern and ancient genetic data. Model 1 is based on BEMs calculated independently ‘within’ each time period, whereas model 2 is based on a conglomerate of ‘pooled’ climatic niche conditions experienced by the species throughout the Quaternary. Model 3 variants are based on models 1 and 2, but include potential bison population declines owing to hunting by humans.
(c) Step 3: genetic data and analysis
Ancient DNA sequence data for bison were originally published and analysed by Shapiro et al.  and Drummond et al. . To estimate the probability of each of the three demographic models, we used 615 bp of control region mitochondrial DNA sequence data for 159 North American, contemporary, historic and radiocarbon-dated bison samples spanning 60 000 years .
We used the software program Bayesian Serial Simcoal (BayeSSC)  to simulate 500 000 iterations of the different demographic scenarios (see the electronic supplementary material, input files for prior distributions) and the ‘abc’ R package  to estimate demographic parameters and determine the best-supported demographic model. We chose an approximate Bayesian computational (ABC) setting  to determine which demographic model was best supported. In general, estimates obtained with full likelihood-based approaches should be more reliable than ABC estimates because they use information from the complete data rather than summarized statistics. However, ABC is more flexible, can compute multi-population demographic models in a reasonable amount of time and computational power, and can be used for direct model selection [47,48]. To model demographic scenarios using BayeSSC, populations were grouped in different statistics groups using age ranges (in generations) that reflected the BEM time frames described in step 1 (i.e. 42, 30, 21, 6, 0 ka ±3000 years) and the time frames between BEMs (before 45 ka, 33–39 ka and 9–18 ka). For time periods with two populations, fossils were either assigned to the northern or southern population depending on their location (see the electronic supplementary material, tables S1 and S3, and input files).
We chose segregating sites, nucleotide diversity and pairwise Fst as summary statistics for the analysis. This adds up to a total of 28 summary statistics. In general, more information can be added by increasing the number of summary statistics; however, too many summary statistics add stochastic noise to the analysis, and thus increase the error estimating the distance between empirical to simulated data during the regression step [47,49]. We thus used an algorithm introduced by Blum & Francois , which is based on nonlinear regression and uses neural networks to optimize the dimensionality. The choice of the tolerance level used in the analysis can have a strong impact on the demographic model results, and thus we used different levels in our analysis. We used tolerance levels of 0.002, 0.004 and 0.008, thereby accepting the 1000, 2000 and 4000 closest values, respectively, for the ABC parameter estimations. Expected deviance according to the deviance information criterion (DIC) , implemented in the R package ‘abc’, was applied to infer the best-supported demographic model. We simulated the respective demographic models with 1000 iterations using the parameters of each iteration separately (sampled from the posterior distribution) as fixed model parameters. The summary statistics were then used to calculate DIC values.
To determine meaningful upper limits for the modern population size, we first performed initial runs using broad uniform priors ranging up to 30 000 000. Excessively broad priors can destabilize ABC parameter estimation. We simulated 1000 datasets and performed an ABC analysis, which showed that the higher values resulted in unreasonably high estimates of genetic diversity and that the ABC analysis clearly favoured smaller modern population size values. Thus, the upper limit for the modern population size was refined to 100 000, which resulted in a much higher effective number of simulations. This estimate was used in the full 500 000-iteration simulations.
An Untangled Problem for Evolution: DNA Topoisomerase.
BackgroundDNA replicates itself. That is Biology 101. The process is quite complex and any biology student required to rattle off the procedure on an exam can confirm. Although there are separate kinds of topoisomerase,1 they perform the same essential function within the process of replication.
What is topoisomerase?In each human cell, there are approximately 2 meters of DNA compacted within a nucleus of the diameter of about 10 micrometers.2 For perspective, imagine that the nucleus is represented by a standard basketball. The length of all the DNA compacted into the ball would go round-trip between the earth and the moon just shy of twice. Topoisomerase has the job of ensuring that the shape of DNA is manageable. Add in the fact that DNA is not just a helix, but a double-helix. As one might imagine, the task of trying to do anything productive with this material presents very serious concerns that need to be dealt with. During replication, the topology (i.e. the shape and structure) of DNA and the movement in unzipping result in issues like tangling and kinking, so the process faces a formidable and complex challenge of bioengineering.3 It’s topoisomerase’s job to alleviate those issues so they don’t halt the replication process.4 If DNA cannot replicate, then an organism cannot grow. Additionally, DNA doesn’t just unzip for replication. It unzips in a process called transcription, which, with the help of RNA, is the first half of producing vital proteins in the cell.
Two examples might help illustrate some of the mechanical concerns that arise. Without even trying, the ordinary use of a wired phone produces coiling in the helical cord. No extraordinary spinning or energy is required to produce the effect, but it doesn’t take long for the cable to go from fresh and neat to knotted spaghetti, which seem to never come out. Even if the tangle is gone, there always seems to be a small bend or a missing turn were the tangle was. The second illustration may serve to demonstrate the problem of kinking just prior to cutting the DNA strand in two. To demonstrate the concept, you can take a piece of twine and try to unwind it by pulling away both its twisted strands away from each other. A "Y"-shape forms: the two strands being pulled apart are each the top segments of the "Y" and the rest of the twine is the bottom. Now at some point, enough tension will build up at the center of the "Y" to where you can no longer pull the strands apart. Fortunately, topoisomerase is there to work out these sorts of problems. Its process is so vital and complex that some have referred to topoisomerase as the "magicians of the DNA world"5 because of their uncanny ability to manipulate the DNA strand and master its topology.
The Problems for evolution
1: Does not fit early evolutionary models for an ancestral topoisomeraseIt is still unclear how exactly topoisomerase originated on an evolutionary model, particularly because all the different types of topoisomerase appear to have originated independently of one another.6 A core tenet of evolutionary biology is that all life can be traced to a common ancestor (more technically, LUCA: last universal common ancestor). It is natural to assume that LUCA would have some master topoisomerase that would have developed into the varying domains7 of natural life: Achaea, Bacteria, and Eukaryota. But, this seems to be an incorrect inference that does not align with the actual data.8 It is speculated that topoisomerase might have emerged with either an ancient viral lineage or a LUCA that began with RNA genetic information or a combination of the two.9
These RNA-infused models would be controversial, because such scenarios involve the already contested10 (among researchers in that field) RNA-world hypothesis for the origin of life. And, the progressive introduction of topoisomerase and other necessary replication proteins into the evolutionary lineup seems an inventive proposition, but requires massive amounts of biological informational leaps that are unfounded, especially within the limited time window for the emergence of life on earth that evolution demands.
2: Conceptual Paradox11Most of all, the origin of topoisomerase presents a chick-and-egg paradox. To replicate itself the DNA molecule needs topoisomerase to unwind it. And, DNA would need to code for a protein to do the unwinding (in order to pass on the genes to code for topoisomerase). It’s like having a car with no gas. You need to drive to the gas station to get gas, but you need gas to get to the gas station. DNA might have obtained topoisomerase from somewhere else, but the instant it makes a copy of itself, the other strand needs topoisomerase also. The scenario requires the leap to the reality, where DNA actually codes for this protein itself. That is not to mention that the actual mechanics of unwinding get even more complicated, because topoisomerase is able to cut and re-join small sections of DNA, which require precise biochemistry to ensure the right pieces get linked back up with each other and at the proper speed. Too slow, and the DNA doesn’t unwind fast enough for cell division (mitosis) and the cell dies before replication can occur.12 Too quickly, and the cell can set off its own self-destruct sequence (apoptosis) or cause irreversible damage to the genes that results in cancer.13
ConclusionLastly, it is also important to keep in mind that this is only a piece of a wider issue for evolutionary biology to explain DNA replication in general: proteins that function as stabilizing clamps, ones that actually unzip the DNA strand, others that error-check the new copies, or the fact that one side of the DNA is duplicated backwards and in sections…
References2. Joseph E. Deweese and Neil Osheroff, "The DNA cleavage reaction of topoisomerase II: wolf in sheep’s clothing," Nucleic Acids Research 37, No. 3 (2009): 738-748.
4. James J. Champoux, "DNA topoisomerases: structure, function, and mechanism," Annu Rev Biochem 70 (2001): 369–413.
5. James C. Wang, "Cellular roles of DNA topoisomerases: a molecular perspective," Nature Reviews Molecular Cell Biology 3 (June 2002): 430-440. Wang also makes use of this phrasing in his conclusion of his ‘Minireview’ in "Topoisomerases: Why So Many?" in The Journal of Biological Chemistry vol. 286, 11 (April 1991): 6659-6662.
6. Patrick Forterre and Daniele Gadelle, "Phylogenomics of DNA topoisomerases: their originand putative roles in the emergence of modern organisms," Nucleic Acids Research (2009), 1. The two coauthored a similar article with Simoneta Gribaldo and Marie-Claude Serre, called "Origin and evolution of DNA topoisomerases," in Biochimie (April 2007): 426-446.
7. A domain is the largest formal classification of living organisms on earth. By contrast, a species is the smallest formal classification.
8. Forterre and Gadelle, 1, 12 (and throughout the article in fact).
10. Patrick Forterre, one of the authors mentioned in this article, notes the RNA-world hypothesis, but references its implausibility on page 146 in the chapter entitled Origin and Evolution of DNA and DNA Replication Machineries of "The Genetic Code and the Origin of Life," published by Kluwer Academic and Plenum Publishers. Jonathan Filée and Hanny Myllykallio.
11. This paradox seems lost on some. The article by Allyn J. Schoeffler and James M. Berger entitled, "DNA topoisomerases: harnessing and constraining energy to govern chromosome topology," Quarterly Reviews of Biophysics 41, 1 (2008): 41–101, puts forth: "Understanding topoisomerase specialization is necessary to illuminate the evolutionary interplay between supercoiling homeostasis and the requirement for multiple topoisomerase activities to shepherd DNA through transcriptional, replication, and recombinational processes." The authors acknowledge that topoisomerase is subject to the mechanisms of evolution. How exactly does a mechanism necessary for evolution to work – even at the basic level – get thrown into the mix if it’s not functional with the first strand of DNA? That is not to say that it cannot improve or adapt with time, but it must work from the onset to begin with. And not just it, but DNA and its other requisite enzymes and processes.
12. Deweese and Osheroff, 741-742.
13. On the other hand, however, the essential nature of topoisomerase to cell replication (which requires DNA replication) has been exploited in anti-cancer drugs, which target the topoisomerases within the unwanted cells.
Image entitled "Replication complex" by Boumphreyfr - Own work. Licensed under CC BY-SA 3.0 via Wikimedia Commons.