Resolving taphonomic and preparation biases in silicified faunas through paired acid residues and X-ray microscopy

Geological setting

The Edinburg Formation is an Upper Ordovician unit of massively bedded black limestone with occasional shale interbeds, representing deep ramp to basinal deposition at the northwest of the Taconic foreland basin system after drowning of shallower ramp carbonate facies including the underlying Lincolnshire Limestone (Holland & Patzkowsky, 1996). K-bentonites derived from Taconic volcanism are sporadically present in the sediments of the foreland basin and have been sampled for radiometric ages, most notably the Millbrig bed, which lies up section of the Edinburg in the Martinsburg Formation and has a U-Pb date of 452.86 Ma (Mitchell et al., 2004; Sell, Ainsaar & Leslie, 2013).

All samples in this study were collected from the Liberty Hall facies of the Edinburg Formation exposed at Strasburg Junction rail cut, a well-studied site in the Shenandoah Valley of Virginia (Fig. 1, Table 1; Cooper & Cooper, 1946; Whittington & Evitt, 1953; Read, 1980; Jacobs & Carlucci, 2019). Hand sample observations revealed two starkly different lithologies (Table 1). Rocks from Horizon 1, approximately 0.5 m above the contact between the Lincolnshire Limestone and Edinburg Formation, were coarse-grained, with a sparry or dismicritic texture likely indicative of recrystallization during early diagenesis. Sparse fossils and possible intraclasts visible on fresh surfaces were weakly aligned with the direction of original bedding, suggesting directional sorting during deposition. Horizons 3, 4, 5, and 6 cropped out over 30 m up section (Table 1) and were fine-grained with only occasional sparry or hematitic grains. Fossils, where visible on fresh or weathered surfaces, were oriented randomly to original bedding and comprised arthropods, brachiopods, and occasional crinoid ossicles. Rocks from these horizons also displayed heavy rinds on weathered surfaces, frequently stained rusty orange to pale yellow with iron oxides/oxyhydroxides.

Fossil extraction

From five horizons with fossils visible on weathered surfaces, 15 samples were prepared as cylindrical cores (Table 1), approximately 2 cm in diameter and ranging from 1 to 4 cm in height, normal to bedding. Each core was imaged via X-ray tomographic microscopy (μCT) at the X-ray Microanalysis Laboratory at the University of Missouri using a Zeiss Xradia 510 Versa X-ray microscope with an isotropic voxel size of approximately 30 μm, producing three-dimensional renders of the core interiors. All μCT scans were processed using Dragonfly software Build 941–v.4.2.2 for Windows, Object Research Systems (ORS) Inc, Montreal, Canada, 2018 (http://www.theobjects.com/dragonfly). Post-image processing within the software was conducted to reduce imaging artifacts using the Ring Removal and Median filters. Regions of high brightness (corresponding to high-density ferrous material) and low brightness (corresponding to low-density siliceous material) relative to the local background were thresholded and manually segmented, excluding regions at the very top and bottom of cores due to persistent boundary artifacts. Due to radial and longitudinal variation in background brightness across the cores, multiple brightness thresholds were used to segment siliceous material in some cores; in such cases, the total siliceous volume is the union of the volumes from segmentation by the various thresholds. Segmented volumes for ferrous and siliceous material were refined by removing small (<100 voxel total volume) islands and were then partitioned into distinct regions-of-interest (ROIs), using six-connectivity (voxels considered to be connected if sharing faces rather than only edges or vertices) for both purposes. Individual fossils and other objects were identified among these ROIs by visual inspection in 3-D representations of voxels and against individual 2-D projections of μCT imagery. In cases where a fragmentary or incompletely resolved fossil fell into multiple non-connected ROIs, those ROIs were merged into one; in cases where a single ROI contained multiple contiguous objects, it was manually partitioned.

After each core had been imaged with μCT, they were macerated in 10% acetic acid to dissolve calcareous material until fully disaggregated and no longer visibly evolving gas bubbles. Acetic acid was chosen over alternatives such as hydrochloric acid or formic acid to allow for slow dissolution and avoidance of fossil damage. Insoluble residues were washed with water and sonicated in 30-s intervals, iterating until the supernatant ran clear. This washing process was then repeated using Calgon solution (0.052 M Na₆(PO₃)₆, 0.286 M NaHCO₃) as a deflocculant, sonicating as before, until all loose clay was removed. Cleaned residues were then washed over a 250-μm sieve, oven dried, and picked for identifiable fossils.

Statistical analyses

Following scanning and segmentation, each object resolved in the μCT data was measured for minimum, maximum, and mean three-dimensional Feret diameter. Fossils that could be confidently identified to at least the phylum level were scored for additional taxon-specific features and anatomical measurements (Table 2).

Fossils recovered from residues were photographed using a GIGAmacro Magnify2 Robotic Imaging System with Canon EOS Rebel T8i DSLR and Nikon T1 1× and 3× objectives. For large specimens, additional photographs were taken under a reflected light microscope (Nikon SMZ1500 tethered to a Nikon D600 DSLR) to record features not visible in top-down view. Photographs were processed and analyzed with FIJI/ImageJ software, using the Trainable Weka Segmentation plugin to isolate fossils against the image background (Andreola et al., 2004; Schindelin et al., 2012). Each fossil was measured for minimum and maximum two-dimensional Feret diameter as well as taxon-specific measurements corresponding to those taken from μCT data (Table 2).

For each core, both the segmented μCT render and acid residues were counted for the total abundance of fossil specimens. Trilobites were identified to the family level based on their general geometry, furrow pattern, and, in residues, their prosopon. Trilobite material not reliably assignable to a single family was treated as a separate category. Prosopon and other textural features were not easily resolvable in the μCT dataset, and so were not generally considered for taxonomic assignment in counts based on μCT results. Bryozoans were morphologically classified by growth form either as thin-branching, thick-branching, or fenestrate. Ostracods, gastropods, and bivalves were not further classified due their generally small sample size and coarse silicification, precluding the reliable identification of taxonomically relevant features.

Using R statistical software, datasets were subjected to NMDS in three dimensions based on the Bray-Curtis dissimilarity index, using the metaMDS function provided in the vegan package (Dixon, 2003; R Core Team, 2021). Sites with no counted fossils were excluded from the analysis. Further paleoecological analysis of taxon abundances within and between samples was performed using PAST statistical software, calculating ecological dominance within samples (Simpson’s D) and assessing compositional similarity between samples using the Bray-Curtis index (Hammer & Harper, 2001, 2022).

Tests of preparation bias

A contingency table was constructed containing μCT identified specimens counts, broken down by taxonomic grouping and by presence or absence in residues. The relationship between taxon and recovery was investigated using Pearson’s ${\chi }^2$ test (Table 4). The only taxonomic group which deviated significantly (p = 0.04) from overall likelihood of recovery was Unidentified Trilobita—material clearly from trilobites but lacking anatomical features in μCT sufficient to assign it to a family—which, despite being identified in μCT, are disproportionately unlikely to be recovered in residues.

Application of Pearson’s ${\chi }^2$ test to total counts of taxa identified in μCT and in residue found strong support (p = 1.89 × 10⁻³⁵) for different proportional abundances in the two preparation types (Table 5). Raphiophorids (p = 1.99 × 10⁻⁵), remopleuridids (p = 0.02), unidentified trilobites (p = 9.92 × 10⁻¹⁷), bivalves (p = 0.02), and gastropods (p = 2.08 × 10⁻⁴) are significantly more abundant in μCT, while thin-branching bryozoans (p = 3.26 × 10⁻¹⁴), thick-branching bryozoans (p = 1.15 × 10⁻³), and ostracods (p = 1.27 × 10⁻¹³) are significantly more abundant in residues.

The presence of size bias in the likelihood of recovery was evaluated using the Mann-Whitney U-test with a null hypothesis of no difference in medians and a standard threshold of significance of α = 0.05 (Mann & Whitney, 1947). In comparisons between fossils identified in both μCT and residue and those identified in μCT but not in residue, there was no statistically significant difference of maximum (p = 0.54), minimum (p = 0.09), or mean (p = 0.24) Feret diameters, nor of elongation factor (p = 0.16) defined as the ratio of maximum to minimum Feret diameters. However, when testing overall distributions of size between all fossils identified in μCT and all those found in residue, the U-test found a dramatic difference in maximum Feret diameter (p = 1.75 × 10⁻⁷⁹) between the medians of the two groups (3.98 mm for μCT, 1.41 mm for residue).

Exploratory paleoecological ordinations

NMDS performed separately on the μCT (Table 6) and residue (Table 7) datasets resulted in ordinations with stress scores of 0.06 and 0.08, respectively (Fig. 4). Taxon loadings in both ordinations consistently grouped bivalves and gastropods nearby each other. Ordinations of both μCT and residues tended to group medium-sized benthic trilobites (raphiophorids, pterygometopids, and cheirurids) together along with unidentified trilobite material (most of which likely derived from one of those families), with pelagic trilobites (remopleuridids) plotting closer to gastropods and bivalves. Asaphids, likely represented here by the extremely large benthic trilobite Isotelus, consistently fall near the latter cluster but are represented only by fragments of cuticle. The ordination position of agnostids, whose life habit remains controversial (Fortey & Owens, 1999), is inconsistent between μCT and residues, falling near the benthic and pelagic trilobite clusters in those analyses, respectively. In both cases, they are closely accompanied by ostracods (which can occupy pelagic or benthic niches and whose life habit was not interpreted in this study). Due to the extreme disparity between identification in μCT and recovery in residues of ostracods, their location within the ordination for μCT data is likely not informative.

Thin-branching and fenestrate bryozoans plotted relatively close together, forming an isolated cluster. Thick-branching bryozoans were not closely associated with the other two morphotypes in the residue ordination and were not detected in μCT.

Samples from individual horizons tended to fall near each other in loose association. Two main clusters are apparent in both residue and μCT ordinations: one characterized by low NMDS1 scores and dominated by bryozoans, and one with higher NMDS1 scores and dominated by arthropods (trilobites and ostracods).

For each sample, dominance was computed for both μCT (Table 6) and residue (Table 7) taxon totals, with 95% confidence intervals based on 9999 bootstrap replicates (Fig. 5). The results of this analysis differed notably between the two; μCT dominance values fell in a tight band between 0.25 and 0.40 (excepting two outliers with n < 3), while residue dominance values had a bimodal distribution with one group of values between 0.40 and 0.65 and another falling between 0.15 and 0.30 (Fig. 5). Rarefaction curves were calculated for both μCT (Table 6) and residue (Table 7) counts of each sample using PAST (Hammer & Harper, 2022; Fig. 6). Curves were visually inspected for the presence of an inflection point as a rough qualitative assessment of sampling completeness; a sample’s curve “leveling off” (sudden decrease in slope) is considered informal evidence that taxa in the true population are well-represented in the sample (Sanders, 1968; Raup, 1975). This inflection point was observed in residue counts from most samples (prominently in 03-C1, 03-C3, 04-C1, and 06-C3), but not in μCT counts. This suggests that the smaller total counts of individuals in μCT data are leading to the non-recovery of rarer taxa. Such under-sampling is a liability of ecological analyses based on those samples.

Cluster analysis

Multivariate cluster analyses were performed in PAST to provide alternate metrics of similarity in the paleoecological composition of the samples, using the Bray-Curtis index as before (Bray & Curtis, 1957; Hammer & Harper, 2022). Dendrograms with branch length scaled to similarity were constructed for μCT counts, residue counts, and a combined dataset treating μCT and residue counts for each sample as two separate sites (Figs. 7 and 8). Cluster analysis of μCT counts grouped samples from Horizon 04, and treated Horizon 05 similarly, but scattered samples from Horizons 03 and 06 across the tree. Using residue data, Horizon 05 formed a cluster as before but Horizon 04 split, with cores 04-C2 and 04-C3 remaining close but 04-C1 further removed; Horizons 03 and 06 were dispersed as in the μCT tree. In both trees, samples from Horizon 01 were far removed from other samples, tending to form a cluster basal to the rest of the samples.

The combined tree retains the general close groupings of the two individual trees due to sharing the same dissimilarity metric (Fig. 8). Notably, even for Horizons 01, 04, and 05 whose samples tended to cluster in both the μCT and residue trees, those clusters are not closely related in the overall analysis; there are few clear patterns in this tree, but μCT samples tend to cluster more closely with other μCT samples, and residue samples with residue samples, than do μCT and residue counts of the same sample or even of samples from the same horizon.

Size and taxon effects on recovery

We hypothesized at the outset of this study that body size would be negatively associated with likelihood of recovery in residues; larger fossils may be more susceptible to breakage, and larger organisms with thicker skeletal elements may not fully silicify, producing a brittle outer husk. However, the results of our analysis do not support this possible effect. While fossils recovered in residue have a smaller median body size than those resolved in μCT, this appears to be due not to preferential breakage of larger fossils during washing but rather disproportionate non-resolution of smaller fossils in μCT datasets. Rather than μCT acting as a baseline to test the biases of acid maceration based on body size, the results of this study suggest the opposite.

Overrepresentation with respect to the abundance of pyritized taxa (gastropods and bivalves) in µCT may be the result of the high contrast between iron-bearing phases and the matrix, making such fossils easier to resolve in μCT than otherwise-comparable siliceous fossils. Alternatively, because most of the pyritized fossil materials represent internal molds, they may be harder to identify outside of the 3D context provided in μCT, where they remain within their shell. If the latter is the case, this likely reflects lower fidelity of moldic pyrite preservation than that of the replacive silicification. Results of the NMDS ordinations consistently plotted molluscan taxa nearby one another, which may reflect either a shared infaunal environment or instead early burial conditions conducive to pyritization in those deposits. However, the low abundance of bivalves (2 in μCT, 1 recovered from residues, always alongside multiple gastropods) makes this apparent association tenuous. It is noteworthy that disseminated pyritization in the form of non-fossil granules was present in several samples (e.g., 05-C4) that contained no mollusks, suggesting that the presence or absence of mollusks in a sample may reflect their abundance in the paleoenvironment and is not solely controlled by whether early burial chemistry allowed pyritization to proceed. Ostracods, meanwhile, are disproportionately overrepresented in residue due to their small size making them difficult to resolve in μCT at the resolution used in this study. Branching bryozoans also suffer from this, but their numbers in residue are likely inflated by fragmentation of large individuals into many smaller ones.

Nonrecovery bias and potential causes

Revealing that the number of unidentified trilobite fossils recognized in μCT is significantly less likely to be recovered in residue is to be expected. Fragmentation and incomplete preservation are common causes for trilobite material to lack identifying characteristics, and both are likely to promote further breakage and degradation during washing by compromising the structural integrity of the sclerite. This may also influence likelihood of recovery by impeding identification of the fossil in residue; a fragment of cuticle lacking identifying characteristics is less likely to be recognized as corresponding to an object observed in μCT, and minor breakage is likely to disrupt recognizable aspects of its outline and other key features.

Diagenetic biases arise from a variety of sources and are heavily dependent on shell microstructure, organic matter, and availability of reactive and replacive ions to the shell. It is well-documented that certain taxa are more susceptible to silicification than others, and, even within taxa, textural differences can make certain skeletal elements more likely to be preserved than others (Daley & Boyd, 1996; Cherns et al., 2011; Butts, 2014). Organic matter content within a shell also influences the likelihood of silicification by providing nucleation sites for the deposition of silica, though the chemistry of this is complex (Wallace, DeYoreo & Dove, 2009; Butts, 2014). Similarly, replacement via pyrite tends to initiate at sites of organic matrix in the shell but is also dependent upon reducing microenvironments and/or microbial zonation within the sediment (Fisher, 1986; Canfield & Raiswell, 1991; Schiffbauer et al., 2014).

Further biases can arise from the sample extraction and preparation process itself. During washing, larger fossils may be tumbled against other grains or container walls, leading to abrasion or breakage, while smaller fossils may be crushed beneath larger objects as grains settle. Elongated grains may be more prone to breakage than spheroidal grains of the same volume due to a lower minimum cross-sectional area, which can lead to taxon or sclerite bias against fossils with rod-like geometry, such as branching bryozoans and the spines of some trilobites. Breakage during preparation can cause taxa to be undercounted due to indiscriminate destruction of individuals or damage sufficient to remove recognizable features. Paradoxically, it can also inflate taxon counts by turning one fossil into many still-recognizable fragments. These overlapping effects can produce a wide variety of biases in size and shape, and the cumulative effect on measurements of abundance cannot reasonably be predicted a priori.

The composition of the host rock can further impose preparation-related biases. For instance, well-cemented rocks frequently require longer maceration periods than poorly consolidated ones, while argillaceous rocks may need more thorough washing to drive off insoluble clays, often including sonication. The longer sediment is washed and manipulated, the more breakage tends to occur, making matrix texture an important factor of preparation bias. While the acid maceration processes can be accelerated by using alternative acids, such as hydrochloric, the violent effervescence caused by the intensity of the reaction can also have a deleterious effect on the extraction of delicate forms.

Artifacts and limitations of µCT

While a powerful tool for visualization, μCT does present some unavoidable sources of potential error: artifacts arising from specimen capture and processing, particularly the “hardening” of the X-ray beam through preferential absorption of lower-energy photons by the sample. Beam hardening is a ubiquitous issue in μCT, but not all hardening artifacts are of equal impact. Since the beam is slightly hardened by passing through the outer surface of the sample, the interior is always slightly darkened relative to the outermost layer; these “cupping artifacts” are produced predictably based on sample geometry and can be corrected relatively easily (Schladitz, 2011; Jung et al., 2011; Abel, Laurini & Richter, 2012).

Heterogeneous samples may experience additional artifacts. Regions of higher-density material (e.g., iron minerals in sediments, bones in biological samples) can drastically harden beams passing through them, causing dark blotches to appear around bright features within a sample and especially in the spaces between multiple bright features. These artifacts are far more difficult to correct due to their irregular shape and remain the subject of ongoing research in economic geology and other materials science fields (Remeysen & Swennen, 2006; Park, Chung & Seo, 2015; Bam et al., 2019). Given the abundant pyrite in the material examined in this study, this stands as a caveat to the results reported here; similar methods applied to rocks with lower density variation may yield clearer μCT data and therefore more complete identification of fossils within the samples. Filtering during pre-processing, or correction applications in post-processing software packages, can reliably improve these artifacts, but often at the cost of introducing noise along the axis of sample rotation (Kyriakou, Prell & Kalender, 2009; Yousuf & Asaduzzaman, 2009).

Relevance for paleoecological interpretation

Broad categories of biofacies are generally consistent between residue- and μCT-derived taxon abundances as interpreted using NMDS. While total counts for μCT are much lower than those of residues, relative abundances in the former largely recapitulate the latter with some notable exceptions; ostracods are almost absent from μCT data, due to resolution constraints imposed from the diameter of the cores, even when highly abundant in residues, and bryozoan counts tend to be much higher in residues than in μCT due to breakage transforming one large fragment into many small fragments. For this reason, the application of μCT data alone to describe paleoecological structure should be undertaken with caution. Ordinations used in paleoecology are sufficiently complex that such biases in recovery can, by introducing error and uncertainty upstream, irreparably taint conclusions drawn downstream.

Overall, the use of μCT counts on their own to characterize silicified remains in these limestone biofacies is of uncertain value. Sample volume constraints make bulk samples unfeasible, and the cost of instrument time can be a roadblock to replication and large-scale sampling. However, μCT may be best suited when chemical preparation is dangerous or impractical due to matrix or fossil mineralogy, or when studying sponges, bryozoans, corals, or other modular organisms for which fragmentation during washing can obscure or inflate the number of individuals present. Samples best suited to the approach outlined in this study will be densely fossiliferous (mitigating the limitation of sample volume), bear taxa with morphological features conducive to identification in μCT (related to shape rather than texture) and contain grains more susceptible to dissolution than the matrix or cement joining them together. Such rocks should also have a contrast in density, whereby the material of interest has a higher X-ray attenuation, for better visualization, whilst keeping in mind the potential effects of beam hardening artifacts. Conodonts and other phosphatic fossils, for example, are occasionally known from silicified sediments, and characterizing taxonomic assemblages using μCT techniques may constitute an alternative to digestion from the toxic and hazardous hydrofluoric acid normally used to extract fossils from siliceous cements (Green, 2001).