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

An emerging approach to ecosystem monitoring involves the use of physiological biomarker analyses in combination with gene transcription assays. For the first time, we employed these tools to evaluate the Pacific razor clam (Siliqua patula), which is important both economically and ecologically, as a bioindicator species in the northeast Pacific. Our objectives were to (1) develop biomarker and gene transcription assays with which to monitor the health of the Pacific razor clam, (2) acquire baseline biomarker and gene transcription reference ranges for razor clams, (3) assess the relationship between physiological and gene transcription assays and (4) determine if site-level differences were present. Pacific razor clams were collected in July 2015 and 2016 at three sites within each of two national parks in southcentral Alaska. In addition to determining reference ranges, we found differences in biomarker assay and gene transcription results between parks and sites which indicate variation in both large-scale and local environmental conditions. Our intent is to employ these methods to evaluate Pacific razor clams as a bioindicator of nearshore ecosystem health. Links between the results of the biomarker and gene transcription assays were observed that support the applicability of both assays in ecosystem monitoring. However, we recognize the need for controlled studies to examine the range of responses in physiology and gene transcripts to different stressors.


INTRODUCTION
Nearshore marine habitats worldwide are facing unprecedented challenges due to expanded commercial and industrial development, resource extraction and environmental 1950s to 1963, Alaska harvested the most Pacific razor clams of any state in the U.S. (National Marine Fisheries Service, 2016).
Pacific razor clams can be found in both intertidal and subtidal areas (to about 55 m depth), from the eastern Aleutian Islands, Alaska, to central California. Populations have been affected by overexploitation (Bishop & Powers, 2003;Szarzi et al., 2010), poor spawning (ADF&G, 2014) and earthquakes (Baxter, 1971). Historically, in Alaska, Pacific razor clams could be found from Southeast Alaska west across the Gulf of Alaska to the east end of the Aleutian chain (Kaiser & Konigsberg, 1977). Now they are found primarily in only a few areas, including the Alaska Peninsula, eastern Prince William Sound and the east and west coasts of Cook Inlet (ADF&G, 2010) and the status of subpopulations varies. For example, even within the relatively small area of Cook Inlet, the personal use fishery harvest of razor clams has been closed on the east side since 2015 due to low abundance estimates (ADF&G, 2018), while in contrast, the west side supports the only active commercial harvest of razor clams in Alaska (350,000-400,000 lbs/year) (ADF&G, 2019, personal communication) and a personal use fishery with no harvest limit (ADF&G, 2019).
Traditional bivalve health assessments used in environmental monitoring include biomarker assays that provide information on how an organism is functioning (Dahlhoff, 2004). A suite of biomarker assays can be used to evaluate various aspects of physiology to detect the effects of environmental and anthropogenic stressors on overall health of the bivalve. Changes in tissue and shell growth due to factors such as increased temperature, pollution, acidification or eutrophication can be detected by measuring the condition factor and shell thickness of the bivalve (Carmichael, Shriver & Valiela, 2004;Gagné et al., 2007;Gaylord et al., 2011). Hemocytes are the primary cell in the bivalve immune system and they phagocytose antigens and release reactive oxygen species as a killing mechanism (Pipe, 1990). Suppression of the immune response by contaminants can lead to disease susceptibility and pathogen exposure can stimulate an immune response (Akaishi et al., 2007;Gagné et al., 2002;Pipe & Coles, 1995). The immune status of bivalves can be evaluated by determining the number of circulating hemocytes and measuring hydrogen peroxide production. The RNA:DNA ratio assesses the metabolic condition of bivalves, because DNA levels remain constant in the cell while RNA varies based on intrinsic status and/or extrinsic perturbations. Therefore, short-term growth or stress can be identified by fluctuations in the RNA:DNA ratio (Caldarone et al., 2001;Lesser et al., 2010). The cytochrome P450 enzyme catalyzes reactions involving xenobiotics, and increased activity can reveal the presence of a pollutant (Rewitz et al., 2006). Changes in heat shock protein levels suggest exposure to a stressor, such as increased temperature or a contaminant (Cruz-Rodriguez & Chu, 2002).
Transcriptomics, the molecular investigation of changes in gene transcription, offers a new approach for assessing mechanisms by which stressors may affect organisms. The earliest observable signs of health impairment can be altered levels of gene transcripts, evident prior to clinical manifestation (Farr & Dunn, 1999;McLoughlin et al., 2006;Poynton & Vulpe, 2009). Gene transcription is the process by which information from the DNA template of a particular gene is transcribed into messenger RNA (mRNA) and eventually translated into a functional protein. The amount of mRNA transcribed from a particular gene is physiologically dictated by a number of intrinsic and extrinsic factors and analysis of mRNA can provide information about changes in the physiological state of an organism. Alterations (including increased or decreased amounts of mRNA) of transcribed genes indicative of physiological processes at the cellular level can identify exposure to stressors and elucidate mechanisms by which these stressors may have deleterious effects (Snell, Brogdon & Morgan, 2003;McLoughlin et al., 2006;Mos et al., 2008;Wu et al., 2008;Poynton & Vulpe, 2009;Miller et al., 2011;Evans & Hofmann, 2012).
Relatively little is known about the physiological or molecular responses of razor clams to environmental stimuli (Niu et al., 2013), as research into these processes has been minimal. Biomarker and gene-based diagnostics provide an opportunity to assess the health of individuals and populations of razor clams. However, a key requirement for the application of diagnostic technology in monitoring programs is the establishment of normal or "reference" ranges of values, to facilitate differentiation of natural variation among individuals from changes associated with altered environmental influences. Reference ranges can be determined via controlled laboratory experiments, where individual stressors are presented to organisms in systems devoid of multiple variables. Due to the absence of background noise, these studies should provide clearly defined reference ranges. However, interpretation of these ranges when applied to natural systems with multiple fluctuating variables can be problematic. Alternatively, reference ranges can be developed in natural systems, using location and point in time as references for comparisons with future monitoring efforts (i.e., longitudinal sampling; McLoughlin et al., 2006).
In this study, we hypothesized that biomarker and gene transcription assays on Pacific razor clams would provide complementary results, with no significant differences detected among sites considered to be relatively pristine. The sites were expected to provide a range of representative natural conditions to generate reference baseline data applicable to Pacific razor clams in this region.
Our specific objectives were to (1) develop biomarker and gene transcription assays with which to monitor the health of the Pacific razor clam, (2) acquire baseline biomarker and gene transcription reference ranges for Pacific razor clams, (3) assess the relationship between physiological and gene transcription assays and (4) determine if site-level differences were present.

MATERIALS AND METHODS
This study was conducted concurrently with a similar study on mussels, recently reported by Counihan et al. (2019). Study areas (Lake Clark National Park and Preserve, LACL and Katmai National Park and Preserve, KATM) were the same although the specific collection substrates within each park varied for mussels (protected mixed-sediment) vs. razor clams (exposed soft sediment). As analytical methods for the two studies were identical, data were collected as previously described in Counihan et al. (2019), but the specific genes targeted in the mRNA transcription varied for each 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 −1 horseradish peroxidase, type II) was added to each well and incubated for another 30 min in the dark. The reaction was stopped by

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 mg ml −1 , 0.2 mg mL −1 , 0.4 µg ml −1 , 0.8 µg ml −1 , 1.6 µg ml −1 , 3.2 µg ml −1 , 6.4 µg ml −1 , 10.0 µg ml −1 ). RNA from bovine pancreas (Sigma-Aldrich) was used for a RNA standard curve (0.4 µg ml −1 , 0.8 µg ml −1 , 1.6 µg ml −1 , 3.0 µg ml −1 , 6.0 µg ml −1 , 8.0 µg ml −1 , 12.0 µg ml −1 and 16.0 µg ml −1 ). 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 −1 ) 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 −1 ) 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 (

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 ml −1 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 mg of RNA template from each razor clam. Reaction conditions included 4 units reverse transcriptase (Omniscript, Qiagen, Valencia, CA, USA), 1 mM 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 Gao et al., 2007).
Peptidylprolyl isomerase A (PPIA) Proinflammatory, increased in response to pathogen stimulus  18S Reference (Niu et al., 2014;Moreira et al., 2012) Elongation Factor Alpha-1 (EF1a) Reference (Araya et al., 2008;Mateo et al., 2010) 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 ml of cDNA was added to a mix containing 12.5 ml of Applied Biosystems Fast SYBR Green Ò Master Mix (5 mM Mg 2+) (Qiagen, Valencia, CA, USA), 0.5 ml each of forward and reverse sequence specific primers (Invitrogen, Carlsbad, CA, USA) and 10.5 ml of RNase-free water; total reaction mixture was 25 ml. 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.

RESULTS
Over the 2 years, approximately 20 Pacific razor clams were collected from each site (Hallo Bay at KATM was not sampled in 2015 due to logistical constraints). A total of 116 razor clams were used for biomarker assays, of which 73 had a full complement of biomarker data (missing data due to logistical constraints of taking morphometric measurements in the field for razor clams from Chinitna Bay and Polly Creek in 2015) and 112 were used for gene transcription. The razor clams ranged in length from 28.5 to 158.3 mm and were generally larger at LACL (86.1-158.3 mm) than at KATM (28.2-103.7 mm). As size clearly varied, and presumably was an indication of age, we estimated age based on annuli counts. Age was tested (lme4 package, R 2.8.1; R Development Core Team, 2012) and no significant effects on razor clam biomarkers or gene transcription were noted; thus neither age nor size were considered in further analyses.

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 2 O 2 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 2 O 2 . 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.

DISCUSSION
This was the first study to assess razor clam health as an indicator of nearshore ecosystem health. In most nearshore ecosystem monitoring studies, the physiological responses of bivalves at a study site are compared with bivalves at a non-impacted reference site (Akaishi et al., 2007;Coray, St. Jean & Bard, 2007;Gagné et al., 2007;Halldórsson, Svavarsson & Granmo, 2005;Halldórsson et al., 2008;Shaw et al., 2011). However, when sampling occurred for this study, none of the sites were suspected of being impacted. Rather, our goal was to evaluate a set of biomarker and gene transcription assays which we developed (Objective 1) to assess change in the nearshore ecosystem, which included establishment of reference values for biomarker and gene transcription assays in Pacific razor clams (Objective 2) at sites in LACL and KATM and description of relationships among both assays (Objective 3). Additionally, we wanted to test our hypothesis that Pacific razor clams at all sites were physiologically similar (Objective 4). Although our primary purpose was to evaluate Pacific razor clams as a bioindicator species for long-term monitoring of nearshore marine areas, a further benefit is the acquisition of data that can be used to assess health of individuals and populations of Pacific razor clams, which may support managers who are grappling with fluctuating populations and closures of fisheries (ADF&G, 2010). Although differences in Pacific razor clams were identified among locations, the variances of physiological and gene transcription values within populations were relatively low when compared to data on bay mussels (Mytilus trossulus) collected from the same general locations in the two parks (Counihan et al., 2019). We suspect this may be due in part to lower variability in physical attributes (such as temperature and wave action) of razor clam habitat. Razor clams are infaunal bivalves and predominantly surrounded by interstitial sea water even at low tide. This would preclude them from many of the environmental fluctuations experienced by epifaunal bivalves such as mussels,  which have the ability to settle on a variety of intertidal surfaces from bedrock to gravel and are exposed to variations in air temperature, solar radiation and wave action (Helmuth et al., 2010).

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 2 O 2 , 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 andduring 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 , anthropogenic sound (Peng et al., 2016) and bacterial challenge (Niu et al., 2014(Niu et al., , 2016Peng et al., 2016), and have provided a foundation for the development of protocols for using razor clams as indicators of ecosystem health.

CONCLUSIONS
We hypothesized that biomarker and gene transcription assays focused on Pacific razor clams would provide complementary results, with no significant differences detected among sites considered to be relatively pristine. However, we found differences in physiological assay and gene transcription results between parks and among sites which indicate variation in both large-scale and local environmental conditions. Changing environmental conditions in coastal ecosystems necessitate methods to assess the health of intertidal communities and facilitate management decisions to sustain these resources. Gene transcription assays will be a valuable technique because they provide early evidence of changes in physiologic status. Biomarker assays are an established approach to assess bivalve health. Links between the results of the biomarker and gene transcription assays were observed that support the applicability of both assays in ecosystem monitoring. Further development of gene transcription methods as well as controlled laboratory exposure studies will improve our capacity to monitor for early signs of physiological impacts in razor clams and their ecosystems due to changing environmental conditions.