Nominal and mean (standard deviation) measured concentrations of SO42- and As in treatments (n = 4 each) and controls (n = 8) during the 28‐day exposure.
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 Venerupis philippinarum ). There are no existing microarrays for freshwater mussels, however, and construction of such screening platforms would require a large body of work developing and characterizing expressed sequence tags for hundreds to thousands of genes. Subtractive hybridization has been applied to other mollusks, such as the invasive zebra mussel Dreissena polymorpha , Pacific oyster Crassostrea gigas , and peppery furrow shell Scrobicularia plana , and has been shown a viable but technically challenging approach to identifying differentially expressed genes. Next‐generation DNA‐sequencing technology has made it possible to cost‐effectively screen expression of all genes transcribed in tissues of interest, even in species for which prior knowledge of the genome is lacking. The approach is tantamount to sequencing all the RNAs produced in that tissue; hence, the approach is referred to as RNAseq. In our context, RNAseq makes quantitative comparison of gene expression in selected tissues among toxin‐challenged and control individuals possible. Against this background, the goal of this study was to relate gene expression end points in a representative freshwater mussel, pheasantshell A. pectorosa, to exposure to arsenate and sulfate, two pollutants resulting from coal combustion and mining, respectively.
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 Corbicula fluminea is 20,740 µg/L As, with no mortality observed at concentrations up to 5000 µg/L for 21 days . In surface waters, the majority of total As will be As(V) under oxygenated conditions. However, investigations of sublethal effects of As on bivalves have generally focused on As(III). Exposures of bivalves to As(III) have demonstrated histological effects including increased damage to digestive gland tissue  and biochemical effects including alteration of levels of adenosine triphosphate , inhibition of the detoxification enzyme catalase , and changes in the activity of glutathione‐S‐transferase (stimulation and inhibition) [14, 16]. Detoxification of As(III) and As(V) involves reduced glutathione and glutathione‐dependent enzymes, and effects on this system in bivalves have been shown to be variable and dependent on both concentration and chemical speciation [14, 16]. The sublethal effects of environmentally relevant concentrations of As(V) on freshwater mussels are currently unknown.
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 Sphaerium simile of 2078 mg/L. However, exposure of C. fluminea to 1500 mg/L SO42- reduced feeding rates over a 4‐week period .
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 (S. simile and C. fluminea; [13, 17, 18, 21]). The majority of toxicity tests with As have been conducted with arsenite (As(III)), the more toxic form of As. However, due to the oxygenated environment in our exposure system, mussels were exposed to As as arsenate (As(V)). Hence, the purpose of this study was to determine whether As (as As(V)) and SO42- caused sublethal changes in freshwater mussels, the nature of these changes, and the relationships between biochemical, histological, and genetic markers.
2.1. Mussel collection and holding
Pheasantshell mussels (A. pectorosa) were collected from the Clinch River, Hancock County, Tennessee, on September 23, 2014. A total of 143 individuals ranging in size from 70 to 90 mm were collected by hand using snorkel and mask. Mussels were held for 14 days of acclimation at the Freshwater Mollusk Conservation Center (FMCC), Virginia Tech, Blacksburg, Virginia, in flow‐through systems continuously replenished with unfiltered water from a man‐made, on‐site pond.
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 (n = 1) were used to test acute toxicity of five SO42- concentrations: 1500, 1000, 500, 250, and 0 mg/L (control). Treatments were prepared via direct addition of appropriate volumes of a 607‐mM (58.3‐g/L) stock solution of SO42-. Five buckets were used to test acute toxicity of five As(V) concentrations: 2000, 985, 515, 170, and 0 µg/L (control). Treatments were prepared via direct addition of the 100‐mM stock solution of As(V). On Day 7, all 30 mussels were alive; there was no mortality in any of the As(V) or SO42- treatments or controls.
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 (n = 4 for each treatment/control). Because of the availability of experimental units and mussels, two separate controls were used (control_1 and control_2); both consisted of 100% unfiltered pond water. The As(V) treatment (hereafter, HA) concentration of 1000 µg/L was achieved by adding appropriate volumes of 100 mM stock solution to unfiltered pond water. The SO42- treatment (hereafter HS) of 1250 mg/L was achieved by adding appropriate concentrations of a 569‐mM (54.65‐g/L) stock solution to unfiltered pond water. The final water volume in all treatment and control buckets was 17.5 L.
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 e‐values greater than 1e-10 or identity less than 50% were discarded. The clean reads were mapped to the reference assembly using Bowtie v. 1.0.0  with parameters set to ‘‐l 25 ‐I 1 ‐X 1000 ‐a ‐m 200.’ RSEM  was used to calculate the gene expression. Differentially expressed genes were calculated using the edgeR  package in R software (http://www.r‐project.org/), and Benjamini‐Hochberg adjusted p‐values less than 0.05 were considered to be significant.
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 p‐values (<9.5 × 10-8) were identified among genes with differential expression between control, HS, and HA. All annotated genes with differential expression between control, HS, and HA females were also selected. For these 48 selected genes, gene expression was correlated with biochemical measurements and histological variables. Digestive glands from nine mussels were included in the analysis of gene expression, but only eight mussels were included in correlations with biochemical measurements of the glutathione antioxidant system; that is, one extract lost during the set of biochemical analyses was one of the mussels randomly selected for analyses of gene expression.
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, p = 0.22). Mussel growth was not quantified in this study due to the use of adult mussels and the short duration of the exposure.
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 (p < 0.05) differential expression (DE) of particular genes was observed among control mussels (CT) and those exposed to either arsenate (HA) or sulfate (HS). Chemical stress was as likely to cause underexpression as it was to cause overexpression of genes relative to levels observed in control mussels. Contrasts also were observed among mussels exposed to the respective pollutants, indicating that pollutants induced up‐ or downexpression of different suites of genes. From 50 to 100 differentially expressed genes were found for each comparison. Figures 1–3 compare levels of expression of particular genes among groups of mussels.
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 Hyriopsis cumingii.
Among annotated underexpressed genes in the arsenate treatment relative to controls were guanylate binding protein 1 (GBP1) and poly [ADP‐ribose] polymerase (BRAFLDRAFT_74879). In mammalian cells, GBP1 has been shown to respond to exposure to multiple stress agents including paclitaxel and doxorubicin . The enzymes poly(ADP‐ribose) (PAR) and polymerase‐1 (PARP‐1) are central for cellular stress responses and for directing cells to specific fates (e.g., DNA repair vs. cell death) . Among annotated overexpressed genes in arsenate relative to control were poly [ADP‐ribose] polymerase 15‐like (LOC101731886), and sodium‐ and chloride‐dependent glycine transporter 2‐like (LOC100899820).
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 (r, p < 0.05) were observed among genetics and histological data from the nine sampled females. Fractions of reproductive acini containing mature or developing oocytes were significantly negatively correlated with gene identification code (gene code) c140332_g1 (gene symbol = LOC101731886) (r = -0.78, p = 0.01). Fractions of digestive diverticula cells with lesions were significantly correlated with gene code c145102_g1 (none) (r = 0.74, p = 0.04) and negatively correlated with gene codes c103784_g1 (CGI_10002926) (r = -0.76, p = 0.03), c149150_g3 (none) (r = -0.78, p = 0.02), and c152424_g6 (none) (r = -0.73, p = 0.04). Fractions of kidney cells containing lipofuscin inclusions were significantly correlated with c138105_g1 (none) (r = 0.76, p = 0.02).
There were strong correlations (r > 0.85, p < 0.01) between the expression of five genes and biochemical measurements. Numerous weaker correlations (0.05 > p > 0.01) are not presented or discussed. The expression of one annotated gene, c150857_g2 (GBP1), was positively correlated with the activity of glutathione‐S‐transferase, but there was little separation between treatments and control_2 for both gene expression and enzyme activity. The expression of one annotated gene (c152797_g1; BRAFLDRAFT_74879) and one unannotated gene (c126914_g1) was positively correlated with the concentration of reduced glutathione, and there was good separation between treatments and control_2. The expression of two unannotated genes (c155139_g1 and c155139_g2) was negatively correlated with the concentration of reduced glutathione and there was good separation between the treatments and control_2. Correlation between the expression of a sixth gene (unannotated; c154720_g4) did not meet the α = 0.01 criterion for strong correlation, but the plot demonstrated a similar pattern as the other two genes negatively correlated with reduced glutathione concentration, with good separation between treatments and control_2. For the eight mussels included in these correlations, concentrations of reduced glutathione were higher in the control and lower in the two treatments, whether these genes are directly related to reduced glutathione production or utilization is unknown. The small sample of mussels utilized for genetic analysis limits the interpretation of the relationships but suggests that these genes warrant further investigation, particularly as related to As(V) exposure.
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) (PARP) . The PARP enzyme family is associated with physiological functions during cell division, transcriptional regulation, repair of DNA damage, and cytoplasmic cell response . Females that provided the data for the correlations had completed oogenesis, and oocytes in acini were mostly atretic or resorbing; it seems logical that nucleus‐related activities controlled by PARP in postgametogenic females would be negatively related to abundances of developing oocytes, since active gametogenesis had ceased. The positive relationship between abundances of lipofuscin in kidney cells and upregulation of gene LOC101731886 (expression of PARP) may be related to oxidative stress. Increased production of reactive oxygen species (ROS) has been related to DNA single‐strand damage incurred during contaminant exposure and inadequate aquaculture holding conditions [43, 44]. Since abundances of intra‐lysosomal lipofuscin granules are related to both contaminant exposure and oxidative stress [45, 46], it is logical that they are positively related to the expression of PARP activity to repair DNA.
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 Elliptio complanata , RNAseq has not been evaluated for its utility for monitoring of the responses of freshwater mussels to coal‐related environmental contaminants. Successful execution of this project provides proof of principle for using RNAseq technology to approach issues of toxicogenomics in freshwater mussels. Among key findings, our results collectively support the view that substantial changes of gene expression occur before dramatic changes in biochemical and histological effects. Our understanding of how changes in gene expression result in physiological and histological effects on the organism is hampered by the current situation in which many freshwater mollusk genes remain unannotated.
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.