Non-technical Summary
Mass extinctions do more than eliminate large numbers of species: they also reshape how life is distributed across the planet. In this study, we examine how the Late Ordovician mass extinction, one of the most severe crises in Earth history, altered the global arrangement of marine communities. Using a large fossil dataset, we tracked changes in taxonomic diversity and how different regions shared their faunas before, during, and after the extinction. We find that the extinction not only reduced global biodiversity but also broke down long-standing biogeographic patterns. Before the crisis, most genera had geographically restricted distributions. After the extinction, however, many genera became widespread, and distinct regional communities collapsed into a more globally uniform biota. This transient homogenization happened because many locally restricted genera died out, while a small number of environmentally tolerant survivors expanded their ranges across previously effective barriers. Although global biodiversity recovered rapidly in the early Silurian, the reestablishment of distinct regional ecosystems took several million years. Our results suggest that short-lived global homogenization followed by a slow rebuilding of regional differences may represent a common ecological response to major environmental upheavals.
Introduction
Mass extinctions are among the most profound disruptions in Earth’s history, reshaping both the taxonomic composition and the spatial organization of life (Jablonski Reference Jablonski2001). The Late Ordovician mass extinction (LOME) was the first of the five commonly accepted major extinction events of the Phanerozoic (Raup and Sepkoski Reference Raup and Sepkoski1982; Marshall Reference Marshall2023). The LOME may have eradicated between 72% (Stanley Reference Stanley2016) and 85% of marine species (Harper et al. Reference Harper, Hammarlund and Rasmussen2014). The crisis occurred in at least two pulses (Harper et al. Reference Harper, Hammarlund and Rasmussen2014; Rasmussen et al. Reference Rasmussen, Kröger, Nielsen and Colmenar2019): the first near the Katian/Hirnantian boundary, broadly attributed to glacial cooling (Finnegan et al. Reference Finnegan, Rasmussen and Harper2016; Wang et al. Reference Wang, Zhan and Percival2019), and the second in the mid-Hirnantian deglaciation, linked to expanded oceanic anoxia (Stockey et al. Reference Stockey, Cole, Planavsky, Loydell, Frýda and Sperling2020).
The LOME triggered profound faunal turnover across multiple clades (Darroch and Wagner Reference Darroch and Wagner2015; Finnegan et al. Reference Finnegan, Rasmussen and Harper2016). Biogeographic patterns of diversification across the Ordovician/Silurian boundary have been well explored (Harper Reference Harper2024), but the consequences of LOME for global biogeographic structure (e.g., the spatial differentiation and connectivity of marine communities) remain poorly constrained.
The dominant hypothesis is that mass extinctions cause biotic homogenization by selectively eliminating endemic taxa and promoting cosmopolitanism through survivor dispersal, as observed in Early Triassic ammonoids and bivalves in the aftermath of the end-Permian mass extinction (Brayard et al. Reference Brayard, Escarguel and Bucher2007; Echevarría and Ros-Franch Reference Echevarría and Ros-Franch2024). This event is known to have led to marine homogenization via high-latitude endemic extinction and the poleward expansion of tropical survivors (Al Aswad et al. Reference Al Aswad, Penn, Monarrez, Bazzi, Deutsch and Payne2025). Loss of provincialism was followed by gradual reestablishment of biogeographic differentiation with increasing endemism, such as in the Middle to Late Triassic (Echevarría and Ros-Franch Reference Echevarría and Ros-Franch2024). In contrast to such observations, Antell et al. (Reference Antell, Kiessling, Aberhan and Saupe2020) have argued that there is limited evidence for ecological release after mass extinction events and claimed that biodiversity and biogeographic distributions are independent at large temporal scales.
Previous work on Ordovician–Silurian biogeography offers mixed insights. One quantitative analysis suggested that major biogeographic reorganizations were limited to the first extinction pulse at the Katian/Hirnantian boundary (Kocsis et al. Reference Kocsis, Reddin and Kiessling2018). Another study documented that changes in beta diversity were mainly caused by the selective loss of taxa with very small geographic ranges, while different extinction pulses (cooling followed by warming) differently affected alpha and beta diversity by altering dispersal barriers (Darroch and Wagner Reference Darroch and Wagner2015). Heterogeneous spatial sampling complicates attempts to compare biogeographic structure through time, underscoring the need for a standardized framework to assess biogeographic changes (Kocsis et al. Reference Kocsis, Reddin, Scotese, Valdes and Kiessling2021).
Here, we combined standardized diversity partitioning with bipartite network analysis to quantify first-order biogeographic change across the Ordovician/Silurian boundary. Specifically, we evaluated whether the LOME drove global biogeographic homogenization (a collapse of provincial differentiation) and whether spatial reorganization lagged taxonomic recovery.
We addressed the following three hypotheses: (1) Marine provincialism decreased in the aftermath of crises, accompanied by fragmentation and subsequent reconfiguration of the global biogeographic network through time. (2) Spatial faunal differentiation (beta diversity) declined during this interval, reflecting increased biotic homogenization (Darroch and Wagner Reference Darroch and Wagner2015). (3) Shifts in the relative dominance of cosmopolitan versus endemic genera strongly influenced large-scale biogeographic patterns (Huang et al. Reference Huang, Rong and Cocks2012, Reference Huang, Jin and Rong2018).
Materials and Methods
Data Compilation and Cleaning
We downloaded all Ordovician–Silurian marine invertebrate occurrence data from the Paleobiology Database (PBDB; accessed 26 February 2023) and analyzed the data with scripts written in the R programming environment (R Core Team 2023). To ensure taxonomic and stratigraphic consistency, we cleaned the data using the fossilbrush R package (Flannery-Sutherland et al. Reference Flannery-Sutherland, Raja, Kocsis and Kiessling2022), which detects and resolves anomalies such as stratigraphic misplacements, misspellings, outdated names, and missing classifications.
Stratigraphic validation followed the Sepkoski (Reference Sepkoski2002) compendium as a reference. Occurrences were flagged when observed stratigraphic ranges overlapped <75% with reference intervals. Anomalies included first appearance older and/or last appearance younger than the reference. Additionally, the Pacmacro algorithm (default parameters: tail density = 5%, tail stratigraphic proportion = 40%) identified and truncated artificially extended ranges. For example, in the conodont genus Wurmiella, several occurrences with broad or dispersed stratigraphic assignments were flagged as tail records by Pacmacro, illustrating how the method reduces the influence of stratigraphically diffuse occurrences on downstream range checking. The final curated dataset comprises 26,532 occurrences, 1149 collections, and 2474 genera spanning the Sandbian–Telychian interval (Late Ordovician–early Silurian).
Temporal and Spatial Binning
To improve temporal resolution beyond standard PBDB stages, we subdivided the Katian into two substages following Rasmussen et al. (Reference Rasmussen, Kröger, Nielsen and Colmenar2019), which resulted in seven time bins: Sandbian, Katian 1, Katian 2, Hirnantian, Rhuddanian, Aeronian, and Telychian. The mean duration of these intervals is 3.57 ± 1.45 Myr, reflecting natural variation in the temporal lengths of the constituent stratigraphic stages. While a finer stratigraphic subdivision is possible (Rasmussen et al. Reference Rasmussen, Kröger, Nielsen and Colmenar2019), we used this subdivision as a compromise between bins of similar duration and sufficient sample size for meaningful biogeographic analysis.
Each collection was assigned to a tectonic plate (geoplate) based on PBDB metadata. A small proportion of collections (n = 354, 2.11% of the original dataset) required correction of plate assignments based on published literature and paleogeographic context. Paleogeographic coordinates were obtained from the PBDB. For each collection, reconstructed paleo-coordinates were taken from the PBDB fields paleolat and paleolng, and tectonic plate assignments from the field geoplate. These paleo-coordinates are derived from global plate tectonic reconstructions implemented in GPlates and are based on the model of Wright et al. (Reference Wright, Zahirovic, Müller and Seton2013).
Collections were aggregated into equal-area hexagonal grid cells using the icosa R package (Kocsis et al. Reference Kocsis, Reddin, Scotese, Valdes and Kiessling2021). We analyzed two spatial resolutions: (1) coarse scale: ~252 cells per interval (2 × 10⁶ km² per cell, 8° edge length); and (2) fine scale: ~812 cells per interval (6.3 × 10⁵ km² per cell). The main text presents results based on the fine spatial grain, while the corresponding analyses at the coarse spatial grain are provided in Supplementary Figure S1. Both spatial resolutions yield qualitatively consistent patterns. Equal-area hexagonal grid cells were used in the network and occupancy analyses, where grid cells served as spatial units for constructing paleobiogeographic networks and calculating network topology metrics.
Sampling Standardization and Diversity Partitioning
Fossil occurrence data are inherently unevenly distributed in space and time (Alroy et al. Reference Alroy, Marshall, Bambach, Bezusko, Foote, Fürsich and Hansen2001). To correct for these effects, we applied sampling standardization. Following Alroy et al. (Reference Alroy, Marshall, Bambach, Bezusko, Foote, Fürsich and Hansen2001), we applied occurrence-weighted subsampling (OW; quota = 160 occurrences). In addition, following the approach of Darroch and Wagner (Reference Darroch and Wagner2015), we implemented a stratigraphically structured subsampling procedure in which occurrences were first partitioned into distinct stratigraphic units (formations), and an equal number of occurrences (6 formations × 10 occurrences per formation) was drawn from each unit. This approach reduces the overrepresentation of heavily sampled formations and ensures more balanced spatial and stratigraphic sampling across intervals. For each subsampling method, we applied a fixed quota across all time bins to standardize sampling intensity, while 500 iterations (trials) were performed to stabilize the estimates. Taxa represented by a single occurrence were omitted from all analyses.
By partitioning diversity into α (alpha, local), β (beta, between-community), and γ (gamma, global) components, we established a framework for disentangling the underlying structure and dynamics of biodiversity patterns. In this study, alpha, beta, and gamma diversity were consistently calculated at the formation level to ensure comparability across diversity components. Beta diversity was quantified using the Sørensen and Jaccard dissimilarity indices and partitioned into turnover and nestedness components after Baselga (Reference Baselga2010). Turnover (βsim) captures taxonomic replacement among regions, while nestedness (βnes) reflects differences, in which taxon-poor assemblages represent subsets of taxon-rich assemblages. Calculations were performed using the betapart R package (Baselga and Orme Reference Baselga and Orme2012). To evaluate the robustness of our results, parallel analyses based on individual fossil collections were also performed and are presented in the Supplementary Material. To assess whether observed patterns of beta diversity deviate from random expectations, we performed null-model analyses using the Raup-Crick metric implemented in the raupcrick() function in the vegan package (Chase et al. Reference Chase, Kraft, Smith, Vellend and Inouye2011). This probabilistic measure accounts for differences in taxonomic richness by comparing observed dissimilarity to a null distribution generated under random assembly. By comparing temporal trajectories of alpha, beta, and gamma diversity, we evaluated whether the LOME was associated with a tendency toward faunal homogenization (as reflected by declines in beta diversity) and how postextinction recovery unfolded. Because beta diversity alone cannot delineate provincialization, these patterns were interpreted as complementary evidence and further examined using the network-based biogeographic analyses described in the “Network Construction and Metrics” section.
Network Construction and Metrics
To quantify biogeographic structure, we constructed bipartite occurrence networks linking genera to geographic units (geoplates or grid cells) in each time bin (Wang and Huang Reference Wang and Huang2020). Each network node represents either a genus or a geographic unit, and edges indicate occurrences. Networks were analyzed using the Louvain community detection algorithm (Blondel et al. Reference Blondel, Guillaume, Lambiotte and Lefebvre2008) to identify cohesive modules and characterize overall topology.
In this study, we distinguished three related but non-equivalent concepts. Modularity describes how strongly genera and geographic units form cohesive, separated modules. Regional differentiation refers to compositional differences among regions, quantified by beta diversity and its turnover and nestedness components. Provinciality denotes the broader pattern of distinct faunal assemblages across regions. High modularity and regional differentiation together indicate stronger provincial structure. We calculated four structural metrics: (1) Degree centrality is the average number of edges per node, indicating network connectivity (Freeman Reference Freeman1979; Dunne et al. Reference Dunne, Williams and Martinez2002). (2) Modularity is the degree to which the network subdivides into cohesive modules (Thébault and Fontaine Reference Thébault and Fontaine2010; Sales-Pardo Reference Sales-Pardo2017). High modularity values denote strong provinciality. (3) Average path length is the mean shortest distance between node pairs, reflecting faunal isolation and connectivity. Higher values indicate greater faunal isolation and weaker connectivity, whereas lower values indicate stronger connectivity among regions (Watts and Strogatz Reference Watts and Strogatz1998; Guimerà and Amaral Reference Guimerà and Amaral2005). Declining modularity and shorter path lengths are interpreted as reduced provincial differentiation and increased cosmopolitan connectivity. (4) Nestedness based on overlap and decreasing fill (NODF) measures the extent to which the genera occurring in less diverse geographic units form subsets of those in more diverse units (Almeida-Neto et al. Reference Almeida-Neto, Guimaraes, Guimaraes, Loyola and Ulrich2008). All metrics were computed using the igraph package (Csárdi and Nepusz Reference Csárdi and Nepusz2006) and averaged over 500 randomized subsampling iterations.
Network analysis has been previously applied to fossil occurrence data (Kiel Reference Kiel2017; Rojas et al. Reference Rojas, Patarroyo, Mao, Bengtson and Kowalewski2017; Kocsis et al. Reference Kocsis, Reddin and Kiessling2018), but our approach extends this methodology in three key directions. We integrate subsampling to control for sampling bias, combine network metrics with diversity partitioning to bridge taxonomic and spatial scales, and explicitly trace the temporal reorganization of networks through crises.
All diversity and network metrics are reported as means and standard deviations across subsampling replicates. Temporal changes are expressed as Δ values, defined as differences in mean values between successive time bins.
Geographic Occupancy
To assess changes in the spatial breadth of genera, we calculated geographic occupancy (Kiessling and Kocsis Reference Kiessling and Kocsis2016) for each time bin as the proportion of sampled geographic equal-area grid cells in which a genus occurred. While this metric is not strictly a measurement of range size, it is less sensitive to sampling bias than, for example, maximum great circle distance (Kiessling and Kocsis Reference Kiessling and Kocsis2016). Changes in occupancy through time were used to track shifts in the spatial breadth of taxa and the degree of biogeographic structure. Pairwise Wilcoxon tests evaluated whether occupancy shifts between intervals were statistically significant. Combining occupancy data with network metrics allowed us to link taxonomic survivorship and ecological breadth to the structural evolution of marine biogeography.
Paleogeographic Mapping of Extinction Intensity
To visualize spatial variation in extinction intensity, occurrence data were aggregated into 10° × 10° paleolatitudinal–paleolongitudinal cells following Raup and Jablonski (Reference Raup and Jablonski1993). Genus-level extinction intensity was calculated for each cell and time interval, following Foote (Reference Foote2000). Only cells containing at least two sampled genera were included. Paleogeographic coordinates were plotted on PALEOMAP reconstructions from Scotese (Reference Scotese2016). Because estimates from sparsely sampled cells are unstable, the number of sampled genera is displayed in each cell, and the maps are interpreted descriptively rather than as sampling-standardized tests.
Analytical Predictions
To explicitly link hypotheses with analytical metrics, we evaluated each hypothesis using corresponding predictions. For Hypothesis 1, which predicts reduced and reconfigured marine provincialism after the crisis, we used network metrics to assess provincial structure, isolation, connectivity, and ordered subset structure. Higher modularity indicates stronger provincial structure, while lower values indicate weaker partitioning. Longer average path length reflects greater isolation, whereas shorter values indicate stronger connectivity. Higher degree centrality reflects greater faunal sharing, and higher nestedness (NODF) indicates stronger ordered subset structure. Thus, decreased modularity and average path length, combined with increased degree centrality and NODF, support reduced provincialism and network reconfiguration.
For Hypothesis 2, predicting declining spatial differentiation and increased biotic homogenization, we used total beta diversity and its turnover and nestedness components. Decreases in total beta diversity and turnover indicate increasing homogenization, whereas increases suggest stronger regional differentiation. The nestedness component assesses whether differences are driven by ordered richness loss rather than taxonomic replacement.
For Hypothesis 3, which predicts that shifts between cosmopolitan and endemic genera shaped large-scale biogeographic patterns, we examined changes in geographic occupancy of surviving genera. Increased occupancy after the extinction interval indicates broader ranges and greater influence of cosmopolitan taxa, especially when accompanied by reduced provincial structure and increased faunal sharing.
Results
Alpha and Gamma Diversity Dynamics
Following sampling standardization, local and global genus richness remained stable through the Sandbian–Katian interval and increased slightly in the late Katian (Fig. 1A). Alpha diversity remained stable in the Hirnantian (9.76 ± 0.22) but declined in the Rhuddanian (9.43 ± 0.32), whereas global diversity had declined already in the Hirnantian (52.82 ± 3.49) and plummeted further in the Rhuddanian (47.56 ± 3.49). A rebound of both metrics occurred in the Aeronian and Telychian.
Sampling-standardized diversity patterns across the studied interval. A, Mean alpha and total gamma diversity. B, Beta diversity (Sørensen dissimilarity and Jaccard dissimilarity; by Baselga Reference Baselga2010). C, Nestedness (Sørensen-based Nestedness component and Jaccard-based Nestedness component; by Baselga Reference Baselga2010). D, Turnover (Sørensen-based Turnover component and Jaccard-based Turnover component; by Baselga Reference Baselga2010). The stratigraphic positions of the two main extinction pulses are indicated following the marking style of Darroch and Wagner (Reference Darroch and Wagner2015). San, Sandbian; Kat1, Katian 1; Kat2, Katian 2; Hir, Hirnantian; Rhud, Rhuddanian; Aero, Aeronian; Tely, Telychian.
Figure 1. Long description
A multi-panel line graph with four panels labeled A through D. All panels share a common X-axis representing Age in M a, ranging from 455 to 435, with stratigraphic stages labeled San, Kat 1, Kat 2, Hir, Rhud, Aero, and Tely. Two vertical dashed lines marked with stars and numbers 1 and 2 indicate extinction pulses during the Hir stage.
Panel A: Mean alpha diversity (orange line) and total gamma diversity (green line). Both lines remain relatively stable until the Hir stage, where they drop sharply to a trough in the Rhud stage before recovering. Gamma diversity ranges from 45 to 60, and Alpha diversity ranges from 9.0 to 10.5.
Panel B: Beta diversity using Jaccard (red) and Sorenson (blue) dissimilarities. Both metrics show a slight decline starting at extinction pulse 1, reaching a minimum in the Rhud stage, followed by a gradual recovery. Jaccard values are higher, near 0.98, while Sorenson values are near 0.96.
Panel C: Nestedness using Jaccard (red) and Sorenson (blue) components. Unlike other metrics, nestedness shows a slight peak during the Rhud stage, rising from approximately 0.008 to 0.011 before leveling off.
Panel D: Turnover using Jaccard (red) and Sorenson (blue) components. The trends mirror Panel B, showing a distinct dip during the extinction interval (Hir to Rhud) with values ranging between 0.94 and 0.98.
Beta Diversity and Its Components
Beta diversity declined markedly across the Hirnantian/Rhuddanian transition (Δ = −0.006) (Fig. 1B). Decomposition of beta diversity reveals that this reduction was driven almost entirely by a sharp drop (Δ = −0.009) in turnover rather than by a change of nestedness, which increased (Δ = +0.003) (Fig. 1C,D). Nestedness contributed little to overall dissimilarity. This decline was also recovered using the richness-corrected Raup–Crick dissimilarity metric, indicating that the observed reduction in faunal differentiation cannot be attributed solely to changes in taxonomic richness (Supplementary Fig. S3). To assess the robustness of these patterns, we applied two independent subsampling approaches (occurrence-weighted and formation-based). Comparable diversity patterns were recovered at the alternative spatial resolution (Supplementary Fig. S2). Both methods yielded consistent temporal trends in beta diversity and its components (see Supplementary Material), indicating that the observed patterns are unlikely to be sampling artifacts.
Spatial Distribution of Extinction Intensity
Extinction intensity showed substantial spatial heterogeneity across the LOME (Fig. 2). No simple global latitudinal gradient was evident in Katian 2. During the Hirnantian, relatively high extinction intensities occurred in several low-latitude cells associated with Laurentia, whereas some midlatitude cells, including parts of Baltica, showed lower values. However, spatial coverage was uneven, and several cells contained relatively few sampled genera. The maps therefore provide descriptive evidence for spatial heterogeneity rather than a formal test of latitudinal selectivity.
Mapping extinction intensity in the Late Ordovician mass extinction. A, Extinction intensity during Katian 2. B, Extinction intensity during the Hirnantian (latest Ordovician). C, Extinction intensity during the Rhuddanian (earliest Silurian). Grid cells are defined at 10° intervals of paleolongitude and paleolatitude. White numbers within each cell denote the number of genera sampled in that geographic region. Cell color represents the genus-level extinction rate, calculated following the method of Foote (Reference Foote2000), with the color scale provided at the bottom. Paleogeographic base maps are from the PALEOMAP raster series (Scotese Reference Scotese2016).
Figure 2. Long description
A three-panel series of global paleogeographic maps. The Y-axis represents Paleolatitude from -50 to 50, and the X-axis represents Paleolongitude from -100 to 100. A color scale at the bottom indicates extinction rates from 0.5 (dark blue) to 2.0 (bright yellow).
* Panel A (Katian 2): Data cells are densely clustered around the central and eastern paleocontinents. Most cells are colored magenta to light purple, indicating moderate extinction rates. High genus counts (e.g., 227, 150, 123) are concentrated in the central tropical regions. A single dark blue cell (low rate) is located in the North-central region with a genus count of 8.
* Panel B (Hirnantian): Data cells are more sparse but show significantly higher extinction rates, indicated by orange and yellow colors. A bright yellow cell (highest rate) is located near the equator at 0 degrees longitude. Orange cells are visible on the western and central paleocontinents. Genus counts are lower than in Panel A, ranging from 9 to 127.
* Panel C (Rhuddanian): Data cells are the most sparse of the three panels. The colors have shifted back toward dark blue and purple, indicating a significant drop in extinction rates following the mass extinction event. Cells are scattered across the central and eastern regions with genus counts ranging from 6 to 58.
Changes in Paleobiogeographic Network Structure
Network metrics corroborate the previous results (Fig. 3, Supplementary Fig. S4): Modularity (Fig. 3C) declined from the Katian 2 onward and reached a minimum in the Hirnantian (0.58 ± 0.071), remaining similarly low in the Rhuddanian (0.60 ± 0.040). It then rose gradually through the Aeronian–Telychian but never returned to pre-extinction levels. Average path length (Fig. 3A) shortened continuously toward the Hirnantian minimum (3.15 ± 0.38), indicating increased connectivity among geographic regions. Degree centrality (number of ties per node) rose concurrently, reflecting enhanced faunal overlap among regions (Fig. 3B). Nestedness (NODF) showed a modest increase (Δ = +0.025) during the Hirnantian–Rhuddanian (Fig. 3D).
Biogeographic metrics across the studied interval. A, Average path length: mean shortest distance between nodes, reflecting faunal isolation and connectivity. B, Degree centrality: average number of edges per node, indicating network connectivity and potential for widespread dispersal. C, Modularity: degree to which the network subdivides into cohesive modules; high values denote strong provinciality. D, NODF (nestedness metric based on overlap and decreasing fill): measures the extent to which less diverse sites contain subsets of genera from more diverse sites, indicating nestedness of communities. The graphs are based on an average partitioning of geographic area (resolution = ~812 cells per interval, 6.3 × 10⁵ km² per cell, data with 26,532 occurrences, 1149 collections, 2474 genera) and sampling standardization (occurrence-weighted subsampling [OW], quota = 160, trials = 500). San, Sandbian; Kat1, Katian 1; Kat2, Katian 2; Hir, Hirnantian; Rhud, Rhuddanian; Aero, Aeronian; Tely, Telychian.
Figure 3. Long description
A four-panel grid of line graphs labeled A through D. All graphs share a common X-axis representing Age in M a, ranging from 455 to 435, with geological intervals labeled San, Kat 1, Kat 2, Hir, Rhud, Aero, and Tely. Two vertical dashed lines marked with stars 1 and 2 indicate the Hirnantian interval.
* Panel A (Top-Left): Average path length in green. The Y-axis ranges from 2.0 to 5.0. The line fluctuates between 3.5 and 4.0, showing a sharp dip to its lowest point during the Hir interval before rising again.
* Panel B (Top-Right): Degree centrality in blue. The Y-axis ranges from 0.0 to 0.8. The trend shows a gradual increase, peaking at approximately 0.5 during the Hir interval, followed by a decline and a final rise in the Tely interval.
* Panel C (Bottom-Left): Modularity in red. The Y-axis ranges from 0.2 to 1.0. The line shows a slight downward trend overall, with a notable trough during the Hir interval and a small peak in the Aero interval.
* Panel D (Bottom-Right): N O D F in purple. The Y-axis ranges from 0.00 to 0.30. The data shows a fluctuating but generally increasing trend, with a peak following the Hir interval and a sharp dip in the Aero interval.
Each line is surrounded by a shaded area representing a confidence interval or error margin.
Changes in Genus-Level Occupancy
Genus-level occupancy patterns confirm these biogeographic changes (Fig. 4). In the Sandbian–Katian interval, most genera exhibited low occupancy, consistent with strong provincial structure. Median occupancy increased slightly from the Katian into the Hirnantian and then rose markedly across the Hirnantian/Rhuddanian transition (Δ = +0.071), indicating a substantial expansion in geographic breadth. Occupancy values remained elevated through the Aeronian before declining in the Telychian. These patterns indicate a transient increase in geographic dispersion following the extinction interval, followed by a partial return toward more restricted distributions. The same temporal pattern is evident when all genera are analyzed (Supplementary Fig. S5) and when the analysis is restricted to genera surviving between consecutive bins (Fig. 4).
Geographic occupancy of surviving genera across the studied intervals. Gray shaded bands indicate the two extinction pulses. Boxes and violin plots of the same color represent the geographic occupancy of the same cohort of surviving genera tracked across two consecutive bins, illustrating changes in their occupancy through time. The black lines indicate medians, and the boxes represent interquartile ranges. Pairwise Wilcoxon signed-rank tests were used, because the same surviving genera were compared between consecutive intervals. Arrows indicate increases and decreases of occupancy, respectively. Asterisks indicate significance levels: *0.01< p < 0.05; **0.001< p < 0.01; ***p < 0.001. Data and treatment as in Fig. 3.
Figure 4. Long description
The graph plots Geographic occupancy on the Y-axis from 0.0 to 0.9 against geological stages on the X-axis.
From left to right, the stages are Sandbian, Katian 1, Katian 2, Hirnantian, Rhuddanian, Aeronian, and Telychian. Two gray vertical bands highlight the mass extinction interval covering the Hirnantian and Rhuddanian stages.
Data is presented in paired cohorts of violin plots and narrow box plots.
- Sandbian to Katian 1: Teal plots show a slight decrease in median occupancy.
- Katian 1 to Katian 2: Orange plots show a further decrease.
- Katian 2 to Hirnantian: Blue plots show a significant increase, marked with three asterisks and an upward arrow.
- Hirnantian to Rhuddanian: Pink plots show a significant increase, marked with three asterisks and an upward arrow.
- Rhuddanian to Aeronian: Green plots show a significant decrease, marked with three asterisks and a downward arrow.
- Aeronian to Telychian: Yellow plots show a significant decrease, marked with three asterisks and a downward arrow.
At the top of the plot, horizontal brackets connect consecutive stages. Statistical significance is indicated by three asterisks (p less than 0.001) for most transitions, accompanied by arrows indicating the direction of change in occupancy. The legend on the right identifies the color-coded stage comparisons.
Discussion
Biogeographic Restructuring during the LOME
The LOME profoundly reorganized the spatial structure of marine biodiversity. Our analyses show that biogeographic change across the LOME was not a simple, monotonic shift toward global homogenization, but a phase-dependent restructuring involving faunal turnover, reduced provinciality, and transient ecological release among survivors. This result extends previous work by combining global-scale PBDB data, taxonomically inclusive marine invertebrate records, beta diversity partitioning, network structure, and survivor occupancy analyses. Our findings partially align with previous studies. Darroch and Wagner (Reference Darroch and Wagner2015) demonstrated that the crisis markedly reduced provincial differentiation and increased faunal connectivity, while Harper (Reference Harper2024) highlighted the preferential loss of environmentally restricted taxa and the rise of cosmopolitan survivors. More broadly, similar patterns of biogeographic homogenization associated with extinction events have been identified at larger temporal scales (Kocsis et al. Reference Kocsis, Reddin and Kiessling2018). The present analysis complements these works by quantitatively linking selective extinction to structural changes in biogeographic organization. This integration of network and occupancy approaches refines our understanding of how the LOME not only reduced diversity but also reshaped the spatial architecture of global marine ecosystems.
Hierarchical Diversity Collapse and Complexity Lag
The observed pattern suggests that the end-Ordovician crisis triggered an early decline in global diversity into the Hirnantian, followed by a collapse of local diversity into the Rhuddanian. Our results are consistent with previous studies showing a prolonged biodiversity decline across the LOME. Rasmussen et al. (Reference Rasmussen, Kröger, Nielsen and Colmenar2019) demonstrated a sustained reduction in global richness extending into the early Silurian, while Harper et al. (Reference Harper, Hammarlund and Rasmussen2014) highlighted the two-phase extinction dynamics and the magnitude of taxonomic loss during the Hirnantian. Extinction selectivity against specialists with narrow thermal, oxygen, or habitat tolerances was pronounced in each pulse, with the first pulse preferentially affecting geographically restricted or depth-restricted taxa, and the second pulse disproportionately impacting local specialists (Finnegan et al. Reference Finnegan, Rasmussen and Harper2016). The persistence of ecological generalists and depth-tolerant taxa supports the “generalist advantage” model (Huang et al. Reference Huang, Harper, Zhou and Rong2020), highlighting that survival depended on ecological breadth and habitat flexibility. Although taxonomic richness rebounded in the Aeronian, the recovery of beta diversity and network modularity lagged behind, remaining suppressed through the Rhuddanian and Aeronian and only returning toward pre-extinction levels in the Telychian (Figs. 1, 3).
Such lagging of taxonomic and biogeographic rebound is consistent with patterns following other mass extinctions. Ecological reorganization typically lags behind taxonomic rebound (Erwin Reference Erwin2001) and is spatially heterogeneous (Lowery et al. Reference Lowery, Bralower, Owens, Rodríguez-Tovar, Jones, Smit and Whalen2018). The delayed recovery observed here likely reflects the time needed to reestablish biogeographic differentiation through evolutionary innovation and renewed environmental gradients.
Beta Diversity Decline and Faunal Homogenization
The sharp reduction in beta diversity during the extinction interval was driven almost entirely by decreased spatial turnover rather than nestedness and was thus mainly driven by reduced taxonomic replacement among regions rather than richness gradients. This result mirrors modern biotic homogenization processes following environmental stress or habitat loss (Olden et al. Reference Olden, Poff, Douglas, Douglas and Fausch2004). The rise of cosmopolitan survivors, such as components of the Hirnantian fauna (Rong et al. Reference Rong, Chen and Harper2002), eroded regional distinctiveness and reduced faunal replacement among provinces. Occupancy trends mirror the modularity and beta diversity results, showing that the LOME replaced regionally distinct faunas with a smaller set of widespread survivors, and that geographic re-provincialization lagged behind taxonomic recovery.
In contrast to our finding of a further reduction of beta diversity, Darroch and Wagner (Reference Darroch and Wagner2015) assessed brachiopod patterns and suggested that beta diversity increased following the second extinction pulse. To assess whether these differences arise from analytical methods or taxonomic scope, we repeated the analyses using brachiopods only (Supplementary Fig. S6). The resulting beta diversity trajectory is consistent with Darroch and Wagner (Reference Darroch and Wagner2015), showing an increase following the second extinction pulse. In addition, network modularity reached a minimum during the Hirnantian (Supplementary Fig. S7), indicating a transient reduction in biogeographic structure. These results indicate that the differences between previous studies and our findings primarily reflect contrasts in taxonomic scope and spatial scale. Future studies need to explore why brachiopods responded differently from the total benthic fauna.
Earlier studies analyzed single clades or regional datasets, whereas the present study integrates multiple marine invertebrate groups across a broader spatial framework. For example, regional analyses from Laurentia indicate that postextinction recovery was rapid and strongly influenced by immigration and faunal turnover, leading to a decoupling between regional and global diversity patterns (Krug et al. Reference Krug, Patzkowsky and Knoll2004). These results highlight that regional diversity can be maintained or quickly restored through dispersal processes even when global diversity remains depressed. Regional recovery dynamics may differ from global patterns, with immigration and faunal turnover contributing to early Silurian differentiation. These factors suggest that postextinction biogeographic restructuring can vary depending on the taxa analyzed and geographic scale.
Network Collapse and Reassembly
Network analysis provides an independent structural perspective of biotic homogenization during the LOME. The transition from a modular, compartmentalized network to one dominated by highly connected nodes suggests global biogeographic homogenization. Analogous to “spatial generalists” in ecological networks, the surviving genera, which attained cosmopolitan distributions in the aftermath of the event, functioned as “hub taxa,” bridging formerly isolated provinces and mediating global dispersal (Stigall Reference Stigall2018).
The gradual reappearance of modularity during the early Silurian may reflect the progressive reestablishment of biogeographic structure linked to the reorganization of environmental barriers and climatic gradients associated with paleogeographic change (Kocsis et al. Reference Kocsis, Reddin, Scotese, Valdes and Kiessling2021; Scotese Reference Scotese2021). Here, biogeographic complexity refers to the degree of spatial structuring, including modularity, regional differentiation, and provinciality. The delayed recovery of these properties relative to taxonomic richness indicates a lag in the reassembly of spatial structure following the extinction. This pattern is consistent with postextinction recovery dynamics, in which ecospace and ecological structure are rebuilt gradually after an initial survival phase (Erwin Reference Erwin2001).
Environmental and Ecological Drivers
The two pulses of the LOME correspond to distinct but complementary stress regimes. The first, at the Katian/Hirnantian boundary, was dominated by glaciation-driven cooling and sea-level fall that compressed thermal niches and disrupted latitudinal diversity gradients (Finnegan et al. Reference Finnegan, Heim, Peters and Fischer2012; Wang et al. Reference Wang, Zhan and Percival2019). The second, near the Hirnantian/Rhuddanian boundary, was marked by widespread anoxia and expansion of euxinic conditions (Stockey et al. Reference Stockey, Cole, Planavsky, Loydell, Frýda and Sperling2020; Kozik et al. Reference Kozik, Young, Newby, Liu, Chen, Hammarlund, Bond, Them and Owens2022).
These contrasting regimes likely imposed different ecological effects. The initial glaciation-driven event primarily affected taxa through habitat contraction and climatic shifts, whereas the subsequent anoxia-driven pulse may have more strongly impacted taxa with high metabolic demand and restricted ecological tolerances, particularly deep-water and stenotopic benthos (Finnegan et al. Reference Finnegan, Rasmussen and Harper2016). However, the extent of ecological selectivity remains debated, with some studies emphasizing a broad ecological collapse, reduced provinciality, and preferential survival of cosmopolitan taxa following the second pulse (Brenchley et al. Reference Brenchley, Marshall and Underwood2001). Eurytopic genera may have preferentially persisted in intermittent oxygenated refugia (Rong and Zhan Reference Rong and Zhan2006), suggesting that survival was at least partly structured by ecological tolerance, even if broader patterns of ecological collapse and reduced provinciality also played an important role.
Mapping extinction intensity in paleogeographic space reveals an indistinct spatial pattern during the LOME (Fig. 2). Extinction in the late Hirnantian was higher in low-latitude regions (e.g., Laurentia) and lower in some midlatitude areas (e.g., Baltica), suggesting midlatitude refugia. This pattern is consistent with the expansion of the Hirnantia fauna (Rong and Zhan Reference Rong and Zhan2006) and increased biogeographic homogenization during the Hirnantian–Rhuddanian interval.
The mass extinction produced a transient biogeographic regime characterized by “crisis cosmopolitanism” sensu Button et al. (Reference Button, Lloyd, Ezcurra and Butler2017): the short-term expansion of tolerant, widespread taxa across eroded environmental and geographic barriers. This pattern resembles the Hirnantia fauna, a broadly distributed, low-diversity assemblage characteristic of the Hirnantian (Rong et al. Reference Rong, Harper, Huang, Li, Zhang and Chen2020), and may represent an early manifestation of crisis cosmopolitanism. However, the Hirnantia fauna primarily reflects environmental filtering during the extinction interval, whereas the main phase of ecological release and range expansion likely occurred in the Rhuddanian. The LOME did not have iconic disaster taxa like Claraia, which experienced a short-term global boom after the end-Permian mass extinction (Petsios and Bottjer Reference Petsios and Bottjer2016). Despite the absence of an iconic disaster taxon, the broader biogeographic signal is consistent with crisis cosmopolitanism. The subsequent reappearance of endemic lineages in the early Silurian marks the onset of re-provincialization, driven by renewed environmental heterogeneity and biotic interactions (Cai et al. Reference Cai, Snyder, Hastings and D’Souza2020; Muscente et al. Reference Muscente, Martindale, Prabhu, Ma, Fox, Hazen and Knoll2021).
Ecological release predicts that survivors may expand geographically or ecologically when competitors, predators, or incumbents are removed (Erwin Reference Erwin2001; Yoder et al. Reference Yoder, Clancey, Des Roches, Eastman, Gentry, Godsoe and Hagey2010). However, this mechanism should not be treated as a universal outcome of mass extinction. Antell et al. (Reference Antell, Kiessling, Aberhan and Saupe2020), for example, found little general evidence that marine bivalves and brachiopods expanded their geographic ranges following intervals of net biodiversity loss, suggesting that ecological release may be difficult to detect at broad temporal and spatial scales and may be outweighed by noncompetitive factors. Our results therefore do not simply invoke ecological release as a general expectation but identify a more specific phase-dependent signal in which survivor occupancy increased most clearly across the Hirnantian/Rhuddanian transition.
Conclusions
Our study provides a rigorous test of biogeographic homogenization during the LOME using standardized diversity partitioning and network topology. The LOME caused severe taxonomic losses and fundamentally restructured global marine biogeography, with strongly phase-dependent patterns. Signals of homogenization and reduced faunal turnover were weak across the Katian/Hirnantian transition but became pronounced across the Hirnantian/Rhuddanian boundary, coincident with a shift from glaciation-driven cooling to widespread anoxia. Integrating diversity partitioning, network topology, and occupancy analyses indicates that this later phase was associated with reduced faunal differentiation and a transient phase of biogeographic homogenization dominated by cosmopolitan survivors. This “crisis cosmopolitanism” reflects a large-scale ecological response to environmental collapse and the breakdown of geographic and ecological barriers.
Although alpha and gamma diversity recovered relatively quickly in the early Silurian, network modularity and beta diversity recovered more slowly, remaining suppressed through the Rhuddanian and Aeronian and only approaching, but not fully returning to, pre-extinction levels by the Telychian (several million years after the Hirnantian extinction). These results highlight that postextinction recovery operates hierarchically: taxonomic diversity rebounds before ecological and spatial complexity.
The LOME exemplifies how extinction-driven reorganization extends beyond taxonomic loss to encompass the spatial and structural dimensions of biodiversity. Comparable patterns of temporary cosmopolitan expansion followed by renewed provincial differentiation are observed in later crises, including the end-Permian mass extinction and other Phanerozoic events (Button et al. Reference Button, Lloyd, Ezcurra and Butler2017), as well as in Quaternary glacial–interglacial transitions (Yasuhara and Danovaro Reference Yasuhara and Danovaro2016; Yasuhara et al. Reference Yasuhara, Hunt, Dowsett, Robinson and Stoll2020). Such parallels suggest that crisis cosmopolitanism is a recurring macroecological response to global stress: an emergent property of metacommunity dynamics under environmental upheaval.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (42372039), Ministry of Science and Technology of China (2023FY100901), and China Scholarship Council (File No. 202504910072). We thank Bing Huang (Nanjing Institute of Geology and Palaeontology, Chinese Academy of Sciences) and Ádám T. Kocsis (Friedrich-Alexander-Universität Erlangen-Nürnberg) for helpful discussions and assistance with the network analyses. We also thank Simon Darroch and Erin Saupe for their editorial handling, and the anonymous reviewers for their constructive comments and suggestions that helped improve the manuscript. This is a contribution to IGCP project 700 (“Palaeozoic Carbonate Build-Ups in Southeast Asia”) and project 735 (“Rocks and the Rise of Ordovician Life”).
Author Contribution
S.L. curated the data, performed the analyses, and drafted the manuscript. Q.L. contributed to study design and interpretation, L.N. to data curation and methodology, and W.K. to methodology and supervision. All authors contributed to revising the manuscript and reviewed and approved the final manuscript.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
All analyses were conducted in R v. 4.3.1. The data supporting this study are archived in the Dryad Digital Repository: https://doi.org/10.5061/dryad.rv15dv4q9. The R scripts are archived in Zenodo: https://doi.org/10.5281/zenodo.22070262. Supplementary Figs. S1–S7 are archived in Zenodo:https://doi.org/10.5281/zenodo.22070264.
Open access
