Differential response to stress in Ostrea lurida as measured by gene expression

Experimental design

Adult, hatchery produced oysters from three wild source populations (Dabob Bay, Fidalgo Bay, and Oyster Bay grown for 19 months at Clam Bay) in WA were used for this experiment. All oysters were held at 8 °C for two weeks at the University of Washington prior to the experiment. Oysters from each population (n = 8 per population) were subjected to acute temperature stress (submerged in 500 mL 38 °C sea water for 1 h), mechanical stress (120 g × 5 min; Sorvall T21, ST-H750 rotor) or served as controls (maintained at 8 °C). After the stress treatments, oysters were returned to 8 °C seawater and sampled at 1 h post stress (n = 72). Ctenidia tissue was resected from each individual and stored separately in 500 µL RNAzol RT (Molecular Research Center, Inc.), frozen on dry ice. All samples were stored at −80 °C for later analysis.

RNA isolation

RNA was isolated using RNAzol RT (Molecular Research Center, Inc., Cincinnati, OH, USA) according to the manufacturer’s protocol for total RNA isolation. Briefly, ctenidia tissue was homogenized in RNAzol RT, volume was brought up to 1 mL with RNAzol RT, vortexed vigorously for 15 s, and incubated at room temperature (RT) for 10 min. 400 µL of 0.1% DEPC-H₂O was added to the homogenized ctenidia tissue, vortexed for 15 s, and incubated at RT for 15 min. The samples were centrifuged for 15 min, 16,000 g, at RT. After centrifugation, 750 µL of the supernatant was transferred to a clean tube, an equal volume of isopropanol added, vortexed for 10 s, and incubated at RT for 15 min. The samples were centrifuged at 12,000 g for 10 min at RT. The supernatant was discarded and the pellets were washed with 500 µL of 75% ethanol (made with 0.1% DEPC-H₂O) and centrifuged at 4,000 g for 3 min at room temperature. This wash step was then repeated. Ethanol was removed and pellets were resuspended in 100 µL of 0.1% DEPC-H₂O. Samples were quantified using a NanoDrop1000 (ThermoFisher, Waltham, MA, USA) and stored at −80 °C.

DNase treatment and reverse transcription

Total RNA was treated with DNase to remove residual genomic DNA (gDNA) using the Turbo DNA-free Kit (Ambion/Life Technologies, Carlsbad, CA, USA). The manufacturer’s rigorous protocol was followed. Briefly, 1.5 µg of total RNA was treated in 0.5 mL tubes in a reaction volume of 50 µL. The samples were incubated with 1 µL of DNase for 30 min at 37 °C. An additional 1 µL of DNase was added to each sample and incubated at 37 °C for an additional 30 min. The DNase was inactivated with 0.2 volumes of the inactivation reagent according to the manufacturer’s protocol. Samples were quantified using a NanoDrop1000 (ThermoFisher, Waltham, MA, USA). Treated RNA was verified to be free of gDNA via qPCR using actin primers (see Primer Design section below) known to amplify gDNA.

Reverse transcription was performed using M-MLV Reverse Transcriptase (Promega) with oligo dT primers (Promega, Madison, WI, USA), using 250 ng of DNased RNA. The RNA was combined with primers (0.25 µg) in a volume of 74.75 uL, incubated at 70 °C for 5 min in a thermal cycler without a heated lid (PTC-200; MJ Research), and immediately placed on ice. A master mix of 5× Reverse Transcriptase Buffer (1× final concentration; Promega), 10 mM each of dNTPs (0.5 mM final concentration of each dNTP; Promega, Madison, WI, USA), and M-MLV Reverse Transcriptase (50 U/reaction) was made and 25.25 µL of the mix was added to each sample (final reaction volume 100 µL). Samples were incubated at 42 °C for 1 h, followed by 95 °C for 3 min in a thermal cycler without a heated lid (PTC-200; MJ Research, Waltham, MA, USA), and then stored at −20 °C.

Quantitative PCR

Quantitative PCR

Quantitative PCR reactions were carried out using Ssofast Evagreen Supermix (BioRad, Hercules, CA, USA). Forward and reverse primers (Integrated DNA Technologies, Coralville, IA, USA) were used at a final concentration of 0.25 µM each. Sample cDNA was diluted (1:9) with molecular-grade water. Nine microliters of diluted cDNA was used as template. Reaction volumes were 20 µL and were run in low-profile, non-skirted, white qPCR plates (USA Scientific) with optically clear lids (USA Scientific, Ocala, FL, USA) in a BioRad CFX Real Time Thermocycler (BioRad, Hercules, CA, USA) and DNA Engine Opticon 2 System (BioRad, Hercules, CA, USA). Cycling conditions were: one cycle of 95 °C for 10 min; 40 cycles of 95 °C for 30 s, 60 °C for 1 min, 72 °C for 30 s. Two qPCR replicates were run for each sample, for each primer set.

Statistical analysis

To calculate relative expression levels for each gene, cycle quantity (Cq) or cycle threshold (Ct) values were calculated using BioRad CFX Manager 3.1 (version 3.1.1517.0823, Windows 8.1) and Opticon Manager 3 (Windows 8.1), respectively. This was accomplished by subtracting global minimum fluorescence from samples and determining the point in the cycle which amplification reached exponential amplification phase. Default settings were accepted for each program to ensure reproducibility. The BioRad CFX Manager used default settings of single threshold for Cq determination and baseline subtracted curved fit for each run. The Opticon Manager used default settings of subtract baseline via global minimum, which estimated the threshold as being between 0.019 and 0.028. Gene expression values were determined as normalized mRNA levels using the following equation (ΔCt): 2^−ΔCt; where ΔCt is: (target Ct—actin Ct) (Schmittgen & Livak, 2008). Actin expression levels were determined to be consistent across all samples and served as an internal amplification control to use for expression normalization. Data from ΔCt did not exhibit normal distributions, so were log transformed (log ΔCt), to establish normal data distributions for statistical analysis. Two-way analysis of variance (ANOVA) followed by Tukey’s Honestly Significant Difference post hoc test (base, R Core Team, 2013) were performed on log ΔCt for each target (p < 0.05).

Gene expression analysis

Without considering separate populations, acute heat shock resulted in statistically significant increases in expression of coactivator-associated arginine methyltransferase 1 (CARM1) (n = 24 oysters per treatment, ANOVA, df = 2, Tukey’s HSD p = 0.00007) (Fig. 1) and Histone 2AV (H2AV) (n = 24 oysters per treatment, ANOVA, df = 2, Tukey’s HSD p = 0.001) (Fig. 2). A statistically significant increase in expression of tumor necrosis factor receptor-associated factor 3 (TRAF3) (n = 24 oysters per treatment, ANOVA, df = 2, Tukey’s HSD p = 0.008) (Fig. 3) occurred upon exposure to mechanical stress.

There was a clear difference in response to mechanical stress in oysters from Oyster Bay as compared to oysters from Dabob and Fidalgo Bays. Specifically, upon heat shock, H2AV expression in oysters from Oyster Bay increased (n = 8 oysters per population, ANOVA, df = 4, Tukey’s HSD = 0.05) (Fig. 2) when compared to the control. When exposed to mechanical stress, bone morphogenic protein 2 (BMP2) (n = 8 oysters per population, ANOVA, df = 4, Tukey’s HSD p = 0.03) (Fig. 4) and growth-factor receptor bound protein 2 (GRB2) (n = 8 oysters per population, ANOVA, df = 4, Tukey’s HSD p = 0.03) (Fig. 5) expression was decreased in the Oyster Bay population, whereas there was no significant differences in responses in the other populations. Additionally, significant interactions were identified between population and treatment in both BMP2 and GRB2 (p < 0.05).

There was no statistical difference in expression in Peptidoglycan recognition protein 1 (PGRP), toll-like receptor 2 type 1 (TLR), and prostaglandin E2 receptor EP4 subtype (PGEEP4) (Figs. 6–8, respectively) within any comparison. Heat shock protein 70 gene expression was significantly different between temperature and mechanical stress (n = 24 oysters per treatment, ANOVA, df = 4, Tukey’s HSD p = 0.006) (Fig. 9).

Response to temperature stress

The response of Ostrea lurida to acute heat stress appears to include an alteration in gene regulatory activity and the innate immune response, as indicated by significant increases of H2AV (Fig. 2) and CARM1 (Fig. 2) gene expression one hour post-temperature stress.

Histone 2AV, H2AV, is a variant of the histone H2A protein. This variant has been shown to act as a transcription promoter agent as well as assist with heterochromatin formation. Truebano et al. (2010) characterized changes in transcription in Antarcticclams, Laternula elliptica, and found that an H2A variant was significantly upregulated under heat stress conditions (3 °C for 12 h). In addition to involvement in the heat stress response, histone H2A has been shown to exhibit antimicrobial properties in three invertebrates: two marine invertebrates (Pacific white shrimp and scallops; Patat et al., 2004; Li et al., 2007a), as well as in a freshwater shrimp (Arockiaraj et al., 2013). In D. melanogaster, H2Av is phosphorylated in response to DNA damage (Madigan, Chotkowski & Glaser, 2002) to inhibit apoptosis, suggesting an additional role in in cellular survival.

Coactivator-associated arginine methyltransferase 1, CARM1, is involved in transcriptional activation via methylation of histones (Chen et al., 1999; Lee et al., 2005). This in turn affects the ability of transcription factors to bind and transcription to proceed. It is possible that increases in CARM1 expression could indicate that overall gene regulatory activity is increased in response to temperature stress. Our results are similar to those of Wang et al. (2011) where researchers described an increase in expression of Histone-arginine methyltransferase in the sea cucumber, Apostichus japonicus, after experiencing 25 °C temperatures for seven days. The authors suggested that this was due to an induced dormancy and lower metabolic rate to provide resources for stress resilience. CARM1 is also a component of the cellular immune response, as it has been identified as a regulator of NF-kB (Covic et al., 2005). Thus another explanation is that acute heat could possible impact the immune response, likely in a negative manner. Future work, that would be relevant to restoration activities, should increase the number of stressors examined in oyster to include pathogens.

Increases in HSPs are often observed in response to stress, but this study only found a significant difference of mRNA expression of HSP70 in the Oyster Bay population between mechanical and heat stresses (Fig. 9). Brown et al. (2004) found the maximum HSP expression in O. lurida occurred 24–48 h post exposure to 39 °C. The absence of a strong response of HSP70, relative to the control group, could be related to temporal changes in expression or an isoform-specific response, as there are many genes in this gene family, particularly in oysters (Clegg et al., 1998; Piano et al., 2005). Mediterranean mussels, Mytilus galloprovincialis, have shown different isoforms of heat shock proteins and cognates that have differential expression patterns caused by heat, mercury exposure, and chromium exposures stressors suggesting that the isoforms have slightly different functions (Franzellitti & Fabbri, 2005). Additionally, there are members of the HSP70 gene family that are constitutively expressed and do not exhibit increases in mRNA in response to heat stress (Sorger & Pelham, 1987; Somji et al., 1999). Without a sequenced genome for Ostrea lurida, combined with utilizing an incomplete transcriptome, it is difficult to ascertain how many isoforms might exist, as well as the number of alternatively spliced products. Upon addition of new genomic resources the entire family of molecular chaperones could be examined and compared across populations.

Response to mechanical stress

Mechanical stress increased expression of inflammation-related target genes. In all populations, there was a significant increase in immune system-related responses seen via the expression of tumor necrosis factor receptor-associated factor 3, TRAF3 (Fig. 3), which is involved in internal tissue damage recognition and apoptosis. The main function of TRAF3 is to assist in cell death initiation caused by stress conditions within tissues (Arch, Gedrich & Thompson, 1998). Upregulation in relation to mechanical stress could be akin to inflammation occurring due to edema from the mechanical stress and used to remove damaged cells as suggested by Roberts, Sunila & Wikfors (2012) when C. virginica were exposed to mechanical stress. Significant differences in expression of other immune system targets such as PGRP, TLR, and PGEEP4 were not detected (Figs. 6–8, respectively), but other studies have found that the time scale for expression may vary (Meistertzheim et al., 2007; Farcy et al., 2009).

Population differences

We suspected that the Dabob Bay population would have demonstrated a more pronounced response to stress as this population is subjected to greater environmental fluctuations with respect to salinity and temperature (Heare et al., 2017). Contrary to our hypothesis, oysters from Oyster Bay were the only population that exhibited a difference in gene expression in response to mechanical or heat stress. Oysters from Oyster Bay parents showed an increase in H2AV expression during heat stress as compared to control (Fig. 2), a decrease in BMP2 and GRB2 upon mechanical stress (Figs. 4 and 5, respectively), and differences in HSP70 expression between heat and mechanical stresses (Fig. 9). Given the putative function of H2AV in transcriptional regulation (Table 1), the increase in expression could be indicative of the role of this protein in controlling the molecular response to stress. Bone morphogenic protein 2, BMP2, and growth-factor receptor bound protein 2, GRB2, were significantly decreased in expression which could be indicative of growth inhibition. Both genes are related to growth and development of tissues, with BMP2 being a pre-cursor to osteoblastic cells that produce shell (Pereira Mouriès et al., 2002) and GRB2 is used for signal transduction between cells during growth phases (Oda et al., 2005). By down-regulating these targets, this may be an effort to reduce energetically costly processes in favor of processes that promote survival during stress events. Organisms faced with stress are often required to reallocate energy resources to homeostasis-related functions in an effort to improve long-term survival of the species (Sokolova et al., 2012). This change in expression coupled with the up-regulation of H2AV (Fig. 2) is in accord with the idea of shifting priorities for stress resilience.

Interactions were identified between population and treatment for both BMP2 and GRB2. Differences between gene expression in control and mechanical stress in the Oyster Bay population are driving this interaction for both genes. Although statistical interactions of this nature are difficult to interpret, it could be related to fact the Oyster Bay population is from a relatively “low-stress” environment (i.e., abundant food and less-pronounced temperature fluctuations).