Introduction
Pleistocene invertebrate fossils from southern California's uplifted coastal marine terraces have been studied for more than 150 years and provide an excellent opportunity to analyze the responses of coastal ecosystems to major episodes of climate change (e.g., Conrad Reference Conrad1855; Valentine Reference Valentine1961; Lindberg and Lipps Reference Lindberg, Lipps, Jablonski, Erwin and Lipps1996; Muhs and Groves Reference Muhs and Groves2018). Up to 17% of species occurring in this exceptionally complete fossil record (Valentine Reference Valentine1989) are no longer present in the California region (Pacific coast from 34.4°N to 27.8°N; Valentine and Jablonski Reference Valentine, Jablonski, Ricklefs and Schluter1993). These so-called “extralimital” species—that is, species that formerly occupied regions outside their present geographic limits—record a dynamic history of species range shifts (Valentine Reference Valentine and Sears1959, Reference Valentine1961; Valentine and Jablonski Reference Valentine, Jablonski, Ricklefs and Schluter1993; Jablonski and Sepkoski Reference Jablonski and Sepkoski1996; Roy et al. Reference Roy, Valentine, Jablonski and Kidwell1996; Hall Reference Hall2002) and have been used to biostatratigraphically correlate deposits (e.g., Wehmiller et al. Reference Wehmiller, Lajoie, Kvenvolden, Peterson, Belknap and Kennedy1977; Kennedy et al. Reference Kennedy, Lajoie and Wehmiller1982), infer paleoenvironments (e.g., Woodring et al. Reference Woodring, Bramlette and Kew1946; Valentine Reference Valentine1961; Addicott Reference Addicott1966; Kennedy et al. Reference Kennedy, Lajoie and Wehmiller1982; Powell Reference Powell2001; Powell et al. Reference Powell, Grant and Conkling2004, Reference Powell, II, Stanton, Vendrasco and Liff-Grief2009), and examine the ecological and evolutionary underpinnings of range shifts (e.g., Roy et al. Reference Roy, Jablonski and Valentine1995, Reference Roy, Jablonski and Valentine2001; Hellburg et al. Reference Hellburg, Deborah and Roy2001).
Extralimital occurrences appear to track climatic cycles, with extralimitals emigrating from the north during colder climate states and from the south during warmer climate states (Valentine Reference Valentine1961; Valentine and Jablonski Reference Valentine, Jablonski, Ricklefs and Schluter1993; Roy et al. Reference Roy, Valentine, Jablonski and Kidwell1996; Powell et al. Reference Powell, Lajoie and Ponti2000; Muhs et al. Reference Muhs, Simmons and Steinke2002, Reference Muhs, Simmons, Schumann, Groves, Mitrovica and Laurel2012, Reference Muhs, Groves and Schumann2014; Muhs and Groves Reference Muhs and Groves2018). Extralimital species from warmer-than-present climate states thus have potential to provide information about which species may be expected to recolonize once regional climates sufficiently change to allow them to track their environmental tolerances back into their paleo-ranges. In this sense, Pleistocene extralimitals can provide pre-anthropogenic baseline expectations for species movements in response to regional climate change.
Although this record has exceptional promise for unraveling the drivers of species range shifts, due to the historical difficulty of independently dating marine terrace deposits (see discussions in Muhs et al. Reference Muhs, Simmons and Steinke2002, Reference Muhs, Simmons, Schumann, Groves, Mitrovica and Laurel2012, Reference Muhs, Groves and Schumann2014; Muhs and Groves Reference Muhs and Groves2018) most previous investigations of range dynamics have been conducted at the composite mid–late Pleistocene level (e.g., Lindberg and Lipps Reference Lindberg, Lipps, Jablonski, Erwin and Lipps1996; Roy et al. Reference Roy, Jablonski and Valentine1995, Reference Roy, Valentine, Jablonski and Kidwell1996, Reference Roy, Jablonski and Valentine2001) and have not been able to resolve range dynamics during discrete climate states. In addition, there are still open questions about which factors are most important in determining species range shifts. The observation that southern extralimital species tend to have broader present-day latitudinal ranges than non-extralimital southern species (Roy et al. Reference Roy, Jablonski and Valentine1995) suggests that environmental niche breadth was an important factor controlling past range shifts, but because latitudinal range is correlated with a variety of environmental factors, it is not clear which specific factors had the greatest influence.
Here we take advantage of advances in late Pleistocene geochronology (Muhs et al. Reference Muhs, Simmons and Steinke2002, Reference Muhs, Simmons, Kennedy, Ludwig and Groves2006, Reference Muhs, Simmons, Schumann, Groves, Mitrovica and Laurel2012, Reference Muhs, Groves and Schumann2014; Muhs and Groves Reference Muhs and Groves2018) and satellite measurements of recent sea-surface temperatures and productivity to examine the modern-day characteristics of extralimital species that expanded their ranges into the California region (Fig. 1) during the warmest climate state of the last 300,000 years: the peak of the last interglacial, Marine Isotope Substage (MIS) 5e (~120,000 years ago; Shackleton and Opdike Reference Shackleton and Opdike1973). We focus on bivalve species, because their biogeographic distributions are exceptionally well studied in the eastern Pacific (e.g., Coan et al. Reference Coan, Scott and Bernard2000; Coan and Valentich-Scott Reference Coan and Valentich-Scott2012). Our specific goals are: (1) to provide an exhaustive compilation of MIS 5e extralimital species from paleo-embayment and paleo–open coast habitats—which differ markedly in their temperature, salinity, nutrient content, substrate, and wave-energy conditions (Ricketts et al. Reference Ricketts, Calvin and Hedgpeth1985; Coan et al. Reference Coan, Scott and Bernard2000; Coan and Valentich-Scott Reference Coan and Valentich-Scott2012)—in southern California; and (2) to determine which factors were most important in controlling species range shifts by comparing the distributional (environmental and geographic) characteristics of MIS 5e extralimital species with those of the broader Northern Panamic fauna (Fig. 1). By investigating the characteristics of extralimital species from the last full interglacial, we further seek to establish a pre-anthropogenic baseline for helping to understand likely range shifts in response to contemporary warming.
California and Northern Panamic region map with Marine Isotope Substage (MIS) 5e fossil localities. Paleo-embayments indicated with asterisks.
Methods
MIS 5e California Region Database
Only fossil collections from marine terrace localities that have been confidently dated to the MIS 5e (or the MIS 5.5 of Martinson et al. Reference Martinson, Pisias, Hays, Imbrie, Moore and Shackleton1987 [130,000–115,000 years ago; Muhs and Groves Reference Muhs and Groves2018]) using uranium-thorium radiometric dating and/or amino acid racemization that show no evidence of substantial temporal mixing were included in the MIS 5e California region database (Supplementary Appendix 1). Corresponding MIS 5e faunal lists were gathered by surveying the published literature and museum collections for the California region (Supplementary Appendix 2). In total, the MIS 5e database contains 142 taxonomically standardized (WoRMS Editorial Board 2018, with minor changes to reflect updates in the published taxonomic literature when necessary) bivalve species from 22 localities spread within and bounding the modern California region that ranges from Punta Eugenia (27.8°N) to Point Conception, California (34.45°N; Supplementary Appendices 1, 2). Because offshore islands can contain faunas distinct from the adjacent mainland coast (Valentine Reference Valentine1966; Lindberg et al. Reference Lindberg, Roth, Kellogg, Hubbs and Powers1980; McLean and Coan Reference McLean and Coan1996), we also compiled exhaustive modern faunal lists for each island in the MIS 5e dataset (i.e., Channel Islands and Guadalupe Island; Supplementary Appendix 3). To identify MIS 5e extralimital species, modern geographic distributions and island occurrences were compared with MIS 5e distributions and island occurrences. Localities were classified as either paleo-embayments or paleo–open coast habitats based on the published literature, wherein paleo-embayments are distinguished on the basis of fine-grained sediments and regional geomorphology and paleo–open coast habitats are characterized by coarser grain sediments and other depositional features characteristic of high-energy, nearshore settings (for further detail, refer to papers listed in Supplementary Appendix 1).
Contemporary Northern Panamic Region Database
To determine what, if any, distributional (i.e., environmental and geographic) characteristics distinguish MIS 5e extralimitals from species that had the potential to occupy the California region during MIS 5e but did not, we compiled a comprehensive database of marine bivalve species that are known to occur within water depths of less than 20 m (a cutoff that is comparable to the maximum inferred depths of marine terrace deposits in the MIS 5e dataset (Kanakoff and Emerson Reference Kanakoff and Emerson1959 [~18 m max.]; Powell et al. Reference Powell, Grant and Conkling2004 [~20 m max.]), in the region most adjacent to southern California, that is, the Northern Panamic region (Baja California Sur and Gulf of California; Fig. 1, Supplementary Appendix 3). Because the MIS 5e dataset contained only one northern extralimital species (see “Results”), we did not compile a northern California/Oregon region fauna. A few species in the dataset have their northernmost Pacific main coast range endpoints within the Northern Panamic region but also occur on Guadalupe Island (29°N); we include such species in the Northern Panamic region fauna, because offshore islands are known to house faunas of disjunct mainland regional affinities (Valentine Reference Valentine1961; Lindberg et al. Reference Lindberg, Roth, Kellogg, Hubbs and Powers1980). Species with maximum body sizes less than 5 mm are excluded from this study, because such small-bodied species are substantially underrepresented in the marine terrace fossil record (Valentine Reference Valentine1989). These constraints help to ensure that only species found within a reasonable distance from the modern California region (and therefore with a reasonably high probability of invading the California region compared with species living farther south) and from depths comparable to the MIS 5e marine terrace record were included in the Northern Panamic region species pool. The Northern Panamic region faunal database contains 343 bivalve species.
Contemporary Species Geographic Distributions
Modern geographic distributions along the Pacific mainland coast and Gulf of California (eastern and western coasts, respectively) were gathered from scientific monographs and databases (Supplementary Appendix 3). We collected occurrence data for eight major offshore Pacific islands off North and South America (Guadalupe, Rocas Alijos, Tres Marías, Revillagigedo, Clipperton Atoll, Coco, Malpelo, and Galapagos Islands) from published faunal surveys, monographs, and scientific papers (Supplementary Appendix 3). Two aspects of geographic range were extracted from these data: latitudinal range and oceanic island occupancy. We focus on these geographic measures, because latitudinal range is the most robust measure of geographic distribution along the north-south-trending mainland Pacific coast and the number of oceanic islands occupied is a measure of dispersal potential and establishment success in offshore coastal habitats (Lindberg and Lipps Reference Lindberg, Lipps, Jablonski, Erwin and Lipps1996). Latitudinal range was taken as the maximum minus the minimum latitude at which each species has been recorded. Oceanic island occupancy is taken as the number of major offshore islands a species is known to occur on (see previous list).
Contemporary Species Environmental Distributions
To determine present-day environmental ranges, we assembled annual and seasonal (maximum and minimum average monthly) chlorophyll α concentration (a proxy for primary productivity) and sea-surface temperature (SST) from Giovanni v. 4.17.2 using the MODIS Aqua 4 km grid cell resolution data from 2003–2015 (Acker and Leptoukh Reference Acker and Leptoukh2007). We focus on SST and productivity, because they have been found to be the most important environmental variables in controlling global shallow-marine biogeographic structure (Belanger et al. Reference Belanger, Jablonski, Roy, Berke, Krug and Valentine2012; Fenberg et al. Reference Fenberg, Menge, Raimondi and Rivandeneira2014) and are readily quantifiable from satellite measurements. Each quarter-degree grid cell was merged with the ETOPO1 Global Relief Model (Amante and Eakins Reference Amante and Eakins2009). Only coastal marine grid cells (containing depths of < 20 m) were included to be comparable to the coastal depths preserved within the marine terrace record.
Maximum and minimum annual and seasonal SST and productivity values were computed by averaging the upper and lower 10% of values observed throughout a species’ entire geographic distribution. Because sufficient and reliable occurrence data are not available for all shallow-marine bivalve species, this “range-through” approach makes the simplifying assumption that species range throughout their geographic ranges along the major coastlines of the eastern Pacific (e.g., Belanger et al. Reference Belanger, Jablonski, Roy, Berke, Krug and Valentine2012). However, this assumption is likely to introduce some distortion for species with patchy geographic distributions or specialized habitats. To address this potential inaccuracy, we also compute present-day “range-limiting” environmental conditions by taking the minimum SST and maximum and minimum productivity conditions observed at a species’ most environmentally extreme geographic range endpoint; this endpoint could be on an oceanic island, within the Gulf of California, or along the Pacific main coast, depending on which of these endpoints had the most extreme temperature and productivity conditions, respectively. For instance, if a species ranges from Punta Eugenia to Cabo San Lucas, with Punta Eugenia representing the cooler of the two range endpoints, then the average annual and seasonal temperatures observed at Punta Eugenia were taken as the species’ annual and seasonal range-limiting temperatures, respectively. This method has the strength of relying exclusively on site-specific occurrences (i.e., by definition, species have verified occurrences at their range endpoints) but makes the assumption that range endpoints represent the extremes of a species’ environmental niche. For rare species that are not commonly collected or reported, this range-limiting method may underestimate environmental niches.
Statistical Analyses
Differences between extralimital occurrences within paleo-embayments and paleo–open coast habitats were assessed using the G-test (Sokal and Rohlf Reference Sokal and Rohlf2012). To evaluate the relative importance of predictor variables (oceanic island occupancy, latitudinal range, chlorophyll α concentration, and SST distributions) in explaining MIS 5e extralimital status, we used multiple logistic regressions as implemented in base R (R Core Team 2019). Two sets of multiple logistic regressions were run: one that contained combinations of geographic and range-through environmental conditions, and another that contained combinations of geographic and range-limiting environmental conditions. Collinear variables (R 2 > 0.5) were identified using the cor.test function as implemented in base R (R Core Team 2019). Support for models with all noncollinear combinations of predictor variables were assessed using a stepwise (both forward and backward) regression framework as implemented in base R. Standardized (beta) coefficients were computed for the best-supported models using the beta function as implemented in the reghelper R package.
Although harder to interpret in a classical statistical framework, classification-based machine-learning approaches generally outperform logistic regression models at prediction, in part because they do not assume a linear relationship between predictor variables and log-odds of the response variable (Elith et al. Reference Elith, Leathwick and Hastie2008). In addition, they allow relationships between predictor and response variables to be readily visualized via partial dependence plots. We therefore used a gradient-boosting machine (GBM) model (Elith et al. Reference Elith, Leathwick and Hastie2008) based on the predictors included in the best-supported logistic regression models to visualize relationships between these predictors and the probability of being an MIS 5e extralimital in California. We then used this model to evaluate which current Northern Panamic species might have the greatest likelihood of expanding their ranges into California as the region warms. The GBM model was R implemented using the gbm.step() function in the R package dismo (Hijmans et al. Reference Hijmans, Phillips, Leathwick and Elith2015). All data, R code, and referenced appendices are available at Dryad: https://doi.org/10.5061/dryad.kprr4xh28.
Results
MIS 5e Extralimital Occurrences
Twenty-one species found as fossils in our MIS 5e database (n = 142) have southern extralimital ranges, that is, they no longer exist in the California region and today occur exclusively to the south. Our MIS 5e dataset contains only one northern extralimital species (Patinopecten caurinus; from Kanakoff and Emerson's [1959] 66-2 fossil locality), which is a robustly shelled temperate-water species that is common in older fossil deposits in California but is not found in any other MIS 5e fossil deposits included herein. While 21 species may appear to be a relatively low number of extralimital species, this amounts to ~15% species difference between the MIS 5e and present-day California region, a percentage that is higher than, though comparable to, previous faunal difference estimates (~12%) for the marine terrace record over the late Pleistocene (Roy et al. Reference Roy, Jablonski and Valentine1995). Sixteen out of 21 southern MIS 5e extralimitals species in our dataset are found within localities interpreted as paleo-embayments (in some San Pedro [Muhs and Groves Reference Muhs and Groves2018], Newport Beach [Kanakoff and Emerson Reference Kanakoff and Emerson1959; Powell et al. Reference Powell, Grant and Conkling2004], and Carmel Valley [Kern Reference Kern1971] localities); the difference between paleo-embayment (n = 9) and paleo–open coast (n = 13) fossil localities in terms of extralimital occurrences is marked (G = 16.27, p < 0.001; Table 1).
Paleo-embayment vs. paleo–open coast Marine Isotope Substage (MIS) 5e extralimital occurrences. G-test results for extralimital presence vs. absence at paleo-embayment (n = 9) and paleo–open coast (n = 13) MIS 5e localities. G, G-test statistic.
Extralimital Characteristics
Of 21 pairwise comparisons of environmental and geographic distribution predictor variables, we find that most SST predictor pairs are correlated (Table 2); these pairs were thus excluded from multiple regression models and assessed separately. Models containing range-limiting environmental condition predictors perform better than models containing range-through environmental conditions (Table 2). Both range-through and range-limiting models are qualitatively comparable in their ranking of each predictor variable's relative importance: overall, models containing mean annual SST perform better than winter and summer SST and all show the expected negative relationship (i.e., cooler temperatures differentiate MIS 5e extralimitals). Seasonal minimum and maximum chlorophyll α concentrations also show their expected positive relationships (i.e., higher chlorophyll α concentrations differentiate extralimitals), with minimum performing better than maximum chlorophyll α concentration (Table 3). Of the geographic predictors, latitudinal range exhibits little predictive value, but the number of oceanic islands is well-supported, with extralimitals tending to occur on more oceanic islands compared with the broader Northern Panamic species pool (Table 2).
Environmental and geographic predictors correlation matrix. Pairs with statistically significant (i.e., p < 0.05) correlations that are greater than |0.5| are shown in bold; these pairs were excluded from regression analyses. SST, sea-surface temperature.
Parameter estimates for the best-supported logistic regression models including range-limiting and range-through environmental conditions, respectively. Standardized beta coefficients shown. SST, sea-surface temperature.
Both forward and backward stepwise regressions prefer models containing annual SST, minimum and maximum chlorophyll α concentration, and number of oceanic islands occupied (Table 3). Range-limiting annual SST is the predictor that most strongly differentiates MIS 5e extralimitals from the broader Northern Panamic fauna (beta = −1.43; Fig. 2, Table 3), followed by number of oceanic islands occupied (beta = 0.62; Fig. 2, Table 3), and minimum chlorophyll α (beta = 0.45; Fig. 2, Table 3); maximum chlorophyll α concentration is also marginally supported (beta = 0.27; p = 0.13). Partial dependence plots based on the GBM model containing these predictors (Fig. 3) show that Northern Panamic species with a modern mean annual temperature range < ~17°C were substantially more likely to occur in California during MIS 5e than species with higher mean annual range-limiting temperatures. Additional marginal increases are associated with currently occupying four or more oceanic islands and (weakly) with currently occupying relatively productive waters (minimum chlorophyll α concentration > ~0.42 and maximum chlorophyll α concentration > ~17). Using this model to predict relative probability of being extralimital for all Northern Panamic species shows the relative performance of the model in discriminating extralimital species from the broader Northern Panamic fauna and also highlights the non-extralimital species with distributional characteristics that are most similar to MIS 5e extralimital species (Fig. 4).
Best-supported Marine Isotope Substage (MIS) 5e extralimital model predictor frequency distributions. Histograms of MIS 5e extralimital species (top; n = 21) compared with frequency distributions of non-extralimital Northern Panamic region species pool (bottom; n = 322) for the three predictors included in the best-supported MIS 5e extralimital model: range-limiting minimum annual sea-surface temperature (left), number of oceanic islands occupied (center), and range-limiting minimum chlorophyll α concentration (right).
Partial dependence plots for all four predictors selected for inclusion in the gradient-boosting machine (GBM) models. Plots show the marginal effect of variation in the predictor variable on the probability of being an Marine Isotope Substage (MIS) 5e southern extralimital. The limited range of the y-axis reflects the fact that no single predictor variable perfectly separates extralimital and non-extralimital species, but the strongest single predictor is mean annual sea-surface temperature (SST) across the modern range of the species—Northern Panamic species with mean annual SST< ~17°C were substantially more likely to occur in California during MIS 5e than species with higher mean annual temperatures.
Predicted relative probability of colonizing California for all species in the present-day Northern Panamic species pool, based on the best-supported gradient-boosting machine (GBM) model (extralimital status ~ range-limiting minimum annual temperature + minimum chlorophyll α concentration + maximum chlorophyll α concentration + number of oceanic islands occupied). Black distinguishes species known to have been actual Marine Isotope Substage (MIS) 5e extralimitals from non-extralimital species (gray). Inset: The 20 species with highest modeled relative probabilities and their extralimital status.
Discussion
Regional Climate Implications of MIS 5e Extralimital Occurrences
The preponderance of southern extralimital species in the MIS 5e record of California (15%, n = 142; Supplementary Appendix 3) appears to support the commonly held view that global warming during MIS 5e resulted in warmer-than-present conditions within coastal Californian marine environments (Muhs et al. Reference Muhs, Simmons and Steinke2002, Reference Muhs, Simmons, Schumann, Groves, Mitrovica and Laurel2012, Reference Muhs, Groves and Schumann2014). However, these findings are also consistent with an alternative MIS 5e climate scenario, the “Bakun hypothesis,” that has received mixed support from some recent oceanographic models (Sydeman et al. Reference Sydeman, Garcia-Reyes, Schoeman, Rykaczewski, Thompson, Black and Bograd2014; Wang et al. Reference Wang, Gouhier, Menge and Ganguly2015). This scenario posits that warm continental temperatures intensified Ekman transport and coastal upwelling during MIS 5e, thus creating relatively cool conditions within exposed outer-coast environments (Muhs and Kyser Reference Muhs and Kyser1987; Muhs et al. Reference Muhs, Miller, Whelan and Kennedy1992; Ortlieb et al. Reference Ortlieb, Diaz and Guzman1996), but warmer-than-present conditions within embayments sheltered from cool upwelling currents (Bakun Reference Bakun1990; Muhs et al. Reference Muhs, Miller, Whelan and Kennedy1992; Bakun et al. Reference Bakun, Black, Bograd, García-Reyes, Miller, Rykaczewski and Sydeman2015). The pattern of MIS 5e southern extralimital occurrences lends some support to this hypothesis (Table 1): 16 of 21 southern extralimital species in the confidently dated MIS 5e database were found within paleo-embayment environments (in some San Pedro localities [Muhs and Groves Reference Muhs and Groves2018], Newport Beach [Kanakoff and Emerson Reference Kanakoff and Emerson1959; Powell et al. Reference Powell, Grant and Conkling2004], and Carmel Valley [Kern Reference Kern1971] localities), implying that outer-coast SSTs were more similar to today even while embayments were warmer.
Range-through versus Range-limiting Environmental Ranges
Our analyses show better performance for models with environmental range predictors based on range-limiting rather than range-through assumptions (Table 3), suggesting that the range-through approach may obscure important distinctions by overestimating ranges. However, the fact that both range-through and range-limiting approaches rank the importance of environmental condition predictors similarly (Table 3) suggests that a common signal is apparent in both approaches.
Applying Contemporary Environmental Ranges to Understand MIS 5e Environmental Conditions
An assumption of the present study is that contemporary environmental ranges of extralimital species can be used to make inferences about the climatic and consequent biogeographical dynamics of the MIS 5e. This assumption could be misleading if environmental conditions in the California region during the MIS 5e were markedly different from modern conditions anywhere in the eastern Pacific (i.e., “no-analogue,” sensu Jackson and Overpeck Reference Jackson and Overpeck2000; Williams and Jackson Reference Williams and Jackson2007). However, the fact that the faunal composition of the California region today differs only modestly from that during MIS 5e (85% faunal similarity between MIS 5e and today; see also Valentine Reference Valentine1989; Valentine and Jablonski Reference Valentine, Jablonski, Ricklefs and Schluter1993) suggests that this is unlikely to be the case. So-called thermally disjunct assemblages containing both northern and southern extralimital species (Woodring Reference Woodring1951; Valentine Reference Valentine and Sears1959, Reference Valentine1973; Valentine and Emerson Reference Valentine and Emerson1961; Lindberg and Lipps Reference Lindberg, Lipps, Jablonski, Erwin and Lipps1996) were once interpreted as representing no-analogue conditions in California during the late Pleistocene, but these assemblages have since been shown to reflect multi-millennial-scale time-averaging due to reflooding of MIS 5e terraces during the cooler MIS 5c sea-level highstand (Muhs et al. Reference Muhs, Simmons, Kennedy, Ludwig and Groves2006, Reference Muhs, Simmons, Schumann, Groves, Mitrovica and Laurel2012, Reference Muhs, Groves and Schumann2014; Muhs and Groves Reference Muhs and Groves2018).
A second consideration with respect to drawing inferences about MIS 5e from contemporary environmental range data is that this approach may be invalid if species’ environmental niches have evolved considerably since MIS 5e. The stability of environmental niches over geologic timescales is an area of much recent research (Pearman et al. Reference Pearman, Guisan, Broennimann and Randin2008; Nogués-Bravo Reference Nogués-Bravo2009; Stigall Reference Stigall2014; Maguire et al. Reference Maguire, Blois, Nieto-Lugilde, Fitzpatrick and Williams2015; Lieberman and Saupe Reference Lieberman and Saupe2016; Qiao et al. Reference Qiao, Saupe, Soberón, Peterson and Myers2016; Waterson et al. Reference Waterson, Edgar, Schmidt and Valdes2017). Most studies have found that environmental niches are relatively stable over Milankovitch scale cycles (Martinez-Meyer and Peterson Reference Martinez-Meyer and Peterson2006; Pearman et al. Reference Pearman, Guisan, Broennimann and Randin2008; Saupe et al. Reference Saupe, Hendricks, Portell, Dowsett, Haywood, Hunter and Lieberman2014, Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015), although realized environmental niches may shift through time (Veloz et al. Reference Veloz, Williams, Blois, He, Otto-Bliesner and Liu2012). These emerging findings support the use of contemporary environmental ranges to help infer the climatic and biogeographic dynamics of the past, at least on Holocene–Pleistocene timescales.
Predictors of MIS 5e Extralimitals
Previous work showed that late Pleistocene southern extralimital species tended to have broader latitudinal ranges than non-extralimitals in the Northern Panamic species pool (Roy et al. Reference Roy, Jablonski and Valentine1995). In our analyses, latitudinal range does not emerge as an important predictor of MIS 5e southern extralimitals when environmental range predictors, which are variably correlated with latitudinal range (R 2 values range from 0.77 to 0.03 for pairwise correlations; Table 2), are also considered. Of the correlates of latitudinal range, cooler annual SST minima and greater productivity distinguish MIS 5e extralimitals from the broader Northern Panamic species pool (Fig. 2, Table 3). These findings suggest that MIS 5e extralimital species have an affinity or tolerance for relatively productive habitats and cool temperatures —both of which characterize the California region's environmental conditions (Belanger et al. Reference Belanger, Jablonski, Roy, Berke, Krug and Valentine2012; Fenberg et al. Reference Fenberg, Menge, Raimondi and Rivandeneira2014). The greater importance of annual SST compared with seasonal SSTs is also consistent with this interpretation, because the California region is substantially less seasonal than the Gulf of California (where most species in the broader Panamic species pool live today). The Pacific coast of southern California and the Gulf of California experience similar winter SSTs but summer SSTs are much lower on the Pacific coast of southern California than in the Gulf of California. This may explain why mean annual temperature is a more important determinant of MIS 5e extralimital status than minimum winter temperature. Thus, our findings support previous inferences about the underlying drivers of range shifts in coastal California (Roy et al. Reference Roy, Jablonski and Valentine1995) and also provide direct paleontological support for global-scale analyses suggesting that temperature plays the dominant role in governing marine biogeographic structure (Belanger et al. Reference Belanger, Jablonski, Roy, Berke, Krug and Valentine2012; Sunday et al. Reference Sunday, Bates and Dulvy2012), with productivity an important but secondary factor (Belanger et al. Reference Belanger, Jablonski, Roy, Berke, Krug and Valentine2012).
We also find that MIS 5e extralimital species occupy significantly more offshore islands than non-extralimitals in the Northern Panamic species pool (Fig. 2, Table 3). Offshore island occupancy in marine mollusks is strongly related to ecological traits such as long-lived planktonic larvae, rafting adult habitats (i.e., some wood-boring or kelp-inhabiting taxa), and generalist feeding ecologies (e.g., Hansen Reference Hansen1980; Jablonski Reference Jablonski1986; Lester et al. Reference Lester, Gaines and Kinlan2007) that influence larval and adult dispersal (Strathmann Reference Strathmann1987; Granthman et al. Reference Grantham, Eckert and Shanks2003). Thus, the emergence of island occupancy as a significant predictor of MIS 5e extralimitals, even while accounting for temperature and productivity ranges, points to a key role for dispersal potential in facilitating climate tracking and consequent species emigrations (see also Lindberg and Lipps [Reference Lindberg, Lipps, Jablonski, Erwin and Lipps1996] and Dobrowski et al. [Reference Dobrowski, Thorne, Greenberg, Safford, Mynsberge, Crimmins and Swanson2011] for similar conclusions).
A Pre-anthropogenic Baseline for Contemporary Warming
Globally, the MIS 5e was a period of climatic warming, with the warmest temperatures of the last 300,000 years (Lisiecki and Raymo Reference Lisiecki and Raymo2005). Thus, analysis of extralimitals during the MIS 5e can serve as a pre-anthropogenic analogue to cast baseline expectations for contemporary climate change–driven range shifts into southern California. Our analysis of MIS 5e extralimitals supports the dual importance of environmental niches and life-history attributes influencing dispersal, thus implying that both characteristics play complementary roles in determining which species responded to MIS 5e pre-anthropogenic warming by successfully tracking their niches into the California region. Predicted relative probabilities from the GBM model can be used to both evaluate overall model performance (extralimital species have relatively high modeled probabilities in comparison to the majority of species within the Northern Panamic species pool) and to highlight species that might be expected to expand their ranges into California in a warming climate. These include MIS 5e extralimitals themselves as well as species that are not known to have occupied California during MIS 5e but have geographic and environmental distributional characteristics similar to MIS 5e extralimitals (Fig. 4).
How accurately models based on MIS 5e range shifts might perform in predicting contemporary and future change is an open question (Jackson and Blois Reference Jackson and Blois2015; Nogués-Bravo et al. Reference Nogués-Bravo, Rodríguez-Sánchez, Orsini, de Boer, Jansson, Morlon, Fordham and Jackson2018). Human impacts such as coastal habitat degradation, pollution, introduced species, and diseases, as well as food web disruption, will invariably alter the pre-Anthropocene rules of climate-mediated geographic range shifts. Our MIS 5e baseline predictions are grounded in a pre-human warming event, and as such, we hold no expectation that they should be exactly replicated in response to anthropogenic global change. However, as the pace of contemporary climate change accelerates, the composition of species that accumulate along California's shore will serve to test our pre-anthropogenic baseline expectations.
Conclusions
The MIS 5e bivalve fauna of southern California is characterized by the presence of southern extralimital species (15%), the majority of which occurred in paleo-embayment, rather than paleo–open coast settings (16 of 21; Table 1). Thus, MIS 5e appears to have coincided with warmer-than-present conditions within southern California's paleo-embayments, but more similar conditions relative to today within open-coast habitats. Extralimital bivalves differ from their contemporary species pool: extralimitals are distinguished by having more numerous offshore island occurrences and cooler, more seasonally productive environmental limits (Fig. 2). These findings suggest that southern extralimitals that occupied the California region during MIS 5e have life-history traits that enable high dispersal rates and environmental niches that include affinities to conditions quite similar to the California's more temperate habitats. Thus, MIS 5e extralimitals, and non-extralimitals with similar characteristics (Fig. 4, Supplementary Appendix 3), may be more likely to colonize the California coast as the contemporary climate continues to warm. More broadly, our analyses provide direct paleontological support for the inference that SSTs play a dominant role in governing marine biogeographic structure through time and bolster efforts to use present and past species distributions to predict future range shifts.
Acknowledgments
C. Powell II, D. Muhs, D. Lindberg, J. Valentine, and J. Lipps provided invaluable advice and expertise on Pleistocene marine geomorphology and paleozoogeography. We gratefully acknowledge D. Muhs and L. Thomas for sharing faunal lists and coordinating fieldwork on San Nicolas Island. E. Clites, P. Holroyd, A. Hendy, and M. Stecheson facilitated the incorporation of museum collections at the University of California at Berkeley, Museum of Paleontology and Natural History Museum of Los Angeles County, Invertebrate Paleontology section. The authors thank three anonymous reviewers and A. Stigall for thoughtful critiques. This work was supported by grants made to E.A.O. from the American Natural History Museum, American Philosophical Society, Conchologists of America, Evolving Earth Foundation, National Geographic Society Young Explorers Program, Natural History Museum of Los Angeles County, Northern California Geological Society, Sigma Xi, University of California Museum of Paleontology, and the Western Society of Malacologists. This study was carried out when E.A.O. was supported by Exploring California Biodiversity GK-12 Fellowship, University of California Museum of Paleontology Fellowship, and the National Science Foundation (EAR 1740214).
