Quantifying ecological impacts of mass extinctions with network analysis of fossil communities.

Mass extinctions documented by the fossil record provide critical benchmarks for assessing changes through time in biodiversity and ecology. Efforts to compare biotic crises of the past and present, however, encounter difficulty because taxonomic and ecological changes are decoupled, and although various metrics exist for describing taxonomic turnover, no methods have yet been proposed to quantify the ecological impacts of extinction events. To address this issue, we apply a network-based approach to exploring the evolution of marine animal communities over the Phanerozoic Eon. Network analysis of fossil co-occurrence data enables us to identify nonrandom associations of interrelated paleocommunities. These associations, or evolutionary paleocommunities, dominated total diversity during successive intervals of relative community stasis. Community turnover occurred largely during mass extinctions and radiations, when ecological reorganization resulted in the decline of one association and the rise of another. Altogether, we identify five evolutionary paleocommunities at the generic and familial levels in addition to three ordinal associations that correspond to Sepkoski's Cambrian, Paleozoic, and Modern evolutionary faunas. In this context, we quantify magnitudes of ecological change by measuring shifts in the representation of evolutionary paleocommunities over geologic time. Our work shows that the Great Ordovician Biodiversification Event had the largest effect on ecology, followed in descending order by the Permian-Triassic, Cretaceous-Paleogene, Devonian, and Triassic-Jurassic mass extinctions. Despite its taxonomic severity, the Ordovician extinction did not strongly affect co-occurrences of taxa, affirming its limited ecological impact. Network paleoecology offers promising approaches to exploring ecological consequences of extinctions and radiations.

Mass extinctions and radiations entail significant changes in biodiversity and community ecology (1–3). Paleontologists have traditionally focused on measuring taxonomic richness, a proxy for biodiversity, as a means of identifying critical transitions in the geological record and evaluating their significance for biosphere evolution (4, 5). However, this approach does not take into account ecological changes, which are generally decoupled from taxonomic turnover and cannot be counted like taxa (6–9). Major paleoecological changes include appearances/disappearances of communities and ecosystems, as exemplified by reefs; increases/decreases in tiering complexity; and contractions/expansions of relative clade abundances and their geographic/environmental distributions. To compare critical transitions in terms of ecological severity, some studies have adopted qualitative approaches, treating paleoecological changes as ranked phenomena and ascribing extinctions/radiations to corresponding levels (2, 3, 6). Ecological changes, however, do not constitute a hierarchy (6), and qualitative approaches can lead to inconsistent results. The lack of quantitative metrics for evaluating ecological disruptions in Earth history hinders our ability to compare biotic crises of the past reliably and understand their implications for twenty-first-century global change.

Here, we apply a network-based approach to study community paleoecology and quantify impacts of extinctions and radiations. Network paleoecology is the study of long-term changes in biosphere structure, function, and evolution using network analysis, computational models, and empirical data from the fossil record. Network analysis of fossil co-occurrence data allows us to detect groups of paleocommunities and paleocommunity types that appeared and disappeared en masse across critical transitions (). For this work, we recognize paleocommunities as recurrent and distinguishable associations of taxa and paleocommunity types as aggregates of communities with similar, but not identical, taxonomic makeups and environmental preferences (10). Our approach builds on the observations that (i) communities generally exhibit stasis in structure and taxonomic makeup over geologic time and (ii) extinctions and radiations lead to turnover of groups of interrelated communities (11, 12). Such groups correspond to ecologic evolutionary units (EEUs), or long (30–140 Ma) intervals of relative paleocommunity stasis separated by relatively short (3–5 Ma) episodes of ecological reorganization preserved in the geologic record (11–13). During the periods of stasis, few communities appeared or disappeared, extinctions and radiations did not significantly affect the relative abundances or interactions of clades, and species only rarely moved beyond their environments of origination (13). The Phanerozoic record contains between 9 and 12 EEUs, depending on how they are defined (11–13). These EEUs may be interpreted as ecological subdivisions of evolutionary faunas (1, 13). However, whereas the evolutionary faunas were identified by quantifying fossil diversity patterns (1), EEUs were defined qualitatively (12, 13). Hence, opinions differ concerning their number, boundaries, and definitions (12, 13). To overcome these issues, our approach provides a means of analyzing fossil data, identifying EEU-like groups of interrelated paleocommunities, and quantifying pulses of change in their representation over geologic time.

Network Analysis

Network theory deals with the study of complex systems of interconnected entities (14–16). A typical network consists of nodes as well as the interactions (links) among them (). For this work, nodes are taxa, and two taxa are linked if they co-occur at any fossil collection sampling point of Phanerozoic age (541–0 Ma). A network-based approach has several advantages over other methodologies for paleocommunity analysis (17, 18). First, it allows us to apply and compare multiple methods of detecting communities (19), which are, by definition, networks (i.e., groups of actually and potentially interacting entities). Second, it supports an assortment of metrics describing local (node specific) and global (whole network) properties of networks and the factors influencing their structures (20). Lastly, network theory accommodates multimodal (pseudohierarchical and nested) data structures with macro-, meso-, and microlevel communities of varying significance (21). In this context, we downloaded data on co-occurrences of marine animal fossils in the fossil collections of the Paleobiology Database (PBDB) and analyzed five related networks (Fig. 1 and and Table S1). Each fossil collection in the PBDB corresponds to a sampling point with a unique geographic and stratigraphic location (). Our links, consequently, represent the smallest co-occurrence units in the PBDB, and signify actual cases of specimen sampling and reporting. Accordingly, our network linkages contain patterns related to geologic age, geographic location, and preservational pathway, all of which dictate co-occurrence potential. Because taxa that lived at different times in the Phanerozoic Eon cannot be linked, the most ancient and least ancient taxa cluster on opposite sides of a given network. This polarity signifies an implicit time axis for assessing ecological change.

Fig. 1.

Graphs of nonweighted networks partitioned with the Louvain algorithm. (A) The o–n network of marine animal orders (n = 621), which fall into three main modules, resembling the Cambrian, Paleozoic, and Modern evolutionary faunas (1). (B) The f–n network of marine animal families (n = 3,729), which fall into five modules, similar to the main clusters in the corresponding generic-level (g1–n and g2–n) networks. Values in parentheses indicate the modules’ proportions of their respective networks. See figures in for additional graphs ().

The five networks differ with regard to their taxonomic level and community specificity (). Three of the networks—the ordinal (o) and familial (f) networks, as well as one of the generic (g1) networks—differ in taxonomic rank, but in each case, represent >95% of valid marine animal taxa with links at each rank ( and Dataset S1). Communities can share all of the same taxa but differ with regard to subtaxa (i.e., one community described at high rank may correspond to two at a lower level), so parallel analyses of the o, f, and g1 networks enable us to assess rank dependency of results. The remaining (g2 and g3) networks represent g1 subsets, and consist entirely of genera recorded within reef facies in the PBDB, specifically, all reef-occurring (g2) and reef-building (g3: bryozoan, cnidaria, rudist, and sponge) genera. Taxa without connections to nodes in these five networks are generally rare (i.e., paleogeographically, taphonomically, and/or sampling limited) taxa, and were omitted from the analyses. Together, the first three networks enable us to explore ecological changes that affected the marine realm as a whole, whereas the other two support focused analyses of reefs, which are biodiversity hotspots that are sensitive to environmental perturbation (6).

To each network, we applied a variety of algorithms for detecting communities. These algorithms facilitate discovery of groups of nodes, which appear as clusters (modules) in network graphs (Fig. 1 and ). In essence, these various algorithms differ with regard to the definition of community and, for a given dataset, may identify somewhat different groups (19). The algorithms also vary in computational complexity, and network size affects analysis speed. Hence, in selecting an algorithm, one must consider both the communities of interest and the global properties of the network. We primarily focus on the Louvain algorithm (22), which is generally recommended for partitioning large networks (n > 6,000 nodes) and those with potentially high mixing parameters (), i.e., great numbers of links between nodes belonging to different communities (19). We believe this algorithm is appropriate, given the properties of our networks.

The Louvain algorithm partitions networks by optimizing for modularity (23), a “global” property related to the strength of division of a network into groups of nodes. For a given set of clusters within a network, modularity is the fraction of links that connect nodes of the same communities minus the corresponding fraction expected in an equivalent network with a random distribution of connections (). In general, modularity scores range between 0 and 1, and values >0.3 are good indicators of community structure (23). Through optimization of this value, the Louvain algorithm attempts to find the best community structure, defined as the set of modules for which the number of links within the communities is maximized () and the number between them is minimized. Weight attributes (e.g., connection strength) also influence Louvain modularity calculation and optimization (22), and can affect network partitioning results. For this reason, we applied the algorithm to nonweighted and weighted versions of each network. In the nonweighted versions (o–n, f–n, g1–n, g2–n, and g3–n), links are binary (present or absent) and equally weighted. Conversely, in weighted versions (o–w, f–w, g1–w, g2–w, and g3–w), link weights are co-occurrence counts from the PBDB. Following this approach, in both weighted and nonweighted networks, modules are separated by regions of relatively sparse connections. These regions are consequences of community turnover, which impacted the range of interactions among taxa in space and time. Partitions in our weighted networks additionally reflect patterns in the frequency and sampling of co-occurring taxa.

Phanerozoic Fossil Networks

The o, f, and g1 networks exhibit variation in the following global properties ( and Tables S4–S6): diameter (maximum degree of separation), edge density (the ratio between the numbers of actual connections and possible links), transitivity (a coefficient measuring the probability that neighbors of a node are connected), and modularity. Edge density and transitivity are greatest in the o network, and diameter and modularity are greatest in the g1 network (). Thus, the g1 network is large and loosely connected (), the o network is small and dense (Fig. 1 and ), and the f network exhibits intermediate properties (Fig. 1 and ). Because most taxa consist of multiple subtaxa, which contribute co-occurrences up the taxonomic hierarchy, taxa will naturally be less numerous and share more connections at higher than lower ranks. Beyond these metrics, the three networks also differ in terms of assortativity coefficients that measure homophily (), or the tendency of nodes to associate with others possessing similar properties (e.g., geologic age, systematic affinities, or geographic ranges). Taxa in the g1 network exhibit the strongest tendency to associate with others of similar age (first absolute age datum and mean age of occurrence), taxonomic placement (phylum or class), and paleogeographic location, followed by those in the f and o networks, respectively. On average, orders have longer stratigraphic ranges (135.5 Ma) and have higher probabilities of co-occurring than families (55.6 Ma) or genera (25.1 Ma). Therefore, orders are most commonly linked to dissimilar taxa and have the greatest potential to diverge from connected taxa with respect to first and mean ages. Regardless, age represents the best predictor of association for taxa in all networks ( and Table S6), manifesting in network polarity (i.e., an implicit time axis).

For each network, the Louvain algorithm facilitated detection of modules with modularity scores >0.3 (Fig. 1 and and Table S4). In comparison with the other suitable methods, the Louvain algorithm consistently returned the highest modularity scores, affirming that it performed best at identifying nonrandom associations of taxa. It also typically returned a low number of modules, making it one of the most conservative approaches. Although some networks (e.g., g1–n and g1–w) contain numerous small clusters (Dataset S1), the majority of taxa in each network are divided among a few (3–8) large clusters (Fig. 2 and ), each containing hundreds to thousands of nodes. As a result, we can interpret the large Louvain modules as relatively inclusive communities containing smaller subcommunities (). For the most part, each large module includes taxa from a variety of marine environments ( and S11, Table S7, and Dataset S1) and consecutive periods of Earth history (Fig. 2). Taking these observations into consideration, we interpret the large modules as macrolevel units encompassing microlevel paleocommunities and mesolevel paleocommunity types (10). Given that geologic age exerts a strong influence on data structure, we surmise that the boundaries between modules demarcate turnover (loss and emergence) of interrelated paleocommunities over time. We propose the term evolutionary paleocommunity for the macrolevel unit, a designation recognizing that associations of interrelated paleocommunities emerge as a result of evolutionary patterns (11–13). Naturally, the macrolevel units differ from each other in terms of taxonomic makeup (Fig. 3), and although they broadly overlap in stratigraphy, they evidently peaked in generic richness and dominated total diversity at different times (Fig. 2). These observations affirm that the evolutionary paleocommunities, like EEUs, succeeded each other through the Phanerozoic, each supplanting its predecessor during a biotic crisis involving ecological reorganization (Fig. 2).

Fig. 2.

Plots show diversity versus geologic age for modules in the nonweighted networks. The value of each bin equals the total number of boundary crossing genera, i.e., total diversity minus singletons (4). (A) The o–n network of marine animal orders (Fig. 1 and ). (B) The f–n network of marine animal families (Fig. 1 and ). (C) The g1–n network of marine animal genera (n= 26,525; ). (D) The g2–n network of reef-occurring animal genera (n = 7,621; ). (E) The g3–n network of reef-building animal genera (n = 2,793; ).

Fig. 3.

Taxonomic loading of modules in the nonweighted network (g1–n) of marine animal genera (Fig. 2 and ). Plot shows proportions of modules for various nonnested clades (classes, subclasses, and infraclasses). These clades account for 95% of the taxa in the network. Detailed breakdowns of each module in this study are presented as interactive sunburst diagrams, which can be found online: https://deeptime.tw.rpi.edu/viz/SunBurst_Fossils_ByCommDetection/Sunburst_f_all.html.

The o–n network contains five modules (Fig. 1 and ) with distinct clusters dominating total diversity in the Cambrian, Ordovician–Permian, and Mesozoic–Cenozoic intervals (Fig. 2). Two other modules, which are dominated by specific groups (Triassic chordates and Paleogene “worms”), are probably consequences of taxonomic homophily induced by sampling and reporting biases. Although these other modules account for ∼12% of orders, their contributions to total diversity at the generic level are relatively minor for all geologic stages, and they never dominate generic diversity (Fig. 2). Our analyses, therefore, demarcate three evolutionary paleocommunities at the ordinal level. These associations broadly correspond to Cambrian, Paleozoic, and Modern evolutionary faunas (1), indicating that ecological restructuring followed major shifts in the balance of animal diversity over time.

Because community structure can vary with taxonomic rank, the ordinal (o) network differs from the corresponding familial (f) and generic (g1) networks in this study (). The f–n network contains five modules (Fig. 1 and ), which peaked in diversity in the Cambrian, Ordovician–Devonian, Carboniferous–Permian, Mesozoic, and Cenozoic intervals (Fig. 2). Evidently, the Devonian and Cretaceous–Paleogene (K–Pg) crises resulted in division of families (but not orders) into distinct associations. Partitioning of the g1–n and g2–n networks returned comparable clusters (Fig. 2 and and S7), demarcating five main modules (>97% of genera) in addition to various small modules that never dominate diversity (). Taken as a whole, our analyses of the f, g1, and g2 networks support the existence of five evolutionary paleocommunities of animal life at the family and genus levels. This result applies both to the entire marine realm (Figs. 1 and 2 and and S6) and to reefs (Fig. 2 and ), affirming that the history of reefs largely parallels that of the fossil record as a whole (24).

The g3–n network contains nine modules (), including six modules that account for the majority (∼98%) of taxa as well as three (“other”) low diversity modules (Fig. 2). Four of the six main modules resemble analogous units in the f–n, g1–n, and g2–n networks (Figs. 1 and 2 and , S6, and S7). In the Cenozoic, two (g3–CzA and g3–CzB) modules () represent roughly comparable numbers of taxa (Fig. 2), with ranges that broadly overlap through stratigraphy. These modules represent consequences of taxonomic homophily, as the g3–CzA and g3–CzB modules are largely comprised of anthozoan corals and bryozoans, respectively. In view of these biases, these two modules are perhaps best interpreted as a single Cz-equivalent evolutionary paleocommunity (Fig. 1 and , S6, and S7).

Inclusion of weights did little to change the results. Partitioning of the weighted o–w, f–w, g1–w, and g2–w networks yielded modules comparable to those in the nonweighted networks, except that in each case (), the algorithm divided Mesozoic and Cenozoic taxa among additional small modules (). Each of these additional modules is dominated by one or several clades (e.g., Anthozoa, Bryozoa, or Chordata), suggesting that they reflect biases in sampling of taxa in those clades. Despite these biases, the analyses corroborate identification of five evolutionary paleocommunities at the family and genus levels as well as the three macrolevel associations of orders (). The g3–w network differs most substantially from its nonweighted counterpart (). This network consists of 10 modules, including 8 representing the majority (∼93%) of taxa in addition to two (other) Mesozoic sponge- and rudist-dominated clusters (). Four of the eight main modules resemble clusters detected in the corresponding g3–n network. However, unlike its counterpart, the g3–w network contains g3–w–O, g3–w–SD, g3–w–Tr, and g3–w–JK modules, which dominated diversity in the Ordovician, Silurian–Devonian, Triassic, and Jurassic–Cretaceous intervals, respectively. Thus, the g3–w network provides evidence that the Ordovician and Triassic–Jurassic mass extinctions strongly affected sessile benthic reef-building animals.

Implications of Network Paleoecology for Mass Extinctions

Our analyses lay the groundwork for several approaches to evaluating ecological impacts of critical transitions. First, the results permit qualitative assessment of their impacts at various taxonomic and ecological levels. For example, our results demonstrate that the Permian–Triassic (P–Tr) mass extinction affected communities at the ordinal, familial, and generic levels, as opposed to the Devonian and K–Pg crises, which significantly affected familial and generic but not ordinal compositions. Our results also show that although marine communities did not change radically during the Ordovician and Triassic–Jurassic (T–J) mass extinctions, ecological restructuring during these events affected reefs more than other marine ecosystems (25). Second, network analysis supports an approach to quantifying magnitudes of ecological change. In essence, we calculate changes in the relative loadings of evolutionary paleocommunities through time (Fig. 4). This approach recognizes that the emergence/disappearance of communities parallels the origination/extinction of taxa and that the ecological severity of an event is related to its selectivity, i.e., the degree that taxonomic turnover diminishes certain types of communities in favor of others (11–13). Our metric of ecological change is “total swing,” which we determined for each geologic stage of the Phanerozoic above the Fortunian (Fig. 4). This metric is calculated from shifts over geologic time in the relative (percentage) contributions of network modules (i.e., evolutionary paleocommunities) to total diversity. For a given geologic age, in which one evolutionary paleocommunity supplants another, the total swing is approximately equal to the relative increase in representation of the former as well as the absolute decrease in representation of the latter (Fig. 4).

Fig. 4.

Ecological severities of biotic crises and critical transitions. Plot shows the total swing in diversity of each Phanerozoic (post-Fortunian) geologic stage for the nonweighted (o–n, f–n, and g1–n) networks. A swing is a relative change in representation (i.e., growth of one module and reduction of another) across an interval, and is expressed in terms of percent total diversity minus singletons (data downloaded from the PBDB using the contain method of time binning). To calculate total swing for a geologic stage, the absolute difference in percent total diversity between the upper and lower boundaries was determined for each module. The sum of these values was then determined and divided by 2. Results did not significantly vary with binning approach ().

By considering total swing over time, we can assess the record of ecological change (Fig. 4 and ). In the Cambrian, ecological changes appear concentrated in the stage 3 and Drumian–Guzhangian intervals (), which encompass the “Botomian” and “Dresbachian” extinctions (26), respectively. However, these patterns may at least in part reflect rock record biases (26) and/or artifacts of taphonomic windows (27), stratigraphic gaps (28), and endemic faunas (29). As expected, and in agreement with previous studies (2, 3, 8, 9), our data show that the major biotic crises of the Phanerozoic Eon (Figs. 2 and 4)—the big five mass extinctions (5)—involved significant ecological reorganization. In each case, ecological change proceeded over an extended period of time, varying in duration (up to tens of millions of years) with the nature of biodiversity loss (elevated extinction versus reduced origination) and the rate and timing of recovery. The Devonian and middle Permian–Triassic stand out as intervals of sustained ecological change that likely encompass multiple pulses of taxonomic loss and recovery (3, 17, 30–32).

The P–Tr mass extinction involved the greatest disruptions to communities, and the Ordovician event entailed the least significant changes in ecology. These results do not vary with taxonomic rank (Fig. 4), and are consistent with interpretations of qualitative data (2, 3, 8, 9). Our total swing metric, however, supports a ranking of the K–Pg, Tr–J, and Devonian crises that differs from previous treatments ( and S9), showing that the K–Pg event entailed the greatest swings at the generic and familial levels, and that the Tr–J event involved the most severe changes at the ordinal level (Fig. 4). Interestingly, comparing all critical transitions in the Phanerozoic record, the greatest ecological changes at the family and genus levels did not result from a mass extinction, but rather from the Great Ordovician Biodiversification Event (GOBE), a major radiation of skeletonized invertebrates (33), possibly driven by ocean cooling and/or oxygenation (34, 35). Like the P–Tr mass extinction, the GOBE left a strong imprint on ecological organization, as evidenced by the turnover of evolutionary faunas (1) and ordinal-level evolutionary paleocommunities (Figs. 1 and 2).

The Tr–J and Ordovician extinctions differ from the other major biotic crises of the Phanerozoic. Both extinctions entailed ecological reorganization of reef communities, including significant changes in the frequencies and interactions of reef-building animals ( and S12), but neither event left a strong imprint in the comprehensive (o, f, and g1) networks in this study or entailed turnover of the evolutionary paleocommunities represented therein. Overall, these results corroborate qualitative studies ( and S9) that rank the severity of Ordovician extinction well below other events in Earth history (2, 3, 8). Indeed, our total swing metric suggests that environmental and biodiversity crises in the Pennsylvanian, mid-Permian, late Triassic, mid-Cretaceous, and Eocene may have induced ecological changes of greater impact than the Ordovician extinction ().

During middle to late Devonian times, a severe reduction in the origination of taxa caused a mass depletion of biodiversity (32, 36). The high swing values of the Givetian, Frasnian, and Famennian stages (Fig. 4) suggest that the significance of this biotic crisis, the Devonian mass extinction, has been underestimated (3, 8, 9). Our results indicate that, in terms of ecological severity, the Devonian ecosystem changes as a whole rank second or third behind the P–Tr mass extinction but on par with the K–Pg event (Fig. 4 and and Tables S8 and S9). Although the P–Tr and K–Pg extinctions were associated with relatively rapid (∼10 Ma) episodes of ecological change, the Devonian crisis entailed protracted replacement of one evolutionary paleocommunity by another over 50 My (Figs. 2 and 4 and and S13). In this light, the Devonian does not conform to the EEU model (11–13), which implies that ecological reorganization typically occurs during short intervals of community turnover. Instead, ecological reorganization occurred throughout the time of mass biodiversity depletion, amounting to secular variation in marine animal communities.

Conclusions

Methods rooted in network theory provide powerful tools for describing and evaluating long-term ecological change, specifically the rise and fall of interrelated communities over time. Identification of evolutionary paleocommunities supports qualitative and quantitative approaches to assessing the severities of mass extinctions. These approaches have great potential, as they utilize fossil co-occurrence data, which are common and easy to compile. Given that network analysis has received limited attention in paleobiology, much work remains to be done. Ultimately, network paleoecology may help to bridge gaps between ancient and modern networks of life, illuminate the threats associated with continuing species losses, and contribute to development of conservation strategies for mitigating the present biotic crisis (37).

Methods

Data on fossil occurrences, taxonomy, and diversity over time were accessed from the Paleobiology Database on various dates in spring 2017 (https://paleobiodb.org/#/). Taxonomic data were revised to correct for inconsistently ranked clades, update classification schemes, and fill in empty fields (). In compiling diversity over time data, we assigned occurrences to geologic stages based on the overlap of age ranges using the “contain” (occurrence age ranges are strictly contained within time bins) and “major” (50% or greater overlap) methods. Our estimates are based on occurrences of regular (body) genera, exclude uncertain taxa, and assume extant taxa range to present. Fossil co-occurrence data were accessed using a custom application programming interface (API) written in Python. The API (https://github.com/zhongh/dtdi-api/tree/b42540b6062e390f31f57aa9259214c33f332ccc) outputs adjacency lists and matrixes containing co-occurrence counts. The built-in functionality of the API allows for varying taxa with respect to age, level (rank), and clade. We processed the API outputs to remove form and trace taxa, invalid taxa (i.e., taxa of uncertain rank and those with subordinate names), and taxa lacking connections or occurring in small isolated clusters.

The datasets were analyzed in RStudio using functions in the igraph, GGally, network, ggplot2, and sunburstR packages. Network graphs were generated using the ggnet2 function of ggplot2 and its default parameters, and nodes of equal size were placed without self-loops according to the Fruchterman-Reingold force-directed algorithm. Measures of whole-network properties were computed using functions of the igraph package. The networks were partitioned into mutually exclusive groups (modules/clusters) based on their nonweighted (equally weighted) and weighted (co-occurrence count) link attributes using community-detection algorithms of the igraph package. For each module, the number of occurrences of taxa was determined for a variety of paleoenvironments, based on data in the “environment” field of the PBDB. For simplicity, we combined similar environments into more inclusive paleoenvironmental categories (). Lastly, sunburst diagrams depicting the taxonomic breakdowns of the network modules were generated for analysis using the sunburstR package.