INTRODUCTION

Ecological niche modeling (ENM) is used commonly in biogeographic studies to predict the geographic ranges of taxa. A series of ENM programs have been developed (see reviews in Elith & others, 2006). Although algorithms and data input differ, all ENM methods reconstruct the niche of a species based on a set of known species occurrence points and the environmental conditions at those locations (Stockwell & Noble, 1992; Stockwell & Peters, 1999; Peterson, 2001). The accuracy of the niche models produced is assessed using internal tests as part of the model creation or by retaining a subset of species occurrence points not used for model creation for model validation (e.g., Guisan & others, 2007; Walls & Stigall, 2011). Much less commonly, models have been validated through secondary field work (e.g., Feria & Peterson, 2002; Raxworthy & others, 2003). Although the accuracy of various ENM methods has been well established for analyses with modern organisms (e.g., Elith & others, 2006; Guisan & others, 2007; Elith & Graham 2009), paleontological analyses have largely relied on the internal tests generated during model creation to determine model validity (e.g., Stigall Rode & Lieberman, 2005; Maguire & Stigall, 2009; Dudei & Stigall, 2010; Walls & Stigall, 2011). Sequestering data for secondary validation is often not possible when models are generated with small sample sizes, such as is typically the case when applying ENM models to the fossil record; and secondary field validation of paleontological ENM studies has not previously been undertaken. Herein, we assess the results of a previously published ENM analysis of Late Ordovician brachiopods (Walls & Stigall, 2011) using a set of species occurrence data acquired from a new donation of museum specimens and secondary field work.

Table 1.

Brachiopod species modeled per time slice and number of species occurrence data points utilized in the creation of original GARP prediction models in Walls and Stigall (2011); †, models that are not statistically significant; Vinlandostrophia = Platystrophia King, 1850, in part (new).

Specifically, this analysis examines the validity of the GARP (Genetic Algorithm for Rule-set Prediction) modeling system (Stockwell & Peters, 1999) for generating paleontological ecological niche models using sedimentological proxies for environmental variables. Although taxon sampling is often comparable between modern and fossil ENM analyses, environmental data sets differ substantially in their generation. Modern environmental data is typically acquired via download from online databases, whereas paleoenvironmental conditions must be inferred from field investigation, analysis of published stratigraphic columns, or derived from climate models (useful for Quaternary analyses only). All paleontological applications have reported results of high statistical significance based on internal model tests, but none of these analyses have validated in the field the method's efficacy in the fossil record. In this study, the ability of GARP to accurately predict species distributions in the fossil record is assessed by testing the Walls and Stigall (2011) model results against data collected post hoc from targeted field and museum collections. These results provide the first direct test of the utility of ENM methods for predicting niches and modeling geographic ranges of fossil taxa.

GARP has been used extensively with modern fauna. Applications include predicting effects of climate change on species distributions, assessing invasive species risks, testing evolutionary hypotheses, and guiding conservation efforts, such as selecting habitat for endangered species (e.g., Peterson, Soberón, & Sánchez-Cordero, 1999; Araújo & Williams, 2000; Peterson, 2003; Thomas & others, 2004; Thuiller & others, 2005; Chen, Wiley, & McNyset, 2007). Although GARP is widely employed, only a handful of prior studies (e.g., Feria & Peterson, 2002; Raxworthy & others, 2003; Stockman, Beamer, & Bond, 2006; Kostelnick & others, 2007; DeVaney & others, 2009) have attempted to validate model results in the field. Field validation for each of these analyses recorded high rates of predictive accuracy, with reported accuracy levels ranging from 75% to 96% prediction success. Notably, omission levels (false negatives) and commission levels (false positives) were generally also low, typically beneath 10 or 20%.

GARP is the most commonly used modeling program in published paleo-ENM studies (e.g., Martínez-Meyer, Peterson, & Hargove, 2004; Stigall Rode & Lieberman, 2005; Martínez-Meyer & Peterson, 2006; Maguire & Stigall, 2009; Dudei & Stigall, 2010; Walls & Stigall, 2011). GARP is designed to analyze species occurrence data sets that are assembled by sampling data that are neither uniform nor designed for statistical tests (such as data collected from museums or randomly selected outcrops) and where environmental data consist of poorly structured domains (Stockwell & Peters, 1999). These features are critical for analyses of fossil species, because geographic occurrence data can rarely be collected within a statistical sampling regime. The GARP system has been tested extensively and has been shown to achieve high accuracy with low numbers of species occurrence data, even when there are as few as five environmental parameters (Peterson & Cohoon, 1999; Peterson, 2001; Anderson, Peterson, & Gómez-Laverde, 2002; Stockwell & Peterson, 2002; Hernandez & others, 2006). Effectiveness with small sample sizes is critical in studies of fossil taxa where species occurrence and environmental data may be limited. Furthermore, GARP requires only presence data (Peterson, Papeş, & Eaton, 2007). Because absence in the fossil record is often due to undersampling or lack of outcrop availability rather than true absence, this is a key attribute required for paleontological applications. Although other modeling systems, notably Maxent (Phillips, Anderson, & Schapire, 2006), are gaining traction in paleontological analyses, their use in deep-time analyses with sedimentary proxies for environmental data has been limited to date (i.e., Malizia & Stigall, 2011; Stigall, 2011). We, therefore, restrict the analyses in this paper to a test of the GARP modeling system.

METHODS

The field validation presented in this paper examines GARP output models published in Walls and Stigall (2011). In that study, GARP was utilized to model the geographic ranges of eight brachiopod species during three time slices within the Late Ordovician C3 depositional sequence (Table 1, Fig. 1) of the Cincinnati, Ohio region (Walls & Stigall, 2011). This stratigraphic interval is characterized by gradual shallowing of sea level of the Cincinnatian shallow marine ramp system, due to the basin infilling with sediments weathered from the Taconic Highlands (Brett & Algeo, 2001). In order to study how species were affected by the relative sea level fall, the C3 sequence was divided into early, middle, and late time slices to provide a framework for assessing how the geographic range of each species changed through time.

The primary results of the biogeographic analysis of Walls and Stigall (2011) indicated that the amount of relative sea level fluctuation impacted the stability of niche parameters of brachiopod taxa differently between time slices and taxa. Data for the original analysis in Walls and Stigall (2011) were obtained from field collection, published literature, and unpublished references, including theses and field guides. From these data, five environmental basemaps were created to model the niche of each species (inferred water depth, limestone bedding style, limestone thickness, biofacies, and lithology). Between 5 and 48 spatially unique points were used to model species niches in each time slice (Table 1). Models were generated using Desktop GARP ( www.nhm.ku.edu/desktopgarp/), as detailed in Walls and Stigall (2011). The 50 best subset models were retained. In the final maps, values of 50 indicate regions of the research area where all of the 50 best models predicted the species occurrence, while values of 0 indicate areas where none of the best models predict species occurrence, and so forth.

Results from Walls and Stigall (2011) indicate a clear link between relative sea level and species ranges. As relative sea level fell from the early to middle time slice, species exhibited high levels of niche stability (average of 75% niche overlap between early to middle time slices). Because the decline in relative sea level decreased the maximum potential geographic extent of specific niches, species inhabiting these niches also experienced a decrease in geographic range. This is demonstrated by a reduction in average geographic range for species in the middle versus the early time slice (75% versus 89% of the map area occupied). From the middle to late time slice, however, the average range of brachiopod species increased (to 80% of the map area) as five of the eight species modeled exhibited greater niche evolution (average of only 49% niche overlap with middle to late time slices) in response to the fall in relative sea level. Subsequent analyses by Malizia and Stigall (2011) that used the methods of Warren, Glor, and Turelli (2008) to examine niche conservatism generated with the Maxent algorithm produced congruent results.

The present study tests the significance of the 22 geographic range models produced by Walls and Stigall's (2011) GARP analysis of the C3 sequence using subsequent field validation. A total of 18 new sites were identified to test the original GARP models (Fig. 2, Table 2, Appendix 1). Species occurrence data for these new sites were based on 10 sites identified from a new donation of specimens to the Ohio University Paleontology Collections and investigations of 8 new field locations, which were selected from regions where the previously generated GARP models suggested species occurrence but that were previously unsampled. Field validation sites were strategically selected to sample portions of the map area where 0 to 5 of the 50 best models predicted species occurrence, where 26 to 45 of the best models predicted occurrence, and where 45 to 50 best models predicted species occurrence. For the 8 new field sites, field work consisted of identification of brachiopod species in situ and construction of stratigraphic columns (see Appendix 2). Stratigraphic columns were constructed to provide a basis for correlation of the new field sites with the stratigraphic framework detailed in Walls and Stigall (2011). Each site was inspected for several hours to exhaustively determine presence or absence of brachiopod species. Additional sites were obtained by canvassing the Ohio University Paleontological Collection for species occurrence information. Absence data were only counted for sites visited during the course of field work for this study. These absences, which persisted after directed searches for generally common species, are considered to be true absences. Absences within the museum data set, however, cannot be distinguished from undersampling and are thus excluded from this study. For each species, a minimum of 5 total sites were investigated (Table 1).

Figure 1.

Generalized stratigraphic column of Upper Ordovician strata in the Cincinnati region; this study focuses on the C3 depositional sequence, highlighted in blue (adapted from Holland & Patzkowsky, 1996).

Accuracy was assessed relative to whether a species was predicted to occur at a site in more than half of the best subset models. This threshold was selected for analysis, because it is neither too broad to dampen the results nor so specific that no pattern can be assessed. For the predicted distribution models examined, there is strong agreement among the majority of the range predictions, resulting in a sharply defined region of predicted occurrence for the 26 or more of the 50 best subset maps (Fig. 3). However, inclusion of all 50 of the best subset maps would result in a substantial expansion of the predicted range that results in much of the predicted range being supported by less than 20% of the best subset models. Because of this pattern, a threshold of greater than 50% of the best subset models is more informative than the more commonly used threshold of all of the best subset models (e.g., Kostelnick & others, 2007; DeVaney & others, 2009). Alternate thresholding methods, such as calculation of the E parameter as proposed by Peterson, Papeş, and Soberón (2008), have been proposed in the literature. Incorporation of these methods, however, required generation of detailed AUC (area under the curve) analyses to determine model validity in concert with the original ENM generation. Since the intent of this study is a post hoc analysis of the accuracy of the original Walls and Stigall (2011) results, these methods are not feasible for the current project.

Model accuracy at field validation sites was analyzed, both statistically and by directly calculating accuracy, omission (underprediction), and commission (overprediction) rates. For statistical analysis, the probability of model accuracy was adjusted for the geographic area of the predicted distribution, which is similar to the method employed by Raxworthy and others (2003). Probability values were calculated by multiplying the probability of random capture of a species at a location, based on percentage of the model region a species was predicted to occupy (for verified species occurrences) or to be absent (for verified species absences) for all accurately predicted species occurrences. This method, therefore, calculates the probability that species occurrence accuracy is due to random distribution within the field area adjusted by the niche breadth of each species. Fischer's Exact Test was employed for data aggregated by time slice. Fischer's Exact Test was not employed for individual models, however, because the sample sizes for individual species are too small for that test to overcome type II errors. Statistical methods specific to ENM validity assessment, such as calculation of partial receiver operator curve AUCs (Peterson, Papeş, & Soberón, 2008) and Pearson's Low Sample Test (Pearson, Papeş, & Eaton, 2007), have been developed. However, these techniques are conducted in concert with the original niche model generation and/or operate on the raw data files from the original modelling and, therefore, are not able to be applied this post hoc analysis. Although these ENM specific methods would provide sensitive tests of model accuracy, the methods employed herein have a history in the literature and been used previously to successfully investigate accuracy of GARP models (e.g., Raxworthy & others, 2003).

Figure 2.

Location of field validation sites; circles, species occurrence localities utilized for original GARP models of Walls and Stigall (2011); triangles, sites utilized for field validation; 1, early time slice; 2, middle time slice; 3, late time slice; IN, Indiana; OH, Ohio; KY, Kentucky (new).

Table 2.

Sites investigated for field validation of GARP models by time slice. *, sites based on Ohio University Paleontology Collection; other locations are sites of field investigation; additional details in Appendices 1 and 2 (new).

RESULTS

A total of 140 species occurrence sites were investigated across the 18 localities to test the original 22 GARP models (Table 3). Of these, 54 sites of predicted species occurrences were verified to be correctly predicted. Overprediction, where the species was predicted to occur but was not found, was identified at 13 sites. For 42 sites, species occurrences were verified where that species was predicted to be absent, which indicates model underprediction. Finally, 31 sites of predicted species absence were confirmed to be correct. In summary, 85 sites (60.8%) correctly verified the original GARP prediction maps, with 13 sites (9.2%) of underprediction and 42 (30%) sites of overprediction. The Fisher's Exact Test result of the combined validation data for all 22 models is statistically significant (p = 0.002), supporting high accuracy of aggregate model results.

Statistical analyses also demonstrated high levels of predictive accuracy within the original GARP models. Statistical significance of the original models was assessed, based on the probability of a species occurring by random chance in the map area versus the probability of locating the species only where it was predicted to occur. Results indicate that 18 out of the 22 species models are significant at the 95% confidence level (Table 4).

Figure 3.

Field validation sites overlain on original GARP prediction models from Walls and Stigall (2011); darker shading, region where 45–50 of the best 50 models predicted species occurrence; light shading, region where 0–5 of the best 50 models predicted species occurrence; circles, validation sites where the species was identified as present; triangles, validation sites where the species was identified as absent. 1, Prediction map for Hebertella occidentalis for early time slice; 2, prediction map for Rafinesquina alternata for middle time slice; 3, prediction map for Vinlandostrophia ponderosa for late time slice; IN, Indiana; OH, Ohio; KY, Kentucky (new).

DISCUSSION

Statistical results of field validation analysis demonstrate the original GARP models were accurate at predicting species occurrences. Eighteen out of the 22 GARP models were shown to be statistically significant at the 95% confidence level (Table 4), with 61% of all the sites investigated correctly predicting species occurrences or absences (Table 3). Some species, however, have higher predictive accuracy than others. Of the 8 species models that were field validated, all models for the species Dalmanella meeki, Hebertella occidentalis, Vinlandostrophia laticosta, V. ponderosa, and Zygospira modesta were demonstrated to be statistically significant. Conversely, one model each for V. cypha, Gnamptorhynchos auburnensis, and two models of Rafinesquina alternata were not statistically different from random. Model accuracy can be explained primarily by 4 parameters: data density, niche breadth, species abundance, and taxonomic accuracy.

The primary control on model accuracy was the number of species occurrence points available for model generation. Within our tests, models generated from larger numbers of spatially distinct species occurrences exhibited higher levels of accuracy in field validation tests (r = 0.42, p = 0.05). Of species with statistically significant models in in all three time slices, the percentage of occurrences accurately predicted was higher for Vinlandostrophia ponderosa (average 80%) and Zygospira modesta (average 63%) than the other three species (averages 45%–60%) (Table 3). The original models for these 2 species were constructed with the greatest amount of species occurrence data of the five most accurately modeled species (Table 1). Vinlandostrophia ponderosa had the highest number of species occurrence points for the original models (16 to 48), whereas Zygospira modesta models were based on 7 to 14 species occurrence points. In comparison, the original model for Vinlandostrophia cypha was constructed with 5 species occurrence points in each time slice and was not statistically significant for the late time slice (Table 1). Computer modeling has shown GARP to be able to achieve high predictive accuracy with as few as 5 species occurrence points (Peterson & Cohoon, 1999), although it performs most reliably with greater than 10 species occurrences (Hernandez & others, 2006). In our field validation, 4 of the 5 niche models constructed with only 5 species occurrence points were statistically significant (Tables 1, 4); however, they are prone to presence/absence accuracy errors (Table 3).

Although models with more data produce more accurate results in general, the minimum number of discrete data points required for successful models is related to niche breadth. Notably, 2 of the 3 models generated for Rafinesquina alternata were not statistically significant, even though they included more discrete species occurrence points than many other models in this study that were statistically significant (Table 1). Rafinesquina alternata is a very common species and was collected at every site visited. This species has previously been interpreted as an ecological generalist with wide environmental tolerances (Holland, 2005; Stigall, 2010, 2011), and R. alternata occupies a larger percentage of the study region in each time slice than any other species (Table 4). This flexibility to inhabit most paleoevironments in the study area makes it difficult for the ENM program to model a constrained set of conditions that encompass the ecological niche of the species. Even though R. alternata is represented by between 17 and 25 species occurrence points per time slice, the inaccurate models included the least data, with 17 and 22 species occurrence points. Even though accurate range models were produced for specialist taxa (e.g., Dalmanella meeki) with as few as 5 species occurrence points, it is possible that generalist species may require at least 20 species occurrence points for the niche model to be accurate within the framework of this analysis. Analyses of modern taxa (e.g., Raxworthy & others, 2003; Hernandez & others, 2006; Guisan & others, 2007) also demonstrated difficulty in modeling broadly distributed species, compared to narrowly adapted species.

Another factor influencing model accuracy is the relative abundance of each species within the ecosystem and, consequently, at the outcrop. With abundant species, the effects of sampling bias are reduced compared to uncommon or rare species (Westrop & Adrain, 2001; Finnegan & Droser, 2005; Tarver, Braddy, & Benton, 2007). If a species is rare, sites of known species occurrence usable for niche modeling may comprise very few locations. These locations could represent chance collections within what may be a much larger range, geographically and ecologically. Therefore, it is likely that analyses pertaining to rare species will underpredict species occurrence due to undersampling of the ecological variation that exists within the full geographic range of the species. In our analysis, underprediction was documented in 30% of species models investigated (Table 3). Species with fewer discrete occurrence points in the original G ARP models more commonly exhibit underprediction (compare Tables 1 and 3). For example, field validation of Dalmanella meeki range models recovered sites of underprediction in all three time slices (Table 3).

Table 3.

Field validation data and error rates; predicted presence or absence is relative to the area where 26 or more of the 30 best subset model predictions of Walls and Stigall (2011) predict species occurrence; additional details in Appendix 3 (new).

Conversely, common species should be expected to have fewer sampling biases and have more accurate range models, which is congruent with our model results. For example, Vinlandostrophia ponderosa and Zygospira modesta are species that are both easily identifiable and are abundant in the stratigraphic layers in which they occur (Holland, 2005). As noted above, these species models exhibit the highest accuracy rates and were verified to occur at all sites of predicted species presence, regardless of time slice (Table 3). Furthermore, all range models for these species were shown to be statistically significant. However, these species also occurred at some sites where they were predicted to be absent in 26 or more of the 50 best models. Of the 40 sites investigated for these 2 species, the species occurred at 34 sites. Twenty-eight sites visited verified the predictions of the original GARP models, whereas 12 sites yielded specimens in locations where they were not predicted to occur. This suggests that for very common species (much like rare species), underprediction may be more of a concern than overprediction of species ranges for ecological niche models.

An additional factor contributing to model validity is the taxonomic accuracy. If species included in this analysis are paraphyletic or polyphyletic groups rather than discrete biological entities, then no coherent niche exists for the set of organisms represented by the species occurrence data set, and poor model results should be expected. This likely contributed to the poor model results of Vinlandostrophia cypha and Gnamptorhynchos auburnensis. Each species exhibited one statistically invalid model; G. auburnensis in the middle time slice and V. cypha in the late time slice (Table 4). Although the poor predictive ability of these models is due partly to low sample size, as discussed above, it is also likely a factor of taxonomic confusion. Vinlandostrophia cypha and G. auburnensis are both problematic species. In fact, the entire genus Vinlandostrophia Zuykov & Harper, 2007 (=Platystrophia King, 1850, in part) is taxonomically unstable. Cummings and Mauck (1902) determined there is no single character or combination of characters that can be utilized to separate specimens referred to G. auburnensis, V. laticosta, and V. dentata into distinct species. Alberstadt (1979) determined that many species within the genus Vinlandostrophia, including V. cypha and V. laticosta, were so similar morphologically that a complete reevaluation of the genus was necessary. Zuykov and Harper (2007) also indicated the need for revision of species assigned to the genus Vinlandostrophia. This taxonomic confusion can obfuscate attempts at creating accurate niche models by polluting the data input to the models.

CONCLUSIONS

The utility of ecological niche modeling, specifically the GARP modeling system, is confirmed as a useful tool for paleobiogeographic studies through field validation of brachiopod species in the Late Ordovician units of the Cincinnati Arch. Eighteen out of the original 22 species range prediction models of Walls and Stigall (2011) yielded statistically significant results based on validation by secondary field work. Over 60% of field sites investigated for species occurrences were correctly predicted by the niche models, and model errors were more commonly underprediction errors (30% of sites) than overprediction errors (9% of sites). Several factors promoted higher accuracy among niche models. First, most of the species utilized in the study are easily identifiable and fairly widespread, such as Zygospira modesta and Vinlandostrophia ponderosa. During secondary field validation, these species were correctly identified at every site where the species were predicted to occur in the GARP models. Another key factor contributing to model success for particular species is the number of species occurrence data points used for constructing the original GARP models. Models produced with a larger amount of original occurrence data are more likely to be statistically significant.

Three species, Vinlandostrophia cypha, Gnamptorhynchos auburnensis, and Rafinesquina alternata, were shown to have significant model error. The poor predictive results of V. cypha and G. auburnensis can be explained by taxonomic confusion and undesampling. The validity of these taxa as biological species has been questioned, and they are likely to be locally prolific subspecies or community variations. Also, results of this study, as well as Raxworthy and others (2003), Hernandez and others (2006), and Guisan and others (2007), demonstrate that species with broad environmental tolerances can yield poorer statistical results than specialist species. Rafinesquina alternata is a widely distributed species for which two of the three models were not statistically significant. This may be explained by too limited a set of initial data. Widely distributed species may require more thorough sampling to yield statistically significant results. In this study, a sample size of 25 localities yielded significant results, whereas the models based on 17 and 22 locations were not significant.

Table 4.

Statistical significance of original GARP models based on field validation. Area of prediction represents percentage of map area covered by 26 or more of the 50 best subset model predictions of Walls and Stigall (2011); bold, values that are not statistically significant (new).

This study is the first direct test of the accuracy of paleontological environmental niche models. Even though ENM methods have been widely used with extant taxa and are increasingly utilized in paleontology, this study is the first to perform a post hoc validation of ENM models that utilized sedimentary data to infer environmental parameters. Results of this analysis validate ecological niche modeling in general and the GARP modeling system in particular as useful tools for paleontological investigation. It clearly demonstrates that utilizing fossil occurrence data (sourced from new field work, museum collections, or databases) with environmental parameters estimated from sedimentary parameters can provide a detailed quantitative assessment of paleobiogeographic patterns. These methods can provide a framework for researchers to study the biogeography of fossil organisms and assess how environmental changes have affected the evolution of life through time.