Introduction
The quality of ecological inference drawn from the fossil record heavily depends on the taphonomic quality and comparability of the samples under investigation (Efremov Reference Efremov1940; Johnson Reference Johnson1960; Fagerstrom Reference Fagerstrom1964; Lawrence Reference Lawrence1968; Raup Reference Raup1976; Behrensmeyer et al. Reference Behrensmeyer, Western and Dechant Boaz1979; Badgley Reference Badgley1986; Smith Reference Smith2001; Blob and Badgley Reference Blob and Badgley2007). This basic tenet holds whether the focus is a single well-preserved specimen or locality, or a suite of well-documented fossil occurrences that permit investigation of geographic patterns and temporal trends (Clyde and Gingerich Reference Clyde and Gingerich1998; Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004). However, in the latter case, the overall comparability of the sites is of primary importance, because paleoecological analyses are more robust when localities share similar taphonomic histories and are sampled and ultimately interrogated in similar fashion.
With this presumably in mind, Behrensmeyer (Reference Behrensmeyer1988a, Reference Behrensmeyer1991) proposed the concept of “isotaphonomy” to promote the practice of identifying and comparing like assemblages in studies of paleoecology and evolution. By design, an isotaphonomic study compares fossil localities that share “an equivalent suite of natural sampling conditions” and represent a “specific habitat or suite of habitats or facies” (Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000: p. 130). Sites are typically deemed isotaphonomic if they derive from the same facies and were exposed to the same potential biasing agents (e.g., Clyde and Gingerich, Reference Clyde and Gingerich1998; Chaney and DiMichele Reference Chaney and DiMichele2003; Aguirre et al. Reference Aguirre, Richiano, Farinati and Fucks2011). The implication is that sites that satisfy this definition will yield comparable samples that should be free from concerns of disproportionate bias.
The concept of isotaphonomy in the vertebrate fossil record has been evaluated using vertebrate microfossil bonebeds (VMBs, a.k.a. microvertebrate assemblages, vertebrate microfossil assemblages, microsites), which are concentrations of predominantly small, disarticulated, and often taxonomically diverse vertebrate hard parts (Wood et al. Reference Wood, Thomas and Visser1988; Eberth et al. Reference Eberth, Rogers and Fiorillo2007; Rogers and Kidwell 2007; Rogers et al. 2007; Rogers and Brady Reference Rogers and Brady2010). VMBs are fairly common in Mesozoic and Cenozoic terrestrial records, where they have been exploited to recover otherwise rarely found small-bodied taxa (e.g., Simpson Reference Simpson1926; Sloan and Van Valen Reference Sloan and Van Valen1965; Estes Reference Estes1969; Sloan Reference Sloan1969; Sahni Reference Sahni1972; Archibald Reference Archibald1982; Lillegraven and McKenna Reference Lillegraven and McKenna1986; Lillegraven and Eberle Reference Lillegraven and Eberle1999; Foster and Heckert Reference Foster and Heckert2011) and studied to estimate relative taxonomic abundance and species richness in ancient vertebrate communities (e.g., Estes Reference Estes1964, Reference Estes1976; Estes and Berberian Reference Estes and Berberian1970; Dodson Reference Dodson1987; Bryant Reference Bryant1989; Brinkman Reference Brinkman1990; Peng et al. Reference Peng, Russell and Brinkman2001; Sankey Reference Sankey2001; Jamniczky et al. Reference Jamniczky, Brinkman and Russell2003; Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004, Reference Brinkman, Eberth and Currie2007; Heckert Reference Heckert2004; Carrano and Velez-Juarbe Reference Carrano and Velez-Juarbe2006; DeMar and Breithaupt Reference DeMar and Breithaupt2006; Heckert et al. Reference Heckert, Mitchell, Schneider and Olsen2012). Workers have applied a variety of statistical tests to VMB assemblages to accomplish the latter. In some cases the findings are consistent with the premise that samples derived from similar facies yield similar—in some cases virtually indistinguishable—taphonomic signatures (e.g., Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004; Rogers and Brady Reference Rogers and Brady2010). Yet in other studies, significant distinctions in taphonomic attributes were discerned in sites that would otherwise be classified as isotaphonomic on the basis of facies context and taphonomic mode (Blob and Fiorillo Reference Blob and Fiorillo1996; Wilson Reference Wilson2008; Peterson et al. Reference Peterson, Scherer and Huffman2011).
In this report we explore the concept of isotaphonomy with a focus on VMBs preserved in the Upper Cretaceous (Campanian) Judith River Formation in north-central Montana (Fig. 1). The Judith River Formation preserves an abundance of VMBs in well-documented facies contexts (Rogers Reference Rogers1998; Rogers and Kidwell Reference Rogers and Kidwell2000; Rogers and Brady Reference Rogers and Brady2010), allowing comparison of taphonomic signatures and thereby providing an ideal opportunity to assess isotaphonomy and evaluate pathways to accumulation. The sites have been thoroughly, repeatedly, and rigorously sampled over more than two decades of fieldwork, and preserve a taxonomically diverse ecosystem that includes aquatic, amphibious, and terrestrial forms ranging in body size from grams to tonnes.
Judith River Formation (JRF) study area and regional stratigraphy. A, The six vertebrate microfossil bonebeds (VMBs) under investigation are located in the Upper Missouri River Breaks National Monument (UMRBNM; shaded gray) in north-central Montana. B, Schematic stratigraphic cross section through the study interval, modified from Gill and Cobban (Reference Gill and Cobban1973) and Rogers et al. (Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016). The Judith River Formation correlates to the west with the terrestrial Two Medicine Formation and is bounded above and below by marine shales of the Bearpaw and Claggett formations. Regional regressive (R7, R8, R9) and transgressive (T8, T9) phases of deposition are from Kauffman (Reference Kauffman1977). The mid-Judith discontinuity (MJD; Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016) separates the McClelland Ferry Member below from the Coal Ridge Member above. The six VMBs featured in this report occur in the Coal Ridge Member, in strata dated between 76.2 Ma and 75.2 Ma (Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016). Positions of Santonian/Campanian and Campanian/Maastrichtian boundaries (Ogg and Hinnov Reference Ogg and Hinnov2012) are approximated.
Here we compare taphonomic attributes among six VMB sites both within and between facies, and evaluate potential biases that relate to sampling strategies and depositional context. We also review and refine scenarios relating to the accumulation of VMBs. Our findings indicate that Judith River VMBs are indeed isotaphonomic and are particularly well suited for the study of Late Cretaceous wetland ecosystems. More broadly, the qualities of these localities are common to many Mesozoic VMBs, opening up the possibility that such sites can be rigorously and robustly compared across greater expanses of time and space.
Geologic Background and Site Descriptions
The Judith River Formation is a classic fossil-bearing unit in the Western Interior, with references to its fossil record dating back to the mid-1800s (Leidy Reference Leidy1856, Reference Leidy1860; Meek and Hayden Reference Meek and Hayden1856; Hayden Reference Hayden1857). In its north-central Montana type area (Fig. 1A) the formation consists of ~180 m of terrestrial and shallow-marine strata of middle to late Campanian age (Goodwin and Deino Reference Goodwin and Deino1989; Rogers Reference Rogers1998; Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016). Judith River strata correlate with terrestrial deposits of the middle and upper Two Medicine Formation to the west (Rogers Reference Rogers1998), and are situated between open-marine deposits of the Claggett and Bearpaw formations (Stanton and Hatcher Reference Stanton and Hatcher1905; Gill and Cobban Reference Gill and Cobban1973; Fig. 1B). The Foremost, Oldman, and Dinosaur Park formations of the Belly River Group represent equivalent strata to the north in Canada (Eberth and Hamblin Reference Eberth and Hamblin1993; Jerzykiewicz and Norris Reference Jerzykiewicz and Norris1994; Hamblin and Abrahamson Reference Hamblin and Abrahamson1996).
A recent reappraisal of the stratigraphy of the Judith River Formation (Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016) advanced three new formal members (Fig. 1B) and linked them to distinct accommodation regimes and shoreline stacking patterns. The VMBs described here occur within the new Coal Ridge Member, a coastal-plain alluvial succession that accumulated landward of backstepping shorelines during transgression of the Bearpaw Sea (Rogers Reference Rogers1994, Reference Rogers1998; Rogers and Kidwell Reference Rogers and Kidwell2000). The Coal Ridge Member spans approximately the upper half of the Judith River Formation and consists of fine-grained fluvial sandstone bodies interstratified with drab gray and olive-green overbank deposits. Fossils of vertebrates, invertebrates, and plants are well represented throughout the formation but are most abundant in terrestrial strata of the Coal Ridge Member in the type area (Sahni Reference Sahni1972; Rogers and Kidwell Reference Rogers and Kidwell2000; Rogers and Brady Reference Rogers and Brady2010). New radioisotopic ages indicate that the Coal Ridge Member accumulated between 76.2 Ma and 75.2 Ma (Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016).
VMBs are a prominent taphonomic (preservational) mode in the Coal Ridge Member, which has yielded several tens of localities thus far. Most VMB localities are preserved in tabular to broadly lenticular beds of mudstone, siltstone, or fine-grained sandstone that can be traced up to 250 m, yielding fossil bone in association with carbonaceous plant debris and abundant freshwater mollusks over the entire extent. Geometry, sedimentology, and fossil content are consistent with deposition in low-energy aqueous settings, such as flood basin ponds or lakes and fluvial backwaters. Rogers and Brady (Reference Rogers and Brady2010) categorized this type of bioclastic accumulation as a pond/lake VMB.
The Coal Ridge Member also preserves VMBs in fluvial sandstone bodies. Bones and teeth in channel-hosted VMBs are admixed with intraclast deposits of claystone pebbles, ironstone pebbles, ironstone steinkerns and burrow casts, carbonaceous debris, silicified wood fragments, and invertebrate shell debris. The thickness of bone-bearing intraclast lag deposits varies among sites, ranging from thin, centimeter-scale pavements draping basal and internal scours to fossiliferous beds 10 or more centimeters thick. The sandstone bodies that host VMBs exhibit an array of sedimentary structures at a variety of scales indicative of flow in active channels, and some show evidence of tidal influence (Rogers and Brady Reference Rogers and Brady2010; Rogers et al. Reference Rogers, Kidwell, Deino, Mitchell, Nelson and Thole2016).
In the present study, three VMBs (UC-8303, UC-8315, UC-914) are of the pond/lake variety and three (UC-8302, UC-8439, UC-913) are channel-hosted assemblages. Sedimentological and general taphonomic characteristics of each site are described in Table 1. Fossil material from these sites is curated at the Science Museum of Minnesota and the National Museum of Natural History, Smithsonian Institution (NMNH).
Characteristics of pond/lake and channel-hosted vertebrate microfossil bonebeds (VMBs) included in this study.
Methods
VMBs investigated in this report were initially documented in two previous studies, one detailing associations of vertebrate skeletal concentrations and discontinuity surfaces (Rogers and Kidwell Reference Rogers and Kidwell2000), and the other investigating the origins of microfossil bonebeds (Rogers and Brady Reference Rogers and Brady2010). In these previous studies the stratigraphy of all six sites was ascertained relative to significant bounding surfaces, and the sedimentology of each site was documented with regard to the host facies and in relation to surrounding beds (Table 1). Specific effort focused on discerning the distribution and packing of bioclasts within the host facies and delimiting the geometry and lateral expanse of each bonebed (Rogers and Brady Reference Rogers and Brady2010).
Fossils for comparative study were recovered via surface collection and bulk sampling. Both methods of collection are commonly employed in studies of VMBs. Surface collections were made on multiple occasions between summer 2002 and summer 2015. Collection teams were charged with recovering all fossil bones on the outcrop surface, including unidentifiable bone fragments (Fig. 2A,B). Appropriate bounds of collection were determined in relation to terrain morphology and expected hillslope processes (Frostick and Reed 1983), and only those elements that could be confidently sourced to the productive bed were sampled. Specimens recovered via surface collection were identified as specifically as possible with regard to element type (anatomical identification, position, and side of body) and taxonomy and were binned into one of nine major groups: osteichthyan fishes, chondrichthyan fishes, amphibians, crocodilians, champsosaurs, turtles, dinosaurs (including birds), squamates, and mammals.
Two distinct methods of fossil collecting in the Judith River Formation. A and B, Surface collecting. Field views of vertebrate fossils on outcrop. All fossil bones were collected, including unidentifiable fragments. C, Bulk sampling and sieving. Schematic progression of a block of fossiliferous matrix reduced to vertebrate and invertebrate bioclasts via sieving and sorting, with vertebrate bioclasts set in frame for image capture and analysis.
Because our goal is taphonomic comparison of bioclast populations among sites rather than the comparison of refined faunal census data, we did not apply corrections to account for breakage and potential overrepresentation of body parts such as scales and scutes (e.g., Blob and Fiorillo Reference Blob and Fiorillo1996). Regardless of identifiability, all specimens were scored for completeness, measured using digital calipers (long axis in millimeters), and characterized by shape (equidimensional/compact, plate-like, elongate, or conical). These four shape categories parallel those used in previous studies of shape representation in microfossil bonebeds (Blob and Fiorillo Reference Blob and Fiorillo1996; Wilson Reference Wilson2008; Peterson et al. Reference Peterson, Scherer and Huffman2011). The scoring of completeness was ordinal, with “1” being complete (or very nearly so), “2” being intact enough to approximate original size (always >50%, and typically >75% represented), and “3” being indecipherable with regard to original bioclast size or shape. Unidentifiable bone fragments relegated to completeness category 3 were included in counts of specimens per site but excluded from analyses of size and, in most cases, shape.
Bulk samples of fossiliferous matrix were collected from each site on multiple occasions between summer 2010 and summer 2012. Material was processed using an automated sieve system designed at Macalester College. Eight nested sieve pairs with openings of 500 μm and 2 mm were washed concurrently, with 1–2 kg of material allocated to each pair. Sieve pairs were slowly submerged and drained in baths of standing water every 10s until free of matrix. The steady movement and slow rate of this process minimized breakage during processing. Remnant bioclasts in the 500 μm and 2 mm sieves were removed and sorted into vertebrate, invertebrate, and plant fractions under light microscopes. Each bioclastic fraction was weighed, and the invertebrate and plant fractions were stored for future study. Counts of vertebrate bioclasts were made during this initial sorting process.
Representative subsamples of approximately 1000 bioclasts per site were studied in detail. The sizes of vertebrate bioclasts in these subsamples, including the omnipresent “unidentifiable” fraction, were documented using microscope-based image analysis (ImageJ; Rasband 1997–Reference Rasband2015) to yield precise measurements of long axes down to the 500 μm recovery limit (Fig. 2C). As with the surface-collected samples, the sieved subsamples were also characterized with regard to identifiability and shape, with identifiable material assigned to the same nine taxonomic categories.
Several different statistical tests were used to compare data sets: chi-square test, Fisher’s exact test, the nonparametric Kolmogorov–Smirnov (K-S) two-sample test, and Spearman’s rank-order correlation. We used Fisher’s exact test in conjunction with the chi-square test because some localities had data categories with counts of less than five. These tests were chosen to facilitate comparison with previous studies of VMB taphonomy (Blob and Fiorillo Reference Blob and Fiorillo1996; Wilson Reference Wilson2008; Rogers and Brady Reference Rogers and Brady2010; Peterson et al. Reference Peterson, Scherer and Huffman2011). Because these tests used multiple comparisons, we applied a Bonferroni correction (e.g., Dunn Reference Dunn1961), dividing the original desired p-value (0.05) by the number of comparisons (15) to achieve a corrected p-value of 0.003. We also faced a “shared zero” problem for Spearman’s rank correlation tests, in which two sites that both lacked the same taxonomic group would thereby appear more similar. To accommodate this, we removed any taxonomic groups that were absent from both sites in a particular pairwise comparison. However, because it is possible that such absences might be genuine, especially at high taxonomic levels and with large samples sizes, we opted to include results without this correction as well.
Results
Fossil Abundance and Identifiability
We investigated a total of 5486 surface-collected specimens, with site-specific yields ranging from 317 specimens (UC-8302) to 2319 specimens (UC-8303) (Table 2). The relative abundance of fossils in a given VMB is difficult to gauge from surface collections because recovery depends partly on extrinsic variables, including (1) the susceptibility of the bone-bearing facies to erosion, (2) the person-hours devoted to collection, and, perhaps most importantly, (3) the extent of the bed itself. For example, UC-8302 spans only a few meters laterally, and relatively few bones were recovered by surface collecting. In contrast, site UC-8303 is by far the most expansive VMB under consideration—extending well over 250 m laterally—and has yielded almost three times as many surface-derived fossils as the next most productive site (Table 2).
Collection data for six VMBs from the Judith River Formation. For the sieve collections, N total refers to all bones recovered from processed matrix, and n subsample refers to the number of vertebrate bioclasts analyzed from each site for this report. Pond/lake assemblages in bold.

The bulk-sampling approach yielded a total of 41,620 vertebrate specimens. Counts of bioclasts recovered via bulk processing are presented in Table 2, as is the mass of bioclastic material recovered from each site, reported both for the total bioclastic population (vertebrate+invertebrate+plant) and for the vertebrate fraction alone. VMBs of the Judith River Formation typically yield hundreds of coarse to very coarse sand-sized vertebrate bioclasts per kilogram of matrix, regardless of facies context. Despite their abundance, these bioclasts represent only a very minor component of the host matrix by percent of total weight (Table 2).
Another measurable aspect of both the sieved samples and related surface collections is the identifiability of the recovered bioclasts, which is arguably a reflection of taphonomic history. Localities that experienced more intense or longer-acting processes that promote bone degradation would be expected to yield fewer identifiable skeletal elements on average. Recovery methods should also impact the overall identifiability of a collection. Surface collecting depends on the identification of fossil bone by collectors, the potential for which could be enhanced by the size and preservational quality of a given bioclast. Thus, larger bioclasts with more readily identifiable features are more likely to be discovered and recovered.
Surface collections vary in the identifiability of their yield of bioclasts from 44% (UC-8303) to 67% (UC-913). On average, ~54% of the material recovered from surface collecting can be identified as a particular skeletal element that belongs to a specific taxonomic group. Sieved collections yield a far lower percentage of identifiable specimens, ranging from ~14% of all bioclasts in the subsample of UC-913 to 26% in UC-8315. On average, ~18% of recovered specimens in sieved collections are identifiable (Table 2).
Fossil Size and Shape Distributions
We determined size distributions of surface collections using only those bioclasts with completeness scores of 1 or 2 (Table 3A, Supplementary Table S1), and excluded fragmentary specimens whose original bioclast size was indeterminate (category 3). Size distributions of sieved collections were ascertained for representative subsets of ~1000 elements per site (Table 3B, Supplementary Table S2).
Size data for collections from six VMBs in the Judith River Formation. Pond/lake assemblages in bold.

The mean long axis of surface-collected specimens is 12.5 mm (Table 3A), and distributions are leptokurtic (to varying degrees) and positively skewed (Fig. 3). The mean long axis of sieved specimens is 1.6 mm (Table 3B), and again distributions are leptokurtic and positively skewed (Fig. 4). Interestingly, in all sieved samples, very coarse sand-sized bioclasts outnumber coarse sand-sized bioclasts by at least 2:1. The tailing off of the finest bioclast fraction was unexpected and suggests that there is biased (incomplete) recovery of the smallest fraction, or a lower threshold to the preservation of bone debris in the VMBs under study. Perhaps bone particles of coarse-sand caliber and finer were more susceptible to chemical dissolution in the burial environment, or perhaps the material properties of bone entail a structural constraint that limits the preservation of particulate debris below a certain size. In any case, we are likely sampling near the true lower size limits of preservation.
Size distributions of vertebrate bioclasts recovered via surface collection at six Judith River Formation VMBs. The mean long axis of surface-collected specimens is 12.5 mm, and distributions are leptokurtic and positively skewed. Pond/lake assemblages in bold. Abbreviations: %, percent of total sample; n, number of specimens.
Size distributions of vertebrate bioclasts recovered via bulk sampling and sieving at six Judith River Formation VMBs. The mean long axis of sieve-collected specimens is 1.6 mm, and distributions are leptokurtic and positively skewed. Pond/lake assemblages in bold. Abbreviations: %, percent of total sample; n, number of specimens.
Cumulative frequency (ogive) plots effectively illustrate both the general similarity of size distributions among samples recovered using the same collection methods and the marked distions in size distributions in samples derived from different collection methods (Fig. 5). Ogive plots also afford the opportunity to assess size distributions from a sedimentological perspective using Folk’s (Reference Folk1966) approach to the estimation of grain size and grain sorting. In surface collections, all sites yield moderately sorted to moderately well-sorted samples of bioclasts of predominantly medium pebble size grade (8–16 mm). In sieved collections, sites yield moderately well-sorted samples of bioclasts of predominantly very coarse sand size grade (1–2 mm).
Ogive plots of surface and sieve collections from six Judith River Formation VMBs. Surface collections from all VMBs consist of moderately sorted to moderately well-sorted samples of bioclasts of medium pebble size grade (8–16 mm). Sieved collections from all bonebeds yield moderately well-sorted samples of bioclasts in the very coarse sand size grade (1–2 mm). Site UC-914 is an outlier with regard to size in the sieved samples. Pond/lake assemblages in gray; channel-hosted assemblages in black.
Despite graphical similarity in size distributions of comparably sampled fossils (surface vs. sieve), most site-to-site comparisons yielded statistically significant differences (Table 4). For example, K-S tests of size distributions among surface collections found significant differences in 9 of 15 possible comparisons (p≤0.003). Similarly, 8 of 15 possible K-S test comparisons for sieve-based collections yielded significant differences at the p≤0.003 level (Table 4). Note that most of the dissimilarities in sieves pertain to comparisons with site UC-914.
K-S tests of size distributions, with significance levels and test statistic D reported for both surface (top row in each cell) and sieve (bottom row in each cell) collections. In surface collections, the majority of comparisons (9 of 15) yield significant differences at the p≤0.003 level. Sieved collections show slightly more similarity in size distributions among sites, with 8 of 15 comparisons yielding significant differences. In sieved collections, all sites are indistinguishable from at least two other sites, with the exception of UC-914, for which all comparisons are highly significant. Comparisons made across facies are enclosed in box. Pond/lake assemblages in bold.
Shape data were recovered from category 1 and 2 (and rarely category 3) bioclasts in surface collections and a representative sample of ~1000 bioclasts in sieve collections (Table 5, Fig. 6). Shape distributions are generally comparable in surface collections, with plate-shaped and equidimensional bioclasts dominating at all sites. Plate-shaped bioclasts are most abundant in four of the six sites, ranging from 33% in UC-8302 and UC-8303 to 82% in UC-913. Equidimensional bioclast frequencies range from 12% in UC-913 to 49% in UC-8302. Elongate and conical elements together represent less than 20% of surface collections from all sites, with conical elements always least abundant (Fig. 6). Shape distributions in sieve-based subsamples show even greater consistency across sites, with equidimensional bioclasts dominating in all localities (Table 5). The only notable outlier in the sieved collections is site UC-914, which has a disproportionate abundance of elongate elements (27%) relative to the other five sites (15% on average).
Shape data were recovered from surface collections and sieved collections from six Judith River Formation VMBs. A, Shape distributions of surface collections. These are generally comparable, with plate-shaped and equidimensional bioclasts dominating in all sites. B, Shape distributions of sieve-based subsamples. These show great consistency across sites, with equidimensional bioclasts dominating in all localities. The only notable outlier is site UC-914, which has a disproportionate abundance of elongate elements relative to the other five sites.
Shape data relating to collections from six VMBs in the Judith River Formation. Pond/lake assemblages in bold.

Chi-square tests and Fisher’s exact tests of shape distributions among surface collections found highly significant differences in most site-to-site comparisons (Table 6). Site UC-913, with its preponderance of plate-shaped bioclasts, was most distinct. Most sieve sites were statistically indistinguishable, or very nearly so, with one exception: UC-914 showed significant differences in all comparisons at the p≤0.003 level (Table 6).
Chi-square tests of shape distributions (df=3), with the chi-square statistic and significance levels reported for both surface (top row in each cell) and sieve (bottom two rows in each cell) collections. Because sample size is small in some sieve-based comparisons, results for the Fisher’s exact test are also reported (in parentheses). In surface collections, UC-8315 and UC-8439 are most similar to other sites (two of five possible comparisons show no significant distinctions), whereas UC-913 is the most distinct: all comparisons are significantly different at the p≤0.003 level. Sieved collections show more similarity overall, with UC-914 as the single notable exception. Comparisons made across facies are enclosed in box. Pond/lake assemblages in bold.
Faunal Composition
We determined faunal composition on a site-by-site basis by placing all identifiable bioclasts recovered in surface and sieved collections (Table 2) into the nine major taxonomic bins (see “Methods”). Surface collections exhibit more richness overall and, with but one notable exception (UC-913), are characterized by appreciable numbers of both terrestrial (dinosaur, squamate, mammal) and amphibious/aquatic (crocodilian, champsosaur, turtle, amphibian, osteichthyan, chondrichthyan) groups (Fig. 7). Particularly abundant groups in surface collections are all characterized by moderate to large body size, and include dinosaurs, crocodilians, turtles, and osteichthyans (specifically gar; Fig. 7, Table 7). In contrast, sieve collections are resoundingly dominated by the skeletal remains of small osteichthyans (minnow scale), with amphibious/aquatic animals representing 95% or more of most collections (Fig. 8).
Faunal composition of surface collections from the six Judith River VMBs. With one exception (UC-913), surface collections are characterized by an abundance of both terrestrial and amphibious/aquatic forms. Particularly abundant groups in surface collections are all characterized by moderate to large body size and include dinosaurs, crocodilians, turtles, and osteichthyans (% gar within Osteichthyes indicated for each site). Pond/lake assemblages in bold.
Faunal composition of sieved subsamples from the six Judith River VMBs. Sieved collections are resoundingly dominated by the skeletal remains of small fishes (% gar within Osteichthyes indicated), with amphibious/aquatic animals representing 95% or more of most collections. Animals with terrestrial affinities are notably less abundant in sieved subsamples. Pond/lake assemblages in bold.
Faunal data recovered from both surface and sieve collections from six vertebrate microfossil bonebeds in the Judith River Formation. Pond/lake assemblages in bold.

Spearman’s rank correlation was employed to test the strength of correlations in relation to taxonomic abundances among faunal assemblages. All possible pairwise comparisons of surface collections yield very strong positive correlation, with r-values ranging from 0.82 to 0.99. Pairwise comparisons of sieve collections are more variable in outcome, with r-values ranging from 0.52 to 0.96 (Table 8).
Spearman rank correlation of faunal assemblages from surface collections (top two rows in each cell) and sieve collections (bottom two rows in each cell). Results with instances of shared zeros culled from comparisons are presented in parentheses. Surface collections yield r-values consistent with very strong positive correlation among data sets. Sieved collections yield r-values consistent with moderate to very strong positive correlation. The significance of correlations varies, and many fail to be significant at the p≤0.003 level, which presumably reflects small sample size. Comparisons made across facies are enclosed in box. Pond/lake assemblages in bold.
Comparative Taphonomy
Studies that address the taphonomy of VMB localities often do so as part of larger considerations of the vertebrate record in specific stratigraphic units (e.g., Dodson Reference Dodson1971; Wood et al. Reference Wood, Thomas and Visser1988; Fiorillo Reference Fiorillo1989; Hunt Reference Hunt1991; DeMar and Breithaupt Reference DeMar and Breithaupt2006; Carrano et al. Reference Carrano, Oreska and Lockwood2016). In these studies, general taphonomic observations are gathered to provide context for associated faunal reconstructions, and related data sets are typically limited in scope. Reports that delve more deeply into the taphonomy of VMBs are relatively few (e.g., Mellett Reference Mellett1974; Maas Reference Maas1985; Eberth Reference Eberth1990; Blob and Fiorillo Reference Blob and Fiorillo1996; Khajuria and Prasad Reference Khajuria and Prasad1998; Brinkman et al. Reference Brinkman, Russell, Eberth and Peng2004; Wilson Reference Wilson2008; Rogers and Brady Reference Rogers and Brady2010; Peterson et al. Reference Peterson, Scherer and Huffman2011) and have mostly focused on questions of how these types of sites originate and how they might express potential bias.
Here we review some of the taphonomic features that serve to characterize and in some cases distinguish our surface and sieve-based collections. Significant biases in size, shape, and faunal representation are evident, consistent with previous studies that have addressed the impact of sampling strategy on fossil collections (Wolff Reference Wolff1975; Kowalewski and Hoffmeister Reference Kowalewski and Hoffmeister2003; Behrensmeyer and Barry Reference Behrensmeyer and Barry2005; Hunter and Donovan Reference Hunter and Donovan2005; Blob and Badgley Reference Blob and Badgley2007; Sessa et al. Reference Sessa, Patzkowsky and Bralower2009; Peterson et al. Reference Peterson, Scherer and Huffman2011; Lyman Reference Lyman2012; Ottens et al. Reference Ottens, Dietl, Kelley and Stanford2012). We follow with a consideration of the taphonomy of the six sites in relation to their depositional setting (pond/lake vs. channel) by comparing taphonomic aspects of size and shape among sites both within and between facies and evaluating potential biases related to depositional context. We conclude with a brief consideration of site UC-914, which presents an interesting outlier among our six VMB localities.
Collection-related Biases and Distinctions
Surface collections from all localities are dominated by vertebrate bioclasts that fall in the medium pebble size grade (8–16 mm long axis). Bulk-sampled and sieved collections from the same sites are dominated by vertebrate bioclasts that approximate coarse to very coarse sand (0.5–2 mm long axis). Vertebrate bioclasts that fall between the medium pebble and very coarse sand size classes are underrepresented in pooled collections (Figs. 3, 4, 5). We interpret the dearth of such bioclasts in the 2–8 mm (granule to “small pebble”) size range as a sampling artifact, because there is no a priori reason to suspect a bimodal size distribution of vertebrate bioclasts in any of the localities under investigation, and, as clarified below, there is good reason to infer some degree of collection bias.
In surface collecting, the likelihood of discovery and recovery of vertebrate fossils should decline with diminishing size (e.g., Payne Reference Payne1972). On an outcrop, small bioclasts (millimeter scale) are less likely to be collected because they are (1) more difficult to find (more prone to fall into surface cracks and crevices and/or to be obscured by other surface detritus, etc.) and (2) more difficult to identify (with regard to element or taxon). This collection bias is compounded by surface processes, including overland runoff after rainstorms, that would preferentially entrain and transport smaller fossils. Data indicate that erosional and transportational processes in shallow overland flows are comparable to those documented for streams (Moss and Walker Reference Moss and Walker1978), especially on relatively shallow slopes of 10° or less (Frostick and Reid Reference Frostick and Reid1983). In badland settings prone to seasonal precipitation, a winnowed surface lag of predominantly coarse bioclastic debris would be expected to mantle an outcrop in close proximity to its source bed(s), and it is this presorted residuum that would remain to be intersected and sampled by collectors. The relatively low recovery rate of the 2–8 mm size fraction—and everything smaller—on the surface can be explained by this combination of collecting bias and selective winnowing.
Bulk sampling, followed by sieving and picking under the light microscope, is a recognized means of managing and reducing bias in fossil collections, especially in relation to the recovery of fossils from small animals (McKenna Reference McKenna1962; Wolff Reference Wolff1973, Reference Wolff1975; Smith et al. Reference Smith, Stearley and Badgley1988; Peterson et al. Reference Peterson, Scherer and Huffman2011). Exhuming and processing fossil-rich matrix from below the weathered carapace should minimize (if not eliminate) the impact of surface processes that might act to preferentially concentrate certain sizes and shapes of fossils. Careful characterization of all bioclasts recovered in sieves should effectively reduce the impacts of biased collecting and provide a more accurate estimation of the actual size distribution.
In sieved subsamples, vertebrate bioclasts in the coarse sand (0.5–1 mm) to very coarse sand (1–2 mm) size range vastly outnumber larger fossils. On average, sand-sized bioclasts represent ~82% of sieve collections, ranging from 67% in UC-914 to 89% in UC-913. Vertebrate fossils with a ≥1 cm long axis are exceedingly rare in our sieve collections, and data suggest that for every few pebble-sized bioclasts recovered on weathered outcrop surfaces, many hundreds of sand-size bioclasts likely remain uncollected and undocumented. Continued bulk sampling and sieving will slowly fill the aforementioned 2–8 mm vacuity in our pooled size distributions (Fig. 9), but as it does, counts of bioclasts in the sand-sized fraction will skyrocket. The apparent bimodality of the pooled size distribution therefore reflects the disjunction between the sieved size distribution and the artificially truncated surface size distribution.
Bimodal size distribution of vertebrate fossils from site UC-8303 recovered via surface collection and sieving. The relative scarcity of vertebrate bioclasts in the 2–8 mm granule to “small pebble” size range is interpreted as an artifact of sampling.
The shape frequencies in our samples of vertebrate fossils also show distinctions that relate to collection history. Equidimensional bioclasts are by far the most abundant shape class in all sieve-based samples, with relative frequencies ranging from 53% to 59%. Plate-shaped bioclasts in turn are notably abundant in surface collections (ranging from 33–82%) and are the most abundant shape category in four of the six sites. In their study of fossil transport and sorting on badland slopes, Frostick and Reid (Reference Frostick and Reid1983) found that plate-shaped elements showed the lowest relative rates of downslope displacement when compared with equidimensional and elongate elements. Surface lags of Plio-Pleistocene fossils in their East African study area were relatively enriched in plate-shaped bioclasts, and the same holds true in our investigation of Late Cretaceous fossils in the badlands of Montana. We interpret the abundance of plate-shaped bioclasts in our surface collections as a potential bias driven by surface slope processes (Frostick and Reid Reference Frostick and Reid1983). The preponderance of equidimensional bioclasts in our sieve collections is a more accurate reflection of original shape representation in our VMBs.
Like size and shape parameters, faunal representation also varies based on collection type and appears to track some of the biases outlined above. Surface collections, which are skewed toward the larger bioclast fraction, are dominated by the fossil remains of animals of medium to large body size, including dinosaurs, crocodilians, champsosaurs, turtles, and gar. Hard parts from these groups represent more than 75% of all identifiable material documented in surface collections. In contrast, the identifiable fractions of bulk-sampled and sieved collections are dominated by fossils derived from small animals. The minute teeth, scales, and vertebrae of minnow-sized fishes are in every case most abundant, followed to a lesser degree by the hard parts of small amphibians (frogs, salamanders, and albanerpetontids) and lizards. Large terrestrial animals, including dinosaurs, are much more rarely documented in our sieved subsamples.
The shape bias described above is also potentially evident in faunal data sets from surface collections, with groups anatomically destined to yield abundant plate-shaped bioclasts, such as gar (ganoid scales) and turtles (dermal scutes), particularly abundant. Others have explored the potential for shape-related bias in faunal data sets derived from VMBs (e.g., Blob and Fiorillo Reference Blob and Fiorillo1996; Wilson Reference Wilson2008; Peterson et al. Reference Peterson, Scherer and Huffman2011), and there is little doubt that the hydrodynamic aspects of shape and size can influence the taxonomic composition of collections both preburial and postexhumation via this anatomical variability.
Finally, it is important to acknowledge that the sieved subsamples described in this study, while less susceptible to size and shape biases driven by surface processes, are still biased from a faunal perspective, because they are far too small to provide an adequate perspective on total faunal composition. This bias can be mitigated with continued processing of matrix coupled with an organized protocol of systematic sieving and rarefaction analysis to confirm adequate sampling (e.g., Jamniczky et al. Reference Jamniczky, Brinkman and Russell2003). Nevertheless, these sieved subsamples do provide a robust indication of the most abundant animals represented in Judith River VMBs and can still provide data sets for comparisons among sites.
Paleoenvironment-related Biases and Distinctions
Given the distinct depositional settings from which our samples were drawn, we would predict taphonomic differences relating to bioclast size and shape. For example, the higher-energy conditions that would have prevailed during accumulation of channel-hosted assemblages (UC-8302, UC-8439, UC-913) should be expressed in size and shape sorting. Unidirectional currents coursing through active channels would have generated shear stresses capable of selectively entraining and transporting particulate matter, including bioclastic debris. In contrast, the lower-energy conditions that prevailed during deposition of sediment and biological detritus in the three quiet-water sites (UC-8303, UC-8315, UC-914) would have been less likely to entrain and sort bioclasts according to their hydrodynamic potentials.
Comparisons of size and shape distributions in Judith River VMBs in relation to facies context yield a variety of outcomes. When viewed in a traditional sedimentological sense using ogive plots (Fig. 5), size frequency distributions of bioclasts from different facies show considerable overlap. Despite compelling graphical similarity in size distributions among sites, most show significant differences when evaluated from a statistical perspective, although there still remains no clear pattern in relation to depositional setting (Table 4). Size frequencies of bioclasts recovered from channel-hosted VMBs are no more similar to one another than they are to those from VMBs representing pond/lake settings.
The same essentially holds true when shape frequencies are compared. In general, surface collections from both channel-hosted and pond/lake localities are dominated by plate-shaped bioclasts, and bulk-sampled and sieved collections are consistently dominated by equidimensional bioclasts. Statistical comparisons of surface collections yield significant differences in most comparisons, but there is no predictable link to facies. Statistical comparisons of sieved collections yield few significant differences and no consistent correlation with facies (Fig. 6, Table 6).
We would also expect faunal representation in our VMBs to vary somewhat given that different depositional settings were sampled, but again there is no indication that assemblages representing channel settings differ substantially from those representing lacustrine settings. Correlations between the rank order of taxa from sites from the same depositional setting were no stronger than correlations from different depositional settings (Table 8). Facies context is not a predictor of faunal composition in this suite of VMBs.
A Preservational Outlier
Site UC-914 is the only notable taphonomic outlier among our VMB localities. It is analogous to the other five sites from a faunal perspective, because dinosaurs, gar, turtles, champsosaurs, amphibians, and small osteichthyan fossils dominate the assemblage. Its surface collection mirrors that of the other sites as well, with an abundance of plate-shaped and equidimensional bioclasts of predominantly medium pebble size—but there the similarities end.
When sieved subsamples are compared, site UC-914 stands apart in both size and shape distributions. With regard to size, UC-914 is significantly different from all other sites (p≤0.003, K-S tests) and is characterized by a consistently coarser but still predominantly sand-sized assortment of bioclasts (Fig. 5, Table 4). UC-914 yields statistically significant distinctions in all site-to-site shape-frequency comparisons (p≤0.003, chi-square and Fisher’s exact tests), preserving an anomalous abundance of elongate bioclasts and far fewer plate-shaped elements (Fig. 6, Table 6).
Site UC-914 is a pond/lake VMB that preserves vertebrate fossils in association with carbonaceous plant debris and freshwater invertebrates (Table 1). However, unlike the other sites, it preserves only molds of mollusk shells and not the original calcite and aragonite material. Diagenetic fluids presumably leached carbonate skeletons from the burial environment after the sedimentary matrix lithified enough to maintain molds. We have not resolved the chemistry that drove this preservational bias, but previous work on the geochemical taphonomy of Judith River VMBs provides some constraints (Rogers et al. Reference Rogers, Fricke, Addona, Canavan, Dwyer, Harwood, Koenig, Murray, Thole and Williams2010). Normalized patterns of rare earth element concentrations collected from fossil bones in UC-914 correspond with those derived from bones in UC-8302, UC-8303, and UC-8439, suggesting that the vertebrate hard parts interacted with similar diagenetic fluids during recrystallization, at least from a trace element perspective. The cerium content of recrystallized bioapatite in UC-914 is also comparable to that found in the other sites (especially UC-8303), indicating similar reduction–oxidation potentials in these burial environments.
Diagenetic fluids capable of dissolving mollusk shell debris may have also impacted the finest vertebrate bioclasts, which, based on basic surface area to volume considerations, would have likely been more susceptible to chemical attack and dissolution than coarser bioclasts. The shape distinctions that render UC-914 unique are more difficult to explain and could be a specific paleoecological signal of this site (see below). That said, the thin teleost scales that bolster counts of plate-shaped elements in all sieved subsamples may have also been susceptible to dissolution, and they do constitute a smaller fraction of UC-914’s collections than at other sites (5.8% vs. 9.6% on average).
Discussion
Isotaphonomy in Concept and Practice
Behrensmeyer (Reference Behrensmeyer1988a, Reference Behrensmeyer1991) introduced the “isotaphonomic approach” in the context of comparing fossil localities with comparable taphonomic attributes and biases through time. The term resurfaced in Behrensmeyer and Hook (Reference Behrensmeyer and Hook1992: Fig. 2.2 and p. 91), who used it to refer to deliberate comparisons of vertebrate localities that derive from similar sedimentary (preservational) contexts, with the goal of managing “taphonomic noise” in paleoecological reconstructions. Behrensmeyer and others (2000) proposed criteria to help identify isotaphonomic assemblages and expanded the scope of isotaphonomic inquiry to include both temporal (as originally proposed) and geographic comparisons.
Most workers have framed their isotaphonomic studies using this procedural approach and have intentionally compared fossil assemblages derived from a specific set of depositional environments. For example, Clyde and Gingerich (Reference Clyde and Gingerich1998) used taphonomic and sedimentological criteria to assemble an isotaphonomic faunal database of Wasatchian mammal localities in the Bighorn Basin of Wyoming. They argued that potential preservational biases were minimized by their focus on surface collections derived exclusively from overbank paleosol facies. Chaney and DiMichele (Reference Chaney and DiMichele2003) took a similar approach in their study of early Permian floras from the Clear Fork Group in north-central Texas by restricting their paleobotanical sampling to mudstone facies that accumulated in quiet-water settings. The implication in these studies, and others of similar ilk (e.g., Gastaldo et al. Reference Gastaldo, Adendorff, Bamford, LaBandeira, Neveling and Sims2005; Aguirre et al. Reference Aguirre, Richiano, Farinati and Fucks2011), is that sites sourced from the same facies, if sampled in similar manner, will yield assemblages that should be free from disproportionate taphonomic bias.
Others have applied a more “strict” definition to the notion of isotaphonomy (e.g., Blob and Fiorillo Reference Blob and Fiorillo1996: p. 433; Moore and Norman Reference Moore and Norman2009), using statistical tests to evaluate taphonomic equivalence among fossil assemblages. For example, Moore and Norman (Reference Moore and Norman2009) compared an assortment of mammal fossils that were surface collected across outcrop belts of the Oligocene Brule Formation in South Dakota. These authors compared the taphonomic attributes of two assemblages using multivariate statistics, arguing that the assemblages were non-isotaphonomic because these tests yielded significant differences. Moore (Reference Moore2012) followed with a similar treatment of data from 36 vertebrate assemblages spanning the terrestrial Mesozoic and Cenozoic. He compiled taphonomic characters from literature review and used multivariate statistical tests to determine whether assemblages were isotaphonomic. Almost all sites under consideration failed Moore’s (Reference Moore2012) “strict” tests of isotaphonomy.
A similar “strict” statistical approach to evaluating isotaphonomy has been applied in several studies of VMB localities. Blob and Fiorillo (Reference Blob and Fiorillo1996) compared two VMB sites preserved in channel sandstones in the Judith River Formation using combined sieved and quarried samples. They approximated fossil size using 1 cm bins and categorized fossil shape using the same sedimentological shape classes that we employ here. Significant differences were reported in both size and shape comparisons. Wilson (Reference Wilson2008) compared size and shape distributions in two VMBs preserved in “floodplain” and “channel-lag” deposits from the Hell Creek Formation. After assembling collections via hand quarrying, measuring specimen long axes with calipers, and assigning shapes using the same suite of descriptors, she reported significant differences in both size and shape distributions. Peterson et al. (Reference Peterson, Scherer and Huffman2011) also compared two VMBs from the Hell Creek Formation representing two different facies (“sandstone” vs. “mudstone”) using material derived from both surface collecting and sieving. These authors did not combine collections in their analysis (contra Blob and Fiorillo Reference Blob and Fiorillo1996) but rather compared them to ascertain collection biases. They grouped fossils into two general size categories, small (1–10 mm) and large (10–20 mm), and described them using the same four shape categories. Significant differences in size and shape were again reported.
The outcomes of these “strict” analyses of isotaphonomy are not particularly encouraging, especially in relation to VMBs. The statistically significant differences reported in most site-to-site comparisons suggest that related faunal data sets may not be reliable records of original vertebrate communities. This in turn implies that the potential for paleoecological reconstruction using VMBs in a broad comparative sense is limited.
However, we argue that these results are based on an unrealistic presumption—that there should be almost no variation in taphonomic attributes across isotaphonomic sites. But every fossil site ultimately represents an individual event of bioclastic accumulation, one that should be expected to exhibit its own idiosyncrasies—and with sufficiently fine methods of sampling and detection, we should expect to record these variations. This should be equally true for VMBs, regardless of whether comparisons are made among sites that accumulated in the same general depositional setting or even among samples collected within a single locality. Indeed, as noted above, the goal of an isotaphonomic approach is precisely to address this naturally occurring “noise” within sedimentary contexts that are known from other data to be similar.
Consequently we find that the expectation of some degree of taphonomic similitude among fossil assemblages collected from similar sedimentary facies is reasonable and follows logically from empirical studies (e.g., Behrensmeyer Reference Behrensmeyer1988b; Wood et al. Reference Wood, Thomas and Visser1988; Eberth Reference Eberth1990), much as one might expect general similarities among members of a species while still being able to observe the distinctions of individual organisms. The critical question, then, is what level of similarity should actually be required to support (or reject) an interpretation of isotaphonomy (e.g., Blob and Fiorillo Reference Blob and Fiorillo1996)? In site-to-site comparisons, how close is “close enough” to reject the likelihood of disproportionate bias? Should some taphonomic characteristics be more heavily weighted than others? Should analyses be focused on subsets of smaller bioclasts that may be more sensitive to chemical and physical taphonomic processes? Will some depositional environments inherently express more taphonomic variability than others? These are but a few of the questions that remain to be resolved.
In its original connotation, isotaphonomy was advanced as a guiding principle that, when applied in a study of many fossil localities, should promote the comparison of samples with roughly equivalent biases. Applying sensitive statistics to large data sets to test for isotaphonomy risks venturing beyond this original definition. Given the general lack of clarity with regard to the definition of “strict” isotaphonomic comparability, and the very reasonable expectation of some degree of taphonomic variability among fossil assemblages, we are not surprised that sensitive statistical analyses detect variations among VMB sites. We agree that fossil assemblages should be evaluated for their taphonomic similitude before bold paleoecologic inferences are advanced, but the “strict” approach, as currently applied, risks focusing too closely on the “noise” instead of the signal. A positive test statistic is certainly reassuring and, as evidenced below, is not unattainable—but its absence is not ipso facto a reason to abandon prospects of meaningful comparison.
Are Judith River VMBs Isotaphonomic?
The Judith River VMB collections detailed in this study are biased. The surface collections have almost certainly been winnowed and sorted by modern slope processes (e.g., Frostick and Reid Reference Frostick and Reid1983) and further modified by collector-imposed bias. The sieved subsamples, while arguably largely immune to these particular pitfalls, are undersampled with regard to total faunal composition (e.g., Jamniczky et al. Reference Jamniczky, Brinkman and Russell2003). Moreover, locality UC-914 has apparently experienced a unique diagenetic history that has altered its fossil content. The good news is that these biases can be managed, especially by additional bulk sampling and sieving, but the question still remains: If collection-related biases are recognized, and adjustments made to compensate, and if diagenetic outliers are identified and excluded, can Judith River VMB localities be compared without major concerns of disproportionate preservational bias? Are they functionally isotaphonomic?
Interestingly, the significant distinctions in size and shape distributions reported in most previous studies (Blob and Fiorillo Reference Blob and Fiorillo1996; Wilson Reference Wilson2008; Peterson et al. Reference Peterson, Scherer and Huffman2011) are not developed in our multisite VMB data set. In graphical comparisons of bioclast size using histograms and ogive plots, our VMB localities show few defining distinctions. Indeed, the site-to-site similarities are compelling when comparisons are made among the sieved subsamples (Figs. 4, 5). In these subsamples, the size distributions of bioclasts from all sites except UC-914 clearly overlap in ogive plots (Fig. 5). From a traditional sedimentological standpoint, the bioclast populations sieved from five of our six sites are essentially the same. The one marginal outlier (UC-914), as detailed above, shows independent evidence of a distinct diagenetic history.
The same basic pattern holds true when shape data are compared among sites. Surface collections, which we acknowledge are biased to some degree, are still generally comparable. Roughly equivalent proportions of shape classes are evident in histograms of the sieved subsamples (Fig. 6). UC-914 is, again, the one site with a unique distribution of shape categories.
The taphonomic comparability of most VMBs is also supported when statistical tests are applied to sieved subsamples (exclusive of site UC-914, for reasons outlined above). When size distributions from the five remaining VMBs are compared (Table 4), 7 of 10 potential site-to-site comparisons yield indistinguishable results at the p≤0.003 level. The case is made even more compelling when tests are applied to shape data from sieved collections. With UC-914 excluded from consideration, all comparisons yield statistically indistinguishable or marginal results at the p≤0.003 level (Table 6). The few site-to-site comparisons that yield significant or marginal test statistics in relation to size distributions have indistinguishable shape profiles.
The case for isotaphonomy among the majority of Judith River VMBs is strong, whether based on graphical comparisons or evaluated in a “strict” statistical manner. This result bodes well for future studies of paleoecology—particularly when it comes to investigations of faunal membership and community structure in Late Cretaceous wetland ecosystems, but perhaps much more broadly. We highlight one potentially important aspect of the Judith River Formation VMB record: localities from different depositional environments (pond/lake vs. channel) are taphonomically indistinguishable in most comparisons of size and shape. This result is counterintuitive, given that hydraulic sorting processes should vary between quiet-water and active-channel settings, but it potentially hints at an even more powerful insight. The taphonomic comparability of VMB localities from different facies may be the signature of broadly similar taphonomic processes across entire alluvial depositional systems. The isotaphonomic equivalence of pond/lake and channel-hosted VMBs could also be the indication of shared history on the path to accumulation and final burial (see below). In either case, there are grounds for arguing the overall comparability of most of the sites studied here.
Refining the VMB Formative Model
Rogers and Brady (Reference Rogers and Brady2010) reviewed existing hypotheses relating to the accumulation of vertebrate skeletal debris in VMBs (e.g., Mellett Reference Mellett1974; Korth Reference Korth1979; Maas Reference Maas1985; Wood et al. Reference Wood, Thomas and Visser1988; Bryant Reference Bryant1989; Hunt Reference Hunt1991; Srivastava and Kumar Reference Srivastava and Kumar1996; Badgley et al. Reference Badgley, Downs and Flynn1998; DeMar and Breithaupt Reference DeMar and Breithaupt2006) and argued for a model that better accommodated the VMB record as represented in the Judith River Formation. In their model, vertebrate hard parts, along with the remains of freshwater invertebrates and plants, initially accumulated in situ as biological detritus in the sedimentary fills of long-lived quiet-water aquatic settings. Over time, durable bioclasts accumulated to concentrated levels to form pond/lake VMBs. With these in place, Rogers and Brady (Reference Rogers and Brady2010) reasoned that channel-hosted VMBs would be an expected outcome of mobile rivers reworking subjacent and laterally disposed bioclast-rich facies. They interpreted channel-hosted VMBs as residual lags that had been reworked from preexisting concentrations but only minimally transported once exhumed from source beds.
The present follow-up study compares some of the same sites using more refined size, shape, and faunal composition data. With this new data set, we briefly revisit and update some of the conclusions and interpretations advanced in Rogers and Brady (Reference Rogers and Brady2010). First, our new data are entirely consistent with previous interpretations of size and shape distributions of vertebrate fossils in VMBs. VMB localities in the Judith River Formation preserve comparable, and in some cases indistinguishable, size and shape distributions of vertebrate bioclasts regardless of facies context (Figs. 3, 4, 5, 6, Tables 4, 6), with the obvious caveat that comparisons are made between samples recovered via the same methods. The final, consolidating point remains the same: pond/lake and channel-hosted VMBs yield similar taphonomic signatures.
Paleofaunal data collected as part of the current study further corroborate the VMB formative hypothesis outlined above. VMB assemblages are characterized by the same general taxonomic groups preserved in comparable proportions, with one conspicuous exception (Figs. 7, 8, Table 8): the UC-913 surface collection preserves a relatively disproportionate abundance of gar scales (Fig. 7). The challenge is whether to interpret this faunal anomaly in UC-913 as an artifact of taphonomic bias related to differential shape sorting (e.g., Blob and Fiorillo Reference Blob and Fiorillo1996) or as true paleoecological signal. Our sieved subsamples provide a baseline to disentangle paleoecological signals from preservational bias because: (1) they are less susceptible to collection bias, and (2) they include an abundance of small bioclasts ranging down to 500 µm that would have been especially sensitive to abiotic sorting processes prior to final burial. If size and shape characteristics in sieved collections are taphonomically equivalent, faunal patterning in corresponding surface collections should not be ascribed to preburial size and shape sorting. Surface processes may still be a significant biasing factor, but it is unlikely that surface processes alone would bias faunal composition to the extent exhibited by UC-913. Moreover, surface-process biases would likely impact multiple sites from a given field area similarly.
From a taphonomic perspective, site UC-913 is similar to most other sites included in this study. Thus, rather than interpret the preponderance of gar scales in UC-913 as a taphonomic bias that reflects differential shape sorting, we contend that it is better construed as a reflection of original abundance in the sampled environment, which following the VMB formative model of Rogers and Brady (Reference Rogers and Brady2010), would extend back to a lacustrine source bed. In the sampling window presented by UC-913 on outcrop today, gar scales are common for reasons that cannot be readily explained by taphonomic biasing. This conclusion is consistent with counts of plate-shaped turtle shell fragments in the UC-913 surface collection, which are in line with the abundance of turtle shell documented in other sites.
The finding that one of our sites preserves an unusual abundance of a particular taxon is consistent with how we envision biological remains accumulating in VMB localities. The Late Cretaceous wetlands that served as initial depositories for the Judith River VMB fauna were complex environments, with microhabitats and species distributions reflecting onshore–offshore patterning linked to water depth, substrate characteristics, vegetation zonation, and potentially facilitative ecological interactions (e.g., Bronstein Reference Bronstein2009; Nell et al. Reference Nell, Frederick, Mazzotti, Vliet and Brandt2016). Animal populations that inhabited these settings would have partitioned habitat just as extant animals do today. Terrestrial animals such as dinosaurs and most mammals and lizards would have frequented the margins of wetlands, foraging along the banks and in the shallows. Their skeletal remains should be more abundant in facies that accumulated nearer to the water’s edge. Nesting birds may have frequented lake margins where proximity to formidable predators such as crocodilians may have limited other nest-predation pressures (Nell et al. Reference Nell, Frederick, Mazzotti, Vliet and Brandt2016). Along these same lines, small fishes and amphibians, along with juveniles of other groups (crocodilians, champsosaurs, turtles) would have likely inhabited the nearshore shallows in response to habitat preference and predation pressures, and their skeletal record should be distributed accordingly. In contrast, the fossil remains of large fishes and other animals that prefer open water would be expected to be more abundant in more offshore facies.
When the ancient depositional environments that host VMB concentrations are viewed in this context, faunal variations are expected, even within the confines of individual localities (cf. Carrano et al. Reference Carrano, Oreska and Lockwood2016). The challenges are establishing a taphonomic framework that permits the management of bias and implementing a sampling strategy that targets potential paleoecological patterns.
Conclusion
In this comparative taphonomic study of vertebrate microfossil bonebeds, we explored potential biases relating to collection methodology and paleoenvironment and evaluated the taphonomic equivalence among sites. Marked distinctions in bioclast size and shape distributions were documented in comparisons of surface-based and sieved collections. These were interpreted to reflect collector-imposed bias on outcrop, compounded by surface processes that winnowed and sorted exposed fossil debris. Sieved collections, which are effectively immune to these biasing effects, afford the best opportunity to compare sites and evaluate taphonomic equivalence. Indeed, very few significant distinctions in size, shape, and faunal composition are evident among sieved subsamples, regardless of facies context.
Our study was isotaphonomic by design (sensu Behrensmeyer Reference Behrensmeyer1988a, Reference Behrensmeyer1991; Behrensmeyer and Hook Reference Behrensmeyer and Hook1992; Behrensmeyer et al. Reference Behrensmeyer, Kidwell and Gastaldo2000), using comparisons of sites from similar preservational contexts to identify and ultimately manage potential taphonomic biases. We also addressed “strict” isotaphonomic comparability (sensu Blob and Fiorillo Reference Blob and Fiorillo1996; Moore and Norman Reference Moore and Norman2009) by evaluating taphonomic equivalence with statistical tests. Despite the various appeals of strict isotaphonomy, a statistically significant difference may not be a sufficient reason to abandon meaningful comparisons in the absence of consensus on how similar sites must be to be considered isotaphonomic (Blob and Fiorillo Reference Blob and Fiorillo1996). Reaching this consensus will be no simple task, because taphonomic attributes vary to some degree in all fossil collections, even within the same general depositional setting or the same locality. Some degree of inherent variability is expected among localities, especially when large data sets are in play, and some taphonomic variables will likely vary more than others.
Finally, the case for taphonomic equivalence among the majority of Judith River VMBs is strong, whether sites are evaluated in a general respect (e.g., Figs. 5, 6) or tested in a “strict” isotaphonomic sense. The finding that most Judith River VMBs are similar from a taphonomic standpoint is consistent with the formative model proposed by Rogers and Brady (Reference Rogers and Brady2010), which argued that channel-hosted VMBs were reworked from preexisting pond/lake VMBs and thus share taphonomic history. The taphonomic comparability of Judith River VMB localities also bodes well for future studies of paleoecology. With a comparative taphonomic framework established and outliers identified and excluded (such as UC-914), any true paleoecological signals (such as the abundance of gar in UC-913) can be disentangled from preservational bias. The rich faunal records archived in Judith River VMBs can now be explored and deciphered in a comparative sense to reveal the diversity and community structure of Late Cretaceous wetland ecosystems. The next logical step is to test the approach advanced here in other VMB-bearing units (e.g., Hell Creek Formation) to ascertain whether the Judith River record provides a broadly applicable model of VMB accumulation or instead represents a unique preservational window into life and death in ancient freshwater ecosystems.
Acknowledgments
This research was supported by grants from the National Science Foundation (EAR-1052673), the Bureau of Land Management (L10AC16281, NLCS), and Macalester College. The Bureau of Land Management provided permits and logistical support in the Upper Missouri River Breaks National Monument (special thanks to L. Anderson, L. Eichhorn, Z. Fulbright, J. Harkson, G. Liggett, J. Mitchell, C. Rye, and G. Smith). K. Moffett (Macalester College) is acknowledged for his design and construction of the automated sieving apparatus that so greatly advanced fossil recovery. We thank H. Boekenheide, B. Faulkner, J. Fuentes, S. Jabo, P. Kroehler, A. Lawrence, M. Marshall, M. Miller, A. Murch, M. Pinsdorf, P. Sullivan, A. Telfer, J. Thole, K. Wayrynen, G. Wesley-Hunt, and Smithsonian NMNH FossiLab volunteers (especially M. Fisher) for their efforts in collecting and sorting vertebrate microfossil bonebeds. We also thank A. K. Behrensmeyer, M. Brady, D. Brinkman, D. Eberth, S. Kidwell, M. Kowalewski, and R. Terry for discussions and insights related to the accumulation and analysis of fossils in alluvial systems. We thank the three anonymous colleagues who provided detailed critical reviews of this article, which allowed us to make significant improvements. The generous hospitality of the people of Winifred, Montana, is also sincerely appreciated.
Supplementary Material
Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.j626h



