Nominal and mean (standard deviation) measured concentrations of SO42- and As in treatments (
Freshwater mussels of the Clinch and Powell rivers of Virginia in the southeastern United States have been heavily impacted by runoff, leachates, or spills of materials related to coal extraction, processing, and use. Assays quantifying sublethal impacts of such wastes are needed. We assessed gene transcriptional markers in a laboratory study under controlled conditions, focusing upon arsenic (arsenate, As(V)) and sulphate, contaminants related to coal mining and processing. Pheasantshells Actinonaias pectorosa collected from the Clinch River were subjected to a 28-day chronic exposure to control or environmentally relevant concentrations of each compound. We compared gene expression in digestive gland among parasite-free, female pheasantshells among control and contaminant-exposed individuals using the Illumina HiSeq platform. Statistically significant differential expression of particular genes was observed among control mussels and those exposed to either arsenate or sulfate. Chemical stress was as likely to cause under-expression as it was to cause over-expression of particular genes. Arsenate and sulfate induced up- or down-expression of different suites of 50-100 genes. Our results provide proof-of-principle for using RNAseq technology to approach issues of toxicogenomics in freshwater mussels. The candidate markers could be validated for quantitative PCR assays for rapidly assessing single-gene responses to exposure to toxic compounds.
- differential gene expression
- transcriptional markers
1.1. Response of aquatic organisms to exposure to toxic substances
Exposure to toxic compounds has not only lethal but also important sublethal effects upon affected individuals. Successful identification of molecular mechanisms underlying response to toxic exposure depends upon development and use of a suitable set of assays, and several approaches are potentially available. While various biomarkers [1, 2] or gene expression [3, 4] assays have been demonstrated for marine mollusks in recent years, there has been no parallel work for freshwater mussels, many of which are of conservation concern. More fundamentally, adaptation of existing assays to freshwater mussels is not entirely attractive because a focus on candidate genes of known function certainly will miss genes or biochemical pathways that are not known. A preferable approach would be to screen global gene expression and thereby identify genes of interest not only in known but also in unknown pathways. Three approaches are available. Microarrays are available for several marine mollusks, and have in some cases been used for purposes of characterizing responses to environmental stressors including toxins or other bioactive compounds (e.g., marine mussels [5–8], Manila clam
1.2. Freshwater mussels
Freshwater mussels (Class Bivalvia: Family Unionidae) have their center of biodiversity in the southeastern United States. However, many regionally important species face a variety of threats, including exposure to toxic compounds. In particular, mussels of the Clinch and Powell River systems of southwest Virginia have been heavily impacted by runoff, leachates, or spills of materials related to coal extraction, processing, and use for electric power generation. Given the continued operation of coal extraction and processing facilities, the shift from deep‐ to surface‐mining practices, and the increase in coal‐bed gas extraction wells in the Clinch and Powell River system, defensible biological assays for assessing the impacts on key components of this aquatic ecosystem must be developed. These assays could provide critical information for assessing the impacts of future toxic events and thereby allow characterization of potential long‐term effects of resource extraction on freshwater mussel populations. This project was conducted to show the response of freshwater mussels to selected physiological stressors. Within this context, the objective of the study was to screen gene‐transcriptional markers in a laboratory study under controlled conditions, focusing upon arsenate and sulfate, two contaminants related to coal mining and processing. Identification of such genetic markers could prove essential for use in future nonlethal determinations of contaminant impacts to federally listed mussels.
1.3. Contaminant selection
The coal industry is a potential source of numerous contaminants to aquatic environments. In several Virginia watersheds, habitat for imperiled freshwater mussel populations exists in close proximity to coal industry‐associated activities. These activities include mining and coal‐fired electric power plants with on‐site storage of coal combustion residue (CCR). Two contaminants were selected for this study based on their association with surface mining (sulfate, SO42-) or CCR storage (arsenate, As(V)). The contaminants As and SO42- were selected due to their environmental relevance and the potential for concentrations to become elevated due to activities related to coal mining and power generation.
Freshwater bivalves appear to be relatively tolerant to acute effects of As; the estimated 96‐h LC50 for Asian clam
Freshwater bivalves also appear to be relatively tolerant to acute effects of SO42-. Soucek and Kennedy  determined a sulfate 96‐h LC50 for grooved fingernail clam
For each contaminant, a relatively ‘high’ concentration was selected based on worst‐case conditions previously measured in heavily impacted environments. The concentration of SO42- (1250 mg/L) was based on low‐flow measurements up to 1200 mg/L in a notoriously polluted river in West Virginia, USA . The high concentration of As (1000 µg/L) was based on measurements in pore‐water downstream of a massive CCR spill (up to 1200 µg/L; ). For both chemicals, the selected concentrations were not expected to cause mortality based on results of acute and chronic toxicity studies conducted with other bivalves (
2.1. Mussel collection and holding
Pheasantshell mussels (
Sodium arsenate heptahydrate (Na2HAsO4*7H2O, American Chemical Society (ACS)‐certified reagent grade) was purchased from J.T. Baker Co. (Avantor Performance Materials, Center Valley, PA). A 100‐mM (7.492 g/L) As (V) stock solution in ultrapure water was prepared at the beginning of the study. ACS‐certified reagent‐grade sodium sulfate (Na2SO4) was purchased from Fisher Scientific (Fair Lawn, NJ). Stock solutions were prepared weekly in ultrapure water.
2.3. Mussel exposure system
Mussels were exposed in 18‐L downweller‐bucket systems  modified to accommodate adult mussels . Buckets were placed into a 757‐L container filled with well water to serve as a temperature control bath. The target exposure temperature of 23°C was maintained in the water baths using aquarium heaters. Each bath held up to five buckets. In the bucket systems, unfiltered water from the FMCC pond was used for acclimation, control water, and make‐up water for each treatment. Mussels were fed daily with a 1:1:1 algal cell ratio from three premixed commercial microalgae diets (Nanno 3600, Shellfish Diet 1800, and TP 1800; Reed Mariculture, Campbell, CA). A concentrated stock solution was prepared in deionized water and 1 mL of stock was added to each bucket to achieve a final concentration of 24,000 cells/mL.
2.4. A 7‐day acute exposure
Thirty mussels were allocated to 10 buckets (three mussels each, arbitrarily chosen) on October 7, 2014 for a preliminary 7‐day test of acute toxicity. This test was conducted to determine the potential for mussel survival in solutions of As(V) and SO42- over 28 days; treatments were not replicated. Five buckets (
2.5. A 28‐day chronic exposure
Sixty‐four mussels were randomly selected on October 16, 2014. Four mussels were allocated to each of 16 buckets, and the length of each mussel was recorded using dial calipers. Mussels were allowed to acclimate in the buckets for 11 days prior to the start of the chronic exposure. During the acclimation period, dissolved oxygen and ammonia levels were monitored using the methods described below to ensure that the feeding rate was appropriate to maintain acceptable water quality.
On Day 0 of the exposure (October 27, 2014), each bucket received a 100% water exchange and was randomly assigned to one of four treatments/controls (
For each bucket, mussel mortality was checked daily and recorded, and a 100% water exchange was conducted weekly. Prior to water exchanges, temperature (°C), specific conductance (µS/cm), dissolved oxygen (DO; mg/L and % saturation), and pH were measured using a calibrated YSI 556 Multi‐Probe Sensor (YSI Inc., Yellow Springs, OH), and water samples were collected. Total ammonia‐nitrogen (TAN) was measured in unfiltered samples using a Hach DR/2400 meter following the manufacturer's methods. Concentrations of elemental As and S were measured in filtered (0.45‐µm) water samples. Elemental quantification was performed by the Virginia Tech Soil Testing Laboratory using an inductively coupled plasma atomic emission spectrometer (ICP‐AES, Spectro Analytical Instrumentation, Kleve, Germany) following the laboratory's standard operating procedures and quality assurance/quality control methods . Sulfate concentrations were calculated from measured total S; all S were assumed to be present as SO42- due to the buckets’ oxygenated environment. Although only total As was measured, all As are assumed to have been present as As (V) due to the oxygenated environment in the buckets. Alkalinity and hardness of the pond water were measured weekly as part of a concurrent study . The 28‐day exposure concluded on November 24, 2014.
2.6. Mussel dissection
At the end of the study on Day 28, all mussels were removed from the buckets. Each mussel was measured using dial calipers and opened using reverse pliers. Each mussel was evaluated for gravidity, that is, the presence of glochidia in the outer marsupial gills. Mussels were removed from their shells by severing the adductor muscles. The mantle tissue was removed and samples of gill, kidney, and gonad were excised and fixed in neutral buffered formalin for histological evaluation. The digestive gland was removed from the remainder of the tissue and divided into three portions, with one portion each preserved separately for histological evaluation and RNA extraction. The remaining portion was further divided into thirds for separate biochemical analyses. Samples for biochemical work were immediately frozen. Samples for screening of gene expression were preserved in RNALater (Qiagen). All dissection tools were sterilized in a concentrated bleach solution between dissections of each mussel.
2.7. Transcriptome analysis
Only parasite‐free, female mussels (determined via histological analysis) were used for transcriptomic analysis. RNA isolated from digestive gland tissue samples was transcribed in vitro to double‐stranded cDNA, yielding a library of transcripts for each individual mussel. Individual‐specific adapters were ligated onto transcripts for each individual so that we could subsequently identify RNA sourced from the respective individuals. Barcoded cDNA samples—three from control and six from contaminant‐exposed individuals (three per treatment)—were sequenced using the Illumina HiSeq platform. Adapters were removed from the raw sequencing reads. Duplicated and low‐quality reads were discarded using FastqMcf  with default parameters. To exclude possible contamination, all reads were aligned to a bacterial database downloaded from NCBI (http://www.ncbi.nlm.nih.gov/), and only unmapped reads were used for assembly. Processed reads from all nine samples were merged together, assembled with Trinity  with parameter—trimmomatic, after duplications were removed. TransDecoder (Trinity package) was used to identify candidate‐coding regions within assembled transcripts, and transcripts with open‐reading frame (ORF) lengths less than 300 (100 amino acids) were removed from the assembly. The final transcript assembly was used as a reference for gene annotation and expression calculation. Transcripts were mapped to the nonredundant protein database [(NR database) from NCBI] using BLAST (v. 2.2.28). Alignments with threshold
2.8. Histological analysis
Sections from gonads, digestive glands, and gills were stained with hematoxylin and eosin for routine histological evaluations, and sections from kidneys were stained with Long Ziehl‐Neelsen stain for elaboration of lipofuscin . Histological evaluations of tissues were conducted by light microscopy using point counting . Evaluations determined fractions of reproductive acini containing mature and/or developing gametes, acini containing atretic and/or resorbing gametes, digestive gland diverticula cells with abnormal cytoplasm, gill filament termini without cilia, and kidney diverticula cells containing lipofuscin inclusions. Genders (female, male, hermaphrodite, and indeterminate) of mussels and incidences of parasitic infestation were recorded. Histological data were analyzed using generalized linear‐mixed models for binomial data with SAS GLIMMIX.
2.9. Biochemical analysis
Digestive gland samples for biochemical analysis were placed in 1.5‐mL tubes, immediately frozen, and held at -80°C. Analyses of biochemical parameters in mussel digestive glands were conducted in duplicate and expressed per mg of protein determined using a bicinchoninic acid method (Pierce BCA Protein Assay Kit, Thermo Scientific, Rockford, IL). Activities of catalase (CAT), glutathione S‐transferase (GST), glutathione peroxidase (GPx), glutathione reductase (GR), malondialdehyde (MDA), total‐ and Na+, K+−ATPase, and lactate dehydrogenase (LDH), and the concentration of reduced glutathione (GSH) were determined.
2.10. Correlations with histological and biochemical markers
Correlations between biochemical measurements, histological variables, and gene expression were evaluated using rank‐transformed data (Spearman's rho) to minimize effects of the different data types and scales for each class of variable. Due to the large number of variables examined, the criterion for statistical significance was reduced (α = 0.01) and plots of all significant relationships were examined to ensure linearity. To select genes for assessment of correlations, genes with the lowest false‐detection rate (<9.67 × 10-5) and
3.1. Water chemistry
Mean‐measured concentrations of SO42- and As were similar to nominal concentrations (Table 1). Background SO42- (∼17 mg/L) was present in the pond water; thus, measured concentrations were slightly above nominal concentrations. There was no measureable background As in the pond water.
Water quality measurements were within acceptable ranges for toxicity tests with freshwater mussels. During the study, the pond had a mean measured hardness of 204 mg/L CaCO3 and a mean measured alkalinity of 196 mg/L CaCO3. Mean measured temperatures were similar between treatments, and temperature was generally within 1°C of the 23°C target temperature (Table 2). Measured pH was similar across all treatments and varied little over the course of the study. Specific conductivity was similar in the controls and As treatments, and was elevated in the SO42- treatments as expected. Dissolved oxygen concentrations remained relatively high for the duration of the study (lowest recorded value: 75% saturation, 6.20 mg/L, on Day 28). Dissolved oxygen concentrations were consistent between treatments but did show variation during the study, with a decrease in overall mean saturation from 86% on Day 7 (range 79.1–88.7%) to 78% on Day 28 (range 75–78.8%). TAN concentrations were also consistent between treatments, but varied over the course of the study. The overall mean TAN concentration was 0.001 mg/L on Day 7 (range 0–0.01 mg/L) and increased to 0.068 mg/L on Day 21 (range 0.05–0.08 mg/L). The maximum measured TAN concentration was 0.080 mg/L in several control and SO42- units, but this value is far below the 2013 acute and chronic water quality standards for ammonia (0.96 and 0.24 mg/L, respectively, at pH 8.6 and 23°C), derived based on the sensitivity of mussels .
|Treatment||Temp. (°C)||pH||Specific conductivity (mS/cm)|
|Arsenic||22.8 (0.9)||8.58 (0.03)||0.436 (0.004)|
|Sulfate||22.6 (0.8)||8.64 (0.05)||2.89 (0.10)|
|Control_1||22.6 (1.1)||8.51 (0.27)||0.431 (0.005)|
|Control_2||22.8 (1.2)||8.61 (0.07)||0.434 (0.004)|
No mussels died during the course of the experiment. The mean length of all exposed mussels was 85.1 mm (range 79.0–99.3 mm). The mean lengths were not significantly different between treatments (GLIMMIX,
3.3. Reference assembly
The analytic program Trinity handled 61,774 transcripts, assembling them into the putative transcripts of 34,019 genes. The average length of these contiguous sequences, or “contigs,” was 1012 base pairs.
3.4. Differential expression analysis
Statistically significant (
Comparing gene expression among control and arsenate‐treated mussels (Figure 1), 18,572 transcripts were not affected by treatment, 28 were underexpressed in controls relative to treated mussels, and 30 were overexpressed.
Comparing gene expression among control and sulfate‐treated mussels (Figure 2), 18,469 transcripts were not affected by treatment, 74 were underexpressed in controls relative to treated mussels, and 86 were overexpressed.
Comparing gene expression among pollutants, 18,444 transcripts were not differentially expressed, 128 were more highly expressed in sulfate‐exposed mussels, and 57 more highly in arsenate‐exposed mussels (Figure 3). Overall, the sulfate treatment was a ∼3× greater stress factor in terms of the number of differentially expressed genes than was the arsenate (HA) treatment. We offer the explanation that in aquatic species, maintaining homeostasis in the face of ionic and osmotic stressors affects many cellular processes, while heavy metal toxicity affects relatively few. Comparing the effects of HS and HA directly, it was evident that HS causes >2 times more overexpression, again potentially pointing to a greater effect of HS. Interestingly, and as would be expected, the expression of different genes was affected. As discussed below, these genes are candidate markers for indicating exposure of mussels to the respective pollutants.
3.5. Annotation of expressed genes
The proportion of genes that we were able to annotate by reference to GenBank was relatively low (63%), as mollusks are nonmodel organisms and hugely underrepresented in genomic databases. We note that the proportions of genes that we annotated are very similar to those reported by Bai et al. [33, 34] for the freshwater pearl mussel
Among annotated underexpressed genes in the arsenate treatment relative to controls were guanylate binding protein 1 (
Among annotated underexpressed genes in the sulfate treatment relative to controls were serine protease inhibitor dipetalogastin precursor and zinc‐binding Sp‐Hypp_8991. Interestingly, the expression of serine protease inhibitor seems to offer stress tolerance via delayed senescence . Among annotated overexpressed genes in sulfate relative to control were NADPH‐dependent alpha‐keto amide reductase (YDL124W), and poly [ADP‐ribose] polymerase 15‐like (LOC101731886). Multiple aldo‐keto reductases fulfill functionally redundant roles in stress response in yeast .
3.6. Correlations with histological, biochemical, and gene expression markers
Significant correlations (
There were strong correlations (
Gene expression profiling via next‐generation sequencing has emerged as a tool to assess the biological impact of environmental pollutants and natural chemical stressors. This approach supports better understanding of whole‐genome expression responses to chemical stress, aids in the identification of ecological and toxicological modes of action, and provides hundreds (or in some cases thousands) of rigorously tested markers for stress responses. Gene expression profiles thus provide a more comprehensive, sensitive, and actionable insight into toxicity than typical toxicological parameters such as morphological changes, altered reproductive capacity, or mortality . For example, it has been demonstrated that chemicals from the same class of compounds give rise to discernible gene expression profiles that bear more similarity to each other than to patterns corresponding to exposure to compounds from a different class . It would be useful for a database for gene expression responses to environmental pollutants to be developed for recognizing compound‐specific expression profiles and molecular signatures of stress responses.
Functional annotations of genes may provide insights into affected regulatory and proliferative and repair/adaptive pathways, with immediate management and conservation implications. To improve the rate of functional annotations, expanded sequencing with more complete sequence information will be required. Recent bioinformatics developments, such as improvements to Trinotate (https://trinotate.github.io/), also hold promise of improved functional annotations. However, we note that complete annotations of affected genes are not necessary. Unannotated transcripts are potentially useful markers not only for compound‐specific expression profiling but also candidates for cost‐effective and rapid individual biomarker assays based on real‐time polymerase chain reaction (PCR), which could be developed for use in the laboratory or even in field settings. Despite the limited sampling, we found dozens of highly significant gene expression markers (after the correction for false discovery rate) for the two stressors tested. As many as 66 genes were expressed only in the presence of the two compounds, with zero transcripts (i.e., no sequences found) in the control group.
We sought quantitative differences between the RNA pools of control and arsenate‐ or sulfate‐exposed mussels. We sought evidence of altered expression of genes known to be linked to toxin response (e.g., heat‐shock proteins, glutathione transferase, metallothioneins, cytochrome P450, etc.), and confirmed that some of the pathways, such as poly(ADP‐ribose), were indeed induced. We also sought to identify genes whose expression is influenced by toxin exposure, but of which we have no existing knowledge. We found many such putative markers of toxin‐induced responses. In this latter case, more complete sequence data will be required to enhance functional annotation, but even if we do not know the functions of all the associated genes, they are still candidate markers for pollutant response.
4.1. Comparisons among biomarker types
Exposure to toxic compounds would be expected to affect gene expression, molecular flux through biochemical pathways, and ultimately histological structures in freshwater mussels. Further, the respective markers would be expected to be causally related to one another. Hence, we sought statistical indicators of such correlations among our data sets.
Correlations among histological and genetics data provide provisional bases for hypotheses regarding links among histological lesions and genes identified during this study. While many of the significant correlations involved genes without known function, a few involved genes of known function. Gene LOC101731886 was significantly negatively correlated with fractions of reproductive acini containing mature or developing gametes. The gene has been linked to the activity of poly [ADP‐ribose] polymerase (15‐like) (
Abundances of lesions in digestive cells were negatively correlated with genes BRAFLDRAFT_118372 (histamine N‐methyltransferase activity) and CGI_10002926 (uncharacterized protein with unknown function) [47–49]. A negative correlation between gene regulation of histamine methylation and lesions in digestive cells seems reasonable, because with increasing abundance of lesions, destruction of cytoplasmic organelles probably occurs. A positive correlation between the expression of Sp‐Hypp_8991 and digestive lesions may be related to activities to repair nuclear damage incurred during treatment and captive holding.
4.2. Implications for future research
Although there has been some limited application in marine mollusks, and one study has identified candidate stress‐response genes in the freshwater mussel
RNAseq‐based assessment of global gene expression identified candidate markers for quantitative polymerase chain reaction (qPCR) assays, which could be developed and validated for rapidly assessing single‐gene responses to exposure to toxic compounds. As noted above, Hamadeh et al.  demonstrated that chemicals from the same class of compounds give rise to discernible gene expression profiles that bear more similarity to each other than to patterns corresponding to exposure to compounds from a different class . It would be useful for a comprehensive database for gene expression responses to environmental pollutants to be developed toward recognizing compound‐specific expression profiles and molecular signatures of stress responses. Additional, well‐controlled laboratory studies would be appropriate for the purpose of defining defensible biomarkers for toxicant‐induced harm in freshwater mussels. Ultimately, practical work would apply (modifications of) the respective approaches to a range of sites in the Clinch, Powell, and other aquatic ecosystems.
This work was supported by a grant from the US Fish and Wildlife Service. The participation of EMH in this project was supported in part by the Virginia Agricultural Experiment Station under the US Department of Agriculture Hatch Program. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the US Fish and Wildlife Service.