Introduction
Determining the correlates of extinction is crucial to understanding macroevolutionary processes operating on geologic timescales (Jablonski Reference Jablonski1986; McKinney Reference McKinney1997; Kiessling and Aberhan Reference Kiessling and Aberhan2007; Payne and Finnegan Reference Payne and Finnegan2007; Meseguer et al. Reference Meseguer, Lobo, Ree, Beerling and Sanmartín2015; Saupe et al. Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015) and for identifying taxa at potential risk of extinction today (Lee and Jetz Reference Lee and Jetz2011; Finnegan et al. Reference Finnegan, Anderson, Harnik, Simpson, Tittensor, Byrnes, Finkel, Lindberg, Liow, Lockwood and Lotze2015; Kiessling et al. Reference Kiessling, Raja, Roden, Turvey and Saupe2019; Smits and Finnegan Reference Smits and Finnegan2019). The International Union for Conservation of Nature Red List, an indicator of the global conservation status of biological species, uses both abundance and geographic range size as determinants of extinction risk for modern species (IUCN 2001). Abundance and geographic range size are often positively correlated (Gaston Reference Gaston1994; Gaston et al. Reference Gaston, Blackburn and Lawton1997; Holt et al. Reference Holt, Lawton, Gaston and Blackburn1997; Harnik Reference Harnik2011) but are not always equally important determinants of extinction risk (e.g., Kiessling and Aberhan Reference Kiessling and Aberhan2007; Payne et al. Reference Payne, Truebe, Nützel and Chang2011; Harnik et al. Reference Harnik, Simpson and Payne2012). For example, examining Neogene bivalves in the Pacific, Stanley (Reference Stanley1986) found that abundance, and not geographic range size, was a strong predictor of extinction selectivity. Similarly, Payne et al. (Reference Payne, Truebe, Nützel and Chang2011) showed an inverse association between abundance and extinction risk that was largely independent of geographic range size. Others, however, have found geographic range size to be a strong predictor of extinction risk, with this relationship not explained by the positive correlation between abundance and geographic range size (e.g., Kiessling and Aberhan Reference Kiessling and Aberhan2007; Harnik et al. Reference Harnik, Simpson and Payne2012).
Univariate analyses cannot glean the relative importance of abundance and geographic range size as predictors of extinction risk, because these variables are measured in different units and are rarely evaluated together in a multivariate framework that adequately establishes their relative importance (Harnik et al. Reference Harnik, Simpson and Payne2012). Harnik et al. (Reference Harnik, Simpson and Payne2012) is an exception, in which the authors analyzed the association between extinction risk and traits related to taxon rarity (geographic range size, habitat breadth, and local abundance) using a multivariate framework and a global database of Phanerozoic marine fossil genera. The authors found that geographic range size was a strong predictor of extinction risk, whereas abundance had little effect. However, the authors were unable to calculate the relationship between geographic range size and extinction risk during the late Carboniferous and early Permian due to low numbers of generic extinctions during this interval. Consequently, here we assemble spatiotemporal fossil occurrence data to fill this knowledge gap and evaluate the relative strength of abundance and geographic range size as correlates of extinction during the late Paleozoic ice age (LPIA) using brachiopods from the midcontinent of the United States.
The late Paleozoic is interesting from a macroevolutionary perspective, because it was a time of “sluggish” evolution when global rates of origination and extinction were low (Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2005; Segessenman and Kammer Reference Segessenman and Kammer2018; Balseiro and Halpern Reference Balseiro and Halpern2019; Kolis and Lieberman Reference Kolis and Lieberman2019) in spite of dramatic, cyclical changes in climate and the environment (Parrish Reference Parrish1993; Olszewski and Patzkowsky Reference Olszewski and Patzkowsky2001a,Reference Olszewski and Patzkowskyb; Raymond and Metz Reference Raymond and Metz2004; Horton et al. Reference Horton, Poulsen, Montañez and DiMichele2012; Balseiro Reference Balseiro2016). The LPIA was the longest-lasting glacial period of the Phanerozoic, which lasted from 320 to 260 Ma. During the LPIA, Southern Hemisphere glacial–interglacial cycles were governed by Milankovitch orbital cycles (Raymond and Metz Reference Raymond and Metz2004). The resulting glacioeustatic sea-level changes produced cyclothems (Montañez and Poulsen Reference Montañez and Poulsen2013) characterized by vacillation of deep-marine, shallow-marine, and terrestrial environments (Heckel et al. Reference Heckel, Dennison and Ettensohn1994; Heckel Reference Heckel, Fielding, Frank and Isbell2008).
Late Paleozoic brachiopods within the midcontinent of the United States show little taxonomic turnover (Olszewski and Patzkowsky Reference Olszewski and Patzkowsky2001a,Reference Olszewski and Patzkowskyb), making them excellent targets for a study of macroevolutionary dynamics during this “quiescent” interval. Brachiopods are also extremely abundant and typically well preserved (Foote and Sepkoski Reference Foote and Sepkoski1999). Recent museum digitization efforts and ongoing initiatives, such as those funded by the Advancing the Digitization of Biodiversity Collections (ADBC) directorate of the U.S. National Science Foundation (NSF), allow for assembly of brachiopod datasets from which geographic range and abundance data can be derived (Page et al. Reference Page, MacFadden, Fortes, Soltis and Riccardi2015). Moreover, extinction and origination rates for brachiopods at a global scale match broadly those observed in a range of marine shelled invertebrate taxa (Powell Reference Powell2005).
The midcontinent of North America is an excellent target region to study the late Paleozoic, because it contains the best-studied record of the greenhouse/icehouse transitions during the Carboniferous and Permian (Heckel Reference Heckel1977, Reference Heckel1986; Heckel et al. Reference Heckel, Dennison and Ettensohn1994; Lupia and Armitage Reference Lupia and Armitage2013). The stratigraphy of the region is well constrained and, in many cases, correlated across state lines (Heckel Reference Heckel1986; Olszewski and Patzkowsky Reference Olszewski and Patzkowsky2003). The North American midcontinent was located near the equator for much of the Carboniferous, where changes associated with glaciation and global cooling were likely pronounced (Algeo and Heckel Reference Algeo and Heckel2008). Previous work has shown that cooling temperatures and increased seasonality, which marked the onset of glaciation during the Mississippian ca. 323 Ma (Parrish Reference Parrish1993), led to preferential loss of tropical brachiopod species with narrow latitudinal ranges (Powell Reference Powell2005). According to Powell (Reference Powell2005), brachiopod survivors had broader latitudinal distributions, potentially reflecting broader thermal tolerances, longer taxonomic durations, and larger populations. The absence of narrowly distributed and quickly evolving genera from the tropics led to a change in macroevolutionary dynamics, in which broadly adapted genera with long stratigraphic durations were common globally, not only in extratropical latitudes (Brezinski Reference Brezinski1988; Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2005; Segassenman and Kammer Reference Segessenman and Kammer2018).
Powell (Reference Powell2005) proposed that the low extinction and origination rates during the Late Mississippian to middle Permian were due, in part, to the higher proportion of brachiopod taxa with larger latitudinal ranges, potentially reflecting broader ecological tolerances. The same characteristics that typically buffer against extinction (large geographic range size; broad ecological tolerance; and abundant, stable populations) also inhibit isolation of populations that would lead to speciation, potentially impeding generation of new taxa, themselves usually characterized by small geographic distributions at higher risk of extinction (Vrba Reference Vrba1980; Jablonski Reference Jablonski1986; Stanley Reference Stanley1986; Foote Reference Foote2007; Antell et al. Reference Antell, Kiessling, Aberhan and Saupe2020). Thus, reduced variation and skew toward larger range sizes may have rendered this variable an ineffective predictor of extinction risk during the Late Mississippian to middle Permian. Under these circumstances, abundance, rather than geographic range size, may be a better predictor of extinction risk and macroevolutionary dynamics during the LPIA. Although geographic range size and abundance frequently covary, they do not exhibit a one-to-one relationship and have been shown to predict extinction risk to varying degrees (e.g., Stanley Reference Stanley1986; Kiessling and Aberhan Reference Kiessling and Aberhan2007; Payne et al. Reference Payne, Truebe, Nützel and Chang2011; Harnik et al. Reference Harnik, Simpson and Payne2012). We hypothesize that abundance is a better predictor of brachiopod extinction risk in the midcontinent and test this hypothesis using generalized linear models and 32,766 spatiotemporal occurrences over nine stages from the Chesterian to Leonardian (encompassing the LPIA).
Materials and Methods
Spatiotemporal Occurrence Data
Individual, global specimen occurrences for brachiopod species present in the Carboniferous–Permian of the midcontinent of the United States were obtained from multiple, spatially explicit databases, including the Division of Invertebrate Paleontology, Biodiversity Institute (KUMIP), the Yale University Peabody Museum of Natural History (YPM), and the Paleobiology Database (PBDB). A total of 32,766 specimen occurrence records were obtained from existing databases and digitization efforts, comprising 4998 records from the PBDB and 27,768 records from the KUMIP and YPM. Only midcontinental species were targeted. However, the entire geographic range of these midcontinental species was reconstructed globally to capture true extinctions rather than local extirpations. Specimen records from KUMIP and YPM were chosen because these institutions have large numbers of brachiopod specimens from the midcontinent, with a high degree of stratigraphic and geographic control. We retained only occurrences with geographic uncertainty radii under 50 km. Records that lacked species-level identifications (those including sp., cf., aff., or ? in the species designation) or those that did not have stratigraphic information to the level of formation were removed. Taxonomy was standardized and updated using Muir-Wood and Copper (Reference Muir-Wood and Cooper1960), Hoare (Reference Hoare1961), Moore (Reference Moore1964), Williams et al. (Reference Williams, Rowell, Muir-Wood, Pitrat, Schmidt, Stehli, Ager, Wright, Elliott, Amsden, Rudwick, Hatai, Biernat, McLaren, Boucot, Johnson, Stanton, Grant and Jope1965), and Carter and Carter (Reference Carter and Carter1970), which represent key references on brachiopods from this region and time interval.
To ensure that distributional data were derived from geologic units of similar ages, a stratigraphic database for the midcontinent was generated from an extensive survey of the primary literature (see Supplementary Material). Informal North American stages (e.g., Kinderhookian–Guadalupian; Heckel and Clayton Reference Heckel and Clayton2006; Menning et al. Reference Menning, Alekseev, Chuvashov, Davydov, Devuyst, Forke, Grunt, Hance, Heckel, Izokh, Jin, Jones, Kotlyar, Kozur, Nemyrovska, Schneider, Wang, Weddige, Weyer and Work2006) were followed to allow for the highest temporal resolution while maintaining the greatest sample size. Occurrences that could not be classified confidently to a North American stage were removed from analysis. Origination and extinction rates in each LPIA North American stage (Chesterian–Leonardian) were calculated using the modified gap-filler method of Alroy (Reference Alroy2015) in R v. 3.4.1 (R Core Team 2017) using the divDyn v. 0.8.0 package (Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019). The modified gap-filler method of Alroy (Reference Alroy2015) was used to minimize potential bias (e.g., those associated with the Signor-Lipps effect) while maintaining precision and accuracy. Rates from stages with <100 occurrences were not recorded but were used in calculations. To facilitate comparison with published species-level Paleozoic rates (e.g., Stigall Reference Stigall2010; Kolis and Lieberman Reference Kolis and Lieberman2019), we additionally calculated extinction and speciation rates as per-million-year rates using the equations of Foote (Reference Foote2000).
We assigned each species a chronostratigraphic duration based on the longest combination of (1) the chronostratigraphic duration reported in the PBDB; (2) the lithostratigraphic units the species occurred in according to Carter and Carter (Reference Carter and Carter1970); and (3) the lithostratigraphic units containing the species within museum databases. The chronostratigraphic ages of the lithostratigraphic units listed in these sources were determined using the compiled stratigraphic literature (see Supplementary Material) and the National Geologic Map Database (or Geolex; U.S. Geological Survey. n.d.). These chronostratigraphic durations were consulted when calculating per-million-year extinction and speciation rates, such that species with gaps in their records were never spuriously counted as extinctions. We tested the correlation between extinction and speciation rates using a generalized linear model (GLM) with a Gaussian distribution in R v. 3.4.1. To test for the effect of differential preservation on macroevolutionary rates, speciation and extinction rates were correlated with sampling intensity (estimated as the total number of brachiopod occurrences per stage; Powell Reference Powell2005) using a GLM with a Gaussian distribution in R v. 3.4.1. A positive correlation between evolutionary rates and sampling intensity could indicate that macroevolutionary rates were driven by differences in sampling intensity rather than biological patterns (Powell Reference Powell2005).
We ensured that paleogeographic analyses were conducted on spatially unique occurrences by culling each species to a single occurrence within an occupied grid cell of 0.025° × 0.025° resolution for each stage (equivalent to ~3 km at the equator). This procedure removed artificial inflation of spatially unique occurrences caused by differences in georeferencing protocols among institutions and individuals; for example, two museum workflows yielding different decimal degree latitude and longitude estimates for the same collection locality. We removed singletons from analyses (i.e., any species by time bin with n = 1 spatially unique occurrence) to limit the effect of poorly sampled taxa. After removal of singletons, there remained 29,720 individual occurrences (4915 of which were spatially unique) belonging to 164 species in 83 genera.
The resulting species by time bin datasets were imported into ArcGIS v. 10.5.1, and the present-day latitude/longitude coordinates were rotated to their paleo-position using the University of Texas Institute for Geophysics (UTIG) plate model in PaleoWeb v. 1.0 (Rothwell Group Inc). Paleo-coordinate reconstructions were performed based on the start of the appropriate stage in millions of years before present (Fig. 1). Due to the equatorial position of Laurasia in the late Paleozoic, we applied the South American Albers Equal Area Conic map projection to reconstructed paleolatitude and paleolongitude occurrences before range-size metrics were calculated. The South American projection was chosen given the preponderance of coordinates located in the region occupied by present-day Brazil once they were rotated to their paleolatitude and paleolongitude.
Paleogeographic reconstruction of spatially unique occurrence locations created using the PALEOMAP Paleo Atlas for GPlates v. 3 (Scotese Reference Scotese2016) for illustration purposes only. We used the University of Texas Institute for Geophysics (UTIG) plate model in PaleoWeb 1.0 for our analyses. Geographic range size was measured as a latitudinal range (yellow line with brackets) and as a convex hull (yellow polygon). Analyses used the median convex-hull area of all possible convex hulls created by jackknifing occurrences, which reduced the impact of geographic outliers on geographic range size. A, Desmoinesia muricatina (n = 118 spatially unique points) from the Desmoinesian. B, Neochonetes transversalis (n = 343 spatially unique points) from the Virgilian.
Geographic Range-Size Metrics
We quantified geographic range size and abundance for each unique species by stage. Geographic range was quantified using two different metrics: minimum convex hull and latitudinal range (for an example, see Fig. 1). A convex hull is a two-dimensional metric that estimates the amount of area occupied by a taxon and was calculated as the median of a series of convex hulls produced by exhaustive jackknifing of individual occurrence points (e.g., Stigall and Lieberman Reference Stigall and Lieberman2006; Hendricks et al. Reference Hendricks, Lieberman and Stigall2008; Myers and Lieberman Reference Myers and Lieberman2011; Darroch and Wagner Reference Darroch and Wagner2015; Saupe et al. Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015; Darroch et al. Reference Darroch, Casey, Antell, Sweeney and Saupe2020). This method has been shown to be especially efficacious for quantifying ranges of fossil taxa (Darroch and Saupe Reference Darroch and Saupe2018; Darroch et al. Reference Darroch, Casey, Antell, Sweeney and Saupe2020). If a species was characterized by only two spatially unique occurrences, a 10 km buffer was applied to each occurrence, and the area of the resulting line substituted for the convex hull (Hendricks et al. Reference Hendricks, Lieberman and Stigall2008; Myers and Lieberman Reference Myers and Lieberman2011). Each convex-hull geographic range estimate was log transformed for normality.
Latitudinal range (maximum observed paleolatitude minus minimum observed paleolatitude) is also used commonly to characterize geographic range sizes (Powell Reference Powell2005, Reference Powell2007; Foote and Miller Reference Foote and Miller2013; Finnegan et al. Reference Finnegan, Rasmussen and Harper2016; Balseiro and Halpern Reference Balseiro and Halpern2019; Darroch et al. Reference Darroch, Casey, Antell, Sweeney and Saupe2020). Unlike a convex hull, latitudinal range is a linear metric and may reflect breadth of thermal tolerance (Jackson Reference Jackson1974; Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2007; Sunday et al. Reference Sunday, Bates and Dulvy2012), with larger latitudinal ranges potentially indicating greater thermal tolerances. Latitudinal range estimates were square-root transformed for normality.
Abundance was quantified as the number of occurrences for a species within a stage. Calculating population size is not easy, even for modern organisms (He and Gaston Reference He and Gaston2000). However, occurrence data have been shown to provide adequate estimates of average abundance for both fossil and modern organisms (Buzas et al. Reference Buzas, Koch, Culver and Sohl1982; Kunin Reference Kunin1998; Alroy Reference Alroy2000; He and Gaston Reference He and Gaston2003). Abundance estimates were square-root transformed for normality.
All three variables (convex hull, latitudinal range, and abundance) were calculated in R v. 3.4.1 (R Core Team 2017) using the packages sp v. 1.3-2 (Pebesma and Bivand Reference Pebesma and Bivand2005; Bivand et al. Reference Bivand, Pebesma and Gomez-Rubio2013) and PBSmapping v. 2.72.1 (Schnute et al. Reference Schnute, Boers, Haigh, Grandin, Johnson, Wessel and Antonio2013).
Mixed-Effects Models
We explored the nature of the relationship between species’ extinction and abundance/range size metrics using generalized linear mixed-effects models (LMM) using both transformed (log or square root) unstandardized variables and standardized variables. For the latter, variables were standardized at stage level by dividing the value for each species by the maximum value of that variable in the stage (Foote et al. Reference Foote, Crampton, Beu and Cooper2008). For example, the log geographic range size for Wellerella truncata during the Virgilian was divided by the maximum log-transformed geographic range size of any species from the Virgilian. Therefore, standardization scaled all three variables to vary between 0 and 1 in each stage. The results for the standardized variable analysis are included in the Supplementary Material (Supplementary Tables S3 and S4). For both unstandardized and standardized analyses, stages were limited to those associated with the LPIA (e.g., Chesterian–Leonardian). The time interval and the Linnaean family rank associated with each record were included as random effects in the LMM. We included these random effects because species from the same time bin and family were expected to pseudoreplicate one another due to shared environmental conditions/spatial sampling structure and evolutionary history, respectively. These random effects obviated the need to correct for multiple testing, because a single coefficient was calculated for each fixed effect over the entire time series.
We estimated the influence of the following three fixed effects in every combination: latitudinal range, abundance, and convex-hull area. This resulted in a total of eight possible models, including a model with no fixed effects added (i.e., a horizontal line). We implemented mixed-effects models with functions in the package lme4 v. 1.1-21 (Venables and Ripley Reference Venables, Ripley, Venables and Ripley2002; Bates et al. Reference Bates, Mächler, Bolker and Walker2015).
Akaike's information criterion (AIC) was estimated instead of corrected AIC, given the large number of species in the dataset, to determine the subset of models that received nontrivial weight (ΔAIC < 10; for derivation and discussion of weights, see Burnham and Anderson Reference Burnham and Anderson2002). Models are listed in order of relative performance. For fixed effects of models, we recorded coefficient estimates on the logit scale and the 95% confidence interval bounds from the likelihood profile. Models were checked for overdispersion and linearity between continuous predictor variables and the logit of the outcome.
Temporal Analysis of Trends in Range Size and Abundance
Foote (Reference Foote2007) suggested that genera tend to exhibit declining geographic range sizes preceding their extinction. Consequently, we evaluated temporal trends in geographic range size, latitudinal range, and abundance throughout a species’ lifetime, including occurrence records from the Carboniferous to Permian (Kinderhookian–Lopingian). Unlike the generalized linear mixed effects framework, which accounted for shared environmental conditions and spatial sampling structure by specifying stage as a random effect, this analysis compared geographic range size and abundance estimates for a single species between stages characterized by differing area of rock outcrop. To account for effects of rock availability on our three metrics (geographic range size, latitudinal range, and abundance), variables were standardized at the stage level as described earlier (Foote et al. Reference Foote, Crampton, Beu and Cooper2008). Temporal trends were evaluated using both standardized and unstandardized variables. Results for unstandardized variables are reported in the Supplementary Material (Supplementary Table S5).
Species present in three or more consecutive stages were determined to decline before extinction if they conformed to the following conditions: (1) the terminal value was lower than the value in the immediately preceding stage; and (2) the terminal value was lower than the mean for the species (Fig. 2). The number of species in decline was tabulated for each metric (convex hull, latitudinal range, and abundance). A one-tailed binomial test on the proportion of species in decline versus the number of stable and increasing species determined whether this proportion was significantly higher than 50%. That is, the test examined whether more species experienced declines in abundance or geographic range size before extinction than species that increased or remained unchanged; effect size was evaluated using Cohen's g.
Standardized geographic range size and abundance through time for the brachiopod species Cancrinella boonensis. The dashed horizontal line indicates the mean value for the species over this interval. A, Standardized geographic range size, measured as convex-hull area (km2) over maximum convex-hull area (km2) for the stage. Cancrinella boonensis specimen IP.008072 Yale Peabody Museum of Natural History; photo by J. Utrup, 2011. B, Standardized geographic range size, measured as species latitudinal range over maximum latitudinal range for the stage. C, Standardized abundance, measured as species abundance over maximum abundance for the stage. Abundance is the only metric that meets the criteria for decline before extinction—i.e., the terminal value is lower than the preceding value and below the species mean. Stage abbreviations: O, Osagean; Me, Meramecian; C, Chesterian; Mo, Morrowan; A, Atokan; D, Desmoinesian; M, Missourian; V, Virgilian; Wol, Wolfcampian.
Results
Speciation and Extinction Rates
Rates of speciation and extinction were low for most of the stages associated with the LPIA (Supplementary Table S1). For the Chesterian to Missourian stages, proportional extinction ranged between 0.014 and 0.154. Highest proportional extinction was in the Wolfcampian at 0.884, preceded by the second-highest proportional extinction of 0.363 in the Virgilian. Proportional extinction was lowest for the Missourian at 0.014. The median proportional extinction rate for the LPIA stages was 0.091. Extinction rates measured as per-million-year rates were lowest in the Morrowan (0.004) and highest in the Virgilian (0.047). As with proportional extinction, the highest per-million-year extinction rates were found in the Virgilian and Wolfcampian, but their relative ranks were reversed. For context, background extinction rates derived from species-level brachiopod data for the Middle Devonian range from 0.2 to 0.6 per million years (Stigall Reference Stigall2010), while late Paleozoic extinction rates derived from species-level cephalopod data range from 0.07 to 0.34 per million years (Kolis and Lieberman Reference Kolis and Lieberman2019). Thus, per-million-year extinction rates for all stages are below the range of Middle Devonian per-million-year background rates and comparable to or lower than other species-level late Paleozoic studies. Note that this dataset could not capture the Serupkhovian mass extinction (Chesterian), because we focused on LPIA taxa and omitted any taxa entirely absent from the Pennsylvanian or Permian.
Proportional speciation ranged from 0.154 to 0.486 and was always higher than proportional extinction, except during the Wolfcampian. The highest proportional speciation, 0.486, occurred during the Virgilian. Median proportional speciation rate for the LPIA was 0.173. Proportional speciation and extinction rates are not significantly correlated (GLM t = 0.66, p = 0.53). Per-million-year speciation rates ranged between 0.008 in the Wolfcampian and 0.054 in the Morrowan and are low relative to per-million-year extinction rates derived from species-level data on late Paleozoic cephalopods (0.03 to 0.41; Kolis and Lieberman Reference Kolis and Lieberman2019). Per-million-year speciation and extinction rates are not significantly correlated (GLM t = 0.11, p = 0.92).
To consider whether variations in extinction and speciation rate might be related to changes in sampling intensity, we examined the relationship between sampling intensity and speciation and extinction rates using GLMs (Powell Reference Powell2005). No significant correlations were found between sampling intensity (estimated as the total number of brachiopod occurrences in each stage) and either extinction rate (GLM t = 0.41, p = 0.69) or speciation rate (GLM t = 1.23, p = 0.25) (Supplementary Table S2), indicating that patterns in origination and extinction likely reflect real biological fluctuations, not simply changes in preservation or sample availability. However, it is worth noting that edge effects are likely to inflate origination and extinction rates at both the start and end of the examined time interval (Foote Reference Foote2000); extending the study interval to account for edge effects would likely reduce the rates for both the Chesterian and Wolfcampian reported herein.
Mixed-Effects Models
The three assessed predictors of extinction risk (convex hull, latitudinal range, and abundance) were not highly collinear, based on the variance inflation factor (VIF). The VIF represents the proportion of variance in one predictor explained by all the other predictors (<2.5, with 1 being no collinearity; Zorro et al. Reference Zuur, Ieno and Elphick2010) and was 2.74, 2.18, and 1.41 for convex hull, latitudinal range, and abundance, respectively.
Mixed-effects models indicate that abundance is the strongest predictor of extinction risk. The three highest ranking mixed-effects models all included abundance as a predictor (Table 1). Moreover, abundance was the only predictor in the top three models to have confidence interval estimates exclusive of zero (Table 1) and was the sole predictor included in the best-supported model (Table 2). By contrast, range size–only models received little support (Table 1). Both latitudinal range and convex-hull area performed best when included as predictors with abundance (Tables 1 and 2, models ranked 2 and 3), but confidence intervals for these predictors were almost always inclusive of zero (Table 1). Of the best-supported models, only one model included measures of geographic range with confidence intervals that were exclusive of zero (model 5; convex hull as the only predictor). Results remain virtually unchanged when conducted using standardized variables (Supplementary Tables S3 and S4).
Coefficient estimates and confidence intervals (CI) for mixed-effects models using unstandardized variables. Coefficient estimates (logit scale) and confidence intervals for each predictor in the eight models, with extinction status as response. Models were built on 377 range size/abundance records for 164 species from nine Phanerozoic stages. Abundance and latitudinal range were square-root transformed; convex hull was log transformed for normality. The cumulative weight of the best-supported models was 99.7%. Confidence intervals, however, are wide, because there were relatively few extinctions even when pooled across intervals (88 extinctions compared with 289 survivals across the nine stages). Confidence intervals inclusive of zero are in bold.
Relative performance of all mixed-effects models using unstandardized variables. Models predicted species extinction status based on measures of abundance and range size (latitudinal range and convex hull). Weights are the ratio of ΔAIC from a given model to the sum of ΔAIC values across all candidate models and are interpreted as the probability that a given model is the “best” (minimizes the Kullback–Leibler discrepancy) of the candidate models, given the data.
Temporal Analysis of Trends in Range Size and Abundance
No statistical support was found for declining geographic range size in species before their extinction. In particular, the proportion of species with range size declines before extinction was no higher than the proportion of species without such declines, regardless of whether latitudinal range or convex hulls were used to quantify geographic range (Table 3). By contrast, species were found to decline in abundance in the interval leading to extinction significantly more than half the time (proportion of species displaying declining abundance = 62%, p = 0.05; Table 3), although Cohen's g suggests the effect size of this difference can be classified as small (g = 0.12). Results remain unchanged when conducted using unstandardized variables (Supplementary Table S5).
One-tailed binomial test results for the probability that species decline in range size and abundance before extinction using standardized variables. Species present in three or more consecutive stages were determined to decline before extinction if they conformed to the following conditions: (1) the terminal value was lower than the value in the immediately preceding stage, and (2) the terminal value was lower than the mean for the species. Values include abundance and geographic range size measured as either convex hull or latitudinal range.
Discussion
Late Paleozoic brachiopods from the North American midcontinent provide a diverse and abundant record of marine life during a distinctive time in Earth history associated with profound climatic oscillations. This time interval is especially noteworthy for the low extinction and speciation rates displayed in marine invertebrates (Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2005; Segessenman and Kammer Reference Segessenman and Kammer2018; Balseiro and Halpern Reference Balseiro and Halpern2019; Kolis and Lieberman Reference Kolis and Lieberman2019). Our analysis of macroevolutionary rates (Supplementary Table S1) is consistent with results from previous studies indicating this was a time of sluggish macroevolution (Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2005; Segessenman and Kammer Reference Segessenman and Kammer2018; Balseiro and Halpern Reference Balseiro and Halpern2019; Kolis and Lieberman Reference Kolis and Lieberman2019). In particular, per-million-year extinction and speciation rates for LPIA intervals are low (Supplementary Table S1), consistent with previously published work (Stanley and Powell Reference Stanley and Powell2003; Powell Reference Powell2005). Per-million-year speciation rates are low throughout the interval, comparable to per-million-year speciation rates from cephalopods during the same interval (Kolis and Lieberman Reference Kolis and Lieberman2019). Similarly, the slight increase in extinction during the Late Pennsylvanian (Virgilian) and early Permian (Wolfcampian) is consistent with patterns of brachiopod extinction found previously by Olzewski and Patzkowsky (Reference Olszewski and Patzkowsky2001a), who attributed the increase to gradual sea-level fall and a shift toward more arid conditions evidenced by the increased appearance of evaporites. The weak to absent correlation between extinction and speciation rates is somewhat unusual (Stanley et al. 1990) and adds further credence to the notion this time period is exceptional with respect to macroevolutionary dynamics.
Species that went extinct during the LPIA had smaller geographic range sizes and lower abundances (Fig. 3) than survivors. Although geographic range size has received robust support as a determinant of extinction risk during many, if not most, geologic intervals (Jablonski Reference Jablonski1986; McKinney Reference McKinney1997; Payne and Finnegan Reference Payne and Finnegan2007; Harnik et al. Reference Harnik, Simpson and Payne2012; Saupe et al. Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015), we found that abundance, rather than geographic range size, was a strong predictor of extinction risk for midcontinental brachiopod species during the late Paleozoic of the North American midcontinent. These patterns receive further support from our analyses, which suggest that species’ geographic range sizes did not decline before extinction but their abundance did (Table 3, Fig. 2). Our results reinforce previous work on LPIA invertebrates from the midcontinent. For example, Kolis and Lieberman (Reference Kolis and Lieberman2019) found that geographic range sizes for cephalopod species did not correlate with extinction rates. Using a multiple linear regression, Powell (Reference Powell2007) demonstrated that abundance and geographic range size contributed equally to genus duration in midcontinental brachiopods. Although Powell's (Reference Powell2007) results for the importance of abundance are similar to those reported here, our results may differ for geographic range because of differences in statistical methodology, taxonomic level, and/or the use of summary metrics across the entire duration of the genus to characterize correlates of duration, rather than extinction risk.
Geographic range size and abundance estimates for species that survive and go extinct at stage boundaries (n = 88 extinct, n = 289 survive). A total of 164 unique species were analyzed across nine stages (Chesterian–Leonardian). Analyses used the median convex-hull area (km2) of all possible convex hulls created by jackknifing occurrences and was log transformed. Latitudinal range and abundance were square-root transformed from original count and degree latitude measurements, respectively. Species trait data were pooled, such that multiple measurements for species were recorded if they occurred in more than one stage.
The abundance–survivorship relationship reported herein could be confounded by ecological or physiological differences among species not controlled for in this analysis. Other potential correlates of extinction risk could include differences in basal metabolic rate (Strotz et al. Reference Strotz, Saupe, Kimmig and Lieberman2018), trophic ecology (Norris Reference Norris1992), factors affecting population structure (Norris Reference Norris1992), variations in dispersal ability (Powell Reference Powell2007; Birand et al. Reference Birand, Vose and Gavrilets2012), or local abundance (i.e., average relative abundance across samples or localities; Stanley Reference Stanley1986). Many of these factors are difficult to ascertain for brachiopods. For example, brachiopods utilize a range of larval development strategies (Thayer Reference Thayer1981; James et al. Reference James, Ansell, Collins, Curry, Peck and Rhodes1992; Peck and Robinson Reference Peck and Robinson1994) that cannot be inferred directly for extinct species in most cases (Valentine and Jablonski Reference Valentine and Jablonski1983; Rowell Reference Rowell1986), which makes dispersal ability difficult to estimate. It is also possible that implicit or explicit biases in the way fossils are collected and/or the way museum collections and databases are built fundamentally impugn the use of species abundance in these types of analyses.
Stanley and Powell (Reference Stanley and Powell2003) and Powell (Reference Powell2005) suggested the low rates of speciation and extinction in the late Paleozoic represented a “new macroevolutionary state” relative to the rest of the Paleozoic. This dampened species turnover has been attributed to the loss of extinction-prone taxa earlier in the Paleozoic, which resulted in a higher proportion of widely distributed and long-duration brachiopod genera in the paleo-tropics, not only at higher latitudes (Brett and Baird Reference Brett, Baird, Erwin and Anstey1995; Powell Reference Powell2005). The physiographic conditions of the late Paleozoic sea in the North American midcontinent—a gently-sloping ramp with few geographic barriers to latitudinal movement—likely aided dispersal by marine organisms and decreased opportunities for both extinction and speciation (Heckel Reference Heckel1986; Olszewski and Patzkowsky Reference Olszewski and Patzkowsky2003). Thus, enhanced abiotic potential for broad ranges in the late Paleozoic sea of the North American midcontinent, combined with the biotic characteristics of most brachiopods present, may have yielded a fauna that lacked a prominent geographic component of extinction susceptibility. Instead, the abundance–extinction relationship observed during the LPIA could reflect the vulnerability of populations with low abundance to the proximate causes of extinction, namely demographic stochasticity, environmental stochasticity, and genetic deterioration (Goodman Reference Goodman1987; Lande Reference Lande1993; Pimm et al. Reference Pimm, Diamond, Reed, Russell and Verner1993; Wissel et al. Reference Wissel, Stephan, Zaschke and Remmert1994; Henle et al. Reference Henle, Davies, Kleyer, Margules and Settele2004a,Reference Henle, Sarre and Wiegandb). Abundance may play an important role in regulating extinction risk during biodiversity crises or extinctions that, like the late Paleozoic, are characterized by low rates of origination (e.g., the Late Devonian biodiversity crisis or the end-Triassic mass extinction; Bambach et al. Reference Bambach, Knoll and Wang2004). Such claims, however, require additional research.
Extinction and speciation rates are frequently correlated with geographic range size (Vrba Reference Vrba1980; Jablonski Reference Jablonski1986; Stanley Reference Stanley1986, Reference Stanley, Ross and Allmon1990b; Eldredge Reference Eldredge1989; Jablonski and Roy Reference Jablonski and Roy2003), and small-ranged taxa are often culled preferentially during extinction events of different scales (Jablonski Reference Jablonski1986; Payne and Finnegan Reference Payne and Finnegan2007; Powell Reference Powell2007; Clapham and Payne Reference Clapham and Payne2011). Therefore, background and mass extinctions could lead to preferential loss of volatile taxa with intrinsically high rates of origination and extinction over time (Gilinksy Reference Gilinsky1994; Lieberman and Melott Reference Lieberman and Melott2013). During “normal” evolutionary times, the extinction-prone, narrowly distributed taxa removed by extinction would be replaced after several million years, as new speciation events tend to generate small-ranged taxa (Foote Reference Foote2007; Liow and Stenseth Reference Liow and Stenseth2007; Antell et al. Reference Antell, Kiessling, Aberhan and Saupe2020). The eventual reappearance of these small-ranged taxa would have therefore reinstated geographic range size as an important determinant of extinction risk. However, the relative dearth of narrowly distributed taxa during the late Paleozoic, combined with the distinctive physiographic nature of the region studied, could have conspired to suppress both extinction and speciation rates and hinder emergence from this distinctive “macroevolutionary state”. This proposed mechanism shares many similarities with Stanley's (Reference Stanley1990a) explanation for the apparent periodicity of mass extinctions, in which the suggested recurrence of extinction events is limited by the length of time necessary to recover extinction-prone or vulnerable taxa.
The abundance–extinction risk relationship uncovered here may be unique to our study system of the North American midcontinent during the LPIA, which focused on brachiopods, tropical latitudes, and epicontinental seaways. For example, study of brachiopod dynamics from midlatitudes and open ocean–facing systems found high regional turnover in response to the LPIA (Balseiro Reference Balseiro2016) and regional extirpation selectivity associated with both genus range size and body size (Balseiro and Halpern Reference Balseiro and Halpern2019). These differing dynamics suggest that brachiopod faunas outside tropical, epicontinental seaways may have responded to the LPIA differently than those in the U.S. midcontinent. Indeed, macroevolutionary dynamics of epicontinental seas may differ from those operating along ocean-facing shelves (Miller and Foote Reference Miller and Foote2009), such that extinction risk patterns could differ in open-shelf communities of the same age or in post-Paleozoic periods when epicontinental seas decrease in prevalence (Peters Reference Peters2007). However, if the abundance–extinction risk relationship found herein is extendable to other regions and taxa during the late Paleozoic, this could add further nuance to Jablonksi's (Reference Jablonski1986) “alternation of macroevolutionary regimes,” which hypothesizes that mass extinctions are unique and not simply intensifications of background extinction dynamics. The possibility exists that macroevolutionary processes shift between intervals of background extinction and of mass extinction, but also between times of background extinction and intervals of sluggish turnover.
Acknowledgments
We thank M. Powell and M. Foote for thoughtful and constructive reviews. We thank museum staff for help with brachiopod collections, including S. Butts (YPM), J. Utrup (YPM), U. Farrell (KUMIP at the time, now N. Lopez Carranza). We thank T. West for assistance with ArcGIS and PaleoWeb. Financial support for this project was provided by NSF DEB 1256993, EF 1206757, and DBI 1602067 to B.S.L., and Leverhulme grant DGR01020 to E.E.S. This material is based upon work supported while working at the National Science Foundation.
Open access
