Understanding the complex patterns of latitudinal diversity gradients (LDGs) in deep time has been hampered by the absence of long-term records of LDGs through multiple climatic changes. We used records of marine invertebrate fossils to generate LDGs in each age bin from the Carboniferous icehouse to the Triassic greenhouse climates. We evaluated LDGs by subsampling, calculated evolutionary rates for different latitudinal zones, and assessed the modularity of the fossil data within latitudinal zones using the Louvain algorithm. Our results suggest that the LDG peaks may be shaped by multiple factors rather than alternating icehouse and greenhouse climates, although icehouse climates usually restrict diversity at high latitudes. In nearly all clades, peak diversity shifted northward during the late Carboniferous and early Permian, reflecting the northward drift of plates and increased habitat area. Changes in the steepness of the LDG was most pronounced at low latitudes and during biotic crises and recovery, while icehouse to greenhouse transitions created more deviation at high latitudes. Our results show a strong historical influence from previous LDG patterns in LDG dynamics, but one that was interrupted by upheavals in composition after dramatic environmental changes.

Latitudinal gradients in species diversity between equatorial and higher latitudinal zones have long been an important focus of ecological and biogeographic research in both terrestrial and marine systems (Willig et al., 2003; Mittelbach et al., 2007; Mannion et al., 2014). The modern latitudinal diversity gradient (LDG) generally decreases from equatorial to polar areas. However, when fossil data are corrected for sampling bias, times where an LDG is absent, as well as extra-tropical diversity peaks, have been identified in deep time (Mannion et al., 2014). Past research has emphasized the importance of energy (rather than simply climate) as a key factor in the control of LDGs, but ecologists and paleontologists have differed over the importance of historical legacy versus current conditions in the generation of LDGs. One way to resolve these differing views over the proximate causes is to assess LDG changes over long intervals with contrasting climates (Jablonski et al., 2017).

During the Carboniferous to Triassic Periods (~150 m.y.), global climate shifted from long-term icehouse (late Paleozoic ice age [LPIA], >60 m.y. duration) to greenhouse conditions (Montañez and Poulsen, 2013), and this has been regarded as the best analog to the modern transition from icehouse to the current global warming (Penn et al., 2018). This ~150 m.y. interval provides an excellent deep-time system for studying cause-effect relationships between paleoclimatic changes and LDG dynamics. This interval includes several critical bioevents associated with drastic environmental changes (McGhee et al., 2012) and post-crisis rebounds, which generated complex shifts in the taxonomic composition of marine communities. Together, these major macroevolutionary events molded LDG patterns.

A late Permian tropical LDG peak and LDG peaks in temperate regions from the Triassic to Cretaceous have been previously identified (Mannion et al., 2014). However, these studies in that review did not use long-term continuous data. Brachiopod LDGs highlighted the effects of invasion and extirpation on the succession of LDG patterns (Qiao and Shen, 2014; Powell et al., 2015), but whether these results apply to other clades during this interval has not been investigated. Studies across the Permian-Triassic boundary showed a flat LDG after the end-Permian mass extinction (EPME) in ammonoids (Brayard et al., 2006) and other invertebrates (Song et al., 2020), revealing the impact of dramatic environmental changes. Despite many studies of fossil and recent LDG patterns in various environments, the dynamics of these changes are still not well understood (Jablonski et al., 2013) due at least in part to a lack of long-term successive LDG patterns and incomplete and heterogeneous sampling (Jones et al., 2021). In this study, we utilized Carboniferous to Triassic stage-level fossil data to construct LDGs and evolutionary trends, and partitioned the latitudinal zones with network analysis to explore LDG dynamics through an icehouse-to-greenhouse transition.

We assembled data from 13,955 Carboniferous to Triassic genera and subgenera from 303,770 occurrences in the Paleobiology Data-base (https://paleobiodb.org/; see Fig. 1 for detailed numbers of taxa). We refined the data by removing all taxa that were not marine invertebrates with stage-level temporal resolution, and paleo-locations were reconstructed by the GPlates web service (https://gws.gplates.org/). Because of their short duration and similar climatic conditions, the Induan and Olenekian (Early Triassic) were combined into a single age bin. The paleolatitudinal zones were set at 20° intervals from equatorial to polar regions, with the last 30° as one zone to compensate for decreasing area. Generic diversities were calculated at each zone for all age bins (Fig. 2; Figs. S2–S5 in the Supplemental Material1) with corrections for possible sampling biases by multiple methods. Our methods, combined with the raw results, revealed a more comprehensive picture of LDGs than results from raw fossil records. Evolutionary rates and the deviation of diversity (the absolute value of the diversity difference with the adjoining closer-to-equator latitudinal zone divided by mean diversity of a zone at the age bin) were calculated (Fig. 3).

In addition, a bipartite network (Fig. S7) was constructed by treating all genera and latitudinal zones in each age bin as nodes and connecting two nodes when a genus was found in a latitudinal zone in an age bin (Fig. 4A), to detect the relationship among LDGs and their diversity composition. These connections through common genera reflect historic succession and faunal similarity of these latitudinal zones. The Louvain algorithm (Blondel et al., 2008), which partitions networks to optimize the modularity of a network (in this case, to identify latitudinal zones with similar generic composition), and the jackknifing resampling technique were used to recognize latitudinal-zone communities and test their robustness (we use communities here in the network sense, not as paleoecological communities). The jackknifing method resampled 90% of edges (a genus in a latitudinal zone in an age bin) and used the subsampled edges and nodes to build a new bipartite net, then reclustered the net by the Louvain algorithm, repeating this process 1000 times. Further details of the methods and sampling issues are provided in the Supplemental Material.

Unimodal LDG patterns have been identified during both icehouse and greenhouse climates, while flattened or multimodal patterns appeared when climate and energy changed or bioevents occurred. Expected patterns include a tropical LDG peak during icehouse climatic regimes with flattened gradients or a temperate diversity peak during greenhouse climates (Mannion et al., 2014; Meseguer and Condamine, 2020). However, evidence supporting this rule (Krug and Patzkowsky, 2007; Benson and Upchurch, 2012) was mostly from LDG studies of restricted geologic intervals rather than studies of long-term secular patterns. Both the raw and subsampled data (Fig. 2; Figs. S2–S5) show an expansion of high-diversity zones but no apparent movement of the LDG peak beyond the tropics through the transition from the LPIA to a warmer climate regime. Both a flat LDG and multiple LDG peaks occurred under both icehouse and greenhouse climates (Fig. 2; Figs. S2–S4; see a description in the Supplemental Material). Similar patterns have recently been reported for other icehouse and greenhouse climate-transition periods (i.e., a unimodal LDG appearing before global cooling at the Eocene-Oligocene boundary [Crame, 2020] and multiple late Ordovician peaks [Zacaï et al., 2021]). However, the Cenozoic greenhouse and icehouse climatic alternation seems to show a consistent picture of tropical diversity decline and midlatitude diversity (Kiessling et al., 2012; Boag et al., 2021). This might result from the temperature-sensitive clades studied (coral and mollusk) and relatively stable plate positions. From our analysis, we conclude that although energy and area set limits to LDGs, the mode and the peak LDG position may be shaped by multiple factors instead of alternations between icehouse or greenhouse climate as claimed by Mannion et al. (2014) and Meseguer and Condamine (2020).

Although a unimodal, equatorial LDG was commonly found, the LDG center shifted northward during the Carboniferous to Triassic, as evident in the LDG heat map (Fig. 2; Figs. S2–S5) and for nearly all clades (Fig. S6). This phenomenon had previously been identified in specific clades (e.g., Fitzgerald and Carlson, 2006; Powell, 2009) and for Early Triassic marine genera (Naimark and Markov, 2011). Our results indicate a unimodal-bimodal-unimodal alternation in LDG during the early Permian (Fig. 2; Fig. S3), after which the LDG center migrated from the Southern to the Northern Hemisphere. In contrast to the explanations of Naimark and Markov (2011) that diversity peaks appeared near the equator during icehouse intervals and in mid-latitudes during the greenhouse interval, the northward shift in LDG is better explained by tectonic reconfiguration and associated geographic shifts in habitat area. Plate motion and changes in sea level over time modulated habitat area and thus carrying capacity. With continental defragmentation from the late Carboniferous to Late Triassic (Fig. 3G) the supercontinent Pangea and many plates of the Cimmerian belts within the Paleotethys Ocean drifted northward with large habitable shelf areas. Previous quantitative analyses based on maps of reconstructed paleogeography revealed that shelf area in the Northern Hemisphere was more than tripled in extent, especially in middle and high latitudes (Powell, 2009; Scotese, 2021).

The LPIA greatly restricted diversity, especially at high latitudes, and produced a steeper gradient at low latitudes. During the LPIA, high-latitude diversities were generally low (Fig. 2), and this strong climatic effect was also reflected in the LDGs of specific clades (Fig. S6), especially for temperature-sensitive animals (Boag et al., 2021). Clades with broader thermal limits, like Gastropoda, show consistent patterns, while clades with more restricted thermal limits, like Anthozoa and Foraminifera (Song et al., 2014), experienced greater loss at high latitudes during the LPIA and at all latitudes during the EPME. However, the diversity deviation among latitudinal zones (Fig. 3F) showed that the greatest changes in steepness usually occurred at low latitudes and especially after the major mid-Carboniferous and EPME bioevents. In addition, at the end of the LPIA, diversity deviations were larger at high latitudes, reflecting the recovery of high-latitude faunas. From these patterns, we infer a strong role for biotic crises and recoveries in modulating changes in the slope of the LDG, with stronger deviations at lower latitudes.

The tropics have been seen as both a cradle and a museum of diversity (Fischer, 1960; Jablonski, 1993), with both high origination and holdover rates allowing the accumulation of diversity. We did find higher origination and holdover rates at lower latitudes through most of the study interval (Figs. 3A and 3E), but this did not always generate a unimodal LDG. Rather, our results emphasize the importance of climate stability. Substantial changes to LDGs and changes in generic composition coincide with climate shifts at the onset and end of the LPIA and the EPME (Fig. 4B; also see the description in the Supplemental Material). We conclude that the persistence of taxa and higher rates of origination in the tropics did not dominate LDG dynamics in the face of drastic changes in environment or habitat area. A flat LDG appeared in the aftermath of a mass extinction caused by rapid warming and a persistent lethally hot climate (Sun et al., 2012), reflecting severely reduced diversity in tropical areas. As shelf area expanded in the Northern Hemisphere temperate regions during the Triassic, an LDG peak was formed.

However, climate stability seems to be crucial in maintaining this historic succession. For a historic LDG to continue to influence biodiversity requires that many taxa survive to the next age bin. As previously discussed, the latitudinal zones across all age bins could be divided into several communities based on the composition of shared taxa. These communities generally had latitudinal continuity (especially under greenhouse climate in the Triassic) but were truncated by global climatic changes during the LPIA and the EPME hypothermal events. Because the Louvain algorithm optimizes partitions within communities, with fewer connections between communities (Fig. 4A), clear boundaries at the onset and end of LPIA and the EPME revealed the great change in generic composition at all latitudes (Fig. 4B; see the Supplemental Material). This analysis also distinguished the middle to low latitudes from high latitudes between the LPIA and the EPME. The similarity of Triassic communities revealed the striking impact of the EPME and the generally harsh greenhouse climate of the Triassic as well as the preponderance of cephalopods in many assemblages. These community partitions also reveal that environmental changes over geologic time played a more important role than geographic differences in latitudinal taxonomic composition. In addition to the temporal community boundaries, the latitudinal boundaries before the Triassic implied a more latitudinally independent evolutionary model, but widespread taxa prevailed in the Early Triassic (e.g., ammonoids and bivalves).

Our study illustrates the value of a long-term, deep-time record for establishing the historical dimensions of LDGs. But both ecological and paleoecological studies have shown that modern LDGs can shift on decadal time scales due to global warming (Yasuhara et al., 2020; Chaudhary et al., 2021), possibly reflecting rapid changes in temperature and oxygen content (Yasuhara and Deutsch, 2022). The stage-level fossil data used in our analysis would miss such dynamics. Examining more rapid changes in LDGs in deep time requires data sets with much higher temporal resolution.

This work is supported by the Second Tibetan Plateau Scientific Expedition and Research (grant 2019QZKK0706), the Frontiers Science Center for Critical Earth Material Cycling Fund (grant DLTD2102), and the Strategic Priority Research Programs of the Chinese Academy of Sciences (grant XDB26000000). We thank Moriaki Yasuhara, Michael Benton, and an anonymous reviewer for their comments on a previous draft of this manuscript.

1Supplemental Material. Supplemental information about methods, results, raw data, and r scripts. Please visit https://doi.org/10.1130/GEOL.S.20164271 to access the supplemental material, and contact editing@geosociety.org with any questions.
Gold Open Access: This paper is published under the terms of the CC-BY license.