We combined the ForCenS planktonic foraminifera assemblages sampled in Atlantic Ocean core-top surface sediment samples (Siccha and Kučera, 2017) with 54 additional samples for the NE Atlantic coastal upwelling systems (Salgueiro et al., 2008; Voelker and Salgueiro, 2017; Salgueiro et al., 2020, represented in the figure as Salgueiro) and 40 samples from the North Atlantic’s subpolar gyre (Barrenechea Angeles et al., 2020; Sahoo et al., 2022). This resulted in a total of 1741 samples ( Figure 1 ), which were collected between the 1960’s and 2016. (For metadata standardization see Supplementary Material).

Biological diversity comprises two main aspects: species richness and evenness or equitability (Magurran, 2004). In ecological studies, both the total number of species (species richness, S), the Shannon-Wiener diversity index (H) and evenness index (E) are commonly used to estimate community diversity. Species richness can be defined as the total number of species encountered in a sample (S). Although species richness may depend on sample size, in this case, the use of census counts (in % of total counts) refutes that possibility as it departs from an equal sample size. Evenness describes the distribution of individuals among the species, i.e. it is “Even” if all species (yet if few) are equally abundant (“even”/more diverse), or “uneven” if the community is dominated by one or a few species only (less diverse). To measure community evenness, the Probability of Interspecific Encounter (E) (Hurlbert, 1971), also known as Simpson’s evenness index and Gini-Simpson index, was used. The E index is unbiased relatively to sample size, and it represents the probability that two randomly chosen individuals from an assemblage will represent two different species. According to Gotelli and Ellison (2004), this index produces easily interpretable units of probability and corresponds intuitively to a diversity measure that is based on the encounter of a novel species, and it varies between 0 and 1 (with 1 representing a community with all species equally abundant (perfectly even, and therefore no species is dominant over another) and 0 representing a highly uneven one, i.e. dominated by one single species). The alpha diversity indices were calculated using the Vegan package (Oksanen et al., 2020) in the R free software (http://www.r-project.org/) and predicted over the sampled area into a regular grid by applying the spatial model estimated variogram, using kriging interpolation, for each index separately (Bivand et al., 2008; Cressie and Wikle, 2011). Kriging has been proven to be amongst the best interpolation methods to represent marine species distributions, when only geographical variables are considered (Rufino et al., 2021).

In the current work, we compare foraminifera assemblages and the calculated diversity with SST, Chl-a concentration (as a surrogate for photosynthetic biomass) and EKE, all variables that can be approximated by satellite synoptic measurements at a large scale ( Figure 2 ).

Annual averaged SST was extracted from the World Ocean Atlas 2013 (WOA13: Garcia et al. (2014)), which is a set of objectively analyzed (0.25° grid) climatological fields. The used variables correspond to the average of five “decadal” climatologies: 1955-1964, 1965-1974, 1975-1984, 1985-1994, and 1995-2006. Although salinity and nutrients in WOA13 (such as phosphate, nitrate, silicate) were initially considered in a preliminary analysis, we ended up not using them due to the strong collinearity with the selected variables (see Supplementary Figure 1 ).

Annual averaged Chl-a concentration (mg/m^-3) is the product generated by the NASA Ocean Biogeochemical Model (NOBM) based on data assimilation of remotely sensed Chl-a. Garver–Siegel–Maritorena model bio-optical products were derived from SeaWiFS and MODIS-Aqua satellite data and processed by ICESS at the University of California – Santa Barbara (https://oceancolor.gsfc.nasa.gov).

The EKE per unit mass are the square of the velocities (cm²/s²) from the SSALTO/DUACS (Segment Sol multi-missions dALTimetrie, d’Orbitographie et de localisation précise/Data Unification and Altimeter Combination System) delayed-time level-4 product distributed by Copernicus Marine Environment Monitoring Service (https://resources.marine.copernicus.eu/product-detail/SEALEVEL_GLO_PHY_CLIMATE_L4_MY_008_057/INFORMATION). The sea level anomaly is the sea surface height above mean sea surface that is included in the absolute dynamic topography, and the mean was calculated over the 1993 – 2012 period. The velocities are geostrophic, estimated from the sea surface height anomaly (for the EKE) based on multi-satellite altimetry and thus should be interpreted as vertically averaged velocities dominated by the upper layer. Pseudo-streamlines are estimated from the geostrophic velocity field based on the absolute dynamic topography to provide a visualization of the mean currents ( Supplementary Figure 2 ). EKE (log₁₀ transformed) was therefore used as a proxy for ocean instability and mixing (Qiu et al., 2014; Uchida et al., 2017) and the streamlines mark the presence of currents and fronts associated to them. In fact, EKE marks the subtropical gyres western boundary currents: the Gulf Stream with onset of the North Atlantic Drift (and may be some interaction with the Labrador Current off Newfoundland); the Malvinas Current/Brazil-Malvinas Confluence Zone and parts of the Antarctic Circumpolar Current; the onset of the South Atlantic Current; the North Equatorial Counter Current and less strong the South Equatorial Current; and the Agulhas “leakage” around South Africa.

The relationship between the alpha diversity indices (H, S and E) and the oceanographic variables SST, Chl-a and EKE used here as explanatory variables or covariates, was statistically analyzed using a generalized additive mixed model (GAMM) approach. EKE and Chl-a were log₁₀ transformed to stabilize the variances, i.e. lx = log₁₀ (x+1) (a base of 10 was used to facilitate the interpretation). For H diversity, ‘Gaussian’ distribution was used, whereas for S (counts of species), the ‘Tweedie’ distribution was considered, since according to Wood et al. (2016), it is suitable for count data (using an estimated power of 1.01) and followed by the checking of the residual’s analyses. For E, a Beta distribution was used, with a logit function family, as recommended for proportions between [0,1]. The assumption of the lack of collinearity was tested (correlation between covariates) by estimating Spearman’s correlation coefficient between the three explanatory variables considered, and it was found that these were always below 0.7, which is the recommended rule of thumb (Zuur et al., 2010).

Spatial autocorrelation implies that samples that are closer to each other are more related than samples that are farther apart and therefore, are not independent, as assumed by most statistical models. To account for spatial autocorrelation, a statistical smoothing term was incorporated in the model as a random effect. In the H and S indices analysis, three spatial auto-correlation structures were evaluated for the random part of the model (exponential, spherical and, Gaussian), and, for each case the most appropriate one was selected using Akaike’s information criterion residuals distribution analysis and ‘acf’ and performance measurement criteria. In the case of Evenness, a term was added to the model, encompassing both latitude and longitude to account for the spatial component, since the use of Beta family did not permit the inclusion of a spatial autocorrelation structure (Wood, 2006).

In all cases, a full model was initially estimated using the three covariates, which was then subject to a model selection procedure, that included the interactions between terms. The best combination determined for each index based on Akaike’s information criterion (Zuur et al., 2010). Model assumptions were evaluated by analyzing the residuals distribution, including residuals auto-correlation (‘acf’) and covariates concurvity (to test collinearity) (Wood, 2006).

Gradient forest analysis (GF) is a method that has been proven to be highly robust and efficient in identifying broad-scale environmental gradients of biodiversity patterns (Leaper et al., 2011; Ellis et al., 2012; Pitcher et al., 2012). This method has been used to: determine thresholds in/as ecological indicators (Large et al., 2015; Waite et al., 2020); identify thresholds of anthropogenic change in marine communities (Couce et al., 2020) in genomics (Fitzpatrick and Keller, 2015); determine the major drivers of phytoplankton distribution in lakes (Roubeix et al., 2016). To consider multiple, co-occurring species in model fitting, a community-level model was used (Nieto-Lugilde et al., 2018). Within these the Gradient Forest Analysis (GF) was chosen because tree-based techniques are among the top performing technique for species distribution analysis that are capable of modeling nonlinear responses and interaction effects, and permit to identify species turnover thresholds along environmental gradients. Using a large set of independent data for evaluation, this technique was found to be highly effective at summarizing spatial variation in both the composition of species assemblages and species turnover (Stephenson et al., 2018). Gradient Forest analysis was therefore applied to the relative abundance of foraminifera and the three oceanographic satellite variables to estimate the location and importance of the community composition thresholds along environmental gradients in the Atlantic Ocean (Ellis et al., 2012). For this analysis, the species relative abundance was natural log transformed, after adding a constant of one (log1p R function). The GF (Pitcher et al., 2011; Ellis et al., 2012) is a machine-learning modeling technique that predicts the distribution of a taxon as a function of relevant environmental variables, on the basis of fitting an ensemble of regression or classification trees. That is, multiple random forest models are run for each taxon separately and then simultaneously trained together to represent the community. Standardized measurements of change along environmental gradients for all taxa are therefore performed simultaneously to analyze the combined results.

This information is used to build ‘response curves’: empirical functions of change in composition for each environmental variable. After accounting for the distribution of the data along the gradient of the predictor, thresholds for which significant changes of community composition take place can be identified for each of the predictive variables, using the frequency histograms of the values at which splits happen in the individual decision trees. Furthermore, GF allows for the quantification of the relative importance of the predictor variables by averaging them across all taxa. Community threshold is defined as the zone, the critical value(s), along an environmental gradient where sharp increases or decreases in the occurrence of several taxa boost a change in community composition.

Results from this analysis were then used to transform the environmental layers to represent biological space, using the oceanographic layers: SST, EKE and Chl-a. A principal component analyses (PCA), was also applied to the predictions, and the first three axis of the PCA scaled into RGB space (Red-Green-Blue) for mapping species turnover. In addition, the main biogeographic regions of planktonic foraminifera were estimated using the k-means cluster method. The number of cluster groups was selected using the Elbow method ( Supplementary Figure 3 ). Further details on the method and applications can be found elsewhere (Ellis et al., 2012; Baker and Hollowed, 2014; Large et al., 2015; Nieto-Lugilde et al., 2018; Stephenson et al., 2018; Couce et al., 2020). All the analysis were done using the R libraries mgcv (e.g., Wood, 2011; Wood et al., 2016; Wood, 2017), the package gradientForest (Ellis et al., 2012), and the package NbClust (Charrad et al., 2014) and factoextra (Kassambara and Mundt, 2020), for the determination of the number of clusters. Visualization of the GAMMs was done with the aid of library ‘mgcViz’ (Fasiolo et al., 2018), in R software.

Spatial auto-correlation on alpha diversity indices was studied for the three indices and the variograms fitted with a spherical model (model results: S: nugget = 4.367, sill = 68.146 and range = 24669.715 (which would be equivalent to a linear model); H: nugget = 0.036, sill = 0.145 and range = 1930.019; E: nugget = 0.002, sill = 0.021 and range = 1378.223; Supplementary Figure 4 ). Species richness (S – Figure 2 ) shows the lowest number of species in the Arctic Ocean (1-4 species), followed by the Nordic Seas (1-11 species). Between 50 °N and 50 °S, the Atlantic shows almost symmetrical bimodal diversity richness south and north of the equator, but peaking values appear between 25 °N – 25 °S (15 to 27 species). Shannon diversity (H – Figure 2 ) is lower than 1 (corresponding to < 10 species) only close to the poles. Evenness (E – Figure 2 ) shows that 90% of the samples have less than 0.5 of encounter probability, that is, the probability of finding the same foraminifera species twice is less than 50% in most of the area. Ninety-eight percent of the samples with an E < 50% contain < 10 species of planktonic foraminifera and were located closer to the poles (> 50 °N and > 50 °S). In summary, the centers of the subtropical gyres show lower diversity (H), a lower number of species (S), and lower E (i.e., more uneven species distribution). In contrast, the equatorial region and the northern and southern boundaries of the hemispheric subtropical gyres and the eastern boundary currents exhibit higher species richness, a more even distribution of species abundance (lower E), and a H distribution pattern closer to E than to S.

The response of the alpha diversity indices to a heterogeneous environment represented by SST, Chl-a, and EKE showed that the best model for all three indices included was given by the following formula: index ~ s(l10EKE, sst.woa13) + s(Chla). Further, adding spatial structure to the analysis improved the performance of the model and drastically decreased residual spatial autocorrelation in S and H. For S, a spherical model was fitted (adj- R² = 82%; Tweedie (1.25)), whereas for the H index, the best model also included a spatial exponential correlation structure (adj- R² = 71%). Results of the GAMMs for H, S and E show similar patterns and can be found in Supplementary Figure 4 . The effect of SST was the strongest one for all indices (S, H, E) and controlled the large diversity patterns. A significant interaction between SST and EKE indicates that the relationship between EKE and diversity changes across the SST ranges (depicted by the interaction panel color shade, from yellow to blue, while the black lines represent species richness variability through EKE and SST, as can be observed in the Supplementary Figure 5 ). All three alpha diversity indices show an almost linear slightly negative relationship with Chl-a, but less critical than the SST/EKE interaction.

In the current work, GF analysis successfully modeled the whole assemblage simultaneously and determined the thresholds associated with continuous environmental gradients, the main species contributing to those changes, and the inference of the principal regions of occurrence along the Atlantic. The model explains a large proportion of the foraminifera assemblage variability (R² = 59%, obtained by adding the R²-values for all predictors). The most critical environmental variables for the model were SST (R² = 29%), Chl-a (R² = 19%) and EKE (R² = 11%). The primary community thresholds for foraminifera species turnover appeared at 5-10 °C and 22-28 °C SST, 0.05-0.15 mg m^-3 Chl-a and, 1.2-2.0 cm² s^-2 log₁₀ EKE energy ( Figure 3 ), respectively.

Out of the 37 foraminifera species considered in our combined dataset, only 28 ( Supplementary Material ) contributed to the species turnover model ( Figure 4 ). The species not considered by the model corresponded to those whose maximum abundance was ≤ 2% in our data set. The response of each species to the different environmental variables presented in Figure 5 , reveals that the most important species highlighted in each environmental gradient have R² > 0.7 ( Figure 4 ). From these, six species are also fundamental for transfer functions (Jonkers and Kučera, 2019), which further reinforces the results obtained and the analysis used. As expected, due to its dominance in polar waters, Neogloboquadrina pachyderma (purple in Figure 5 ; most relevant species in GF) is within the species that show marked thresholds relative to SST, with a sharp, pronounced transition between 7 and 9 °C. In contrast, the tropical species Trilobatus sacculifer (red; second most relevant species in GF) and Globigerinoides ruber (formerly Globigerinoides ruber pink; orange) mark the extreme warm threshold between 22 and 27°C. Neogloboquadrina incompta (blue), the third most important species (4^th in GF), signals a more abrupt threshold at 7°C followed by a stepwise transition between 15 and 23°C. Globoconella inflata (green) shows thresholds at 10°C and 22°C.

The most sensitive species to changes in Chl-a ( Figure 5 ) are Globigerinoides white (formerly Globigerinoides ruber white) (green), Globigerinoides conglobatus (red), Globorotalia truncatulinoides (green-blue), Globigerinoides tenellus (olive green), and N. incompta (blue), all of these species with a R² ≥ 0.6 for the GF model ( Figure 4 ). Although all of them show thresholds below 0.1 mg m^-3, subtropical species G. ruber, G. conglobatus and G. tenellus present the lowest (0.03 mg m^-3) and more abrupt ones, followed by G. truncatulinoides with a more gradual transition between 0.04 and 0.11 mg m^-3, whereas N. incompta shows a break at 0.09 followed by another, more gradual one at 0.4 mg m^-3.

With respect to EKE, Pulleniatina obliquiloculata (dark pink), Globorotalia hirsuta (blue-green), Globigerinella siphonifera (green-brown), G. tenellus (olive green), and Orbulina universa (purple) show S-type relationships with steep thresholds but narrow intermediate plateaus ( Figure 5 ). Their importance to the GF model (R²) exceeds 0.58, with the exception of O. universa with a value of 0.45 ( Figure 4 ). G. tenellus, G. hirsuta and G. siphonifera behave similarly up to 0.9 cm² s^-2 followed by a clear divergence. G. tenellus shows steep slopes at 1 cm² s^-2 and 1.5 cm² s^-2, while G. hirsuta has an abrupt transition at 1.2 cm² s^-2 and then gradually increases until reaching a plateau at 2.4 cm² s^-2. G. siphonifera, on the other hand, continues the low slope increase up to 1.5 cm² s^-2, and then exhibits steep slopes at 1.7 cm² s^-2 and 2.1 cm² s^-2. P. obliquiloculata responds to EKE only above 1 cm² s^-2, and reveals a threshold at 1.4 cm² s^-2, followed by a less steep rise until 1.8 cm² s^-2 ( Figure 5 ). O. universa reveals a steep increase until 0.9 cm² s^-2, followed by a slowly rising “plateau” between 0.9 and 1.6 cm² s^-2 and then another steep increase until 2.2 cm² s^-2.

When the species turnover map is transformed into the environmental space of the oceanographic variables and the bioregions are identified by cluster analysis ( Figures 6A, B ), not only the key latitudinal regions are visible, but the main current features are also superimposed; the eastern boundary upwelling systems (EBUS) off Namibia (Benguela Current) and Northwest Africa (Canary Current) as well as the Equatorial divergence, and the northern and southern subtropical fronts, i.e., the Gulf Stream-North Atlantic Drift path in the north, and the South Atlantic current in the south.

The importance of SST as a driver for foraminifera’ diversity has been well shown in previous studies (e.g., Morey et al., 2005; Fenton et al., 2016; Rillo et al., 2022). Some authors, such as Fenton et al. (2016), attempted to evaluate the effect on the diversity of several environmental variables generally highly correlated with SST, namely nutrients and salinity. These studies did not provide additional information on the influence of ocean circulation nor studied SST thresholds for species turnover along environmental gradients. The present work reveals critical community thresholds at 5-10 °C and 22-28 °C SST, in agreement with previous observations, namely in transfer functions, and mainly related to Neogloboquadrina pachyderma and Neogloboquadrina incompta. Other than the preponderance of SST, the effect of a broader range of environmental factors have been hypothesized to explain local core top planktonic foraminifera diversity (e.g., Morey et al., 2005; Fenton et al., 2016; Rillo et al., 2022), but none of these authors considered proxies for ocean circulation. Only Richter et al. (2022) consider that plankton DNA biogeography at the global scale is linked to transport and large-scale ocean currents. Following this idea, we used EKE as a proxy of ocean fronts and circulation. SST, EKE, and Chl-a, which are all intrinsically connected with species physiology and ecology, and are likely to be the main drivers of foraminifera diversity. For example (1) SST, an index of solar energy input, influences the kinetics of biochemical reactions, metabolic rates and, thus, the probability of species mutation (Rhode, 1992; Allen et al., 2006; Eme et al., 2022), (2) Chl-a, as an index of productivity and food supply, determines population abundances, which are also directly related to the probability of mutation and species adaptation (Roy et al., 1998; Morey et al., 2005), and (3) kinetic energy, represented by EKE is a proxy for water column instability and ocean fronts, therefore it influences the rate of immigration and the probability of local extinction. These three variables reflect the main changes of energy between the individuals and the environment.

The maps of species turnover transformed to the environmental space of these oceanographic variables, and the resulting cluster groups ( Figure 6 ), clearly show the major latitudinal regions known since the work of Bé (1959) and Bé and Tolderlund (1971). Since then, the defined bioregions and the location of the respective transition zones, i.e. boundaries with great faunistic contrast between zones, have been adapted and cited with very little changes, by many authors (e.g. Schmidt et al., 2004; Kučera, 2007; Rillo et al., 2022). In our study, only four provinces were identified, unlike previous works that describe five zones, probably because a large part of the fifth province, the polar zone (like in Schmidt et al. (2004) and Bé and Tolderlund, 1971), is in our study combined with the subpolar zone, potentially due to the exclusion of samples in the Antarctic polar region. Additionally, the upwelling areas in our study are not within the contiguous provinces (tropical and subtropical), but rather correspond to zones which are farther apart from the equator (i.e. temperate zone in Figure 6 ), unlike in Schmidt et al. (2004) where it is considered as a sixth province. Otherwise, the zones and location of transition areas found in our work purely using statistics are in good agreement with the planktonic foraminifera provinces already established by Bé and Tolderlund (1971) for the modern ocean, despite slight differences in the position of the borders of each province.

The central area around the equator corresponds to the tropical zone ( Figure 6B ). The north west corner coincides with previous works, being located near 38°N (as in Bé and Tolderlund, 1971; Schmidt et al., 2004 and refs. therein), but then diverges from those studies by the diagonal border down to 10°N off NW Africa, i.e., in our study the tropical zone in the eastern North Atlantic is greatly diminished. A similar pattern is observed for the boundary in the South Atlantic that also runs diagonally from about 15°S in the East to about 31°S in the West. Not only the diagonal shape, but also the latitude differs from previous works, with around 15-18°S in Schmidt et al. (2004) and Bé and Tolderlund (1971), whereas Rillo et al. (2022) place it up to 37°S. In fact, the south-west extension of the tropical zone in the South Atlantic in the current work largely overlaps with the subtropical zone ( Figure 6B ) in previous works. To sum up, in our study the subtropical zone gains area in the eastern North Atlantic and loses area in the western South Atlantic. The temperate zone (or transitional zone as termed in Bé and Tolderlund, 1971; Figure 6B ), also marks the upwelling regions in our study, especially off SW Africa and in the Cape Blanc upwelling system (NW Africa). The south-west transition of the temperate zone is always located at the tip of South America (54°S), in agreement with the previous studies. The north-west temperate zone limit reaches up to 52°N in the southern Labrador Sea, whereas in previous definitions it was located near ~40°N, representing the beginning of the sub-polar zone ( Figure 6B ). The expansion of the temperate zone observed in this study in the North Atlantic can clearly be related to the influence of the North Atlantic Current and its branches (e.g., Daniault et al., 2016). Thus, the northern transition between sub-polar and temperate zone is further north in the current work, i.e., at 58°N, instead of 50°N in the previous works, whereas the transition between the temperate and subtropical zones is located further to the south (37°N vs. 42°N). The differences can be justified by the new samples included in the current study and by the statistical approach considered. Further work (and data, if possible) is required to improve the definition of the transition zones, potentially also integrating data from other taxonomic groups. Overall, there is a good correspondence between our and Rillo et al. (2022) borders of the polar and temperate bioregions, but not in the tropical and equatorial region, where our analysis clearly marks a subdivision of the two, aside from highlighting the significant circulation systems. Additionally, the provinces based on our species turnover also represent large circulation patterns and EBUS.

We believe that the location of the transitions of our bioregions improve the previous models due to higher accuracy in some aspects (e.g., reflecting main circulation pathways in the subtropical gyres). For example, in the Iberian margin where we added new samples by Salgueiro et al. (2008), the location of the transition between subtropical and temperate regions is found to be at 37°N (Cape São Vicente) and not at ~44°N as in previous works, which is in agreement with the local study of Salgueiro et al. (2008).

Bé and Tolderlund (1971) already concluded that the planktonic foraminifera characterize the major current systems of the world oceans, but quantitative analysis was not possible by then. For assessing the impacts of large-scale frontal systems on plankton foraminifera diversity we used EKE as it is the best currently available proxy at a large scale. However, further work is required in developing better proxies for this important variable to better verify their influence on marine diversity.

Differences in diversity are also depicted for EBUS. These systems, linked to the eastern boundary currents, lead to the occurrence of jets, filaments, and eddies that form along the Ekman divergence. However, their mesoscale variability is not depicted by altimetry, which is indeed known to be used to determine large-scale circulation (Cheney and Steel, 2001; Figure 2 ). Nevertheless, the Atlantic EBUS and the Equatorial Current system are clearly identified as biogeographic regions marked in the species turnover map ( Figure 6 ), revealing the proliferation of opportunistic phytoplankton in response to the upwelling of nutrient-rich waters (Margalef, 1978; Raymont, 1980). In this scenario, ecological succession is followed by opportunist omnivorous and herbivorous planktonic foraminifera species, as reported in sediment trap studies (Abrantes et al., 2002; Schiebel et al., 2004) and reflected in species turnover.

When assemblage variability is accessed by GF, diverse morphologically defined species emerge with different importance and thresholds along the gradients of the environmental parameters identified as the main diversity drivers by our GAMM. For SST the most important species identified by GF ( Figure 5 - N. pachyderma, T. sacculifer, N. incompta, G. inflata and G. ruber) are the same and appear in the same order of relevance as the ones driving transfer function results in the North and South Atlantic (Jonkers and Kučera, 2019). Furthermore, the thresholds revealed by GF for those species agree with the SST ranges used by Imbrie and Kipp (1971) to define the main planktonic foraminifera provinces: Tropical: ≥24 °C; Subtropical: 18 - 24 °C; Temperate/Transitional: 13 - 18 °C; Subpolar: 7 - 13 °C; and Polar: ≤ 7 °C. The fact that GF correctly identified the critical species and their temperature ranges, provides confidence to the results obtained by the same GF for Chl-a and EKE. Furthermore, it encourages us to validate if the compositional change marked by the individual species along each variable gradient can be explained by our current knowledge about the ecology of those species or can otherwise contribute to better understanding of their ecology. First, one can observe the species distribution pattern ( Figure 7 ). N. pachyderma shows a polar to subpolar distribution justifying its low SST threshold. Similarly, T. sacculifer and G. ruber show a direct relation to tropical waters, also previously documented (e.g., Schiebel et al., 2002; Kučera, 2007). However, according to Schiebel and Hemleben (2017) and references therein, they spread over a wide range of SST and salinity, as suggested by the GF outcome ( Figures 5 , 7 ). N. incompta and G. inflata appear allied to intermediate SSTs ( Figures 5 , 7 ). The broad SST threshold from 15 to 22°C of N. incompta, a species with a single genotype in the Atlantic Ocean (Darling and Wade, 2008), is likely to reflect the diverse areas in which it thrives, from the mixing zone between the Gulf Stream and Labrador Current in the NW Atlantic, to the Canary and Benguela EBUS, where G. inflata also contributes considerably (Meggers et al., 2002; Lončarić et al., 2007; Salgueiro et al., 2008; Lessa et al., 2020). GF spotted thresholds for G. inflata fit with the broad range of subtropical to transitional waters inhabited by the species genotype I (Morard et al., 2011). However, the ramp between 5 and 10°C may reflect the association of its genotype II to the Antarctic subpolar waters (> 40°S) (Morard et al., 2011), in particular to the Malvinas Current in the SE Atlantic ( Figure 7 ).

Considering Chl-a, the species identified by GF as the ones with a more substantial effect on community change along the Chl-a gradient are G. white, G. tenellus, G. conglobatus, G. truncatulinoides, and N. incompta ( Figure 5 ). The three species with higher importance, G. white, G. conglobatus and G. tenellus, live in oligotrophic, subtropical to tropical waters (Rebotim et al., 2017; Schiebel and Hemleben, 2017; Jentzen et al., 2018; Lessa et al., 2020). Their distribution ( Figure 7 ) displays an important presence along the Canary Current EBUS. In particular G. white, as reported by several authors (e.g., Meggers et al., 2002; Salgueiro et al., 2008), and that can explain the species higher Chl-a threshold relative to the other two species ( Figure 5 ). The presence of G. white in both oligo- and eutrophic environments reflects the existence of two different genotypes identified in the Arabian Sea associated with eutrophic and oligotrophic conditions (Seears et al., 2012). The genotype preferring oligotrophic conditions also occurs in the subtropical North Atlantic and the Mediterranean Sea.

Based on the work of Morard et al. (2019) G. white, G. conglobatus and G. tenellus are genetically linked and that might be one of the reasons explaining their similar response to environmental variability in the GF analysis. As spinose species, they prefer to prey on zooplankton like copepods and ciliates, keeping them between their spines during feeding (Schiebel and Hemleben, 2017, and references therein). Contrarily, in laboratory cultures, neanic and juvenile specimens preferentially digest phytoplankton (Hemleben et al., 1989). In addition, all three species have dinoflagellates as photosymbionts (Hemleben et al., 1989; Takagi et al., 2019) and show a positive relationship between their test size and symbionts abundance (Takagi et al., 2019). Carrying photosynthetic symbionts enables them to develop better under higher nutrient conditions. Nevertheless, the fact that juveniles feed on phytoplankton may also suggest a relation of the most important spinose species with Chl-a and explain the lack of correlation between primary production and shell flux to the seafloor (Jonkers and Kučera, 2015).

The fourth species is the non-spinose species G. truncatulinoides, which, like other non-spinose species, is assumed to be omnivorous with a preference for herbivorous food (e.g., coccolithophores, diatoms) (Schiebel and Hemleben, 2017 and references therein). Although a deep-dwelling species, living below the thermocline in waters down to 600 m (e.g., Lončarić et al., 2006; Wilke et al., 2009; Ujiié et al., 2010; Jentzen et al., 2018), even reaching 2000 m water depths at the Azores Front (Schiebel et al., 2002), it migrated vertically up to the mixed layer, when it reproduces during winter (Hemleben et al., 1989; Wilke et al., 2009; Rebotim et al., 2017; Lessa et al., 2020). Furthermore, the G. truncatulinoides morphospecies incorporate four genotypes in the Atlantic Ocean with ecological niches extending from the polar-subantarctic waters in the South Atlantic to the North Atlantic subtropical gyre (de Vargas et al., 2001; Ujiié et al., 2010). This wide ecological range is picked up by the GF analysis ( Figure 6B ) and visible in Figure 7 , with populations in the oligotrophic subtropical gyres and the Caribbean Sea likely associated with the lower Chl-a range up to 0.1 mg m^-3 (see also Figure 2 ). Populations related to Chl-a between 0.1 and 0.3 mg m^-3, on the other hand, likely relate to those thriving in the subantarctic Malvinas Current and the highly productive Brazil-Malvinas confluence zone in the Argentine basin as well as those profiting from enhanced productivity along the subtropical front in the South Atlantic, the Azores front in the North Atlantic and upwelling associated with the equatorial currents. Therefore, the association between subtropical-tropical G. truncatulinoides populations and Chl-a likely reflects the species reproduction strategy and is determined by the winter primary production maximum at those latitudes (e.g., Lévy et al., 2005).

The distribution of N. incompta (the last of the five-top species) includes the mid-latitudinal Atlantic, the subarctic/subantarctic domains, and the Atlantic province of the Nordic Seas ( Figure 7 ), as previously reported (Johannesen et al., 1994; Boltovskoy et al., 1996; Schröder-Ritzrau et al., 2001; Schiebel and Hemleben, 2017). N. incompta is a species that tends to live below the Chl-a maximum (Field, 2004; Rebotim et al., 2017) and is assumed to be herbivorous as other non-spinose species (Hemleben et al., 1989). However, nothing is known about the food preferences of its genotype I. In the California Current EBUS, genotype II appears to graze on phytoplankton and prey on bacteria (Bird et al., 2018). The relationship of N. incompta with Chl-a ( Figure 5 ) reveals breaks at 0.09 and 0.4 mg m^-3. Compared with Figure 2 , the populations of the Norwegian Sea and the sub-Antarctic domain are probably the ones responsible for the lower Chl-a range since higher fluxes to the seafloor are reported to occur only after the blooms in the hemispheric summer to early fall (Schröder-Ritzrau et al., 2001; King and Howard, 2003). As a result, populations inhabiting regions with Clh-a concentrations above 0.4 mg m^-3, most probably will be those linked to the EBUS in the North and South Atlantic (Meggers et al., 2002; Salgueiro et al., 2008; Lessa et al., 2020).

The species that react the most to EKE variations are P. obliquiloculata, G. hirsuta, G. tenellus, G. siphonifera, and O. universa, reaching their highest importance above 2.4 cm² s^-2, whereas O. universa already starts to level out at 2.2 cm² s^-2 ( Figure 5 ). These species distribution ( Figure 7 ) at the latitudes of the equatorial current, counter- and undercurrent system and within the Caribbean Sea confirms that these species are strongly influenced by ocean currents and EKE. The occurrence of P. obliquiloculata and G. tenellus mainly in the area of the Gulf Stream ( Figure 7 ) is particularly interesting as it justifies the species relationship to EKE values above 2.1 cm² s^-2.

Although with different importance, P. obliquiloculata and G. hirsuta, show EKE thresholds at 1.5 and 1.8 cm² s^-2 and 1.3 and 2.3 cm² s^-2, respectively. Both species are non-spinose, tropical to temperate and live in the intermediate-depth range (100-300 m; Schmuker and Schiebel, 2002; Wilke et al., 2009; Rebotim et al., 2017; Jentzen et al., 2018). P. obliquiloculata is associated with the deep Chl-a maximum (Schiebel and Hemleben, 2017 and references therein), and its higher fluxes to the seafloor occur in winter (Salmon et al., 2015). The fact that the equatorial currents and their branches that feed into the Caribbean Sea and the Gulf Stream extend into subsurface levels and thus into the depth range of this species, may explain its reaction to EKE ( Figure 7 ). Although rare in the living fauna, this species has a high preservation potential in the sediment, which may result in an increase of its relative abundance in the sediment assemblage, and contribute to the GF identification over other tropical species like Globorotalia menardii and Neogloboquadrina dutertrei that are also known to be associated with the equatorial current system and the Gulf Stream (Fairbanks et al., 1980; Ufkes et al., 1998). G. hirsuta, on the other hand, feeds on diatoms and has an annual reproduction cycle similar to G. truncatulinoides, moving into subsurface waters in summer (Schiebel and Hemleben, 2017; Rebotim et al., 2017). In the North Atlantic, the species is associated with the southern side of the Azores front and Azores current (Storz et al., 2009; Rebotim et al., 2017), where its average living depth extends down to 400 m. Likewise, it exhibits higher abundances in the latitudes of the subtropical front in the South Atlantic ( Figure 7 ). So, the gradual rise in importance to EKE values between 1.3 and 2.3 cm² s^-2 can clearly be linked to its affinities to those major frontal systems in both hemispheres.

Spinose species G. tenellus lives in the upper mixed surface layer (e.g., Rebotim et al., 2017; Jentzen et al., 2018) and contributes significantly to the assemblage of the subtropical gyre in the South Atlantic (Lessa et al., 2020). The lower range EKE thresholds ( Figure 2 ) might reflect the high abundances observed between 15-20°S in the South Atlantic’s subtropical gyre ( Figure 7 ). This species also occurs with similar abundances, in the Canary Current EBUS and along the major currents in the southwestern quadrant, which might explain the relation to higher EKE values ( Figure 2 ).

Tropical to subtropical spinose, symbiont bearing, and omnivorous species G. siphonifera includes four genotypes with several subtypes and two morphotypes (Schiebel and Hemleben, 2017 and references therein), which are not distinguished in ForCenS. Not all of those subtypes occur in the Atlantic Ocean, and some of them are associated with particular oceanic regions or currents (Schiebel and Hemleben, 2017). That species genotypic variety seems to be reflected by the GF results. The linear response to 1.5 cm² s^-2 EKE appears to be related to the patch within the South Atlantic’s subtropical gyre or the subtropical gyres in general. Therefore, this group would include both cosmopolitan genotypes and those linked to the North Atlantic Drift, Azores and Canary Currents. The populations occurring after the break at 1.5 cm² s^-2 and the one at 2.1 cm² s^-2 are likely to encompass those living within the equatorial current and counter-current region, i.e., the same as observed for P. obliquiloculata, and those within the central and southern Caribbean Sea (Figure 7), i.e., mainly the tropical genotype associated with SST > 25°C (Schiebel and Hemleben, 2017). The third group, related to the highest EKE values, are composed of populations thriving within the Loop Current in the Gulf of Mexico and the Gulf Stream and potentially include the cosmopolitan genotype associated with SST > 17.5°C. The distribution of this group conforms to the high abundances observed by Fairbanks et al. (1980) in the slope waters and a cold-core ring, and the one solely observed in the North Atlantic Drift. Previous work has found that small cumulative changes in natural environmental parameters (e.g., shear stress), climate variables (e.g., temperature) and anthropogenic activities (e.g., bottom trawling) can lead to much larger responses than those predicted from linear effects on benthic biodiversity (Couce et al., 2020). Our results suggest a similar pattern also for plankton biodiversity.

Lastly, O. universa is a surface dwelling, mostly carnivorous species (e.g., Rebotim et al., 2017; Jentzen et al., 2018; Lessa et al., 2020). Three genotypes morphologically distinguishable based on shell thickness and pore size have been defined for this species (Morard et al., 2009; André et al., 2014; Marshall et al., 2015). Those geno-/morphotypes seem to have specific ecological preferences (Morard et al., 2009), which can be related to intervals in the cumulative importance curve ( Figure 4 ). Genotype III (Mediterranean Sea) is the one more widely distributed in the Atlantic Ocean and linked to nutrient-rich waters (Morard et al., 2009). This variant contributes to EKE response throughout the wide range of 0.3 to 1.6 cm² s^-2, with the lower end reflecting the central part of the South Atlantic’s subtropical gyre ( Figures 2 , 5 , 7 ). Due to its preference by nutrient-rich waters, we relate the higher O. universa abundances in the EBUS regions off western Iberia, NW and SW Africa ( Figure 7 ) to the occurrence of this genotype and an EKE range from 1 to 1.6 cm² s^-2. Like for other species, the response at the upper end of the EKE spectrum, i.e., values above 1.6 cm² s^-2, is likely driven by the occurrences of genotypes I and II in the Caribbean Sea (genotype I), the equatorial current system and the Gulf Stream, where both genotype I and III are found (Morard et al., 2009). Genotype II is associated with the oligotrophic waters of the Sargasso Sea, and probably also contributes to the EKE range between 1 and 1.6 cm² s^-2 ( Figures 2 , 7 ).

In summary, our study reveals that planktonic foraminifera diversity in the Atlantic Ocean, and on a large scale, is not only determined by SST but also by EKE (currents) and Chl-a (photosynthetic biomass). Furthermore, the GF analysis allowed the redefinition of the location of the biogeographic transition zones (e.g., Bé and Tolderlund, 1971; Kučera, 2007), and identifies the species that most react to the variation of the different parameters. These results prompt the interest to use diversity indexes and the parameter-related species to attempt the definition of functions to reconstruct past EKE and Chl-a. Besides, the fact that large-scale ocean circulation is projected to change in response to climate warming calls for the use of diversity indexes to monitor such variations. Such monitoring could be achieved by combining the sampling of microplankton at specific locations along the encountered diversity gradients. Changes in latitude and/or longitude of those diversity gradients will reflect changes in the prevailing oceanographic circulation.