Scale dependence of temporal biodiversity change in modern and fossil marine plankton

Aim: Biodiversity dynamics comprise evolutionary and ecological changes on multiple temporal scales from millions of years to decades, but they are often interpreted within a single time frame. Planktonic foraminifera communities offer a unique op-portunity for analysing the dynamics of marine biodiversity over different temporal scales. Our study aims to provide a baseline for assessments of biodiversity patterns over multiple time-scales,

Clearly, a proper baseline is needed to make accurate predictions on how climate change and anthropogenic activities will influence biodiversity in the future (Cardinale et al., 2018;Jackson, 2008).
Defining a clear baseline for biodiversity before industrialization is difficult, because most monitoring programs post-date the onset of the industrial revolution (Dornelas et al., 2014;Elahi et al., 2015;Lotze et al., 2006). As a result, we lack information on the variation in species numbers and community composition (as one aspect of biodiversity) before human impact. However, for some taxonomic groups, the fossil record provides the necessary temporal context encompassing the range of ecological variation on time-scales from decades to millions of years and, therefore, there is increasing interest in using palaeontological records to answer ecological questions (Fritz et al., 2013;Harnik et al., 2012;Nogués-Bravo et al., 2018;Yasuhara et al., 2017). Recent findings by Jonkers, Hillebrand, and Kucera (2019) confirmed that composition of modern foraminifera assemblages differ from those in pre-industrial times as a result of temperature change, highlighting the importance of historical baselines for the assessment of ecological communities.
Another critical aspect of studies reporting anthropogenic impact on biodiversity is the focus on local species richness (Chase & Knight, 2013;Hillebrand, Blasius, et al., 2018). Recent analyses reveal that some communities are undergoing major compositional change without change in species number (Hillebrand, Blasius, et al., 2018;Magurran et al., 2018). This is because colonization by new species sometimes equals or exceeds rates of local extinctions (Sax & Gaines, 2003), so that local richness fluctuates around a steady state, despite changes in community composition and species abundance.
Reordering of species dominance patterns can occur along with changes in species identity (species appearance and disappearance) or can be independent of identity change. If there is a change in dominance, but no change in species presence/absence, biodiversity change is driven by reordering the relative abundances of species already present in the local community (Hillebrand, Blasius, et al., 2018).
Therefore, an analysis of baseline biodiversity change that is relevant for biodiversity conservation has to consider both richness (presence/ absence) and dominance (relative abundance) (Barnosky et al., 2017;Cardinale et al., 2018). Accordingly, we here focus on two metrics of biodiversity change: (a) change in species presence/absence (species turnover measured by Jaccard's dissimilarity), and (b) reordering of species dominance patterns (Wishart's dissimilarity) as described by Hillebrand, Blasius, et al. (2018). For completeness, we also discuss changes in species richness (as it is the most common measure of biodiversity), and the effective number of species (ENS), which has been suggested to be a more robust diversity metric than richness regarding different sampling schemes (Chase & Knight, 2013). By comparing the identity-based metrics (richness, Jaccard) with dominance-based metrics (ENS, Wishart), we provide a comprehensive assessment of the magnitude of biodiversity change across time-scales.
We use planktonic foraminifera as a model system to study scale dependence of temporal biodiversity patterns in the marine realm.
Planktonic foraminifera are open-ocean protists that have the most complete taxonomic record of any clade, which allows for direct comparisons between recent and geologically distant communities (Fenton, Pearson, Dunkley Jones, Farnsworth, et al., 2016). They are widely sampled both in the recent past and throughout deeptime, so that rich datasets on foraminiferal assemblages are available across large spatial and temporal scales Yasuhara et al., 2017). Planktonic foraminifera are sensitive to environmental change, such as temperature (Ezard, Quental, & Benton, 2016;Jonkers et al., 2019), and are therefore widely used for palaeoclimate reconstructions (Kucera, 2007). In this study, we analysed data on community composition of planktonic foraminifera over the past 7 million years to track dynamics of biodiversity change in the ocean. Rather than defining the drivers of biodiversity change, of which there are many, we assessed the magnitude of biodiversity fluctuations across five distinct temporal scales from decades to millions of years. Benefitting from access to modern monitoring time series as well as to fossil data, we analysed community turnover in planktonic foraminifera at these time-scales, specifically asking how the magnitude of change detected in modern time series of species turnover compares with observations in the fossil record. In particular, we tested the following hypotheses: Hypothesis 1: In the multi-million-year time series, an increasing number of speciation and extinction events increases temporal turnover based on presence-absence, and might increase the overall probably forced by anthropogenic effects, than those observed over much longer durations.

K E Y W O R D S
biodiversity, community turnover, global change, planktonic foraminifera number of species observed at locations where the speciation rates exceed the extinction rates. For dominance-based metrics (ENS, Wishart's dissimilarity), we expect similar patterns but even higher turnover as -in addition to speciation and extinction -changes in relative abundance will accumulate over time.
Hypothesis 2: In the more recent records (million-year to decadal scales), where speciation and extinction do not affect temporal biodiversity patterns (Yasuhara, Hunt, Dowsett, Robinson, & Stoll, 2012;Yasuhara et al., 2017), the number of species will increase and turnover will decrease in longer time series that are intrinsically characterized by higher averaging than shorter time series.
In addition, we compared the observed changes in biodiversity to changes in temperature over the same time, testing the Hypothesis 3 that global temperature cycles are associated with higher rates of turnover during warming events. We expect that these effects are most pronounced during ice ages and after the last glaciation (million-year and multi-millennial scales) characterized by the strongest temperature fluctuations.

| Data sources
We compiled data on foraminifera community composition (species' relative abundances) from 20 sediment cores and 6 sediment traps with a minimum length of 5 years distributed over the world's oceans ( Figure 1). Only cores that contained the full taxonomic record were included in the analysis. Cores were selected to explore different scales of temporal variability in foraminifera community composition, thus datasets that had too low sampling resolution to cover this variability (e.g. too short records) were omitted. The datasets (Supporting Information Appendix S1) were grouped into five categories depending on how far back in time they extend: multi-million-year scale (0 to 7 Myr, eight cores), million-year scale (0 to 0.5 Myr, four cores), multi-millennial scale (0 to 15 thousand years, three cores), millennial scale (0 to 1,100 years, five cores) and decadal scale (0 to 32 years, six sediment trap records). The longest datasets were characterized by the highest sample integration (highest time averaging, Supporting Information Appendix S1: Table   S1). All data from sediment traps were analysed within the decadal scale category. From the sediment traps we extracted yearly averages of the foraminifera fluxes to eliminate seasonal dynamics (Jonkers & Kučera, 2015). The year 2016 was set as 'present day' for all datasets (cores and traps) to facilitate a direct comparison of the results.
We also compiled estimates of global temperature for each time frame (Supporting Information Appendix S1), except for the 0-7 Myr period for which no global compilation is available. Instead, we used a tropical (almost exclusively Pacific) temperature stack , which is appropriate as all data that span this period are from this basin.
F I G U R E 1 Geographical distribution of time series analysed in this study. Names of the sites and their desription are listed in the Supporting Information Appendix S1: Table S1 [Colour figure can be viewed at wileyonlinelibrary.com]

| Scaling in a neutral community model
To create reference time series for different scales, we designed a stochastic process simulating a neutral community model (for more detailed description, please see Supporting Information Appendix S2). In this model stochastic speciation events are introduced through a constant speciation rate that is balanced by the extinction rate in order to generate a species richness that fluctuates around the level of 20. When a speciation event occurs, a new species is released at the detection limit (zero). When species abundance declines crossing the detection limit, it is recorded with the value zero, but left in the simulation with a chance of reoccurrence. Eventually, when it hits the extinction threshold (at We used this model to simulate time series and infer frequency distributions of diversity metrics at different time-scales (time-series lengths 0.1, 1, 10, 1,000 and 10,000 kyr, corresponding to decadal, millennial, multi-millennial, million-year and multi-million-year timescale categories, respectively). We used six, five, three, four and eight simulation runs for decadal, millennial, multi-millennial, million-year and multi-million-year time-scale categories, respectively.
In contrast to the empirical data, time resolution was adapted to render 100 equidistant time points per simulation run with a sampling interval of 0.001, 0.01, 0.1, 10 and 100 kyr, respectively. This means that going from the shortest to the longest scale the simulated time series were uniformly downsampled by a factor of 10 -5 and time-series length is positively correlated with sampling interval.
The initial foraminifera community contained 20 species, except for the decadal scale records, for which we started with 6, 6, 11, 11, 12 and 13 species in the related six runs.

| Statistical analysis
We calculated species richness as the number of morphospecies (fossil species can only be defined using the morphological species concept) identified in each sample. The ENS was calculated as described by Chase and Knight (2013) (see also Supporting Information Appendix S2). The advantage of using ENS over richness is that ENS is less sensitive to sampling effort and does not overestimate the impact of rare species on the overall diversity.
Jaccard's dissimilarity and Wishart's dissimilarity were calculated as described by Hillebrand, Blasius, et al. (2018). Jaccard's dissimilarity is a presence-absence-based species exchange ratio between two samples and represents changes in community composition in terms of species identities, whereas Wishart's dissimilarity is a measure of dominance change (turnover in species proportional abundances). The latter is closely related to the Simpson index and thus related to robust measures of diversity, such as the ENS (Chase & Knight, 2013). For each core and trap series we calculated both dissimilarity metrics comparing foraminifera assemblages between each unique pair of samples (the diagonal of the dissimilarity matrix was always omitted to exclude zero dissimilarities between each sample and itself). Then, we calculated a median of pairwise dissimilarities between a sample of a given age and all other samples in the core or the sediment trap ( Figure 2). The advantage of this procedure is that the temporal trend does not depend on a single reference community at the beginning or at the end of the time series. Finally, for each site we calculated a trend using generalized additive models (GAMs) with cubic regression splines. The error covariance was assumed to follow a time-dependent autoregressive process to account for temporal correlation.
General trends within each five categories were derived by applying GAMs to the fitted values for individual series and site ID was included as a random effect (random intercept) to account for differences between sampling sites (different location and sample integration). We used the GAM approach, because it enables us to account for statistical irregularities in the data, such as non-uniform sampling frequency and interpolating over missing observations. GAMs also allow the response to be estimated as a linear or nonlinear function of the predictors, giving a better chance of detecting non-monotonic patterns, for example, glacial-interglacial fluctuations. Model assumptions, including normality, independence and homogeneity of variance, were verified by visual analysis of residuals. All analyses were performed in R (R Development Core Team, 2017, version 3.3.3) using package 'mgcv' for GAMs (Wood, 2017). More details on statistical procedures are included in Supporting Information Appendix S2.
To further explore the effects of temporal scaling (Tomašových & Kidwell, 2010) we made a simple comparison of the mean diversity values (measured using our four indices) with time-series length, sampling frequency (time separation between the samples) and sample integration (inverse of the sedimentation rate measured as a x n = + A(x n−1 − ) + n first derivative of the depth-age model for each analysed sediment core).
To address whether observed biodiversity patterns were associated with temperature change, we performed cross-correlation analyses for each time category. Results of these analyses are presented in detail in Supporting Information Appendix S2.

| Frequency distribution of diversity metrics
The frequency distribution of diversity measures shows substantial changes in standing diversity, identity and dominance change over is too short to record consistent extirpation or invasion, thus suggesting lower turnover. A distinctly higher turnover is seen on time-scales up to 10 6 years (Figures 3b and 4b), where we also see the strongest response to climatic forcing (Supporting Information Appendix S2: Figure S5).

F I G U R E 2
Calculation of temporal trends in Jaccard's or Wishart's dissimilarities at one site. Community composition between each pair of samples was compared and temporal trends were derived from median dissimilarity at every time step (detailed description in text). For better visualization only the first 9 of the total 78 time steps are presented

| Dominance-based biodiversity change over time
Based on the observed changes in the mode and range of dominance-based diversity metrics across time-scales, we focus on temporal trends of ENS and Wishart's dissimilarity here, but show similar graphs for richness, and Jaccard turnover in Supporting Information Appendix S3. Analysis of foraminifera ENS over the past 7 million years shows nonlinear pattern in single cores, but highlights a consistent diversity maximum around 3 million years ago (Figure 5a), close to the time of major intensification of the Northern Hemisphere glaciation (Figure 5c) (Raymo, 1994), corresponding to the aftermath of the last diversification of planktonic foraminifera after which only extinctions occurred (Fritz et al., 2013). Dominance-based turnover is generally high over these time-scales (Figure 5b, see above), but shows a decline towards the time of the Northern Hemisphere glaciation 3 Ma (Figure 5b). After the glaciation started (Figure 5c), ENS decreased and Wishart's dissimilarity did not change, but remained at a high overall level.

Over glacial-interglacial cycles (million-year time-scale), ENS and
Wishart's dissimilarity show cycles with a frequency corresponding to the glacial cycles and peaks during interglacial periods (Figure 5d-f and Supporting Information Appendix S2: Figure S5). Details differ between cores, but maximum ENS and turnover in species composition is often associated with warming periods. The overall magnitude of turnover is much smaller than on the multi-million-year time-scale (cf. Figure 5b and e, respectively).
In the multi-millennial datasets, temporal turnover dynamics are related to climate instability at the end of the last glacial period c.
12,000 years ago (Figure 5g-i). However, ENS peaked significantly later (c. 10 kyr ago) than the turnover metric (c. 13.5 kyr ago), indicating that reordering of community composition can take place without gaining or losing effective species. In fact, turnover (Figure 5h) is high during the warming phase, when it is in the same range as at In contrast to the significant changes in ENS and Wishart's dissimilarity over longer time-scales, both remained constant over the last hundreds to thousand years (Figure 5j-l). Although small sample size should be acknowledged, differences between sites are obvious, as both metrics are higher and more variable at lower latitudes (sites C90, SO90 and SBB2 compared to MD99_3 and MD99_5). However, the overall trend of ENS and Wishart's dissimilarity remains neutral, likely reflecting the stable climate during this period (Figure 5l and Supporting Information Appendix S2: Figure S2) and no effect of extinction and speciation on this temporal scale.
In the decadal records, ENS increased with time and the shortest trap series showed the highest ENS (Figure 5m). Notably, however, Wishart's dissimilarity remained very high (comparable to glacial-interglacial turnover; Figure 5n).

| Effects of temporal resolution
The duration of the time series had no effect on the mean ENS in the empirical data (Figure 6b). Highest species richness and turnover values (Jaccard's and Wishart's dissimilarity) were observed in the longest records (multi-million-year time series), because of a high number of speciation and extinction events that could be captured on this time-scale. As the longest time series had also significantly longer temporal separation between sampling events, the same pattern was observed for the relationships between biodiversity measures and sampling frequency (Figure 6e-h), that is, time series with extremely low sampling frequency were characterized by highest richness and turnover values.
Fossil assemblages are averaged over a certain amount of time (sample coarsening), which can influence estimates of species turnover (Tomašových & Kidwell, 2010), but such sample integration does not affect biodiversity patterns in our study (Figure 6i-l). This can be probably explained by counteracting effects of duration and sample integration time (Tomašových & Kidwell, 2010) with the duration effects overriding the effects of sample resolution.

| D ISCUSS I ON
Our analysis provides a first evaluation of marine biodiversity patterns across temporal scales from millions of years to decades. It indicates that trends in community turnover are key to understanding these patterns. The magnitude of temporal turnover observed in the decadal time series is smaller in terms of richness, but compositionally larger than what is seen on the millennial time-scale and the range of the compositional turnover is consistent with the climatically driven pattern of glacial and post-glacial time series (Figure 5 and Supporting Information Appendix S2: Figure S5). When considering the records from millennial to multi-million-year time-scales, we find a significant correlation between climatic changes and the turnover in composition and standing diversity of foraminifera in presence-absencebased estimates (richness, Jaccard), as well as in dominance-affected measures (ENS, Wishart) supporting Hypothesis 3 that global warming events are associated with higher rates of turnover. This is visible in Figure 5 and Supporting Information Figure S1 (Appendix S3), and supported by the results of cross-correlation analyses between global temperature records and observed biodiversity patterns (Supporting Information Appendix S2: Figure S5). Here, we focus our discussion on ENS and Wishart's dissimilarity, as these are more likely to be informative on the short time-scales we are interested in for predicting future diversity responses to climate change.
On the geological time-scale of millions of years, ENS peaked at the onset of the Northern Hemisphere glaciation that occurred 3.5 million years ago, and turnover is decreasing in a cooling world. These changes are smaller in magnitude, but consistent at the million-year time-scale: turnover and ENS are higher in warming periods ( Figure 5 and Supporting Information Appendix S2: Figure S5). The magnitude of change decreases when considering multi-millennial time frames and is very small to absent at millennial time-scales. On long time-scales (millions of years), changes in community turnover are in part driven by slowly occurring species extinctions (global or geographically selective) and speciation (evolution), whereas on shorter time-scales climate change plays the dominant role and the turnover dynamics reflect geographical redistribution of species, rather than extinction (Jackson, 2008;Yasuhara, Hunt, Dowsett, et al., 2012). Geographical differences between different sites are most pronounced in millennial records that are characterized by higher values of ENS and Wishart's dissimilarity at low latitudes ( Figure 5, sites C90, SBB2, SO90) F I G U R E 5 Temporal trends of foraminifera biodiversity over different time periods (as in Figure 2) as effective number of species (upper panels) and Wishart's dissimilarity (mid panels) in relation to temperature (bottom panels). Colours represent different sites, which are listed in the Supporting Information Appendix S1: Table S1. Black solid lines are general trends (± 95% confidence intervals) calculated using generalized additive models (GAMs) [Colour figure can be viewed at wileyonlinelibrary.com] compared to high latitudes ( Figure 5, sites MD99_3 and MD99_5), which is in line with the spatial distribution of fossil foraminifera diversity (high tropical and low polar biodiversity, Yasuhara, Hunt, Dowsett, et al., 2012). Climate-driven changes in species community structure have been described for many organism types in single locations or regions (small mammals: McGill, Hadly, & Maurer, 2005;dinoflagellates: de Vernal et al., 2005;corals: Pandolfi & Jackson, 2006;ostracodes: Yasuhara, Cronin, Hunt, & Hodell, 2009;benthic foraminifera: Reymond, Bode, Renema, & Pandolfi, 2011;planktonic foraminifera: Flower & Kennett, 1995;diatoms: Ren, Gersonde, Esper, & Sancetta, 2014;molluscs: Tomašových, Dominici, Zuschin, & Merle, 2014) and recently also for planktonic foraminifera before and after industrialization (Jonkers et al., 2019). However, most studies take a reverse approach and use changes in species assemblages to reconstruct past climate, rather than focusing on patterns of species turnover associated with climate change (Yasuhara, Hunt, Breitburg, et al., 2012). Our study goes beyond earlier work by analysing dominance-driven biodiversity change worldwide on different time-scales.
Consistent with the lack of substantial climatic change over the last millennium, ENS remained stable and turnover was small ( Figure 5). This observation contrasts with the analysis of decadal records. Here, we see strong fluctuations of the dominance structure (Wishart's dissimilarity) independent of changes in ENS, with the magnitude of change comparable to turnover observed over millions of years. This is in line with recent literature showing that community composition is rapidly changing, whereas the number of species remains stable over time (Dornelas et al., 2014;Hillebrand, Blasius, et al., 2018;Magurran et al., 2018). This high modern turnover of assemblages may reflect anthropogenic forcing and the comparison with our neutral scaling model indicates that stochastic processes cannot explain observed variability in foraminifera biodiversity on decadal time-scales ( Figure 4). Thus, the lower turnover on the millennial time-scale is not a result of higher temporal averaging, but reflects ecological processes (Supporting Information Appendix S2: Figure   S5). Since the decadal time series were annually averaged, seasonal variation in modern time series does not explain these differences as they are not reflected in our turnover estimates. These data thus suggest that the high turnover in modern assemblages is a response to fast environmental changes, such as increasing mean and variance in temperature.
Our analysis of foraminifera biodiversity patterns may be influenced by the effects of temporal resolution (time-series length, F I G U R E 6 Effects of temporal resolution on foraminifera biodiversity: (a-d) for a temporal extent of the study; (e-h) for a sampling frequency (temporal distance between sampling events); (i-l) for a sample integration (calculated as the inverse of sedimentation rate). ENS = effective number of species sampling frequency and sample integration) on species turnover (Tomašových & Kidwell, 2010). In line with Hypothesis 1, we observed the highest number of species and highest values of turnover in the multi-million-year time series, but these results are partly due to scaling ( Figure 6). However, scaling effects do not explain either the patterns in ENS (Figure 6) or the high variability in Wishart's dissimilarity in the decadal records ( Figure 3). Thus, Hypothesis 2 (decreasing turnover in modern times as a result of time averaging) can be supported for the presence-absence-based turnover, but not for the dominance-based turnover.
In conclusion, our analysis shows how biodiversity patterns are shaped at different time-scales, which is an important step towards understanding the processes driving diversity dynamics across scales.
On the multi-million-year time-scale, slow processes, such as speciation and global extinction, can be addressed. On the shorter timescales (thousands to hundreds of years), climatic forcing most likely prevails, whereas biodiversity patterns observed on the decadal scale are possibly a result of fast occurring processes, such as species dispersal. Temporal scaling effects (Tomašových & Kidwell, 2010) can explain higher richness and turnover values on the multi-million-year time-scale ( Figure 6), but not the high variation in the dominance structure in modern times. We reinforce the importance of historical baselines of biodiversity change in interpretations of current dynamics and reveal that restructuring of communities occurs at all times, but changes in the modern dominance structure may reach higher magnitudes than expected from Earth's history. Given the rapid change in climate, the critical question in ecology is whether species will be able to adapt or move, and thus communities to persist. Analysis of species richness alone is inconclusive and does not allow predictions on biodiversity dynamics (Cardinale et al., 2018), whereas turnover estimates prove to be more sensitive and consistent in their response to past climate change, with higher turnover during warming periods.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are available in the Dryad data repository (https://doi.org/ 10.5061/dryad.2z34t mphf) and the sources are listed in the Appendix.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.