The Fezouata Shale Formation biota is typical for the high latitudes of the Early Ordovician—a quantitative approach

Abstract

The Fezouata Shale Formation has dramatically impacted our understanding of Early Ordovician marine ecosystems before the great Ordovician biodiversification event (GOBE), thanks to the abundance and quality of exceptionally preserved animals within it. Systematic work has noted that the shelly fossil subassemblages of the Fezouata Shale biota are typical of open-marine deposits from the Lower Ordovician, but no studies have tested the quantitative validity of this statement. We extracted 491 occurrences of recalcitrant fossil genera from the Paleobiology Database to reconstruct 31 subassemblages to explore the paleoecology of the Fezouata Shale and other contemporary, high-latitude (66°S–90°S) deposits from the Lower Ordovician (485.4–470 Ma) and test the interpretation that the Fezouata Shale biota is typical for an Ordovician open-marine environment. Sørensen's dissimilarity metrics and Wilcoxon tests indicate that the subassemblages of the Tremadocian-aged lower Fezouata Shale are approximately 20% more heterogenous than the Floian-aged upper Fezouata Shale. Dissimilarity metrics and visualization suggest that while the lower Fezouata and upper Fezouata share faunal components, the two sections have distinct faunas. We find that the faunal composition of the lower Fezouata Shale is comparable with other Tremadocian-aged subassemblages from high latitudes, suggesting that it is typical for an Early Ordovician open-marine environment. We also find differences in faunal composition between Tremadocian- and Floian-aged deposits. Our results corroborate previous field-based and qualitative systematic studies that concluded that the shelly assemblages of the Fezouata Shale are comparable with those of other Lower Ordovician deposits from high latitudes. This establishes the first quantitative baseline for examining the composition and variability within the assemblages of the Fezouata Shale and will be key to future studies attempting to discern the degree to which it can inform our understanding of marine ecosystems just before the start of the GOBE.

Non-technical Summary

The Fezouata Shale Formation in present-day Morocco is a site of exceptional fossil preservation from the Lower Ordovician that provides a unique view of animal life before one of the most important radiation events in Earth's history, the great Ordovician biodiversification event (GOBE). Previous work on the fossil diversity of the Fezouata Shale has suggested that there are faunistic differences between the two major intervals with exceptional preservation and that the overall shelly biota of the Fezouata Shale is comparable to other Lower Ordovician sites that reflect open-marine conditions. In this study, we make the first comprehensive quantitative comparison between the Fezouata Shale Formation and other high-latitude Early Ordovician sites based on their shelly fossil biotas with publicly available fossil occurrence information from the Paleobiology Database. We find that the fossil subassemblages of the stratigraphically older lower Fezouata Shale are more heterogeneous than those of the younger upper Fezouata Shale. The fossil biota preserved in the lower Fezouata Shale is most similar to those found in other high-latitude deposits from the Lower Ordovician. We also find that there are differences in faunal composition between Tremadocian- and Floian-aged deposits. Our work provides the first quantitative support for faunistic differences between the lower and upper Fezouata Shale Formation and indicate that the lower Fezouata Shale conventional fossil biota is typical for the Tremadocian, further contextualizing the ecology of the polar regions before the GOBE as informed by this major site of exceptional fossil preservation.

Introduction

The Fezouata Shale Formation contains one of the most abundant and diverse marine fossil biotas from the Lower Ordovician known to date (Van Roy et al. 2010, 2015a; Lefebvre et al. 2016c). The first exceptionally preserved specimens from this formation were discovered by the fossil collector Mohamed “Ou Said” Ben Moula approximately two decades ago (Van Roy et al. 2015a; Lefebvre et al. 2016b). Today, multiple Fezouata Shale Formation subassemblages that preserve fossils of soft-bodied organisms are known throughout the Draa Valley, an area encompassing a few hundred square kilometers in the Anti-Atlas region of southeastern Morocco (Van Roy et al. 2010, 2015a; Martin et al. 2016a; Saleh et al. 2021b, 2022a,b). Remarkable fossil preservation in the Fezouata Shale Formation is exposed in two biostratigraphically distinct fossiliferous intervals that can be readily distinguished based on the abundance and composition of preserved graptolites. The upper Tremadocian lower Fezouata Shale Formation has been linked to the Araneograptus murrayi and Hunnegraptus copiosus biozones, while the middle Floian upper Fezouata Shale Formation is correlated with the Baltograptus jacksoni biozone (Gutiérrez-Marco and Martin 2016; Lefebvre et al. 2018; Saleh et al. 2018).

Much of the interest in the Fezouata Shale Formation stems from its highly fossiliferous faunal assemblages. More than 160 marine animal genera have been identified from the two exceptional intervals of the Fezouata Shale Formation (Van Roy et al. 2015a; Lefebvre et al. 2016b). The biota contains major animal groups that can be tracked back to the Cambrian radiations (e.g., radiodonts, aglaspidids, marrellomorphs, nektaspids, and sachitids) alongside newer representatives of the Paleozoic evolutionary faunas (e.g., asterozoans, cephalopods, and chelicerates), thus extending the known temporal range for multiple animal groups (Van Roy et al. 2015b; Legg 2016; Ortega-Hernández et al. 2016; Vinther et al. 2017; Hunter and Ortega-Hernández 2021; Pérez-Peris et al. 2021). Given the ongoing discussion around the timing and duration of the evolutionary radiations that occurred during the Middle Ordovician within a variety of animal clades, known collectively as the great Ordovician biodiversification event (GOBE) (Droser 2003; Servais et al. 2010, 2023; Servais and Harper 2018; Stigall 2018; Stigall et al. 2019), the significance of the Fezouata Shale biota goes beyond simply being an Ordovician Konservat-Lagerstätte. Instead, the Fezouata Shale is critical to our understanding of the biodiversity of animal-dominated ecosystems before the initial stages of the GOBE (Servais et al. 2023).

It is necessary to characterize the taphonomy, paleoenvironment, and paleoecology of the Fezouata Shale Formation biota to fully understand the contribution of its fossils to the broader macroecological trends in the early Paleozoic. Critically, our understanding of these aspects of the Fezouata Shale Formation has advanced dramatically in the last decade. Taphonomic investigations suggest that the organisms preserved in the Fezouata Shale biota were buried in situ through storm-associated deposits (Saleh et al. 2020, 2021b) below storm-wave base on an open siliciclastic ramp roughly 70° south of the equator (Saleh et al. 2020, 2021b; Martin et al. 2016a). Data using brachiopods and bivalves suggest that the communities of the Fezouata Shale biota were ephemeral (Saleh et al. 2018). Work focused on trilobites, graptolites, and cephalopods has provided information on the distribution of these clades within the Fezouata Shale (Kröger and Lefebvre 2012; Gutiérrez-Marco and Martin 2016; Martin et al. 2016b). This, alongside publications summarizing faunal lists for the deposit (Van Roy et al. 2015a; Lefebvre et al. 2016b), has laid the fundamental groundwork for our understanding of the paleoecology of the Fezouata Shale. More recent studies have also revealed well-preserved instances of interspecific ecological interactions within the Fezouata Shale biota, including the use of shells as hard substrates for the growth of epibenthic organisms (Nanglu et al. 2023). Although it has been shown that the Fezouata Shale Formation is dominated by echinoderms, arthropods, graptolites, and brachiopods (Kouraiss et al. 2019; Saleh et al. 2022b), there has been little to no community-level work examining temporal or spatial patterns across the entire formation.

A high-resolution analysis of the Fezouata Shale Formation's ecology in its totality, as well as how the distinct stratigraphic assemblages within the formation compare with other Ordovician localities, has never been performed. Initial qualitative assessments of the Fezouata Shale biota have suggested that the shelly communities are broadly comparable to those of other open-marine deposits from the Lower Ordovician and typified by trilobites, brachiopods, cephalopod and bivalve mollusks, pelagic graptolites, and echinoderms (Van Roy et al. 2010, 2015a; Kröger and Lefebvre 2012; Lefebvre et al. 2016c, 2018). A quantitative evaluation of some of the main shelly fossils (e.g., trilobites, and echinoderms) preserved in the central Anti-Atlas suggests differences in the assembly of Tremadocian- and Floian-aged communities (Saleh et al. 2022a). These claims are noteworthy, because, if accurate, they better contextualize the Fezouata Shale biota within the broader ecological changes taking place before the GOBE (Servais et al. 2010; Servais and Harper 2018) and allow this deposit to serve as a key point of comparison for understanding the Cambrian to Ordovician evolutionary transition. However, the fundamental question of whether the Fezouata Shale biota is representative of a typical Ordovician fauna, whether measured by faunal composition, the magnitude of diversity, or ecological structure, remains unanswered.

In this study, we downloaded occurrence data from the Paleobiology Database (PDBD) to quantitatively assess the faunal similarity between the Fezouata Shale relative to other high-latitude (66°S–90°S) marine deposits from the Lower Ordovician (485–470 Ma). This approach allows us to further address two major questions: (1) Do the Tremadocian-aged lower Fezouata Shale Formation and the Floian-aged upper Fezouata Shale Formation intervals differ in terms of their faunal composition? (2) Is the Fezouata Shale Formation biota quantitatively similar to other high-latitude Lower Ordovician fossil assemblages?

Materials and Methods

Tools

We performed data manipulation using Excel and in RStudio 4.1.0 with the packages dplyr, tidyr, readr, reshape2, janitor, and plyr (Wickham 2011, 2020; Wickham et al. 2017, 2023a,b; Firke et al. 2020; R Core Team 2022). We computed ecological metrics in RStudio using the package vegan (Oksanen et al. 2019). We carried out statistics and visualization in RStudio using the packages vegan, stats, dunn.test and funfuns, ggpol, ggplot2, and colorRamps (Wickham 2009; Keitt 2012; Dinno 2017; Oksanen et al. 2019; Tiedemann 2020; R Core Team 2022; Trachsel 2023). We downloaded occurrence data from the PDBD on June 21, 2022 (Supplementary Dataset 1).

Data Selection

Fossil occurrences of genera (485–470 Ma) include information related to age (Tremadocian vs. Floian), paleolatitudes, and geological formations of origin. For the purposes of this study, we only considered metazoan occurrences, so all prokaryotes and protists (e.g., cyanobacteria, radiolarians, and foraminiferans) were excluded. We further restricted our dataset to fossil occurrences from high latitudes (66°S–90°S), corresponding to the broad paleolatitude of the Fezouata Shale biota for comparative purposes. We then standardized these data to include only genera found in conventional fossiliferous deposits, including biomineralizing clades such as mollusks, brachiopods, trilobites, and echinoderms (Fig. 1A–J), as well as non-biomineralizing genera that have highly recalcitrant structures (Supplementary Dataset 2). Highly recalcitrant genera are mainly represented by graptolites, whose tubaria are presumably collagenous and are a substantial component of many Ordovician strata without soft-tissue preservation (Maletz et al. 2016) (Fig. 1K–M). This sampling strategy, alongside the use of Sørensen's index, reduces the impact that rare outlier taxa and/or those that are only able to fossilize under exceptional conditions found in the Fezouata Shale have on our results. Our data selection ultimately allows direct comparisons between a variety of fossil deposits, regardless of the presence or degree of soft-tissue preservation.

Figure 1. Representative biomineralizing and recalcitrant fossil taxa from the Lower Ordovician Fezouata Shale Formation of Morocco. A, The sponge Choia sp. (YPM.IP.226567). B, The mollusk Pelecyogyra (YPM.IP.518716). C, The stylophoran echinoderm Thoralicystis (MCZ.IP.201057). D, Lingulid brachiopod (MCZ.IP.201058). E, The trilobite Bavarilla (MCZ.IP.201059). F, The conulariid Eoconularia sp. attached to a brachiopod (YPM.IP.530270). G, The hyolith Pauxillites (MCZ.IP.201060). H, The cnidarian Sphenothallus (YPM.IP.226558). I, The eocrinoid Balantiocystis (MCZ.IP.201061). J, The stylophoran Chauvelicystis (MCZ.IP.202501). K, The planktic graptolite Tetragraptus (MCZ.IP.201062). L, The planktic graptolite Clonograptus (MCZ.IP.201063). M, The benthic graptolite Dictyonema (MCZ.IP.201064). Institutional abbreviations: MCZ.IP: Invertebrate Paleontology Collection at the Museum of Comparative Zoology, Harvard University, Cambridge, Mass., U.S.A.; YPM.IP. Invertebrate Paleontology Collection at the Yale Peabody Museum, New Haven, Conn., U.S.A. Scale bars, (A) 5 cm; (B, C, E–H, K–M) 1 cm; (D, I, J) 5 mm.

Fossil occurrences that originate from the same formation may have different sets of paleo-coordinates because they belong to different time-averaged subassemblages, enabling us to analyze patterns of spatial variation between individual time-averaged subassemblages. For each geological formation and its corresponding subassemblages (i.e., with a unique set of paleo-coordinates), fossil occurrences were pooled to reconstruct time-averaged subassemblages. We converted the dataset to genus-level presence–absence occurrence data to reduce the impact that collector's bias and poor sampling may have on subsequent analysis (Chao et al. 2005; Whitaker and Kimmig 2020; Nanglu and Cullen 2023); we also discarded time-averaged subassemblages with fewer than five fossil occurrences.

The total dataset includes 491 fossil occurrences consisting of 227 different genera originating from 31 localities and 11 geological formations. The geological formations represented include: the Borrachon, the Dere, and the Santed of Spain (Álvaro and Martínez-Benítez 2023), the Couches de Barroubio, the Couches de la Maurerie, and the Saint-Chinian of France (Tortello et al. 2006), the Trenice and Milina of the Czech Republic (Kraft et al. 2014), the Wysockzi of Poland (Modliński and Szymański 2001), and the lower and upper Fezouata Shale of Morocco (Van Roy et al. 2010, 2015a). We treat the lower Fezouata Shale (Tremadocian) and the upper Fezouata Shale (Floian) as distinct geological formations to facilitate the discussion of their temporal and taxonomic differences. These data were arranged into a presence/absence matrix (Supplementary Dataset 3) that was used in all subsequent analyses.

Analysis of Biomineralizing and Conventional Fossil Assemblages

The removal of soft-bodied taxa that are unlikely to be preserved except under exceptional conditions such as those found at Fezouata Shale Formation allows us to mitigate the effects of taphonomic bias from our comparisons between fossil biotas. However, even within this highly conservative dataset, there are still differences in the fossilization potential between taxa. This is most notable when considering completely biomineralizing taxa (e.g., brachiopods, trilobites, and echinoderms; Fig. 1A–J) relative to organisms with recalcitrant but non-biomineralized structures, particularly graptolites (Fig. 1K–M). To examine the effect of non-biomineralizing recalcitrant taxa in comparisons between high-latitude Lower Ordovician marine fossil biotas, we performed two iterations of our set of analyses. First, we investigated the faunal structure and dissimilarity between communities for the initial dataset that combines both the biomineralizing taxa as well as fossils with recalcitrant structures of analyses, the latter of which are represented by graptolites and one problematic taxon known as Marcusodictyon (Vinn 2016); we refer to this as the “conventional” fossil assemblage in our “Results” and “Discussion.” Second, we investigated the same questions as before, but using a more exclusive dataset consisting of strictly biomineralizing taxa; we refer to this as the “biomineralized” fossil assemblage dataset in our “Results” and “Discussion.” The entire set of standards for dataset partitioning is detailed in Figure 2, and the relative abundance of conventional and biomineralized fossil taxa used in the analyses are shown graphically in Figure 3.

Figure 2. Schematic of data partitioning including conventional fossils (top) and strictly biomineralizing taxa (bottom). Analyses were repeated in both scenarios to explore the effects of common but not biomineralizing taxa (e.g., graptolites) on the inference of broadscale diversity patterns.

Figure 3. Relative abundance of conventional fossil genera (A) and biomineralizing fossil genera (B) in the studied dataset. Diagram depicts the percentage of major animal groups (e.g., arthropods, mollusks, echinoderms) alongside the percentages of lower-level, but taxonomically important, clades within those larger groups (e.g., asaphid trilobites, cephalopod mollusks, asterozoan echinoderms).

Time-Averaged Assemblage Dissimilarity, Data Visualization, and Statistical Analysis

First, we compared the time-averaged subassemblages of the same formation to one another to provide a measure of intra-formational variability using Sørensen's dissimilarity index. Similar in function to the Jaccard index, Sørensen's index was specifically chosen due to its resilience to outlier-induced inflations that can impact the interpretation of beta-diversity comparisons between communities (Koleff et al. 2003; Hao et al. 2019). At least three pairwise comparisons between subassemblages are needed to create a basic distribution and box plot, so this analysis was restricted to geological formations represented in the PBDB by at least three different fossiliferous subassemblages. From the total dataset, the only geological formations from the high latitudes of the Lower Ordovician that fulfill these requirements are the lower Fezouata Shale Formation (Tremadocian), the upper Fezouata Shale Formation (Floian), the Trenice Formation (Tremadocian), and the Saint-Chinian Formation (Tremadocian). We used box plots to visualize the intra-formational variability for each geological formation. We used Kruskal-Wallis testing and post hoc Dunn tests with Bonferroni corrections to determine whether the variability within these high-latitude formations from the Lower Ordovician were statistically different.

Next, we compared all the subassemblages (N _conventional = 31; N _{biomineralized} = 30) from the total dataset (N = 11 formations) with each other in terms of their composition using Sørensen's dissimilarity index. We used dissimilarity metrics to produce a nonmetric dimensional scaling (NMDS) ordination wherein subassemblages were labeled as either lower Fezouata Shale Formation (N = 1), upper Fezouata Shale Formation (N = 1), or pooled into a non-Fezouata category (N = 9) to explore the degree of compositional similarity between these fossil biotas preserved in these formations. We preformed tests of homogeneity of multivariate spread (betadisper) and permutational analysis of variance (PERMANOVA) with Bonferroni corrections to determine whether the allocated groupings and their spread had any statistical significance. We also used these dissimilarity metrics to produce an NMDS in which all subassemblages were partitioned based on their age instead of formation, categorized as either originating from the Tremadocian or the Floian. Here, we used tests of homogeneity of multivariate spread and PERMANOVA to determine whether the allocated aged-based groupings (Tremadocian vs. Floian) and their spread had any statistical significance.

Results

Intraformational Variability and Fezouata Shale Diversity Metrics

Only four formations (lower Fezouata Shale, upper Fezouata Shale, Trenice, and Saint-Chinian) originating from the high latitudes of the Lower Ordovician contained enough subassemblages (at least three) needed for analyzing the intra-formational variability (Fig. 4A,B). The subassemblages of the lower Fezouata Shale Formation are less like one another than those of the upper Fezouata Shale Formation, as the median dissimilarity between lower Fezouata Shale subassemblages is 0.96–0.95, while the median dissimilarity between subassemblages of the upper Fezouata Shale Formation is 0.79–0.76 (Table 1). The subassemblages of the lower Fezouata Shale Formation exhibited the lowest spread in dissimilarity values among all the analyzed formations (interquartile range [IQR] = 0.1–0.09), while those of the upper Fezouata Shale Formation had the highest (IQR = 0.23–0.22). We did not detect statistical differences in variability between any of the four formations for the analysis of conventional fossil taxa (pairwise Dunn test, p > 0.05). By contrast, the analysis restricted to biomineralized taxa shows that the variability of the upper Fezouata Shale Formation was statistically different from that of both the lower Fezouata Shale Formation (Dunn test, p = 0.0298) and the Saint-Chinian Formation (Dunn test, p = 0.0234) but not the Trenice Formation (Dunn test, p = 1.00). The failure of the analysis to find significance between the Trenice Formation and other formations is likely due to the limited number of reconstructed subassemblages (N = 3).

Figure 4. Box plots showing intra-formational variability as measured by Sørensen's index. A, Analysis of conventional fossil taxa. B, Analysis restricted to strictly biomineralizing fossil taxa. Red lines and asterisks represent a statistical difference between the two formations.

Table 1. Data distribution of intra-formational variability as measured by Sørensen's index. Includes median dissimilarity and information regarding the spread (IQR) for analysis of conventional and biomineralizing taxa. IQR, interquartile range.

The median α-diversity of the lower Fezouata Shale formation was 11.5 and 9.5 for analysis of conventional and biomineralizing taxa, while the median α-diversity of the Upper Fezouata Shale was 11 for both analysis of conventional and biomineralizing taxa (Table 2). Given that dissimilarity is a proxy for β-diversity (Beck et al. 2013; Ricotta 2017), the median β-diversity of the lower Fezouata Shale formation ranged from 0.96 to 0.95 (for analysis of conventional and biomineralizing taxa, respectively), while the median β-diversity of the upper Fezouata Shale formation ranged from 0.79 to 0.76 (for conventional and biomineralizing taxa, respectively) (Tables 1, 2). The γ-diversity of the lower Fezouata Shale formation was 45 for analysis of conventional taxa and dropped to 40 for analysis of biomineralizing taxa (Table 2). The γ-diversity was higher for the upper Fezouata Shale formation and ranges from 58 for analysis of conventional taxa to 52 when restricted to biomineralizing taxa (Table 2).

Table 2. Measures of diversity (α, β, γ) within the lower and upper Fezouata Formation for analysis of conventional and biomineralizing taxa. Includes median β-diversity (dissimilarity) as measured by the Sørensen's index.

Ordination by Fezouata Membership and Geological Stage

To make comparisons among high-latitude Lower Ordovician fossil biotas more broadly, we employed a three-dimensional NMDS (stress ~ 0.058) outlined by membership in the lower Fezouata Shale Formation, the upper Fezouata Shale Formation, or non-Fezouata formations (N = 9), including analyses considering conventional and biomineralizing fossils (Fig. 5). We find significant differences in spread between the lower Fezouata, upper Fezouata, and non-Fezouata clusters for the analysis of conventional taxa (F-test = 6.30, p = 0.00553) and biomineralizing taxa (F-test = 6.49, p = 0.00499). There is overlap between the hulls encompassing subassemblages of the lower Fezouata Shale and the non-Fezouata sites for the analysis of conventional (pairwise PERMANOVA, F-test =1.20, R ² = 0.0592, p = 0.645) and biomineralizing taxa (pairwise PERMANOVA, F-test = 1.23, R ² = 0.0608, p = 0.609) along all three axes, with much of the lower Fezouata Shale hull space falling within the hull space of the non-Fezouata sites (Fig. 5A,B, Supplementary Figs. S1, S3, S4). There is little overlap along the second ordinal axis between the convex polygon hulls of the upper Fezouata Shale subassemblages and those of the lower Fezouata Shale for the analysis of conventional (pairwise PERMANOVA, F-test = 1. 63, R ² = 0.120, p = 0.162) and biomineralizing taxa (pairwise PERMANOVA, F-test = 1.71, R ² = 0.134, p = 0.171) (Fig. 5A,B, Supplementary Figs. S1, S3, S4). There is limited overlap along the second ordinal axis between the convex polygon hulls of the upper Fezouata Shale subassemblages and the non-Fezouata sites for the analysis of conventional (pairwise PERMANOVA, F-test = 3.733, R ² = 0.130, p = 0.003) and biomineralizing taxa (pairwise PERMANOVA, F-test = 3.84, R ² = 0.138, p = 0.003) (Fig. 5A,B, Supplementary Figs. S1, S3, S4). These patterns are consistent between both our conventional and biomineralized datasets. Failure of the analysis to find significant differences between lower Fezouata Shale and other formations can be attributed to the fact that we only recovered four lower Fezouata Shale subassemblages, as most statistical tests gain statistical power as sample size increases (Banerjee et al. 2009; Haas 2012).

Figure 5. Nonmetric dimensional scaling (NMDS) for the high latitudes of the Lower Ordovician showing convex hulls for subassemblages from the lower Fezouata Shale Formation (orange), upper Fezouata Shale formation (yellow), and non-Fezouata (green) along NMDS 1 and NMDS 2. A, Analysis of conventional fossil taxa. B, Analysis of strictly biomineralized fossil taxa. Reduction of the number of dimensions leads to a stress of approximately 0.058. Additional two-dimensional visualizations for the ordination are included as Supplementary Fig. 1.

There are significant differences in spread between Tremadocian and Floian clusters for the analyses of conventional (F-test = 5.72, p = 0.0235) and biomineralizing taxa (F-test = 5.92, p = 0.0216). Fossil subassemblages corresponding to the Tremadocian (e.g., lower Fezouata Shale Formation, Saint-Chinian, and Trenice) and Floian (e.g., upper Fezouata Shale Formation) are separated along the second axis in ordinal space for the analyses involving conventional taxa (PERMANOVA, F-test = 3.10, R ² = 0.0965, p = 0.001) and those restricted to biomineralizing taxa (PERMANOVA, F-test = 3.07, R ² = 0.0987, p = 0.001) (Fig. 6A,B, Supplementary Figs. S2, S3, S4).

Figure 6. Nonmetric dimensional scaling (NMDS) the high latitudes of the Lower Ordovician showing convex hulls for subassemblages from the Tremadocian (purple) and Floian (pink) along NMDS 1 and NMDS 2. A, Analysis of conventional fossil taxa. B, Analysis of strictly biomineralized taxa. Reduction of the number of dimensions leads to a stress of roughly 0.058. Additional two-dimensional visualizations for the ordination are included as Supplementary Fig. 2.

Discussion

Bias in the PBDB for the High Latitudes of the Early Ordovician

It is crucial to emphasize that despite our efforts to exhaustively sample information in the PBDB from the high latitudes of the Early Ordovician ecosystems, the subassemblages with the necessary occurrence information to be included are restricted geographically to Gondwana and peri-Gondwana (Iberia and Armorica) (Modliński and Szymański 2001; Tortello et al. 2006; Kraft et al. 2014; Van Roy et al. 2015a; Álvaro and Martínez-Benítez 2023). This pattern likely derives from the impact of colonization and neo-colonization on global inequality (e.g., material resources, stability, and economic and social development), which ultimately affects occurrence data uploaded into online repositories (Amano and Sutherland 2013; Trisos et al. 2021; Raja et al. 2022). More specifically, the bias we uncover points to both an overrepresentation of sampling effort in sites from the Global North (United States, Canada, Europe, etc.) coupled with the overrepresentation of science conducted by researchers from those regions (Raja et al. 2022). Though this study encompasses as much information as possible from the high latitudes of the Early Ordovician, it is important to recognize that this limited geographic range alongside the biases in global sampling likely impacts the results of our study (Whitaker and Kimmig 2020; Benson et al. 2022; Nanglu and Cullen 2023).

Heterogeneity in High-Latitude Fossil Biotas from the Lower Ordovician

All the formations with enough fossiliferous subassemblages to allow for intra-formational comparisons (Table 1) showed evidence of moderate to high heterogeneity, with median intra-formational dissimilarity ranging from 75–79% (i.e., upper Fezouata Shale Formation) to 95–96% (i.e., lower Fezouata Shale Formation) (Fig. 4). The high β-diversity values in our results concur with recent studies showing high dissimilarity throughout the Early Ordovician on a global scale (Penny and Kröger 2019). It is important to note that the Early Ordovician is marked by the beginning of environmental changes and increases in animal dispersal that would eventually become more pronounced during the Middle Ordovician (Miller 1997; Servais et al. 2014; Lee and Riding 2018; Edwards 2019; Lam et al. 2021). In this context, the patterns we recovered may suggest that the ecosystems of the Early Ordovician in the higher latitudes were still dictated by processes reminiscent of those taking place during the Cambrian to some extent (Servais et al. 2014; Rasmussen et al. 2016; Lee and Riding 2018; Saleh et al. 2022a). It is also relevant to consider that, despite the appearance of bioturbators defining the start of the Phanerozoic Eon (Buatois et al. 2020), the Early Ordovician precedes the major radiations of bioturbating organisms that would come to define the modern marine biosphere (Servais et al. 2010; Tarhan 2018). The high β-diversity recovered may also reflect a lack of time averaging introduced via bioturbation that is widespread further into the Phanerozoic (Kidwell 1997).

Homogenization within the Fezouata Shale Formation

The differences in the raw median variability observed between the lower Fezouata Shale and the upper Fezouata Shale indicate a clear pattern of increasing homogenization through time (Fig. 4, Table 1). The raw median variability in the upper Fezouata Shale (0.76) is 19% more homogenous relative to the lower Fezouata Shale (0.95), based on comparisons involving strictly biomineralized taxa (Dunn test, p = 0.0298).

It is notable that statistical differences in intra-formational variability are found in the analysis of the biomineralized dataset, but not in the analysis of the conventional dataset (Fig. 4). This likely results from the increase in intra-formational variability seen in the upper Fezouata Shale when including conventional taxa (median intra-formational variability = 0.793) versus only considering biomineralizing taxa (median intra-formational variability = 0.755) (Table 1). Given that recalcitrant non-biomineralized taxa make up approximately 5% of all analyzed fossil occurrences (Fig. 3), 96% of which correspond to graptolites, this result indicates that graptolites are likely undersampled and disproportionally skew the data so that the upper Fezouata appears more heterogenous when they are included. This is further supported by similar p-values recovered for comparisons in intra-formational variability between the upper and lower Fezouata in the conventional (Dunn test, p = 0.0524) and biomineralized analyses (Dunn test, p = 0.0298).

A previous study focusing on bivalves and brachiopods reported that the lower Fezouata Shale is characterized by low α-diversity and high γ-diversity (Saleh et al. 2018). Our results support this conclusion, as we find that the median α-diversity for a conventional fossil assemblage in the lower Fezouata Shale was 11.5 for analysis including conventional taxa and 9.5 for analysis including biomineralizing taxa (Table 2). The γ-diversity for all four subassemblages of the lower Fezouata Shale Formation summed together for analysis of conventional taxa and biomineralizing taxa are 45 and 40 taxa, respectively (Table 2). A similar pattern of a low α- to γ-diversity ratio exists for both treatments of the upper Fezouata dataset (Table 2). Considering this information, the interpretation for the presence of dominant and opportunistic short-lived populations of mollusks and brachiopods could be extended to the entire recalcitrant fauna, also including taxa such as trilobites, echinoderms, and graptolites (Saleh et al. 2018).

Although it is difficult to establish a single cause, the increase in homogenization observed in the upper Fezouata Shale Formation (Fig. 4) most likely results from a combination of changes in ocean currents, temperature, chemistry, and tectonic movements that would have increased the dispersal potential of many animal groups in the early Paleozoic (Pruss et al. 2010; Saltzman et al. 2015; Algeo et al. 2016; Liljeroth et al. 2017; Lam et al. 2018). It has also been proposed that high nutrient runoff rates may have led to radiations among phytoplankton, providing additional sources of food for marine heterotrophs attempting to avoid benthic predators (Servais et al. 2010, 2016). Likewise, the growing colonization of the water column around the same time of increased ocean circulatory activity most likely allowed organisms to disperse far and wide (Poussart et al. 1999; Servais et al. 2010; Pohl et al. 2016; Rasmussen et al. 2016; Cocks and Torsvik 2021), contributing to more homogenous ecosystems (Saleh et al. 2022a).

The homogenization of the Fezouata Shale Formation can be partially explained by examining the changes in the lifestyles of the trilobites in our dataset (Chatterton and Speyer 1989; Speyer and Chatterton 1989; Signor and Vermeij 1994; Hughes et al. 2006; Fortey 2014; Laibl et al. 2023), given the known colonization of the nektonic and planktonic realm by many animal groups throughout the late Cambrian and Early Ordovician (Servais et al. 2010; Liu et al. 2022). We find that 37% of the trilobite genera recovered from the subassemblages of the lower Fezouata Shale Formation belong to families known to have planktonic postembryonic larval stages (e.g., Asaphidae, Calymenidae, Nileidae), whereas 47% of the trilobite genera recovered from the upper Fezouata Shale Formation belong to families known to have planktonic postembryonic larval stages (e.g., Asaphidae, Calymenidae, Raphiophoridae). In fact, three of the four most widely distributed trilobites (present in at least 50% subassemblages) of the upper Fezouata Shale Formation—Ampyx, Colpocoryphe, and Asaphellus—belong to families with planktonic larval stages (Laibl et al. 2023). While dispersal influenced the movement and distribution of trilobites throughout the Ordovician (Chatterton and Speyer 1989; Perrier et al. 2015; Lam et al. 2018), the communities of the central Anti-Atlas region appear controlled by these dispersal forces to a lesser degree than many other neighboring regions (Saleh et al. 2022a), as reflected in the slight percentage increase in taxa with planktonic larvae we observe between the Fezouata Shale sections (+10%). Regardless of the magnitude of this increase, this pattern readily ties into recent work showing that trilobites and other animal clades with planktonic larval stages became a larger component of marine ecosystems during the Ordovician (Peterson 2005; Nützel et al. 2006; Laibl et al. 2023) and indicate that ecological turnover can be detected between high-resolution timescales, such as the transition between the Tremadocian and Floian.

A Lagerstätte View of Ordovician Marine Systems

The results of the ordination show little overlap along the second axis of the ordination between the lower and upper Fezouata Shale Formation (Fig. 5A,B); however, we find no statistical differences in faunal composition between the lower and upper Fezouata Shale Formation for either the analysis including all conventional taxa or restricted to biomineralizing taxa (pairwise PERMANOVA, R ² ~ 0.127, p ~ 0.167). This discrepancy between the graphical and statistical analyses likely results from the small number of subassemblages recovered from the lower Fezouata Shale Formation (n = 4) available in the PBDB, which would limit the test's power to detect differences between it and other fossil subassemblages (Banerjee et al. 2009; Haas 2012).

While there are numerous shared taxa between the two sections (e.g., the homalozoan echinoderm Anatifopsis, the trilobite Asaphellus, and tergomyan mollusk Carcassonnella), the lack of overlap along the second axis (Fig. 5A,B, Supplementary Fig. S1) and therefore total three-dimensional ordinal space, likely reflects the faunal distinctions between the lower Fezouata Shale and upper Fezouata Shale. These results quantitatively corroborate assessments based on systematic and biostratigraphic work suggesting some faunal differences in the two main intervals with exceptional preservation in the Fezouata Shale Formation (Van Roy et al. 2015a; Gutiérrez-Marco and Martin 2016; Lefebvre et al. 2016b). Studies point to either environmental changes or rapid colonization of the seafloor to explain differences in assemblage composition between the sections of the Fezouata Shale Formation (Lefebvre et al. 2016a; Martin et al. 2016b). Faunal differences were likely further emphasized by changes occurring during the Early Ordovician, such as tectonic activity, increased ecological complexity, and changes in seawater temperature and chemistry (Servais et al. 2008, 2016; Lee and Riding 2018; Servais and Harper 2018; Cocks and Torsvik 2021).

The difference in faunal composition between the Tremadocian and Floian subassemblages recovered in the analysis (PERMANOVA, R ² ~ 0.098, p < 0.05; Fig. 6A,B, Supplementary Fig. S2) may reflect a genuine biological signal, as previous works have noted some turnover and changes in the structure and composition of brachiopod, echinoderm, and trilobite communities in the Early Ordovician (Hansen and Holmer 2010; Balseiro and Waisfeld 2013, 2014; Lam et al. 2021; Laibl et al. 2023). This pattern may also be the result of differential sampling of the Tremadocian and Floian stages in our PBDB data, which is still a ubiquitous issue even in the age of large-scale paleontological datasets (Nanglu and Cullen 2023). For example, from a total of 12 Floian-aged fossil subassemblages included in our dataset, 10 belong to the upper Fezouata Shale Formation; indicating that Floian fossil subassemblages are not currently well represented in the PBDB. Thus, whether the fossil subassemblages of the upper Fezouata Shale Formation are truly representative of the Floian remains to be confirmed. The paucity of non-Fezouata Floian subassemblages that met the minimum diversity, temporal, and latitudinal thresholds for inclusion in our dataset demonstrates that the upper Fezouata Shale provides ones of the most complete snapshots into the community structure of high-latitude animal communities during the Floian available to date.

If the broader differences between Tremadocian- and Floian-aged deposits are ecological in origin and not entirely a result of bias, this pattern further emphasizes the differences in faunal composition between the upper Fezouata Shale (Floian) and both the lower Fezouata Shale (Tremadocian) and other Tremadocian-aged formations. Likewise, the similarities in assemblage composition between the lower Fezouata Shale Formation and other Tremadocian-aged deposits (pairwise PERMANOVA, R ² ~ 0.06, p ~ 0.627; Fig. 5A,B, Supplementary Fig. S1) provide the first assemblage-wide quantitative evidence that the conventional fossil subassemblages of the Fezouata Shale Formation are comparable with other high-latitude marine deposits of the same stratigraphic age (Van Roy et al. 2010, 2015a; Kröger and Lefebvre 2012; Lefebvre et al. 2016b,c). More broadly, our results also suggest that the conventional fossil subassemblages preserved in the Fezouata Shale Formation, at least in the lower sections (Tremadocian in age), are representative of high-latitude Early Ordovician deposits. Considering that the Fezouata Shale biota reflects all major ecological tiering guilds (Saleh et al. 2021a), we hypothesize that it is highly likely that the exceptionally preserved organisms from the Fezouata Shale biota excluded from our analyses might also be similarly representative of high-latitude marine life at the time.

Conclusion

In summary, we find robust patterns for increased homogeneity in the Fezouata Shale Formation through time. Our results represent the first quantitative evidence that the shelly faunal subassemblages of the Fezouata Shale Formation are typical for a high-latitude Early Ordovician marine deposit. This study provides a baseline for future studies planning to include rarer soft-bodied taxa in examinations of the biodiversity and ecology of the Early Ordovician polar biosphere in the critical time frame before the GOBE.