Monitoring nearshore ecosystem health using Pacific razor clams (Siliqua patula) as an indicator species

Study organisms

Pacific razor clams were collected in July 2015 and 2016 at intertidal sites within each of two national parks in southcentral Alaska: Lake Clark National Park and Preserve (sites: Polly Creek, Silver Salmon and Chinitna Bay) and Katmai National Park and Preserve (sites: Hallo Bay, Swikshak and Kashvik Bay) (Fig. 1) (National Park Service permit numbers CF-15-088 and CF-16-089). In 2015, Pacific razor clams were collected from 2 sites at KATM and 3 sites at LACL and in 2016, from 3 sites at each park. Due to logistical constraints, in 2015 we were unable to collect razor clams from Hallo Bay in KATM or take morphometrics in the field on razor clams from Chinitna Bay and Polly Creek in LACL. Approximately 20 Pacific razor clams were collected from each site, with 10 designated for biomarker assays and 10 for gene transcription (Tables S1 and S2). In LACL, beaches sampled were those known to be used for recreational clamming, ensuring presence of razor clams. In KATM, several sites were identified a priori based on historical razor clam distribution (Kaiser & Konigsberg, 1977) and initial site visits were made to determine presence of razor clams. All sites sampled had razor clams present and all razor clams were collected on the morning rising tide.

Pacific razor clams were processed as soon as possible following collection, generally within 1–3 h. For biomarker assays, hemolymph was removed by inserting a tuberculin needle into the anterior adductor muscle, aspirating the hemolymph and transferring it to a cryovial. The hemolymph and razor clam were frozen in liquid nitrogen for transport to the lab where they were transferred to a −80 °C freezer until processing. For molecular assays (gene transcription), gill tissue was removed and placed in cryovials with RNAlater^®. The samples were kept at –20 °C until used in experiments. The length of collected razor clams ranged from 28.5 to 158.3 mm and age, as determined by counting annuli, ranged from 2 to 9 years. The maximum recorded length for S. patula in Alaska is 304 mm (Foster, 1991) and age is 18 years (ADF&G, 2020)

Invertebrate biomarkers

The seven invertebrate biomarker assays detailed below were selected to assess the physiological status of the razor clams as an indicator of nearshore ecosystem health (Table 1) and have been used previously in our lab (Counihan, 2018; Counihan et al., 2019). The methods are based on published research on Mytilus spp.; however, these methods were validated for use with razor clams in our laboratory prior to conducting this project (K. Counihan, 2014, unpublished data).

Morphometrics

Frozen razor clams were thawed on ice and their length, width, height and total wet weight were measured prior to dissection. The posterior adductor muscle and digestive gland were excised for assays. The remaining soft tissue was discarded and the shell weighed. The condition factor was calculated by dividing the total razor clam weight by the shell length. Soft tissue weight ratio was calculated by dividing the soft tissue weight by the total razor clam weight (Gagné et al., 2008).

Shell thickness

Razor clam shells were dried in a 37 °C incubator overnight and then a micrometer was used to measure the shell thickness at five regularly spaced points around the shell, approximately 1 mm from the edge (Versteegh & Hansson, 2012).

Hemocyte count

A 0.01 ml sample of hemolymph was diluted 1:1 with tris-buffered saline (TBS) and the number of cells were counted using a hemocytometer in three replicates (Akaishi et al., 2007).

Hydrogen peroxide production

Hemolymph samples were diluted 1:1 with TBS solution. Samples were tested in triplicate, if enough hemolymph was available, in a spectrophotometric assay by pipetting 0.05 ml of the 1:1 hemolymph:TBS solution into a 96-well plate. After incubating the plate in the dark for 1 h, 0.05 ml of phenol red solution (phosphate buffered saline pH 7.4, 5.5 mM dextrose, 0.56 mM phenol red, 8.5 U ml⁻¹ horseradish peroxidase, type II) was added to each well and incubated for another 30 min in the dark. The reaction was stopped by adding 0.01 ml of 1 N NaOH, and the plate was read on a Molecular Devices SpectraMax Plus microplate reader (Sunnyvale, CA, USA) at 620 nm (Akaishi et al., 2007).

RNA:DNA ratio

Half of the posterior adductor muscle was homogenized with a mortar and pestle and 0.15 ml of 1% sarcosyl tris-EDTA (STEB) was added to the homogenate. The mixture was vortexed for 60 min to completely dissolve the tissue. Each sample had 1.35 ml of tris-EDTA (TE) buffer added. The sample was centrifuged 15 min at 14,000×g at room temperature and the supernatant saved for testing. The samples were diluted 1:20 and 0.075 ml of each sample was added to a 96-well plate in duplicate. Genomic, unsheared DNA from calf thymus (Sigma–Aldrich, St. Louis, MO, USA) was used to prepare a DNA standard curve (0.1 μg ml⁻¹, 0.2 μg mL⁻¹, 0.4 µg ml⁻¹, 0.8 µg ml⁻¹, 1.6 µg ml⁻¹, 3.2 µg ml⁻¹, 6.4 µg ml⁻¹, 10.0 µg ml⁻¹). RNA from bovine pancreas (Sigma–Aldrich) was used for a RNA standard curve (0.4 µg ml⁻¹, 0.8 µg ml⁻¹, 1.6 µg ml⁻¹, 3.0 µg ml⁻¹, 6.0 µg ml⁻¹, 8.0 µg ml⁻¹, 12.0 µg ml⁻¹ and 16.0 µg ml⁻¹). The stock solutions of DNA and RNA were diluted in 0.1% STEB to the appropriate concentrations and 0.075 ml of each standard was added to the plate. The wells had 0.075 ml of ethidium bromide solution (2 µg ml⁻¹) added to them and the microplate was shaken for 15 min. The plate was read in a SpectraMax Gemini EM fluorescent microplate reader (Molecular Devices) with 525 nm excitation, 600 nm emission to determine the total nucleic acid reading. Each well had 0.0075 ml of RNase solution (20 U ml⁻¹) added to it and the plate was shaken for 20 min and read on the microplate reader again with the same settings. The second reading was the DNA only reading. The RNA content was determined by subtracting the second reading from the first (Caldarone et al., 2001).

Cytochrome P450 activity

The digestive gland was homogenized in buffer (25 mM Hepes, 125 mM NaCl, 0.1 mM EDTA, 0.1 mM dithiothreitol) at a 1:5 weight:volume ratio. The mixture was centrifuged at 1,500×g for 10 min at 2 °C and the supernatant transferred to a clean tube. The supernatant was centrifuged at 10,000×g for 20 min at 2 °C. The supernatant was discarded and the pellet resuspended in 0.15 ml of microsome buffer (25 mM Hepes, 140 mM NaCl, 1 mM KH₂PO₄). A microsome sample was incubated with 30 uM rifampin for 48 h to be used as a positive control. A 96-well plate was inoculated with 0.05 ml of each microsome sample in triplicate, and a media only blank was included. Each well had 0.05 ml of 50 µM BFC solution (7-benzyl-4(trifluoromethyl) coumarin in phenol red-free Dulbecco’s Modified Eagle Medium) added and the plate was incubated for 4 h at room temperature. The plate had 0.04 ml of stop solution (80% CH₃CN, 20% 0.5 M Tris base) added to each well and it was read on a SpectraMax Gemini EM fluorescent microplate reader (Molecular Devices) at 410 excitation, 530 emission (Mensah-Osman et al., 2007).

Heat shock protein

Half of the posterior adductor muscle was homogenized and lysed in buffer (150 mM NaCl, 1% Triton-X 100, 0.5% sodium deoxycholate, 0.1% sodium dodecyl sulphate, 50 mM Tris, 1 mM phenylmethylsulfonyl fluoride). The homogenate was centrifuged at 12,000 rpm for 20 min at 4 °C and the supernatant collected. The protein concentration was determined using a Bradford assay. The sample was mixed 1:1 with 2X Laemmli buffer, boiled at 100 °C for 5 min and loaded on a SDS-PAGE gel. Each gel was run with a positive heat shock protein control and a molecular weight marker. The gel was electrophoretically transferred onto a polyvinyl difluoride (PVDF) membrane. Total protein was visualized by staining the gel with Pierce reversible protein stain (Thermo Scientific, Pittsburgh, PA USA). The stain was removed and the membrane was probed with a mouse anti-HSP40 primary antibody (Abcam, Cambridge, MA USA). The membrane was washed and incubated with a secondary alkaline phosphatase labeled anti-mouse antibody (Abcam). The membrane was washed again and alkaline phosphatase substrate added. The membrane was photographed to document bands and analyzed with Image Studio Lite software, version 5.2.5 (Li-Core, Lincoln, NE, USA).

Invertebrate gene transcription

Five genes of interest and two reference genes were quantified using the quantitative polymerase chain reaction (qPCR) technique to assess the status of the razor clams as an indicator of nearshore ecosystem health (Table 2). The transcript panel utilized here was developed based on studies of clams subjected to single or multiple stressors.

Genes were selected for the transcript panel based on their known responses to stressors, including pathogens, xenobiotics, ocean acidification, heat stress and dissolved oxygen levels (Choi et al., 2006; Choi, Jo & Choi, 2008; Cellura, Toubiana & Roch, 2007; Gao et al., 2007; Perrigault, Tanguy & Allam, 2009; Li et al., 2011; Chen et al., 2012; Niu et al., 2013; Leite et al., 2013; Wang et al., 2013).

Tissue collection and RNA extraction

Gill tissue was collected from each razor clam and placed into RNAlater^® (Ambion/Life Technologies, Grand Island, NY, USA). All tissue samples were stored at −80 °C. Total RNA was extracted from pulverized gill tissue using the RNeasy Lipid Tissue Mini Kit (Qiagen; www.qiagen.com). To remove contaminating genomic (g)DNA, the spin columns were treated with 10 U μl⁻¹ of RNase-free DNase I (DNase, Amersham Pharmacia Biotech Inc.; www.apbiotech.com) at 20 °C for 15 min. RNA was then stored at −80 °C pending further analyses.

cDNA synthesis

A standard cDNA synthesis was performed on 2 μg of RNA template from each razor clam. Reaction conditions included 4 units reverse transcriptase (Omniscript, Qiagen, Valencia, CA, USA), 1 μM random hexamers, 0.5 mM each dNTP, and 10 units RNase inhibitor, in RT buffer (Qiagen, Valencia, CA, USA). Reactions were incubated for 60 min at 37 °C, followed by an enzyme inactivation step of 5 min at 93 °C and then stored at −20 °C until further analysis.

Primer design

Degenerate primers were designed based upon multi-species alignments (GenBank). Briefly, degenerate primer pairs developed for the razor clam were used on cDNA from three randomly selected razor clam samples. Degenerate primer pairs were designed to amplify five genes of interest and two ribosomal housekeeping genes (Tables 2 and 3). The PCR amplifications using these primers were performed on 20 ng of each cDNA sample in 50 μl volumes containing 20–60 pmol of each primer, 40 mM Tris-KOH (pH 8.3), 15 mM KOAc, 3.5 mM Mg (OAc)2, 3.75 μg ml⁻¹ bovine serum albumin (BSA), 0.005% Tween-20, 0.005% Nonidet-P40, 200 μM each dNTP, and 5U of Advantage^® 2 Taq polymerase (Clontech, Palo Alto, CA, USA). The PCR was performed on an MJ Research PTC-200 thermal cycler (MJ Research, Watertown, MA, USA) and consisted of 1 cycle at 94 °C for 3 min, and then 40 cycles at 94 °C for 30 s, at 60 °C for 30 s, and 72 °C for 2 min, with a final extension step of 72 °C for 10 min. The products of these reactions were electrophoresed on 1.5% agarose gels and resulting bands visualized by ethidium bromide staining. Definitive bands representing PCR products of a predicted base pair size of the targeted gene were excised from the gel, and extracted and purified using a commercially available nucleic acid-binding resin (QIAEX II Gel extraction kit; Qiagen, Valencia, CA, USA).

Real-time PCR

Real-time PCR reactions for the individual, Pacific razor clam-specific housekeeping genes (18S and EF1-a) and genes of interest were run in separate wells (Tables 2 and 3). Briefly, 1 μl of cDNA was added to a mix containing 12.5 μl of Applied Biosystems Fast SYBR Green^® Master Mix (5 mM Mg 2+) (Qiagen, Valencia, CA, USA), 0.5 μl each of forward and reverse sequence specific primers (Invitrogen, Carlsbad, CA, USA) and 10.5 μl of RNase-free water; total reaction mixture was 25 μl. The reaction mixture cDNA samples for each gene of interest and reference genes were loaded into Fast 96 well plates in duplicate and sealed with optical sealing tape (Applied Biosystems, Foster City, CA, USA). Reaction mixtures that contained water but no cDNA were used as negative controls.

Statistical analysis

Both reference genes were evaluated for stability and ranked using the web-based analysis tool RefFinder (https://www.heartcure.com.au/for-researchers/) (Chen et al., 2015). Due to the differential stability of the reference genes, cycle threshold crossing values (C_T) for the genes of interest were normalized to the most stable reference gene (18S).

Analysis of qPCR data was conducted using normalized values (housekeeping gene threshold crossing subtracted from the gene of interest threshold crossing); the lower the normalized value, the more transcripts are present. A change in normalized value of 2 is approximately equivalent to a 4-fold change in the amount of the transcript.

For all data (gene transcription and biomarker assays), medians, 2.5% and 97.5% percentiles and ranges were calculated (NCSS, Statistical and Power Analysis Software, Kaysville, UT, USA). We used Gaussian linear random effects models to estimate site means for each gene transcription factor and biomarker parameter. We included sampling year as a random effect to account for variances introduced by collecting sample units over the course of 2 years. We fit separate models for each gene with maximum likelihood estimation using the lme4 package in R 3.5.0 (R Development Core Team, 2018), which accounts for the unbalanced data using Satterthwaite’s method. We conducted post-hoc Tukey tests of site level differences. To obtain site means for each year, we fit a Gaussian linear model with a site by year interaction. The resulting site means were used to test for correlations between genes and biomarker parameters. Relationships between gene transcript and biomarker data, as well as within gene transcript and biomarker data, were assessed in R using Pearson correlations.

We conducted two-dimensional non-parametric multidimensional scaling of the Bray–Curtis dissimilarity from gene transcripts (N = 112) and biomarkers (N = 73) using the Vegan package in R version 3.5.0. The graphical representations show individual razor clams clustered by similarity in transcription and biomarker values and not by pre-defined groups such as location. We obtained vectors describing the strength of each gene and biomarker contribution to the two non-metric multidimensional scaling (NMDS) axes for graphical display. We evaluated goodness of fit for NMDS models using stress plots.

Reference ranges and correlations

Biomarker and gene transcript medians, 2.5% and 97.5% percentiles, and ranges were created for use in future monitoring efforts (Table 4). Site-specific medians and ranges for the biomarker and gene transcription assays during each year are provided as Supplemental Material (Tables S1 and S2).

Statistically significant differences among sites were identified for five of the seven biomarkers (condition factor, shell thickness, hemocyte count, H₂O₂ and P450; Fig. 2), and four of the five genes (CaM, Ferr, HSP70, and PPIA; Fig. 3). One statistically significant positive correlation within the biomarker assay was present between condition factor and shell thickness (Fig. 4). Within the gene transcription assays, statistically significant positive correlations were identified between: CaM and Ferr, CaM and PPIA, and Ferr and PPIA (Fig. 5). Statistically significant negative correlations within the gene transcription assays were identified between: CaM and HSP70 and PPIA and HSP70 (Fig. 5). Interestingly, both populations are segregated by PPIA transcript levels. No negative correlations within the biomarkers were identified. Significant correlations (P ≤ 0.05) between the biomarker and gene transcription results also were identified (Table 5). For interpretation of results, we considered correlations with P values that were >0.30 (or <−0.30) to be of possible biological significance. For these calculations, the mean values for each site-year were used, as biomarkers and gene transcription were not measured on the same individual razor clams (and thus the sample size was lower, at n = 9). Regarding the gene transcription-biomarker correlations only (Table 5), as a larger normalized value is indicative of lower transcription levels, a negative R value actually indicates a positive relationship, and a positive R value indicates a negative relationship. A positive correlation was observed between Ferr transcription and hemocyte count and negative correlations between PPIA transcription and condition factor, as well as PPIA transcription and shell thickness.

NMDS results

Non-parametric analysis of biomarker data by site yielded results similar to the Pearson correlations analyses (2-dimensional stress = 0.222) (Fig. 6). Results show that NMDS 1 is primarily defined by condition factor, shell thickness, RNA:DNA, and hemocyte count. NMDS 2 is associated with HSP40, P450 and H₂O₂. The plot does not clearly show clusters, although there is some separation between two KATM sites (Kashvik and Swikshak) and the rest of the sites on NMDS2.

Non-parametric multivariate analysis of gene transcription data by site were similar to the Pearson correlation analyses (2-dimensional stress = 0.092) (Fig. 7). Our results show that NMDS 1 is heavily influenced by PPIA and HSP70. NMDS 2 is most strongly associated with HSP 90. The plot clearly shows two clusters almost entirely divided by park. Some individuals from KATM are found in the LACL cluster, but none of the individuals from LACL are found in the KATM cluster.

Nutrients

Higher quality and/or quantity of nutrients in the marine intertidal have been associated with increased condition and thicker shells in bivalves (Carmichael, Shriver & Valiela, 2004). The condition factor and shell thickness of Pacific razor clams were higher at LACL than KATM. On average, the condition factor of Pacific razor clams from LACL was 4.7 times higher and shells 1.6 times thicker than Pacific razor clams from KATM. In a concurrent study, LACL mussels had higher condition factors and thicker shells than KATM mussels (Counihan et al., 2019). The condition factor of Pacific razor clams at all KATM sites, except Swikshak in 2015, was below the range observed in Pacific razor clams held in a controlled environment with sufficient nutrients (K. Counihan, 2014, unpublished data). These results suggest that there were differences in nutrient availability and/or quality between parks.

Pathogen presence

Pathogen exposure generally results in stimulation of the immune system (i.e., increased hemocytes, hydrogen peroxide production, transcription of PPIA and Ferr) (Pipe, 1990; Galloway & Depledge, 2001; Song et al., 2009). The primary immune cells in bivalves are hemocytes and they release reactive oxygen species in response to foreign antigens (Pipe, 1992). Increased immune activity was detected at KATM in 2015, where hemocyte count, hydrogen peroxide production and PPIA and Ferr transcription all were elevated compared to LACL, suggesting an immune response was occurring at KATM. PPIA transcription was consistently elevated at KATM and Pacific razor clams also were in poorer condition compared to LACL; however, causal mechanisms underlying this relationship are not clear. The hemocyte count in Pacific razor clams from Polly Creek in LACL nearly doubled from 2015 to 2016, but there were no concurrent increases in hydrogen peroxide production or PPIA and Ferr transcription. A two-fold or higher increase in hemocyte count has been associated with immune stimulation (Coray, St. Jean & Bard, 2007; Duchemin et al., 2008). Exposure to contaminants has been shown to result in increased hemocyte production without the stimulation of a cell-mediated response (Galloway & Depledge, 2001). However, in other studies, decreased hemocyte production has resulted from exposure to contaminants (Akaishi et al., 2007).

Although we did not attempt to quantify pathogens or algal toxins in this study, both could potentially influence the status of the razor clam populations that we sampled. For example, nuclear inclusion X (NIX) is a disease of serious concern for Pacific razor clams and has decimated populations in Washington state in the past (Elston, 1986). However, NIX is thought to be rare or absent in Alaskan razor clam populations (Meyers & Burton, 2009) and likely was not a factor in this study. Harmful algal blooms, which are increasing at northern latitudes as ocean temperatures warm (Gobler et al., 2017), could also be an issue. Of particular concern are the toxic dinoflagellates in the genus Alexandrium, the causative agents of paralytic shellfish poisoning (Vandersea et al., 2018; Tobin et al., 2019). Elevated water temperatures were documented in the Gulf of Alaska in 2015–2016 (Danielson et al., 2019) as well as increased abundance of Alexandrium in lower Cook Inlet (Vandersea et al., 2018). Algal toxins are generally recognized as an additional stressor in warming waters (Griffith & Gobler, 2019) and could be a factor in the elevated immune activity detected in razor clams from KATM in 2015 (hemocyte count, H₂O₂, and PPIA and Ferr transcription were all elevated). However, further study is needed to ascertain Pacific razor clam responses to algal toxins, and to measure relative toxin exposures in clam populations sampled in different areas.

Contaminant exposure

Natural oil seeps have been reported throughout Cook Inlet and several are in the vicinity of our sampling sites (Becker & Manen, 1988). A study conducted in 2002 detected PAHs in sediments at multiple coastal sites in southcentral Alaska and the source was attributed to natural seeps (Saupe, Gendron & Dasher, 2005). The cytochrome P450 enzyme catalyzes reactions involving xenobiotics and endogenous compounds (Rewitz et al., 2006; Kingtong & Janvilisri, 2011). Cytochrome P450 activity and HSP40 levels were elevated at Swikshak in KATM and Silver Salmon and Polly Creek in LACL during 2016. HSP40 increases in response to stressors such as contaminants (Wang et al., 2012). The amount of oil seepage can vary over time and could result in variations in P450 activity among sites. However, as we did not measure PAHs at individual sites, we can only speculate about the relationship between P450 and HSP40 levels and PAH presence at Swikshak and other sites.

Metabolic activity

Metabolic activity is evaluated by quantifying the RNA:DNA ratio, which detects short-term fluctuations in RNA concentration over the course of 1–3 days prior to measurement (Buckley, Caldarone & Ong, 1999; Dahlhoff, 2004). RNA production was 1.5 times higher in KATM Pacific razor clams as compared to LACL Pacific razor clams and was higher at all sites in 2016 compared to 2015. Mussels collected from both parks in 2015 and 2016 also had elevated RNA levels at KATM and during 2016 (Counihan et al., 2019). These results indicate that metabolic activity was heightened at KATM and at both locations during 2016. Protein synthesis may have been induced by multiple factors such as food availability or a stressor (Dahlhoff, 2004). However, shell thickness and condition factor data suggest higher nutrients at LACL; thus increased protein synthesis at KATM is not likely caused by greater food availability. Additionally, all of the sites with higher RNA:DNA ratios in 2016 also had elevated P450 and/or HSP40 levels, suggesting the presence of stressor-stimulated protein production.

Correlations

Correlations within the biomarker assays, within the gene transcription panel, and among the biomarkers and genes were determined. In the biomarker assays, a positive correlation was found between condition factor and shell thickness. Condition factor indicates the nutritional status of razor clams, and it was not surprising that shell thickness would positively relate to overall condition, as razor clams with more nutrients might be expected to allocate more nutrients to shell formation. In other bivalve species, nutrient deficits have been shown to initiate shell metabolization (Masthanamma, Purushotham & Ramamurthi, 1984).

Given the multiple functions of the genes in our transcript panel, as well as the interconnectedness of genes in general, we expected and found numerous correlations among the genes (Fig. 5). For example, CaM, PPIA and Ferr can all be indicators of pathogen exposure. Although traditionally thought of as a response to thermal stress, HSP70 can also indicate general physiological stress. Thus, as expected, we found positive correlations between CaM and PPIA, CaM and Ferr and PPIA and Ferr. Negative correlations were found between HSP70 and PPIA, as well as HSP70 and Ferr. One potential explanation may be shifts in resource allocation away from one response towards another (Ernande et al., 2004).

As we describe significant differences in physiological and gene transcription results between parks and among sites, it is important to acknowledge the existence of synergistic or antagonistic effects of multiple stressors (Côté, Darling & Brown, 2016; Griffith & Gobler, 2019). Due to the numerous stressors present in any natural environment, it is not possible to correlate every causative influence with an effect on an organism’s physiology. In addition, physiological and genetic pathways are a web of often interconnecting paths; a single gene or physiological response may be initiated by multiple stressors. With multiple assay types (i.e., biomarker and gene transcription), our ability to understand the potential input of both large and fine scale processes improves. For example, our results point, in general, to a higher quantity or quality of nutrients at LACL in comparison with KATM, a higher pathogen response at KATM in comparison with LACL and a higher response to contaminants at Chinitna Bay and Swikshak than at other sites. However, it is also important to note that while we found differences in gene transcription and biomarker assay results between parks and among sites, most or all observations should still be in the range of “normal”.

An objective of our study was to acquire baseline biomarker and gene transcription reference ranges for razor clams. The establishment of baseline or reference ranges in any species is problematic. One option is to use data from captive clinically normal individuals as baseline or reference values. However, captivity comes with its own set of stressors; there is no population free of all external or internal stimuli. Longitudinal data, with the initial timepoint acting as baseline or reference, is an alternative approach. However, the initial measurements, which form the baseline, are simply that: a place from which to make relative comparisons, with the findings in one group compared with other groups. It is important to note in this case that varying levels of transcription within or among groups or individuals may be normal responses to stimuli and an indication of a properly functioning system.

Univariate correlations do not adequately describe the multiple interacting pathways represented by the gene transcripts and biomarkers described above. The vectors generated from our NMDS analysis can be used to evaluate which transcripts and biomarkers give similar information about the sampled population and which, if any, drive group separation in multidimensional space. For example, in our gene transcript analysis, CaM and Ferr seem to describe similar pathways in our populations, while PPIA and HSP90, perpendicular to each other, describe potentially different pathways. Additionally, vectors that are opposite each other on the same axis can be thought of as negative correlations. The size and position of the PPIA vector in relation to the group separation along the NMDS 1 implicate PPIA as a main driver of separation between parks (the differences in PPIA are also evident in the correlation plots, Fig. 5). Although the biomarker analysis failed to separate populations based on park, the vectors indicate that biomarker pairs: condition factor & shell thickness, HSP40 & P450 and RNA:DNA & hemocyte count may describe similar pathways. The gene transcript and biomarker ordinations provided different, but complementary information, suggesting that the combination of gene transcripts and biomarkers in a single ordination model should be a focus of future studies.

Another objective was the comparison of the two methodologies, gene transcription and biomarker assays. Several of the biomarkers and genes were associated with similar physiological functions and we anticipated correlations would arise. Although we did not expect complete agreement or duplication between the methods, we expected the methods to support one another. Ferritin transcription was positively correlated with hemocyte count. This relationship makes biological sense, as an increase in Ferr transcription could indicate a parasite pressure, which could be associated with elevated hemocyte count, an indication of immune stimulus. Negative correlations were found between PPIA and condition factor, as well as PPIA and shell thickness.

A limitation of this study was the low number of genes in the transcript panel. This stems from the fact that relatively little is known about the genetic makeup of invertebrates in comparison with vertebrates (Niu et al., 2013). GenBank^®, the National Institutes of Health genetic sequence database, an annotated collection of all publicly available DNA sequences, contains only three genetic sequences for the Pacific razor clam, all representing one gene (Clark et al., 2015). The implications of this lack of information are that each “new” target gene must be identified and sequenced using degenerate primers designed from closely-related species (if available) prior to amplification using standard PCR. PCR products are then sequenced and compared with known sequences in GenBank^®. Herein lies the problem: if GenBank^® has a dearth of invertebrate sequences, then comparing a sequence amplified from an equally unexplored invertebrate is likely to yield few matches, if any (Hou et al., 2011). Many of our sequences, when submitted to GenBank^®, yielded no matches. The results of our considerable efforts were five target genes and two reference genes in which we had complete confidence. Similar difficulties were found in transcriptomic studies of the Chinese razor clam, in which only 9.9% of the 147,669 transcribed sequences had significant matches in GenBank^® (Niu et al., 2013).

The limited number of genetic studies on razor clams is incongruent with their value as a commercial species and popularity as a sport fishery. Studies to date have focused on the Chinese razor clam, identifying genes responding to heavy metals (Wang et al., 2016), anthropogenic sound (Peng et al., 2016) and bacterial challenge (Niu et al., 2014, 2016; Peng et al., 2016), and have provided a foundation for the development of protocols for using razor clams as indicators of ecosystem health.