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 spatio-temporal 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 prior to 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.

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 4,998 records from the PBDB and 27,768 records from the KUMIP and YPM. Only midcontinent species were targeted. However, the entire geographic range of these midcontinent 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 extensive survey of the primary literature (see Suppl. 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 whilst 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’ v0.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) whilst maintaining precision and accuracy. Rates from stages with < 100 occurrences were not recorded but used in calculations. To facilitate comparison with published species-level Paleozoic rates (e.g., Stigall 2010; Kolis 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 the above sources were determined using the compiled stratigraphic literature (see Suppl. Material) and the National Geologic Map Database (or Geolex https://ngmdb.usgs.gov/Geolex/search). 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° x 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 (4,915 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.

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; Powell 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).