The effects of geographic range size and abundance on extinction during a time of “sluggish”’ evolution

Abstract

Geographic range size and abundance are important determinants of extinction risk in fossil and extant taxa. However, the relationship between these variables and extinction risk has not been tested extensively during evolutionarily “quiescent” times of low extinction and speciation in the fossil record. Here we examine the influence of geographic range size and abundance on extinction risk during the late Paleozoic (Mississippian–Permian), a time of “sluggish” evolution when global rates of origination and extinction were roughly half those of other Paleozoic intervals. Analyses used spatiotemporal occurrences for 164 brachiopod species from the North American midcontinent. We found abundance to be a better predictor of extinction risk than measures of geographic range size. Moreover, species exhibited reductions in abundance before their extinction but did not display contractions in geographic range size. The weak relationship between geographic range size and extinction in this time and place may reflect the relative preponderance of larger-ranged taxa combined with the physiographic conditions of the region that allowed for easy habitat tracking that dampened both extinction and speciation. These conditions led to a prolonged period (19–25 Myr) during which standard macroevolutionary rules did not apply.

Introduction

Determining the correlates of extinction is crucial to understanding macroevolutionary processes operating on geologic timescales (Jablonski 1986; McKinney 1997; Kiessling and Aberhan 2007; Payne and Finnegan 2007; Meseguer et al. 2015; Saupe et al. 2015) and for identifying taxa at potential risk of extinction today (Lee and Jetz 2011; Finnegan et al. 2015; Kiessling et al. 2019; Smits and Finnegan 2019). 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 1994; Gaston et al. 1997; Holt et al. 1997; Harnik 2011) but are not always equally important determinants of extinction risk (e.g., Kiessling and Aberhan 2007; Payne et al. 2011; Harnik et al. 2012). For example, examining Neogene bivalves in the Pacific, Stanley (1986) found that abundance, and not geographic range size, was a strong predictor of extinction selectivity. Similarly, Payne et al. (2011) 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 2007; Harnik et al. 2012).

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. 2012). Harnik et al. (2012) 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 2003; Powell 2005; Segessenman and Kammer 2018; Balseiro and Halpern 2019; Kolis and Lieberman 2019) in spite of dramatic, cyclical changes in climate and the environment (Parrish 1993; Olszewski and Patzkowsky 2001a,b; Raymond and Metz 2004; Horton et al. 2012; Balseiro 2016). 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 2004). The resulting glacioeustatic sea-level changes produced cyclothems (Montañez and Poulsen 2013) characterized by vacillation of deep-marine, shallow-marine, and terrestrial environments (Heckel et al. 1994; Heckel 2008).

Late Paleozoic brachiopods within the midcontinent of the United States show little taxonomic turnover (Olszewski and Patzkowsky 2001a,b), 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 1999). 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. 2015). Moreover, extinction and origination rates for brachiopods at a global scale match broadly those observed in a range of marine shelled invertebrate taxa (Powell 2005).

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 1977, 1986; Heckel et al. 1994; Lupia and Armitage 2013). The stratigraphy of the region is well constrained and, in many cases, correlated across state lines (Heckel 1986; Olszewski and Patzkowsky 2003). 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 2008). Previous work has shown that cooling temperatures and increased seasonality, which marked the onset of glaciation during the Mississippian ca. 323 Ma (Parrish 1993), led to preferential loss of tropical brachiopod species with narrow latitudinal ranges (Powell 2005). According to Powell (2005), 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 1988; Stanley and Powell 2003; Powell 2005; Segassenman and Kammer 2018).

Powell (2005) 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 1980; Jablonski 1986; Stanley 1986; Foote 2007; Antell et al. 2020). 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 1986; Kiessling and Aberhan 2007; Payne et al. 2011; Harnik et al. 2012). 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 (1960), Hoare (1961), Moore (1964), Williams et al. (1965), and Carter and Carter (1970), 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 2006; Menning et al. 2006) 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 (2015) in R v. 3.4.1 (R Core Team 2017) using the divDyn v. 0.8.0 package (Kocsis et al. 2019). The modified gap-filler method of Alroy (2015) 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 2010; Kolis and Lieberman 2019), we additionally calculated extinction and speciation rates as per-million-year rates using the equations of Foote (2000).

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 (1970); 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 2005) 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 2005).

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.

Figure 1. Paleogeographic reconstruction of spatially unique occurrence locations created using the PALEOMAP Paleo Atlas for GPlates v. 3 (Scotese 2016) 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 2006; Hendricks et al. 2008; Myers and Lieberman 2011; Darroch and Wagner 2015; Saupe et al. 2015; Darroch et al. 2020). This method has been shown to be especially efficacious for quantifying ranges of fossil taxa (Darroch and Saupe 2018; Darroch et al. 2020). 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. 2008; Myers and Lieberman 2011). 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 2005, 2007; Foote and Miller 2013; Finnegan et al. 2016; Balseiro and Halpern 2019; Darroch et al. 2020). Unlike a convex hull, latitudinal range is a linear metric and may reflect breadth of thermal tolerance (Jackson 1974; Stanley and Powell 2003; Powell 2007; Sunday et al. 2012), 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 2000). However, occurrence data have been shown to provide adequate estimates of average abundance for both fossil and modern organisms (Buzas et al. 1982; Kunin 1998; Alroy 2000; He and Gaston 2003). 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 2005; Bivand et al. 2013) and PBSmapping v. 2.72.1 (Schnute et al. 2013).

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. 2008). 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 2002; Bates et al. 2015).

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 2002). 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 (2007) 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. 2008). 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.

Figure 2. 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 (km²) over maximum convex-hull area (km²) 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 2010), while late Paleozoic extinction rates derived from species-level cephalopod data range from 0.07 to 0.34 per million years (Kolis and Lieberman 2019). 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 2019). 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 2005). 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 2000); 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. 2010) 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).

Table 1. 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.

Table 2. 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).

Table 3. 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 2003; Powell 2005; Segessenman and Kammer 2018; Balseiro and Halpern 2019; Kolis and Lieberman 2019). 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 2003; Powell 2005; Segessenman and Kammer 2018; Balseiro and Halpern 2019; Kolis and Lieberman 2019). 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 2003; Powell 2005). 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 2019). 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 (2001a), 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 1986; McKinney 1997; Payne and Finnegan 2007; Harnik et al. 2012; Saupe et al. 2015), 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 (2019) found that geographic range sizes for cephalopod species did not correlate with extinction rates. Using a multiple linear regression, Powell (2007) demonstrated that abundance and geographic range size contributed equally to genus duration in midcontinental brachiopods. Although Powell's (2007) 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.

Figure 3. 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 (km²) 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. 2018), trophic ecology (Norris 1992), factors affecting population structure (Norris 1992), variations in dispersal ability (Powell 2007; Birand et al. 2012), or local abundance (i.e., average relative abundance across samples or localities; Stanley 1986). Many of these factors are difficult to ascertain for brachiopods. For example, brachiopods utilize a range of larval development strategies (Thayer 1981; James et al. 1992; Peck and Robinson 1994) that cannot be inferred directly for extinct species in most cases (Valentine and Jablonski 1983; Rowell 1986), 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 (2003) and Powell (2005) 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 1995; Powell 2005). 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 1986; Olszewski and Patzkowsky 2003). 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 1987; Lande 1993; Pimm et al. 1993; Wissel et al. 1994; Henle et al. 2004a,b). 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. 2004). Such claims, however, require additional research.

Extinction and speciation rates are frequently correlated with geographic range size (Vrba 1980; Jablonski 1986; Stanley 1986, 1990b; Eldredge 1989; Jablonski and Roy 2003), and small-ranged taxa are often culled preferentially during extinction events of different scales (Jablonski 1986; Payne and Finnegan 2007; Powell 2007; Clapham and Payne 2011). Therefore, background and mass extinctions could lead to preferential loss of volatile taxa with intrinsically high rates of origination and extinction over time (Gilinksy 1994; Lieberman and Melott 2013). 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 2007; Liow and Stenseth 2007; Antell et al. 2020). 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 (1990a) 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 2016) and regional extirpation selectivity associated with both genus range size and body size (Balseiro and Halpern 2019). 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 2009), 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 2007). 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 (1986) “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.