Test-Size Evolution of the planktonic Foraminifera Globorotalia menardii in the Eastern Tropical Atlantic since the Late Miocene

The mean test size of planktonic foraminifera (PF) is known to have increased especially during the last 12 Ma, probably in terms of an adaptive response to an intensification of the surface-water stratification. On geologically short timescales, the test size in PF is related to environmental conditions. In an optimal species-specific environment, individuals 10 exhibit a greater maximum and average test size, while the size decreases the more unfavourable the environment becomes. An interesting case was observed in the late Neogene and Quaternary size evolution of Globorotalia menardii, which seems to be too extreme to be only explained by changes in environmental conditions. In the western tropical Atlantic Ocean (WTAO) and the Caribbean Sea, the test size more than doubles from 2.6 Ma to 1.95 Ma and 1.7 Ma, respectively, following an almost uninterrupted and successive phase of test size decrease from 4 Ma. Two hypotheses have been suggested to explain the sudden 15 occurrence of a giant G. menardii form: it was triggered by either (1) a punctuated, regional evolutionary event or (2) the immigration of specimens from the Indian Ocean via the Agulhas Leakage. Morphometric measurements of tests from sediment samples of the Ocean Drilling Program (ODP) Leg 108 Hole 667A in the eastern tropical Atlantic Ocean (ETAO), show that the giant type already appears 0.1 Ma earlier at this location than in the WTAO, which indicates that the extreme size increase in the early Pleistocene was a tropical-Atlantic-Ocean-wide event. A 20 coinciding change in the predominant coiling direction suggests that probably a new morphotype occurred. If the giant size and the uniform change in the predominant coiling direction are an indicator for this new type, the form already occurred in the eastern tropical Pacific Ocean at the Pliocene/Pleistocene boundary at 2.58 Ma. This finding supports the Agulhas Leakage hypothesis. However, the hypothesis of a regional, punctuated evolutionary event cannot be dismissed due to missing data from the Indian Ocean. 25 This paper presents the AMOC/thermocline hypothesis, which not only suggests an alternative explanation for the sudden testsize increase in the early Pleistocene, but also for the test size evolution within the whole tropical Atlantic Ocean and the Caribbean Sea for the last 8 Ma. The test-size evolution shows a similar trend with indicators for changes in the Atlantic Meridional Overturning Circulation (AMOC) strength. The mechanism behind that might be that changes in the AMOC strength have a major influence on the thermal stratification of the upper water column, which is known to be the habitat of 30 G. menardii. Thore Friesenhagen1,2 https://doi.org/10.5194/bg-2021-67 Preprint. Discussion started: 15 March 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
While short-term changes in the test size of planktonic foraminifera (PF) are thought to be related to changes in environmental conditions (e.g. Poore, 1981;Keller, 1985;Ravelo et al., 1990;Wolff et al., 1999;Chaisson and Ravelo, 1997; (1) This area is located within tropical waters, which are known to be the habitat for G. menardii (Caley et al., 2012). Surface 75 sediments show that G. menardii has a high Holocene occurrence at this site. Throughout the studied interval, sediments are used to contain an adequate number of G. menardii and related forms (Manivit, 1989). This location is outside or only peripheral reached by the NW African upwelling system (Fig. 1), so that it was marginally affected by this upwelling system for the investigated time interval of the last 8 Ma (Weaver and Raymo, 1989). Thus, it is supposed to show a relatively and long-term water-column stability on the geological timescale.
(2) This area is within the range of water masses which are 80 affected by the Agulhas Leakage (Biastoch et al., 2009;Rühs et al., 2013), so that biota originating in the Indian Ocean is transported by currents up to this location. (3) The preservation of the fossils is good to moderate (Manivit, 1989), which is partly attributed to a sediment deposition depth (present: 3529.3 m water depth below sea level) above the carbonate compensation depth. (4) For the studied interval from 8 Ma until present, the sedimentation has most likely been continuous.
The sediment sequence is only disturbed by a small slump (Shipboard Scientific Party, 1988), which was avoided for sampling. 85

Sample Selection
The samples were chosen from interglacial periods with a similar age as the investigated samples of the studies of Knappertsbusch (2007; (Table 1). The working hypothesis presumes G. menardii to reach its maximum test size during interglacials, inferred from the observation of overall decrease of population size or even complete absence during glacial intervals in the Atlantic Ocean (Ericson and Wollin, 1956;Sexton and Norris, 2011) and references therein; (Portilho-Ramos 90 et al., 2014).
Due to the lack of stable isotopic data for this site, the Age-Depth plot uses biostratigraphic data of PF and nannoplankton (Weaver and Raymo, 1989) as well as magnetostratigraphic events (Shipboard Scientific Party, 1988;Fig. 2). The Age Depth Plot program by Lazarus (1992) was used to manually draw a line of correlation (loc) through recognised bio-and magnetostratigraphic events (Fig. 2). Using the loc's control points numerical ages were computed by linear interpolation with 95 the help of the Age Maker (Lazarus, 1992) (NEPTUNE Age Model; see supplementary materials file "667A.loc95.txt"). The Age-Depth plot is based on published core-depth information from Hole 667A (Shipboard Scientific Party, 1988) and biostratigraphic occurrences of first and last occurrence dates of nannofossils, planktonic and benthonic foraminifera and magnetostratigraphic polarity reversals given in the initial reports and scientific results of that Leg. The time chronology of Berggren et al. (1995) was applied to allow direct comparison to previous studies of Knappertsbusch (2007;. 100 Table 1: Studied samples, their depths in meter below seafloor (mbsf) and age (Ma) of Hole 667A, following the age-depth plot of

Sample Preparation and Parameter Measurement
The procedure for the treatment of the samples follows that of Knappertsbusch (2016).
Approximately 2-3 cm³ bulk sediment per sample were dried at 40°C over night and weighted. In a following step, the samples were gently boiled with water, containing soda as an additive, and wet-sieved with a 63 µm net. The fraction <63 µm was 110 decanted, dried and preserved. The >63 µm fraction was dried at 40°C for 24 h and weighted afterwards. A microsplitter was used to split the >63 µm fraction until at least 200 menardiform specimens could be picked from the sample. This number of specimens was judged to be a reasonable compromise between efforts for picking and manual mounting, imaging, analytical steps and statistical and the limited amount of time for this project. The specimens were mounted on standard faunal Plummer  Weaver et al. (1989) and from Weaver and Raymo (1989). Magnetostratigraphic data were taken from Shipboard Scientific Party (1988). The vertical bars within the symbols illustrate the depth range in which this event took place. The data for the palaeomagnetic reversals below the x-axis are taken from Berggren et al. (1995). The red bars on the right side indicate cores and core recovery.  (Malmgren et al., 1983). Intact specimens 120 showing a menardiform morphology were picked from the sample splits. They include the whole G. menardii lineage as well as members of the G. tumida lineage. In total, 4482 G. menardii, 764 G. limbata and 228 G. multicamerata specimens were picked from samples at 33 stratigraphic levels back to 8 Ma (Table 1). All study 125 material is stored in the collections of the Natural History Museum Basel.
Digital images of the menardiforms were collected with the Automated Measurement System for Shell Morphology (AMOR), software version 3.28. This system automatically orientates and 130 photographs tests in keel view to achieve orientation for outline analysis (Knappertsbusch et al., 2009). The free software ImageJ 1.52i of the National Institute of Health was used to clean and preprocess images for outline coordinate extraction. Processing steps include removal of adhering particles, smoothing, enhancement of 135 contrast, binarization, closing of single pixel embayments before storing the processed pictures as 640 x 480 pixel and 8 bit greylevel Tiff files. Adapted MorphCol software programmed in Fortran 77 from Absoft by Knappertsbusch (2007; were used (Appendix Fig. A1) to extract cartesian outline coordinates and to 140 derive morphometric measurements. These applications were converted to Fortran 95 versions and adapted for usage on Windows operating systems. The adapted MorphCol programs and codes are deposited at the PANGAEA data repository.
These programs considerably accelerate the process of measuring 145 several different morphometric parameters from the images. Derived parameters include the spiral height ( X) and the axial length ( Y), their ratio (R = X / Y), the area of the specimen in keel view (Ar), the convexities of the spiral (A) and the umbilical (B) side, the ratio of the convexities (RA/B), the upper (φ1) and the lower (φ2) keel angles, the angle at the apex (φ3) as well as the radii of the osculating circles in the upper (Rup) and the lower (Rlo) keel region (see Fig 3). This study focuses on the test-size parameter X, Y and Ar. In order to compare 150 specimens with dextral and sinistral coiling, dextral specimens were vertically mirrored using the adapted "DexFlip_win" program, modified from Knappertsbusch (2016).

Globorotalia menardii Lineage -Species Discrimination
Globorotalia menardii is discriminated from its extinct descendants Globorotalia limbata and Globorotalia multicamerata by 155 the number of chambers in the last whorl (Fig. 4). Menardiform specimens with six or less chambers were determined as G. menardii. Specimens at 2.39 Ma (Wade et al., 2011). Globorotalia multicamerata has more than seven chambers in its last whorl and became extinct in the early Pleistocene at 3.09 Ma (Berggren et al., 1995). Images for figure 4 were made with a Keyence VHX-6000 digital microscope.
The identification of menardiform specimens is based on illustrations in Kennett and Srinivasan (1983), Bolli et al. (1985) and 160 comparison with the reference collection to "49 Cenozoic planktonic foraminiferal zones and subzones prepared by Bolli in 1985Bolli in -1987 which is deposited at the Natural History εuseum Basel. Knappertsbusch (2016) refers to the disappearance of G. limbata as a possible pseudoextinction because of the occurrence of singular specimens of menardiforms with seven chambers in the last whorl after 2.39 Ma. with seven chambers are accounted as G. limbata, a form which became extinct during the early Pleistocene. 165

Univariate and Contoured Frequency Diagrams
Statistical analysis and univariate parameter-versus-age plots were prepared with RStudio (V. 3.5.3; RStudio Team, 2020), using the packages psych (Revelle, 2018), readxl (Wickham and Bryan, 2019), ggplot2 (Wickham, 2016), pacman (Rinker and Kurkiewicz, 2018) and rio (Chan et al., 2018). For the generation of contour frequency diagrams (CFD) the commercial software Origin 2018 and Origin 2019 by OriginLab Corporation was used. CFDs per species help to detect shifts in the 170 dominant test size of populations throughout time. The same method was applied in Knappertsbusch (2007; and enables a direct comparison of evolutionary change in Hole 667A with previous studies. Emergence and divergence of new frequency peaks between subsequent samples may help to empirically identify signs of cladogenetic splitting or anagenetic evolution in the lineage of G. menardii-G. limbata-G. multicamerata. The CFDs were constructed from so-called "gridded files": Basically, these gridded files were obtained by plotting X versus Y, superposing a grid with grid-cell sizes of ΔX = 50 µm 175 and ΔY = 100 µm (see Knappertsbusch 2007; and then counting the number of specimens per grid cell. This gridding procedure was performed with program "Grid2.2_win" (adapted εorphCol software by Knappertsbsuch, 2007;, and the result was a two-dimensional matrix of absolute frequencies of specimens per grid-cell. No smoothing of frequencies was applied, because experiments revealed an increasing loss of frequency variation with increasing size of bin-width. However, in contrast to Knappertsbusch and Mary (2012) and Knappertsbusch (2016), local absolute specimen frequencies were used 180 throughout instead of relative frequencies.
Different contour intervals were used for the CFDs, because the number of G. menardii specimens per sample varies from one (667A-5H-2, 105-106 cm) to 273 (667A-4H-3, 120-121 cm). This approach increases the legibility of the single CFDs, because setting a high contour interval in a sample with few specimens would have levelled out the CFD. Conversely, choosing a low contour interval would lead to exaggerated contour line densities in CFDs when the number of specimens is high.

Volume Density Diagrams
Volume density diagrams (VDD) were made with the commercial software Voxler 4 by Golden Software. This method was shown to be useful to illustrate and visualise evolutionary tendencies in coccolithophores, but also in menardiform globorotaliids (Knappertsbusch and Mary, 2012). Conceptually, they are constructed by stacking the contour frequency 190 diagrams from different time levels. This way, the grid cells of plane bivariate contour frequency diagrams expand to include time as the third dimension, e.g. spiral height, axial length and time. The local frequency is the fourth dimension (F). In this manner, a four-dimensional unit (X, Y, time, F) called "voxel" is generated. The component F of a voxel (local frequency) can then be represented as iso-surface, which is done using Voxler. In other words, the iso-surface of the VDD represents the distribution of a constant local frequency through time (Knappertsbusch, 2016). High iso-values form the core of a VDD and 195 represent abundant specimens. They allow the investigation of the main evolutionary path through time. Low iso-values illustrate rare specimens and show the extremes of test size. They are often related to innovation caused by evolution or represent extreme forms introduced by dispersal.
The protocol for constructing a VDD developed by Knappertsbusch (2007; and Knappertsbusch and Mary (2012) was modified to improve the level of coincidence between the plane CFDs and VDDs. The most important changes concern (1) 200 the usage of absolute instead of normalised frequencies in the input files, (2) a different setup in the gridder option and (3) the modification of the iso-value. A detailed list of the used adjustments is given in the supplementary material (File "VDD_setups.txt").
The commercial software PDF3D ReportGen by Visual Technology Services Ltd was used to create the 3D model from the Open Inventor (.iv) file format of a VDD when exported from Voxler. 205

Results
In a first step of analysis, the test-size evolution of G. menardii at Hole 667A was investigated by plotting X and Y versus the age. This is the simplest type of analysis for evolutionary change and allows a direct comparison with previous data from Knappertsbusch (2007;. At Hole 667A, this test size variation shows different phases of evolution through time: As will be demonstrated in the following section, these two parameters serve as a primary measure for the intraspecific variability of 210 the G. menardii lineage.

Morphological Parameters through Time
The comparison of the test size of G. menardii during times of co-occurrence with its sister taxa G. limbata and G. multicamerata and the size after the extinction of G. limbata and G. multicamerata may give evidence about possible shifts in the ecology of G. menardii. Major changes in the size of G. menardii before and after the extinction of its sister taxa probably 215 point to an adaption to a different, new niche, e.g. in terms of "incumbency replacement" (Rosenzweig and McCord, 1991).
Between 7.96-2.58 εa, the evolution of X in G. menardii shows three peaks at 7.11 Ma, 5.78 Ma and 3.99 Ma in the mean and median values (Fig. 5). Except for the sample 667A-10H-1, 97-98 cm at 5.26 Ma, at which the maximum size of G. menardii does not decrease as the mean and median do, the maxima of X follow the trends of the corresponding mean and median values. The maximum values exhibit one peak at 7.11 εa and two "peak plateaus" from 5.78 εa to 5.26 εa and 4.35 X peaks (Fig. 5d) as G. menardii at 7.11 Ma, 5.78 Ma and 4.14 Ma. In average, populations of G. limbata are slightly larger in size than those of G. menardii. Specimens with seven chambers in the last whorl, which are considered as G. limbata, still 225 occur after 2.58 Ma, but only sporadically and in low numbers and no statistically significant statements are possible for those times.
Globorotalia multicamerata attains the largest size of the three species at times before 3 Ma (Fig. 5c). It surpasses G. menardii and G. limbata in test size mean and maximum values in all samples in which it occurs (Fig. 5d). Exceptions are the samples at 6.07 Ma, in which it has the same mean value as G. limbata, and at 2.057 Ma. No specimen was found at 5.78 Ma. Thus, 230 G. multicamerata only exhibit one major peak in the maxima values at 3.69 Ma and in the mean values at 3.99 Ma. Similar to X, the mean and median values of Y also show three major peaks (7.11 Ma, 5.78 Ma, 4.14 Ma) for G. menardii and G. limbata between 7.96 Ma and 2.58 Ma (Fig. 6, 7). Maxima of Y exhibit similar peaks, but note a fourth peak in Y at 3.204 Ma for G. menardii (Fig. 6a). 235 Measurements of Ar are shown in figure 7. Between 7.96 Ma and 2.58 Ma Ar of G. menardii reveals three peaks at 7.11 Ma, 5.78 Ma and 3.204 Ma and a plateau from 4.35 Ma to 3.99 Ma. The data also show a peak in Ar for G. limbata at 4.14 Ma.
For G. multicamerata, the maximum and mean Y and Ar values show a similar pattern as X (Fig. 6c, 7c), but with a major 240 peak at 4.14 Ma. This species exhibits the largest size in these two parameters in comparison to the other two species.
The three parameters show a high degree of overlap between the three species. However, morphological overlap between these species point to strong interspecific size variation. Globorotalia multicamerata exhibited the largest mean population test size, G. menardii the smallest mean size, while G. limbata was intermediate.

Contour Frequency Diagrams of Spiral Height and Axial Length 245
As already mentioned in the methods section (chapter 2.5.), CFDs may help to detect patterns of cladogenetic splitting or anagenetic evolution by identifying shifts in the dominant test size of populations through time. The underlying grid-cell size for CFDs (and VDD in the next section) is 50 µm in X direction and 100 µm in Y direction.
In general, the contour frequency plots of G. menardii (Fig. 8) show that size measurements vary almost linearly by a diagonal semicontinuous morphocline in the X and Y morphospace. This trend is due to a flattening of the test during the ontogenetic 250 growth of the individuals (Caromel et al., 2016). As was already recognised in the univariate parameter-vs-time diagrams, two different phases of shell size development can be distinguished in the CFDs. The first phase ranges from 7.96 Ma until about 2.88 Ma and is characterised by populations with a dominant test size smaller than 300 µm in X and smaller than 600 µm for Y for G. menardii, and predominantly unimodal distributions of the population size (Fig. 8)    0.003 Ma) and the sample from 4.35 Ma visually display a bimodal distribution, in which the peaks are separated either at ca. X = 200 µm or at X = 300 µm (Fig. 8). Whether or not this pattern indicates speciation within G. menardii menardii will be investigated in the following. In case of a speciation, modal centres would connect into continuous branches that diverge for the last 2 Ma. Populations need to be closer inspected, which is done in a vertical section of stacked CFDs via a so-called Volume Density Diagrams (see next chapter 3.4). 270 Bimodal pattern may also be caused by seasonality. It is known that the annual shift of the trade winds and of the Intertropical Convergence Zone (ITCZ) influence the thermocline depth on both sides of the tropical Atlantic Ocean (Merle, 1983;Chaisson and Ravelo, 1997). Different seasonal expression in the depth of the thermocline may have caused a different reaction to growth in vertically separated populations leading to different modal distributions (Chaisson, 2003).

Volume Density Diagrams 275
The iso-surface of figure 9 illustrates the test size of rare, often innovative specimens, which either evolved within the Atlantic Ocean or intruded by dispersal. As the VDD is basically a stacking of the individual CFDs, it shows the same peaks at 7.11 Ma, 5.78 Ma and 4.14 Ma for G. menardii. The VDD clearly illustrates the size decrease during the interval from 4.14 Ma until 2.58 Ma, and the striking size increase from 2.58 Ma to 2.057 Ma (Fig. 9a). The size reached at 2.057 Ma is unprecedented.
Of special note is the aberrant steeper slope of the youngest CFD (0.003 Ma; Fig. 8), which is displayed with respect to rest 280 of the VDD towards elongated and flattened specimens. Such a trend to flat specimens was also observed in the uppermost Quaternary of DSDP Site 502 (Knappertsbusch, 2007). In the present case specimens have developed a strong keel and so are presumably not classified as G. menardii cultrata.
An interactive version of the VDD can be found in the 3D-PDF file "VDD_3D_PDF.pdf" (see supplementary materials). clades is indicated around 1.735 Ma. This sample was already mentioned to develop bimodality in CFDs (Fig. 8). In the youngest part of the core this bifurcation is no longer observed, despite the presence of distinct modal centres in individual 295 CFDs, and G. menardii tends to gradually increase its test size.
The complexity of the size evolution of G. menardii through time is further illustrated in two parallel sections in 45° orientation with different offsets and three orthogonal sections at 135° (Appendix Fig. A2-A7). The different perspectives of the VDD show other density peak trends. An "ideal" description of maximal evolutionary trends would require a flexural vertical section plain at 45°. 300

Changes in Coiling Direction in G. menardii 305
The data also show changes in the coiling direction of G. menardii, which may be related to understand evolutionary changes (see for example Bolli, 1950). In the ETAO, three different periods in the predominant coiling direction of G. menardii were observed (Fig. 11a). In the first period from 7.96 Ma until 5.268 Ma, coiling seems to frequently swing between sinistral and dextral. During the second period from 5.268 Ma to 2.057 Ma dextrally coiled specimens dominated (>90 %, except at 2.58 Ma with 78.5 %). In the youngest 310 period, lasting from 2.057 Ma to present, sinistral coiling prevailed strongly (>95 %). These periods are in agreement with Bolli and Saunders (1985) and references therein (Bolli, 1950;Bermúdez and Bolli, 1969;Robinson, 1969;Bolli, 1970;Lamb and Beard, 1972;Bolli and Premoli Silva, 1973). It is interesting that sites from the WTAO (925B), the Caribbean Sea (502) and the eastern tropical Pacific Ocean (503) exhibit a similar history of changes in the coiling direction in menardiforms (Fig.   11), although phase 1 extents in these sites until ca. 4.15 Ma and the stratigraphic resolution for trans-oceanic correlation 315 remains rather low.
Nevertheless, the reversal in the preferential coiling direction from dextral (phase 2) to sinistral (phase 3) at ca. 2 Ma is nearly synchronous in all of the above-mentioned sites and coincides with the stratigraphic entry of giant G. menardii forms in the Atlantic Ocean.
Combined, these observations may point to the establishment of a new Atlantic G. menardii clade past 2 Ma, that is also seen 320 in the main clade from the VDD at Hole 667A.
Interestingly, the giant sinistral specimens ( Y = >1000 µm) occurred already in the eastern tropical Pacific Ocean Site 503 2.58 Ma ago, ca. 0.5 Myr earlier than in the Atlantic Ocean (Fig. 11b).

Size Variation in Globorotalia menardii 325
A striking test-size increase of G. menardii is observed at Hole 667A. Within the short time interval from 2.58 to 2.057 Ma, the size more than doubles (Fig. 5, 6, 7, 8). Knappertsbusch (2007; observed a similar expansion in test size evolution in western Atlantic ODP Hole 925B and the Caribbean Sea DSDP Site 502 between 2.58 Ma and 1.95 Ma and 1.7 Ma, respectively. He considered two hypothesis which could explain this observation: a rapid faunal immigration via Agulhas Leakage or rapid evolutionary test-size increase by punctuated evolution. 330 The new data from Hole 667A are discussed in context of these two hypotheses. A third hypothesis is introduced which explains the sudden test-size increase by a rapid development of the thermocline strength.

Agulhas Leakage Hypothesis
In the Agulhas Leakage hypothesis, G. menardii is assumed to have been entrained from the subtropical Indian Ocean into the tropical Atlantic Ocean by episodic and especially strong Agulhas Faunal Leakage events (Knappertsbusch, 2016). 335 The Agulhas Leakage is known to disperse Indian Ocean biota into the Atlantic Ocean on a large scale via giant eddies (Peeters et al., 2004;Caley et al., 2012;André et al., 2013;Villar et al., 2015). These eddies form when watermasses of the Agulhas Current separate from the retroflection-point off South Africa (e.g. Lutjeharms and Van Ballegooyen, 1988;Norris, 1999 Fig. 1). At ODP Site 1087, which is located in the southern Benguela region, the Agulhas Leakage was found to exist since 1.3 Ma by presence/absence of G. menardii (Caley et al., 340 2012).
Globorotalia menardii is a well-known tropical dweller (Caley et al., 2012;Schiebel and Hemleben, 2017  An alternative idea is proposed by Norris (1999), according to which unfavourable environmental conditions in the WTAO prevented G. menardii to stabilise viable populations, which could explain the size differences during 2.58-1.95 Ma at Hole 360 925B. The Indian Ocean-influenced watermasses were perhaps further transported to the ETAO via the North Equatorial Countercurrent (Fig. 1), where more favourable conditions prevailed, allowing G. menardii to thrive. A similar hypothesis of presence and absence of suitable environmental conditions was already considered to explain a distinct short pulse of Globorotalia truncatulinoides in the southern Atlantic Ocean at 2.54 Ma (Spencer-Cervato and Thierstein, 1997;Sexton and Norris, 2008). 365 According to Chaisson and Ravelo (1997), a tradewind seesaw between the ETAO and the WTAO prevailed, which possibly let to unfavourable environmental conditions for G. menardii at Site 925B between 2.5-1.95 Ma. These authors argue that tradewinds influence the thermocline depth at each side of the equatorial Atlantic Ocean in a reverse way: Increased trade winds in the WTAO pile up warm surface waters, leading to a massive thermocline layer and a deeper thermocline. At the same time in the ETAO, increased trade winds shoal the thermocline by inducing upwelling and hence cooling the sea surface 370 temperature. This is in agreement with reconstructions (Billups et al., 1999), observations (Niemitz and Billups, 2005) and models (Merle, 1983;Ravelo et al., 1990) about seasonal latitudinal shifts in the position of the trade winds and the ITCZ. Globorotalia menardii is a typical thermocline dweller (Fairbanks et al., 1982;Curry et al., 1983;Thunell and Reynolds, 1984;375 Keller, 1985;Savin et al., 1985;Ravelo et al., 1990;Schweitzer and Lohmann, 1991;Gasperi and Kennett, 1992;Ravelo and Fairbanks, 1992;Gasperi and Kennett, 1993;Hilbrecht and Thierstein, 1996;Stewart, 2003;Steph et al., 2006;Mohtadi et al., 2009;Regenberg et al., 2010;Wejnert et al., 2010;Sexton and Norris, 2011;Davis et al., 2019) and may react sensitively in reproduction, abundance and morphology to vertical shifts of the regional thermal surface-water stratification.
The observed changes in the predominant coiling direction (Fig. 10) also support the AL hypothesis. The giant and sinistrally 380 coiling G. menardii form was first observed in the eastern tropical Pacific Ocean Site 503 at 2.58 Ma, while it occurred in the Atlantic Ocean Site 667 ca. 0.5 Myr later. Since the final closure of the Panamanian Isthmus between 4 Ma and 2.8 Ma (Chaisson, 2003;Bartoli et al., 2005) prohibited a direct water exchange between the tropical Pacific Ocean and the tropical Atlantic Ocean, the coiling evidence would rather call for spreading of the giant type from the Pacific Ocean into the Atlantic Ocean via the Indian Ocean and the Agulhas Leakage route within 500 Kyr. A study within the Indian Ocean is currently under 385 progress to test the Agulhas Leakage hypothesis.

Punctuated Gradualism by Local Evolution and/or Environmental Adaptation
A regional, more punctuated evolution of G. menardii into giant forms is another possible process to explain the observed test size pattern in the tropical Atlantic at ca. 2 Ma.
In PF and other planktonic microfossils, speciation is sometimes observed to happen within short time. Examples include the 390 classic case of fast speciation of Globorotalia tumida from Globorotalia plesiotumida within only 600 Kyr during the late Miocene/early Pliocene at DSDP Site 214 in the southern Indian Ocean (Malmgren et al., 1983). An even more rapid speciation for the same group is supposed for the western tropical Pacific (ODP Hole 806C), where G. tumida evolved from its ancestor G. plesiotumida in the late Miocene and early Pliocene within 44 Kyr (Hull and Norris, 2009). Pearson and Coxall (2014) observed transitions in the Hantkenina genus from a normal-spined to a tubulospined form within only 300 Kyr. In case of the 395 Pliocene radiolarian Pterocanium prismatium cladogenetic speciation from its ancestor P. charybdeum was reported to occur within 50 Kyr (Lazarus, 1986 A persistent question remains, however: Why did such rapid evolutionary change take place especially and only at the time between 2.3 Ma to 2.057 Ma, and not earlier or later? Answers may be sought in the final closure of the Central American Seaway from ca. 4 Ma and until 2.58Ma-2.057 Ma (Chaisson, 2003). Perhaps mainly the establishment of northern hemisphere ice sheets (Raymo, 1994;Tiedemann et al., 1994;Bartoli et al., 2005)  such rapid evolutionary events. The global climate cooling caused fundamental changes in the stratification of the upper water column (Chapman, 2000) and undoubtedly led to unfavourable environmental conditions for species like menardiform globorotallids in the Atlantic Ocean (see chapter 4.1.3). An ongoing deterioration in viability under environmental pressure of NHG presumably caused first the extinction of G. multicamerata after 2.88 Ma and then the (pseudo-)extinction of G. limbata after 2.58 Ma at Site 667. Isotopic measurements (Keller, 1985;Gasperi and Kennett, 1993;Pfuhl and Shackleton, 2004) 410 suggest that also G. limbata and G. multicamerata were thermocline dwellers, with G. multicamerata living at the top, G. limbata in the centre and G. menardii at the bottom of the thermocline.
These ecological niches were occupied during relatively rapid adaption and evolution from the ancestral G. menardii sensu the "incumbency replacement" process of (Rosenzweig and McCord, 1991). Support for such a process is also given by consideration of the maximal test growth values attained by the involved species. After extinction of G. multicamerata and 415 G. limbata in the course of the NHG, their niches in the upper to middle thermocline became liberated and could be re-occupied by G. menardii. The settlement of latter species at higher levels in the watercolumn may have led to optimum growth and development of larger tests.
Unfortunately, the temporal sampling resolution of this study is too coarse to prove the hypothesis of a punctuated or gradual evolutionary event, but could be resolved as soon as higher temporal and spatial sampling intervals are investigated at Hole 420 667A, 925B and Site 502 in the period between 2.3 and 2.057 Ma.

Possible Influence of the AMOC Strength
Unexpectedly, the measured variations of test-size maxima of G. menardii show in phases a rough parallel trend with the dissolved radiogenic isotope composition of Neodymium ( Nd), which is an indicator for the relative long-term strength of the AMOC (Dausmann et al., 2017;see Fig. 10, 12). During time intervals of increasing test size, Nd, and thus the strength of 425 the AMOC, appeared to generally increase as well. In contrast, a decrease of the test size is on average accompanied by a decreasing trend in Nd, the latter suggesting a weak AεOC. The AεOC is the Atlantic part of the global ocean conveyor belt, which causes a redistribution of heat within the global oceans. At the surface, warm and salty water is transported from the South Atlantic Ocean via the Caribbean Sea into the North Atlantic. There, it sinks down, caused by a loss of buoyancy due to the release of heat, and flows southward at depth as the North Atlantic Deep Water. The release of heat in the North 430 Atlantic influences the climate of northeastern Europe, leading to relatively mild winter temperatures. (McCarthy et al., 2017) Nd is used as a tracer for ocean circulation (Dausmann et al., 2017;Blaser et al., 2019). Erosion and weathering of continental crust, which displays characteristic, isotopic signatures from the Samarium-Neodymium decay system for different continents, is the source of dissolved Nd in the ocean water. After entry to the sea, convection of the characteristic Nd signature to deep waters allows this tracer to reconstruct large-scale ocean circulations (Blaser et al., 2019). 435 Water originating in the Northern Atlantic is known to develop more negative Nd values in comparison to waters of other origin (Dausmann et al., 2017;Blaser et al., 2019). In the study of Dausmann et al. (2017) a continuous high-resolution record for Nd at ODP Site 1088 in the Southern Atlantic was (Fig. 1)  interval and larger-regional settings of the present study.
Although very preliminary, the present empirical observation of a possible relationship between G. menardii size trends and Nd suggests that a connection between menardiform test size and AMOC strength exists. However, for Hole 667A the correlation between the maximum shell size and the linear interpolation of Nd values by Dausmann et al. (2017) remains poor (R² = 0.1477, Appendix Fig. A8). The poor correlation arises by several "outliers" in the time interval from 2.057 Ma to 0.73 445 Ma. In this interval, six out of seven samples show an out-of-phase trend of the maximum test size with Nd. No explanation can be delivered for this observation at the moment.   (Haarsma et al., 2008;dos Santos et al., 2010) of the upper ocean water column in the tropical Atlantic Ocean. 450 Thus, changes in the strength of the AMOC come into mind, which shifted the position of the ITCZ and associated trade winds (Billups et al., 1999;Timmermann et al., 2007), and which in turn affect the thermocline strength (Merle, 1983;Chaisson and Ravelo, 1997;Wolff et al., 1999). It is for example known that the ETAO thermocline reacts sensitively to variations in the AMOC strength (Haarsma et al., 2008;dos Santos et al., 2010). In this manner, the habitat of G. menardii would have been altered as well. In this context, a model for the response of test size of G. menardii under a changing thermocline is presented 455 in the next section.

A Thermocline Model for Size Variation in G. menardii
A number of stable isotopic studies (Curry et al., 1983;Keller, 1985;Savin et al., 1985;Schweitzer and Lohmann, 1991;Gasperi and Kennett, 1992;Ravelo and Fairbanks, 1992;Gasperi and Kennett, 1993;Stewart, 2003;Steph et al., 2006;Mohtadi et al., 2009;Regenberg et al., 2010;Wejnert et al., 2010;Davis et al., 2019), plankton tows (Fairbanks et al., 1982;Thunell 460 and Reynolds, 1984;Ravelo et al., 1990), census data from sediments (Sexton and Norris, 2011), and in situ observation (Hilbrecht and Thierstein, 1996) showed that G. menardii preferably dwells in the thermocline. According to Sexton and Norris (2011) and references therein, this coincides often with vertical habitats of increasing organic particle concentration and segregation, a zone in the thermocline where oxygen consumption due to particle degradation is high and where oxygen content becomes lowered. 465 Changes in the test size of PF are thought to be related to changes in the environmental conditions (Hecht, 1976;Malmgren and Kennett, 1976;Naidu and Malmgren, 1995;Schmidt et al., 2004;André et al., 2018): Under optimum conditions, test size of species increases to its maximum, while under non-optimum conditions, the size is reduced, although detailed physiological processes at individual levels are still not entirely understood. thermocline leads to a stronger density gradient between the surface and the subsurface layer. Often the chlorophyll maximum zone is located at this boundary (Fairbanks et al., 1982;Ravelo and Fairbanks, 1992;Steph et al., 2006), where marine snow accumulates (Möller et al., 2012;Prairie et al., 2015). The increased concentration of degrading particulate organic matter hances nutritional conditions and favours the test growth of G. menardii. It is, for example, known that nutrient-rich conditions facilitate test-size increase in the PF species Globigerinoides sacculifer (Bé et al., 1981), Globigerinoides ruber, Globigerinita 475 glutinata, Globigerina bulloides and Neogloboquadrina dutertrei (Naidu and Malmgren, 1995).
The thermocline may play a crucial role in other aspects of G. menardii's lifecycle as well. A strong thermocline and the corresponding high density contrast is thought to concentrate its gametes and food particles at a narrower zone and thus increases their chance to survive (Norris, 1999;Broecker and Pena, 2014). This model of ecological factors within the regional thermocline influencing the phenotypic expression of G. menardii fits with Sexton and Norris' (2011) deglaciation proliferation model modulating the stratigraphic distribution of G. menardii. They suggest that G. menardii tracks thermoclines in areas with a moderately low oxygen concentration of ~50-100 µmol kg -1 , probably reduced by the degradation of organic matter.
Furthermore, Sexton and Norris (2011) postulate the reduction or vanishing of G. menardii populations during glacial times 485 due to better ventilated surface water masses, i.e. a weaker thermocline. Weakening of the AMOC during glacial times (Broecker, 1991;Buizert and Schmittner, 2015) and associated changes in the position of the ITCZ led to a weakening and/or re-positioning of the thermocline, so that ambient conditions became less suitable for growth and proliferation of G. menardii.
The proposed thermocline hypothesis (Fig. 13) offers a possible way to explain the striking test size increase in the Atlantic Ocean after severe cold chills and between 2.58 and 2.057 Ma. Perhaps it can also explain several periods of test-size increase 490 and decrease in the G. menardii lineage in the Pacific and the Atlantic Ocean within the last 8 Ma, assuming analogous conditions.
A causal chain of physiological processes in order to explain the empirical similarity between the AMOC strength and the evolution of test size of G. menardii, however, still remains elusive and need further investigation.

Conclusions 495
1. Test-size measurements of the planktonic foraminifer Globorotalia menardii from the ETAO ODP Site 667A show a striking size increase in the early Pleistocene and a test-size evolution during the past 8 Ma, which was similar to observations observed in the tropical Atlantic and Caribbean Sea (Knappertsbusch, 2007; 5. At present, the alternative hypothesis of a regional and punctuated evolutionary event cannot be dismissed until more paleobiogeographic data are available at higher geographic resolution, especially from the Indian Ocean realm.
6. The results of this paper show that for an improvement of a taxonomic distinction between closely related species with high 520 morphological overlap, it is strongly necessary to better include temporal measurements of morphological divergence.

Code 535
The modified MorphCol programs, which were used to process the raw data, as well as their codes will be available at PANGAEA (www.pangaea.de).

Data availability
The full set of derived and raw data and images will be deposited at PANGAEA repository (www.pangaea.de).
The supplied zip archive Supplementary_Material.zip is an extract of all data and contains the necessary data to reproduce the 540 illustrated figures.

Sample availability
The sample material is deposited in the collections of the Natural History Museum Basel, Switzerland.

Competing interests.
The author declare that he has no conflict of interest.