Extinctions in marine plankton preceded by stabilizing selection

Unless they adapt, populations facing persistent stress are threatened by extinction. Theoretically, populations facing stress can react by either disruption, increasing trait variation, or stabilisation, decreasing trait variation. In the short term, the more economical response is stabilisation, because it quickly transfers a large part of the population closer to a new ecological optimum. However, canalisation is deleterious in the face of persistently increasing stress because it reduces variability and thus decreases the ability to react to further change in stress. Understanding how natural populations react to intensifying stress reaching terminal levels is key to assessing their resilience to environmental change such as that caused by global warming. Because extinctions are hard to predict, observational data on the adaptive reaction of populations facing extinction are rare. In this study, we make use of the glacial salinity rise in the Red Sea as a natural experiment allowing us to analyse the reaction of planktonic Foraminifera to stress escalation in the geological past. We analyse morphological trait state and variance in two species across a salinity rise leading to their local extinction. One species reacted by stabilisation in shape and size, detectable several thousand years prior to extinction. The second species reacted by trait divergence, but each of the two divergent populations remains stable or reacted by further stabilisation. These observations indicate that the default reaction of the studied Foraminifera is stabilisation and that stress escalation did not lead to the local emergence of adapted forms. Inability to breach the global adaptive threshold would explain why communities of Foraminifera, and many other groups of marine plankton, reacted to Quaternary climate change by faithfully tracking their zonally shifting environments. It also means that populations of marine species adapted to response by migration, when exposed to stress outside of the adaptive range, will be at risk of extinction.


Introduction
Extinction events are remarkably difficult to study due to their unpredictability in recent environments (Moritz and Agudo 2013). While large efforts have been made to predict extinctions via population dynamics analyses (Tuljapurkar and Hecht Orzack 1980, Hanski 1998, Drake and Griffen 2010, Guzetta et al. 2017 this approach is still suffering from the large variation of population sizes (Ludwig 1999, Rand et al. 2017). Morphometrics could serve as a versatile alternative tool for predicting stress levels and identifying impending extinctions.
Shape and size of organisms have long been hypothesized to reflect the influence of environmental stress on the physiology of an organism during its lifetime (Furlow et al. 1997, Lens et al. 2002, Klingenberg 2003, Osterauer et al. 2010, Marschner et al. 2013. Therefore, a characterisation of shape and size and their variation should in principle allow an assessment of the severity of stress exposure. Under this assumption, stress exposure that leads to extinction, can be expected to leave a discernible imprint on morphology in pre-extinction populations (West-Eberhard 2003).
Developmental stability (compare Debat and David 2001) is an especially promising morphometric parameter for this purpose. The developmental stability over the ontogeny, as inferred from morphology, has been shown to be influenced by environmental stress (Furlow et al. 1997, Leung et al. 2000, Willmore et al. 2005, i.e. the sum of all biotic and abiotic factors that deviate from the ecological optimum of a species. Mechanistically, environmental stress influences population morphology by decreasing the fitness of specimens which show a high degree of developmental instability (Lens et al. 2002, Hendrickx et al. 2003. On the population level, this can lead to either stabilizing selection when a certain morphotype is the preferred survivor toward the stress (e.g. Van Valen 1965) or disruptive selection when higher variation better guarantees the survival of the population in a rapidly changing environment (e.g. Bull 1987). Both stabilizing and disruptive selection are detectable by assessing the phenotypic variation (i.e. the range of realized phenotypes, Schmalhauzen 1949) of the population, which can therefore be used as a measure for developmental stability.
While studies on developmental stability are frequent (e.g. Willmore et al. 2005, Osterauer et al. 2010, Sánchez-Chardi et al. 2013, O'Keefe et al. 2014, Smith et al. 2014 there is still much controversy as to what extent population morphology truly reflects the stress levels a community is exposed to (Hoffmann andWoods 2001, Santos et al. 2005). This is because it is difficult to study stress reactions over ecologically relevant time-scales, which cannot be simulated in the laboratory, or because the severity and ultimate outcome of the stress reaction is hard to predict.
In this regard, the sedimentary record offers a unique opportunity to study the effects of terminal stress levels (i.e. stress leading to local extinction), because here the outcome can be directly observed. This benefit comes at the cost that the sedimentary record always comprises a temporally integrated sequence, and that the environmental change that caused the local extinction can be difficult to reconstruct in some settings (compare Weinkauf et al. 2014). Shell bearing protists, such as planktonic Foraminifera, are perfect model systems for such studies, because they are preserved in great numbers in the fossil record (Schiebel 2002, Kučera 2007) and thus allow a robust analysis of their variation on the population level.
Understanding the morphological reaction of Foraminifera toward environmental stress could serve as a proxy for environmental reconstructions and further our understanding of evolvability of this organismal group. Past studies have shown that morphological deviations in benthic and planktonic Foraminifera can be caused by environmental forcing (Malmgren and Kennett 1976, Malmgren 1984, Baumfalk et al. 1987, Mary and Knappertsbusch 2013, Moller et al. 2013, Weinkauf et al. 2014, Knappertsbusch 2016, Brombacher et al. 2017. But since these studies quantified morphology in very different ways, the results obtained are scarce and controversial. In light of earlier results, it is reasonable to assume to see a morphological trend in the assemblage of planktonic Foraminifera associated with terminal stress levels. We hypothesize that a correlation between morphology and either environmental proxies or species abundance (as biotic indicator for the stress level the community is exposed to) exists. To test this hypothesis, we use Pleistocene sediments from the Red Sea, where several species of planktonic Foraminifera regionally disappeared from the fossil record (hereafter called local extinction) within aplanktonic zones resulting from environmental change (Fenton et al. 2000). Shortly before the onset of each aplanktonic zone, a distinct sequential local extinction pattern of different planktonic Foraminifera species can be observed, until virtually all planktonic Foraminifera are absent from the fossil record. In contrast to comparable cases where fossil material is used (e.g. Weinkauf et al. 2014, Brombacher et al. 2017), we are here in the unique situation that the extinction events can nearly exclusively be linked to salinity increase (Rohling et al. 2009), allowing a qualification of the stressor which acted on the assemblage.
The study will thus evaluate the value of morphology as proxy for stress and help to disentangle environmental forcing of morphology from pure stress-reactions. Should the environmental factor play a dominant role in organism morphology, regardless of the amount of stress it introduces, then morphology should mainly change as a function of salinity. Should the stress-level a community is exposed to be the dominant factor inducing morphological deviations, we assume to see a closer relationship between morphology and abundance patterns, under the assumption that abundance is a good indicator for ecologically optimal conditions.

Sample material
For the present study we used material from piston core Geo-TÜ KL 09 (19.804 • N, 38.103 • E, Suppl. S1) taken in the Red Sea during the RV Meteor cruise M5-2 (Nellen et al. 1996). We chose material from marine isotope stage 12 (MIS 12), specifically the range from 479.8 to 455.9 kyrs bp (age model by Rohling et al. 2009), to investigate planktonic Foraminifera morphology under terminal environmental stress. The interval covers an aplanktonic zone, which occurred during the Pleistocene in the Red Sea as a result of extremely high salinities (>49 on the psu scale) induced by a changed circulation pattern in the Red Sea basin (Fenton et al. 2000). The MIS 12 aplanktonic zone has been chosen, because it is the most prominent and longest one preserved in KL 09. The spatial resolution for our sampling varies between 0.5 cm and 2 cm such that a homogenous temporal resolution of 200-280 years/sample is achieved. Since salinity in the Red Sea is tightly coupled with the relative sea level (Thunell et al. 1988), and high resolution sea level reconstructions from the Red Sea exist (Rohling et al. 2009, Fig . 1), we can approximate past sea water salinity to test for its influence on foraminiferal shell morpho-logy. While MIS 12 represents a glacial period, it was shown that the sea surface temperature in the Red Sea area rarely dropped below 24 • C (Rostek et al. 1997), which is well within the temperature tolerance levels of all species common in the Red Sea (Prell et al. 1999a, Prell et al. 1999b, so that sea surface temperature should have only played a minor role on environmental stress levels. We investigated two abundant, symbiontbearing species of planktonic Foraminifera (Suppl. 1). Both species react sensitively to salinity changes, as shown by their consistently early position within the extinction sequence of the aplanktonic zones (Fenton et al. 2000). Orbulina universa d'Orbigny, 1839 is characterized by a trochospiral juvenile shell that is overgrown by a spherical terminal chamber. Trilobatus sacculifer (Brady, 1877) Spezzaferri et al., 2015 shows a trochospiral shell, in which the terminal chamber sometimes develops a saclike shape. Both species are surface dwellers, partly due to their symbionts, and thus occur in comparable environments (compare Rebotim et al. 2017, Schiebel andHemleben 2017). Differences in their reactions can therefore not be the result of the exposure to considerably different environmental forcing.

Figure 1|
Accumulation rates and relative abundances in relation to other species of planktonic Foraminifera of Orbulina universa and Trilobatus sacculifer during marine isotope stage 12 in piston core Geo-TÜ KL 09 from the Red Sea. The incidence of abnormal specimens in O. universa and sacculifer-morphotypes in T. sacculifer is indicated alongside its 95 % confidence interval (shaded area). The aplanktonic zone begins at approximately 465 kyrs BP. The two intervals of dropping abundance for O. universa (Intervals 1 and 2) and the three phases defined for T. sacculifer (Phases 1-3) based on relative abundances are indicated. The three-point moving average relative sea level in the Red Sea, as calculated by Rohling et al. (2009) based on δ 18 O values of Globigerinoides ruber (white), is also shown. Since salinity in the Red Sea is mainly driven by sea level (Thunell et al. 1988), this curve can be considered the inverse of sea surface salinity.

Sample preparation and data acquisition
Sediment material was prepared using standard procedures (compare Suppl. S1), and planktonic foraminiferal specimens of the >150 µm size fraction were analysed morphologically. Specimens from representative aliquots (split with a microsplitter) were picked with a needle and transferred onto glass slides, were they were fixed in position using permanent glue. We were striving for a sample size of at least 50 specimens per sample and for O. universa we always used the entire sample aliquot. Since the manual data extraction in T. sacculifer was much more time-consuming we only used a randomly chosen subsample of 50 specimens/sample from the aliquot in that species. Images for morphological analyses were taken with a Canon EOS 500 D digital mirror reflex camera attached to a Zeiss Stereo.V8 binocular microscope under constant magnification.
Morphological data for O. universa (2774 specimens) were semi-automatically extracted from high-contrast transmitted light images using the software FIJI (ImageJ v. 1.48s, Schindelin et al. 2012). We used the shell size (Feret diameter) and shell roundness (ratio between longest and shortest axis of a fitted ellipse) (Suppl. S1), as well as the incidence of the abnormal ecophenotypes 'Orbulina suturalis ' and 'Biorbulina bilobata' (Caron et al. 1987b). For T. sacculifer (2228 specimens), images were taken under reflected light with specimens oriented such that the apertural plane was lying horizontally, perpendicular to the direction of view (apertural standard view). In these images, 12 landmarks (Suppl. S1) were manually digitized in R v. 3.5.1 (R Core Team 2018). In 68 specimens, parts of the structures were not visible clearly enough to extract all landmarks and they were excluded from all analyses using the landmark data, leaving a dataset of 2160 specimens. Furthermore, the attribution to one of the three morphotypes of the species (trilobus, quadrilobatus, sacculifer; compare André et al. 2013) was recorded. Attribution to a morphotype was done by the same author (MFGW) for all specimens to avoid errors due to varying species concepts.

Morphometric data analysis
All statistical analyses have been performed in R v. 3.5.1 (R Core Team 2018). The normality of data distribution was tested with a Shapiro-Wilk test (Shapiro and Wilk 1965) and the homoscedasticity by a Fligner-Killeen test (Fligner and Killeen 1976) wherever necessary. Confidence intervals of morphological parameters were calculated via bootstrapping with the R-package 'boot' v. 1.3-20 (Davison and Hinkley 1997). Confidence intervals for morphotype occurrences were calculated using multinomial equations (Heslop et al. 2011) and for standard deviations of morphological parameters following equations in Sheskin (2011). In all cases of multiple testing (e.g. pairwise tests between more than two groups), p-values were corrected for the false discovery rate after Benjamini and Yekutieli (2001).
To investigate trends in the incidence of abnormal morphotypes of O. universa and sacculifer-morphotypes of T. sacculifer over time, we used generalized linear models (GLM, Nelder and Wedderburn 1972) on the binomial distribution with logit as link-function. The coefficient of determination was determined according to equations by Nagelkerke (1991) (Nagelkerke-R 2 ). For T. sacculifer we additionally used pairwise two-proportions z-tests to investigate the influence of stress levels on the incidence of sacculifer-morphotypes.
For O. universa, the extracted morphological parameters were subjected to traditional morphometric analyses. For T. sacculifer we used traditional morphometrics as well as geometric morphometric analytical methods as described in Claude (2008) and Zelditch et al. (2012). The landmark coordinates were fully Procrustes fitted using the R-package 'shapes' v. 1.2.4. The centroid sizes were used as size parameters for T. sacculifer shells. All size data were log e -transformed prior to analyses.
Traditional morphometric group-differences were investigated by a Kruskal-Wallis test (Kruskal and Wallis 1952), under circumstances followed by pairwise Mann-Whitney U tests (Mann and Whitney 1947). Bimodality was tested using Hartigan's dip test (Hartigan and Hartigan 1985) as implemented in the Rpackage 'diptest' v. 0.75-7, and the bimodality coefficient after Ellison (1987) was calculated with the R-package 'modes' v. 0.7.0. Trends in variation were analysed by calculating the coefficient of variation with the 95 % confidence interval after Vangel (1996). Trends over time were analysed using a Kendall-Theil robust line fitting model III linear regression (Kendall 1938, Theil 1950, Sen 1968). To test for correlations between size and roundness in O. universa we applied a Kendall rank-order correlation (Kendall 1938).
The shape of T. sacculifer specimens, defined as the Riemannian shape distance (Kendall 1984) of an individual landmark configuration from the grand mean shape, was investigated using geometric morphometrics. The variation of the T. sacculifer populations was approximated as variance of individual Riemannian shape distances within the population. Its confidence intervals and standard errors were derived by bootstrapping, and pairwise comparisons were performed by Student's t-tests with adjusted degrees of freedom (number of specimens in both groups minus two, Zelditch et al. 2012). Superimposed landmark data were analysed for differences between groups using non-parametric multivariate analysis of variance (NPMANOVA, Anderson 2001) on the Euclidean distances with 999 permutation as implemented in the R-package 'vegan' v. 2.5-2. Ensuing pairwise comparisons used the 'testmeanshapes' function of the R-package 'shapes' v. 1.2.4 with 999 permutations. Shape changes were visualized using thin-plate splines (Bookstein 1989), and canonical variates analysis (CVA Fisher 1936) from the R-package ' MASS' v. 7.3-50 (Venables and Ripley 2002) was used to analyse the shape changes between predefined groups.
We further tested all observed morphological trends against three potential models of phyletic evolution using the R-package 'paleoTS' v. 0.5-1 (Hunt 2006). This allows to distinguish between directional selection (general random walk), a directional pattern due to the accumulation of random change (unbiased random walk), and a system that does not change over time (stasis). The corrected Akaike information criterion (AIC c , Akaike 1974) in combination with Akaike weights (Wagenmakers and Farrell 2004) was used to decide, which model best describes the data.

Species abundances
Orbulina universa occurred at generally low abundances (5.6 % on average) that never exceeded 14.5 % of the total assemblage (Fig. 1). The species shows two abundance peaks in the studied interval, in case of the second event leading to local extinction. The first abundance peak occurred around 476 kyrs bp, and was followed by a rapid decline in abundance until 474 kyrs bp. Abundances then rose again until 473 kyrs bp, after which a second decline was observed, which culminated in the local extinction of the species between 466.6 and 467.4 kyrs bp. Based on the abundance, we could separate the O. universa population into two subsets (indicated as Intervals 1 and 2 respectively in Fig. 1) and treat them as a replication of the same general process.
Trilobatus sacculifer generally occurred at higher abundances of up to 38.4 % (on average 12.6 %) of the entire planktonic Foraminifera assemblage. From c.478.5 kyrs bp the abundance of the species decreased gradually until a local extinction between 465.8 and 465.0 kyrs bp, approximately 1000 yrs later than the comparable event in O. universa. Based on the species' relative abundance we defined three phases of population size leading to the extinction (Phases 1-3 in Fig. 1). Phase 1 with high abundances (24 %) spans between 478.2 and 475.1 kyrs bp, Phase 2 with medium abundances (10 %) between 474.9 and 472.6 kyrs bp, and Phase 3 with low abundances (4 %) between 472.4 and 465.7 kyrs bp. Those phases could be explicitly investigated for morphological developments within the community.

Morphotype abundances
In O. universa, two traditionally distinguished morphotypes ('B. bilobata' and 'O. suturalis') occurred with a mean abundance of 2.4 % ( Fig. 1). 'Biorbulina bilobata' occurs in marginally higher abundances (on average 1.4 %) than 'O. suturalis' (on average 0.9 %). The abundance of abnormal morphotypes decreases over time, as shown by a GLM over the entire time period (p < 0.001, Nagelkerke-R 2 = 0.21, Suppl. S1). When testing both intervals separately we observe that the decrease in abnormal morphotypes over time is significant for the second interval which leads to local extinction (p = 0.007, Nagelkerke-R 2 = 0.17) but not for the first interval (p = 0.773, Nagelkerke-R 2 = 0.01). In the T. sacculifer complex, we find higher abundances of the sacculifer-morphotype of up to 24 % broadly coinciding with Phase 1 of the abundance of the species (Fig. 1). During early Phase 2, this morphotype decreased rapidly in abundance and, although it never vanished completely within the limits of confidence, never comprised more than 10 % of the population from Phase 2 onwards. A constant decline over time is indicated by the GLM (p < 0.001, Nagelkerke-R 2 = 0.58). This is confirmed by a z-test: The abundance of the sacculifer-morphotype differs significantly (p < 0.001) between all three Phases, with decreasing mean incidence of that morphotype from Phase 1 (14.7 %) over Phase 2 (5.8 %) to Phase 3 (1.9 %).

Morphology of Orbulina universa across replicated drops in abundance
Morphological parameters of O. universa are presented in Fig. 2. The shell size of O. universa specimen indicates the successive establishment of two populations, as indicated by Hartigan's dip test and the coefficient of bimodality (Suppl. S1). Until 471.8 kyrs bp only one population with large shells (on average 399.2 µm) is present. After that, the size distribution reveals the existence of two populations with different sizes. One population shows an average size of 367.9 µm and can be considered a continuation of the population that was present before the split. While it differs in size from the larger population before the split on average (Mann-Whitney U test, p < 0.001) this can be explained by the gradual decrease in size of the larger population over time revealed by robust line fitting (R 2 = 0.08, p < 0.001, Fig. 2, Suppl. S1). The second population is significantly smaller (p < 0.001), only 154.7 µm on average, and only occurs with few specimens before the split from c.475 kyrs bp. This population shows the same trend of decreasing size toward the local extinction (R 2 = 0.13, p < 0.001, Suppl. S1). The large population became only marginally rarer, being more common before the split, but never disappeared from the sedimentary record until the local extinction. Noteworthy, the coefficient of variation of the larger population is significantly dropping at the splitting point from 0.21 to 0.18. However, the smaller group shows an overall strongly reduced coefficient of variation of only 0.11. The split occurs within Interval 2, and is unrelated to the replicated abundance pattern in O. universa.
Overall, shell size of the assemblage decreases toward the local extinction, while the population wide shell size variation remains constant due to the establishment of the small population ( Fig. 2e). When testing for a relationship between shell size and species abundance, a robust line fitting reveals a significant correlation between accumulation rates shell size for both the large (p = 0.023) and small (p < 0.001) populations (Suppl. S1). However, for the large population the coefficient of determination is negative (R 2 = −0.02), pointing toward a spurious effect, while being comparable to the time effect for the small population (R 2 = 0.11).
To investigate shell roundness in O. universa we excluded 'B. bilobata' and 'O. suturalis' specimens (59 specimens) from the analyses because both deviate significantly from the normal morphology of a terminal shell of that species. Shell roundness shows no indication of bimodality and is constantly decreasing over time according to a robust line fitting (R 2 = 0.06, p < 0.001, Suppl. S1). A secondary trend of shell roundness being correlated with species abundance exists, but it is much less pronounced (R 2 = 0.01, p < 0.001, Suppl. S1). Shell roundness shows an increase in variation contemporaneously with the decrease in shell roundness. Noteworthy, the major decrease in shell roundness and increase in its variation is unrelated with the replicated abundance drops, but occurs well within Interval 2 and no comparable trend is observable during Interval 1 (Fig. 2d,f).

Morphology of Trilobatus sacculifer during a long, continuous extinction event
Morphological parameters of T. sacculifer are presented in Fig. 3. The shell size shows a comparable pattern to the incidence of the sacculifer-morphotype, with larger shells during Phase 1, a rapid size decrease during Phase 2, and small shells during Phase 3. Comparing the values within the phases reveals a decrease in mean shell size from Phase 1 (411 µm) over Phase 2 (280 µm) to Phase 3 (237 µm). The differences in size are significant between all groups (p < 0.001 for a Kruskal-Wallis test, with all p < 0.001 in pairwise Mann-Whitney U tests). Moreover, the variation of shell size decreased   significantly during Phase 2 (Fig. 3a,b). This decrease in shell size and variation is nearly exclusively caused by the lack of large specimens after Phase 2, while the size of the smallest specimens remained rather constant during the entire interval (Fig. 3a).
Using geometric morphometrics allows a more detailed analysis of the shape of T. sacculifer shells (Fig. 3c,d). Since the sacculifermorphotype is a considerably derived morphology, we left it out of the dataset for most analyses (n = 136 specimens), but results including the sacculifer-morphotype specimens are presented in Suppl. S1.
The shape of T. sacculifer specimens between abundance phases differs significantly (NPMAN-OVA, p < 0.001). A post-hoc pairwise comparison of shape differences reveals, that the shape of T. sacculifer specimens differs between all three phases at p = 0.002, and the shape changed especially during the last c.3000 yrs before local extinction (Fig. 3c). Shape variation shows no trend over time according to a robust line fitting (R 2 = −0.002, p = 0.687), but would be significant when including the sacculifer-morphotype (R 2 = 0.35, p < 0.001). This implies that a large part of the shape variation is linked to the disappearance of the sacculifer-morphotype in the population (Suppl. S1). When comparing integrated values across the three phases defined by abundance patterns, the data reveal that variation is highest in Phase 2, mediocre in Phase 1, and lowest in Phase 3 (Fig. 3d). A t-test reveals that differences between Phases 1 and 2 are insignificant (p = 0.944), but specimens from Phase 3 differ significantly from specimens from both Phase 1 (p = 0.010) and Phase 2 (p < 0.001). This implies a significant reduction of phenotypic plasticity in the community when exposed to higher stress levels and impeding local extinction.
When analysing the shape change from phase to phase (Fig. 4) it is revealed that from Phase 1 to Phase 2 the terminal chamber became flatter and the aperture was inflated. At the same time the lower part of the shell containing the older chambers became more voluminous, so that there is a general trend towards a smaller terminal chamber in relation to the older chambers. This trend is reversed in Phase 3, with the terminal chamber becoming more inflated and the aperture flattening again. This led to a strong increase in the size of the terminal chamber in regard to the older shell.

Error analysis
Morphometric data are especially prone to errors, because many things (like the orientation of specimens and digitization of landmarks) are done manually. We therefore analysed our data for the influence of potential error sources, following methods described by Yezerinac et al. (1992). We could not detect any errors of severe size in our data and show a full error discussion in Suppl. S1.

Stabilisation and disruption within populations of planktonic Foraminifera
The concept of phenotypic plasticity (Schmalhauzen 1949) describes the observed morphological variation among specimens of a population, and is thus contrasted to variability, which is the potential of a population to vary within the borders of the genetic encoding (Wagner and Altenberg 1996). This is especially important in planktonic Foraminifera, where the existence of (pseudo-)cryptic species can lead to an underestimation of genetic diversity and thus variability of the population. In T. sacculifer we can rule out that possibility, because despite its large morphospace it only contains one biological species with large variation (André et al. 2013).
In O. universa, the morphospecies encompasses several biospecies (de Vargas et al. 1999), but at present only one genotype is known to occur in the Red Sea (Darling and Wade 2008). However, the multimodality in shell size of O. universa we observe in the upper part of the profile points towards the existence of more than one population. While at the moment we do not know of any instance were pseudo-cryptic species in (c-d) Shell shape (Riemannian shape distance from the mean shape) excluding sacculifer-morphotypes (compare Suppl. S1). (c) Raw values (grey dots) are plotted alongside the sample mean and standard deviation (solid lines). (d) The shape variation is significantly lower in Phase 3 than it was before.
O. universa can be differentiated by shell size (Morard et al. 2009), we must consider it a possibility that two O. universa species where prevalent in the Red Sea during MIS 12. Importantly, by controlling for cryptic diversity we could assume that the populations investigated in our study share a common variability which does not change dramatically over the time of this natural experiment. This implies that the morphological changes we observe in the studied interval are results of a changing variation of a homogeneous community, and therefore the effect of environmental forcing associated with the onset of the aplanktonic zone.
In O. universa we observe bimodality in shell size and significant increase in variation of shell roundness, indicating disruptive selection and decanalization in the population in the upper part of the profile, coinciding with increasing salinity. Especially under random mating models, disruptive selection can introduce an increase in variation, without inducing bimodality (Doebeli 1996), and the degree of disruption can vary in different traits (Wagner and Altenberg

Figure 4|
Thin plate splines of the deformation of Trilobatus sacculifer specimens from marine isotope stage 12 in the Red Sea between three phases defined by relative abundances (compare Fig. 1).

1996, Klingenberg 2008).
A decrease in size and mean roundness when the species approached its local extinction in the aplanktonic zone is apparent. Both parameters show signs of an unbiased random walk pattern (Table 1), implying that the observed morphological change is the result of accumulation of random changes rather than directional selection. When performing the analysis separately for both size groups in O. universa, it is revealed that the trend of unbiased random walk is inherent to the larger population, while the smaller population shows stasis. A reduction in shell size has been proposed to indicate suboptimal environmental conditions (Ortiz et al. 1995, Schmidt et al. 2004) but also a buoyancy adaptation to increasing salinity in order to be able to stay at an optimum depth in denser water (Haenel 1987). Since the strong salinity increase in the Red Sea associated with the onset of the aplanktonic zone shifted the local habitat away from the optimum requirements of the incumbent species, it is therefore hard to say whether the trends we observe are the result of abiotic or biotic forcing. However, if the observed morphology change would be a purely biotic stress response, assuming that abundance is a useful proxy for the suitability of the environment (Drake and Griffen 2010), we would principally expect to see a comparable development in Interval 1 as in Interval 2. In contrast to that assumption, neither shell size nor shell roundness show any signs of a deviation associated with the first abundance drop at the end of Interval 1. We further observe that shell size and shell roundness are significantly different between both intervals at p < 0.001, which is not what one would expect if they would show a comparable internal pattern.
The observed pattern would be consistent with two possible end-member scenarios: (1) a monospecific O. universa population present at the time in the Red Sea shows signs of disruptive selection as a result of exposure to a suboptimal environment (Bull 1987), or (2) a new population comprising a different biospecies with different morphology is introduced into the Red Sea once environmental conditions become more suitable.
To investigate those scenarios, we applied a Kendall rank-order correlation between individual shell size and shell roundness and found that they are significantly correlated (ρ = −0.321, p < 0.001, Fig. 5). This implies that smaller individuals also tend to have less round shells. To test this, we divided the population into two subsets at the local minimum of the size distribution (238.6 µm) and tested the shell roundness in both subgroups against each other with a Mann-Whitney U test. The test confirms that the larger population produced shells that are on average significantly (p < 0.001) rounder (mean = 1.03) than the smaller subpopulation (mean = 1.07). The fact that the maximum shell roundness is practically 1.00 in both populations shows, that the observed trend is neither cause by a biological inability of shells to grow very round below a certain size nor the result of a lower precision of shell roundness determination in smaller shells. Comparing the 95 % confidence interval of the coefficient of variation shows a significantly higher variation of shell roundness in the population with smaller shells (0.036-0.039) than in the population with larger shells (0.020-0.022). Since the smaller population becomes increasingly dominant toward the local extinction, this can explain the observed increase in shell roundness variation toward the local extinction (Fig. 2). Noteworthy, shell roundness variation did not meaningfully increase over time in either population, which shows the variation as inherent to the respective population (Suppl. S1). It is also possible that shell roundness and shell size are biologically integrated in O. universa, so that a change in one parameter necessitates a certain change in the other value (Klingenberg 2008). This would be an alternative explanation for the observed correlation between individual shell size and shell round-ness, but the very low covariance between both parameters (−0.008) together with the fact that small specimens realize the complete range of possible deviations from sphericity principally argues against that hypothesis.
Our observations do not allow us to decide for one of the two potential hypotheses. The Red Sea was possibly invaded by a population with smaller, more variably shaped shells from approximately 471.8 kyrs bp onwards, which increasingly established its presence at the expense of the incumbent population. Should this be true, we would find here the first example where different O. universa biospecies could be distinguished on the basis of relatively easily obtainable morphological characters (compare Morard et al. 2009) and also the first example of more than one O. universa biospecies occurring in the Red Sea (compare Darling and Wade 2008). Furthermore, it would provide evidence for different ecological preferences among those biospecies, which facilitates competitive exclu-sion due to increasing stress levels (Sultan and Stearns 2005). However, we would argue that this explanation is unlikely on the basis of other assumptions. (1) The potentially invading species would have had to be of Indian Ocean origin. Within an environmental setting of increasing salinity, it is hard to perceive that any species from the open-marine Indian Ocean would have a selective advantage over a native species from the Red Sea, that should be better adapted to high salinities (Sofianos and Johns 2003). (2) A trend of decreasing shell size with increasing salinity in O. universa was already observed by Haenel (1987), and was there attributed to buoyancy requirements, although Spero (1988) saw a closer relation to nutrient availability due to the energy requirements to build larger shells. This makes it reasonable that the observed size change is indeed an adaptive response by part of the population. (3) The splitting in two size populations as well as the increased variation in roundness in the smaller population both imply disruptive selection (Badyaev 2005). We thus favour the explanation that the incumbent O. universa population underwent disruptive selection as a result of environmental stress, associated with drastically increasing salinity levels in the Red Sea. This would replicate results obtained on a Mediterranean population that was exposed to higher stress levels in relationship with the onset of sapropel deposition (Weinkauf et al. 2014), and would lend valuable evidence to the assumption, that certain morphological changes can be universally interpreted as indicator of environmental stress, at least when the same species are considered. The small subspecies would in that scenario be the better adapted result of accumulated random change (Hunt 2007), and would fall into a state of evolutionary stasis once optimal adaptation is reached. This conforms with the proven development of stasis in extreme environments (Parsons 1994), while both populations can coexist as long as an equilibrium point exists along the Lotka-Volterra isoclines (Sultan and Stearns 2005). Divergence without speciation is often observed in nature, and new character traits that are encompassed in the variability of the population can be fixed, once a diverging subpopulation occupies a new niche (West-Eberhard 2003). In the present case, the size difference in the two O. universa populations points toward a different in the depth habitat due to differential buoyancy (Haenel 1987).

Figure 5|
Correlation between shell size (as Feret diameter) and shell roundness in Orbulina universa during marine isotope stage 12 in the Red Sea. The population shows clear bimodality in shell size. Smaller specimens show a significantly higher variation of roundness than larger specimens.
While thus disruptive selective pressure in O. universa may have led to the development of two populations, each individual population shows signs of stabilizing selection in its shell size parameter. Trilobatus sacculifer draws a very similar picture in the form of decreasing variation in its shell size and shape when the species approached local extinction. In shell size this decrease complies with the unbiased random walk model, as it does for the shell shape (both with and without sacculifer-morphotypes), expressed as distance to the mean shape.
While the size decrease may itself be a signal for decreasing environmental suitability for the species (Ortiz et al. 1995, Schmidt et al. 2004, the decrease in variation is a clear signal for stabilizing selection inducing microenvironmental canalization (Waddington 1942, Schmalhauzen 1949. It is interesting, that stabilization occurs toward Phase 3, which must have been clearly environmentally suboptimal for T. sacculifer (drop in abundance, drastic increase in salinity). This could indicate either a stabilization of the regional environment, allowing selection for an optimal trait (Van Valen 1965), or a rapidly changing environment enforcing fluctuating selection that benefits a stable phenotype in the long run (Kawecki 2000, Pélabon et al. 2010. The slight increase of disparity in Phase 2 may be an indicator for the second explanation, and indicate a Baldwin effect (Baldwin 1902), i.e. a phase of increased plasticity that allows adaptation to a new environment and evolutionary fixation afterwards. The reduction in the abundance of the sacculifer-morphotype after Phase 1 could be another manifestation of the stabilizing selection, but it could also result from an inability of T. sacculifer shells to build an asymmetric chamber below a certain chamber size threshold. In the latter case, the drop in the abundance of the sacculifer-morphotype would result from the shell size decrease and not necessarily be a signal for stabilizing selection.
We observe certain morphological trends in T. sacculifer over time, as shown in Fig. 6, which can be interpreted as adaptive for a changing environment that gradually becomes less suitable for the population. During Phase 1, shells show more inflated and larger terminal chambers. This is reminiscent of the sacculifermorphotype, which was left out of the CVA due to its strongly derived morphology. This phenotype is often assumed to be more abundant under lower stress levels (Caron et al. 1987a) but has already earlier been shown to be correlated with larger shell sizes (Hemleben et al. 1987). During Phase 2, coinciding with a first strong drop in sea level and thus salinity increase ( Fig. 1), the abundance of the sacculifermorphotype and the shell size decreased, while the size of the terminal chamber decreased as well, indicating a trend towards Kummerforms that are often associated with unfavourable environmental conditions. Finally, in Phase 3 the shell size decreased further while the terminal chamber became larger again in relative terms. The latter trend could be necessitated by the small shell size, so that the terminal chamber had to become relatively larger again to provide enough space for gametogenesis. This phase is also characterized by a canalization peak, probably induced by an environment that was so unstable and unfavourable for the T. sacculifer population that any deviation from a very narrow morphotype would drastically decrease survival rates and fitness. The entire investigated timespan can possibly be interpreted as follows: During Phase 1 an incumbent population of T. sacculifer was well adapted to the high but not excessive salinity levels in the Red Sea during MIS 12. During Phase 2, salinity levels started to rise, and induced a selection that via accumulation of random change led to a modification of the phenotype (Hunt 2007). In this phase, a Baldwin effect set in, increasing the variation of shape and thus the buying the population time until further adaptation could settle in (Baldwin 1902). During Phase 3, the adaptation was as complete as it could possibly become, and the successful morphologies became fixed via canalization (Schmalhauzen 1949, Gibson and Wagner 2000, West-Eberhard 2003. It is interesting to note that the trend towards smaller shells that was observed in O. universa is also present in T. sacculifer. This indicates that probably, buoyancy problems due to the increasing water salinity were responsible for this adaptive trend.

Planktonic Foraminifera morphology as environmental proxy
Morphology has often been shown to be influenced by the environment (e.g. Hendrickx et al. 2003, Willmore et al. 2005, Weinkauf et al. 2014 but recently this hypothesis has been challenged. For example, it has been theorized that morphological change may be under circumstances largely decoupled from the environment, because the observed changes in morphology have no impact on fitness (Hoffmann andWoods 2001, Santos et al. 2005) and morphology can therefore not serve as basis for selection. Alternatively, it was hypothesized that several developmental traits are buffered by the same genes and are therefore biologically integrated (Breuker et al. 2006, Klingenberg 2008, Breno et al. 2011, Klingenberg et al. 2012). This would imply that an observed change in any morphological trait could be the by-product of a change in another (superficially unrelated) parameter that is the actual basis for selection. With our data, however, it seems that the environmental change left a strong imprint in the morphology of the foraminiferal assemblage, which can at least partly be explained by selection for particular traits and thus adaptive processes working on the phenotype.

Figure 6|
Canonical variates analysis (CVA) of the shape of Trilobatus sacculifer (excluding sacculifer-morphotype) from marine isotope stage 12 in the Red Sea. Points indicate specimens in the three pre-defined abundance phases (compare Fig. 1), ellipses of the same colour indicate the 95 % confidence interval of the standard deviation on the centroid of the respective group, black silhouettes depict the morphology at the extremal points of the canonical variate (CV) 1 and 2 axes, respectively. The CVA is significant (T 2 = 0.631, F approx = 31.585, p < 0.001) with 82.2 % of the variance explained by CV 1 and the remaining 17.8 % by CV 2, and a correct classification rate of the discriminant function after leave-one-out cross validation of 69.27 %.
Qualitatively, it is obvious that the abundance decrease in T. sacculifer began much earlier (around 477 kyrs bp) than in O. universa (approximately 473 kyrs bp). However, while the abundance decreased steadily in T. sacculifer, it is more variable in O. universa, with one abundance lowpoint before the study period and another one at the end of Interval 1. This indicates that abundance in O. universa might not be strictly coupled to salinity but also influenced by other factors such as competition with other species. It is noteworthy in this regard, that the emergence of a second population of smaller specimens in O. universa coincided with the renewed increase in abundance. This could point toward the fact, that the new population was better adapted to the changing environment and had better capabilities to compete with other foraminiferal species.
In T. sacculifer, the abundance is likely to mirror salinity changes and the species became extinct c.1000 yrs later than O. universa, but shows first morphological reactions already during early Phase 2, i.e. c.2000 yrs earlier than the deviations in shell size and roundness in O. universa occur. Nevertheless, this offset is within the limits of what could be expected if both species have different salinity threshold tolerances and the morphological change was triggered by salinity increase. Conversely, under the assumption that abundance is a good indicator for stress, it is unlikely that stress and individual fitness alone influenced the population morphology, because during the first abundance drop in O. universa at the end of Interval 1 no morphological deviation could be observed. Orbulina universa seems to be able to resist salinity changes rather long, but after a certain threshold is reached shows heavy reactions, as evidenced by the later onset of morphological change and earlier local extinction. In contrast, Trilobatus sacculifer seems to react earlier to the environmental stress levels, but due to its higher evolvability being better able to adapt to the changing environment and survive longer under equally unfavourable conditions. We thus hypothesize that the reactions seen in both species are a morphological reaction to stress induced by salinity changes, and that the abiotic and biotic component are so closely related that they cannot be disentangled. Conversely, species abundance seems to be a less valuable indicator of environmental stress, and might be more reliable in some species than in others.
While the O. universa community reacted towards the environmental stress by decanalization, we see convincing evidence for canalization in T. sacculifer within several traits, which are either all selected for or which are integrated and cannot vary independently (Klingenberg 2008, Brombacher et al. 2017. While results for O. universa shell size and roundness are in line with observations from a Mediterranean sapropel (Weinkauf et al. 2014), the conservative bet-hedging observed there and evidenced by increasing abundances of abnormal morphotypes shortly before local extinctions could not be replicated here. This, however, can result from the different stress pattern in the Red Sea in comparison to the study by Weinkauf et al. (2014). Bet-hedging is hypothesized to be especially beneficial for a population that is exposed to a very variable environment which changes unpredictably (Philippi andSeger 1989, Beaumont et al. 2009). While this assumption is reasonable for the study presented in Weinkauf et al. (2014), our study here was explicitly chosen to investigate the effect of one parameter that changes gradually. There is no reason to assume that the salinity in the Red Sea at that time changed unpredictably (Thunell et al. 1988), rather than increasing continuously. Since the environment thus developed toward an unfavourable but not unstable state, there is no reason to expect bet-hedging in the community.
It is, however, interesting that a reduction in shell roundness in O. universa was already observed by Weinkauf et al. (2014), but was there (in contrast to this study) associated with a decrease in salinity. This lends evidence to the hypothesis that some morphological parameters are indeed a neutral proxy for environmental stress, widely independent of the exact type of stressor or environmental forcing causing that stress. Another interesting observation is, that all observed morphological trends follow patterns of either unbiased random walk or stasis, but we could not observe general random walk patterns in any parameter. This implies that directional evolution plays no role in the selective response of the planktonic Fo-raminifera investigated here, at least on the timescales and environmental changes covered by our study. Both species showed different adaptive responses, and both were ultimately not successful. This does not invalidate the adaptive nature of these processes, however, but merely means that salinity levels became that high that the epigenetic landscape (within which adaptation of a species is theoretically possible) of the Foraminifera could not accommodate them (Waddington 1974). Similar patterns, where adaptive trends in planktonic Foraminifera resulted ultimately in extinction due to an excessive amplitude of the environmental change, have been observed before by Brombacher et al. (2017).
While we thus have to admit that the direct correlation between environmental and biotic parameters and morphological deviation are hard to grasp and need further research, we could show that both species showed clear morphological trends with increasing stress levels. It could not be clearly shown whether salinity changes as abiotic parameter or increasing stress levels were the main inductor for the observed morphological deviations, because both parameters are still too strictly coupled. The fact that morphological deviations in both species occurred relatively synchronous, however, makes it reasonable that salinity changes were a major influential parameter. Especially when assuming that O. universa showed no morphological deviation at the end of Interval 1, when the abundance already dropped significantly a first time. The results strongly suggest that morphology is a versatile tool to reconstruct past stress levels on the foraminiferal assemblage, and can lead to better estimates about the timing of stress induction than species abundance (and thus population dynamics) can do.