Truncated bimodal latitudinal diversity gradient in early Paleozoic phytoplankton.

The latitudinal diversity gradient (LDG)-the decline in species richness from the equator to the poles-is classically considered as the most pervasive macroecological pattern on Earth, but the timing of its establishment, its ubiquity in the geological past, and explanatory mechanisms remain uncertain. By combining empirical and modeling approaches, we show that the first representatives of marine phytoplankton exhibited an LDG from the beginning of the Cambrian, when most major phyla appeared. However, this LDG showed a single peak of diversity centered on the Southern Hemisphere, in contrast to the equatorial peak classically observed for most modern taxa. We find that this LDG most likely corresponds to a truncated bimodal gradient, which probably results from an uneven sediment preservation, smaller sampling effort, and/or lower initial diversity in the Northern Hemisphere. Variation of the documented LDG through time resulted primarily from fluctuations in annual sea-surface temperature and long-term climate changes.

INTRODUCTION

The latitudinal diversity gradient (LDG), classically described as the increase in taxonomic diversity from the poles to the equator, is probably the most striking and pervasive biogeographical pattern on Earth (–). Despite being a ubiquitous phenomenon documented for a large variety of functional and taxonomic groups, its origin, frequency among fossil groups, and its explanatory mechanisms remain widely debated. Proposed explanations for this pattern can be grouped into six categories: geographic area, geometric constraints, energy availability, habitat heterogeneity, evolutionary mechanisms [() and references therein], and niche theory (). While some of these explanations can be tested through neontological approaches, the effect of long-term phenomena such as climate changes and biological evolutionary events are inextricable using only extant taxa (–). The fossil record thus offers a unique and complementary perspective for studying the LDG, as it can provide insights into its appearance, stability through time, and long-term control mechanisms (, ). Studies of the LDG on fossil data have shown that it extends back to the Paleozoic (, –), but—regardless of its shape—the strength of the gradient has varied through time (, , , ) and the pattern even disappeared during intervals of mass extinctions (, , ). For the early Paleozoic (i.e., Cambrian-Ordovician, ca. 541 to 444 Ma), the presence of an LDG has been documented for three marine groups—brachiopods (, ), bryozoans (), and chitinozoans ()—as well as for the whole marine diversity of the Paleobiology Database (). Their LDGs consist of unimodal gradients (i.e., one diversity peak) centered on low to midlatitudes of the Southern Hemisphere (i.e., between 0° and 45°S), whereas, for most extant organisms, the LDG is either unimodal with a peak of diversity centered at the equator or bimodal with two peaks of diversity centered on the midlatitudes of each hemisphere (, –). In the present study, we expand our perspective on the LDG to the main representatives of phytoplankton in the early Paleozoic, i.e., the acritarchs. Defined as organic walled microfossils of uncertain biological affinity, most of them are now considered to be cysts of marine, generally phytoplanktonic unicellular algae (), and some may have a close biological affinity with dinoflagellates, as suggested by biomarkers and morphological similarities (). So far, acritarchs have been mainly used as tools for biostratigraphical, paleogeographical, and paleoenvironmental reconstructions, but only few studies considered them for the main source of information on phytoplankton in the early Paleozoic and studied them as such through, e.g., macroevolutionary and macroecological approaches (, ). Yet, acritarchs have one of the most densely and continuously sampled fossil record of any group at the global scale over the early Paleozoic and therefore constitute a unique material to study deep-time macroecological dynamics during this crucial period of the history of life. Our approach was to (i) determine the potential timing of establishment of an LDG for early Paleozoic acritarchs and its evolution through time according to long-term climate changes, (ii) analyze the relationship between the latitudinal distribution of acritarch diversity and abiotic factors whose values were reconstructed using paleoclimate simulations, and (iii) compare the observed shape of acritarch LDGs to the patterns predicted by a modeling approach based on the same paleoclimatic and paleogeographic context to estimate the portion of the empirical LGDs that could be due to preservation/sampling bias and distinguish it from the biological signal.

We conducted our analyses over 11 time bins covering the Cambrian-Ordovician interval with a new database of more than 9000 occurrences of acritarch genera from localities with known paleogeographic coordinates (table S1). Latitudinal paleobiodiversity curves were reconstructed by means of “range-through” indices and rarefaction methods. The relationships between acritarch diversity and abiotic factors [i.e., continental shelf area (CSA), sea-surface temperature (SST), and sea-surface salinity (SSS)] were tested for three key stages, i.e., the Cambrian stage 4, the Darriwilian, and the Hirnantian, by means of (i) a generalized linear model (GLM) with a response variable (i.e., diversity) following a Poisson distribution and (ii) the nonparametric Kendall’s rank correlation coefficient applied between the LDG and each abiotic factor separately. These three stages are representative of both the major macroevolutionary events recorded for the early Paleozoic—the Cambrian Explosion (), the Great Ordovician Biodiversification Event (GOBE) (), and the onset of the Late Ordovician mass extinction (), respectively—and the great climate variation that characterize this interval, i.e., greenhouse for the Cambrian, modern-like temperatures for the Middle Ordovician, and a glaciation for the end Ordovician (–). Last, we modeled species diversity for these three intervals using a model based on the MacroEcological Theory on the Arrangement of Life [METAL; (, , )] forced with the SST values simulated using the climate model on available paleogeographical reconstructions ().

RESULTS AND DISCUSSION

We find a unimodal latitudinal gradient of acritarch genus diversity for each time bin (Fig. 1 and fig. S1). The shallow gradient observed at the beginning of the Cambrian (i.e., Terreneuvian, ca. 531 Ma) may be explained not only by a poor sampling effort or preservation bias (fig. S2) but also by the end-Ediacaran mass extinction of acritarchs (, ). Contrary to the classical unimodal LDG observed for many modern taxa (, ), the acritarch LDG does not exhibit a diversity maximum at the equator or near the tropics but rather over mid- to high latitudes of the Southern Hemisphere throughout the early Paleozoic. Starting at 35°S to 40°S during the Cambrian series 2 (ca. 515 Ma; Fig. 1 and fig. S1), the diversity peak then slowly drifts southward. It reaches 45°S during the Furongian (ca. 491.2 Ma) and its highest latitudinal position (i.e., 50°S to 60°S) during the Tremadocian-Dapingian interval (ca. 485.4 to 467.3 Ma), before splitting in two peaks corresponding to the drift of highly diverse areas (i.e., Baltica and the Armorican Terrane Assemblage; Fig. 2) during the Middle Ordovician (Fig. 1 and fig. S1) [see (, ) for links between biodiversity and plate tectonics]. Only the peak located at the lowest latitude and corresponding to Baltica persists until the end of the Ordovician (Fig. 1 and fig. S1). Similar LDGs are also recovered from subsampled diversity when using the method of Marcot et al. () (except for the Miaolingian; figs. S3 and S4), although their latitudinal extent is often reduced compared to the nonrarefied diversity curves because of the small sample size of some latitudinal bands (fig. S1).

Fig. 1

Interpolation map of acritarch genus diversity calculated per 5° latitudinal bands.

Diversity was calculated with the standard range-through approach excluding single-interval taxa [RTexS; see ()] for each series of the Cambrian and each stage of the Ordovician. Note the slight southward drift of peak diversity through the Cambrian, the colonization of intertropical latitudes during the Early Ordovician, and the subsequent duplication of the diversity peak. The null diversity observed above ~30°N is a consequence of the particular paleogeography of the early Paleozoic, with a Northern Hemisphere mainly constituted by the Panthalassa Ocean. Data are lacking for this large area, because the vast majority of its oceanic crust has long been recycled in subduction zones.

Fig. 2

Paleogeographical context of the early Paleozoic.

Paleogeographical reconstruction for the mid-early Paleozoic (i.e., Jiangshanian, 491.75 Ma) computed with GPlates 2.0.0 () using the rotation file of Scotese () and paleolocation of the 328 fossil localities used in this study (black dots) for this time interval. Striped parts correspond to CSAs after Scotese ().

In addition to the position of peak diversity, the steepness of the LDG also changes through time. It markedly increases during the Floian in the Southern Hemisphere (ca. 473.85 Ma), paralleling an increase in diversity in the equatorial and low latitudes of the Northern Hemisphere (Fig. 1 and fig. S1). A similar tropical diversification is known for marine organisms considered as a whole (). This transition from a poor Cambrian to a highly diverse Ordovician intertropical zone is probably related to a global climate change (, , ). While the Cambrian corresponded to a greenhouse climate with equatorial SST probably too high to sustain highly diverse communities (fig. S6) (), a steady cooling trend occurred through the Early Ordovician to reach the range of modern equatorial SST by the Middle Ordovician (fig. S6) (), allowing all major marine groups to colonize low latitudes (), before ending in a glaciation associated to a massive extinction during the Hirnantian (). Temporal variations in the position of the LDG diversity peak of early Paleozoic acritarchs and the general shape of this LDG thus appear to have been related to the climatic regimes of the Cambrian-Ordovician interval and, to a lesser extent, to the position of continental masses. Brayard et al. () have shown that the position of peak diversity and steepness of LDGs in marine environments depend on the shape and magnitude of the latitudinal SST gradient: Sigmoid-like SST gradients generate a bimodal LDG with peaks of diversity located at midlatitudes, while a reduced steepness of the SST gradient leads to a weakened unimodal LDG with a plateau of equal diversity between the tropics. This implies that bimodal LDGs with diversity peaks at intermediate latitudes are a characteristic feature of icehouse periods, whereas unimodal LDGs with an equatorial peak or a plateau of equal diversity are characteristic of greenhouse periods (, ). Alternatively, supporting the observations of Mannion et al. (), Beaugrand et al. () have shown that, in case of a global cooling, the LDG steepness increases because of the movement of species toward the equator, leading to a unimodal gradient centered at the equator, whereas, in case of global warming, it is attenuated because of the poleward shift of species, leading to a bimodal gradient with peaks of diversity at intermediate latitudes (, ). The LDG of early Paleozoic acritarchs is steeper, and its peak of diversity is closer to the equator during the colder Middle-Late Ordovician than during the Cambrian, and no equatorial peak or plateau of equal diversity is observed under the Cambrian greenhouse climatic conditions (fig. S1). This corresponds more to the model of Beaugrand et al. () and its deviation from the model of Brayard et al. () might be explained by a threshold temperature effect: Peak diversity could not occur at the equator during the Cambrian, because equatorial temperatures were too high to allow most organisms to survive, and this peak slowly shifted toward more equatorial latitudes during the Ordovician cooling, with a rapid first step of diversification at equatorial latitudes during the Early Ordovician (i.e., Floian; fig. S1). Hence, it seems that—as for the majority of marine groups of organisms—long-term climate trends had an important influence on the distribution of diversity of early Paleozoic phytoplankton.

The results of our analyses of correlation between SST and acritarch diversity support this hypothesis. For the Cambrian stage 4, seasonal variation in SST (SSTdiff; i.e., the difference between summer and winter SST) appears as a positive predictor of acritarch diversity in the three best GLMs [163 < AICc (corrected Akaike information criterion) < 164, coefficient of determination (R2) = 0.66; Table 1]. For the Darriwilian, which corresponds to the major diversity peak of the GOBE and during which global climate was cooler than during the Cambrian (fig. S6) (), SSTdiff is also significantly correlated to diversity as revealed by Kendall’s test (τ = 0.61, P < 0.001; table S2), and mean annual SST (SSTan) and squared SSTan (SST2) appear alternatively as negative predictors of diversity in the two best GLMs (125 < AICc < 126, R2 = 0.96; Table 1). Such negative relationships between temperature and diversity may arise when temperatures are beyond thermal niches of most taxa under study (, , , ) or when taxa are more diverse at cooler temperatures (). Eventually, during the colder Hirnantian stage (fig. S6), the best GLM shows that SSTan and SSTdiff are both predictors of acritarch diversity (121 < AICc < 123, 0.91 < R2 < 0.92; Table 1), a result supported by Kendall’s test (−0.45 < τ < −0.40, P < 0.05; table S2). Together, these results clearly show the existence of a relationship between SST, especially its seasonal variations, and the spatial distribution of acritarch diversity throughout the early Paleozoic. Temperature has often been cited as one of the main explanatory factors for large-scale distribution patterns of organisms (, , , –). Many studies, in both marine and terrestrial environments, have shown that present-day temperature, along with some other environmental factors such as precipitation on land and productivity in the sea, is a powerful first-order predictor of LDG patterns [() and references therein]. In the marine realm, Tittensor et al. () have shown that SST was the only environmental predictor highly correlated with diversity across 13 taxonomic groups, a finding confirmed for phytoplankton in particular (). This positive relationship has been observed for a large variety of environments, scales, and taxa [e.g., (, , , )]. Our results show that acritarchs are no exception to this rule and that, as far back as the early Paleozoic, the spatial distribution of phytoplankton diversity was controlled, at least in part, by this abiotic variable.

Table 1

GLM results for acritarch genus diversity per 5° latitudinal band for the Cambrian stage 4, the Darriwilian, and the Hirnantian.

Only the two or three best models (i.e., with the lowest AICc value and highest Akaike weight) from all possible combination of the variables considered are shown here for each time interval. Empty cells indicate a variable not included in the corresponding model (i.e., a variable that does not increase the predictability of taxonomic richness). Neg. bin., negative binomial; CSA, continental shelf area (106 km2); SSTan, mean annual SST; SST2, squared SST; SSTdiff, difference between summer and winter SST; SSSan, mean annual SSS; SSS2, squared SSS; SSSdiff, difference between summer and winter SSS; and ∅, not applicable (AICc not being computable for negative binomial GLMs). Dispersion = residual deviance/degrees of freedom, overdispersion when >1.

Temporal interval	Family	Dispersion	R2	CSA	SSTan	SST2	SSTdiff	SSSan	SSS2	SSSdiff	AICc	Akaike weight
Cambrian stage 4	Poisson	1.477	0.66				0.056		0.003		163	0.114
	Poisson	1.489	0.66				0.055	0.236			163	0.101
	Neg. bin.	1.092	0.66				0.053		0.003		164	∅
Darriwilian	Poisson	0.301	0.96		−0.071			16.815	−0.242		125	0.31
	Poisson	0.322	0.96			−0.001		12.758	−0.184		126	0.262
Hirnantian	Poisson	2.276	0.91	0.172	−0.244		0.271		0.017		123	0.303
	Poisson	2.279	0.91	0.175	−0.244		0.268	1.139			123	0.297
	Neg. bin.	1.202	0.92	0.199	−0.264		0.286		0.018		121	∅

Salinity is known to affect the distribution patterns of marine organisms at local and sometimes at regional scales, but its influence on global distribution patterns of diversity seems to be weak (, ). Regarding early Paleozoic acritarchs, GLMs reveal a positive correlation between diversity and either mean annual SSS (SSSan) or squared SSSan (SSS2) during the Cambrian stage 4 (163 < AICc <164, R2 = 0.66; Table 1; supported by Kendall’s test, τ = 0.52, P < 0.001; table S2). For the Darriwilian, SSSan and SSS2 appear as positive and negative predictors of diversity, respectively, using a GLM (125 < AICc < 126, R2 = 0.96; Table 1), while SSSdiff is negatively correlated to diversity using Kendall’s test, τ = −0.33, P < 0.05; table S2). For the Hirnantian, GLM shows that SSSan and SSS2 both positively predict the latitudinal distribution of acritarch diversity (121 < AICc < 123, 0.91 < R2 < 0.92; Table 1; not supported by Kendall’s test; table S2). A relationship thus seems to exist between the early Paleozoic acritarch diversity and salinity, although not very stable according to the GLM results (i.e., the same salinity variable is not systematically found as a predictor of diversity in the best models; Table 1). In addition, given that SSS values remain quite stable in the intertropical zone from the Cambrian stage 4 to the Hirnantian (fig. S6D) while a major shift of diversity is recorded in this area during the Early Ordovician, we argue that temperature was a more limiting factor of acritarch diversity, because the colonization of the intertropical zone by acritarchs during the Floian is not associated with a change in SSS but much more likely with a decline in SST ().

The third abiotic factor investigated here, i.e., CSA, has often been found to be related to marine diversity (–), but it should not be considered a priori as a bias. Indeed, during transgression events, flooded continental areas increase, offering more niche space and opportunities for diversification and affecting the production of sediments, which provides a general explanation for the covariation between fossil marine biodiversity and the amount of marine sedimentary rocks (). However, this signal varies among groups of organisms, spatial scales, and time intervals (). Regarding acritarchs, CSA appears as a positive predictor of diversity for the Hirnantian when analyzed through GLM (121 < AICc < 123; 0.91 < R2 < 0.92; Table 1) and for the Cambrian stage 4 when using Kendall’s test (τ = 0.454, P < 0.001; table S2). We interpret this result as evidence that, although present, the sampling/preservation bias (fig. S2) is probably sufficiently evenly distributed in latitude so as not to have an impact on the relationship between the diversity of acritarchs and the surface area of their living environment.

Assuming that the diversity of early Paleozoic acritarchs is concentrated on continental shelves, we compared their empirical LDGs with the LDGs obtained by modeling species distribution on continental shelves only using the METAL model (, , ) to estimate the magnitude of the preservation/sampling bias and separate its effects from the biological signal. Validated for modern plankton, fish, and cetaceans (), this model shows that the primary cause of LDGs in the oceans lies in the interaction between species thermal tolerance and both seasonal and year-to-year fluctuations in temperature, together with the generation of a mid-domain effect. For the Cambrian stage 4 and the Darriwilian, modeled LDGs are bimodal with two peaks of diversity located respectively at 37.5°S and 42.5°N, and 32.5°S and 32.5°N, when averaged by 5° latitudinal bands (Fig. 3). The Northern peak of diversity is smaller than the Southern one for both stages (Fig. 3). For the Hirnantian, the gradient is steeper and unimodal with a plateau of maximum diversity ranging from ~22.5°S to 27.5°N (Fig. 3). The three modeled LDGs end around 30°N for each stage, because no continental shelf existed above this latitude during the early Paleozoic, the Northern Hemisphere being essentially covered by the Panthalassa Ocean. Modeled and empirical LDGs are significantly correlated for the Darriwilian (Kendall’s τ = 0.40, P < 0.05) but not for the Cambrian stage 4 and the Hirnantian (τ = 0.05, P = 0.72 and τ = 0.04, P = 0.83, respectively). Although the Southern Hemisphere peaks of diversity of the modeled and empirical LDGs are located at close latitudes during the first two stages, some important differences can be observed: Modeled LDGs are bimodal during the first two stages while empirical LDGs are unimodal, and the Hirnantian modeled LDG shows a plateau of equal diversity centered at the equator, while the empirical LDG shows a diversity peak located around 27.5°S. This change from a bimodal to a unimodal modeled LDG between the Darriwilian and the Hirnantian does not match our observations, but the fact that the peak of diversity moves toward lower latitudes during this interval in both cases (i.e., modeled and empirical; Fig. 3) shows the influence of the climate cooling on the latitudinal distribution of diversity, allowing taxa to colonize lower latitudes. The greater divergence in diversity patterns between modeled and empirical data for the Hirnantian (Fig. 3) may also be related to the rapid and intense climate changes that characterize this time interval. The Hirnantian is marked by cold temperatures, a major glaciation, and a severe drop in sea level that drained epicontinental seaways, producing harsh climate in low and midlatitudes (). Considered responsible for the mass extinction that occurred during this time interval (), these abrupt changes probably make the spatial dynamics of biodiversity less easily modelable than that of other early Paleozoic stages.

Fig. 3

Comparison of latitudinal gradient of acritarch generic diversity and modeled species diversity.

(A and B) Cambrian stage 4 (510 Ma, 32 PAL, i.e., preindustrial atmospheric level of CO2, 1 PAL = 280 parts per million), (C and D) Darriwilian (460 Ma, 12 PAL), and (E and F) Hirnantian (445 Ma, 5 PAL). Modeled spatial distribution of species diversity corresponds to mean results of 1000 METAL simulations for continental shelves (>200 m, bright colors) and ocean surface (<200 m, light colors), based on Scotese and Wright’s () maps and SST simulated using the FOAM model, projected on a 5° × 5° spatial grid (A, C, and E) and calculated by 5° latitudinal bands [gray dots on (B), (D), and (F) plots]. Acritarch generic diversity calculated by 5° latitudinal bands [white squares on (B), (D), and (F)] using the range-through counting approach (, ). AGD, acritarch generic diversity; MSD, modeled species diversity.

Bimodal gradients are very common for extant marine organisms (–). They have been documented for marine taxa as varied as bivalves (, ), brachiopods (), bryozoans (), foraminifers (, , ), seaweeds (), crustaceans (), and copepods (, ). The LDG of modern phytoplankton is not homogeneous and varies among groups. The diversity of dinoflagellates—probably the closest relatives of acritarchs—peaks at both mid-southern (~30°S) and mid-northern latitudes (~20° and 45°N) [(), appendix S2], and coccolithophores show a bimodal gradient with highest diversity at ~35°N and 35°S [(), app. S2]. The only other phytoplankton groups that were subject to such analyses, i.e., calcareous nannoplankton and diatoms, show respectively a classical unimodal gradient centered on the equator () and various patterns depending on the data and method chosen (–). Although very classic in the modern oceans, bimodal LDGs have never been observed for Paleozoic organisms (, , ). This marked difference may result either from the fact that modern LDGs are controlled by different factors than early Paleozoic LDGs or that the early Paleozoic fossil record is biased and the LDGs one reconstructs are distorted images of the original ones. The fact that LDGs modeled here are bimodal and that empirical unimodal LDGs centered on the midlatitudes of the Southern Hemisphere are characteristic of Paleozoic organisms support the second hypothesis. This seems all the more likely for acritarchs, because the LDGs of extant pelagic organisms () and dinoflagellates in particular (), as well as modeled LDGs of virtual planktonic species (, , ), are also bimodal. Therefore, we argue that rather than being an atypical unimodal gradient centered on the midlatitudes of the Southern Hemisphere, the early Paleozoic acritarch LDG probably corresponds to a truncated bimodal gradient whose Northern peak of diversity cannot be recovered or never existed. Either it cannot be recovered because of a reduced sediment and/or specimen preservation, and/or a smaller sampling effort in the Northern Hemisphere; or it never existed, because some habitat-related factors essential for survival of acritarchs but not accounted for in the model were lacking in most of the Northern Hemisphere because of the small extent of landmasses.

This first quantitative analysis of the latitudinal distribution of early Paleozoic acritarch diversity shows that (i) the earliest representatives of marine phytoplankton were characterized by an LDG already since the Cambrian Explosion; (ii) this LDG most likely corresponded to a truncated bimodal rather than a unimodal gradient; (iii) variation in the shape of this LDG through the early Paleozoic can be attributed to both long-term climate cooling and plate tectonics; and (iv) at the shorter temporal scale of the stage, the latitudinal distribution of acritarch diversity was influenced by SST and its seasonal fluctuations and also, to a lesser extent, by salinity and the available area of continental shelf. These results strongly contrast with the unimodal nature of other LDGs found for Paleozoic organisms, such as brachiopods and bryozoans (, ). It seems that no single model of LDG prevails, but rather that several models, which vary in amplitude and shape in space and time, as well as according to the ecological characteristics and life traits of organisms [such as size; (–)], have coexisted at least since the early Paleozoic.

MATERIALS AND METHODS

Fossil occurrence data

The dataset that we used, initially created as part of the PhytoPal project by Mullins et al. (, ), consists of 1503 species of acritarchs belonging to 287 genera and distributed among 328 fossil localities, at the global scale (Fig. 2 and table S1). It is a compilation of all occurrence data of Paleozoic organic walled microphytoplankton published before 2007: It synthesizes the presence of a species in a locality during a time interval (i.e., chronostratigraphic unit) as reported by a reference (table S1). The data were filtered and cleaned by removing or reassigning illegitimate, questionable, and synonymous taxa and converting local to global chronostratigraphic units. Temporal resolution is at the series level for the Cambrian and at the stage level for the Ordovician. Geographic coordinates of fossil localities were recovered from the literature or using fossil locality names, and their paleocoordinates were calculated by rotating their present-day coordinates with the software GPlates () version 2.0.0, using the rotation file supplied by Scotese () for each Cambrian series and Ordovician stage (table S3). Fossil localities with only one taxon per time interval were discarded from the analyses. A curve illustrating the variations through time of the global diversity of acritarch genera is available in the Supplementary Materials (fig. S7).

Latitudinal diversity gradient

We tested for the existence of an LDG by calculating the paleobiodiversity of acritarchs by 5° latitudinal bands at the genus level, because it is considered more robust than the species level for such analyses (see Supplementary Text). We estimated the taxonomic richness by means of two range-through indices [see review of Alroy ()] including or excluding single-interval taxa: RTinS and RTexS, respectively, but we focused on the RTexS values for the diversity (Figs. 1 and 3, fig. S1, and table S4) and correlation analyses (Table 1 and table S2). Range-through indices account for every taxon known to occur in a latitudinal band (i.e., real occurrence) as well as every taxon inferred to be present (i.e., recorded south and north of the considered latitudinal band). This range-through approach assumes that studied organisms tend to fill the latitudinal extent between their range end points, and it ignores minor within-range absences that could be the result of geographically and temporally biased sampling (), a common practice for LDG studies in the past (, –). All time intervals considered, the average extent of this filled latitudinal range between two real occurrences corresponds to two latitudinal bands (i.e., 10°). Because of uneven sampling in the fossil record, we additionally calculated rarefied genus richness through shareholder quorum subsampling (, ). However, as this standard method appears not well adapted to our dataset (see Supplementary Text), we also used the less stringent rarefaction method proposed by Marcot et al. () (see Supplementary Text), which offers the advantage of not distorting the pattern of genus diversity. The number of occurrences randomly drawn (without replacement) in each latitudinal band corresponds to the maximum genus diversity observed within any single band of the time interval to maintain the possibility of recovering the face-value richness, which would not be possible with any smaller quota. We conducted the diversity analyses at the series level for the Cambrian and at the stage level for the Ordovician.

Correlation with abiotic factors

We tested the relationships between acritarch diversity and abiotic factors (CSA, SST, and SSS) for three key stages, i.e., the Cambrian stage 4, the Darriwilian, and the Hirnantian. These three time intervals are among the best sampled (see fig. S2); therefore, more likely to produce reliable diversity estimates. We simulated paleoenvironmental conditions using the Fast Ocean Atmosphere Model (FOAM) version 1.5 () (see Supplementary Text) with the same continental reconstructions than that used for the analyses of acritarch LDG (). CSA was derived from the paleogeographical reconstructions, and values of SST and SSS were extracted from the climatic simulations (see Supplementary Text). We tested the relationships between these abiotic factors (i.e., mean values per latitudinal band for SST and SSS) and the latitudinal distribution of acritarch diversity (RTexS; i.e., range-through diversity without single-interval taxa) using two approaches (see Supplementary Text for details): (i) A GLM (see Supplementary Text) with a response variable (i.e., diversity) following a Poisson distribution by means of a log link function, applied to all abiotic factors together, using the R environment () and the packages “MuMIn” () and “MASS” (), and (ii) the nonparametric Kendall’s rank correlation coefficient applied between the LDG and each abiotic factor separately.

Species distribution model

We modeled theoretical species diversity for the Cambrian stage 4, the Darriwilian, and the Hirnantian using a model based on the METAL [see Supplementary Text and (, , )] forced with the SST values simulated using the climate model on the paleogeographical reconstructions of Scotese and Wright (). METAL uses a number of models that are specifically designed to be used at the species or community organizational level (, , , ). Here, the model implements a set of basic ecological/climatic principles to consider the interaction between the thermal niches of species and spatiotemporal fluctuations in temperature at an annual time scale. The model has been fully described and tested in Beaugrand et al. (, , ). The principle of the model is simple: It starts by creating a large number of niches on the basis of temperature (here, all possible niches between −1.8° and 44°C) (). We therefore created pseudospecies having each a unique thermal niche with distinct degrees of eurythermy and thermophily (, ). We allowed pseudospecies to colonize any given region of the global ocean, provided that they could withstand the local annual thermal regime. In a given area, each pseudospecies has a unique niche after the principle of competitive exclusion of Gause () while considering niche overlapping (). By reconstructing pseudocommunities, the model reproduces the spatial arrangement of biodiversity. METAL models have been tested at annual to millennial time scales for different marine plankton groups and fishes for which changes in community composition and biodiversity are the result of the (realized) niche-environment interaction. Correlations between observed community reorganization/biodiversity at these time scales and changes assessed from METAL were highly significant (, ). To estimate pseudospecies diversity, we randomly selected 10,000 among 101,397 possible pseudospecies and repeated the procedure 1000 times. Because our main goal was to examine large-scale diversity patterns, we increased the number of pseudospecies to improve our estimations of diversity. However, large-scale diversity patterns using the maximum number of genera of the acritarch dataset (i.e., 188, 654, and 410 for the Cambrian stage 4, the Darriwilian, and the Hirnantian, respectively) and 10,000 gave almost identical patterns.

We integrated the number of pseudospecies for each 5° latitudinal bands and compared this modeled species diversity to the empirical genus acritarch diversity that we observed. Our model cannot realistically implement biotic interactions in the construction of pseudocommunities on a global scale, although multiple examples suggest that they can be quite important in some ecosystems (, ). There are currently insufficient data to include these potentially important top-down and competitive interactions over the scales needed. Therefore, our overall approach was to consider the fundamental niche-temperature interaction to reconstruct global diversity patterns, leaving scope for future refinement. We are aware, however, that, at a more local scale, biotic interactions may precipitate or alleviate the climatic effects on pseudocommunities ().

We assume here that genus diversity is directly proportional to species diversity. Some studies have shown that this assumption holds at least for contemporary large-scale diversity patterns of some plankton protists [e.g., diatoms; (, )].