In bsaul/elktoeChemistry: Data and code for Elktoe Chemistry project

\doublespacing

library(elktoeChemistry)
library(gt)

Introduction

Laser ablation coupled mass spectrometry (LA-ICP-MS) provides precise elemental analysis of shell and other tissues [@pozebon2014review]. Trace element concentrations obtained from LA-ICP-MS are actively being used to study biological processes [@pozebon2014review]. Studies often compare elemental concentrations between groups to examine physiological responses to an experiment or in situ environmental factors. Ideally, to characterize as much variation in tissue concentration each specimen would have both pre- and post-experiment LA-ICP-MS samples. However, LA-ICP-MS is destructive and anatomical regions of interest are often internal. Consequently, accurately identifying sources of variation within an individual may be impossible in most organisms, even when an experimenter manipulates an instrument designed to elicit a biological response. In addition, natural variation between individuals may be poorly accommodated by the study design.

Adequately blocking sources of variation can be a challenge in any experimental design or analysis. Indeed, after conducting a study intended to block many sources of variation from LA-ICP-MS samples, we found multiple post-experiment obstacles in the goal of interpreting differences in elemental concentrations. In this paper, we discuss methodological considerations in designing and analyzing studies of biological processes using LA-ICP-MS data. We motivate our discussion by presenting results from a completely randomized experiment designed to assess factors associated with health status of the federally endangered freshwater mussel Appalachian Elktoe (Alasmidonta raveneliana).

Freshwater mussels (Bivalvia: Unionidae, Magaritifera) are predominately suspension feeding bivalves that inhabit streams, rivers and lakes. Unionids provide a range of ecosystem services [@vaughn2017ecosystem]. They serve as food for raccoons, muskrats and other species that inhabit riparian habitat, and when abundant help stabilize streambeds. While feeding unionids filter the water column removing algae, bacteria and detritus and process and sequester some contaminants [@gagne2007changes; @vaughn2017ecosystem]. Unionid shells support a microecosytem of bacteria, protozoa and fungi that play a prominent role in stream ecology. Their shells are similar to tree rings, archiving some aspects of their environmental history, and have been used as bioindicators of pollution [@markich2002freshwater; @brown2005freshwater; @shoults2014resolution; @wilson2018freshwater].

North American freshwater mussel species are imperiled throughout their home range. The majority are federally or state listed as threatened, endangered or of special concern [@williams1993conservation]. Environmentally degrading riparian land-use practices have been attributed to their decline. Exposure to certain heavy metals from municipal effluents and other sources have been shown to increase morbidity and mortality in freshwater bivalves [@naimo1995a-review]. @hartmut2007declining attributed declining populations of the freshwater pearl mussel (Margaritifera margaritifera) to heavy metals and pesticides. Exposure to heavy metals have also been associated with genotoxic effects [@black1996dna; @khan2019bioaccumulation] enzymatic function, and behavioral changes [@naimo1995a-review].

The Appalachian Elktoe (Alasmidonta raveneliana) (Bivalvia: Unionidae) is native to the streams of eastern Tennessee and western North Carolina. The Elktoe is listed as federally endangered and remaining populations are found in the Nolichucky River in eastern Tennessee and the Upper Tennessee River system in North Carolina. Although once abundant in the Little Tennessee River, precipitous declines were noted in A. raveneliana populations after major hurricane-related flooding associated with Hurricane Frances in 2004. However, within the same watershed’s A. raveneliana Tuckasegee River populations appeared stable. @miller2013concentrations evaluated sediment concentrations of Cu, Cr, Ni, and Zn in the Little Tennessee and Tuckasegee. When considered with the chemical concentrations of A. raveneliana shells measured by @jarvis2011water, they found these elements (Cu in particular) unlikely to be the primary culprit in the A. raveneliana decline. Reasons for the decline remain unclear [@jarvis2011water; @miller2013concentrations; @miller2016potential; @pandolfi2016effects].

To determine how river-specific factors may be associated with health and physiological outcomes, we conducted a completely randomized experiment by assigning eighty-two A. raveneliana along with an equal number of hatchery-raised Wavy-Rayed Lampmussel (Lampsilis fasciola) to sites within each river. We hypothesized that the distributions of one or more elements accrued in valve layers may differ between experimental sites. Despite obtaining concentrations from 20 trace elements along multiple LA-ICP-MS transects from each specimen and carefully aligning chemistry data with anatomical features on the valves, we found hurdles to rigorously analyzing the data. The considerations we address include alignment of chemistry data and anatomic features, noisy and skewed data where extreme values may be important, censoring of physiological processes due to mortality, stability of results to data pre-processing and denoising, and lack in knowledge of biological mechanisms by which to explain observed differences.

Materials and methods

Study Specimens

Eighty-two A. raveneliana specimens were collected from a single location in the Tuckasegee river. Ten were randomly selected and sacrificed for baseline comparisons, and the remaining 72 were randomly assigned and relocated in cages to six experimental sites, three sites each in the Tuckasegee and Little Tennessee Rivers (see appendix for site locations). The same number of young adult L. faciola specimens, obtained from the NC Wildlife Resources' Marian hatchery, were randomly assigned to a separate cages at the six sites. The plastic mesh cages were anchored in the bottom substrate. All animals were individually identified with randomly selected numeric identification tags glued to the shells. Weight and other metrics were obtained before and after deployment in the cages. After the experiment, a small (approximately 0.5 inch) section of the right valve from each animal was removed and embedded in epoxy for laser ablation analyses. Four specimens were removed from analysis due to anomalies in their valve anatomy or in the laser transects.

LA-ICP-MS

Since the study was initiated to identify factors potentially contributing to the decline of the Elktoe populations, our analysis focused on the sections of shell that biomineralized during the period of the sentinel experiment. Parallel transects from the outer margins of valves, away from the umbo were targeted for analysis. For thicker valves we also analyzed additional offset scans to evaluate reproducibility and include older growth in the prismatic layers. In most cases, the entire preserved growth history was recorded by making multiple offset traverses (Figure \@ref(fig:alignment)). We used a 25$\mu$m round aperture moving at 5$\mu$m/sec. The quadruple sampling period (0.5762 s) enabled a moving compositions average for the 20 surveyed analytes every 2.88 $\mu$m (8.7 readings per 25$\mu$m footprint). For valve cross-sectional thicknesses on the order of 1 mm, transects included roughly 350 measurements. The laser trajectories were drawn perpendicular to annular growth rings in a "step-scanning" approach [@shoults2014resolution], moving from the inner to the outer shell, beginning and ending in epoxy. Scans varied depending on location, residing solely with the outer prismatic layer and periostracum, or being compound scans initiating in the nacreous layer.

Each transect resulted in a chemical time-series for each analyte (concentration vs. distance; a matrix of about 7000 concentration measurements per mm). For data reduction^[Nathan/Heather - can you explain what this means?] we assumed a stoichiometric CaCO3 composition of 39.547 wt\% Ca in nacreous tissue and 40.0432 wt\% Ca in other tissue and did not adjust for mineralogy (aragonite vs. calcite){>>Include EPMA work done to evaluate Ca concentrations of the nacreous and prismatic layers. The nacreous layer compositions would need a slight correction for the difference in Ca concentration (used as the internal standard).<<}. We used USGS MACS-3 (synthetic aragonite) as the primary calibration standard, and analyzed NIST 614 (trace elements at ~1 ppm) and NIST 612 (trace elements at ~35 ppm). For drift compensation, each hour of data collection on unknowns was bracketed by analyses of each standard (in duplicate). We surveyed the shells for twenty industrial metals and bio-metals: Mg, Ca, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Kr (for Se interference correction if detected), Sr, Mo, Cd, Sn, Ba, Hg, Pb and U. Analyses were conducted at the University of Texas at Austin using an ESI New Wave Research NWR193 excimer laser coupled to an Agilent 7500ce ICP-MS.

Anatomical Alignment

Review of high resolution cross section images and annotation of the initial shell/epoxy boundary, annuli, nacre/prismatic layer boundary and terminal shell/epoxy boundary were made using Adobe Illustrator (Adobe Inc, San Jose CA). Annuli were identified as any linear growth feature which could be traced from its origin in the umbo to its emergence through the prismatic and periostracum layers. Non-continuous linear features along the transect were not annotated and assumed to be "false annuli" caused by environmental disturbances [@veinott1996identification].

To align anatomical features with chemistry data, several algorithms were used to first identify a reference point such as the epoxy/nacre boundary in the chemistry data. These methods included detecting changes in mean, variance [@killick2014changepoint; @killick2016changepoint] or the extrema [@borchers2018pracma] of Ca43_CPS, Pb208_CPS, or the ratio of the two (Figure \@ref(fig:alignment)). The difference of the distance between the detected reference point and the distance between the valves edges measured from high resolution scans was then computed. For each valve, the method that minimized this distance was then used to locate anatomical locations in the chemical data. The project website includes the code and complete details of the alignment steps.

knitr::include_graphics(here::here("manuscript", "figures", "figure_alignment.pdf"))

Chemical data pre-processing and de-noising

Valves were grouped into seven batches for LA-ICP-MS processing. During each batch process, the lower limit of detection (LOD) was assessed multiple times using the standards noted above. Concentration values were left-censored based on the mean LOD within a batch for each element. For analysis, parts-per-million (ppm) values were rescaled to $\text{mol/Ca mmol}$ where $\text{Ca mmol}$ was based on Ca\% of ~39.54\% for nacral observations and 40.078\% for all other layers. Chemical signals were denoised by a variety of methods including (1) no filtering, (2) a generalized additive model, and (3) a rolling median absolute deviation filter with a bandwidth of ten observations.

Study Questions

We hypothesized that the distributions of one or more elements accrued in valve layers would differ between experimental sites. We examined differences between sites within the following valve layers: (1) the youngest annuli in the nacre, (2) older nacre, (3) prismatic layer, and (4) the periostracum. The epoxy layers were also included in the analysis as a negative control. No differences should be detected in the epoxy layers.

Statistical Analyses

Specimens' weight and size measurements at the start of the experiment were summarized by site with means and standard deviations. Within each valve, observations were grouped across laser transects into the following six layers: inner epoxy the youngest nacral annuli, older nacral annuli, prismatic layer, periostracum, and outer epoxy. The distribution of elements within each layer were summarized by its central tendency. Due to the amount of values below detection levels in some elements, we estimated the first L-moment based on probability weighted moments for left censored data [@lmomco2018]. In the case that all values for a specimen were below detection within a layer, the detection limit was used as the layer's summary statistic. The L-moments are computed from order statistics and are less sensitive to outliers compared to standard moment estimates such as the arithmetic mean.

Median and interquartile ranges of the moment estimates were tabulated by species, layer, and site. Within each species/layer combination, the null hypothesis of no difference between any site was tested by the nonparametric Kruskal-Wallis method [@hollander2013nonparametric]. To account for possible batch effects in LA-ICP-MS processing, tests used the residuals from a no-intercept linear model including a term for the LA-ICP-MS drawer. The tests were done both including and excluding the baseline site. No control for type 1 error was done. P-values are presented as descriptive measures, and no statistical significance threshold was set. Analyses were done using R version r paste(R.Version()$major, R.Version()$minor, sep = ".") [@rcore2020].

Results

Baseline characteristics of experimental specimens are shown in Table \@ref(tab:table1). In both species, specimens were reasonably balanced in terms of size.

x <- readLines(here::here("manuscript", "figures", "table1.tex"))
x <- x[!grepl("\\\\small ", x)]
x <- c(x[1], "\\label{tab:table1}", x[2:length(x)])
cat(x)

knitr::include_graphics(here::here("manuscript", "figures", "pval_summaries_1.pdf"))

Differences between rivers (and sites) within layer, species, and element are summarized by Kruskal-Wallis test p-values as a heat map in Figure \@ref(fig:pvalues). Smaller p-values are shown as darker colors, while p-values greater than 0.1 are shown as white. Tests were based on residuals from a simple linear model adjusting the first L-moment for LA-ICP-MS drawer to account for possible batch effects. Figure \@ref(fig:pvalues) show results for test of four related null hypotheses of the distributions of the adjusted L$_1$ moments: the two rivers and the baseline site are the same, any site and the baseline are the same, and the previous two tests but excluding the baseline group. In general, tests including the baseline site have smaller p-values, suggesting more discrepancies from a null hypothesis when the baseline group is included as an experimental group.

The A. raveneliana L$_1$ moment values used in the Kruskal-Wallis tests are displayed in Figure \ref(fig:Arav_summary). To aid in display, the values have been standardized within each layer/element. A common pattern in the young nacre is for the distribution of baseline specimen to be lower and have a smaller spread. Tabular summaries corresponding to Figure \ref(fig:Arav_summary) of the 1960 (2 species by 20 elements by 7 sites by 6 layers) element summary statistics by layer and site are available at the project website.

knitr::include_graphics(here::here("manuscript", "figures", "element_summaries_Arav.pdf"))

Discussion

Our study has one of the largest sample sizes for completely randomized experiment in freshwater bivalves and compared to other studies measured more trace elements. In this way, our study may useful database for future research. The study also serves to highlight several methodological challenges to drawing inferences about biological processes from LA-ICP-MS obtained element concentrations. We discuss these considerations in light of our results.

Laboratory Protocols

Figure \@ref(fig:pvalues) shows a large number of possible differences (smaller p-values) between site elemental distributions is in the outer epoxy layer. Figure \@ref(fig:Arav_summary) reveals the underlying patterns for the A. raveneliana specimens. For example, the Little Tennessee sites as group have higher median concentrations of Fe, Mn and Pb in the outer epoxy compared to the Tuckasegee sites (p-values are 0.059, 0.003 and <0.001, respectively). Since the all valves were set in the same epoxy at the same time, concentrations in the epoxy should serve as negative control. We should see no association.

Misalignment of anatomy and chemistry is one explanation for the observed associations in the outer epoxy. Some observations, for example, could be incorrectly labeled as outer epoxy but are actually periostracum. Spillover of LA-ICP-MS measurements from one sample point to the next could also play a role. Since the laser ran from inner edge to outer, spillover may result inner observations looking more like the previous observations. We do see some evidence that the composition of the inner and outer epoxy do differ; however, other causes that cannot be ruled out are chemical interactions with periostracum or potentially leaching from the valves into the epoxy or vice versa.

The periostracum is the outer valve layer exposed to the ambient environment, hence differences in the this layer could be biological relevant. Since our alignment procedure (Section \ref{}) used the inner edge of the valve as the registration point, measurement or scaling errors may be exacerbated at the outer edge of the valve. If we only saw differences in the outer epoxy, we could assume that the differences were due to mis-classification of valve layers. However, of the three elements exemplified above, Mn and Pb also have indications of potential differences in the inner epoxy. The problem of multiple statistical tests aside (see Section \ref{}), the possible mislabeling of anatomical features compounds an already difficult task of parsing the logical consequents resulting from the multitude of queries generated by comparisons of elements within layer within (river/site) within species.

Recomendation 1: Plan anatomical and chemical alignment before the experiment.

We designed the alignment procedure described in Section \ref{} after it become clear during the analysis phase our study questions could only be addressed only after resolving the anatomical locations of each LA-ICP-MS data point in order to align experimental time with anatomy. Ideally, we would have marked the specimens in some way to better identify the portion of valve that grew during the experiment (see Section \ref{biological-mechanisms} for further discussion on this point).

For our study, a discrete sequence of pre-registered ablation points (rather than continuous ablation transects) may have made alignment easier. Software such as the R sclero package [@vihtakari2016sclero] is available to aid processing for this type of data. In the case that continuous ablation transects are used and the transect begin and/or end in epoxy, spiking the epoxy with a trace element known to be rare in the biological specimen can help identify edges. This is especially important where the biological processes of interest are, as in our case, at or near the boundaries between anatomical features.

Recomendation 2: Include a negative control.

Negative controls [@arnold2016negative] are common in bench science. Analyzing the trace elements in the epoxy, which should have no association with our experimental treatment, highlighted several of the consideration we discuss here. Having a negative control for each specimen helped us in considering whether observed differences were due to the experiment or due to LA-ICP-MS or non-experimental artifacts.

Recomendation 3: Block randomize experimental treatments within the LA-ICP-MS array.

In many elements and layers the median L$_1$-moments of A. raveneliana baseline specimens are lower than the concentrations of specimens in either river (see Ni and Zn$_66$ in the young nacre and outer epoxy in Figure \@ref(fig:Arav_summary) for clear examples). The baseline specimens were sacrificed in April with the experimental specimens six months later. Hence, it could be that what we observe is simply seasonal changes in nacral annuli. However, review of the array in which specimens where processed by LA-ICP-MS shows that all the A. raveneliana baseline specimens were done in sequence, all of them within the same drawer and set of standards. The experimental specimens, on the other hand, were ordered randomly in the LA-ICP-MS array. The ordering was a mistake in how the specimens were set in epoxy, but the mistake means that we can't completely rule out that observed results are not a measurement artifact. Hence, we recommend that researchers prepare their specimens in such a way that treatments are evenly distributed across the LA-ICP-MS. A simple method would to randomize specimens within LA-ICP-MS "blocks" (i.e. between standards), including at least one specimen from each treatment strata within a block.

Biological mechanisms

In the previous section, we considered our results in light of the study design and lab protocols. In this section, we consider challenges to interpretation due to our understanding (or lack thereof) of the physiological processes.

If the experimental exposures stop or slow biological processes, then statistics from observed data may be deceiving without additional information that we were unable to collect. Consider the following scenario. An experimenter has three bivalve specimens and exposes each one to a different concentration of mercury: none, low, and high. Presume that Hg is known to be toxic and that nacral Hg concentrations have been shown to positively associate with ambient Hg concentration. The experimenter returns six months later and finds the specimen treated with high amounts of Hg dead and the other two survived. When the LA-ICP-MS concentrations of the youngest nacre are analyzed, the mean concentration of Hg are highest in the specimen exposed to low concentrations. The high-exposed and the not exposed appear to have the same mean concentration of Hg. Unbeknownst to the researcher, the high-exposed specimen died after only two weeks of the experiment. Death stopped valve accretion, hence the high Hg did not appear to be recorded in the valve. Without knowing how long the specimens survived this cannot be accounted for in the analysis.

In our experiment, A. raveneliana mortality and morbidity was increased in the Little Tennessee and increased also increased in a downstream direction in both rivers. Relevant to the scenario above, we do not know when specimens died during the experiment. Consider Mn as an example, which we were particularly interested, as it is an essential trace element with important roles in numerous enzymatic pathways that could have harmful effects on freshwater bivalves. Interestingly, surviving valves from the Little Tennessee had films on the gills that were apparent by scanning electron microscopy, which upon x-ray diffraction analysis consisted largely of Mn. Manganese valve concentrations have been shown to correlate with seasonal changes [@siegele2001manganese; @shoults2014resolution]. The reasons for the seasonal variation are unknown, but in our study, different survival times could have resulted in different accrual patterns. Those that survived longer would have kept accumulating Mn and other elements for later in the season. It could be that those specimens that survived to the end of the study (Tuckasegee 1 site had the lowest mortality) were accruing shell for longer while the trajectories for other specimens were censored at their time of death.

Recomendation 4: If survival outcomes are of interest, find ways to noninvasively monitor the health of specimens during the experiment.

We assumed that the youngest annuli in the nacre accrued during the experiment, but the data provide no way of testing this assumption. In Little Tennessee specimens, where high morbidity and mortality was observed in A. raveneliana it could be that little to no nacre accrued during the experiment. Marking the start of the experiment shells in some way may have provided appropriate data, but @haag2008testing points out that notching shells does harm individuals. Though we carefully aligned valve layers with chemistry data, in some cases layers were most certainly misclassified particularly at layer boundaries. Future studies should consider ways of safely marking the start of the experiment on the valves in order to measure and classify the valve that accrued during the experiment.

Recomendation 5: Find ways to mark the start of the experiment on the anatomy before the experiment.

Based on current understanding of bivalve shell growth [@checa2000a-new-model], the laser transects nearest the mantle edge that included nacre should have captured the accretion that occurred during the experiment. Shell layers are secreted along the mantle edge of each valve, so the valve grow out from the umbo (Figure \ref{}). By aligning of anatomic images and chemical data obtained by LA-ICP-MS, we sought to compare elemental concentrations across experimental assignments in each bivalve layer. The nacral annuli classified as the youngest which should have accumulated during the experiment, while the older nacral annuli should be pre-experiment. While the angle of laser transects cut across the recorded time in the nacre, the same angle cut a more cross-sectional sample of the prismatic layer. Although, prismatic layer observations from transects closer to the mantle edge are younger than those closer to the umbo, resolving the chronological time of the prismatic layer observations relative to the chronological time of the experiment would be challenging. Visual inspection suggests most of the prismatic layer sampled valve accrued prior to the experiment.

The periostracum is exposed to the ambient environment, hence it could be that the observed differences in the outer epoxy shown in Section \ref{} are actually represent differences in the periostracum. To our knowledge, unlike the prismatic and nacral layers, there is no existing literature on the plasticity of periostracum (how it may change after initially created) or on how the periostracum records water chemistry. This is worthy of further study, it could be possible to non-destructively sample periostracum.

Statistical Considerations

• noisy and skewed data where extreme values may be important
• stability of results to data pre-processing and denoising

The causes for A. raveneliana declines remain a mystery, as our results do not show any clear association between contaminants and experiment sites or health outcomes.

The project website bsaul.github.io/elktoeChemistry/ contains supplementary analyses and project details. Data and code for this project are available from the project's github repository github.com/bsaul/elktoeChemistry.

Funding

This work was supported by funds provided by the North Carolina Department of Transportation [NCDOT Project 2013-2012-37]; and the US Fish and Wildlife Service [F18AC00237].

Acknowledgment

The authors thank [TODO]

References

bsaul/elktoeChemistry documentation built on Nov. 17, 2022, 8:10 a.m.