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. Reference Van Roy, Orr, Botting, Muir, Vinther, Lefebvre, Hariri and Briggs2010, Reference Van Roy, Briggs and Gaines2015a; Lefebvre et al. Reference Lefebvre, Lerosey-Aubril, Servais and Van Roy2016c). 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. Reference Van Roy, Briggs and Gaines2015a; Lefebvre et al. Reference Lefebvre, El Hariri, Lerosey-Aubril, Servais and Van Roy2016b). 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. Reference Van Roy, Orr, Botting, Muir, Vinther, Lefebvre, Hariri and Briggs2010, Reference Van Roy, Briggs and Gaines2015a; Martin et al. Reference Martin, Pittet, Gutiérrez-Marco, Vannier, El Hariri, Lerosey-Aubril and Masrour2016a; Saleh et al. Reference Saleh, Vaucher, Antcliffe, Daley, El Hariri, Kouraiss and Lefebvre2021b, Reference Saleh, Guenser, Gibert, Balseiro, Serra, Waisfeld and Antcliffe2022a,Reference Saleh, Vaucher, Vidal, Hariri, Laibl, Daley and Gutiérrez-Marcob). 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 Reference Gutiérrez-Marco and Martin2016; Lefebvre et al. Reference Lefebvre, Gutiérrez-Marco, Lehnert, Martin, Nowak, Akodad, El Hariri and Servais2018; Saleh et al. Reference Saleh, Candela, Harper, Polechová, Lefebvre and Pittet2018).
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. Reference Van Roy, Briggs and Gaines2015a; Lefebvre et al. Reference Lefebvre, El Hariri, Lerosey-Aubril, Servais and Van Roy2016b). 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. Reference Van Roy, Daley and Briggs2015b; Legg Reference Legg2016; Ortega-Hernández et al. Reference Ortega-Hernández, Van Roy and Lerosey-Aubril2016; Vinther et al. Reference Vinther, Parry, Briggs and Van Roy2017; Hunter and Ortega-Hernández Reference Hunter and Ortega-Hernández2021; Pérez-Peris et al. Reference Pérez-Peris, Laibl, Lustri, Gueriau, Antcliffe, Bath Enright and Daley2021). 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 Reference Droser2003; Servais et al. Reference Servais, Owen, Harper, Kröger and Munnecke2010, Reference Servais, Cascales-Miñana, Harper, Lefebvre, Munnecke, Wang and Zhang2023; Servais and Harper Reference Servais and Harper2018; Stigall Reference Stigall2018; Stigall et al. Reference Stigall, Edwards, Freeman and Rasmussen2019), 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. Reference Servais, Cascales-Miñana, Harper, Lefebvre, Munnecke, Wang and Zhang2023).
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. Reference Saleh, Pittet, Sansjofre, Guériau, Lalonde, Perrillat and Vidal2020, Reference Saleh, Vaucher, Antcliffe, Daley, El Hariri, Kouraiss and Lefebvre2021b) below storm-wave base on an open siliciclastic ramp roughly 70° south of the equator (Saleh et al. Reference Saleh, Pittet, Sansjofre, Guériau, Lalonde, Perrillat and Vidal2020, Reference Saleh, Vaucher, Antcliffe, Daley, El Hariri, Kouraiss and Lefebvre2021b; Martin et al. Reference Martin, Pittet, Gutiérrez-Marco, Vannier, El Hariri, Lerosey-Aubril and Masrour2016a). Data using brachiopods and bivalves suggest that the communities of the Fezouata Shale biota were ephemeral (Saleh et al. Reference Saleh, Candela, Harper, Polechová, Lefebvre and Pittet2018). Work focused on trilobites, graptolites, and cephalopods has provided information on the distribution of these clades within the Fezouata Shale (Kröger and Lefebvre Reference Kröger and Lefebvre2012; Gutiérrez-Marco and Martin Reference Gutiérrez-Marco and Martin2016; Martin et al. Reference Martin, Vidal, Vizcaïno, Vaucher, Sansjofre, Lefebvre and Destombes2016b). This, alongside publications summarizing faunal lists for the deposit (Van Roy et al. Reference Van Roy, Briggs and Gaines2015a; Lefebvre et al. Reference Lefebvre, El Hariri, Lerosey-Aubril, Servais and Van Roy2016b), 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. Reference Nanglu, Waskom, Richards and Ortega-Hernández2023). Although it has been shown that the Fezouata Shale Formation is dominated by echinoderms, arthropods, graptolites, and brachiopods (Kouraiss et al. Reference Kouraiss, El Hariri, El Albani, Azizi, Mazurier and Lefebvre2019; Saleh et al. Reference Saleh, Vaucher, Vidal, Hariri, Laibl, Daley and Gutiérrez-Marco2022b), 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. Reference Van Roy, Orr, Botting, Muir, Vinther, Lefebvre, Hariri and Briggs2010, Reference Van Roy, Briggs and Gaines2015a; Kröger and Lefebvre Reference Kröger and Lefebvre2012; Lefebvre et al. Reference Lefebvre, Lerosey-Aubril, Servais and Van Roy2016c, Reference Lefebvre, Gutiérrez-Marco, Lehnert, Martin, Nowak, Akodad, El Hariri and Servais2018). 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. Reference Saleh, Guenser, Gibert, Balseiro, Serra, Waisfeld and Antcliffe2022a). 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. Reference Servais, Owen, Harper, Kröger and Munnecke2010; Servais and Harper Reference Servais and Harper2018) 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 Reference Wickham2011, Reference Wickham2020; Wickham et al. Reference Wickham, Hester and Francois2017, Reference Wickham, François, Henry, Müller and Vaughan2023a,Reference Wickham, Vaughan and Girlichb; Firke et al. Reference Firke, Denney, Haid, Knight, Grosser and Zadra2020; R Core Team 2022). We computed ecological metrics in RStudio using the package vegan (Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn and Minchin2019). We carried out statistics and visualization in RStudio using the packages vegan, stats, dunn.test and funfuns, ggpol, ggplot2, and colorRamps (Wickham Reference Wickham2009; Keitt Reference Keitt2012; Dinno Reference Dinno2017; Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn and Minchin2019; Tiedemann Reference Tiedemann2020; R Core Team 2022; Trachsel Reference Trachsel2023). 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. Reference Maletz, Lenz and Bates2016) (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.
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. Reference Chao, Chazdon, Colwell and Shen2005; Whitaker and Kimmig Reference Whitaker and Kimmig2020; Nanglu and Cullen Reference Nanglu and Cullen2023); 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 Reference Álvaro and Martínez-Benítez2023), the Couches de Barroubio, the Couches de la Maurerie, and the Saint-Chinian of France (Tortello et al. Reference Tortello, Vizcaïno and Álvaro2006), the Trenice and Milina of the Czech Republic (Kraft et al. Reference Kraft, Mergl, Hroch and Kraft2014), the Wysockzi of Poland (Modliński and Szymański Reference Modliński and Szymański2001), and the lower and upper Fezouata Shale of Morocco (Van Roy et al. Reference Van Roy, Orr, Botting, Muir, Vinther, Lefebvre, Hariri and Briggs2010, Reference Van Roy, Briggs and Gaines2015a). 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 Reference Vinn2016); 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.
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.
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. Reference Koleff, Gaston and Lennon2003; Hao et al. Reference Hao, Corral-Rivas, González-Elizondo, Ganeshaiah, Nava-Miranda, Zhang, Zhao and Von Gadow2019). 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).
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.
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. Reference Beck, Holloway and Schwanghart2013; Ricotta Reference Ricotta2017), 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).
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 2 = 0.0592, p = 0.645) and biomineralizing taxa (pairwise PERMANOVA, F-test = 1.23, R 2 = 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 2 = 0.120, p = 0.162) and biomineralizing taxa (pairwise PERMANOVA, F-test = 1.71, R 2 = 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 2 = 0.130, p = 0.003) and biomineralizing taxa (pairwise PERMANOVA, F-test = 3.84, R 2 = 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. Reference Banerjee, Chitnis, Jadhav, Bhawalkar and Chaudhury2009; Haas Reference Haas2012).
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 2 = 0.0965, p = 0.001) and those restricted to biomineralizing taxa (PERMANOVA, F-test = 3.07, R 2 = 0.0987, p = 0.001) (Fig. 6A,B, Supplementary Figs. S2, S3, S4).
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 Reference Modliński and Szymański2001; Tortello et al. Reference Tortello, Vizcaïno and Álvaro2006; Kraft et al. Reference Kraft, Mergl, Hroch and Kraft2014; Van Roy et al. Reference Van Roy, Briggs and Gaines2015a; Álvaro and Martínez-Benítez Reference Álvaro and Martínez-Benítez2023). 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 Reference Amano and Sutherland2013; Trisos et al. Reference Trisos, Auerbach and Katti2021; Raja et al. Reference Raja, Dunne, Matiwane, Khan, Nätscher, Ghilardi and Chattopadhyay2022). 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. Reference Raja, Dunne, Matiwane, Khan, Nätscher, Ghilardi and Chattopadhyay2022). 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 Reference Whitaker and Kimmig2020; Benson et al. Reference Benson, Brown, Campione, Cullen, Evans and Zanno2022; Nanglu and Cullen Reference Nanglu and Cullen2023).
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 Reference Penny and Kröger2019). 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 Reference Miller1997; Servais et al. Reference Servais, Danelian, Harper and Munnecke2014; Lee and Riding Reference Lee and Riding2018; Edwards Reference Edwards2019; Lam et al. Reference Lam, Sheffield and Matzke2021). 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. Reference Servais, Danelian, Harper and Munnecke2014; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen and Eriksson2016; Lee and Riding Reference Lee and Riding2018; Saleh et al. Reference Saleh, Guenser, Gibert, Balseiro, Serra, Waisfeld and Antcliffe2022a). It is also relevant to consider that, despite the appearance of bioturbators defining the start of the Phanerozoic Eon (Buatois et al. Reference Buatois, Mángano, Minter, Zhou, Wisshak, Wilson and Olea2020), the Early Ordovician precedes the major radiations of bioturbating organisms that would come to define the modern marine biosphere (Servais et al. Reference Servais, Owen, Harper, Kröger and Munnecke2010; Tarhan Reference Tarhan2018). The high β-diversity recovered may also reflect a lack of time averaging introduced via bioturbation that is widespread further into the Phanerozoic (Kidwell Reference Kidwell1997).
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. Reference Saleh, Candela, Harper, Polechová, Lefebvre and Pittet2018). 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. Reference Saleh, Candela, Harper, Polechová, Lefebvre and Pittet2018).
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. Reference Pruss, Finnegan, Fischer and Knoll2010; Saltzman et al. Reference Saltzman, Edwards, Adrain and Westrop2015; Algeo et al. Reference Algeo, Marenco and Saltzman2016; Liljeroth et al. Reference Liljeroth, Harper, Carlisle and Nielsen2017; Lam et al. Reference Lam, Stigall and Matzke2018). 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. Reference Servais, Owen, Harper, Kröger and Munnecke2010, Reference Servais, Perrier, Danelian, Klug, Martin, Munnecke and Nowak2016). 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. Reference Poussart, Weaver and Barnes1999; Servais et al. Reference Servais, Owen, Harper, Kröger and Munnecke2010; Pohl et al. Reference Pohl, Nardin, Vandenbroucke and Donnadieu2016; Rasmussen et al. Reference Rasmussen, Ullmann, Jakobsen, Lindskog, Hansen, Hansen and Eriksson2016; Cocks and Torsvik Reference Cocks and Torsvik2021), contributing to more homogenous ecosystems (Saleh et al. Reference Saleh, Guenser, Gibert, Balseiro, Serra, Waisfeld and Antcliffe2022a).
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 Reference Chatterton and Speyer1989; Speyer and Chatterton Reference Speyer and Chatterton1989; Signor and Vermeij Reference Signor and Vermeij1994; Hughes et al. Reference Hughes, Minelli and Fusco2006; Fortey Reference Fortey2014; Laibl et al. Reference Laibl, Saleh and Pérez-Peris2023), given the known colonization of the nektonic and planktonic realm by many animal groups throughout the late Cambrian and Early Ordovician (Servais et al. Reference Servais, Owen, Harper, Kröger and Munnecke2010; Liu et al. Reference Liu, Fan, Zong and Gong2022). 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. Reference Laibl, Saleh and Pérez-Peris2023). While dispersal influenced the movement and distribution of trilobites throughout the Ordovician (Chatterton and Speyer Reference Chatterton and Speyer1989; Perrier et al. Reference Perrier, Williams and Siveter2015; Lam et al. Reference Lam, Stigall and Matzke2018), 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. Reference Saleh, Guenser, Gibert, Balseiro, Serra, Waisfeld and Antcliffe2022a), 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 Reference Peterson2005; Nützel et al. Reference Nützel, Lehnert and Frýda2006; Laibl et al. Reference Laibl, Saleh and Pérez-Peris2023) 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 2 ~ 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. Reference Banerjee, Chitnis, Jadhav, Bhawalkar and Chaudhury2009; Haas Reference Haas2012).
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. Reference Van Roy, Briggs and Gaines2015a; Gutiérrez-Marco and Martin Reference Gutiérrez-Marco and Martin2016; Lefebvre et al. Reference Lefebvre, El Hariri, Lerosey-Aubril, Servais and Van Roy2016b). 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. Reference Lefebvre, Allaire, Guensburg, Hunter, Kouraïss, Martin and Nardin2016a; Martin et al. Reference Martin, Vidal, Vizcaïno, Vaucher, Sansjofre, Lefebvre and Destombes2016b). 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. Reference Servais, Lehnert, Li, Mullins, Munnecke, Nützel and Vecoli2008, Reference Servais, Perrier, Danelian, Klug, Martin, Munnecke and Nowak2016; Lee and Riding Reference Lee and Riding2018; Servais and Harper Reference Servais and Harper2018; Cocks and Torsvik Reference Cocks and Torsvik2021).
The difference in faunal composition between the Tremadocian and Floian subassemblages recovered in the analysis (PERMANOVA, R 2 ~ 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 Reference Hansen and Holmer2010; Balseiro and Waisfeld Reference Balseiro and Waisfeld2013, Reference Balseiro and Waisfeld2014; Lam et al. Reference Lam, Sheffield and Matzke2021; Laibl et al. Reference Laibl, Saleh and Pérez-Peris2023). 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 Reference Nanglu and Cullen2023). 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 2 ~ 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. Reference Van Roy, Orr, Botting, Muir, Vinther, Lefebvre, Hariri and Briggs2010, Reference Van Roy, Briggs and Gaines2015a; Kröger and Lefebvre Reference Kröger and Lefebvre2012; Lefebvre et al. Reference Lefebvre, El Hariri, Lerosey-Aubril, Servais and Van Roy2016b,Reference Lefebvre, Lerosey-Aubril, Servais and Van Royc). 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. Reference Saleh, Bath-Enright, Daley, Lefebvre, Pittet, Vite, Ma, Mángano, Buatois and Antcliffe2021a), 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.
Acknowledgments
This work is funded by NSF CAREER award no. 2047192 “Ecological Turnover at the Dawn of the Great Ordovician Biodiversification Event—Quantifying the Cambro-Ordovician Transition through the Lens of Exceptional Preservation,” and NSF GRFP award no. 2140743. We also acknowledge support from the Department of Organismic and Evolutionary Biology and the Harvard Museum of Comparative Zoology. We thank the Invertebrate Paleontology Collection Staff of both the Museum of Comparative Zoology (J. Cundiff, C. Green, and M. Renczkowski) and the Yale Peabody Museum (J. Utrup) for help in curating and retrieving the fossil specimens illustrated in this study.
Competing Interest
The authors declare no conflict of interest.
Data Availability Statement
The major datasets produced and used for analysis are included as Supplementary Materials (Supplementary Datasets 1, 2, and 3) alongside the R Script created to run code and complete the analysis on Dryad (https://doi.org/10.5061/dryad.1jwstqk26). Additional two-dimensional visualizations alongside Shepard's (i.e., stress) plots for the ordination segregated by formation and age are included as Supplementary Materials (https://doi.org/10.5281/zenodo.10654895).
Open access
