Introduction
The ca. 506-Myr Burgess Shale, located in the Canadian Rockies of British Columbia, is among the most celebrated Konservat-Lagerstätten in the world thanks to the exceptional preservation of abundant and diverse assemblages of soft-bodied animals. For more than 100 yr, the Burgess Shale fauna has provided critical insights into the morphology and evolutionary relationships of some of the oldest members of many metazoan groups within the broader context of the Cambrian radiation of animal life (e.g., Erwin et al. Reference Erwin, Laflamme, Tweedt, Sperling, Pisani and Peterson2011). In addition to these contributions to evolutionary studies, the Burgess Shale and other similar deposits like the Chengjiang biota of South China (Zhao et al. Reference Zhao, Caron, Bottjer, Hu, Yin and Zhu2013) offer a far more complete view of Cambrian marine life than do normal fossil deposits of the same age (Conway Morris Reference Conway Morris1986).
This view has been expanded by the discovery of new Cambrian localities with soft-tissue preservation, in particular in the southern Canadian Rockies (Collins et al. Reference Collins, Briggs and Morris1983; Fletcher and Collins Reference Fletcher and Collins1998, Reference Fletcher and Collins2003; Johnston et al. Reference Johnston, Johnston and Powell2009; Caron et al. Reference Caron, Gaines, Mángano, Streng and Daley2010), underscoring that the Burgess Shale paleocommunity is a mosaic of multiple local communities. However, rigorous, quantitative paleocommunity analyses that incorporate both abundances and stratigraphic data have been largely focused on the Walcott Quarry, the best-sampled and most intensively studied site (however, see Zhao et al. [Reference Zhao, Hu, Caron, Zhu, Yin and Lu2012, Reference Zhao, Caron, Bottjer, Hu, Yin and Zhu2013] for two exceptions from China). These studies, spurred by large amounts of additional fossil material collected by the Royal Ontario Museum in the 1990s, have not only shown how varied the Burgess Shale is in both taxonomic composition and diversity of life habits, but also how the paleocommunity changed throughout the Walcott Quarry interval (Caron and Jackson Reference Caron and Jackson2006, Reference Caron and Jackson2008). Subsequent comparisons with the Tulip Beds, located approximately 4 km southwest from the Walcott Quarry on Mount Stephen, revealed that only half of the genera of the Tulip Beds are shared with the Walcott Quarry (O'Brien and Caron Reference O'Brien and Caron2015). Despite these differences, however, the major taxonomic groups (Arthropoda, Porifera, Priapulida, and Brachiopoda) and ecological modes (epibenthic suspension feeders and nektobenthic predators) were broadly similar at the two localities (O'Brien and Caron Reference O'Brien and Caron2015). Taken together, these results suggested that while environmental conditions may have caused some change in the ecological structure of the Burgess Shale paleocommunity over short intervals, at a broader scale (including older sites in China) Cambrian faunas were relatively stable in structure across large geographic and temporal scales (O'Brien and Caron Reference O'Brien and Caron2015).
These recent studies have begun to elucidate a broader picture of Cambrian paleocommunity ecology. However, to fully investigate local and regional variations of the Burgess Shale paleocommunity composition, as well as the ecological and environmental factors responsible for those variations, analyses including additional localities that provide a greater breadth of both time and spatial extent are required. The lack of a unified data set incorporating fine-scale temporal variations in composition between multiple Burgess Shale localities has also precluded a broader quantitative assessment of the Burgess Shale paleocommunity as a whole. Consequently, the utility of the Burgess Shale as a point of comparison with modern ecological dynamics (faunal turnover, spatial and temporal heterogeneity, etc.) has remained limited.
Here we provide a detailed overview of the paleocommunity diversity and structure of the Burgess Shale at the Marble Canyon fossil locality (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). Discovered in 2012, Marble Canyon is the second most intensively sampled Burgess Shale locality after the Walcott Quarry, with major quarrying operations led by the Royal Ontario Museum in 2014 and 2016. Marble Canyon has already revealed key insights into a number of major animal groups such as Arthropoda (Aria and Caron Reference Aria and Caron2015, Reference Aria and Caron2017, Reference Aria and Caron2019; Aria et al. Reference Aria, Caron and Gaines2015; Moysiuk and Caron Reference Moysiuk and Caron2019) Chordata (Conway Morris and Caron Reference Conway Morris and Caron2014), Lophophorata (Moysiuk et al. Reference Moysiuk, Smith and Caron2017), Hemichordata (Nanglu et al. Reference Nanglu, Caron, Conway Morris and Cameron2016), and Annelida (Nanglu and Caron Reference Nanglu and Caron2018).
We then integrate these data with a census of fossils collected by the Royal Ontario Museum in the 1990s from a roughly 2-m-thick interval from the Raymond Quarry (Fletcher and Collins Reference Fletcher and Collins1998). With the exception of an unpublished M.S. thesis (Devereux Reference Devereux2001), the Raymond Quarry, which is one of the most extensively sampled Burgess Shale sites, has yet to be fully described in a quantitative paleocommunity framework. The collecting strategy employed at both Marble Canyon and Raymond Quarry—centimetric stratigraphic sampling—allows for uncommonly high-resolution consideration of temporal change within the paleocommunity, a scale broadly comparable to the Walcott Quarry. Integrating the Marble Canyon, the Walcott Quarry (Caron and Jackson Reference Caron and Jackson2008; O'Brien and Caron Reference O'Brien and Caron2015), the Tulip Beds (O'Brien and Caron Reference O'Brien and Caron2015), and the Raymond Quarry (Fig. 1A), represents the largest Cambrian Burgess Shale–type paleocommunity analysis to date.
Geologic Setting, Temporal Scales, and Taphonomic Considerations
The Marble Canyon fossil site (MC) is located in Kootenay National Park, British Columbia, Canada, approximately 40 km southeast of the Walcott Quarry in Yoho National Park, and is inferred to be the youngest Burgess Shale paleocommunity based on biostratigraphic evidence and stratigraphic position (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). The Marble Canyon quarry is located near the top of the “thick” Stephen Formation (Fig. 1A), immediately underneath the contact with the overlying Eldon Formation. This high position within the “thick” Stephen Formation is in sharp contrast to other Burgess Shale localities at Fossil Ridge, on Mount Stephen, and on Odaray Mountain, which lie in the lower and middle portions of the formation (Collins et al. Reference Collins, Briggs and Morris1983; Fletcher and Collins Reference Fletcher and Collins1998, Reference Fletcher and Collins2003; O'Brien et al. Reference O'Brien, Caron and Gaines2014).
Although the Cathedral Escarpment is not exposed in the Marble Canyon area, as it is in the type area (Aitken Reference Aitken1997; Fletcher and Collins Reference Fletcher and Collins1998, Reference Fletcher and Collins2003), several factors indicate that the “thick” Stephen Formation at Marble Canyon was deposited on a sloping surface immediately offshore of a pronounced break in submarine topography. The “platformal” or “thin” Stephen Formation (e.g., Aitken Reference Aitken1997), which is immediately adjacent to the Marble Canyon locality (<300 m), is organized into a series of seven depositional parasequences, reflecting the control of sea-level variations over its stratigraphic expression (Caron et al. Reference Caron, Gaines, Mángano, Streng and Daley2010; Gaines Reference Gaines2014). No such depositional cyclicity is evident in the “thick” or “basinal” Stephen Formation at Marble Canyon, indicating deposition in a comparably deep-water setting that was not influenced by changes in sea level (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). Furthermore, the section at Marble Canyon is characterized by prevalence of slumps, slide masses, and other manifestations of soft-sediment deformation (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). Together, these factors point to a sharp slope break, which may have taken the form of a near-vertical escarpment, as in the type area, or an oversteepened slope that descended to the basin. Regardless, these features demonstrate that an overall similar biostratinomic setting was shared among all Burgess Shale sites at Marble Canyon and in Yoho National Park. Fossil assemblages were rapidly entombed in fine-grained claystones in settings that lay near the foot of the Cathedral Escarpment, or in the case of Marble Canyon, at its distal sloping surface, in relatively close proximity to their living environments (Caron and Jackson Reference Caron and Jackson2006; Gaines Reference Gaines2014).
At Marble Canyon and the three other localities considered here, fossils were preserved via a globally distributed taphonomic pathway referred to as “Burgess Shale–type preservation” (Butterfield Reference Butterfield1995; Gaines Reference Gaines2014), by which whole assemblages of soft-bodied organisms were preserved as primary organic remains. This mode of fossilization, which resulted from early suppression of microbial decomposition in the sediments of the early burial environment (Gaines et al. Reference Gaines, Hammarlund, Hou, Qi, Gabbott, Zhao, Peng and Canfield2012), is widespread in early Phanerozoic strata and is much less common thereafter (Butterfield Reference Butterfield1995), a pattern that has been linked to changes in seawater chemistry (Gaines et al. Reference Gaines, Hammarlund, Hou, Qi, Gabbott, Zhao, Peng and Canfield2012). The four assemblages considered here reflect highly consistent biostratinomic and early diagenetic conditions. In addition, the presence of abundant and pristinely preserved vermiform enteropneusts (Oesia disjuncta) across most of the Marble Canyon quarry interval suggests little pre- or postmortem transport or decay (Nanglu et al. Reference Nanglu, Caron and Cameron2015; Beli et al. Reference Beli, Piraino and Cameron2017), similar to Walcott Quarry (Caron and Jackson Reference Caron and Jackson2008), the Tulip Beds (O'Brien et al. Reference O'Brien, Caron and Gaines2014), and likely the Raymond Quarry, based on disparate suites of soft-bodied organisms of low preservation potential.
Although the “thick” Stephen Formation at Marble Canyon exhibits major differences in thickness and stratigraphy from those localities of the type area in Yoho National Park (Fletcher and Collins Reference Fletcher and Collins1998; Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014), the upper and lower contacts of the “thick” Stephen Formation provide excellent markers of regional environmental changes (e.g., Aitken Reference Aitken1997) and allow for ordering of the four localities in relative time, with Marble Canyon being the youngest site, followed by Raymond Quarry, Walcott Quarry, and Tulip Beds (Fig. 1A). Biostratigraphically, the “thick” Stephen Formation lies between the overlying trilobite Ptychagnostus gibbus Zone and the underlying Glossopleura Zone (Sundberg Reference Sundberg1994). According to published time-calibrated biostratigraphic schemes, this interval could potentially represent less than 200,000 yr (Peng et al. Reference Peng, Babcock and Cooper2012). At Marble Canyon, the mudstones of the “thick” Stephen Formation correspond to approximately 35 m of strata (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014), yielding a very crude time estimate of ca. 1.75 m/10 kyr, assuming that sediment accumulation rates were relatively constant at the meter scale. However, the 4-m-thick Marble Canyon quarry lies at the top of the “thick” Stephen Formation within the upper part of the Ehmaniella burgessensis faunule of the Ehmaniella biozone (Sundberg Reference Sundberg1994; Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). In this quarry interval (and all the other quarry intervals studied herein), the presence of several centimeter scale–event deposited claystones (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014), suggests a shorter timescale than the one calculated based on constant sedimentation rates (i.e., at Marble Canyon <24,000 yr). In the Fossil Ridge area, the “thick” Stephen Formation encompasses 140 m of strata (Fletcher and Collins Reference Fletcher and Collins1998), yielding a comparably crude accumulation rate of 7 m/10 kyr, making the same assumptions as above. This approach yields a maximum estimate of ca. 10,000 yr for the deposition of the 7-m “Greater Phyllopod Bed,” which is consistent with the lower boundary estimate already calculated for this interval (e.g., Caron and Jackson Reference Caron and Jackson2008). With the same accumulation rate, the 2-m-thick interval of the Raymond Quarry studied herein would roughly represent 3000 yr, and the interval separating the Walcott and Raymond quarries (approximately 20 m; Fletcher and Collins Reference Fletcher and Collins1998), ca. 30,000 yr. The span of time between the Tulip Beds and the Walcott Quarry or the Raymond Quarry and Marble Canyon are less constrained owing to lateral stratigraphic variations between outcrops. It is important to underscore that these estimates are very poorly constrained, but nevertheless they are conservative and suggest that each assemblage represents no more than a few thousands of years and that assemblages are separated by intervals of a few tens of thousands of years.
Materials and Methods
Field collections at Marble Canyon (MC) were conducted in 2012, 2014, and 2016, yielding a total of 21,687 individual specimens (field observations and collected specimens) representing 73 species. The contact point between the Stephen Formation and the overlying Eldon Formation was used as a reference point for stratigraphic measurements, with specimens measured in centimeters below this contact (Fig. 1B). Stratigraphic level was recorded for all specimens included in subsequent quantitative analyses. Due to prominent lateral variation in bedding throughout the quarry interval, fossils were organized into 10-cm-thick stratigraphic subunits henceforth referred to as sample intervals (SI), which represent multiple discrete burial events. SIs are composed of tens to roughly 100 claystone beds, ranging from 1 to 37 mm in thickness. Thus the sum of all fossils collected within a 10-cm-thick SI, typically representing dozens of burial events, constitutes an induced time-average assemblage of roughly 600 yr in duration (see previous section).
SIs with a total specimen count under 299 were removed to acknowledge possible undersampling effects; a comparable threshold number was used for the Walcott Quarry in previous analyses (Caron and Jackson Reference Caron and Jackson2008). The resultant species abundance matrix included 18 SIs from Marble Canyon, the uppermost of which ends 230 cm below the Eldon–Stephen contact (Fig. 1C). In total, 16,438 specimens were included in all analyses of faunal composition, diversity, and turnover at MC.
Specimens were identified to the species level whenever possible, and abundance was tabulated using complete censuses of slab surfaces following methodologies used in previous quantitative studies of the Burgess Shale (Caron and Jackson Reference Caron and Jackson2006, Reference Caron and Jackson2008; O'Brien and Caron Reference O'Brien and Caron2015). In particular, some specimens were identified based on disarticulated or dissociated material, including carapaces of bivalved arthropods such as Tokummia katalepsis, hurdiid carapaces such as Hurdia sp., or the myomeres of Metaspriggina sp. In rare instances, the number of specimens from surfaces composed exclusively of shelly taxa such as Haplophrentis carinatus or Ptychagnostus sp., with individuals typically overlapping one another, were estimated using the density of a 2 cm2 square and extrapolating to the total surface area covered by those specimens.
Other Burgess Shale Localities
The Walcott Quarry (WQ) is located on Fossil Ridge between Mount Wapta and Mount Field and is part of the Walcott Quarry Member of the “thick” Stephen Formation (Fletcher and Collins Reference Fletcher and Collins1998; Caron and Jackson Reference Caron and Jackson2006). We used the data set originally published by Caron and Jackson (Reference Caron and Jackson2008) and updated by O'Brien and Caron (Reference O'Brien and Caron2015) for our quantitative analyses. This data set comprised 43,296 specimens, divided stratigraphically into 26 SIs with more than 300 collected specimens each.
The Raymond Quarry (RQ) is located approximately 20 m above the base of the Greater Phyllopod Bed of the Walcott Quarry (Raymond Quarry Member; Fletcher and Collins Reference Fletcher and Collins1998, Reference Fletcher and Collins2003). A total of 9539 observations were made systematically from a 2.3-m-thick stratigraphic section using the same methodology as for the MC, using 10 cm bins to create SIs roughly comparable with those of the MC. We identified 114 taxa from this relatively restricted stratigraphic section of the RQ, with dominant taxa being the arthropod Leanchoilia, the sponge Choia, and the priapulid Ottoia.
The Tulip Beds (TB) locality on Mount Stephen, which occurs within the Campsite Cliff Shale Member of the “thick” Stephen Formation, was described in O'Brien et al. (Reference O'Brien, Caron and Gaines2014) and O'Brien and Caron (Reference O'Brien and Caron2015). Although the middle and upper parts of the “thick” Stephen Formation are not present at the TB locality, extrapolation from Mount Stephen to the nearby (<2 km) complete section on Mount Field reveals that the Campsite Cliff Shale Member lies below the Walcott Quarry Member (Fletcher and Collins Reference Fletcher and Collins1998, Reference Fletcher and Collins2003). This locality shares roughly 50% of its genera with the Walcott Quarry and is largely dominated by an indeterminate egg-shaped taxon, the tulip animal Siphusauctum gregarium, the alga Marpolia spissa, and the predator Anomalocaris canadensis. A total of 7906 specimens (the data set from O'Brien and Caron [Reference O'Brien and Caron2015]) from the Tulip Beds were included in our analyses. These specimens were analyzed as a bulk assemblage (O'Brien et al. Reference O'Brien, Caron and Gaines2014).
With all four localities combined, our ecological data matrix includes abundance patterns for 77,179 specimens of 234 taxa (Supplementary Tables 1–3).
Diversity Patterns and Multivariate Analyses
Rarefaction curves were produced in PAST. All other quantitative data was prepared in Microsoft Excel and analyzed in R using the following packages: vegan, MASS, analogue, dendextend, dplyr, indicspecies, and ggplot2.
Rarefaction curves were used to determine the extent to which the MC locality had been adequately sampled to recover the majority of the true total species richness. Whittaker plots (also called rank abundance curves) were also plotted for the three richest and poorest SIs to observe how paleocommunity structure may be influenced by patterns of evenness/dominance of species abundances.
Cluster analysis (performed in R using the function hclust with complete linkage) was used to identify major groups of subunits sharing similar species or ecological compositions. We used the Morisita-Horn overlap index to generate the ecological distance matrices, which is a common metric that has been demonstrated to be robust to variation in sample sizes (Magurran Reference Magurran2004; Barwell et al. Reference Barwell, Isaac and Kunin2015).
Correspondence analysis (CA), a method of ordination (Gauch Reference Gauch1982; Legendre and Legendre Reference Legendre and Legendre2012) was chosen to visualize the scores of species and sampling intervals into two dimensions. CA is a commonly used method for exploring paleocommunity data sets (Hammer and Harper Reference Hammer and Harper2006), including previous studies of the Burgess Shale (Caron and Jackson Reference Caron and Jackson2006, Reference Caron and Jackson2008) and the Chengjiang biota (Zhao et al. Reference Zhao, Caron, Bottjer, Hu, Yin and Zhu2013). CA does not make presuppositions regarding the presence or absence of an environmental gradient in the data set, although this method expects a unimodal distribution of the species (Gauch Reference Gauch1982; Legendre and Legendre Reference Legendre and Legendre2012). Eigenvalues returned by CA, while less straightforward to interpret than in an ordination technique such as principal components analysis, represent a proportion of the variance among the data set explained by their respective axes (Legendre and Legendre Reference Legendre and Legendre2012). CA of the Burgess Shale–wide data set showed an arch effect, which is a common artifact resulting from the driving factors influencing the plotting of the first correspondence axis also affecting the second axis. To accommodate for this effect, detrended correspondence analysis (DCA) was used, which is a method specifically designed to remove distorting properties of arch effects (Legendre and Legendre Reference Legendre and Legendre2012).
When performing CA and DCA, the raw data matrix of species abundances by site was first standardized by column and row weights, which are calculated by dividing the sum of all the values in the respective column or row by the sum of the values of the entire table (Legendre and Legendre Reference Legendre and Legendre2012; Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, Mcglinn, Minchin, Hara, Simpson, Solymos, Stevens, Szoecs and Wagner2019). In both CA and DCA, (1) sites similar in species composition should appear close to each other along at least one axis; (2) species should appear close to sites in which they are abundant; (3) species should appear close to other species with similar distributions across sites; and (4) species that are abundant across multiple sites should appear toward the center of the ordination, while those that are rare across most sites are plotted at the edges (Hammer and Harper Reference Hammer and Harper2006; Legendre and Legendre Reference Legendre and Legendre2012).
Finally, species indicator analysis was used to identify species that may be uniquely characteristic of a given locality. By providing a priori knowledge of how particular sites are grouped, this method calculates an indicator statistic for each species in the paleocommunity data matrix that takes into account: (1) specificity, that is, the likelihood of the sampled site (in this case, SI) belonging to a target site group (in this case, a locality of WQ, TB, RQ, or MC) given that the species has been found; and (2) fidelity, that is, the probability of finding the species at any sites belonging to the target site group (De Cáceres et al. Reference De Cáceres, Legendre, Wiser and Brotons2012). An indicator taxon is one that is, therefore, both highly specific to and relatively ubiquitous at a given locality or combination of localities.
Ecological Groups and Functional Diversity
The ecological mode was described using a modified version of a previously published three-axis system (Bambach et al. Reference Bambach, Bush and Erwin2007). This system posits that the ecological mode of the vast majority of marine taxa can be described using three variables: (1) the vertical position within the water column they occupied, (2) their motility, and (3) their trophic strategy. This framework has been used in other analyses of Cambrian fauna (Conway Morris Reference Conway Morris1986; Caron and Jackson Reference Caron and Jackson2008; Zhao et al. Reference Zhao, Caron, Bottjer, Hu, Yin and Zhu2013; O'Brien and Caron Reference O'Brien and Caron2015) and is particularly useful when considering animal communities that incorporate a wide variety of phyla. While admittedly overgeneralized in some cases, this approach can accommodate the wide variety of ecologies employed by the considerable disparity of body plans represented by a well-preserved marine locality. We then constructed an ecological-group matrix whereby each taxon was assigned an ecological mode using the three-axis system. This data set was then analyzed using cluster analyses, CA and DCA (as described earlier).
Results
Diversity Patterns at Marble Canyon
In total, 73 species-level taxa are present at MC, most coming from discrete bedding intervals appearing black in outcrop (Fig. 1B,C). Abundant species tend to recur across the studied interval (e.g., Liangshanella, Haplophrentis, Oesia, and Peronopsis), while others have patchier temporal distributions (e.g., Yawunik; Fig. 1D). Kootenayscolex and Primicaris seem particularly abundant in SIs characterized by thicker beds (between SI 420 and SI 350; Fig. 1D). As is the case in all major studied Burgess Shale sites (Conway Morris Reference Conway Morris1986; Caron and Jackson Reference Caron and Jackson2008; O'Brien and Caron Reference O'Brien and Caron2012), arthropods are the most abundant and diverse major taxonomic group (Fig. 1D). However, in no other Burgess Shale sites are hemichordates, annelids, or chordates such significant components of the paleocommunity, at 16%, 5%, and 1.3% by abundance, respectively. When considering diversity irrespective of abundance, we found a substantial number of radiodonts, which constitute roughly 11% of the total species richness, many of which are currently undescribed forms (Fig. 1D).
The slope of the rarefaction curve for the Marble Canyon as a bulk assemblage, compared with the Walcott Quarry, Raymond Quarry, and Tulip Beds bulk assemblages, suggests that while all localities have been sufficiently sampled to capture most of the species richness, Marble Canyon has the greatest potential (the curve is furthest from reaching an asymptote) for new species discovery with continued sampling (Fig. 2A). SIs 350, 380, 400, and 410 have the highest observed species richness, while SIs 300, 310, and 320 have the lowest (Figs. 2B, 3). The original description of the Marble Canyon assemblage suggested a high diversity of arthropods comparable with Walcott Quarry (Caron et al. Reference Caron, Gaines, Aria, Mángano and Streng2014). This present study confirms this trend with a much larger sample size (Supplementary Fig. 1), although other taxonomic groups are not as diverse (i.e., sponges), which explains the relatively lower diversity level overall.
Rank abundance curves illustrate two main types of assemblages, the most extreme examples of which are illustrated in Figure 4. Most SIs have relatively steep curves, indicating that their respective temporal subcommunities were numerically dominated by relatively few taxa (Fig. 4). In most cases, these were the four most abundant and stratigraphically wide-ranging taxa, namely Haplophrentis and Oesia, which are stratigraphically widespread but most dominant between SI 330 and SI 230, and Liangshanella and Peronopsis, particularly below SI 330. SIs with less steep rank abundance curves are indicative of assemblages that are less dominated by a small group of taxa. In addition to having the greatest evenness of species composition, SIs 380, 410, and 350 are also among the most species-rich subunits (Figs. 3, 4).
All SIs at Marble Canyon are dominated by one of three groups: Arthropoda, Hemichordata, or Lophophorata (Fig. 5). Hemichordate and lophophorate dominances are almost entirely due to Oesia and Haplophrentis, respectively. The abundance of other species of hemichordate and lophophorates make up only 1.5% of the paleocommunity (0.2% and 1.3%, respectively). Together, these three dominant groups (Arthropoda, Hemichordata, and Lophophorata) constitute between 81% (SI 390) and 99% (SI 230 and SI 500) of the total species abundance, with the upper levels of the quarry (SI 320 and above) dominated by hemichordates and lophophorates and the lower levels dominated by arthropods. This change is due to both the increase in abundance of large predatory arthropods (e.g., Yawunik and Sidneyia) and, more significantly, the increase in abundance of the arthropod Liangshanella, the most abundant taxon in SI 420 to SI 350.
Taxonomic cluster analysis of the Marble Canyon distinguishes two principal clusters (Fig. 6) that we informally call the upper and lower quarry intervals. The lower quarry (LQ) consists of SIs 500 to 480. The upper quarry (UQ) consists of all other SIs. The UQ is further subdivided stratigraphically into the uppermost levels of the UQ (SIs 340 to 230; UQ1) and the lower levels of the UQ (SIs 420 to 350; UQ2).
The CA corroborates the division of SIs detected by cluster analysis (Fig. 7). The first two CA axes represent roughly 67% of the variation among sites. Axis 1 largely defines two major types of assemblages (LQ vs. UQ) and likely represents an environmental factor fluctuating primarily as a function of time. Ptychagnostus and Spartobranchus are strongly associated with the LQ SIs. Axis 2 further discriminates between the UQ1 and UQ2 SIs. Taxa that are highly recurrent throughout the quarry are found close to the center of this axis, such as Yawunik, Sidneyia, Metaspriggina, and Oesia (Fig. 7). Taxa that are abundant but have more stratigraphically structured distributions such as Haplophrentis and Kootenayscolex are plotted farther from the center of the ordination.
Patterns of Ecological Change at Marble Canyon
In SI 370 and from SI 350 to SI 230, epibenthic-sessile-suspension feeding (ESSU) is the most common ecological mode, generally constituting greater than 50% of the entire ecological mode occupancy (Fig. 8). In SI 360 as well as SI 420 to SI 380, nektonic-vagrant-hunting/scavenging (NKHS) is the most common ecological mode. In SI 500 to SI 480, nektobenthic-vagrant-deposit feeding (NKDE) is the most common ecological mode. Together, these three ecological modes (ESSU, NKHS, and NKDE) constitute at least 79% of the total ecological mode occupancy of their respective SIs (SI 390) and at most 97% (SI 500).
The results of cluster analysis and CA recapitulate the SI-level quantitative patterns of dominance described earlier (Figs. 9, 10). Ecological group 1 (EG1) consists of the SIs dominated by the ESSU ecological mode; ecological group 2 (EG2) consists of the SIs dominated by the NKHS ecological mode; and ecological group 3 (EG3) consists of the SIs dominated by the NKDE ecological mode (Fig. 9). These clusters were used to build the convex hulls in the CA of the Marble Canyon SI-level ecological compositions (Fig. 10). Compared with the other modes, the NKDE-dominated SIs (EG3) are far from the center of the ordination, reflecting a clear shift in ecological modes between older (EG3) and younger (EG1 and EG2) intervals.
Overall, we observe long periods of compositional stasis interrupted twice by rapid changes in species assemblages: once from the LQ (below SI 420) fauna to UQ2, and again from UQ2 to UQ1 (approximately between SI 350 and SI 340). These patterns are closely mirrored using both ecological and taxonomic groups, suggesting that the factors that caused major changes in faunal composition similarly caused a change in the environmental conditions of Marble Canyon, favoring a new faunal assemblage with a new ecological structure than what had existed previously.
Patterns of Ecological Change across the “Thick” Stephen Formation
Intuitively, most SIs cluster with SIs from the same locality. However, the UQ2 SIs from MC cluster with a subset of the WQ sites, and the LQ SIs cluster with the TB (Fig. 11). The clustering of UQ2 MC and some WQ sites is likely due to the great abundance of Liangshanella in these SIs; the clustering of LQ SIs and TB is explained by the presence of Peronopsis, a statistically significant indicator species for an MC + TB-type SI (Table 1).
The first two DCA axes represent 67% of the total variation among sites (Fig. 12). Four non-overlapping groups of SIs represent MC, WQ, RQ, and TB (Fig. 12), each locality also having unique indicator species, which taken together suggest significant compositional differences between them (Table 1). These groups are easily distinguished on axis 1 (which summarizes 39% of the variability in the data set), suggesting that the WQ SIs are most representative of the Burgess Shale fauna overall, as they cluster near the center of the ordination. This position might in part reflect the fact that WQ represents the most heavily sampled site with the greatest number of specimens collected. On axis 1, MC and TB are most distant from each other, suggesting these faunas are the most highly differentiated within our data set. Being both the most stratigraphically and geographically disparate pair of localities (MC and TB), axis 1 may then represent a combination of environmental gradients operating at sufficiently large geographic and temporal scales such that the intermediate localities (WQ and RQ) cannot be clearly discriminated from each other on this axis.
Ecologically, four major clusters describe the SIs of our four-locality Burgess Shale data set (Fig. 13); however, each includes stratigraphically disparate SIs. These clusters can be defined by the dominant trophic mode present: the largest cluster represents suspension feeding–dominated SIs and contains SIs from all four localities, including the TB bulk assemblage; the second-largest cluster represents SIs dominated by hunter/scavenger taxa and includes SIs from WQ, RQ, and MC; the last two smaller clusters are dominated by deposit feeders and include SIs from WQ and MC (specifically, the LQ SIs from MC). When the results of the cluster analysis are used to plot convex hulls on the CA of the Burgess Shale ecological matrix, there is minimal overlap among these clusters (Fig. 14). The first two CA axes represent greater than 52% of the total variation, suggesting that the relative proportions of represented trophic modes is of principal importance when defining ecological structure among SIs.
Discussion
Short Timescale Turnover Patterns of the Burgess Shale Paleocommunity
Broadly, we can expect the compositional change of ecological communities to conform to one of two sets of patterns. The first pattern is that species or individuals are distributed entirely randomly; this view can be described as the null model of paleocommunity assembly, with the unified neutral theory of biodiversity being its most influential incarnation (Hubbell Reference Hubbell2005; Rosindell et al. Reference Rosindell, Hubbell and Etienne2011). Alternatively, paleocommunity composition and change may be structured by a combination of both abiotic and biotic factors, resulting in nonrandom patterns of species distributions over time. There is still considerable debate regarding the relative importance of these factors on community assembly, part of which is due to the temporal constraints of neontological data sets. This is particularly true of experimental studies analyzing trends in diversity change (Simberloff and Wilson Reference Simberloff and Wilson1969; Findlay and Kasian Reference Findlay and Kasian1986; Condit et al. Reference Condit, Ashton, Manokaran, LaFrankie, Hubbell and Foster1999), but even observational studies rarely span more than a few decades in length (Grant and Grant Reference Grant and Grant2002; Azaele et al. Reference Azaele, Pigolotti, Banavar and Maritan2006). Paleontological data sets are unique in this regard, providing a much broader temporal sampling than is otherwise possible.
Our data from MC suggest that even at the level of 10-cm-thick SIs, an uncommonly short timescale for a paleontological data set (ca. 600 yr), the fauna alternated between periods of relative stability and rapid species turnover (Figs. 5, 6). At Walcott Quarry, this was taken as an indication of periodic environmental disturbance followed by rapid recolonization (Caron and Jackson Reference Caron and Jackson2008), often resulting in a very similar fauna between successive SIs. Additionally, most SIs of both quarries are numerically dominated by taxa with broad stratigraphic distributions (Fig. 1). At Marble Canyon, these taxa are Oesia, Haplophrentis, and Liangshanella; at Walcott Quarry, Liangshanella, Hazelia, and Marrella are the most temporally recurrent taxa. These taxa also have wide geographic ranges, further suggesting that long periods of compositional stasis were maintained by a subset of taxa with long survivorship, likely due to high dispersal and/or colonization ability (Foote Reference Foote2003; Payne and Finnegan Reference Payne and Finnegan2007; Hopkins Reference Hopkins2011).
The nature of the perturbations that might have resulted in major compositional changes at Marble Canyon as well as among other Burgess Shale sites (see following section) must remain speculative. A variety of abiotic factors are known to result in major changes in modern marine community diversity and structure (Sousa Reference Sousa1984; Conversi et al. Reference Conversi, Dakos, Gårdmark, Ling, Folke, Mumby, Greene, Edwards, Blenckner, Casini, Pershing and Möllmann2015) such as shifting temperatures and perturbations of ocean chemistry at ecologically relevant timescales (Harley et al. Reference Harley, Hughes, Hultgren, Miner, Sorte, Thornber, Rodriguez, Tomanek and Williams2006; Bornette and Puijalon Reference Bornette and Puijalon2011). The near-total absence of sponges from Marble Canyon could be due to either turbid conditions impairing their ability to feed (Pineda et al. Reference Pineda, Strehlow, Sternel, Duckworth, Den Haan, Jones and Webster2017) or to a lack of an appropriate substrate for larval recruitment (Maldonado and Young Reference Maldonado and Young1996; Whalan et al. Reference Whalan, Abdul Wahab, Sprungala, Poole and De Nys2015). While we are unable at this time to provide a mechanistic correlate for changes in faunal composition we observed, the most significant of these changes (after which community structure remained relatively stable for thousands of years) may have been caused by similarly broad-scale environmental fluctuations. Testing this hypothesis would require additional explanatory proxies, including the potential availability of geochemical data, in the future.
Macroecological Patterns of Paleocommunity Structure across the Burgess Shale
Major patterns of taxonomic compositional change observed within and between localities are at odds with a random distribution of taxa from a single species pool. One of the most striking examples of this temporal variation in species composition among localities is the abundance of sponges. They are the second most diverse phylum at the Tulip Beds and Walcott Quarry, and the second most abundant phylum at the Walcott Quarry and Raymond Quarry, but at Marble Canyon, only 5 specimens out of 16,438 are sponges. The priapulids are similarly abundant throughout most well-studied Cambrian localities, particularly at the Raymond Quarry, where Ottoia alone makes up 10.1% of specimens from our sampled interval. At Marble Canyon, however, only a few specimens of Selkirkia and an as-of-yet indeterminate priapulid represent this phylum. Other taxa show a systematic increase in abundance from older to younger localities. While hemichordates are absent from the Tulip Beds, each of the younger localities has a different hemichordate taxon that is prevalent enough to be an indicator species (Spartobranchus, an undescribed enteropneust, and Oesia for Walcott Quarry, Raymond Quarry, and Marble Canyon, respectively). The only annelid at the Tulip Beds is a single specimen of Burgessochaeta. By contrast, the younger Walcott Quarry and Marble Canyon possess significant numbers of annelids (Burgessochaeta and Kootenayscolex, respectively).
A number of possibilities exist to explain the high degree of taxonomic heterogeneity among our sampled localities. The first are ecological interactions mediating species distributions. Interspecific competition leading to competitive exclusion may have sequentially excluded species from younger localities (Hardin Reference Hardin1960). For example, sponges may have been systematically outcompeted in the epifaunal suspension-feeding niche by taxa such as Haplophrentis (Moysiuk et al. Reference Moysiuk, Smith and Caron2017) and Oesia (Nanglu et al. Reference Nanglu, Caron, Conway Morris and Cameron2016). The near-complete absence of Burgessochaeta may be due to its niche being filled by Kootenayscolex, which is thought to have shared a similar mode of life (Caron and Jackson Reference Caron and Jackson2008; Nanglu and Caron Reference Nanglu and Caron2018). The same could be argued of other taxa that presumably shared a significant degree of ecological overlap and show a similar pattern of distributional turnover, such as Yawunik (Aria et al. Reference Aria, Caron and Gaines2015) and Tokummia (Aria and Caron Reference Aria and Caron2017) at Marble Canyon, possibly filling the ecospace of Leanchoilia and Branchiocaris from older localities. However, biotic drivers of community assembly, such as competition (Connell Reference Connell1961a,b; McCook et al. Reference McCook, Jompa and Diaz-Pulido2001) or dispersal ability (Rosindell et al. Reference Rosindell, Hubbell and Etienne2011; Klompmaker and Finnegan Reference Klompmaker and Finnegan2018), are difficult to test in a paleontological setting and thus remain speculative.
Our view of Cambrian trophic ecology is also impacted through our quantitative SI-level analysis of multiple localities. The view that Cambrian paleocommunities were relatively stable in ecological structure has been largely based on limited comparisons between sites at the level of broadly sampled bulk assemblages without quantitative comparisons of how within-locality ecological fluctuations compared across localities (Zhao et al. Reference Zhao, Caron, Bottjer, Hu, Yin and Zhu2013; O'Brien and Caron Reference O'Brien and Caron2015). For example, a bulk assemblage comparison of the Tulip Beds and Walcott Quarry indicated that there was a broadly similar diversity of ecological strategies (and that this pattern was also shared with the Chengjiang localities; O'Brien and Caron Reference O'Brien and Caron2015). Our data suggest that if we consider fine-scale stratigraphic variation rather than combining the ecologies of each major locality into a single pool, the dominant ecological modes were susceptible to considerable variation. SIs may be more similar in ecological structure to temporally disjunct SIs than to their stratigraphically adjacent neighbors. In this context, ecological structure throughout the Burgess Shale was highly heterogeneous both at the level of 10 cm assemblages and localities.
Predation is often suggested as one of the principal forces structuring communities during the Cambrian (Bengtson Reference Bengtson2002; Hu et al. Reference Hu, Steiner, Zhu, Erdtmann, Luo, Chen and Weber2007; Vannier et al. Reference Vannier, Steiner, Renvoisé, Hu and Casanova2007; Vannier Reference Vannier2012; Laflamme et al. Reference Laflamme, Darroch, Tweedt, Peterson and Erwin2013), as well as one of the factors most distinguishing Cambrian communities from the preceding Ediacaran fauna (Laflamme and Narbonne Reference Laflamme and Narbonne2008; Xiao and Laflamme Reference Xiao and Laflamme2009). In our analyses, the epifaunal suspension feeding–dominated SIs are the largest cluster, suggesting that this mode of life was at least as important as predation for defining paleocommunity ecological structure. One way of interpreting these data would be to suggest periodic alternation between top-down (predator and competition-mediated) versus bottom-up (prey and environment-mediated) control of community ecology (Roff et al. Reference Roff, Doropoulos, Rogers, Bozec, Krueck, Aurellado, Priest, Birrell and Mumby2016). However, it has been noted that this dichotomy may not be appropriate in marine settings, as these mechanisms are not mutually exclusive (Ruppert et al. Reference Ruppert, Travers, Smith, Fortin and Meekan2013; Conversi et al. Reference Conversi, Dakos, Gårdmark, Ling, Folke, Mumby, Greene, Edwards, Blenckner, Casini, Pershing and Möllmann2015). For example, many predatory species have larval forms (which we are unlikely to see in the fossil record, except in the most exceptional of circumstances; Maas et al. Reference Maas, Braun, Dong, Donoghue, Müller, Olempska, Repetski, Siveter, Stein and Waloszek2006; Yin et al. Reference Yin, Zhu, Knoll, Yuan, Zhang and Hu2007) that are highly sensitive to environmental conditions, while sessile organisms are frequently in direct competition for substrate (McCook et al. Reference McCook, Jompa and Diaz-Pulido2001). It is, therefore, most appropriate to view the structure of Cambrian paleocommunities as we now view modern benthic communities: the result of multiple interspecific drivers such as predation (Connell Reference Connell1961a), competition (Connell Reference Connell1961b; Paine Reference Paine1974; McCook et al. Reference McCook, Jompa and Diaz-Pulido2001; Hixon and Jones Reference Hixon and Jones2005), niche partitioning (MacArthur Reference MacArthur1958), and dispersal (Hubbell Reference Hubbell2001) acting simultaneously on the backdrop of abiotic factors such as environmental perturbation (Krebs Reference Krebs2009; Conversi et al. Reference Conversi, Dakos, Gårdmark, Ling, Folke, Mumby, Greene, Edwards, Blenckner, Casini, Pershing and Möllmann2015; Pershing et al. Reference Pershing, Mills, Record, Stamieszkin, Wurtzell, Byron, Fitzpatrick, Golet and Koob2015).
Conclusions
While our understanding of the Cambrian explosion continues to progress, the discovery of the Marble Canyon fossil site clearly demonstrates that major discoveries at the level of entire faunas remain to be made. Such discoveries have the potential to shed light not only on the early origins of nearly all major metazoan phyla, but also on the structure and diversity of Cambrian ecosystems. Our study represents the first quantitative analysis of the four best-sampled Burgess Shale localities and the largest analysis of a Cambrian Burgess Shale–type paleocommunity to date. The integration of fine-resolution stratigraphic data possesses significant promise for disentangling the potential mechanisms behind the assembly of some of the earliest complex animal paleocommunities. SI-level changes in both taxonomic and ecological structure show that fluctuations between periods of relative stasis and rapid compositional change are characteristic among Cambrian localities at short timescales. Widening the scope of our observations to stratigraphically disparate localities reveals that the dominance of different ecological modes was patchy throughout the Burgess Shale and that major taxonomic changes occurred within and between localities. These data indicate that the Burgess Shale was highly heterogeneous with respect to ecology and species composition at both ecological and geologic timescales.
Acknowledgments
We thank P. Fenton and M. Akrami for collections assistance at the Royal Ontario Museum, E. Welsh for help with data entry (Raymond Quarry), and T. M. Cullen and M. Orobko for discussions and help with R coding. H. Cyr, D. Currie, M. Aberhan, M. Patzkowsky, and two anonymous reviewers provided helpful comments on an earlier draft of this article. Material for this study was collected under several Parks Canada Research and Collections permits (to J.-B.C.). Major funding support for fieldwork comes from Barbara Polk Milstein and several Royal Ontario Museum grants (to J.-B.C.), the National Geographic Society (2014 research grant to J.-B.C.), the National Science Foundation (2011 EAR-1046233 and 2016 EAR-1556226 awards to R.R.G.–Pomona College). K.N.'s doctoral research was supported by fellowships from the University of Toronto (Department of Ecology and Evolutionary Biology) and J.-B.C.'s NSERC Discovery Grant (number 341944). This is Royal Ontario Museum Burgess Shale project number 81.