Predictability of intraspecific size variation in extant planktonic foraminifera

Planktonic foraminifera (PF) size varies greatly both among and within species. This variation affects our understanding of PF ecology and evolution as well as reconstructions of the ocean-climate system. It is currently accepted that PF species are largest under optimum environmental conditions, where abundance is maximised. This idea is based on observations from marine sediment assemblages; however, these observations either had limited intraspecific resolution or focused on a restricted part of each species’ biogeographical range. Here we compile a new global PF shell size dataset to investigate the relationship between intraspecific size variation and abundance and sea surface temperature (SST). Our dataset contains 3817 individual size measurements on nine PF species in 53 surface sediments around the world. For each species, we fitted a generalised linear model of population shell size as function of local abundance (as an indicator of optimum environmental conditions) and SST. We support previous results that species maximum size and maximum abundance rank along SST; however, this relationship is not supported within species. Only two species out of nine revealed a significant positive relationship between size and abundance, suggesting shell size is not maximised at the species environmental optimum. SST significantly explained variation in shell size for four species out of nine. By incorporating intraspecific variation and sampling broader geographical ranges compared to previous studies, we conclude that the relationships between PF shell size and abundance or SST are either absent or weaker than previously reported.

. From these sea-floor sediment samples, we selected those that contained 124 only modern species (Table S2), were collected within the upper 15 cm of sediment, and each of these five population metrics of the Buckley Collection against the re-sampled data, 166 and calculated the residuals based on the identity function (1:1 relationship). The residu-167 als of the regressions were predominantly positive (Fig. S2), indicating that the Buckley 168 Collection has a consistent collector bias towards large specimens. 169 The mean squared error was lowest for the 95th percentile (Fig. S2), meaning that this 170 metric is the most representative population metric of the Buckley Collection. The robust-171 ness of the size distribution's 95th percentile has also been documented by Schmidt et al. 172 (2004), as it is less sensitive to single outliers than the distribution's maximum value, and 173 to representative sampling at the lower end of the size range than the distribution's mean 174 and median values. Accordingly, in our analyses, we used the 95th percentiles of the pop-

191
The spatial arrangement of dead PF on the sea floor is affected during settling by sub-192 surface currents (Berger and Piper, 1972). Recent models estimate that dead foraminiferal  To test for the effect of retrieving relative abundance data of samples 300 km apart, 203 we ran the same analysis using solely the nearest neighbour of the ForCenS database rela-204 tive to each morphometric sample. The median distance between the morphometric samples 205 and their nearest ForCenS neighbours was 106 km. The analyses using the single nearest 206 ForCenS sample produced consistent results when compared to the analyses using all sam-207 ples within a 300 km distance (Section S5). We present results using the more conservative 208 300 km median relative abundance. (1 • is approximately 111 km at the equator). Again, the distances between the datasets were 214 calculated using the WGS 84 system (Hijmans, 2015). We used SST data from the earliest 215 decade available in the WOA13 database, resulting in SST data averaged for the years be-216 tween 1955 and 1964. We chose this time period because the last historical expedition that 217 we used for our morphometric dataset sailed in 1965 (Table S1)

247
In general, intraspecific size variation is high among populations (Fig. 2) and within popu-248 lations (Fig. S3). Among the nine PF species studied, only T. sacculifer and G. truncatuli-249 noides show a statistically significant positive relationship between shell size and relative abundance. Relative abundance never explains more than 7% of population shell size vari-251 ation (R 2 ad j in Fig. 2a). Regarding mean annual SST, T. sacculifer, G. siphonifera and P.
obliquiloculata increase in size significantly with linear increase of SST (Fig. 2b) while G.

253
truncatulinoides intraspecific shell size variation is significantly explained by a quadratic 254 function of SST. Shell size in the other five species did not covary significantly with SST.

255
No GLM with relative abundance as the sole explanatory variable was the best-supported 256 model (Table 2). Although relative abundance alone significantly explains shell size varia-257 tion within T. sacculifer and G. truncatulinoides (Fig. 2a), the best supported model for T. 258 sacculifer and G. truncatulinoides includes only SST, and adding abundance data has no 259 impact or decreases the amount of intraspecific size variation explained by the SST model 260 (R 2 ad j in Table 2). G. menardii's best supported model was the full GLM of both variables 261 (abundance and quadratic SST) plus their interaction term (Table 2), with ∆ AICc > 2 and 262 high model weight (Table 2). G. ruber and G. conglobatus show equal or similar weights 263 between the null and the relative abundance models; however, relative abundance does not 264 significantly explain shell size variation in these two species when tested alone (Fig. 2a). In 265 N. dutertrei and G. inflata, intraspecific variation was best explained by the null (intercept-266 only) model with R 2 ad j below 3% (Table 2). Visual inspection of the residual plots did not 267 reveal any obvious deviations from homoscedasticity, except for G. inflata (Fig. S4i).

268
The LMER shows that relative abundance and linear SST are both significant fixed ef-269 fects explaining PF population shell size variation (  was not the best supported model for most of the species analysed (Table 2).

295
Sea surface temperature explained more PF shell size variation than relative abundance ( Fig. 2b, Table 2). T. sacculifer, G. siphonifera and P. obliquiloculata, which are tropical-297 subtropical species (Kucera, 2007), showed a positive linear relationship between SST and 298 shell size. Moreover, the transitional G. truncatulinoides showed a quadratic relationship 299 between shell size and SST (Table 2). These results support the idea that PF species are 300 largest at their environmental temperature optimum (Hecht, 1976;Schmidt et al., 2004, 301 2006). However, the other analysed species (namely G. ruber, G. conglobatus, G. menardii,

302
N. dutertrei and G. inflata) showed neither a linear nor quadratic relationship between shell 303 size and SST (Fig. 2b, Table 2), contrary to the expectation of the ecological optimum hy-304 pothesis. The definition of optimal temperature range for a species is based on their relative 305 abundances in the marine sediments, with higher relative abundances indicating more opti-306 mal temperatures (Kucera, 2007). Thus, although SST could explain more intraspecific shell 307 size variation than local abundance, a positive monotonic relationship between shell size and 308 relative abundance of a species would still be expected under the ecological optimum hy-309 pothesis, regardless of the species' biogeography.

310
When increasing model power by analysing all the species together under a LMER 311 framework, relative abundance is a significant explanatory variable of PF intraspecific shell 312 size variation (Table 3). A linear positive relationship between shell size and SST is also 313 significant ( It might be that we did not find a strong relationship between size and abundance within 320 species because of the collector biases found in the NHMUK Henry Buckley Collection of 321 Planktonic Foraminifera (Fig. S2). Another concern regarding our analyses is that we used 322 relative abundance data from the ForCenS database (Siccha and Kucera, 2017) instead of the 323 abundance data estimated from the sediment samples used in the shell size data. As a result, 324 sometimes the ForCenS database yielded 0% of relative abundance of a species in the same 325 region that we had size data for the given species (Fig. 2a). Considering these two issues, we 326 assessed the robustness of our results by testing the same hypothesis on a more uniform, but 327 smaller, dataset. We re-sampled ten original sediment samples used by Buckley to amass 328 his collection (same samples used in the shell size bias analysis, Fig. 1a). We identified, 329 counted and measured the size of all PF individuals in each of the ten samples (Section 330 S3), minimising therefore any potential collector bias. Relative abundances of species were 331 calculated from each re-sampled assemblage itself, meaning that the same specimens were 332 used to extract abundance and size data. We then tested if population shell size could be 333 predicted by relative abundance in this re-sampled dataset using a linear-mixed effect model 334 with species as random effects. The re-sampled dataset included 20 species, summing 65 335 populations from the ten sites. The results showed no significant relationship between size 336 variation and relative abundance (Chi-square test, χ 2 = 2.18, P = 0.14, Table S4), supporting 337 our previous findings using the global Buckley Collection data and our statistical models.

338
Another source of bias in the Buckley Collection is that the samples come from different 339 expeditions using different sediment sampling strategies (Table S1). This source of bias is 340 inherent to this historical collection, as it includes samples from pioneering marine expedi-tions such as the HMS Challenger (1872−76) which lay on the foundation of oceanography 342 and ocean-floor sampling. In a previous study (Rillo et al., 2018), we showed that the PF 343 assemblages estimated from these historical samples are representative of Holocene assem-344 blages and can, therefore, be used in macroecological studies.

345
Ten of the 53 samples in our dataset come from sediments prone to dissolution (i.e., 346 waters deeper than 4000 meters for newly sedimented foraminifera, Berger and Piper 1972).

347
Dissolution may affect species size distributions, as smaller individuals are more prone to 348 dissolution (Kennett, 1976;Be and Hutson, 1977). We tested if water depth could explain 349 population shell size variation using a linear-mixed effects model with species as random 350 effects and found that water depth is not significantly related to PF size variation in our 351 dataset (Chi-square test, χ 2 = 1.83, P = 0.18, Table S5).

389
Another way of looking at species-level patterns is to plot the median population shell 390 size against median relative abundance. PF species showed a negative relationship between 391 size and relative abundance (Fig. 4). Abundant species such as G. ruber and G. inflata reach for a given amount of resources, have slower growth rates and obtain lower population den-397 sities than smaller organisms (Fenchel, 1974;Muller and Geller, 1993;Savage et al., 2004;398 Huete -Ortega et al., 2012). It remains to be tested whether smaller PF species have indeed 399 faster population growth rates.

400
At population-level, the mechanism that would lead to simultaneous increase of cell size would result in higher abundance in the sediment, but relative to other populations of the 414 same species, and not relative to the local assemblage (as the usual PF relative abundance 415 data). In the local assemblage, resource availability is the same for all co-occurring species.

416
As smaller species are generally more abundant in the sediment (Fig. 4), relative abundance 417 data regarding the local assemblage potentially blur within-species ecological patterns.

419
Our results caution against using the relative abundance of a species or SST to predict plank-420 tonic foraminifera intraspecific shell size variation. Regarding the understanding of PF ecol-421 ogy and evolution, maximum shell size might not indicate that a species is at its ecological 422 optimum, and/or the highest relative abundance of a species in the sediment might not co-423 incide with its ecological optimum. The low predictability of PF intraspecific size variation 424 found in our study also has implications for PF biomass estimation. If shell size is pre-        Table 3: Linear mixed-effects models ANOVA, using population size variation as response variable, species as random effects and fixed effects as sea surface temperature annual mean (sst linear effect, sst 2 quadratic effect), relative abundance of species (median within 300 km distance) (abund), plus the interaction between these two explanatory variables (sst:abund    choice also depended on the availability of bulk sediment samples in the OBD Collection.

677
Half of the amount available in the OBD containers was further split into two equal parts, 678 leaving an archive sample and a sample to be processed. The sample processing consisted 679 of weighing, wet washing over a 63µm sieve and drying in a 60 • C oven. The residues were 680 further dry sieved over a 150µm sieve and the coarser fraction was split with a microsplitter 681 as many times as needed to produce a representative aliquot containing around 300 PF shells.

682
All PF specimens in each of the nine final splits were identified by MCR and MK under a 683 stereomicroscope to species level, resulting in a total of 2,611 individuals belonging to 31 684 species (see also Rillo et al. 2018). This way, we calculated the relative abundance of each 685 species in each sample.

686
We then mounted species-specific slides from the re-sampled samples and extracted 687 shell size data in the same way as for the slides of the Buckley Collection (section 2.2).

688
Only species also present in the Buckley Collection samples were measured, resulting in 689 1824 specimens from 20 species (Table SI S3 Table SI S3).

696
The mean squared error was lowest for the 95th percentile ( Fig. SI S2), meaning that this 697 metric is the least biased measurement of the Buckley Collection when considering log-698 transformed shell area.  Using the re-sampled data described above, we tested whether relative abundance variation 701 significantly explains population shell size variation. Since the re-sampled data includes 702 only ten samples (Fig. 1a), there were not enough populations within each species to use species-specific GLM. Instead, we used linear mixed-effect models. The log-transformed 704 95th percentile of the population shell size distributions was modelled as the response vari-705 able, and the independent fixed effect was the local relative abundance (median within 300 706 km distance). Species were modelled as random effects, allowing for random intercepts and 707 slopes (i.e., the intercept and slope of the relationship between shell size and the relative 708 abundance may vary among species). We used the Likelihood Ratio Test (LRT) to com-709 pare the likelihood of the fixed effect. We calculated the LRT between the models with and 710 without the effect. Significance of each fixed effect was given through the LRT. Marginal 711 R 2 (R 2 m ), which is associated with the fixed effects, was calculated for each LMER model 712 (Barton, 2017).  We ran the same generalised linear models (GLM) analysis as described in section 2, but 725 using the species relative abundance retrieved from the nearest neighbouring sample of the 726 ForCenS database (instead of the median relative abundance of the samples within 300 km 727 distance). Although the order of the best-supported models changed for some species, the 728 models weights and the ∆ AICc are still consistent when compared to the model using the 729 median relative abundance within 300 km radius (Table 2). One example is G. ruber: here 730 the best supported model is the abundance one (Table SI S6) whereas for the median relative 731 abundance within 300 km the best supported model is the null model (abund, Table 2).

732
However, the ∆ AICc between these two models of G. ruber is close to zero (0.02) as well 733 as the difference in the models weights (0.01), consistent with Table 2.
734 Table S6: Model selection of the generalised linear models (with the Gamma logarithmic error function) testing if planktonic foraminifera shell size (represented by the 95th percentile of each population size distribution) can be predicted by sea surface temperature annual mean (sst) and relative abundance of species (nearest neighbouring ForCenS sample) (abund), plus the interaction between these two explanatory variables (sst:abund). Columns: explanatory variables, degrees of freedom, log-likelihood, Akaike Information Criterion corrected for small sample size, difference in AICc, model weight, adjusted R squared. response variable was the 95th percentile of the population size distribution whereas the 739 independent explanatory variables were the local relative abundances (median within 300 740 km distance) and mean annual sea surface temperature.