Marine biodiversity and the chessboard of life

Species richness is greater in places where the number of potential niches is high. Consequently, the niche may be fundamental for understanding the arrangement of life and especially, the establishment and maintenance of the well-known Latitudinal Biodiversity Gradient (LBG). However, not all potential niches may be occupied fully in a habitat, as measured by niche vacancy/saturation. Here, we theoretically reconstruct oceanic biodiversity and analyse modeled and observed data together to examine patterns in niche saturation (i.e. the ratio between observed and theoretical biodiversity of a given taxon) for several taxonomic groups. Our results led us to hypothesize that the arrangement of marine life is constrained by the distribution of the maximal number of species’ niches available, which represents a fundamental mathematical limit to the number of species that can co-exist locally. We liken this arrangement to a type of chessboard where each square on the board is a geographic area, itself comprising a distinct number of sub-squares (species’ niches). Each sub-square on the chessboard can accept a unique species of a given ecological guild, whose occurrence is determined by speciation/extinction. Because of the interaction between the thermal niche and changes in temperature, our study shows that the chessboard has more sub-squares at mid-latitudes and we suggest that many clades should exhibit a LBG because their probability of emergence should be higher in the tropics where more niches are available. Our work reveals that each taxonomic group has its own unique chessboard and that global niche saturation increases when organismal complexity decreases. As a result, the mathematical influence of the chessboard is likely to be more prominent for taxonomic groups with low (e.g. plankton) than great (e.g. mammals) biocomplexity. Our study therefore reveals the complex interplay between a fundamental mathematical constraint on biodiversity resulting from the interaction between the species’ ecological niche and fluctuations in the environmental regime (here, temperature), which has a predictable component and a stochastic-like biological influence (diversification rates, origination and clade age) that may alter or blur the former.


Introduction
Understanding how life is arranged has occupied scientists for centuries [1][2][3]. Despite many hypotheses or theories [1,3,4] however, there is still no consensus on the causes of some marine ecogeographic patterns such as the Latitudinal Biodiversity Gradient (LBG). Our recent MacroEcological Theory on the Arrangement of Life (METAL) proposes that current patterns in marine biodiversity is driven by the interaction between the ecological niche of a species (sensu Hutchinson [5]) and the local environmental regime (mean and variability), which subsequently propagates at the community level [6][7][8][9].
The concept of the niche sensu Hutchinson [5] represents the combination of environmental tolerances and resources required by an organism. Hutchinson [5] conceptualised this notion with the so-called n-dimensional hypervolume, in which n ideally corresponds to all environmental factors. Using the ecological niche allows consideration of underlying processes (genetic and physiological) that are difficult to identify and parametrise for a large number of species. Interaction between the niche and climatic and environmental changes propagate from the species to the community and ecosystem levels, and are detectable from the smallest ecosystems to the whole ecosphere [6].
The METAL theory (hereafter referred to as METAL) offers a way to make testable ecological and biogeographical predictions to understand how life is organised and how it responds to global environmental changes, including climate change [6,10]. At the organismal level, METAL explains (i) how species are distributed in space (species range) and time (phenology) [9], (ii) local changes in species' abundance [10][11][12][13], (iii) species responses to climate change in time and space, including the inconsistency in the relationships between climate and species' abundance [12], (iv) phenologic and (v) biogeographic shifts [9]. At the community level, METAL explains (i) some biogeographic laws or ecogeographic patterns (e.g. Rapoport's pattern, infrequency law) [7], (ii) large-scale patterns in biodiversity [7], (iii) changes in biodiversity in space and time [8], (iv) long-term community shifts, including environmentallyinduced abrupt community shifts (also called regime shift) [14] and (v) species and community sensitivity and vulnerability to climate change [8]. In addition, METAL unifies all these different biological and ecological responses to environmental changes in space and time by a unique process: the interaction between the species ecological niche and fluctuations in the environmental regime (S1 and S2 Text).
When we applied METAL at a global scale and community level, it provides an assessment of the ecogeographic patterns in biodiversity [7,8]. Not all potential niches are fully occupied in real communities however, a phenomenon called niche vacancy or saturation [15]. Therefore, METAL can be seen as providing an estimation of the species richness at saturation. Because niche saturation is the ratio of observed species richness of a given taxonomic group on the modeled (from METAL) theoretical number of species at saturation, the procedure enables large-scale patterns in niche saturation to be assessed when applied in conjunction with data for observed species richness of different taxonomic groups [4]. Investigating the LBG and ecogeographic patterns in niche saturation jointly in this way may provide important clues to the key factors and processes that influence the arrangement of marine life.
Here, using METAL [6][7][8][9] and global-scale observations [4], we reconstruct global oceanic biodiversity and examine large-scale ecogeographic patterns in niche saturation for six taxonomic groups (foraminifers, euphausiids, oceanic sharks, tuna/billfish, cetaceans, pinnipeds). We investigate what factors or processes may be at the origin of those ecogeographic patterns to help identify the respective influence of ecological and evolutionary factors. The influence of diversification processes such as speciation and extinction is studied indirectly by the examination of niche saturation, which integrates (theoretical) expected richness at saturation and observed species richness of the six taxonomic groups considered in this study. We investigate whether the mean degree of niche saturation of a taxonomic group depends upon organismal complexity and discuss briefly, how climate change may restructure current patterns of marine biodiversity. Based on those findings, we propose a new hypothesis using the analogy of a chessboard that integrates both ecological and evolutionary processes in a single theoretical framework to explain the arrangement of life in the superficial (epipelagic) part of the ocean.

Physical data
Annual SSTs originated from the dataset ERSST_V3 (1850-2014). This dataset is derived from a reanalysis based on the most recently available International Comprehensive Ocean-Atmosphere Data Set (ICOADS). Improved statistical methods have been applied to produce a stable monthly reconstruction, on a 1˚x 1˚spatial grid, based on sparse data [16]. Both mean annual SST and SST variability (i.e. the coefficient of variation of SSTs) were calculated on the basis of monthly SSTs for the period 1850-2014; a total of 165 years x 12 months = 1980 values were therefore used to assess mean annual SST and SST variability in each geographical cell. SST variability was subsequently scaled between 0 and 1 as follows: With g Ã the scaled variable and g SST variability. We used the coefficient of variation of SSTs instead of the standard deviation to remove the positive influence of the mean of SSTs on the standard deviation. We prefer this parameter to the range of SSTs, which may be affected too much by extreme values at a global scale. The coefficient of variation of SSTs have been frequently used in the past as an index of temporal variability of SSTs [17][18][19].

Biological data
The dataset of observed biodiversity was provided by Dr Derek Tittensor, Dalhousie University [4]. The data were compiled from empirical sampling data (foraminifers and bony fish) or from expert-verified range maps (euphausiids, oceanic sharks and marine mammals) encompassing many decades of records. The data was originally gridded on a 880-km equal-area resolution grid [4]. The taxonomic groups we chose had a species richness that reflected well its biodiversity (i.e. species richness greater than 50% of known species; see their Table 1): foraminifers (Chromalveolata Rhizaria, 39 species; 88% of known species), euphausiids (Crustacea Malacostraca, 86 species; 100% of known species), oceanic sharks (Chondrichthyes, 27 species; 100% of known species), tuna/billfish (Osteichthyes, 12 species; 63% of known species), cetaceans (oceanic marine mammals, 81 species; 96% of known species) and pinnipeds (sea/land marine mammals, 36 species; 100% of known species); this last taxonomic group shows a wellknown inversed latitudinal biodiversity pattern [4], which is explained when biodiversity is decomposed into groups with different levels of steno/eurythermy [7]. Spatial representations of these data are available in Tittensor and colleagues (their Fig 1) [4].
scale biodiversity patterns using a model based on the MacroEcological Theory on the Arrangement of Life (METAL) [6,8,9,14](S1 Text). This theory explains how life is arranged in the sea and how changing environmental conditions alter biological arrangements in space and time at different organisational levels (e.g. species, community, ecosystem), allowing precise predictions to be tested. In particular, we use the model recently developed in Beaugrand and colleagues [8] that estimates the implications of climate change on biodiversity for the past, present and future. This model is based on the use of thermal niches and temperature. While temperature is the most frequently used environmental parameter in theories or empirical studies we are nevertheless aware that the niche is a multi-dimensional concept [20] and that other environmental parameters affect species [2](S2 Text). While additional environmental parameters are probably necessary to provide better estimations on smaller spatial scales we assume that temperature alone provides a correct estimation of theoretical biodiversity patterns on a global scale; this assumption was tested in our previous studies that showed temperature explains a large proportion of variance of both past and contemporaneous patterns in biodiversity for several taxonomic groups [7,8,10]. More information on the METAL and its applications can be found in http://metaltheory.weebly.com/. Calculation of potential niches at saturation (variable 1 in Fig 1): a fundamental assumption of METAL is that each species has a unique ecological niche [21], which determines to a large extent its spatial (biogeography) and temporal (phenology and year-to-year to long-term changes) distribution [6][7][8][9][10]12]. Providing this assumption is true, we can create a pool of pseudo-species, characterized by all thermal optima from low to high and thermal breadth from strict stenotherms to universal eurytherms. Each species, whether ectothermic or endothermic, has evolved anatomical structures and life history traits that allow them to live within a certain range of thermal conditions. In other words, all species have their own thermal niche and so we expect the model will apply to both ectotherms and endotherms. All pseudo-species colonise a given region of the global ocean providing they can withstand the local thermal (annual SST) regime (Fig 2). By this way, we reconstruct the arrangement of life in the oceans and seas and the number of maximal species' niches at saturation is approached by using the total number of pseudo-species per geographical square [7,8].
More specifically, pseudo-species richness S t was assessed by using the METAL model of Beaugrand and colleagues [8]. The response curve of the abundance L of a pseudo-species s in a given site i and time j to change in annual SSTs (1960-2013, S2 Text) was modeled by the following Gaussian function [22]: With L i,j,s the expected abundance of a pseudo-species s at location i and time j (years of the period 1960-2013); c s the maximum value of abundance for species s fixed to one; λ i,j the value of temperature at location i and time j; ϕ s the thermal optimum and h s the thermal breadth for species s (this parameter corresponds to the standard deviation in the normal distribution). Once the niche was modeled, the expected abundance of pseudo-species was determined in space or time by linear interpolation from annual SST for a given year or time period and geographical square. A large number of pseudo-species was created with ϕ s varying between -1.8˚C and 40˚C by 0.1˚C increments and h s varying between 1.1˚C and 10˚C by increments of 0.05˚C. As we showed previously [7], this implementation generates a Mid-Domain Effect (MDE) [23] in the Euclidean space of the niche, which influences the geographical patterns in modeled species richness and niche saturation because more thermal niches overlap in warmer than in colder environments [7]. Our model assumes that no species of the same ecological guild can have exactly the same thermal niche, although relaxing this constraint would not alter our conclusions significantly. In contrast to our framework, some biodiversity models (empirical or geometrical) do not assume niche exclusion [24,25].
After removing randomly half of the created pseudo-species to reduce calculation time [8], a total of 19609 pseudo-species was able to colonise the whole ocean [8]; we showed previously that this does not affect the ecogeographic biodiversity patterns (Beaugrand and colleagues [7], see Scenarios 18-19 and 27 in their S1 Table). The map of biodiversity, which can also be interpreted as a map of the maximal number of species' niche at saturation is displayed in Fig 3. Data are available S1 Data.
Degree of eurythermy and thermophily (variables 2 and 3 in Fig 1): from the METAL model [8], we calculated the degree of eurythermy of a pseudo-community (i.e. weighted mean of the degree of eurythermy for all pseudo-species) and the degree of thermophily of each pseudo-community (i.e. weighted mean of the degree of thermophily for all pseudo-species) per geographic square. These two properties were calculated to better understand ecogeographic patterns in the degree of niche saturation of the six taxonomic groups. We represented the degree of eurythermy using 10 categories, each representing a decile; for example, the first and last categories corresponded to the first and last deciles.
Differences in biodiversity between the Last Glacial Maximum/mid-Pliocene and today (variables 4 and 5 in Fig 1): data on quantitative changes in biodiversity expected between the Last Glacial Maximum (LGM,~20,000 years ago) and today, and between the mid-Pliocene (mid-Piacenzian, 3.264-3.025 Ma) and today, originated from the model of Beaugrand and colleagues [8]. The index of quantitative change was the sum of the differences between theoretical abundance (standardized between 0 and 1) at two time periods weighted by species richness common for the two time periods. This index was calculated for both periods for pseudo-species having an abundance higher than 0.1 (a threshold below which the species was considered to be absent) for both periods. For a given geographical cell, quantitative index F was calculated as follows: With p species richness common to the two time periods, a i the abundance of species i at the second time period and b i , the abundance of species i at the first time period. In this model, the maximum abundance of all species reaches 1. Those gridded data were used as a proxy for environmental perturbations induced by long-term climate change.
All biological variables assessed from METAL were calculated on a spatial grid of 2˚longitude x2˚latitude from 178˚W to 180˚E and from 88˚S to 88˚N.
Step 2: Spatial interpolation of mean and variability in SSTs (variables 6 and 7 in Fig  1). Mean and variability in SSTs were spatially interpolated on a grid of 2˚longitude x2T he chessboard of life latitude from 178˚W to 180˚E and from 88˚S to 88˚N using the inverse squared distance method [26].
Step 3: Estimation of mass-corrected rate of evolution from mean SSTs (variable 8 in Fig 1). Although the Metabolic Theory of Ecology (MTE) remains debated [27], we used a model originating from this framework to calculate the mass-corrected evolutionary rate [28]. The mass-corrected evolutionary rate (here the rate of DNA substitution; mitochondrial DNA) D was calculated from annual SST [29]: With α = -0.74eV, k the Boltzmann's constant (8.62 x 10 −5 eV.K -1 ), V the mean annual SST (K) and μ = 30.28 [29]. Values of α and μ originated from the linear model of Gillooly and colleagues based on mitochondrial DNA data [29]. Those data were utilized to examine the relationships between niche saturation and a proxy for evolutionary rates.
Step 4: Spatial interpolation of observed species richness (variables 9-14 in Fig 1). Observed species richness for foraminifers, euphausiids, tuna/billfish, oceanic shark, cetaceans and Pinnipeds were also spatially interpolated on a grid of 2˚longitude x2˚latitude from 178˚W to 180˚E and from 88˚S to 88˚N using the inverse squared distance method [26].
Step 5: Correlation between expected and observed large-scale biodiversity patterns. Linear correlations between expected and observed biodiversity patterns were calculated for each taxonomic group and between expected biodiversity and both mean and variability in SSTs (Fig 1). To account for spatial autocorrelation in the geographical pattern of species richness (two dimensions), the degrees of freedom were recalculated to indicate the minimum number of samples (n Ã ) needed to maintain a significant relationship at p = 0.05 [19,30,31]. The smaller n Ã , the less likely is the effect of spatial autocorrelation on the probability of significance [7]. We preferred this technique to others (e.g. a technique based on the calculation of the Moran's index or classical semi-variograms) that are based on the assumption of isotropy, which is often violated, as shown on the diversity of North Atlantic calanoid copepods calculated by using (local) Point Cumulative Semi-Variograms [32].
Step 6: Degree of niche saturation. The degree of niche saturation can be defined as the percentage of occupied niches (sensu Hutchinson [5]) in a particular community (see Fig 4 for an example). For a given taxonomic group or an ecological guild, niche saturation can be assessed by calculating the number of species observed in a given location out of the maximum number of species that can coexist in a site [15]; this way of calculating saturation is similar to what was proposed by MacArthur and Wilson in their equilibrium theory of insular zoogeography [33] (page 373). For each of the six taxonomic groups, the degree of niche saturation Q x of a taxonomic group x for a given geographical cell was determined by calculating the ratio of the observed species richness S o x on the maximum species' niche at saturation S t y (theoretical biodiversity assessed from METAL in each geographical cell and seen in Fig 3) The total number of niches (or species) at saturation considered in our model was 19609. The oceans contain between~200000 and~250000 already inventoried species [34,35], which represents between 7.8% and 9.8% of all described species and between 0.9 and 3.9% of all estimated species (between 0.5 and 2.2 million species [35,36]). The taxonomic groups we considered represented between 0.06% (12 species of tuna/billfish) and 0.44% (86 species of euphausiids) of the modeled biodiversity (19609). Results from the calculation of Q for the 6 taxonomic groups (foraminifers, euphausiids, oceanic sharks, tuna/billfish, cetaceans and pinnipeds) are shown in Fig 5. To ensure that maps of niche saturation were comparable, we replaced niche saturation values by rank and scaled them between 0 (complete unsaturation) and 1 (full saturation). From Eq (5), we see that S t y did not vary among taxonomic groups, only among geographical cells. In this analysis, we therefore examined ecogeographic patterns in relative saturation among groups and not directly their values of saturation (see Step 13 for the comparison of the values of niche saturation).
Step 7: Examination of the ecogeographic patterns in the degree of niche saturation of the six taxonomic groups (variables 15-20 in Fig 1) according to annual SSTs and SST variability (variables 6-7 in Fig 1). The degree of niche saturation of the six taxonomic groups ( Fig 5) was assessed as a function of mean SSTs and SST variability (i.e. the coefficient of variation of SSTs). For this analysis, we used the matrices constructed to map the niche saturation of the six taxonomic groups in Fig 5 (Step 6) and the corresponding matrices of mean and variability in SSTs interpolated in Step 2 (Fig 1). A total of 8 matrices were therefore analysed and each corresponded to a grid of 2˚longitude (from 178˚W to 180˚E) x 2˚latitude (from 88˚S to 88˚N). The median of the degree of niche saturation was calculated for each interval of mean annual SST between -1˚C and 30˚C by increment of 0.25˚C and each interval of SST variability between 0 and 1 by increment of 0.05; values of mean annual SST below -1˚C, not shown, were under-represented and most were highly saturated so that they appear as outliers in our subsequent analyses (principal component analyses). The median was only calculated when the minimum number of values of niche saturation inside a given range of mean and variability in SSTs was higher than 10. To minimize the influence of outliers on the graphics, we calculated the 5 th and 95 th percentiles and fixed values below or above to these thresholds. Data were then standardized between 0 (lowest value of niche saturation) and 1 (highest value of niche saturation). The last two transformations were not applied when data were analysed (e.g. polynomial regression, correlation and principal component analyses; see below). The results of this analysis were represented on a scatterplot for each taxonomic group, with bullets that had a colour and a size proportional to their value ( Fig 6); small blue bullets indicated small median values of niche saturation for a given increment of mean and variability in SSTs and conversely, for large red bullets. Data are available in S1 Data.
Step 8: Examination of patterns in explanatory variables (variables 1-5 and 8 in Fig 1) according to annual SSTs and SST variability (variables 6-7 in Fig 1). A similar procedure was applied for explanatory variables for a mean SST between -1˚C and 30˚C (Fig 7): (1) degree of eurythermy of each modeled (METAL) pseudo-community, (2) degree of thermophily of each modeled (METAL) pseudo-community, magnitude of changes in modeled (METAL) pseudo-species richness between (3) the LGM or (4) the mid-Pliocene and today, (5) the maximal expected (METAL) number of niches (or species) at saturation and (6) the mass-corrected evolutionary rate assessed from the MTE. Note that the degree of eurythermy and thermophily was not assessed from observed data because we have only information on observed species richness and not individual species range (see Step 12). This peculiar chessboard has a number of squares that corresponds to the number of geographic squares that divides the oceans. In each square, a sub-square represents a specific ecological niche sensu Hutchinson [20] and their number in a square illustrates maximum species' niches at saturation. Note that the area of a sub-square should not be not interpreted as related to the size of a niche. Future research may show that it is related to the number of individuals within a species, although this also depends upon resource availability and environmental conditions. Sub-squares are not referenced in the geographic space and represent a potential niche that can be taken by a species (not an individual). In this study, the number of sub-squares is fixed by the MacroEcological Theory on the Arrangement of Life (METAL). Q means niche saturation of a given community. S 0 and S T represent observed species richness and maximum species' niches at saturation (expected from METAL), respectively. The species, the chess piece on the chessboard, can appear and disappear by speciation (and immigration) and (local or global) extinction, respectively. The different chess pieces illustrate the differences in species' life history traits such as dispersal (e.g. rooks and bishops move in a different way) and degree of importance or trophic status (e.g. a queen has not the same importance than a pawn). Although marine biodiversity should increase from poles to the tropics, clade origination and life history traits of a specific taxonomic group may blur this pattern. The increase in diversification rates from cold to warm regions or the increased probability of clade origination in the tropics where the number of niches are greater explain an augmentation in niche saturation between temperate and permanently stratified regions because the temperate-tropic difference in maximum species' niches at saturation is small nowadays. Biogeographic shifts can rearrange the sub-squares in a geographic square and when occupied the chess pieces (species), providing that there is no geographical barrier or that dispersal is sufficiently elevated. In case of a global cooling, the Latitudinal Biodiversity Gradient (LBG) reinforces because of the movement of the sub-squares and associated species towards the equator. In case of global warming, the LBG is attenuated because of the movement of sub-squares and corresponding chess pieces polewards. https://doi.org/10.1371/journal.pone.0194006.g004 The chessboard of life Standardization between 0 and 1 was not needed. For this analysis, we used the matrices interpolated in Steps 1-3 (Fig 1). Therefore a total of 8 matrices (the 6 explanatory variables above and mean and variability in SSTs) was analysed, each having a size corresponding to a grid of 2˚longitude (from 178˚W to 180˚E) x 2˚latitude (from 88˚S to 88˚N). As in step 7, results of the procedure were represented for each explanatory variable by means of a scatterplot (Fig 7). Bullets had a colour and a size proportional to the value of the variable; small blue bullets indicated small median values of a given explanatory variable for a given increment of mean and variability in SSTs and conversely for large red bullets (Fig 7). The chessboard of life Step 9: Examination of the patterns in niche saturation within the 2-dimensional Euclidean space determined by mean SST and SST variability. A 2-order polynomial regression was calculated to test whether patterns in niche saturation were detected within the 2-dimensional Euclidean space determined by mean SST and SST variability. R 2 , adjusted R 2 and RMSE (Root Mean Squared Error) were examined. The basic polynomial model used for The chessboard of life all regression analyses was: With β the estimation of the degree of niche saturation of a taxonomic group (foraminifers, euphausiids, oceanic sharks, tuna/billfish, cetaceans, pinnipeds) or the value of an explanatory variable (degree of eurythermy, degree of thermophily, LGM/today quantitative changes, mid-Pliocene/today quantitative changes, theoretical number of niches, mass-corrected rate of evolution), m mean SST and g SST variability. The parameters δ 1 to δ 5 and U were constants assessed from the regression analysis. A total of 12 polynomial regression analyses was performed, one for each variable considered in Figs 6 and 7 (S1 and S2 Tables).
Step 10: Correlations between changes in niche saturation and explanatory variables according to mean SST and SST variability. We calculated linear correlation coefficients between changes in niche saturation for the six taxonomic groups and explanatory variables in Figs 6 and 7 (S3 Table). Those correlations were based on (i) mean SSTs ranging from -1˚C to 30˚C (185 degrees of freedom for all variables but LGM (157)) and (ii) mean SSTs from 2.75˚C to 30˚C (149 degrees of freedom for all variables).
Step 11: Principal Components Analyses (PCAs). Principal Components Analyses (PCAs) were performed to investigate globally the relationships between the degree of niche saturation of the six taxonomic groups and their relationships with explanatory variables (Fig  8). The first PCA (Fig 8A-8C) was applied on the table [category of mean SST and SST variability] x 6 variables, using all categories of mean SST (i.e. mean SST between -1˚C and 30˚C by increment of 0.25˚C) and SST variability (i.e. SST variability between 0 to 1 by increment of 0.05)(see Figs 6 and 7). The second PCA (Fig 8D-8F) was applied on the table [category of mean SST and SST variability] x 6 variables, for categories of mean SST higher than 2.75˚C. The second PCA was performed to examine more closely the relationships between niche saturation of the six taxonomic groups and explanatory variables once some outliers were removed. Eigenvectors were normalised as follows [37]: With U Ã being the matrix of normalised eigenvectors, U the matrix of eigenvectors and Λ the diagonal matrix of eigenvalues. Thus, the variables in the space of eigenvectors represented the linear correlation with the first and the second principal components [38]. We used this mathematical property to add supplementary variables (all explanatory variables), simply by calculating the linear correlation (Pearson correlation coefficient) between these variables and the first two principal components. Thus, explanatory variables had no weight in the calculation of the principal components. The first two normalized eigenvectors and corresponding principal components were represented (Fig 8).
Step 12: Observed species richness as a function of temperature. A bar diagram was performed for each of the 6 taxonomic groups (foraminifers, euphausiids, oceanic sharks, tuna/billfish, cetaceans and pinnipeds), to examine the mean observed species richness per annual SST, between -2˚C and 30˚C and by 0.5˚C increment (S1 Fig and S4 Text). This analysis was performed to estimate the mean thermal optimum and thermal range of each taxonomic group. The highest values in a bar diagram indicated mean thermal optimum T o , which was assessed by calculating the average of the first 5 highest values of species richness in the diagram. The degree of equity in a bar diagram was used as an indication of the thermal range E o of a taxonomic group; here this was assessed by calculating the coefficient of variation of observed species richness.  Fig 3), pseudo-community thermophily (Thermo), LGM/today biodiversity changes (LGM), mid-Pliocene/today biodiversity changes (Plio), theoretical number of niche/species at saturation (theoS) and mass-corrected evolutionary rate (Evol). https://doi.org/10.1371/journal.pone.0194006.g008 The chessboard of life Step 13: Global average degree of niche saturation for four ecological guilds. First estimations based on described species: we used the taxonomic groups as representative of some ecological guilds (or functional types) to assess their mean niche saturation as an absolute value (and not relative as in Fig 4); foraminifers (39 species) were assumed to be representative of protozooplankton, euphausiids (86 species) of metazooplankton, both tuna/billfish (12 species) and oceanic sharks (27 species) of pelagic oceanic fish and both cetaceans (81 species) and pinnipeds (36 species) of marine mammals (Fig 9).
Let S o x be the species richness of a given geographical cell for a given taxonomic group x and S t y the maximum (theoretical) species richness (or total number of niche at saturation) of a given geographical cell determined by the METAL model (o for observed species richness and t for theoretical pseudo-species richness). S o w can be scaled for a given functional group or ecological guild w by dividing the total species richness of a taxonomic group x and by multiplying by the total species richness of the functional group or ecological guild w, as follows: x is the species richness of a given taxonomic group; x = 39 species for foraminifers, x = 86 species for euphausiids, x = 12 species for tuna/billfish, x = 27 species for oceanic sharks, x = 81 species for cetaceans and x = 36 species for pinnipeds. The parameter w was first fixed to 1350 for protozooplankton and 5500 for metazooplankton [39], 1068 for pelagic species [40] and 117 for marine mammals [4]. And similarly, S t y can be scaled to correspond to the total species richness z of the whole species richness of all functional groups considered in the analysis. S t y was divided by y (here 19609) and multiplied by z: z was first equal to 1350 (protozooplankton) + 5500 (metazooplankton) + 1068 (pelagic fish) + 117 (marine mammals) = 8035 species. The corrected degree of niche saturation Q w of a trophic guild or functional type w is given by extrapolating species richness to a functional group and scaling maximum species richness to the total number of species considered in the analysis: With Q x the degree of niche saturation of a given taxonomic group x calculated from Eq (5). Uncertainties related to undescribed and unknown species: the estimations above are highly dependent upon the true species richness of an ecological guild. If estimations of global mean niche saturation of a given ecological guild are based on an underestimated number of species, this will lead to an underestimation in niche saturation. The conclusion could be that communities are highly undersaturated whereas they are not. Therefore, we also corrected our first estimations (see previous section) in different ways, depending upon the ecological guild (S3 Text).
For all ecological guilds, we corrected our estimations on mean niche saturation by adding species remaining to be discovered (S3 Text). Two corrections were applied for plankton, one based on the work of Appeltans and colleagues [36] and the other based on the work carried out as part of TARA expedition [41]. For protozooplankton, w (i.e. the total number of species in the ecological guild or functional type; here protozooplankton) was fixed to 2363 (first estimation [36]) and 30192 species (second estimation [41]) and z (i.e. the total number of species in all ecological guilds or functional types) to 9048 (first estimation [36]) and 36877 species (second estimation [41]) in Eqs (8)- (10). For metazooplankton, w was fixed to 8945 (first estimation [36]) and 22500 (second estimation [41]) species and z to 11480 (first estimation [36]) and 25035 (second estimation [41]) in Eqs (8)- (10). For pelagic fish and marine mammals, our estimations were based exclusively on the work of Appeltans and colleagues [36]. For fish, w was fixed to 1314 species and z to 8281 and for mammals w was fixed to 125 species and z to 8043 in Eqs (8)- (10). In all estimations, x (i.e. the species richness of a given taxonomic group) = 39 species for foraminifers, x = 86 species for euphausiids, x = 12 species for tuna/billfish, x = 27 species for oceanic sharks, x = 81 species for cetaceans and x = 36 species for pinnipeds. And y = 19609 pseudo-species richness for all estimations. The chessboard of life Histograms of the values assessed for each geographical square was then performed for each guild (Fig 9): (i) Protozooplankton, (ii) Metazooplankton, (iii) pelagic oceanic fish and (iv) marine mammals. The mode of each histogram, as well as the median of all values of niche saturation, were assessed to provide a global mean niche saturation value (Table 1).

Results and discussion
Using METAL (S1 Text) we can model ecogeographic biodiversity patterns by creating a pool of theoretical pseudo-species that each has a distinct thermal optimum and breadth. Each theoretical pseudo-species can occupy a given region of the ocean so long as they can withstand local thermal fluctuations (Methods, Fig 2 and S1 Text). In each geographical cell of a gridded ocean, pseudo-species form a pseudo-community that is characterized by its richness. We found that large-scale biodiversity patterns modeled by Gaussian niches (Fig 3) were very similar to biodiversity patterns modeled by rectangular niches (r = 0.99; p<0.01; n = 9927, n Ã = 4) [7]; this result suggests that the shape of the niche does not alter biodiversity patterns significantly. Spatial patterns in expected biodiversity were correlated positively and strongly with mean SSTs (r = 0.89,p<0.01,n = 10578,n Ã = 5) and negatively and strongly with long-term (1850-2014) monthly variability in SSTs (r = -0.74,p<0.01,n = 10578,n Ã = 8). This joint influence of mean and variability in SST was also detected in a previous model performed at a weekly temporal resolution where both eurythermic and thermophilic species were important in generating the maximum of biodiversity that occurs in the mid-latitude regions [7,8]. Although the Mid-Domain Effect (MDE), as defined originally in the geographical space [23], did not influence the ecogeographic biodiversity patterns in our study because the peak of biodiversity did not take place at the equator (i.e. the mid-range of the domain), we showed previously that a MDE occurring in the Euclidean space of the thermal niche affects ecogeographic biodiversity patterns [7].
We found positive correlations between expected and observed ecogeographic biodiversity patterns for foraminifers (r = 0.87,p<0.01,n = 9853,n Ã = 5), euphausiids (r = 0.79,p<0.01, n = 9853,n Ã = 7), oceanic sharks (r = 0.76,p<0.01,n = 9853,n Ã = 7), tuna/billfish (r = 0.82, p<0.01,n = 9853,n Ã = 6), and cetaceans (r = 0.68,p<0.01,n = 9853,n Ã = 9), while pinnipeds were negatively correlated (r = -0.82,p<0.01,n = 9853,n Ã = 6) [7]. The positive correlations suggest the existence of a deterministic component to how marine life is arranged [6]. As a consequence Fig 3 can be seen as a type of chessboard (Fig 4) where each geographic square has a distinct number of sub-squares, determined by the maximum number of species' niches at saturation, and where each sub-square can hold only one species of the same trophic guild, following Gause's principle of competitive exclusion [21]. In other words, the chessboard constrains the number of species that can co-occur in each geographical square; note that the size of a sub-square in Fig 4 is arbitrary (i.e. fixed by the number of sub-squares within a geographical square in the example) and should not be related to the size of the corresponding ecological niche . Fig 3 shows that the number of sub-squares is higher in the tropics and lower polewards. The number of actual species out of the number of sub-squares determines niche saturation in each square for a taxonomic group (Fig 4) [33].
Large-scale ecogeographic patterns in niche saturation have not been shown on a global scale before. The removal of the deterministic effect fixed by the underlying chessboard (here assessed by calculating expected biodiversity from METAL) may provide interesting clues on other mechanisms that affect the arrangement of life. We show that different ecogeographic patterns of niche saturation exist among the 6 taxonomic groups, which suggests the existence of a particular chessboard for each group that might originate from taxon-specific diversification history [42,43] (Fig 5). Although it is commonly assumed that niche saturation increases towards the equator (i) because evolutionary rates have been hypothesised to increase from cold to warm regions [29,44] and (ii) because of the presence of strong climate-induced environmental perturbations in extra-tropical regions that limit species richness [15,45], we found that niche saturation was highest towards the poles for all groups but fish (Fig 5). Those results show that the few sub-squares available in polar regions are frequently occupied, which also suggests that low polar biodiversity should not always be attributed to low diversification rates (origination minus extinction) [46,47], but rather to smaller maximum number of species' niches at saturation that limits local biodiversity fundamentally. When we excluded polar regions, niche saturation of all groups but mammals was higher over permanently stratified regions.
To better understand the origin of ecogeographic patterns in niche saturation, we represented niche saturation and key ecological properties as a function of mean SST above -1˚C and long-term (1850-2014) monthly variability in SSTs; these last two variables are important in explaining marine biodiversity patterns (S2 Text) [4,17]. All patterns in niche saturation and explanatory variables were significant (S1 and S2 Tables). Fig 6 shows patterns in niche saturation as a function of mean SST and long-term (1850-2014) monthly variability in SSTs. Higher niche saturation (red bullets in Fig 6) was observed towards colder regions for all taxonomic groups but oceanic sharks (Fig 6) and it coincided with (i) higher modeled pseudocommunity eurythermy that provides a better resistance to climate [8], (ii) lower modeled pseudo-community thermophily, (iii) higher influence of long-term climate change on modeled biodiversity, (iv) smaller maximum expected number of species'niches and (v) low evolutionary rates as modeled from the MTE (Fig 7). For long-term (1850-2014) monthly variability in SSTs above 0.7, many taxa exhibited low niche saturation, possibly because of conditions that are too harsh [6]. A focus on saturation for SSTs ranging from 2.75˚C to 30˚C (excluding Polar regions) revealed a second pattern, an increase in saturation from cold-temperate to tropical waters for plankton and fish (see small panels in Fig 6). Despite having more sub-squares in the tropics, the density of occupation of the chessboard is higher in tropical than temperate zones for plankton and fish; this is probably because of greater net tropical diversification rates [46,47] or faster species turnover in extratropical regions [48] (Fig 7C, 7D and 7F; small panels).
We performed two Principal Components Analyses (PCAs) on the matrix of niche saturation [categories of mean SST and its long-term monthly variability in SSTs x 6 taxonomic groups] for mean SSTs (i) above -1˚C and (ii) 2.75˚C and with ecological parameters described in Fig 7 as supplementary variables (Fig 8). Plankton and marine mammals were correlated negatively with maximum species' niches at saturation and fish were either uncorrelated (tuna/billfish) or correlated positively (Fig 8C; S3 Table). Plankton and mammals were significantly positively correlated with pseudo-community eurythermy, which suggests they may be resistant to strong thermal fluctuations (Fig 8C; S3 Table); while this result was confirmed by the positive correlation with modeled biodiversity changes between the Last Glacial Maximum (LGM) and today it contrasts with the small thermal range covered by pinnipeds (S1 Fig and  S4 Text). Oceanic sharks were significantly correlated positively with pseudo-community thermophily and corresponding evolutionary rates assessed from the MTE, in contrast to tuna/billfish (Fig 8C; S3 Table). Lack of congruence among taxonomic groups may be explained by their respective diversification histories. For example, pinnipeds originate from Arctoid carnivores [25][26][27] Ma that were present in the cold regions of the North Pacific [49], which may also explain the contrast with other groups and the negative correlation between predicted (pseudo-species richness at saturation) and observed biodiversity (r = -0.82,p<0.01; inversed LBG). Life history strategies that are often taxon-specific may also play a role. For example, marine mammals are endotherms and well-insulated, and so they can resist thermal perturbations. Plankton are eurythermic ectotherms and large-scale dispersers. Oceanic sharks and tuna/billfish are less eurythermic and more thermophilic than plankton (S1 Fig and S4 Text) and can therefore be affected more negatively by climate change and more positively by temperature that may influence evolutionary rates [28]. All groups were monophyletic and their patterns of niche saturation were not attributable to their taxonomic level; with the exception of the Foraminifera (subphylum), most groups were orders (S4 Table).
The second PCA performed on annual SST above 2.75˚C confirmed that mammals were more representative of a eurythermic community and unaffected by climate change (Fig 8D-8F and S3 Table). Plankton and fish are both more thermophilic than mammals (S1 Fig and S4 Text) and are also negatively affected by climate change and positively by evolutionary rates as assessed from the MTE. Plankton was less influenced by maximum species' niches at saturation between 2.75˚C and 30˚C, but correlations remained unchanged with mammals. It therefore appears that changes in net diversification rates influence patterns in niche saturation between tropical and temperate regions where the effect of maximum species' niche is weaker in contrast to the poles (Fig 3). Our results therefore agree with the time hypothesis [50], which stipulates that the tropics assemble more species over a longer time period because they are more stable; this is interesting because studies suggest that the tropics are both a species' cradle (higher origination rates) and a museum (more long-term climatic stability) [46].
Finally, we used the taxonomic groups as representative of ecological guilds to assess their niche saturation in terms of absolute value rather than relative patterns; foraminifers for protozooplankton, euphausiids for metazooplankton, both tuna/billfish and oceanic sharks for pelagic fish, and both cetaceans and pinnipeds for mammals (Fig 9). Niche saturation of mammals was low (~0.5%; median = 0.67%; Table 1) in comparison to other groups (Fig 9); fish saturation exhibited a mode at~10% (median = 7.62%), protozooplankton at~15% (median = 14.47%) and metazooplankton at~35-40% (median = 37.69%). A correction was applied to account for species that have not yet been discovered using two estimates: the first was mainly based upon the work of Appeltans and colleagues [36] and was applied to all guilds. The second, only used for plankton, was based upon the work of TARA [41]. While the first correction did not affect our results substantially, the second increased values of saturation for protozooplankton (Fig 9 and Table 1). After correction, our results showed an increase in saturation with diminishing organism complexity; the highest values were found for protozooplankton (or only plankton when Appeltans' correction was applied) and the smallest was found for mammals. Our results suggest that communities can be far from saturation [15,51] and that the total number of species on the chessboard diminishes with organismal complexity. The decrease with organism complexity can be explained by basic ecological and evolutionary processes. Endosomatic energy decreases from primary producers to higher trophic levels as a consequence of the second law of thermodynamics, diminishing the number of individuals and therefore species richness and niche saturation from producers to top predators [2,6]; a positive relationships between number of individuals and species richness has often been proposed in theories or hypotheses developed to explain the cause of the LBG, for examples, the productivity theory [52], the Area Hypothesis [53], and the Unified Neutral Theory of Biodiversity and Biogeography [54]. In addition to diminishing the number of individuals [28], larger body size also increases generation time, which slows down evolution [28,29]. We propose that the likelihood a taxon exhibits an ecogeographic biodiversity pattern different to the one imposed by the life chessboard is higher when its mean niche saturation is weaker; this proposition is valid in this work because correlations between expected and observed patterns are smaller (cetaceans) or even opposite (pinnipeds) for marine mammals.
Hutchinson [20] recognized that it is impossible to use all niche dimensions and it is therefore important to select those that control a large part of the spatial distribution of species. Temperature is a parameter frequently identified as being instrumental in explaining species, community, ecosystem and biome distribution in global-scale studies, influencing the biology and ecology of all species [4,[55][56][57]. We therefore used a one-dimensional (thermal) niche to reconstruct biodiversity [7,8] at a macro-scale. While we consider our approach to be valid at a macro-scale, we also anticipate that future versions of our framework will incorporate more ecological dimensions to provide more precise estimations of biodiversity and niche saturation at regional or local scales. At smaller scales, environmental parameters such as bathymetry, nutrients (e.g. nitrate, phosphate, silicate, iron), and oxygen concentration are likely to play an important role in shaping the arrangement of some ecological guilds or taxonomic groups [57][58][59]. Resource availability is also expected to be important as it controls the number of individuals, which in turn should promote more species [54,60,61]. However at a large scale, primary production (the unit of endosomatic energy produced per unit of time and space), is not correlated positively with marine biodiversity [6,18], suggesting that resource may not be the primary cause of the LBG.

Conclusions
Many hypotheses have been proposed to explain the primary cause of the latitudinal patterns in biodiversity [2,3,62]. While some authors have proposed that the LBG is due to the larger area of the tropical belts [61], the neutral model of biodiversity and biogeography [54], the middomain effect [63], or evolutionary explanations [64] are particularly popular. Perhaps the most compelling hypotheses have been those that invoke an environmental control of biodiversity such as environmental stability or energy availability [4,65,66]. As we have indicated however, temperature has been often suggested to explain large-scale patterns in the distribution of marine organisms [4,17], although the mechanisms by which this parameter acts to create latitudinal gradients in marine biodiversity has remained elusive [4]. The joint examination of biodiversity and niche saturation on a global scale in this study has provided interesting clues to the mechanisms underlying this pattern. Our results show that oceanic biodiversity is spatially constrained by an underlying chessboard-like structure that exerts a strong deterministic geographic effect on the macro-scale organization of marine eukaryotes by fixing the maximum number of species that can establish regionally; we exclude for the moment prokaryotes for which our framework has not been tested. We use the analogy of a chessboard to characterize this mathematical effect (Fig 4). Each square on the chessboard represents a geographical cell and the potential species that may occupy a square at saturation determine the sub-squares, which may be filled by speciation and immigration or purged by extinction (Fig 4). The deterministic effect imposed by the chessboard is strongest currently between the poles and other oceanic regions and weakest between temperate and permanently stratified biomes because the latter have smaller differences in biodiversity (Fig 3).
We propose that the presence of an underlying chessboard creates the establishment of the LBG for two reasons. First, it constraints large-scale oceanic patterns of biodiversity by determining a maximum number of species that can establish by speciation or immigration in a given region. The constraints the chessboard imposes on biodiversity may be rapidly perceptible because clade diversification takes place relatively quickly initially on a geological time scale [67]. Plankton and oceanic fish, having virtually unlimited dispersal, may rapidly conform to the mathematical structure imposed by the chessboard. Second, many clades should exhibit a LBG because their probability of emergence should be higher in the tropics where there are more available niches and palaeontological data have provided compelling evidence of greater rates of origination for tropical clades [68]. The higher degree of saturation found in the tropics for plankton and fish suggest that higher diversification rates or slower species turnover may originate from higher speciation rates, lower extinction rates or both [46,48,64].
Our study also suggests that the diversification processes (e.g. place of origination and time of emergence) may significantly blur the patterns of biodiversity imposed by the chessboard of life. Pinnipeds for example, that originated in extratropical regions of the North Pacific 25-27 Ma ago [49], do not exhibit an expected LGB. Nevertheless, pinniped biodiversity remains probably constrained in term of species number. Indeed, cetaceans, which are older but also originated in warmer regions [69] where the number of niches is higher, have a greater biodiversity than pinnipeds. Providing that pinnipeds have not evolved life history traits or anatomical structures that prevent them to be eurygraph, the diversification of this clade may continue and eventually lead to a more global-scale distribution that could be constrained by the chessboard of life. Because saturation is more important when biocomplexity diminishes, it follows necessarily, that the deterministic effect induced by the chessboard of life is likely to be more prominent for simpler (e.g. plankton) than for more complex taxonomic groups (e.g. mammals) of eukaryotes.
The rate of net diversification is a fundamental parameter because it determines the degree of niche occupancy in a square and the emerging life history traits and strategies make the chessboard specific to each clade; this explains the lack of universality in ecogeographic patterns (i.e. lack of congruence of the biodiversity among taxonomic groups). In a previous study, we have shown that climate-induced biogeographical shifts may alter chessboard's subsquares and therefore the species associated with them [9,66]. When a sub-square is occupied by a species they are both expected to move equatorwards in the case of global cooling, reinforcing the LBG slope, and to move polewards in the case of global warming, diminishing the LBG slope (Fig 4). This movement makes the chessboard highly dynamic on time scales from decades to millions of years. By merging evolutionary, ecological and mathematical mechanisms, our theory lies at the intersection of ecological and historical biogeography. Our theory, which we consider should also apply to the benthic domain and the terrestrial realm, explains why many taxonomic groups in the oceans exhibit a LBG, and it provides a comprehensive explanation why some clades may have different ecogeographic patterns.  Table. Summary of the results of the 2-order polynomial regression between the degree of niche saturation of each taxonomic group and both mean annual SST (between -1˚C and 30˚C) and annual SST variability. See Fig 6. All probabilities were significant, i.e. lower than 0.01. For all regressions, degrees of freedom were 179. Adjusted R 2 consider the reduction of the degree of freedom when a variable is added to the regression. RMSE: Root Mean Squared Error; the smaller the value, the better the fit of the regression. (DOCX) S2 Table. Summary of the results of the 2-order polynomial regression between the value of each explanative variable and both mean annual SST (between -1˚C and 30˚C) and annual SST variability. See Fig 7. All probabilities were significant, i.e. lower than 0.01. RMSE: Root Mean Squared Error. Higher RMSE are observed for LGM/today and mid-Pliocene/today changes because RMSE has the same unit as the response variable, which has an interval of variation higher than other variables. (DOCX) S3 Table. Table showing the correlations between changes in niche saturation of the six taxonomic groups and the explanatory variables according to mean SST and SST variability. Correlations in the lower triangle are based on mean SST ranging from -1˚C to 30˚C (185 degrees of freedom for all variables but LGM (157)) and those in the upper triangle are from 2.75˚C to 30˚C (149 degrees of freedom for all variables). Bold values are significant at p<0.01. Foram: foraminifers, Eupha: euphausiids, Shark: oceanic sharks, Tuna: tuna/billfish, Ceta: Cetaceans, Pinni:pinnipeds. Supplementary variables, indicated by a black bullet, are degree of eurythermy of the pseudo-community for each unit of mean SST and coefficient of variation of SST (Eury; see Fig 3) and similarly for the degree of pseudo-community thermophily (Thermo), LGM/today quantitative biodiversity changes (LGM), mid-Pliocene/today quantitative biodiversity changes (Plio),theoretical number of niche/species at saturation (theoS) and mass-corrected rate of evolution (Evol).