Anthropogenic climate change is altering global biogeographical patterns. However, it remains difficult to quantify how bioregions are changing because pre-industrial records of species distributions are rare. Marine microfossils, such as planktonic foraminifera, are preserved in seafloor sediments and allow the quantification of bioregions in the past. Using a recently compiled data set of pre-industrial species composition of planktonic foraminifera in 3802 worldwide seafloor sediments, we employed multivariate and statistical model-based approaches to study spatial turnover in order to 1) quantify planktonic foraminifera bioregions and 2) understand the environmental drivers of species turnover. Four latitudinally banded bioregions emerge from the global assemblage data. The polar and temperate bioregions are bi-hemispheric, supporting the idea that planktonic foraminifera species are not limited by dispersal. The equatorial bioregion shows complex longitudinal patterns and overlaps in sea surface temperature (SST) range with the tropical bioregion. Compositional-turnover models (Bayesian bootstrap generalised dissimilarity models) identify SST as the strongest driver of species turnover. The turnover rate is constant across most of the SST gradient, showing no SST threshold values with rapid shifts in species composition, but decelerates above 25°C, suggesting SST is less predictive of species composition in warmer waters. Other environmental predictors affect species turnover non-linearly, and their importance differs across regions. In the Pacific ocean, net primary productivity below 500 mgC m⁻² day⁻¹ drives fast compositional change. Water depth values below 3000 m (which affect calcareous microfossil preservation) increasingly drive changes in species composition among death assemblages in the Pacific and Indian oceans. Together, our results suggest that the dynamics of planktonic foraminifera bioregions are expected to be highly responsive to climate change; however, at lower latitudes, environmental drivers other than SST may affect these dynamics.

All data used in this paper have been published before elsewhere, please cite them appropriately.

For planktonic foraminifera relative abundance data (ForCenS), see Siccha, M., & Kucera, M. 2017 PANGAEA https://doi.org/10.1594/PANGAEA.873570.

For environmental data, see:

• World Ocean Data (WOA) 2018: https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/
• Historical sea-surface temperature: https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.v5.html
• Ocean Productivity: http://sites.science.oregonstate.edu/ocean.productivity/

Code "0_main.R" calls (source) the other codes, so start with this code. Most of the analyses require high-performance computing, so I suggest decreasing the size of the data or the length of analyses (e.g. bootstraps) when running codes on your personal computer.

The code provided here (Zenodo) downloads and analyses these data, except for the Ocean Productivity data that is already provided as "forcens_pp.csv". The forcens_pp dataset contains net primary productivity (NPP) data for each of the ForCenS samples used. Please contact the first author if you would like the R code to obtain annual mean NPP data. Explanation of the forcens_pp columns:

• Sample_ID: identification code of the sample from the ForCenS dataset (see above)
• Latitude and Longitude: coordinates of the sample/site
• site: unique identification number of each ForCenS site given by us
• pp: annual mean values of NPP from the Ocean Productivity project (see above), using the standard VGPM algorithm for each site