Introduction
Comparing vertebrate microfossil assemblages to assess marginal paleocommunity changes, including differences in taxon abundance among sites, may allow paleoecologists to test hypotheses about the drivers of macroevolution (cf. Jablonski and Sepkoski Reference Jablonski and Sepkoski1996); however, multiple biases confound direct site comparisons (Blob and Fiorillo Reference Blob and Fiorillo1996; Cutler et al. Reference Cutler, Behrensmeyer and Chapman1999; Gates et al. Reference Gates, Sampson, Zanno, Roberts, Eaton, Nydam, Hutchinson, Smith, Loewen and Getty2010). Vertebrate microfossil bonebeds (VMBs or “microsites”) record genuine paleocommunity signals, including species composition and relative abundance differences possibly related to paleoenvironmental, stratigraphic, and/or geographic differences among sites, but VMB assemblages also reflect different taphonomic biases (Eberth Reference Eberth1990; Rogers et al. Reference Rogers, Carrano, Curry Rogers, Perez and Regan2017). Assemblage specimen compositions may differ, in part, because different species have different transport and preservation potentials and/or because different amounts of spatial and temporal averaging result in variable specimen mixing at sites (Miller Reference Miller1988; Blob and Fiorillo Reference Blob and Fiorillo1996; Tomasovych and Kidwell Reference Tomasovych and Kidwell2009a, Reference Tomasovych and Kidwell2010; Bürkli and Wilson Reference Bürkli and Wilson2017; Rogers et al. Reference Rogers, Carrano, Curry Rogers, Perez and Regan2017). Fossil specimen concentrations in terrestrial formations can result from attritional, lag, and allochthonous depositions (Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000; Brinkman et al. Reference Brinkman, Russell, Peng, Currie and Koppelhus2005b; Rogers and Brady Reference Rogers and Brady2010; Rogers et al. Reference Rogers, Carrano, Curry Rogers, Perez and Regan2017). These processes facilitate rare species accumulation but often homogenize deposits, obscuring formation beta diversity (Kidwell and Flessa Reference Kidwell and Flessa1995; Tomasovych and Kidwell Reference Tomasovych and Kidwell2009a; Kidwell and Tomosovych Reference Kidwell and Tomasovych2013). This is especially true for tidally influenced channel deposits, which appear particularly conducive to VMB formation due to lag accumulation (Rogers and Kidwell Reference Rogers and Kidwell2000). Many VMBs yield both obligate terrestrial and obligate aquatic taxa (e.g., Carrano et al. Reference Carrano, Oreska and Lockwood2016; Cullen et al. Reference Cullen, Fanti, Capobianco, Ryan and Evans2016; Rogers et al. Reference Rogers, Carrano, Curry Rogers, Perez and Regan2017), and this postmortem mixing of specimens from different habitat types complicates the use of raw VMB specimen counts as paleocommunity proxies.
Significant habitat source mixing in fossil assemblages likely results primarily from temporal, rather than spatial, averaging (Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000; Kidwell and Holland Reference Kidwell and Holland2002). Time averaging can combine specimens in a given fossil assemblage from paleocommunities associated with different habitat types that occupied the paleo-catchment successively (Kidwell and Flessa Reference Kidwell and Flessa1995; Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000; Kidwell and Holland Reference Kidwell and Holland2002). Studies have measured time averaging in recent subfossil assemblages (cf. Kosnik et al. Reference Kosnik, Hua, Kaufman and Wüst2009; Kidwell Reference Kidwell2013), and a few studies have noted that the amount of time averaging in a given subfossil assemblage may increase with deposit age (Terry and Novak Reference Terry and Novak2015; Kowalewski et al. Reference Kowalewski, Casebolt, Hua, Whitacre, Kaufman and Kosnik2017). Terrestrial fossil assemblages are more likely to record habitat source mixing from time averaging than marine subfossil assemblages, but studies have not yet attempted to quantify temporal mixing in older, terrestrial deposits.
In comparison, spatial averaging may increase compositional heterogeneity among otherwise similar assemblages representing a given “paleocommunity type” (sensu Bennington and Bambach Reference Bennington and Bambach1996) but is unlikely to completely obscure differences among living communities contributing to deposits. Both allochthonous specimen transport (Miller and Cummings Reference Miller and Cummings1990; Bürkli and Wilson Reference Bürkli and Wilson2017) and differences in microhabitat density within the paleo-catchments surrounding depositional sites—for example, ephemeral ponds on an otherwise terrestrial landscape—give rise to spatial averaging. Recent death assemblages often show some degree of spatial mixing but preserve evidence of spatial and environmental gradients (Tomasovych and Kidwell Reference Tomasovych and Kidwell2009b; Albano and Sabelli Reference Albano and Sabelli2011; Tyler and Kowalewski Reference Tyler and Kowalewski2017). Most work has focused on marine invertebrate communities, but vertebrate death assemblages can also preserve real spatial distinctions (Miller Reference Miller2012). These patterns may be discernible in much older vertebrate fossil assemblages (Western and Behrensmeyer Reference Western and Behrensmeyer2009), provided that this paleoecologically significant information can be differentiated from habitat source mixing resulting from time averaging.
Studies have used various methods to constrain mixing due to time averaging in recent subfossil assemblages, with approaches ranging from absolute dating and taphonomic categorization to modeling faunal turnover patterns (e.g., Kosnik et al. Reference Kosnik, Hua, Kaufman and Wüst2009; Tomasovych and Kidwell Reference Tomasovych and Kidwell2010; Belanger Reference Belanger2011). However, applying these techniques to terrestrial vertebrate assemblages from deep time remains a challenge. In some cases, paleoecologists can quantify habitat mixing using the presumed life habits of preserved taxa, but this approach is tenuous for vertebrates with uncertain affinities in assemblages that lack modern analogues (e.g., Noto and Grossman Reference Noto and Grossman2010). Most Mesozoic terrestrial assemblages fall into this category, given the abundance of nonavian dinosaurs and uncertainty surrounding the life-habitats of other groups, for example, possibly amphibious crocodilians and mammals (e.g., Luo Reference Luo2007).
Many terrestrial paleocommunity studies, therefore, rely on exploratory approaches to compare assemblages, including cluster analysis (e.g., Miller and Cummings Reference Miller and Cummings1990; Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004; Cullen and Evans Reference Cullen and Evans2016; Araújo-Júnior et al. Reference Araújo-Júnior, Oliveira Porpino and Bergqvist2017) and ordination, often nonmetric multidimensional scaling (e.g., Zambito et al. Reference Zambito, Mitchell and Sheets2008; Moore Reference Moore2012; McMullen et al. Reference McMullen, Holland and O'Keefe2014), but these methods do not allow us to quantify the likelihood of paleoecological interpretations. Bayesian methods offer a promising alternative for comparing presumably biased, community-level samples (e.g., Ellison Reference Ellison2004; Clark Reference Clark2005; McCarthy Reference McCarthy2011). By combining prior information with new data and a model, Bayesian methods allow paleoecologists to quantify the likelihood of specific hypotheses, given available data (cf. McCarthy Reference McCarthy2011; Silvestro et al. Reference Silvestro, Schnitzler, Liow, Antonelli and Salamin2014; Pimiento et al. Reference Pimiento, Griffin, Clements, Silvestro, Varela, Uhen and Jaramillo2017). Although there are precedents for applying prior information to qualitatively infer paleoenvironment in vertebrate microfossil assemblages (e.g., Brinkman et al. Reference Brinkman, Braman, Neuman, Ralrick, Sato, Currie and Koppelhus2005a), Bayesian statistics have not been formally applied to investigate mixing bias in VMBs.
Bayesian mixing models developed to address source-mixing problems in other fields may provide a formal framework for quantifying mixing bias in microfossil assemblages. Bayesian mixing models are most commonly used to quantify consumer diet proportions using prior information about the stable isotope composition of diet constituents (e.g., Hopkins and Ferguson Reference Hopkins and Ferguson2012; Wilkinson et al. Reference Wilkinson, Carpenter, Cole, Pace and Yang2013; Matsubayashi et al. Reference Matsubayashi, Morimoto, Tayasu, Mano, Nakajima, Takahashi, Kobayashi and Nakamura2015); however, geochemists have adapted the approach to quantify different organic-matter source contributions to sediment carbon pools (e.g., Bergamino et al. Reference Bergamino, Dalu and Richoux2014; Larsen et al. Reference Larsen, Bach, Salvatteci, Wang, Andersen, Ventura and McCarthy2015; Craven et al. Reference Craven, Edwards and Flood2017) and to identify abiotic sediment contributions to suspended sediment loads and fluvial deposits (e.g., Koiter et al. Reference Koiter, Lobb, Owens, Petticrew, Tiessen and Li2013). In these examples, source pools are characterized using “tracer” distributions, which are often—but not exclusively—stable isotope ratios. In the latter application, commonly referred to as sediment “fingerprinting,” geochemists may use a variety of geochemical and radionuclide data to define unique source signatures (Koiter et al. Reference Koiter, Lobb, Owens, Petticrew, Tiessen and Li2013; Cooper et al. Reference Cooper, Krueger, Hiscock and Rawlins2014). The tracer composition of a mixed assemblage reflects the magnitude of the contribution from each “end-member” source, which can be back-calculated using a mixing model equation (Parnell et al. Reference Parnell, Phillips, Bearhop, Semmens, Ward, Moore, Jackson, Grey, Kelly and Inger2013). In many cases, the end-member and mixture compositions exhibit considerable variability (Parnell et al. Reference Parnell, Inger, Bearhop and Jackson2010). By sampling the distributions representing each source and the mixture using a Markov chain Monte Carlo (MCMC) approach and iteratively solving the mixing model equation, Bayesian mixing models generate posterior distributions that account for this variability (Parnell et al. Reference Parnell, Inger, Bearhop and Jackson2010, Reference Parnell, Phillips, Bearhop, Semmens, Ward, Moore, Jackson, Grey, Kelly and Inger2013; Phillips et al. Reference Phillips, Inger, Bearhop, Jackson, Moore, Parnell, Semmens and Ward2014). The posterior mean (or median) provides an estimate of the fractional contribution from a given source to the mixture; the distribution spread relates the estimate uncertainty. A discrete solution can be obtained for all sources, provided the source fractions sum to 1 and the number of sources considered in the analysis exceeds the number of tracers by n + 1 (Hopkins and Ferguson Reference Hopkins and Ferguson2012; Parnell et al. Reference Parnell, Phillips, Bearhop, Semmens, Ward, Moore, Jackson, Grey, Kelly and Inger2013). The mixing model finds a geometric solution to describe the mixture based on its distance from the n end-members in n − 1 dimensions.
Mixed fossil deposits are analogous in many respects to the mixed sedimentary deposits studied by geochemists. The percent abundance of particular fossil taxa in VMBs provide a metric analogous to the isotope ratios used by geochemists to characterize different sediment deposits. A Bayesian mixing model can, potentially, be used to back-calculate source paleocommunity contributions to mixed fossil deposits in terrestrial formations using the relative abundance of taxonomic groups as tracers, provided we can define the source paleocommunity end-members. Terry (Reference Terry2010) used a frequentist mixing model to quantify habitat source mixing in Holocene owl pellet assemblages, using live community data to define the source communities. For the purposes of this study, we define a source paleocommunity as a “paleocommunity type” that can be associated with a single habitat type. Studies in the Upper Cretaceous (Campanian) Belly River Group in Alberta, Canada, have identified particular vertebrate microfossil assemblages that appear to be derive from “paleocommunity types” associated with a predominant habitat type, including several with independent sedimentological evidence for the inferred paleoenvironment (Brinkman Reference Brinkman1990; Eberth and Brinkman Reference Eberth and Brinkman1997; Brinkman et al. Reference Brinkman, Braman, Neuman, Ralrick, Sato, Currie and Koppelhus2005a). In comparison, well-studied channel deposit VMBs in the Dinosaur Park Formation (DPF) appear to exhibit considerably more habitat source mixing (Fig. 1), which may be attributable to landscape change associated with marine transgression. These sites, both the presumably mixed DPF sites and the possible “end-member” sites, provide an ideal test case for quantifying habitat source mixing bias in different VMBs. Some of the “end-member habitat” sites may also be mixed to some degree. Note, for example, chondrichthyans at terrestrial Oldman Formation sites, including BB100, BB103, and BB107, and dinosaur specimens at the marine site BB96 (Brinkman et al. Reference Brinkman, Braman, Neuman, Ralrick, Sato, Currie and Koppelhus2005a; Cullen and Evans Reference Cullen and Evans2016). Beavan and Russell (Reference Beavan and Russell1999) interpreted the terrestrial and freshwater specimens at marine site BB96 as reworked on account of their worn preservation. We account for this mixing and other potential biases among assemblages representing the end-member “paleocommunity types” by generating distributions for each end-member and incorporating these distributions into the model using MCMC sampling. We hypothesize that (1) assemblages representing distinct paleocommunity sources can be differentiated using taxon percent-abundance distributions and (2) specimen contributions from terrestrial habitat sources should generally decrease at stratigraphically younger DPF channel deposit sites on account of marine transgression. Note that the second hypothesis does not necessarily assume a constant, continuous transition over time, because the timing and pattern of marine transgression during DPF deposition was likely complex, involving multiple landward incursions (Brinkman Reference Brinkman1990; Beavan and Russell Reference Beavan and Russell1999).
Belly River Group stratigraphic column, showing placement of vertebrate microfossil sites used in this study (modified from Cullen and Evans Reference Cullen and Evans2016). Mixed Dinosaur Park Formation channel deposit sites shown in italics; Oldman Fm. sites represent the terrestrial end-member, Onefour mud-filled incised valley (OMFIV) sites represent the freshwater end-member, and BB96, L2377, and the Foremost Fm. sites represent the marine end-member. See Eberth and Brinkman (Reference Eberth and Brinkman1997) and Cullen and Evans (Reference Cullen and Evans2016) for maps of site locations. Abbreviations: Fm., Formation; LCZ, Lethbridge Coal Zone; Ma, million years ago; ss, sandstone.
Geologic Setting
The DPF overlies the Oldman Formation, which records a well-drained, terrestrial floodplain paleoenvironment, and underlies the marine Bearpaw Formation (Brinkman et al. Reference Brinkman, Russell, Peng, Currie and Koppelhus2005b; Eberth Reference Eberth, Currie and Koppelhus2005). A sequence of interbedded sandstones, mudstones, and shale laminae in the uppermost DPF Lethbridge Coal Zone (LCZ) suggests episodic Bearpaw Sea incursions (Beavan and Russell Reference Beavan and Russell1999). DPF microfossil deposits at opposite ends of the transgressive sequence appear to record predominantly terrestrial and predominantly marine assemblages, respectively; whereas stratigraphically intermediate assemblages appear to exhibit compositional overlap (Brinkman Reference Brinkman1990). Cullen and Evans (Reference Cullen and Evans2016) specifically characterize the youngest DPF sites, BB119 and BB115, as “transitional” with overlying marine assemblages (Beavan and Russell Reference Beavan and Russell1999; Cullen and Evans Reference Cullen and Evans2016; Cullen et al. Reference Cullen, Fanti, Capobianco, Ryan and Evans2016). However, we cannot unambiguously categorize DPF taxa as marine or terrestrial, because some may have been amphibious, whereas others may have inhabited freshwater habitats, such as those represented by the Onefour mud-filled incised valley (OMFIV) sites (Eberth and Brinkman Reference Eberth and Brinkman1997). This potential complexity warrants the use of a three-source, two-tracer Bayesian mixing model framework.
Methods
Data
As with past studies (Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004; Evans et al. Reference Evans, McGarrity, Ryan, Currie and Koppelhus2014; Cullen and Evans Reference Cullen and Evans2016), we ranked the possibly mixed, DPF microfossil channel deposit sites in or near Dinosaur Provincial Park stratigraphically, from oldest to youngest: BB51, BB86, BB97, BB78, BB117, BB25, BB106, BB120, BB54, BB75, BB94, BB102, BB108, BB119, and BB115 (Fig. 1). We excluded three DPF bonebed sites identified as crevasse-splay deposits, BB98, BB31, and BB104 (Cullen and Evans Reference Cullen and Evans2016) from our analysis to limit possible confounding effects from differences in deposit formation among the mixed assemblages.
Source Paleocommunity End-Members (Prior Information)
We compared the DPF channel deposit sites with the separate fossil assemblages associated with predominantly marine, freshwater, or terrestrial paleoenvironments. The taxonomic compositions of these other, habitat-specific assemblages provided the prior information necessary for parsing specimen contributions from these possible source habitats to the mixed deposits. Brinkman et al. (Reference Brinkman, Braman, Neuman, Ralrick, Sato, Currie and Koppelhus2005a) identified BB96 and L2377 in the DPF LCZ as marine shoreface deposits, on account of preserved dinoflagellate cysts. We also used three marine shoreface deposits in the underlying Foremost Formation, PHR1, PHR2, and PHRN (Cullen and Evans Reference Cullen and Evans2016), to define the marine paleocommunity. The 12 OMFIV sites in the upper DPF LCZ represented the freshwater paleocommunity, an interpretation supported by freshwater microfossils and gastropods (but note speculation that the five eastern OMFIV sites may record a more brackish paleoenvironment: Eberth and Brinkman [1997]). We used seven Oldman Formation sites that occur in Dinosaur Provincial Park, BB105, BB118, BB71, BB121, BB103, BB100, and BB107, to define our terrestrial paleocommunity (Cullen and Evans Reference Cullen and Evans2016). The Oldman Formation has been interpreted as a more inland, seasonally arid, terrestrial paleoenvironment (Brinkman Reference Brinkman1990; Eberth and Hamblin Reference Eberth and Hamblin1993; Cullen and Evans Reference Cullen and Evans2016). All source habitat sites used in the analysis occurred in either the same section as the mixed channel deposits (the Oldman and Foremost sites) or in the DPF (the OMFIV sites) to control for geographic confounding factors. These habitat-specific assemblages provided relatively distinct habitat “end-members” for comparison with the compositionally intermediate DPF channel deposit sites.
Bayesian Mixing Model
We characterized each of the mixed DPF channel site assemblages and the site assemblages representing the different terrestrial, freshwater, and marine habitat source end-members according to the percent abundance of different class/order-level vertebrate taxa—Chondrichthyes, Osteichthyes (here restricted to actinopterygians), Amphibia, Testudinata, Choristodera, Crocodylia, and Dinosauria—the model tracers (Fig. 2). Lepidosauromorpha, Plesiosauria, Aves, and Mammalia occurrences affected the percent abundance of other taxa in assemblages but were not specifically evaluated as tracers, because they occurred infrequently at sites (Cullen and Evans Reference Cullen and Evans2016). We relied on taxa that occurred at both the mixed channel sites and in our end-member assemblages.
Vertebrate percent-abundance ranges at end-member sites by end-member habitat type: marine (n = 5 sites), freshwater (n = 11 sites), and terrestrial (n = 7 sites). Box plots show group median (bar) and mean (point); see Table 1 for taxon percent-abundance mean and SD values by end-member habitat type. Abbreviations: Mar, marine; Fresh, freshwater; Terr, terrestrial.
Mean percent abundance of vertebrate taxa across rarefied sites (site specimen n = 100; subsampled 1000×) representing source habitat end-members.
A Bayesian framework allowed us to account for observed compositional overlap among the end-member groups (e.g., Eberth and Brinkman Reference Eberth and Brinkman1997), which limited our ability to define the end-members based solely on taxon presence/absence (Fig. 2). We used an MCMC approach to sample the percent-abundance ranges for each habitat source (the prior distributions) and each mixed assemblage (the data) and iteratively solve for the fractional contribution from each end-member source to each mixed assemblage using the mixing model. The fractional results from each iteration were aggregated to generate a posterior distribution for each source habitat. The posterior distribution mean provided an estimate for the average fractional contribution from that habitat source to the mixed assemblage. The posterior distribution standard deviation allowed us to quantify our confidence in that estimate.
We used the observed percent abundances of two different vertebrate taxon tracers, designated A and B, to parse the relative contributions from these three habitat types to each mixed assemblage. As mentioned in the “Introduction,” the number of tracers in the analysis must equal the number of sources minus 1 to yield discrete solutions. Only one tracer would be needed to solve for the proportional similarity of an intermediate composition relative to two end-members (i.e., an intermediate location on a line between two points). The biplots in Figure 3 show the locations of the mixed sites in two dimensions (taxonomic compositions) relative to three end-members (habitat sources). Our Bayesian mixing model statement for a given DPF mixed assemblage, sitei, was as follows:
$$\eqalign{\% A_{{\rm sit}{\rm e}_{\rm i}} = &({\rm \phi }_{\rm T}x\% A_{\rm T}) + ({\rm \phi }_{\rm F}x\% A_{\rm F}) \cr & + ({\rm \phi }_{\rm M}x\% A_{\rm M})} $$
$$\eqalign{\% B_{{\rm sit}{\rm e}_{\rm i}} = &({\rm \phi }_{\rm T}x\% B_{\rm T}) + ({\rm \phi }_{\rm F}x\% B_{\rm F}) \cr & + ({\rm \phi }_{\rm M}x\% B_{\rm M})} $$
$$1 = {\rm \phi} _{\rm T} + {\rm \phi} _{\rm F} + {\rm \phi} _{\rm M}$$where T, F, and M denoted the terrestrial, freshwater, and marine end-member habitat types, respectively, and ϕ represented each end-member habitat's fractional contribution to the mixed assemblage. The fractional contributions from the three habitat types had to sum to 1. This framework only allowed us to consider two taxa at a time, so we reran the model on each site using different combinations of the taxon tracers.
Biplots showing the percent abundance of particular taxa at 15 Dinosaur Park Formation (DPF) channel deposit sites (numbers = BB localities) and at terrestrial, freshwater, and marine end-member sites (ranges = group percent-abundance means ± 1 SD). The Bayesian mixing model can solve for the fractional contribution from each source to DPF channel deposit sites that fall within the “mixing polygon” (dashed lines).
We employed the standard Bayesian mixing model SIAR (‘SIAR’ package v. 4.2; Parnell and Jackson Reference Parnell and Jackson2015) to run these analyses in R (v. 3.4.2; R Development Core Team 2017). SIAR uses replicate observations to generate a distribution for the mixed assemblage for comparison with the end-member distributions. We generated 100-specimen replicate samples for each site by randomly subsampling the assemblages without replacement. A rarified sample size of 100 specimens allowed us to control for differences in assemblage sample size and include the least productive channel site, BB71, which only yielded 184 specimens (Cullen and Evans Reference Cullen and Evans2016). We randomly drew 100 specimens from each site 1000 times in R using the rrarefy() function in the ‘vegan’ package (v. 2.4-3) and calculated the percent abundance of each vertebrate taxon in each of the 1000 subsamples generated for each site (see Oksanen [Reference Oksanen2018] for the code to run rrarefy()). We also subsampled the assemblages representing the three source habitat types and calculated taxon percent abundances for each site at the rarefied specimen count n = 100, which were then averaged to produce taxon percent-abundance means and standard deviations for each habitat end-member. We used the siarmcmcdirichletv4() function in the ‘SIAR’ package to run the Bayesian mixing model (see Parnell and Jackson [Reference Parnell and Jackson2015] for the model code). SIAR generated Dirichlet-distributed priors for each habitat from the mean and standard deviation estimates for each taxon used as a tracer in the model. We then ran the model using 500,000 MCMC iterations for each run. Posterior distributions were generated for each source after discarding the first 50,000 iterations (the model “burn-in”).
We considered combinations of taxa that allowed us to differentiate source habitats relative to the mixed assemblages, such that at least 7 of the 15 mixed assemblages yielded intermediate compositions relative to the marine, freshwater, and terrestrial end-members. We determined which combination of taxa generated the best model for each mixed assemblage by calculating the Bayesian information criterion (BIC) for each model run using the following equation (cf. Burnham and Anderson Reference Burnham and Anderson1994; Hondula and Pace Reference Hondula and Pace2014: Supplement):
$${\rm BIC}{\rm} = {\rm} -2{^\ast}{\rm ln}\lsqb {{\lpar {{\rm RSS}/{\rm n}} \rpar }^{(n/-2)}} \rsqb {\rm} + K{^\ast}{\rm ln}\lpar n \rpar $$where n is the number of samples for each site, and K is the number of habitat sources. Each assemblage was represented by 1000 subsamples (n = 1000), and each model run included the same three source habitats (K = 3). Model BIC, therefore, only varied with respect to the residual sum of squares (RSS), the mean residual error for each taxon tracer. The RSS was found by averaging the SD output for both tracers and summing the squared means (cf. Hondula and Pace Reference Hondula and Pace2014: Supplement). In cases in which multiple models yielded results for a given site, we selected the model with the lowest RSS.
We determined whether the stratigraphic position of a DPF mixed assemblage was significantly correlated with the source fraction from each habitat type using Spearman rank correlation (rcor.test in R ‘ltm’ package v. 1.1-0).
Results
All of the vertebrate taxa, with the exception of Plesiosauria, occurred in all three end-member habitat types, preventing habitat determination using simple taxon presence/absence (Table 1, Fig. 2). However, the average percent abundance of many of these taxa differed by habitat type (Table 1, Fig. 2). Chondrichthyan specimens accounted for 53.47 ± 4.5 (SD) % of the marine site assemblages but only 7.22 ± 2.2 (SD) % of the freshwater sites and 3.95 ± 1.6 (SD) % of the terrestrial sites. Despite being generally rare, mammals were more abundant at terrestrial sites (0.61 ± 0.7 [SD] %) than at marine (0.06 ± 0.2 [SD] %) or freshwater sites (0.11 ± 0.2 [SD] %). Dinosaurs were also considerably more abundant at terrestrial sites, 20.67 ± 3.6 (SD) %; whereas osteichthyan specimens occurred with the highest average percent abundance at freshwater sites (86.69 ± 2.8 [SD] %). Exceptions included Testudinata and Choristodera, which occurred with a low and fairly consistent average percent abundance in all three habitat types (Table 1).
Simultaneously differentiating the three habitat types from one another required considering pairs of taxa with unique occurrence patterns (e.g., Chondrichthyes at marine sites and Dinosauria at terrestrial sites). Of the 21 taxon pair combinations, Chondrichthyes–Amphibia (C-A) and Chondrichthyes–Dinosauria (C-D) facilitated the clearest habitat end-member separation (Fig. 3). The percent-abundance range (mean ± 1 SD) for each taxon in each habitat type showed either no or minimal overlap, and 14 of the 15 mixed assemblages plotted between the different end-member ranges (Fig. 3). The percent-abundance ranges for other taxa overlapped one another for at least two of the three habitat types, limiting their utility as tracers. For example, the Osteichthyes percent-abundance distributions overlapped at marine and terrestrial sites, and the Crocodylia distribution overlapped at marine and freshwater sites (Fig. 3). Choristodera and Testudinata percent-abundance ranges showed considerable overlap for all three habitat types (Fig. 3).
The percent-abundance patterns for Chondrichthyes, Amphibia, and Dinosauria allowed us to solve the mixing model and calculate fractional specimen contributions from each habitat source. We only applied the Bayesian mixing model to the 14 DPF mixed sites with taxon percent abundances that plotted within end-member mixing polygons, which encompassed the percent-abundance ranges for the three habitat sources with respect to both taxon tracers (Fig. 3). Seven of the 15 DPF channel deposit sites exhibited C-A percent abundances that were intermediate between the three source habitat ranges (BB117, BB54, BB75, BB94, BB102, BB108, and BB119); whereas 9 of the DPF sites exhibited intermediate C-D percent abundances (BB51, BB86, BB97, BB25, BB120, BB75, BB94, BB102, and BB119). Three DPF sites did not plot within either of these mixing polygons (BB78, BB106, and BB115), although BB78 and BB106 plotted within the mixing polygon for the Amphibia–Dinosauria comparison (A-D), along with BB86. The stratigraphically youngest site, BB115, did not plot within the mixing polygon for any of the 21 combinations of taxa, due in part to its unusually high percent abundance of chondrichthyan specimens (61.9%) (Table 2).
Mean percent abundance of vertebrate taxa at rarefied Dinosaur Park Formation sites (site specimen n = 100; subsampled 1000×); sites listed from oldest to youngest. Chond, Chondrichthyes; Ostei, Osteichthyes; Amph, Amphibia; Lepido, Lepidosauromorpha; Plesio, Plesiosauria; Testud, Testudinata; Choris, Choristodera; Croc, Crocodylia; Dino, Dinosauria; Mamm, Mammalia.
We ran the Bayesian mixing model on all of the DPF sites that fell within the mixing polygons defined by the C-A, C-D, and A-D taxon tracer combinations (n = 14) and calculated an RSS for each model output. In cases in which we were able to use different tracer combinations to parse source contributions at a single site, we selected the mixing model output with the lowest RSS to represent that site (Table 3). The mixing model generated well-constrained habitat source posterior distributions (defined by SD) for each of the DPF sites that we were able to model (Fig. 4, Table 4). At most sites, all three posterior distributions showed clear separation. The posterior distributions for the terrestrial and freshwater habitats overlapped one another at BB108. The marine and terrestrial distributions overlapped slightly at BB94, and the marine and freshwater distributions overlapped at sites BB78, BB86, and BB51 (Fig. 4).
Bayesian mixing model posterior distributions for each Dinosaur Park Formation channel deposit site; the model with the lowest residual sum of squares is shown (Chondrichthyes–Amphibia [C-A], Chondrichthyes–Dinosauria [C-D], or Amphibia–Dinosauria [A-D]); the posterior distribution means provide an estimate of the specimen contribution from each source at each site. Note: site BB115 plotted outside the mixing polygon for all 21 models.
Mixing model residual sum of squares for Dinosaur Park Formation channel deposit sites modeled using multiple combinations of taxa. C-D, Chondrichthyes–Dinosauria; C-A, Chondrichthyes–Amphibia; A-D, Amphibia–Dinosauria.
Source fraction output (posterior means and SD) from the mixing model with the lowest residual sum of squares for each Dinosaur Park Formation channel deposit site; sites listed from oldest to youngest. C-D, Chondrichthyes–Dinosauria; C-A, Chondrichthyes–Amphibia; A-D, Amphibia–Dinosauria.
We quantified the fractional contribution from each habitat type to the 14 modeled DPF assemblages using the posterior distribution means. The associated SD provided an indication of our confidence in each contribution estimate (Table 4). Ranking the sites stratigraphically and comparing their habitat source fractions revealed several source contribution changes over time. Terrestrial source contributions were significantly correlated with site age rank (ρ = 0.701, p = 0.007). The seven oldest DPF sites, from BB51 to BB106 in Table 4, were primarily “terrestrial,” with contributions averaging 0.86 ± 0.05 (SE). In comparison, terrestrial contributions averaged 0.22 ± 0.09 (SE) at the seven youngest sites and 0.07 ± 0.03 (SE) at the three youngest sites. Marine contributions generally decreased with site age, but the correlation was not statistically significant (ρ = −0.429, p = 0.128). Three older sites (BB86, BB78, and BB106) had low marine source fractions (<0.1; Table 4). However, two of the seven youngest sites also had marine fractions <0.1, given freshwater fractions of 0.99 and 0.31 at sites BB54 and BB75, respectively (Table 4). Sites BB120 and BB94 showed the most habitat mixing, with fractional contributions from all three habitat sources that exceeded 0.2. The correlation between freshwater contributions and age rank was not statistically significant at the Bonferroni corrected α = 0.008 (ρ = −0.596, p = 0.028)—the freshwater fraction fluctuated between 0.10 and 0.99 at the seven youngest sites; however, we observed a significant negative correlation between freshwater and terrestrial contributions (ρ = −0.851, p < 0.001).
Discussion
Paleocommunity Changes during DPF Deposition
Unique percent-abundance signatures for Chondrichthyes, Amphibia, and Dinosauria allowed us to successfully differentiate likely terrestrial, freshwater, and marine source paleocommunities (hypothesis 1) and quantify their relative contributions to the compositionally mixed, DPF channel deposit assemblages. These three taxa were each proportionally more abundant in a single habitat type, consistent with expectations about their life habits. We were not able to successfully use occurrence patterns for other taxa, because they occurred at either very low or similar percent abundances in multiple habitats.
The mixing model results confirmed increased paleocommunity mixing in younger vertebrate microfossil sites within this stratigraphic sequence (Fig. 5). The second hypothesis, that terrestrial source contributions would decrease at younger sites, was also supported. However, most of the observed paleocommunity mixing was not attributable to direct marine incursion but rather to a significant, inverse relationship between freshwater and terrestrial specimen contributions. The rise in base level over time apparently increased the prevalence of freshwater habitats on the landscape, beginning at site BB120 and continuing through younger sites until BB102 (Fig. 5, Table 4). Studies have documented similar changes in aquatic community prevalence and complexity in response to other marine transgressive events (e.g., Freitas et al. Reference Freitas, Andrade and Cruces2002; Rowe Reference Rowe2007; Raabe and Stumpf Reference Raabe and Stumpf2016). Most of the terrestrial specimens at the oldest DPF sites (BB51 through BB75 in Fig. 5 and Table 4) likely originated in contemporaneous, terrestrial paleocommunities living at or near the depositional site, possibly reflecting lateral channel migration. The more modest terrestrial fractions in the youngest deposits (BB94 through BB119 in Fig. 5 and Table 4) likely accumulated through lag deposition as fluvial channel density increased and scoured older sediments. This reconstruction broadly supports Matson's (Reference Matson2010) sequence of DPF paleoenvironmental change based on palynological, macrofloral, and sedimentological data. Our interpretation is also consistent with Eberth's (Reference Eberth1990) stratigraphic assessment that channels in the lower DPF aggraded vertically, whereas the upper DPF records more lateral channel aggradation. Eberth (Reference Eberth, Currie and Koppelhus2005) identified a lower “sandy unit” and an upper “muddy unit” in the DPF (see Fig. 1). Our analysis indicates that the VMB assemblages in Eberth's (Reference Eberth, Currie and Koppelhus2005) “muddy unit” are considerably more mixed (Fig. 5).
Habitat source fractions (posterior means) by Dinosaur Park Formation channel deposit site, arranged stratigraphically from oldest to youngest. Vertical error bars relate posterior SD.
According to our results, only the youngest channel deposit sites (BB108 and BB119) recorded a predominantly estuarine/marine paleocommunity. Note that we were unable to parse habitat mixing at BB115 due to its essentially marine composition. This relatively late appearance of a predominantly estuarine/marine habitat signal—and, consequently, the lack of a significant relationship between marine source fraction and DPF site age rank—is not surprising. Brinkman (Reference Brinkman1990) identified transitional assemblages in this sequence as “coastal” communities. These sites would not necessarily receive abundant specimens from the encroaching marine environment until they were inundated. The landscape experienced multiple small marine incursions during DPF deposition (Beavan and Russell Reference Beavan and Russell1999; Matson Reference Matson2010), but we lacked the chronostratigraphic resolution to associate our marine source fractions with specific, local transgressive events. Some palynological and sedimentological data exist for individual DPF sites (e.g., Eberth Reference Eberth1990; Brinkman et al. Reference Brinkman, Braman, Neuman, Ralrick, Sato, Currie and Koppelhus2005a), but additional stratigraphic work is needed to associate the assemblage mixing patterns we identified with specific events at local and subregional scales.
The results of the Bayesian mixing model analysis can, nevertheless, be used to control for paleocommunity mixing and compare the DPF sites several different ways. We can scale the number of specimens at each site by the fractional contribution from each habitat type and then recalculate the percent abundance of taxa with known life habits using the appropriate habitat-specific specimen counts. All dinosaur specimens, for example, likely derived from terrestrial sources. Dinosaur percent abundance ranged between 8.12% and 24.7% relative to the total number of specimens in the DPF channel deposits, irrespective of site stratigraphic position. However, dinosaur relative abundance increased as a percentage of just the terrestrial specimens at younger sites (Fig. 6). Ceratopsian and hadrosaur terrestrial percent abundances also increased relative to the terrestrial specimen fraction. The ceratopsian terrestrial abundance doubled from 4.36% to 8.65%, and the hadrosaur abundance increased from 57.9% to 64.5%. This result supports stable isotope work suggesting that ceratopsians inhabited open, coastal environments, whereas hadrosaurs marginally preferred inland, forested settings (Fricke and Pearson Reference Fricke and Pearson2008). Note, however, that the number of dinosaur specimens at the youngest sites shown in Figure 6, BB94 and BB102, approached and then exceeded the total expected number of terrestrial specimens at these predominantly freshwater sites (Fig. 5). This suggests that an increasing percentage of dinosaur specimens were either reworked or washed into the freshwater deposits (allochthonous deposition).
Dinosaur, ceratopsian, and hadrosaur relative abundances as a percentage of terrestrially derived specimens at sites, ranked from oldest to youngest. The dashed line shows the dinosaur percent abundance relative to the total number of specimens at each site. Note: the entirely freshwater site BB54 and the estuarine sites BB108 and BB119 were excluded. Error bars relate scaled terrestrial posterior SD.
After accounting for differential mixing, direct assemblage comparisons may also prove instructive. For example, both “Holostean A” and caudates (salamanders) occur at the freshwater-dominated site BB54 and the terrestrial-dominated BB106 (Table 4), but the holostean is proportionally more abundant at the former (39.2 vs. 9.72%), whereas caudates represent a larger proportion of the latter (18.4 vs. 13.4%) (Cullen and Evans Reference Cullen and Evans2016). Caudates were an important constituent of terrestrial paleocommunities during DPF deposition. The occurrence of “Holostean A” at the terrestrial site likely reflects the occurrence of aquatic microhabitats within the otherwise terrestrial paleo-catchment surrounding BB106. We can also compare assemblages with different stratigraphic positions but roughly equivalent habitat source mixtures, for example BB51 and BB120 (Table 4). The higher percent abundance of Lepisosteus (13.1 vs. 0.9%) and the relative decline of “Holostean A” (from 19.6% to 8.9%) at the younger site, BB120 (Cullen and Evans Reference Cullen and Evans2016), likely represents a real freshwater paleocommunity change resulting from landscape evolution over time.
A Bayesian Framework for Vertebrate Paleoecology
The Bayesian mixing model analysis used in this study facilitates robust, repeatable comparisons among fossil assemblages, but the specific results depend on how analyses are framed—a caveat affecting both Bayesian and non-Bayesian mixing models (Cooper et al. Reference Cooper, Krueger, Hiscock and Rawlins2014; Phillips et al. Reference Phillips, Inger, Bearhop, Jackson, Moore, Parnell, Semmens and Ward2014). The mixing model framework requires that we hold beta diversity constant, which entails making prior assumptions about likely paleocommunity sources contributing to deposits. Incorporating separate prior information avoids circularity concerns but introduces questions about how to select end-member data for comparison. Incorporating additional data may, or may not, improve prior distributions. For example, we considered incorporating assemblages from the penecontemporaneous Judith River Formation to define the terrestrial end-member but ultimately excluded these sites because they record a terrestrial paleocommunity with proportionally more amphibian and squamate specimens (e.g., UC-8439 and UC-914 in Rogers and Brady [Reference Rogers and Brady2010]). This additional variability would increase posterior distribution uncertainty. The compositional differences between the Judith and Oldman microsites may be due to a marginal difference in paleolatitude, an additional confounding variable. However, we decided to include the possibly estuarine OMFIV sites in the freshwater end-member (Eberth and Brinkman Reference Eberth and Brinkman1997) in the absence of a clear reason to exclude them. These sites yielded more Hybodus specimens; however, their chondrichthyan percent abundances were lower than those of the other freshwater sites, due to more specimens of the holostean Atractosteus (Eberth and Brinkman Reference Eberth and Brinkman1997). Incorporating these “estuarine” OMFIV sites actually increased the apparent distinction between the freshwater and marine source distributions. Chondrichthyan specimens at these “freshwater” sites strongly suggest euryhaline tolerances of the particular chondrichthyan species in question, given the large number of chondrichthyan specimens at these sites. However, we acknowledge that postmortem specimen mixing may have augmented the chondrichthyan specimen occurrences.
Some compositional homogenization among end-member sites does not preclude their use in the analysis, because model solutions are based on relative, rather than absolute, distances between each mixed assemblage and the end-members. The mixing model analysis is not affected by differences in sample size among sites, because it compares relative percentages. However, biases that affect specific end-member or mixed assemblage distributions in inconsistent ways would affect the results. Site specimen productivity may affect the results if taxon percent abundances are inaccurate at less productive microfossil sites due to biases affecting particular taxonomic groups (Carrano et al. Reference Carrano, Oreska and Lockwood2016). We assumed that any biases affecting the percent abundance of a particular taxon at end-member sites affected the percent abundance of that taxon in the mixed assemblages the same way. In such cases, the bias would affect the absolute location of the mixed sites and the mixing polygon in the Figure 3 biplots but not the relative location of a given mixed assemblage or end-member. However, if a mixed assemblage experienced unusual specimen loss due to specimen transport, for example, percent abundances at that site may be biased in ways that confound that mixing model result. Isotopic studies that employ mixing models occasionally encounter a similar problem, stable isotope ratio depletion or enrichment in particular samples due to diagenesis, which can be addressed by incorporating a diagenetic and/or a temporal uptake correction factor into the model (e.g., Larsen et al. Reference Larsen, Bach, Salvatteci, Wang, Andersen, Ventura and McCarthy2015; Matsubayashi et al. Reference Matsubayashi, Morimoto, Tayasu, Mano, Nakajima, Takahashi, Kobayashi and Nakamura2015). Paleoecologists can run the model with a similar “taphonomic correction factor,” provided they can quantify expected taxon percent-abundance effects. Taphonomic information, including size sorting, weathering class, and element absences (Blob and Fiorillo Reference Blob and Fiorillo1996; Moore and Norman Reference Moore and Norman2009; Belanger Reference Belanger2011; Behrensmeyer and Miller Reference Behrensmeyer, Miller and Louys2012), can provide independent evidence to assess whether a correction is needed to account for differential preservation or other biases affecting individual deposits, separate from habitat mixing.
Although we successfully applied a Bayesian mixing model framework to DPF sites, some data sets may not be conducive to mixing model analysis. We considered applying this analysis to Cloverly Formation VMBs to resolve possible paleoenvironmental differences among the four most productive sites; however, all Cloverly sites exhibited unusually high crocodilian and turtle abundances (Carrano et al. Reference Carrano, Oreska and Lockwood2016). We cannot parse their specimen compositions without a more extreme end-member condition for these sites. The Cloverly sites may represent a unique end-member paleocommunity relative to the DPF sites, which we could model as a fourth habitat type possibly contributing specimens to the mixed DPF channel deposits. We would need a third taxon tracer to run this four-source model, which we could compare with the three-source model using the BIC. However, regardless of whether the BIC suggests improved model-fit, incorporating the “Cloverly habitat” end-member may introduce both temporal and geographic confounding variables that would undermine the results. Using end-member assemblages from different time intervals may, nevertheless, be defensible, provided the end-member distribution accounts for sources of variability, including paleocommunity evolution over time. All studies employing mixing models must carefully consider context and justify applying particular models as appropriate to study questions (Phillips et al. Reference Phillips, Inger, Bearhop, Jackson, Moore, Parnell, Semmens and Ward2014).
Bayesian mixing models offer promising new avenues for deep-time terrestrial paleoecology. Although Bayesian methods may not be appropriate or insightful in all situations, they provide a flexible, repeatable framework for microsite comparison with distinct advantages over more traditional approaches, which cannot account for the multiple, often unknowable, sources of variability that potentially bias microfossil assemblages. We offer the following recommendations for paleoecologists interested in adapting this technique:
1. Studies should collect replicate samples from microfossil sites. We subsampled existing, aggregate collections to generate the mixed assemblage distributions used in this analysis. The resulting, well-constrained posterior distributions implied a high degree of certainty about assemblage composition (see Fig. 5 error bars). Analyzing genuine, replicate samples obtained from sites in the field would provide a better measure of the compositional variability in fossil deposits.
2. We employed Dirichlet distributed priors in our model, a standard approach taken by mixing model studies (Parnell et al. Reference Parnell, Phillips, Bearhop, Semmens, Ward, Moore, Jackson, Grey, Kelly and Inger2013). Future studies may want to evaluate alternate model structures (see Parnell and Inger [Reference Parnell and Inger2016] for information on employing the ‘simmr’ package, a SIAR upgrade that allows for additional model flexibility).
3. The selection of prior information should be justified as appropriate to the questions at hand on a case-by-case basis, regardless of whether the data yield a mathematical solution for the Bayesian mixing model.
4. Results from Bayesian analyses can inform subsequent analyses; however, more primary information, including additional microfossil sites and independent paleoenvironmental evidence, will help generate better-informed models.
5. This analysis can be conducted at different taxonomic levels, but we note that analyses based on higher-level taxa provide a species-independent frame of reference for comparing assemblages with species absences, minimal species overlap, or potential taxonomic biases.
Conclusions
Bayesian mixing models can be used to quantify different habitat source contributions to individual vertebrate microfossil deposits. The results for the mixed DPF sites confirm a decrease in terrestrial source contributions over time, likely attributable to increased habitat complexity on the paleolandscape. The seven older sites appear predominantly “terrestrial” and are directly comparable to underlying Oldman Formation assemblages; however, the seven younger sites show considerable habitat source mixing, indicative of increased landscape change over time. Comparisons based on these sites must account for this mixing by either quantifying habitat-specific percent abundances or limiting comparisons to deposits with similar source fractions. An increase in habitat source mixing may be expected at sites in other formations that record relatively abrupt, landscape-scale paleoenvironmental change. Paleoecologists should account for this potential source of bias when trying to assess paleocommunity responses to these episodes, which may drive major macroevolutionary changes. A Bayesian mixing model framework, therefore, offers significant utility for terrestrial paleoecology, because this approach can help paleoecologists use vertebrate microfossil assemblages to address macroevolutionary questions.
Acknowledgments
We wish to thank G. Wilkinson for advice on Bayesian mixing model analysis in R and P. Wagner for relevant advice on a related project. We are also grateful to the editor and the two anonymous reviewers for their feedback, which helped improve the final manuscript. This work was supported by NSF EAR-1052673 to MTC.
