Coupling of geographic range and provincialism in Cambrian marine invertebrates

The Cambrian witnessed a noteworthy radiation in biodiversity (Sepkoski 1997; Maloof et al. 2010; Na and Kiessling 2015) and a remarkable divergence in Bauplans (Erwin et al. 2011) of marine metazoans. Although the temporal and spatial patterns of Cambrian biodiversity have been widely documented (e.g., Zhuravlev and Riding 2001; Jago et al. 2006; Li et al. 2007; Peng et al. 2012; Na and Kiessling 2015; Kröger et al. 2019; Rasmussen et al. 2019; Fan et al. 2020), less attention has been paid to the fundamental biogeographic structure and its changes through time. Biogeographic patterns may provide insights into evolutionary processes during large-scale diversifications at different temporal scales (Valentine et al. 1978; Fortey et al. 1996; Crame and Owen 2002; Huang et al. 2018).

Cambrian provinciality has been recorded among various fossil groups at different taxonomic levels (Theokritoff 1979; Lu 1981; Zhuravlev 1986; Steiner et al. 2007; Yang et al. 2015). A two-part or three-part provincial scheme has been suggested to represent early Cambrian trilobite biogeography (Lu et al. 1974; Palmer and Repina 1993; Fortey et al. 1996): the Redlichiid Province in Gondwana; the Olenellid Province in Laurentia, Baltica, and Siberia; and sometimes an intermediate province representing overlap domains (Lu 1981; Pillola 1989). The identification of trilobite biogeography relied largely on distinctive endemic trilobites among major continents and has proved unreliable in light of phylogenetic analyses (Lieberman 1998; Meert and Lieberman 2004). Taking phylogenetic history of trilobites into account, more complex biogeographic schemes among trilobites were advocated for the early and the late Cambrian (Meert and Lieberman 2004; Álvaro et al. 2013), and many attempts have been made to compare phylogenetic units among tectonic blocks to reconstruct plate tectonic configurations (Jell 1974; Lieberman 1997, 2005; Álvaro et al. 2013; Steiner et al. 2021). However, traditional approaches based on taxon- and time slice–specific analysis (e.g., Theokritoff 1979; Yang et al. 2015) can hardly provide a panoramic perspective on secular trends in biogeography and their link to diversification patterns (Kocsis et al. 2018a,b). Moreover, the heterogenous fossil coverage distorts our understanding of deep-time biogeographic patterns (Alroy et al. 2008; Close et al. 2020), which has necessitated multi-taxon approaches and a more reproducible quantitative framework in delineating biogeographic units (Kocsis et al. 2021).

Here, we derive the first-order biogeographic patterns of macroinvertebrates across the entire Cambrian period. We apply a network analysis of species occurrence data while accounting for uneven sampling. The network-based approach has been widely applied to fossils to reconstruct biogeographic patterns in deep time (Sidor et al. 2013; Vilhena et al. 2013; Rojas et al. 2017; Dunne et al. 2018; Huang et al. 2018; Rong et al. 2020), but no study has focused on the Cambrian period. Our methodology traces community structures and outlines time-traceable biogeographic units (bioregions; Kiel 2017; Kocsis et al. 2018a,b). We employ a multi-fauna approach, but for comparison with earlier works, we also apply our method to only trilobite occurrences. By assessing the relationship between biogeographic changes and species’ geographic ranges, we further explore the timing of provincial divergence/convergence and derive the potential drivers of biogeographic changes during the Cambrian radiation.

We downloaded occurrence data of Cambrian marine invertebrate species from the Paleobiology Database (https://paleobiodb.org) on September 3, 2021. The raw data comprised 9214 collections and 41,912 fossil occurrences. For temporal binning, we updated the framework for Cambrian biostratigraphic correlation following Na and Kiessling (2015) (Supplementary File F1), with time binning created with the divDyn package (Kocsis et al. 2019) in R (R Core Team 2021). Ages refer to the Geologic Time Scale 2020 (Gradstein et al. 2020). Only collections that could be assigned to a single stage were retained. We also removed all occurrences that could not be assigned to a species with confidence as well as form species and ichnospecies. Species occurrences with the qualifiers “cf.”, “?”, “aff.”, “n. sp.”, “n. gen.”, or “indet.” were disregarded. The final dataset comprised 4397 collections and 15,622 species occurrences of 2644 species (Supplementary Table S1; also see Supplementary Material available through Zenodo). Nearly half of the original collections were dropped due to the uncertainty in stratigraphic binning.

Paleo-coordinates were obtained by rotating the present-day coordinates of fossil collections with the GPlates v. 2.2.0 software (Müller et al. 2018). We used the rotation file supplied by C. R. Scotese in the PALEOMAP PaleoAtlas for GPlates package (Scotese 2016). Accurate Cambrian reconstructions of the configuration of Earth's continents come with substantial spatial and temporal uncertainties (Landing et al. 2013). With a different paleogeographic model that reconstructs Cambrian geography differently, our results might be somewhat different. We rotated all collections to their positions at 510 Ma (the middle of the Cambrian) to (1) make the variation among cells comparable through time without largely influenced by tectonic variations, (2) better determine how provincial assignments evolved in a particular region, and (3) better relate the source data to our paleobiogeographic reconstructions.

We aggregated the collections based on their paleo-coordinates to equal-area hexagonal grid cells using a tessellated icosahedral geographic grid from the icosa package in R (Kocsis 2017). To decrease the mismatch between the areas covered by a geographic cell in two subsequent time slices due to plate tectonic movements while retaining within-continent variation as much as we could, we used a rather coarse spatial resolution that gives roughly 252 cells per time interval with one grid cell being roughly 2 million km² in area (edge length = 8°). In this way, 2520 spatiotemporal cells (10 geological stages × 252 geographic cells) were outlined and used in the network analysis to construct a time-uniform biogeographic partitioning scheme. Following a previous analysis (Kocsis et al. 2018b), we employed a minimum quota of at least 10 occurrences in each spatiotemporal unit to be included in the bipartite network.

The collections in a time interval were assigned to equal-area geographic grid cells, and species were tabulated in every spatiotemporal cell. From that matrix, a bipartite occurrence network (Vilhena and Antonelli 2015) was constructed based on shared species between the 2520 spatiotemporal units, wherein every spatiotemporal cell was represented with a vertex, and units were connected if the spatiotemporal cells shared species. The strength of connections was determined by the number of shared species and sampling. We performed a weighted projection from the bipartite network onto the node subset (Rojas et al. 2017) and applied the Infomap community detection algorithm (Rosvall and Bergstrom 2008) to partition the spatiotemporal grid cells within the projected network (Edler et al. 2016; Rojas et al. 2017). The defined modules (i.e., bioregions) were then projected into geographic space.

To investigate the biogeographic dynamics, we used two indices (Kocsis et al. 2018a): (1) by-cell turnover and (2) by-bioregion turnover. By-cell turnover represents total biogeographic turnover and is described as the changes of bioregion assignments of each grid cell between consecutive time bins. This metric ranges from 0 (no change) to 1 (complete turnover) and is insensitive to changes in sampling intensity (Kocsis et al. 2018a). By-bioregion turnover estimates the proportions of emergence and disintegration of bioregions to the total number of outlined bioregions in the time interval. In other words, this index traces the origination and extinction of bioregions through time. By-bioregion turnover series are based on subsampling (classical rarefaction) to 10 geographic cells in every age with 100 iterations.

Although the dataset includes a large number of Cambrian data (more than 40,000 fossil occurrences), it remains incomplete. When spatial sampling is sparse, as in our case, the biogeographic membership of a cell can be sensitive to the exact location of a cell boundary when there is considerable turnover over distances that one cell covers. To overcome the biasing effect of this phenomenon, we report time series based on the biogeographic partitioning as averages of 400 replicates with the grid rotated randomly in each iteration. To further determine whether biogeographic patterns change with spatial resolution, the analyses were repeated with different geographic resolutions (edge lengths of 10°, 6.7°, 5.7°, 5°, and 4.4°, respectively).

Following Kocsis et al. (2021), we employed Hurlbert's probability of interspecific encounter (PIE) to estimate provinciality. PIE measures the probability that two randomly sampled grid cells belong to different provinces. This provinciality metric ranges from values of 0 (perfectly uniform) to 1 (perfectly uneven). Although direct calculations of species’ ranges (e.g., Stigall and Lieberman 2006; Hendricks et al. 2008; Na and Kiessling 2015) can provide information about species’ maximum possible geographic distribution, they can be biased. To circumvent the problem of heterogenous spatial sampling, we employed proportional occupancy as a proxy of species’ geographic ranges. Geographic occupancy is here defined as the proportion of all sampled grid cells in which a particular species occurs (Kiessling and Kocsis 2016). We also calculated species’ stratigraphic durations by the number of stratigraphic intervals during which a species existed. All calculations were implemented in the R programming language and environment (R Core Team 2021).

The first three Cambrian stages, that is, the main phase of diversification during the Cambrian radiation (Na and Kiessling 2015), saw a prominent increase in the number of bioregions (Fig. 1). The number of bioregions peaked in Stage 4. The total biogeographic turnover increased substantially from Stage 2 until the Wuliuan (Fig. 2). Biogeographic turnover was highest in the Wuliuan and the Paibian stages. Patterns are robust to changes of spatial resolution and minimum occurrence quota (Supplementary Figs. S1, S2). The rate of emergence of bioregions decreased through the first four stages and fluctuated at a lower level afterward, while the rate of disintegration of bioregions increased slightly from Stage 3 to the Guzhangian (Fig. 2).

The extent of biogeographic changes varied geographically (Fig. 1). The increase in bioregions during the Cambrian radiation was prominent in tropical areas of Laurentia, Siberia, and Gondwana. Tropical areas were also dominated by single-interval species (singletons; Table 1). The distribution of species’ stratigraphic durations was right-skewed in both tropical and nontropical areas, but the mean species’ stratigraphic longevity was slightly shorter in tropical areas than in nontropical areas (Table 1). Similar patterns are also evident in trilobite biogeography (Table 1, Fig. 3). Overall bioregions spanned across continents throughout the Cambrian, but most trilobite bioregions were constrained to a single continent from Stage 3 to the Wuliuan (Fig. 3).

The disintegration and subsequent reorganization of bioregions since Stage 3 showed no distinct pattern, but most bioregions seemed to follow expansion–contraction trajectories with respect to their geographic distribution. For example, the fauna of bioregion 1 (red grids in Fig. 1) that persisted from the Fortunian to Stage 4 achieved worldwide distribution (from the tropics to the southern polar region) during Stages 2 and 3 until its disintegration during Stage 4. Only the Wuliuan, the Guzhangian, and the Jiangshanian featured some short-lived (single-interval) bioregions (Fig. 1).

The PIE provinciality metric indicates a profound increase in provinciality from Stage 2 to Stage 3 (Fig. 4). Although the increase of provinciality at the beginning of the Cambrian is concordant with the main diversification phase (Fig. 5), we found no correlation between changes of sampling-standardized species diversity and provinciality (rho = 0.18, p = 0.64), suggesting that the changes in provinciality were not related to changes in species diversity. On the contrary, we found a strong correlation between the provinciality and the number of singletons (rho = 0.78, p = 0.01) over the whole period studied, suggesting that the singletons exerted a potential role in driving Cambrian provincialism.

Species’ geographic occupancy decreased visibly from the Fortunian to Stage 3 (Fig. 6), despite the profound increase in the diversity of mobile species during the same period (Fig. 7). This contraction of geographic occupancy during the critical diversification phase is also evident in subsets of mobile/non-mobile and singleton/non-singleton species (Supplementary Fig. S3).

Our results confirm a profound increase in provinciality during the Cambrian radiation. This increase of provinciality is paralleled by a decrease of geographic ranges, which is consistent with earlier findings (e.g., Hendricks et al. 2008). A rise of global diversity can be fueled by an increase in provincialism (Valentine et al. 1978; Lieberman 2003; Stigall et al. 2017). However, we found no correlation between provinciality and diversity over the entire Cambrian.

Our trilobite biogeographic partitioning is largely concordant with traditional trilobite demarcations between Laurentia and Gondwana, but failed to delineate subrealms among tectonic blocks and microcontinents, such as Australia, Avalonia, India, Kazakhstan, and Mongolia (Fig. 3). Moreover, a biogeographic boundary between northern Laurentia and southern Laurentia, as suggested by earlier studies (Meert and Lieberman 2004; Álvaro et al. 2013), is not discernible during Series 2 and the Miaolingian in our trilobite partitioning outcome, although it is evident in our overall partitioning (Figs. 1, 3). This suggests that the dissimilarity in multi-taxonomic composition within Laurentia was not strongly associated with faunal changeover in local olenelloid trilobites (Pates et al. 2019). The differences between the traditional and network-based biogeographic models may reflect limitations of our approach, which is ignorant of phylogenetic histories.

Abrupt biogeographic differentiation during the Cambrian radiation has been suggested by vicariance associated with the breakup of supercontinents (Lieberman 2005; Meert and Lieberman 2008; Na and Kiessling 2015) such as Pannotia (Scotese 2009). The within-continent variation in bioregions was prevalent since Cambrian Stage 3 (Fig. 1), suggesting that the vicariance process may have occurred not only at continental scales but also at regional and local scales. Within-continent bioregionalization is also evident with different spatial resolution of grid cells (Supplementary Fig. S2) and supported by previous analyses (Lieberman 1998; Brock et al. 2000; Meert and Lieberman 2004; Skovsted 2004; Paterson et al. 2016).

Across the entire Phanerozoic, both climate and continental configuration have been important drivers of marine provinciality. There is evidence that climate might have been a stronger driver (Kocsis et al. 2021), with cooling leading to increased provincialism. Determining the relationship between climate and Cambrian provinciality is compromised by limited climate proxy data in the Cambrian. An early Cambrian greenhouse and subsequent cooling is suggested by a few mid-Cambrian geochemical data (Hearing et al. 2018; Wotte et al. 2019) and modeling (Hearing et al. 2021). Goldberg et al. (2021) demonstrated a late Cambrian warming, but there were no data points below the Miaolingian. It is thus unclear whether the increasing provinciality across the Cambrian was linked to a long-term cooling trend. The reduced temperature gradients in a warmer world, however, might in part explain the slight decrease in provincialism after the Wuliuan (Fig. 4).

The increase in the number of bioregions during the Cambrian radiation may have been facilitated by morphological and ecological innovations such as dispersal and predation, but underlying mechanisms may vary with time (Nürnberg and Aberhan 2013). The global distribution of a single bioregion (Fig. 1) and maximum geographic ranges of species in the Terreneuvian (Fig. 6) could result from a low-competition system (Na and Kiessling 2015), but could also be due to the lack of diagnostic characters in the earliest Cambrian faunas (Steiner et al. 2007). The recent discovery of a Cambrian bryozoan in traditional small shelly beds (Ernst and Wilson 2021; Zhang et al. 2021) proves a hidden diversity in the earliest Cambrian, which could also account for higher geographic occupancy and lower provinciality in the first two Cambrian stages than afterward. The development of novelties in trophic networks later in Stage 3, such as pervasive predation, would have caused niche overlap and contraction (Na and Kiessling 2015). This process would ultimately have decreased geographic ranges of species and increased biogeographic differentiation.

The metazoan larvae reported from abundant phosphatized microfossils in strata at the basal Cambrian (e.g., Chen et al. 2000; Steiner et al. 2004) indicate a larval revolution (Raff 2009). Cosmopolitanism in the Fortunian and Stage 2 could have been achieved by larval dispersal through passive transport (De Bie et al. 2012). Short larval durations under warm temperature (Álvarez-Noriega et al. 2020) might also be associated with the pattern that biogeographic boundaries occurred principally in low latitudes during Stages 3 and 4 (Figs. 1, 3, Supplementary Fig. S2).

Despite an increase in average mobility and perhaps in dispersibility of marine invertebrate species from the Fortunian to Stage 3 (Fig. 7), biogeographic patterns became more pronounced. This may suggest that geographic distribution of species is not dependent on their dispersal ability but on taxonomic duration during which they could persist (Kiessling and Aberhan 2007). Finding a direct link between dispersibility and geographic distribution could be compromised by our long-term timescale, which is in millions of years. But the biogeographic structure could be partly regulated by some widespread species. For example, the uniform biogeographic structures that were featured since the Wuliuan were likely related to the expanded distribution of the cosmopolitan agnostoid trilobites (Babcock et al. 2017), such as Ptychagnostus atavus and Glyptagnostus reticulatus (Zhu et al. 2006).

Another possible explanation for the increase in provinciality might be the large number of species with short stratigraphic duration (e.g., singletons) in the Cambrian (Supplementary Fig. S4), as species might not have attained their potential geographic range simply because they had short life spans (Kiessling and Aberhan 2007; Foote et al. 2008). The high proportion of singletons in low latitudes and shorter mean durations (Table 1) suggest a more rapid evolutionary turnover of tropical species, which would reduce geographic ranges and accelerate biogeographic differentiation.

Extinction events strongly affect biogeographic patterns (Kocsis et al. 2018a), and our results suggest that the Cambrian is no exception. Several biological crises are evident in the Cambrian. The best known is the extinction of archaeocyath sponges in Stage 4 (Sinsk event), the end of the Marjumiid Biomere at the Guzhangian/Paibian boundary, and the end of the Pterocephaiid Biomere at the Paibian/Jiangshanian boundary (e.g., Babcock et al. 2015). In the aftermath of the extinction of archaeocyath sponges (Zhuravlev and Wood 1996), a major biogeographic turnover occurred in the Wuliuan (Fig. 2), suggesting that the collapse of reefs exerted a strong influence on biogeographic structures. The other two events paved the way for a more biologically uniform world in the subsequent stages (Paibian, Jiangshanian, and Stage 10), when expansions of oxygen minimum zones, with their corresponding habitat disruption, were probably recurrent across the shallow shelf (Kröger et al. 2019).

Although our network-based approaches have removed some biases arising from an incomplete fossil record and uneven spatial coverage, a series of problems still remain, such as taxonomic issues and stratigraphic resolution. For example, the long stage duration (e.g., Fortunian) might lead to an underestimation of bioregions and an overestimation of the frequency of singletons compared with short stage duration (e.g., Paibian). If evolutionary rates were constant through time, differences in stage duration would severely distort biogeographic structures. However, in our case, given that the Terreneuvian is characterized by low taxonomic turnover rates and high genus survivorships (Kröger et al. 2019), more time slices cut into the Terreneuvian are unlikely to change our biogeographic patterns. On the other hand, Cambrian fossil groups may have been taxonomically oversplit, such as in small shelly faunas (Steiner et al. 2007), trilobites (Geyer 2019; Zhao et al. 2020), archaeocyath sponges. (Zhuravlev and Wood 2020). To what extent this taxonomic issue would impact our biogeographic patterns through time will rely on the variation in proportion of the oversplit taxa among stages, which is worth exploring in the future.

In summary, an increase in regional differences of faunal composition occurred during the Cambrian radiation and was linked with a decrease in by-species geographic distribution. We show that the various ways in which faunal composition was partitioned among regions may have been determined by geography, biotic interaction, climatic changes, and tectonics, providing further evidence for the evolution of biogeographic patterns during and after the Cambrian radiation.

We thank M. Hopkins, B. S. Lieberman, and C. M. Ø. Rasmussen for their insightful and helpful reviews, which greatly improved this article. L.N. and Q.L. acknowledge funding from the Youth Innovation Promotion Association of the Chinese Academy of Sciences (CAS; 2019310), and the State Key Laboratory of Paleobiology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, CAS) (20202103), and the CAS (XDB26000000). A.T.K. acknowledges funding from the Deutsche Forschungsgemeinschaft (Ko 5382/1-2 and Ko 5382/2-1). Q.L. is grateful for the invitation to be a guest researcher at GeoZentrum Nordbayern. This is a contribution to IGCP Project 668: “The Stratigraphic and Magmatic History of Early Palaeozoic Equatorial Gondwana and Its Associated Evolutionary Dynamics and IGCP 735 “Rocks and the Rise of Ordovician Life: Filling Knowledge Gaps in the Early Palaeozoic Biodiversification.” This is Paleobiology Database publication no. 436.

The authors declare no competing interests.

Raw data and scripts required to replicate the analysis are accessible through Zenodo: https://doi.org/10.5281/zenodo.6898357. Supplementary File F1: Cambrian biostratigraphic correlation; Supplementary File F2: Supplementary Figs. S1–S4 and Supplementary Table S1.